Home > classes > @ao > computeperiodogram.m

computeperiodogram

PURPOSE ^

COMPUTEPERIODOGRAM Periodogram spectral estimation.

SYNOPSIS ^

function [P,f] = computeperiodogram(x,win,nfft,esttype,varargin)

DESCRIPTION ^

 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 $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % COMPUTEPERIODOGRAM   Periodogram spectral estimation.
0002 %   This function is used to calculate the Power Spectrum Sxx, and the
0003 %   Cross Power Spectrum Sxy.
0004 %
0005 %   Sxx = COMPUTEPERIODOGRAM(X,WIN,NFFT) where x is a vector returns the
0006 %   Power Spectrum over the whole Nyquist interval, [0, 2pi).
0007 %
0008 %   Sxy = COMPUTEPERIODOGRAM({X,Y},WIN,NFFT) returns the Cross Power
0009 %   Spectrum over the whole Nyquist interval, [0, 2pi).
0010 %
0011 %   Inputs:
0012 %    X           - Signal vector or a cell array of two elements containing
0013 %                  two signal vectors.
0014 %    WIN         - Window
0015 %    NFFT        - Number of frequency points (FFT) or vector of
0016 %    frequencies at whicch periodogram is desired (Goertzel)
0017 %    WINCOMPFLAG - A string indicating the type of window compensation to
0018 %                  be done. The choices are:
0019 %                  'ms'  - compensate for Mean-square (Power) Spectrum;
0020 %                          maintain the correct power peak heights.
0021 %                  'psd' - compensate for Power Spectral Denstiy (PSD);
0022 %                          maintain correct area under the PSD curve.
0023 %
0024 %   Output:
0025 %    Sxx         - Power spectrum [Power] over the whole Nyquist interval.
0026 %      or
0027 %    Sxy         - Cross power spectrum [Power] over the whole Nyquist
0028 %                  interval.
0029 %
0030 % A direct copy of MATLAB's function for LTPDA
0031 %
0032 % M Hewitson 08-05-08
0033 %
0034 % $Id: computeperiodogram.m,v 1.2 2008/08/01 13:19:42 ingo Exp $
0035 %
0036 
0037 %   Author(s): P. Pacheco
0038 %   Copyright 1988-2006 The MathWorks, Inc.
0039 %   $Revision: 1.2 $  $Date: 2008/08/01 13:19:42 $
0040 
0041 function [P,f] = computeperiodogram(x,win,nfft,esttype,varargin)
0042 
0043   error(nargchk(3,5,nargin,'struct'));
0044   if nargin < 4,
0045     esttype = 'psd'; % Default, compenstate for window's power.
0046   end
0047 
0048   if nargin < 5 || isempty(varargin{1}),
0049     Fs = 2*pi;
0050   else
0051     Fs = varargin{1};
0052   end
0053 
0054   % Validate inputs and convert row vectors to column vectors.
0055   [x,Lx,y,is2sig,win,errid,errmsg] = validateinputs(x,win,nfft);
0056   if ~isempty(errmsg), error(errid, errmsg); end
0057 
0058   % Determine if FFT or Goertzel
0059   freqVectorSpecified = false;
0060   if (length(nfft)> 1), freqVectorSpecified = true; end
0061 
0062   % Window the data
0063   xw = x.*win;
0064   if is2sig, yw = y.*win; end
0065 
0066   % Evaluate the window normalization constant.  A 1/N factor has been
0067   % omitted since it will cancel below.
0068   if strcmpi(esttype,'ms'),
0069     % The window is convolved with every power spectrum peak, therefore
0070     % compensate for the DC value squared to obtain correct peak heights.
0071     U = sum(win)^2;
0072   else
0073     U = win'*win;  % compensates for the power of the window.
0074   end
0075 
0076   % Compute the periodogram power spectrum [Power] estimate
0077   % A 1/N factor has been omitted since it cancels
0078 
0079   [Xx,f] = ao.computeDFT(xw,nfft,Fs);
0080   if is2sig, [Yy,f] = ao.computeDFT(yw,nfft,Fs); end
0081 
0082   P = Xx.*conj(Xx)/U;      % Auto spectrum.
0083   if is2sig,
0084     P = Xx.*conj(Yy)/U;  % Cross spectrum.
0085   end
0086 
0087 end
0088 
0089 %--------------------------------------------------------------------------
0090 function [x,Lx,y,is2sig,win,errid,errmsg] = validateinputs(x,win,nfft)
0091   % Validate the inputs to computexperiodogram and convert the input signal
0092   % vectors into column vectors.
0093 
0094   errid = '';
0095   errmsg = ''; % in case of early return.
0096 
0097   % Set defaults and convert to row vectors to columns.
0098   y     = [];
0099   Ly    = 0;
0100   is2sig= false;
0101   win   = win(:);
0102   Lw    = length(win);
0103 
0104   % Determine if one or two signal vectors was specified.
0105   Lx = length(x);
0106   if iscell(x),
0107     if length(x) > 1,
0108       y = x{2};
0109       y = y(:);
0110       is2sig = true;
0111     end
0112     x = x{1};
0113     Lx = length(x);
0114   end
0115   x = x(:);
0116 
0117   if is2sig,
0118     Ly  = length(y);
0119     if Lx ~= Ly,
0120       errid = generatemsgid('invalidInputSignalLength');
0121       errmsg = 'The length of the two input vectors must be equal to calculate the cross spectral density.';
0122       return;
0123     end
0124   end
0125 
0126   if Lx ~= Lw,
0127     errid = generatemsgid('invalidWindow');
0128     errmsg = 'The WINDOW must be a vector of the same length as the signal.';
0129     return;
0130   end
0131 
0132 end
0133 
0134

Generated on Mon 08-Sep-2008 13:18:47 by m2html © 2003