Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/computeperiodogram.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % COMPUTEPERIODOGRAM Periodogram spectral estimation.
2 % This function is used to calculate the Power Spectrum Sxx, and the
3 % Cross Power Spectrum Sxy.
4 %
5 % Sxx = COMPUTEPERIODOGRAM(X,WIN,NFFT) where x is a vector returns the
6 % Power Spectrum over the whole Nyquist interval, [0, 2pi).
7 %
8 % Sxy = COMPUTEPERIODOGRAM({X,Y},WIN,NFFT) returns the Cross Power
9 % Spectrum over the whole Nyquist interval, [0, 2pi).
10 %
11 % Inputs:
12 % X - Signal vector or a cell array of two elements containing
13 % two signal vectors.
14 % WIN - Window
15 % NFFT - Number of frequency points (FFT) or vector of
16 % frequencies at whicch periodogram is desired (Goertzel)
17 % WINCOMPFLAG - A string indicating the type of window compensation to
18 % be done. The choices are:
19 % 'ms' - compensate for Mean-square (Power) Spectrum;
20 % maintain the correct power peak heights.
21 % 'psd' - compensate for Power Spectral Denstiy (PSD);
22 % maintain correct area under the PSD curve.
23 %
24 % Output:
25 % Sxx - Power spectrum [Power] over the whole Nyquist interval.
26 % or
27 % Sxy - Cross power spectrum [Power] over the whole Nyquist
28 % interval.
29 %
30 % A direct copy of MATLAB's function for LTPDA
31 %
32 % M Hewitson 08-05-08
33 %
34 % $Id: computeperiodogram.m,v 1.2 2008/08/01 13:19:42 ingo Exp $
35 %
36
37 % Author(s): P. Pacheco
38 % Copyright 1988-2006 The MathWorks, Inc.
39 % $Revision: 1.2 $ $Date: 2008/08/01 13:19:42 $
40
41 function [P,f] = computeperiodogram(x,win,nfft,esttype,varargin)
42
43 error(nargchk(3,5,nargin,'struct'));
44 if nargin < 4,
45 esttype = 'psd'; % Default, compenstate for window's power.
46 end
47
48 if nargin < 5 || isempty(varargin{1}),
49 Fs = 2*pi;
50 else
51 Fs = varargin{1};
52 end
53
54 % Validate inputs and convert row vectors to column vectors.
55 [x,Lx,y,is2sig,win,errid,errmsg] = validateinputs(x,win,nfft);
56 if ~isempty(errmsg), error(errid, errmsg); end
57
58 % Determine if FFT or Goertzel
59 freqVectorSpecified = false;
60 if (length(nfft)> 1), freqVectorSpecified = true; end
61
62 % Window the data
63 xw = x.*win;
64 if is2sig, yw = y.*win; end
65
66 % Evaluate the window normalization constant. A 1/N factor has been
67 % omitted since it will cancel below.
68 if strcmpi(esttype,'ms'),
69 % The window is convolved with every power spectrum peak, therefore
70 % compensate for the DC value squared to obtain correct peak heights.
71 U = sum(win)^2;
72 else
73 U = win'*win; % compensates for the power of the window.
74 end
75
76 % Compute the periodogram power spectrum [Power] estimate
77 % A 1/N factor has been omitted since it cancels
78
79 [Xx,f] = ao.computeDFT(xw,nfft,Fs);
80 if is2sig, [Yy,f] = ao.computeDFT(yw,nfft,Fs); end
81
82 P = Xx.*conj(Xx)/U; % Auto spectrum.
83 if is2sig,
84 P = Xx.*conj(Yy)/U; % Cross spectrum.
85 end
86
87 end
88
89 %--------------------------------------------------------------------------
90 function [x,Lx,y,is2sig,win,errid,errmsg] = validateinputs(x,win,nfft)
91 % Validate the inputs to computexperiodogram and convert the input signal
92 % vectors into column vectors.
93
94 errid = '';
95 errmsg = ''; % in case of early return.
96
97 % Set defaults and convert to row vectors to columns.
98 y = [];
99 Ly = 0;
100 is2sig= false;
101 win = win(:);
102 Lw = length(win);
103
104 % Determine if one or two signal vectors was specified.
105 Lx = length(x);
106 if iscell(x),
107 if length(x) > 1,
108 y = x{2};
109 y = y(:);
110 is2sig = true;
111 end
112 x = x{1};
113 Lx = length(x);
114 end
115 x = x(:);
116
117 if is2sig,
118 Ly = length(y);
119 if Lx ~= Ly,
120 errid = generatemsgid('invalidInputSignalLength');
121 errmsg = 'The length of the two input vectors must be equal to calculate the cross spectral density.';
122 return;
123 end
124 end
125
126 if Lx ~= Lw,
127 errid = generatemsgid('invalidWindow');
128 errmsg = 'The WINDOW must be a vector of the same length as the signal.';
129 return;
130 end
131
132 end
133
134