Mercurial > hg > ltpda

comparison m-toolbox/test/utils/test_ndeigcsd.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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % TEST ndeigcsd
3 %
4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 % HISTORY: 23-04-2009 L Ferraioli
6 % Creation
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8
9
10
11 %% 2 Dim Test
12
13 %%% Define startinf TFs
14
15 fs = 10;
16 f = [logspace(-6,log10(fs/2),2000)]';
17
18 % Model Stefano TF11
19 dRes11 = [2.44554138162509e-011 - 1.79482547894083e-011i;
20 2.44554138162509e-011 + 1.79482547894083e-011i;
21 2.66402334803101e-009 + 1.1025122049153e-009i;
22 2.66402334803101e-009 - 1.1025122049153e-009i;
23 -7.3560293387644e-009;
24 -1.82811618589835e-009 - 1.21803627800855e-009i;
25 -1.82811618589835e-009 + 1.21803627800855e-009i;
26 1.16258677367555e-009;
27 1.65216557639319e-016;
28 -1.78092396888606e-016;
29 -2.80420398962379e-017;
30 9.21305973049041e-013 - 8.24686706827269e-014i;
31 9.21305973049041e-013 + 8.24686706827269e-014i;
32 5.10730060739905e-010 - 3.76571756625722e-011i;
33 5.10730060739905e-010 + 3.76571756625722e-011i;
34 3.45893698149735e-009;
35 3.98139182134446e-014 - 8.25503935419059e-014i;
36 3.98139182134446e-014 + 8.25503935419059e-014i;
37 -1.40595719147164e-011];
38
39 dPoles11 = [0.843464045655194 - 0.0959986292915475i;
40 0.843464045655194 + 0.0959986292915475i;
41 0.953187595424927 - 0.0190043625473383i;
42 0.953187595424927 + 0.0190043625473383i;
43 0.967176277937188;
44 0.995012027005247 - 0.00268322602801729i;
45 0.995012027005247 + 0.00268322602801729i;
46 0.996564761885673;
47 0.999999366165445;
48 0.999981722418555;
49 0.999921882627659;
50 0.999624431675213 - 0.000813407848742761i;
51 0.999624431675213 + 0.000813407848742761i;
52 0.997312006278751 - 0.00265611346834941i;
53 0.997312006278751 + 0.00265611346834941i;
54 0.990516544257531;
55 0.477796923118318 - 0.311064085401834i;
56 0.477796923118318 + 0.311064085401834i;
57 0];
58
59 dDTerms11 = 0;
60
61 % Model Stefano TF12
62 dRes12 = [1.44258422208796e-017 + 7.07359428613009e-019i;
63 1.44258422208796e-017 - 7.07359428613009e-019i;
64 -3.4918408053655e-021 - 1.05662874569329e-021i;
65 -3.4918408053655e-021 + 1.05662874569329e-021i;
66 -7.61773292876976e-021;
67 4.84357724603939e-020 + 2.38824204294595e-019i;
68 4.84357724603939e-020 - 2.38824204294595e-019i;
69 -4.07088520945753e-020 - 2.31474543846105e-019i;
70 -4.07088520945753e-020 + 2.31474543846105e-019i;
71 8.73316588658882e-023;
72 -5.21840635377469e-020;
73 1.8461911504859e-023;
74 5.20105247464461e-020;
75 -4.68960092394415e-022;
76 -1.44261407664171e-017 + 6.8922564526833e-019i;
77 -1.44261407664171e-017 - 6.8922564526833e-019i;
78 3.13688133935426e-022];
79
80 dPoles12 = [0.477546340377332 - 0.310830571032376i;
81 0.477546340377332 + 0.310830571032376i;
82 0.99790715414307 - 0.0028490561287024i;
83 0.99790715414307 + 0.0028490561287024i;
84 0.998014205354671 ;
85 0.999585354543332 - 0.000780408757425194i;
86 0.999585354543332 + 0.000780408757425194i;
87 0.99966003029931 - 0.000830944038363768i;
88 0.99966003029931 + 0.000830944038363768i;
89 0.999962770401331 ;
90 0.999981881865521 ;
91 0.999999365763457 ;
92 0.999981706320212 ;
93 0.99992421574188 ;
94 0.477898460791003 + 0.311001926610074i;
95 0.477898460791003 - 0.311001926610074i;
96 0];
97 dDTerms12 = 0;
98
99 % Model Stefano Tf21
100 dRes21 = [-1.80035241582968e-016 + 1.99543917791863e-015i;
101 -1.80035241582968e-016 - 1.99543917791863e-015i;
102 -1.85590889333759e-013 - 1.23844418827409e-014i;
103 -1.85590889333759e-013 + 1.23844418827409e-014i;
104 5.03656596876842e-013 ;
105 -2.62470963499904e-013 + 2.30024232938878e-012i;
106 -2.62470963499904e-013 - 2.30024232938878e-012i;
107 -9.83780507870955e-018 ;
108 3.40426735130194e-021 ;
109 9.78322351492755e-018 ;
110 -1.65010934542937e-020 ;
111 2.60918565203438e-015 + 1.0546609464659e-015i;
112 2.60918565203438e-015 - 1.0546609464659e-015i;
113 3.18105585405455e-014 + 2.48839990780042e-013i;
114 3.18105585405455e-014 - 2.48839990780042e-013i;
115 3.23021641947666e-013 ;
116 4.81265000078114e-016 - 3.18269170053848e-017i;
117 4.81265000078114e-016 + 3.18269170053848e-017i;
118 5.16260024128201e-018];
119
120 dPoles21 = [0.872004077421604 - 0.110344282822693i;
121 0.872004077421604 + 0.110344282822693i;
122 0.956884129232757 - 0.0225532091775074i;
123 0.956884129232757 + 0.0225532091775074i;
124 0.966514825697177 ;
125 0.995140550419744 - 0.000755127639524413i;
126 0.995140550419744 + 0.000755127639524413i;
127 0.999981802194393 ;
128 0.99999936576546 ;
129 0.999981722418555 ;
130 0.999921882627659 ;
131 0.999624431675213 + 0.000813407848742761i;
132 0.999624431675213 - 0.000813407848742761i;
133 0.997312006278751 + 0.00265611346834941i;
134 0.997312006278751 - 0.00265611346834941i;
135 0.990516544257531 ;
136 0.477796923118318 + 0.311064085401834i;
137 0.477796923118318 - 0.311064085401834i;
138 0];
139
140 dDTerms21 = 0;
141
142 % Model Stefano Tf22
143 dRes22 = [1.1284521501259e-014;
144 -3.72133611555879e-014 - 2.08232683444075e-014i;
145 -3.72133611555879e-014 + 2.08232683444075e-014i;
146 9.84930639106637e-014 - 1.46640810672565e-013i;
147 9.84930639106637e-014 + 1.46640810672565e-013i;
148 2.51323684013671e-014 ;
149 -5.64078525288305e-014 ;
150 -1.62476406586366e-014 ;
151 -1.08424815979566e-011 + 8.32328079357669e-012i;
152 -1.08424815979566e-011 - 8.32328079357669e-012i;
153 2.41831559776112e-011];
154
155 dPoles22 = [0.988511243978897;
156 0.997305870640646 + 0.00211760900132725i;
157 0.997305870640646 - 0.00211760900132725i;
158 0.999626453270255 + 0.0008125673525946i;
159 0.999626453270255 - 0.0008125673525946i;
160 0.999999366366222 ;
161 0.999981706320212 ;
162 0.99992421574188 ;
163 0.477898460791003 - 0.311001926610074i;
164 0.477898460791003 + 0.311001926610074i;
165 0];
166
167 dDTerms22 = 0;
168
169 % response calculation
170
171 pfparams.type = 'disc';
172 pfparams.freq = f;
173 pfparams.fs = 10;
174
175 % TF11
176 pfparams.res = dRes11;
177 pfparams.pol = dPoles11;
178 pfparams.dterm = dDTerms11;
179 pfr = utils.math.pfresp(pfparams);
180 mtf11 = pfr.resp;
181
182 % TF12
183 pfparams.res = dRes12;
184 pfparams.pol = dPoles12;
185 pfparams.dterm = dDTerms12;
186 pfr = utils.math.pfresp(pfparams);
187 mtf12 = pfr.resp;
188
189 % TF21
190 pfparams.res = dRes21;
191 pfparams.pol = dPoles21;
192 pfparams.dterm = dDTerms21;
193 pfr = utils.math.pfresp(pfparams);
194 mtf21 = pfr.resp;
195
196 % TF22
197 pfparams.res = dRes22;
198 pfparams.pol = dPoles22;
199 pfparams.dterm = dDTerms22;
200 pfr = utils.math.pfresp(pfparams);
201 mtf22 = pfr.resp;
202
203 % CSD calculation with Papoulis Method Style
204 % csd11 = mtf11.*conj(mtf11)+mtf12.*conj(mtf12);
205 % csd12 = mtf11.*conj(mtf21)+mtf12.*conj(mtf22);
206 % csd22 = mtf22.*conj(mtf22)+mtf21.*conj(mtf21);
207 % csd21 = conj(csd12);
208
209 % CSD = zeros(2,2,length(f));
210 % for hh = 1:length(f)
211 % CSD(:,:,hh) = [csd11(hh) csd12(hh);csd21(hh) csd22(hh)];
212 % end
213
214 % CSD = zeros(2,2,length(f));
215 % for hh = 1:length(f)
216 % CSD(:,:,hh) = [csd11(hh) 0;0 csd22(hh)];
217 % end
218
219
220 %% Get CSD Model
221
222 fs = 10;
223 f = [logspace(-6,log10(fs/2),2000)]';
224
225 % currPath = cd;
226 % cd MDC1 % assuming to start from ..\ltpda\software\m-toolbox\test
227 % [TF,CSD] = mdc1_tf_models(plist('f',f,'fs',fs));
228 % cd(currPath)
229
230 % get noise psd shapes
231 % parameters names
232 parnames = {'DH','S11','S1D','SD1','SDD','Dh1','Dh2','w1','w2',...
233 'TH','Th1','Th2'};
234
235 % nominal params values
236 parvalues = {1,1,0,-1e-4,1,1,1,-1.3e-6,-2.0e-6,0.35,0.25,0.28};
237
238 % Noise model
239 N = matrix(plist('built-in','mdc3_ifo2ifo_noise'));
240
241 % evalueat model
242 [rw,cl]=size(N.objs);
243 for ii=1:rw
244 for jj=1:cl
245 N.objs(ii,jj).setParams(parnames,parvalues);
246 N.objs(ii,jj).setXvals(f);
247 sp(ii,jj)=eval(N.objs(ii,jj));
248 if ii==jj
249 sp(ii,jj)= abs(sp(ii,jj));
250 end
251 end
252 end
253
254 csd11 = sp(1,1).y;
255 csd12 = sp(1,2).y;
256 csd21 = sp(2,1).y;
257 csd22 = sp(2,2).y;
258
259 CSD = zeros(2,2,length(f));
260 for hh = 1:length(f)
261 CSD(:,:,hh) = [csd11(hh) csd12(hh);csd21(hh) csd22(hh)];
262 end
263
264 %% Eigcsd 2-dim
265
266 h = utils.math.ndeigcsd(CSD,'MTD','PAP');
267
268 %% Test eigenshuffle
269
270 [l,m,npts] = size(CSD);
271
272 % Finding suppression
273 k = min(sqrt(csd11./csd22));
274 if k>=1
275 suppr = floor(k);
276 else
277 n=0;
278 while k<1
279 k=k*10;
280 n=n+1;
281 end
282 k = floor(k);
283 suppr = k*10^(-n);
284 end
285
286 supmat = [1 0;0 suppr];
287 % isupmat = [1 0;0 1/suppr];
288 isupmat = inv(supmat);
289
290 PP = CSD;
291 for ii=1:npts
292 PP(:,:,ii)=supmat*CSD(:,:,npts)*supmat;
293 end
294
295 [V,D] = eigenshuffle(CSD);
296 h = ones(l,m,npts);
297
298 for ii=1:npts
299 HH = isupmat*V(:,:,ii)*sqrt(diag(D(:,ii)));
300 h(:,:,ii) = HH;
301 end
302
303 %% Plot CSD
304
305 % CSD calculation
306 ncsd11 = squeeze(h(1,1,:).*conj(h(1,1,:))+h(1,2,:).*conj(h(1,2,:)));
307 ncsd12 = squeeze(h(1,1,:).*conj(h(2,1,:))+h(1,2,:).*conj(h(2,2,:)));
308 ncsd22 = squeeze(h(2,2,:).*conj(h(2,2,:))+h(2,1,:).*conj(h(2,1,:)));
309 ncsd21 = conj(csd12);
310
311 figure()
312 loglog(f,abs(csd11),'k')
313 hold on
314 grid on
315 loglog(f,abs(ncsd11),'r')
316
317 figure()
318 loglog(f,abs(csd12),'k')
319 hold on
320 grid on
321 loglog(f,abs(ncsd12),'r')
322
323 figure()
324 semilogx(f,angle(csd12),'k')
325 hold on
326 grid on
327 semilogx(f,angle(ncsd12),'r')
328
329 figure()
330 loglog(f,abs(csd22),'k')
331 hold on
332 grid on
333 loglog(f,abs(ncsd22),'r')
334
335 %% plot TFs
336
337 figure()
338 loglog(f,abs(mtf11),'k')
339 hold on
340 grid on
341 loglog(f,abs(squeeze(h(1,1,:))),'r')
342
343 figure()
344 loglog(f,abs(mtf12),'k')
345 hold on
346 grid on
347 loglog(f,abs(squeeze(h(1,2,:))),'r')
348
349 figure()
350 loglog(f,abs(mtf21),'k')
351 hold on
352 grid on
353 loglog(f,abs(squeeze(h(2,1,:))),'r')
354
355 figure()
356 loglog(f,abs(mtf22),'k')
357 hold on
358 grid on
359 loglog(f,abs(squeeze(h(2,2,:))),'r')
360
361 %% plot TFs phase
362
363 figure()
364 semilogx(f,angle(mtf11),'k')
365 hold on
366 grid on
367 semilogx(f,angle(squeeze(h(1,1,:))),'r')
368
369 figure()
370 semilogx(f,angle(mtf12),'k')
371 hold on
372 grid on
373 semilogx(f,angle(squeeze(h(1,2,:))),'r')
374
375 figure()
376 semilogx(f,angle(mtf21),'k')
377 hold on
378 grid on
379 semilogx(f,angle(squeeze(h(2,1,:))),'r')
380
381 figure()
382 semilogx(f,angle(mtf22),'k')
383 hold on
384 grid on
385 semilogx(f,angle(squeeze(h(2,2,:))),'r')
386
387
388
389 %% 3 Dim Test
390
391
392 % sampling frequency
393 fs = 1; % Hz
394 % frequency vector
395 f = logspace(-6,log10(0.5),300);
396 f = f.';
397
398 %%% Set the 3 channel model
399
400 % d = [1 -1.1 0.5 0.02];
401 %
402 % h11 = miir(1e-1.*[1 -0.5],d,fs);
403 % h12 = miir(1e-4.*[0 -0.5],d,fs);
404 % h13 = miir(1e-2.*[1 0.2],d,fs);
405 %
406 % h21 = miir(1e-3.*[0 0.4],d,fs);
407 % h22 = miir(1e-5.*[1 -0.6],d,fs);
408 % h23 = miir(1e-6.*[1 0.3],d,fs);
409 %
410 % h31 = miir(1e-4.*[0 0.1 -0.2],d,fs);
411 % h32 = miir(1e-7.*[0 -0.6],d,fs);
412 % h33 = miir(1e-5.*[1 0.05 0.3],d,fs);
413
414
415 h11 = pzmodel(plist('GAIN', [0.01], 'POLES', [pz(9.9999999999999995e-007,NaN) pz(0.01,NaN) pz(0.02,NaN)], 'ZEROS', [pz(0.0001,NaN) pz(0.002,NaN)]));
416 h12 = pzmodel(plist('GAIN', [0.0001], 'POLES', [pz(0.10000000000000001,NaN) pz(0.0030000000000000001,NaN) pz(0.0050000000000000001,NaN) pz(1.0000000000000001e-005,NaN)], 'ZEROS', [pz(0.00050000000000000001,NaN) pz(0.001,NaN)]));
417 h13 = pzmodel(plist('GAIN', [0.00001], 'POLES', [pz(5.0000000000000004e-006,NaN) pz(0.02,NaN)], 'ZEROS', pz(0.00040000000000000002,NaN)));
418
419 h21 = pzmodel(plist('GAIN', [0.00001], 'POLES', [pz(5.0000000000000004e-006,NaN) pz(0.10000000000000001,NaN) pz(0.050000000000000003,NaN)], 'ZEROS', [pz(0.002,NaN) pz(0.00020000000000000001,NaN)]));
420 h22 = pzmodel(plist('GAIN', [0.001], 'POLES', [pz(0.01,NaN) pz(9.9999999999999995e-008,NaN) pz(0.002,NaN) pz(0.001,NaN)], 'ZEROS', [pz(1.0000000000000001e-005,NaN) pz(2.0000000000000002e-005,NaN) pz(5.0000000000000002e-005,NaN) pz(0.20000000000000001,NaN)]));
421 h23 = pzmodel(plist('GAIN', [0.0001], 'POLES', [pz(0.01,NaN) pz(0.029999999999999999,NaN) pz(0.10000000000000001,NaN) pz(5.0000000000000002e-005,NaN)], 'ZEROS', [pz(0.0001,NaN) pz(0.0050000000000000001,NaN)]));
422
423 h31 = pzmodel(plist('GAIN', [0.001], 'POLES', [pz(0.01,NaN) pz(0.029999999999999999,NaN) pz(0.10000000000000001,NaN) pz(0.00050000000000000001,NaN)], 'ZEROS', [pz(0.0001,NaN) pz(0.0050000000000000001,NaN)]));
424 h32 = pzmodel(plist('GAIN', [0.00001], 'POLES', [pz(0.01,NaN) pz(0.029999999999999999,NaN) pz(0.10000000000000001,NaN) pz(0.00050000000000000001,NaN) pz(0.00029999999999999997,NaN)], 'ZEROS', [pz(0.002,NaN) pz(0.059999999999999998,NaN)]));
425 h33 = pzmodel(plist('GAIN', [0.05], 'POLES', [pz(0.0001,NaN) pz(9.9999999999999995e-008,NaN) pz(0.002,NaN)], 'ZEROS', [pz(3.0000000000000001e-005,NaN) pz(5.0000000000000002e-005,NaN) pz(0.10000000000000001,NaN)]));
426
427 %%% TFs resps
428
429 % filters response
430 plresp = plist('f',f);
431
432 rh11 = resp(h11,plresp);
433 rh12 = resp(h12,plresp);
434 rh13 = resp(h13,plresp);
435
436 rh21 = resp(h21,plresp);
437 rh22 = resp(h22,plresp);
438 rh23 = resp(h23,plresp);
439
440 rh31 = resp(h31,plresp);
441 rh32 = resp(h32,plresp);
442 rh33 = resp(h33,plresp);
443
444 % iplot(rh11,rh12,rh13)
445 % iplot(rh21,rh22,rh23)
446 % iplot(rh31,rh32,rh33)
447
448 % iplot(rh11,rh22,rh33)
449
450 %%% Spectra calculation with Kay style method
451
452 % Note that the definition of cross-spectral matrix follows that of s M
453 % Kay, Modern Spectral Estimation
454 fs = 1; % Hz
455 % csd11
456 G11 = rh11.y.*conj(rh11.y) + rh12.y.*conj(rh12.y) + rh13.y.*conj(rh13.y);
457 G11 = ao(plist('xvals', f, 'yvals', G11, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
458 G11.setName;
459
460 % csd12
461 G12 = conj(rh11.y).*rh21.y + conj(rh12.y).*rh22.y + conj(rh13.y).*rh23.y;
462 G12 = ao(plist('xvals', f, 'yvals', G12, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
463 G12.setName;
464
465 % csd13
466 G13 = conj(rh11.y).*rh31.y + conj(rh12.y).*rh32.y + conj(rh13.y).*rh33.y;
467 G13 = ao(plist('xvals', f, 'yvals', G13, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
468 G13.setName;
469
470 % csd21
471 G21 = conj(G12);
472 G21.setName;
473
474 % csd22
475 G22 = rh21.y.*conj(rh21.y) + rh22.y.*conj(rh22.y) + rh23.y.*conj(rh23.y);
476 G22 = ao(plist('xvals', f, 'yvals', G22, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
477 G22.setName;
478
479 % csd23
480 G23 = conj(rh21.y).*rh31.y + conj(rh22.y).*rh32.y + conj(rh23.y).*rh33.y;
481 G23 = ao(plist('xvals', f, 'yvals', G23, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
482 G23.setName;
483
484 % csd31
485 G31 = conj(G13);
486 G31.setName;
487
488 % csd32
489 G32 = conj(G23);
490 G32.setName;
491
492 % csd33
493 G33 = conj(rh31.y).*rh31.y + conj(rh32.y).*rh32.y + conj(rh33.y).*rh33.y;
494 G33 = ao(plist('xvals', f, 'yvals', G33, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
495 G33.setName;
496
497 % iplot(G11,G12,G13,G22,G23,G33)
498
499 % iplot(G11,G22,G33)
500
501 %%% Build CSD matrix
502
503 CSD = zeros(3,3,length(f));
504 for hh = 1:length(f)
505 CSD(:,:,hh) = [G11.y(hh) G12.y(hh) G13.y(hh);G21.y(hh) G22.y(hh) G23.y(hh);G31.y(hh) G32.y(hh) G33.y(hh)];
506 end
507
508 % for hh = 1:length(f)
509 % CSD(:,:,hh) = [G11.y(hh) 0 0;0 G22.y(hh) 0;0 0 G33.y(hh)];
510 % end
511
512 %%% 3-dim Eigcsd with kay method stule
513 h = utils.math.ndeigcsd(CSD,'MTD','KAY');
514
515 % %%% Eigenshuffle
516 %
517 % [l,m,kk] = size(CSD);
518 % suppr = ones(l,l);
519 % for ii = 2:l
520 % k = ones(l,1);
521 %
522 % for jj = ii-1:-1:1
523 % k(jj) = min(sqrt(CSD(jj,jj,:)./CSD(ii,ii,:)));
524 % if k(jj)>=1
525 % suppr(jj,ii) = floor(k(jj));
526 % else
527 % n=0;
528 % while k(jj)<1
529 % k(jj)=k(jj)*10;
530 % n=n+1;
531 % end
532 % k(jj) = floor(k(jj));
533 % suppr(jj,ii) = k(jj)*10^(-n);
534 % % suppr(ii) = suppr(ii)*suppr(ii-1);
535 % end
536 % end
537 % % csuppr(ii) = prod(suppr(:,ii));
538 % end
539 % csuppr = prod(suppr,2);
540 % supmat = diag(csuppr);
541 % % ssup = rot90(rot90(supmat));
542 % ssup = supmat;
543 % % ssup = eye(3);
544 % % ssup = supmat*supmat.';
545 % issup = inv(ssup);
546 %
547 % for jj = 1:kk
548 % nCSD(:,:,jj) = ssup*CSD(:,:,jj)*ssup;
549 % end
550 %
551 % [Vseq,Dseq] = eigenshuffle(nCSD);
552 %
553 % % get h = V*sqrt(D);
554 % [nn,mm,kk] = size(Vseq);
555 % h = zeros(nn,mm,kk);
556 % for dd = 1:kk
557 % % h(:,:,dd) = rot90(issup*Vseq(:,:,dd)*sqrt(diag(Dseq(:,dd))));
558 % h(:,:,dd) = issup*Vseq(:,:,dd)*sqrt(diag(Dseq(:,dd)));
559 % end
560
561 %% Build TF AOs
562
563 nh11 = ao(plist('xvals', f, 'yvals', squeeze(h(1,1,:)), 'fs', fs, 'dtype', 'fsdata'));
564 nh11.setName;
565
566 nh12 = ao(plist('xvals', f, 'yvals', squeeze(h(1,2,:)), 'fs', fs, 'dtype', 'fsdata'));
567 nh12.setName;
568
569 nh13 = ao(plist('xvals', f, 'yvals', squeeze(h(1,3,:)), 'fs', fs, 'dtype', 'fsdata'));
570 nh13.setName;
571
572 nh21 = ao(plist('xvals', f, 'yvals', squeeze(h(2,1,:)), 'fs', fs, 'dtype', 'fsdata'));
573 nh21.setName;
574
575 nh22 = ao(plist('xvals', f, 'yvals', squeeze(h(2,2,:)), 'fs', fs, 'dtype', 'fsdata'));
576 nh22.setName;
577
578 nh23 = ao(plist('xvals', f, 'yvals', squeeze(h(2,3,:)), 'fs', fs, 'dtype', 'fsdata'));
579 nh23.setName;
580
581 nh31 = ao(plist('xvals', f, 'yvals', squeeze(h(3,1,:)), 'fs', fs, 'dtype', 'fsdata'));
582 nh31.setName;
583
584 nh32 = ao(plist('xvals', f, 'yvals', squeeze(h(3,2,:)), 'fs', fs, 'dtype', 'fsdata'));
585 nh32.setName;
586
587 nh33 = ao(plist('xvals', f, 'yvals', squeeze(h(3,3,:)), 'fs', fs, 'dtype', 'fsdata'));
588 nh33.setName;
589
590 iplot(rh11,nh11)
591 iplot(rh12,nh12)
592 iplot(rh13,nh13)
593 iplot(rh21,nh21)
594 iplot(rh22,nh22)
595 iplot(rh23,nh23)
596 iplot(rh31,nh31)
597 iplot(rh32,nh32)
598 iplot(rh33,nh33)
599
600 %% Build CSD AOs
601
602 % Note that the definition of cross-spectral matrix follows that of s M
603 % Kay, Modern Spectral Estimation
604
605 % csd11
606 H11 = nh11.y.*conj(nh11.y) + nh12.y.*conj(nh12.y) + nh13.y.*conj(nh13.y);
607 H11 = ao(plist('xvals', f, 'yvals', H11, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
608 H11.setName;
609
610 % csd12
611 H12 = conj(nh11.y).*nh21.y + conj(nh12.y).*nh22.y + conj(nh13.y).*nh23.y;
612 H12 = ao(plist('xvals', f, 'yvals', H12, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
613 H12.setName;
614
615 % csd13
616 H13 = conj(nh11.y).*nh31.y + conj(nh12.y).*nh32.y + conj(nh13.y).*nh33.y;
617 H13 = ao(plist('xvals', f, 'yvals', H13, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
618 H13.setName;
619
620 % csd21
621 H21 = conj(H12);
622 H21.setName;
623
624 % csd22
625 H22 = nh21.y.*conj(nh21.y) + nh22.y.*conj(nh22.y) + nh23.y.*conj(nh23.y);
626 H22 = ao(plist('xvals', f, 'yvals', H22, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
627 H22.setName;
628
629 % csd23
630 H23 = conj(nh21.y).*nh31.y + conj(nh22.y).*nh32.y + conj(nh23.y).*nh33.y;
631 H23 = ao(plist('xvals', f, 'yvals', H23, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
632 H23.setName;
633
634 % csd31
635 H31 = conj(H13);
636 H31.setName;
637
638 % csd32
639 H32 = conj(H23);
640 H32.setName;
641
642 % csd33
643 H33 = conj(nh31.y).*nh31.y + conj(nh32.y).*nh32.y + conj(nh33.y).*nh33.y;
644 H33 = ao(plist('xvals', f, 'yvals', H33, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
645 H33.setName;
646
647 iplot(G11,H11)
648 iplot(G12,H12)
649 iplot(G13,H13)
650 iplot(G22,H22)
651 iplot(G23,H23)
652 iplot(G33,H33)
653
654
655
656
657
658
659
660