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 $
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