Mercurial > hg > ltpda
view m-toolbox/classes/@ao/computeperiodogram.m @ 33:5e7477b94d94 database-connection-manager
Add known repositories list to LTPDAPreferences
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
line source
% COMPUTEPERIODOGRAM Periodogram spectral estimation. % This function is used to calculate the Power Spectrum Sxx, and the % Cross Power Spectrum Sxy. % % Sxx = COMPUTEPERIODOGRAM(X,WIN,NFFT) where x is a vector returns the % Power Spectrum over the whole Nyquist interval, [0, 2pi). % % Sxy = COMPUTEPERIODOGRAM({X,Y},WIN,NFFT) returns the Cross Power % Spectrum over the whole Nyquist interval, [0, 2pi). % % Inputs: % X - Signal vector or a cell array of two elements containing % two signal vectors. % WIN - Window % NFFT - Number of frequency points (FFT) or vector of % frequencies at whicch periodogram is desired (Goertzel) % WINCOMPFLAG - A string indicating the type of window compensation to % be done. The choices are: % 'ms' - compensate for Mean-square (Power) Spectrum; % maintain the correct power peak heights. % 'psd' - compensate for Power Spectral Denstiy (PSD); % maintain correct area under the PSD curve. % % Output: % Sxx - Power spectrum [Power] over the whole Nyquist interval. % or % Sxy - Cross power spectrum [Power] over the whole Nyquist % interval. % % A direct copy of MATLAB's function for LTPDA % % M Hewitson 08-05-08 % % $Id: computeperiodogram.m,v 1.2 2008/08/01 13:19:42 ingo Exp $ % % Author(s): P. Pacheco % Copyright 1988-2006 The MathWorks, Inc. % $Revision: 1.2 $ $Date: 2008/08/01 13:19:42 $ function [P,f] = computeperiodogram(x,win,nfft,esttype,varargin) error(nargchk(3,5,nargin,'struct')); if nargin < 4, esttype = 'psd'; % Default, compenstate for window's power. end if nargin < 5 || isempty(varargin{1}), Fs = 2*pi; else Fs = varargin{1}; end % Validate inputs and convert row vectors to column vectors. [x,Lx,y,is2sig,win,errid,errmsg] = validateinputs(x,win,nfft); if ~isempty(errmsg), error(errid, errmsg); end % Determine if FFT or Goertzel freqVectorSpecified = false; if (length(nfft)> 1), freqVectorSpecified = true; end % Window the data xw = x.*win; if is2sig, yw = y.*win; end % Evaluate the window normalization constant. A 1/N factor has been % omitted since it will cancel below. if strcmpi(esttype,'ms'), % The window is convolved with every power spectrum peak, therefore % compensate for the DC value squared to obtain correct peak heights. U = sum(win)^2; else U = win'*win; % compensates for the power of the window. end % Compute the periodogram power spectrum [Power] estimate % A 1/N factor has been omitted since it cancels [Xx,f] = ao.computeDFT(xw,nfft,Fs); if is2sig, [Yy,f] = ao.computeDFT(yw,nfft,Fs); end P = Xx.*conj(Xx)/U; % Auto spectrum. if is2sig, P = Xx.*conj(Yy)/U; % Cross spectrum. end end %-------------------------------------------------------------------------- function [x,Lx,y,is2sig,win,errid,errmsg] = validateinputs(x,win,nfft) % Validate the inputs to computexperiodogram and convert the input signal % vectors into column vectors. errid = ''; errmsg = ''; % in case of early return. % Set defaults and convert to row vectors to columns. y = []; Ly = 0; is2sig= false; win = win(:); Lw = length(win); % Determine if one or two signal vectors was specified. Lx = length(x); if iscell(x), if length(x) > 1, y = x{2}; y = y(:); is2sig = true; end x = x{1}; Lx = length(x); end x = x(:); if is2sig, Ly = length(y); if Lx ~= Ly, errid = generatemsgid('invalidInputSignalLength'); errmsg = 'The length of the two input vectors must be equal to calculate the cross spectral density.'; return; end end if Lx ~= Lw, errid = generatemsgid('invalidWindow'); errmsg = 'The WINDOW must be a vector of the same length as the signal.'; return; end end