Mercurial > octave-dspies
view scripts/signal/periodogram.m @ 18995:8ac4ab4ae5f4
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.
author | Drew Abbot <drewabbot@gmail.com> and Rik <rik@octave.org> |
---|---|
date | Thu, 07 Aug 2014 10:13:30 -0700 |
parents | d63878346099 |
children | 956fc864c39f |
line wrap: on
line source
## 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 ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @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{}) ## ## Return the periodogram (Power Spectral Density) of @var{x}. ## ## The possible inputs are: ## ## @table @var ## @item x ## ## 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. ## ## @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 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 Fs ## sampling rate. The default is 1. ## ## @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 ## ## 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}. ## ## 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 <Friedrich.Leisch@ci.tuwien.ac.at> ## 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 (! 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 (isempty (nfft)) nfft = max (256, 2.^nextpow2 (n)); elseif (! isscalar (nfft)) error ("periodogram: NFFT must be a scalar"); endif if (strcmpi (range, "onesided")) range = 1; 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 (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 ## FIXME: Need some functional tests %% Test input validation %!error periodogram () %!error periodogram (1,2,3,4,5,6) %!error <X must be a real or complex vector> periodogram (ones (2,2)) %!error <WIN must be a vector.*same length> periodogram (1:5, ones (2,2)) %!error <WIN must be a vector.*same length> periodogram (1:5, 1:6) %!error <NFFT must be a scalar> periodogram (1:5, 1:5, 1:5) %!error <"centered" range type is not implemented> periodogram (1:5, "centered")