Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/wosa.m @ 43:bc767aaa99a8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
CVS Update
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Tue, 06 Dec 2011 11:09:25 +0100
parents
children
comparison
equal deleted inserted replaced
0:f0afece42f48 43:bc767aaa99a8
1 % WOSA implements Welch's overlaped segmented averaging algorithm with
2 % segment detrending and variance estimation.
3 %
4 % [pxx, f, info] = wosa(x,type,pl)
5 % [pxx, f, info] = wosa(x,y,type,pl)
6 %
7 % INPUTS: x - input analysis objects
8 % y - input analysis objects
9 % type - type of estimation:
10 % 'psd' - compute Power Spectral Denstiy (PSD)
11 % 'cpsd' - compute cross-spectral density
12 % 'tfe' - estimate transfer function between inputs
13 % 'mscohere' - estimate magnitude-squared cross-coherence
14 % 'cohere' - estimate complex cross-coherence
15 % pl - input parameter list
16 %
17 % PARAMETERS: 'Win' - a specwin window object [default: Kaiser -200dB psll]
18 % 'Olap' - segment percent overlap [default: taken from window function]
19 % 'Nfft' - number of samples in each fft [default: length of input data]
20 % 'Scale' - one of
21 % 'ASD' - amplitude spectral density
22 % 'PSD' - power spectral density [default]
23 % 'AS' - amplitude spectrum
24 % 'PS' - power spectrum
25 % * applies only to spectrum 'Type' 'psd'
26 % 'Order' - order of segment detrending
27 % -1 - no detrending
28 % 0 - subtract mean [default]
29 % 1 - subtract linear fit
30 % N - subtract fit of polynomial, order N
31 %
32 % Version: $Id: wosa.m,v 1.5 2011/12/02 07:08:11 hewitson Exp $
33 %
34
35 function varargout = wosa(varargin)
36 import utils.const.*
37
38 % Process inputs
39 if nargin == 3
40 a = varargin{1};
41 esttype = varargin{2};
42 pl = varargin{3};
43 inunits = a.data.yunits;
44 L = a.len;
45 else
46 a = varargin{1};
47 b = varargin{2};
48 esttype = varargin{3};
49 pl = varargin{4};
50 if a.data.fs ~= b.data.fs
51 error('The two time-series have different sample rates.');
52 end
53 inunits = b.data.yunits / a.data.yunits;
54 L = min(a.len, b.len);
55 end
56
57 % Parse inputs
58 win = find(pl, 'Win');
59 nfft = find(pl, 'Nfft');
60 percentOlap = find(pl, 'Olap')/100;
61 scale = find(pl, 'scale');
62 xOlap = round(percentOlap*nfft);
63 detrendOrder = find(pl, 'order');
64 fs = a.fs;
65 winVals = win.win.'; % because we always get a column from ao.y
66
67 % Compute segment details
68
69 nSegments = fix((L - xOlap)./(nfft - xOlap));
70 utils.helper.msg(msg.PROC3, 'N segment: %d', nfft);
71 utils.helper.msg(msg.PROC3, 'N overlap: %d', xOlap);
72 utils.helper.msg(msg.PROC3, 'N segments: %d', nSegments);
73
74 % Compute start and end indices of each segment
75 segmentStep = nfft-xOlap;
76 segmentStarts = 1:segmentStep:nSegments*segmentStep;
77 segmentEnds = segmentStarts+nfft-1;
78
79 % Estimate the averaged periodogram for the desired quantity
80 switch lower(esttype)
81 case 'psd'
82 % Compute averaged periodogram
83 [Sxx, Svxx] = psdPeriodogram(a, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder);
84 case 'cpsd'
85 [Sxx, Svxx] = cpsdPeriodogram(a, b, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder);
86 case 'tfe'
87 [Sxx, Sxy, Syy] = tfePeriodogram(a, b, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder);
88 case {'mscohere','cohere'}
89 [Sxx, Sxy, Syy] = tfePeriodogram(a, b, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder);
90 otherwise
91 error('Unknown estimation type %s', esttype);
92 end
93
94 % Scale to PSD
95 switch lower(esttype)
96 case {'psd','cpsd'}
97 [P, Pvxx] = scaleToPSD(Sxx, Svxx, nfft, fs);
98 % the 1/nSegments factor should come after welchscale if we don't
99 % want to apply sqrt() to it.
100 % We correct for that here. It is only needed for 'asd','as' in
101 % psd/cpsd, the other cases go always through 'PSD'.
102 if (strcmpi(scale,'PSD') || strcmpi(scale,'PS'))
103 dP = Pvxx;
104 elseif (strcmpi(scale,'ASD') || strcmpi(scale,'AS'))
105 dP = Pvxx/nSegments;
106 else
107 error('### Unknown scale')
108 end
109 case 'tfe'
110 % Compute the 1-sided or 2-sided PSD [Power/freq] or mean-square [Power].
111 % Also, corresponding freq vector and freq units.
112 % In the Cross PSD, the frequency vector and xunits are not used.
113 Pxx = scaleToPSD(Sxx, [], nfft, fs);
114 Pxy = scaleToPSD(Sxy, [], nfft, fs);
115 Pyy = scaleToPSD(Syy, [], nfft, fs);
116 % mean and std
117 P = Pxy ./ Pxx; % Txy
118 if nSegments == 1
119 dP =[];
120 else
121 dP = (nSegments/(nSegments-1)^2)*(Pyy./Pxx).*(1 - (abs(Pxy).^2)./(Pxx.*Pyy));
122 end
123 case 'mscohere'
124 % Magnitude Square Coherence estimate.
125 % Auto PSD for 2nd input vector. The freq vector & xunits are not
126 % used.
127 Pxx = scaleToPSD(Sxx, [], nfft, fs);
128 Pxy = scaleToPSD(Sxy, [], nfft, fs);
129 Pyy = scaleToPSD(Syy, [], nfft, fs);
130 % mean and std
131 P = (abs(Pxy).^2)./(Pxx.*Pyy); % Magnitude-squared coherence
132 dP = (2*P/nSegments).*(1-P).^2;
133 case 'cohere'
134 % Complex Coherence estimate.
135 % Auto PSD for 2nd input vector. The freq vector & xunits are not
136 % used.
137 Pxx = scaleToPSD(Sxx, [], nfft, fs);
138 Pxy = scaleToPSD(Sxy, [], nfft, fs);
139 Pyy = scaleToPSD(Syy, [], nfft, fs);
140 P = Pxy./sqrt(Pxx.*Pyy); % Complex coherence
141 dP = (2*abs(P)/nSegments).*(1-abs(P)).^2;
142
143 end
144
145 % Compute frequencies
146 freqs = psdfreqvec('npts', nfft,'Fs', fs, 'Range', 'half').';
147
148 % Scale to required units
149 [Pxx, dP, info] = utils.math.welchscale(P, dP, winVals, fs, scale, inunits);
150 info.navs = nSegments;
151
152 if nSegments ==1
153 dev = [];
154 else
155 dev = sqrt(dP);
156 end
157
158 % Set outputs
159 varargout = {Pxx, freqs, info, dev};
160
161 end
162
163 % scale averaged periodogram to PSD
164 function [Pxx, Pvxx] = scaleToPSD(Sxx, Svxx, nfft, fs)
165
166 % Take 1-sided spectrum which means we double the power in the
167 % appropriate bins
168 if rem(nfft,2),
169 indices = 1:(nfft+1)/2; % ODD
170 Sxx1sided = Sxx(indices,:);
171 % double the power except for the DC bin
172 Sxx = [Sxx1sided(1,:); 2*Sxx1sided(2:end,:)];
173 if ~isempty(Svxx)
174 Svxx1sided = Svxx(indices,:);
175 Svxx = [Svxx1sided(1,:); 4*Svxx1sided(2:end,:)];
176 end
177 else
178 indices = 1:nfft/2+1; % EVEN
179 Sxx1sided = Sxx(indices,:);
180 % Double power except the DC bin and the Nyquist bin
181 Sxx = [Sxx1sided(1,:); 2*Sxx1sided(2:end-1,:); Sxx1sided(end,:)];
182 if ~isempty(Svxx)
183 Svxx1sided = Svxx(indices,:); % Take only [0,pi] or [0,pi)
184 Svxx = [Svxx1sided(1,:); 4*Svxx1sided(2:end-1,:); Svxx1sided(end,:)];
185 end
186 end
187
188 % Now scale to PSD
189 Pxx = Sxx./fs;
190 Pvxx = Svxx./fs^2;
191
192 end
193
194 % compute tfe
195 function [Sxx, Sxy, Syy] = tfePeriodogram(x, y, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder)
196
197 nfft = segmentEnds(1);
198 Sxx = zeros(nfft,1); % Initialize Sxx
199 Sxy = zeros(nfft,1); % Initialize Sxy
200 Syy = zeros(nfft,1); % Initialize Syy
201 % loop over segments
202 for ii = 1:nSegments
203 if detrendOrder < 0
204 xseg = x.y(segmentStarts(ii):segmentEnds(ii));
205 yseg = y.y(segmentStarts(ii):segmentEnds(ii));
206 else
207 [xseg,coeffs] = ltpda_polyreg(x.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
208 [yseg,coeffs] = ltpda_polyreg(y.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
209 end
210
211 % Compute periodograms
212 Sxxk = wosa_periodogram(xseg, [], winVals, nfft);
213 Sxyk = wosa_periodogram(yseg, xseg, winVals, nfft);
214 Syyk = wosa_periodogram(yseg, [], winVals, nfft);
215
216 Sxx = Sxx + Sxxk;
217 Sxy = Sxy + Sxyk;
218 Syy = Syy + Syyk;
219 % don't need to be divided by k because only rations are used here
220 end
221
222 end
223
224 % compute cpsd
225 function [Sxx, Svxx] = cpsdPeriodogram(x, y, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder)
226
227 Mnxx = 0;
228 Mn2xx = 0;
229 nfft = segmentEnds(1);
230 for ii = 1:nSegments
231 if detrendOrder < 0
232 xseg = x.y(segmentStarts(ii):segmentEnds(ii));
233 yseg = y.y(segmentStarts(ii):segmentEnds(ii));
234 else
235 [xseg,coeffs] = ltpda_polyreg(x.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
236 [yseg,coeffs] = ltpda_polyreg(y.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
237 end
238
239 % Compute periodogram
240 Sxxk = wosa_periodogram(xseg, yseg, winVals, nfft);
241
242 % Welford's algorithm to update mean and variance
243 Qxx = Sxxk - Mnxx;
244 Mnxx = Mnxx +Qxx/ii;
245 Mn2xx = Mn2xx + abs(Qxx.*conj(Sxxk - Mnxx));
246 end
247 Sxx = Mnxx;
248 if nSegments ==1
249 Svxx = [];
250 else
251 Svxx = Mn2xx/(nSegments-1)/nSegments;
252 end
253
254
255 end
256
257 % compute psd
258 function [Sxx, Svxx] = psdPeriodogram(x, winVals, nSegments, segmentStarts, segmentEnds, detrendOrder)
259
260 Mnxx = 0;
261 Mn2xx = 0;
262 nfft = segmentEnds(1);
263 % Loop over the segments
264 for ii = 1:nSegments
265 % Detrend if desired
266 if detrendOrder < 0
267 seg = x.y(segmentStarts(ii):segmentEnds(ii));
268 else
269 [seg,coeffs] = ltpda_polyreg(x.y(segmentStarts(ii):segmentEnds(ii)), detrendOrder);
270 end
271 % Compute periodogram
272 Sxxk = wosa_periodogram(seg, [], winVals,nfft);
273 % Welford's algorithm for updating mean and variance
274 if ii == 1
275 Mnxx = Sxxk;
276 else
277 Qxx = Sxxk - Mnxx;
278 Mnxx = Mnxx + Qxx/ii;
279 Mn2xx = Mn2xx + Qxx.*(Sxxk - Mnxx);
280 end
281 end
282 Sxx = Mnxx;
283 if nSegments == 1
284 Svxx = [];
285 else
286 Svxx = Mn2xx/(nSegments-1)/nSegments;
287 end
288
289 end
290
291 % Scaled periodogram of one or two input signals
292 function Sxx = wosa_periodogram(x, y, win, nfft)
293
294 % window data
295 xwin = x.*win;
296 isCross = false;
297 if ~isempty(y)
298 ywin = y.*win;
299 isCross = true;
300 end
301
302 % take fft
303 X = fft(xwin, nfft);
304 if isCross
305 Y = fft(ywin, nfft);
306 end
307
308 % Compute scale factor to compensate for the window power
309 K = win'*win;
310
311 % Compute scaled power
312 Sxx = X.*conj(X)/K;
313 if isCross,
314 Sxx = X.*conj(Y)/K;
315 end
316
317 end