# HG changeset patch # User Drew Abbot and Rik # Date 1407431610 25200 # Node ID 8ac4ab4ae5f40ff8074cd0df43b0bae260575466 # Parent 79d2dd9351cca15a658d2d68652a5c81b375fc11 periodogram.m: Overhaul function (bug #39279, bug #42859). * contributors.in: Add Drew Abbot to list of contributors. * periodogram.m: Rewrite documentation. Simplify input parsing of arguments. Accept both row and column inputs for X. Correct onesided computation when NFFT is odd. Add an error message about unrecognized range specification "centered". Add input validation tests. diff -r 79d2dd9351cc -r 8ac4ab4ae5f4 doc/interpreter/contributors.in --- a/doc/interpreter/contributors.in Thu Aug 07 11:46:12 2014 -0400 +++ b/doc/interpreter/contributors.in Thu Aug 07 10:13:30 2014 -0700 @@ -1,4 +1,5 @@ Ben Abbott +Drew Abbot Andy Adler Adam H. Aitkenhead Giles Anderson diff -r 79d2dd9351cc -r 8ac4ab4ae5f4 scripts/signal/periodogram.m --- a/scripts/signal/periodogram.m Thu Aug 07 11:46:12 2014 -0400 +++ b/scripts/signal/periodogram.m Thu Aug 07 10:13:30 2014 -0700 @@ -1,5 +1,6 @@ ## Copyright (C) 1995-2013 Friedrich Leisch ## Copyright (C) 2010 Alois Schloegl +## Copyright (C) 2014 Drew Abbot ## ## This file is part of Octave. ## @@ -18,43 +19,60 @@ ## . ## -*- 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}. +## @deftypefn {Function File} {[@var{Pxx}, @var{w}] =} periodogram (@var{x}) +## @deftypefnx {Function File} {[@var{Pxx}, @var{w}] =} periodogram (@var{x}, @var{win}) +## @deftypefnx {Function File} {[@var{Pxx}, @var{w}] =} periodogram (@var{x}, @var{win}, @var{nfft}) +## @deftypefnx {Function File} {[@var{Pxx}, @var{f}] =} periodogram (@var{x}, @var{win}, @var{nfft}, @var{Fs}) +## @deftypefnx {Function File} {[@var{Pxx}, @var{f}] =} periodogram (@dots{}, "@var{range}") +## @deftypefnx {Function File} {} periodogram (@dots{}) ## -## [Pxx,w] = periodogram (@var{x}). +## Return the periodogram (Power Spectral Density) of @var{x}. ## -## [Pxx,w] = periodogram (@var{x},win). +## The possible inputs are: ## -## [Pxx,w] = periodogram (@var{x},win,nfft). -## -## [Pxx,f] = periodogram (@var{x},win,nfft,Fs). +## @table @var +## @item x ## -## [Pxx,f] = periodogram (@var{x},win,nfft,Fs,"range"). +## data vector. If @var{x} is real-valued a one-sided spectrum is estimated. +## If @var{x} is complex-valued, or @qcode{"@var{range}"} specifies +## @qcode{"@nospell{twosided}"}, the full spectrum is estimated. ## -## @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 +## window weight data. If window is empty or unspecified a default rectangular +## window is used. Otherwise, the window is applied to the signal +## (@code{@var{x} .* @var{win}}) before computing the periodogram. The window +## data must be a vector of the same length as @var{x}. ## -## @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. The default is 256 or the next higher power of +## 2 greater than the length of @var{x} +## (@code{max (256, 2.^nextpow2 (length (x)))}). If @var{nfft} is greater +## than the length of the input then @var{x} will be zero-padded to the length +## of @var{nfft}. ## -## @item nfft: number of frequency bins, default max (256, 2.^ceil (log2 (length (x)))). -## -## @item Fs: sampling rate, default 1. +## @item Fs +## sampling rate. The default is 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 range +## range of spectrum. @qcode{"@nospell{onesided}"} computes spectrum from +## [0..nfft/2+1]. @qcode{"@nospell{twosided}"} computes spectrum from +## [0..nfft-1]. +## @end table ## -## @item @nospell{Pxx}: one-, or two-sided power spectrum. +## The optional second output @var{w} are the normalized angular frequencies. +## For a onesided calculation @var{w} is in the range [0, pi] if @var{nfft} +## is even and [0, pi) if @var{nfft} is odd. Similarly, for a twosided +## calculation @var{w} is in the range [0, 2*pi] or [0, 2*pi) depending on +## @var{nfft}. ## -## @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 +## If a sampling frequency is specified, @var{Fs}, then the output frequencies +## @var{f} will be in the range [0, @var{Fs}/2] or [0, @var{Fs}/2) for +## onesided calculations. For twosided calculations the range will be +## [0, @var{Fs}). +## +## When called with no outputs the periodogram is immediately plotted in the +## current figure window. +## @seealso{fft} ## @end deftypefn ## Author: FL @@ -63,12 +81,11 @@ function [pxx, f] = periodogram (x, varargin) ## check input arguments - if (nargin < 1 || nargin > 5) print_usage (); endif - nfft = []; fs = []; range = []; window = []; + nfft = fs = range = window = []; j = 1; for k = 1:length (varargin) if (ischar (varargin{k})) @@ -88,66 +105,54 @@ endif endfor - [r, c] = size (x); - if (r == 1) - r = c; + if (! isvector (x)) + error ("periodogram: X must be a real or complex vector"); + endif + x = x(:); # Use column vectors from now on + + n = rows (x); + + if (! isempty (window)) + if (! isvector (window) || length (window) != n) + error ("periodogram: WIN must be a vector of the same length as X"); + endif + window = window(:); + x .*= window; endif - if (ischar (window)) - range = window; - window = []; - endif; - if (ischar (nfft)) - range = nfft; - nfft = []; - endif; - if (ischar (fs)) - range = fs; - fs = []; - endif; - - if (! isempty (window)) - if (all (size (x) == size (window))) - x .*= window; - elseif (rows (x) == rows (window) && columns (window) == 1) - x .*= window (:,ones (1,c)); - endif; + if (isempty (nfft)) + nfft = max (256, 2.^nextpow2 (n)); + elseif (! isscalar (nfft)) + error ("periodogram: NFFT must be a scalar"); endif - if (numel (nfft)>1) - error ("nfft must be scalar"); - endif - if (isempty (nfft)) - nfft = max (256, 2.^ceil (log2 (r))); - endif - - if (strcmp (range, "onesided")) + if (strcmpi (range, "onesided")) range = 1; - elseif (strcmp (range, "twosided")) + elseif (strcmpi (range, "twosided")) range = 2; + elseif (strcmpi (range, "centered")) + error ('periodogram: "centered" range type is not implemented'); else range = 2-isreal (x); endif ## compute periodogram - if (r>nfft) + if (n > nfft) Pxx = 0; rr = rem (length (x), nfft); if (rr) - x = [x(:); (zeros (nfft-rr, 1))]; + x = [x(:); zeros(nfft-rr, 1)]; endif x = sum (reshape (x, nfft, []), 2); endif - if (isempty (window)) - n = r; - else + if (! isempty (window)) n = sumsq (window); - end; - Pxx = (abs (fft (x, nfft))) .^ 2 / n ; + endif; + Pxx = (abs (fft (x, nfft))) .^ 2 / n; - if (nargin<4) + if (nargin < 4) Pxx /= 2*pi; elseif (! isempty (fs)) Pxx /= fs; @@ -156,24 +161,30 @@ ## generate output arguments if (range == 1) # onesided - Pxx = Pxx(1:nfft/2+1) + [0; Pxx(end:-1:(nfft/2+2)); 0]; + 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; + f = (0:nfft/2)' / nfft; elseif (range == 2) - f = (0:nfft-1)'/nfft; + f = (0:nfft-1)' / nfft; endif - if (nargin<4) - f *= 2*pi; # generate w=2*pi*f + if (nargin < 4) + f *= 2*pi; # generate w=2*pi*f elseif (! isempty (fs)) f *= fs; endif endif if (nargout == 0) - if (nargin<4) + if (nargin < 4) plot (f/(2*pi), 10*log10 (Pxx)); xlabel ("normalized frequency [x pi rad]"); ylabel ("Power density [dB/rad/sample]"); @@ -190,3 +201,16 @@ endfunction + +## FIXME: Need some functional tests + + +%% Test input validation +%!error periodogram () +%!error periodogram (1,2,3,4,5,6) +%!error periodogram (ones (2,2)) +%!error periodogram (1:5, ones (2,2)) +%!error periodogram (1:5, 1:6) +%!error periodogram (1:5, 1:5, 1:5) +%!error <"centered" range type is not implemented> periodogram (1:5, "centered") +