## Copyright (C) 1995-2013 Friedrich Leisch
## Copyright (C) 2010 Alois Schloegl
## Copyright (C) 2014 Drew Abbot
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING. If not, see
## .
## -*- texinfo -*-
## @deftypefn {Function File} {[Pxx, @var{w}] =} periodogram (@var{x})
## For a data matrix @var{x} from a sample of size @var{n}, return the
## periodogram. The angular frequency is returned in @var{w}.
##
## [Pxx,w] = periodogram (@var{x}).
##
## [Pxx,w] = periodogram (@var{x},win).
##
## [Pxx,w] = periodogram (@var{x},win,nfft).
##
## [Pxx,f] = periodogram (@var{x},win,nfft,Fs).
##
## [Pxx,f] = periodogram (@var{x},win,nfft,Fs,"range").
##
## @itemize
## @item x: data; if real-valued a one-sided spectrum is estimated,
## if complex-valued or range indicates @qcode{"@nospell{twosided}"}, the full
## spectrum is estimated.
##
## @item win: weight data with window, x.*win is used for further computation,
## if window is empty, a rectangular window is used.
##
## @item nfft: number of frequency bins, default max (256, 2.^ceil (log2 (length (x)))).
##
## @item Fs: sampling rate, default 1.
##
## @item range: @qcode{"@nospell{onesided}"} computes spectrum from [0..nfft/2+1].
## @qcode{"@nospell{twosided}"} computes spectrum from [0..nfft-1]. These
## strings can appear at any position in the list input arguments after
## window.
##
## @item @nospell{Pxx}: one-, or two-sided power spectrum.
##
## @item w: angular frequency [0..2*pi) (two-sided) or [0..pi] one-sided.
##
## @item f: frequency [0..Fs) (two-sided) or [0..Fs/2] one-sided.
## @end itemize
## @end deftypefn
## Author: FL
## Description: Compute the periodogram
function [pxx, f] = periodogram (x, varargin)
## check input arguments
if (nargin < 1 || nargin > 5)
print_usage ();
endif
nfft = []; fs = []; range = []; window = [];
j = 1;
for k = 1:length (varargin)
if (ischar (varargin{k}))
range = varargin{k};
else
switch (j)
case 1
window = varargin{k};
case 2
nfft = varargin{k};
case 3
fs = varargin{k};
case 4
range = varargin{k};
endswitch
j++;
endif
endfor
if (ischar (window))
range = window;
window = [];
endif;
if (ischar (nfft))
range = nfft;
nfft = [];
endif;
if (ischar (fs))
range = fs;
fs = [];
endif;
if (isrow (x))
x = x.';
endif
n = size (x,1);
if (! isempty (window))
if (isrow (window))
window = window.';
endif
if ( size(window,1) != n || ~isvector(window) )
error('periodogram: the window must be a vector of the same length as the signal.')
endif
x .*= window;
endif
if (numel (nfft)>1)
error ("nfft must be scalar");
endif
if (isempty (nfft))
nfft = max (256, 2.^ceil (log2 (n)));
endif
if (strcmp (range, "onesided"))
range = 1;
elseif (strcmp (range, "twosided"))
range = 2;
else
range = 2-isreal (x);
endif
## compute periodogram
if (n>nfft)
Pxx = 0;
rr = rem (length (x), nfft);
if (rr)
x = [x(:); (zeros (nfft-rr, 1))];
endif
x = sum (reshape (x, nfft, []), 2);
endif
if (! isempty (window))
n = sumsq (window);
endif;
Pxx = (abs (fft (x, nfft))) .^ 2 / n ;
if (nargin<4)
Pxx /= 2*pi;
elseif (! isempty (fs))
Pxx /= fs;
endif
## generate output arguments
if (range == 1) # onesided
if ( ~ rem(nfft,2) ) # nfft is even
psd_len = nfft/2+1;
Pxx = Pxx(1:psd_len) + [0; Pxx(nfft:-1:psd_len+1); 0];
else # nfft is odd
psd_len = (nfft+1)/2;
Pxx = Pxx(1:psd_len) + [0; Pxx(nfft:-1:psd_len+1)];
endif
endif
if (nargout != 1)
if (range == 1)
f = (0:nfft/2)'/nfft;
elseif (range == 2)
f = (0:nfft-1)'/nfft;
endif
if (nargin<4)
f *= 2*pi; # generate w=2*pi*f
elseif (! isempty (fs))
f *= fs;
endif
endif
if (nargout == 0)
if (nargin<4)
plot (f/(2*pi), 10*log10 (Pxx));
xlabel ("normalized frequency [x pi rad]");
ylabel ("Power density [dB/rad/sample]");
else
plot (f, 10*log10 (Pxx));
xlabel ("frequency [Hz]");
ylabel ("Power density [dB/Hz]");
endif
grid on;
title ("Periodogram Power Spectral Density Estimate");
else
pxx = Pxx;
endif
endfunction