--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/utils/test_ndeigcsd.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,660 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% TEST ndeigcsd
+% 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% HISTORY:     23-04-2009 L Ferraioli
+%                 Creation
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+%% 2 Dim Test
+
+%%% Define startinf TFs
+
+fs = 10;
+f = [logspace(-6,log10(fs/2),2000)]';
+
+% Model Stefano TF11
+dRes11 = [2.44554138162509e-011 - 1.79482547894083e-011i;
+2.44554138162509e-011 + 1.79482547894083e-011i;
+2.66402334803101e-009 +  1.1025122049153e-009i;
+2.66402334803101e-009 -  1.1025122049153e-009i;
+-7.3560293387644e-009;
+-1.82811618589835e-009 - 1.21803627800855e-009i;
+-1.82811618589835e-009 + 1.21803627800855e-009i;
+1.16258677367555e-009;
+1.65216557639319e-016;                         
+-1.78092396888606e-016;
+-2.80420398962379e-017;
+9.21305973049041e-013 - 8.24686706827269e-014i;
+9.21305973049041e-013 + 8.24686706827269e-014i;
+5.10730060739905e-010 - 3.76571756625722e-011i;
+5.10730060739905e-010 + 3.76571756625722e-011i;
+3.45893698149735e-009;
+3.98139182134446e-014 - 8.25503935419059e-014i;
+3.98139182134446e-014 + 8.25503935419059e-014i;
+-1.40595719147164e-011];
+
+dPoles11 = [0.843464045655194 - 0.0959986292915475i;
+0.843464045655194 + 0.0959986292915475i;
+0.953187595424927 - 0.0190043625473383i;
+0.953187595424927 + 0.0190043625473383i;
+0.967176277937188;
+0.995012027005247 - 0.00268322602801729i;
+0.995012027005247 + 0.00268322602801729i;
+0.996564761885673;
+0.999999366165445;
+0.999981722418555;
+0.999921882627659;
+0.999624431675213 - 0.000813407848742761i;
+0.999624431675213 + 0.000813407848742761i;
+0.997312006278751 - 0.00265611346834941i;
+0.997312006278751 + 0.00265611346834941i;
+0.990516544257531;
+0.477796923118318 - 0.311064085401834i;
+0.477796923118318 + 0.311064085401834i;
+0];
+
+dDTerms11 = 0;
+
+% Model Stefano TF12
+dRes12 = [1.44258422208796e-017 + 7.07359428613009e-019i;
+1.44258422208796e-017 - 7.07359428613009e-019i;
+-3.4918408053655e-021 - 1.05662874569329e-021i;
+-3.4918408053655e-021 + 1.05662874569329e-021i;
+-7.61773292876976e-021;
+4.84357724603939e-020 + 2.38824204294595e-019i;
+4.84357724603939e-020 - 2.38824204294595e-019i;
+-4.07088520945753e-020 - 2.31474543846105e-019i;
+-4.07088520945753e-020 + 2.31474543846105e-019i;
+8.73316588658882e-023;
+-5.21840635377469e-020;
+1.8461911504859e-023;                         
+5.20105247464461e-020;
+-4.68960092394415e-022;
+-1.44261407664171e-017 +  6.8922564526833e-019i;
+-1.44261407664171e-017 -  6.8922564526833e-019i;
+3.13688133935426e-022];
+
+dPoles12 = [0.477546340377332 -     0.310830571032376i;
+          0.477546340377332 +     0.310830571032376i;
+           0.99790715414307 -    0.0028490561287024i;
+           0.99790715414307 +    0.0028490561287024i;
+          0.998014205354671                         ;
+          0.999585354543332 -  0.000780408757425194i;
+          0.999585354543332 +  0.000780408757425194i;
+           0.99966003029931 -  0.000830944038363768i;
+           0.99966003029931 +  0.000830944038363768i;
+          0.999962770401331                         ;
+          0.999981881865521                         ;
+          0.999999365763457                         ;
+          0.999981706320212                         ;
+           0.99992421574188                         ;
+          0.477898460791003 +     0.311001926610074i;
+          0.477898460791003 -     0.311001926610074i;
+                          0];
+dDTerms12 = 0;
+
+% Model Stefano Tf21
+dRes21 = [-1.80035241582968e-016 + 1.99543917791863e-015i;
+     -1.80035241582968e-016 - 1.99543917791863e-015i;
+     -1.85590889333759e-013 - 1.23844418827409e-014i;
+     -1.85590889333759e-013 + 1.23844418827409e-014i;
+      5.03656596876842e-013                         ;
+     -2.62470963499904e-013 + 2.30024232938878e-012i;
+     -2.62470963499904e-013 - 2.30024232938878e-012i;
+     -9.83780507870955e-018                         ;
+      3.40426735130194e-021                         ;
+      9.78322351492755e-018                         ;
+     -1.65010934542937e-020                         ;
+      2.60918565203438e-015 +  1.0546609464659e-015i;
+      2.60918565203438e-015 -  1.0546609464659e-015i;
+      3.18105585405455e-014 + 2.48839990780042e-013i;
+      3.18105585405455e-014 - 2.48839990780042e-013i;
+      3.23021641947666e-013                         ;
+      4.81265000078114e-016 - 3.18269170053848e-017i;
+      4.81265000078114e-016 + 3.18269170053848e-017i;
+      5.16260024128201e-018];
+
+dPoles21 = [0.872004077421604 -     0.110344282822693i;
+          0.872004077421604 +     0.110344282822693i;
+          0.956884129232757 -    0.0225532091775074i;
+          0.956884129232757 +    0.0225532091775074i;
+          0.966514825697177                         ;
+          0.995140550419744 -  0.000755127639524413i;
+          0.995140550419744 +  0.000755127639524413i;
+          0.999981802194393                         ;
+           0.99999936576546                         ;
+          0.999981722418555                         ;
+          0.999921882627659                         ;
+          0.999624431675213 +  0.000813407848742761i;
+          0.999624431675213 -  0.000813407848742761i;
+          0.997312006278751 +   0.00265611346834941i;
+          0.997312006278751 -   0.00265611346834941i;
+          0.990516544257531                         ;
+          0.477796923118318 +     0.311064085401834i;
+          0.477796923118318 -     0.311064085401834i;
+                          0];
+                        
+dDTerms21 = 0;
+
+% Model Stefano Tf22
+dRes22 = [1.1284521501259e-014;
+     -3.72133611555879e-014 - 2.08232683444075e-014i;
+     -3.72133611555879e-014 + 2.08232683444075e-014i;
+      9.84930639106637e-014 - 1.46640810672565e-013i;
+      9.84930639106637e-014 + 1.46640810672565e-013i;
+      2.51323684013671e-014                         ;
+     -5.64078525288305e-014                         ;
+     -1.62476406586366e-014                         ;
+     -1.08424815979566e-011 + 8.32328079357669e-012i;
+     -1.08424815979566e-011 - 8.32328079357669e-012i;
+      2.41831559776112e-011];
+    
+dPoles22 = [0.988511243978897;
+          0.997305870640646 +   0.00211760900132725i;
+          0.997305870640646 -   0.00211760900132725i;
+          0.999626453270255 +    0.0008125673525946i;
+          0.999626453270255 -    0.0008125673525946i;
+          0.999999366366222                         ;
+          0.999981706320212                         ;
+           0.99992421574188                         ;
+          0.477898460791003 -     0.311001926610074i;
+          0.477898460791003 +     0.311001926610074i;
+                          0];
+                        
+dDTerms22 = 0;
+
+% response calculation
+
+pfparams.type = 'disc';
+pfparams.freq = f;
+pfparams.fs = 10;
+
+% TF11
+pfparams.res = dRes11;
+pfparams.pol = dPoles11;
+pfparams.dterm = dDTerms11;
+pfr = utils.math.pfresp(pfparams);
+mtf11 = pfr.resp;
+
+% TF12
+pfparams.res = dRes12;
+pfparams.pol = dPoles12;
+pfparams.dterm = dDTerms12;
+pfr = utils.math.pfresp(pfparams);
+mtf12 = pfr.resp;
+
+% TF21
+pfparams.res = dRes21;
+pfparams.pol = dPoles21;
+pfparams.dterm = dDTerms21;
+pfr = utils.math.pfresp(pfparams);
+mtf21 = pfr.resp;
+
+% TF22
+pfparams.res = dRes22;
+pfparams.pol = dPoles22;
+pfparams.dterm = dDTerms22;
+pfr = utils.math.pfresp(pfparams);
+mtf22 = pfr.resp;
+
+% CSD calculation with Papoulis Method Style
+% csd11 = mtf11.*conj(mtf11)+mtf12.*conj(mtf12);
+% csd12 = mtf11.*conj(mtf21)+mtf12.*conj(mtf22);
+% csd22 = mtf22.*conj(mtf22)+mtf21.*conj(mtf21);
+% csd21 = conj(csd12);
+
+% CSD = zeros(2,2,length(f));
+% for hh = 1:length(f)
+%   CSD(:,:,hh) = [csd11(hh) csd12(hh);csd21(hh) csd22(hh)];
+% end
+
+% CSD = zeros(2,2,length(f));
+% for hh = 1:length(f)
+%   CSD(:,:,hh) = [csd11(hh) 0;0 csd22(hh)];
+% end
+
+
+%% Get CSD Model
+
+fs = 10;
+f = [logspace(-6,log10(fs/2),2000)]';
+
+% currPath = cd;
+% cd MDC1 % assuming to start from ..\ltpda\software\m-toolbox\test
+% [TF,CSD] = mdc1_tf_models(plist('f',f,'fs',fs));
+% cd(currPath)
+
+% get noise psd shapes
+% parameters names
+parnames = {'DH','S11','S1D','SD1','SDD','Dh1','Dh2','w1','w2',...
+  'TH','Th1','Th2'};
+
+% nominal params values
+parvalues = {1,1,0,-1e-4,1,1,1,-1.3e-6,-2.0e-6,0.35,0.25,0.28};
+
+% Noise model
+N = matrix(plist('built-in','mdc3_ifo2ifo_noise'));
+
+% evalueat model
+[rw,cl]=size(N.objs);
+for ii=1:rw
+  for jj=1:cl
+    N.objs(ii,jj).setParams(parnames,parvalues);
+    N.objs(ii,jj).setXvals(f);
+    sp(ii,jj)=eval(N.objs(ii,jj));
+    if ii==jj
+      sp(ii,jj)= abs(sp(ii,jj));
+    end
+  end
+end
+
+csd11 = sp(1,1).y;
+csd12 = sp(1,2).y;
+csd21 = sp(2,1).y;
+csd22 = sp(2,2).y;
+
+CSD = zeros(2,2,length(f));
+for hh = 1:length(f)
+  CSD(:,:,hh) = [csd11(hh) csd12(hh);csd21(hh) csd22(hh)];
+end
+
+%% Eigcsd 2-dim
+
+h = utils.math.ndeigcsd(CSD,'MTD','PAP');
+
+%% Test eigenshuffle
+
+[l,m,npts] = size(CSD);
+
+% Finding suppression  
+k = min(sqrt(csd11./csd22));
+if k>=1
+  suppr = floor(k);
+else
+  n=0;
+  while k<1
+    k=k*10;
+    n=n+1;
+  end
+  k = floor(k);
+  suppr = k*10^(-n);
+end
+
+supmat = [1 0;0 suppr];
+% isupmat = [1 0;0 1/suppr];
+isupmat = inv(supmat);
+
+PP = CSD;
+for ii=1:npts
+  PP(:,:,ii)=supmat*CSD(:,:,npts)*supmat;
+end
+
+[V,D] = eigenshuffle(CSD);
+h = ones(l,m,npts);
+
+for ii=1:npts
+  HH = isupmat*V(:,:,ii)*sqrt(diag(D(:,ii)));
+  h(:,:,ii) = HH;
+end
+
+%% Plot CSD
+
+% CSD calculation
+ncsd11 = squeeze(h(1,1,:).*conj(h(1,1,:))+h(1,2,:).*conj(h(1,2,:)));
+ncsd12 = squeeze(h(1,1,:).*conj(h(2,1,:))+h(1,2,:).*conj(h(2,2,:)));
+ncsd22 = squeeze(h(2,2,:).*conj(h(2,2,:))+h(2,1,:).*conj(h(2,1,:)));
+ncsd21 = conj(csd12);
+
+figure()
+loglog(f,abs(csd11),'k')
+hold on
+grid on
+loglog(f,abs(ncsd11),'r')
+
+figure()
+loglog(f,abs(csd12),'k')
+hold on
+grid on
+loglog(f,abs(ncsd12),'r')
+
+figure()
+semilogx(f,angle(csd12),'k')
+hold on
+grid on
+semilogx(f,angle(ncsd12),'r')
+
+figure()
+loglog(f,abs(csd22),'k')
+hold on
+grid on
+loglog(f,abs(ncsd22),'r')
+
+%% plot TFs
+
+figure()
+loglog(f,abs(mtf11),'k')
+hold on
+grid on
+loglog(f,abs(squeeze(h(1,1,:))),'r')
+
+figure()
+loglog(f,abs(mtf12),'k')
+hold on
+grid on
+loglog(f,abs(squeeze(h(1,2,:))),'r')
+
+figure()
+loglog(f,abs(mtf21),'k')
+hold on
+grid on
+loglog(f,abs(squeeze(h(2,1,:))),'r')
+
+figure()
+loglog(f,abs(mtf22),'k')
+hold on
+grid on
+loglog(f,abs(squeeze(h(2,2,:))),'r')
+
+%% plot TFs phase
+
+figure()
+semilogx(f,angle(mtf11),'k')
+hold on
+grid on
+semilogx(f,angle(squeeze(h(1,1,:))),'r')
+
+figure()
+semilogx(f,angle(mtf12),'k')
+hold on
+grid on
+semilogx(f,angle(squeeze(h(1,2,:))),'r')
+
+figure()
+semilogx(f,angle(mtf21),'k')
+hold on
+grid on
+semilogx(f,angle(squeeze(h(2,1,:))),'r')
+
+figure()
+semilogx(f,angle(mtf22),'k')
+hold on
+grid on
+semilogx(f,angle(squeeze(h(2,2,:))),'r')
+
+
+
+%% 3 Dim Test
+
+
+% sampling frequency
+fs = 1; % Hz
+% frequency vector
+f = logspace(-6,log10(0.5),300);
+f = f.';
+
+%%% Set the 3 channel model
+
+% d = [1 -1.1 0.5 0.02];
+% 
+% h11 = miir(1e-1.*[1 -0.5],d,fs);
+% h12 = miir(1e-4.*[0 -0.5],d,fs);
+% h13 = miir(1e-2.*[1 0.2],d,fs);
+% 
+% h21 = miir(1e-3.*[0 0.4],d,fs);
+% h22 = miir(1e-5.*[1 -0.6],d,fs);
+% h23 = miir(1e-6.*[1 0.3],d,fs);
+% 
+% h31 = miir(1e-4.*[0 0.1 -0.2],d,fs);
+% h32 = miir(1e-7.*[0 -0.6],d,fs);
+% h33 = miir(1e-5.*[1 0.05 0.3],d,fs);
+
+
+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)]));
+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)]));
+h13 = pzmodel(plist('GAIN', [0.00001], 'POLES', [pz(5.0000000000000004e-006,NaN) pz(0.02,NaN)], 'ZEROS', pz(0.00040000000000000002,NaN)));
+
+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)]));
+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)]));
+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)]));
+
+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)]));
+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)]));
+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)]));
+
+%%% TFs resps
+
+% filters response
+plresp = plist('f',f);
+
+rh11 = resp(h11,plresp);
+rh12 = resp(h12,plresp);
+rh13 = resp(h13,plresp);
+
+rh21 = resp(h21,plresp);
+rh22 = resp(h22,plresp);
+rh23 = resp(h23,plresp);
+
+rh31 = resp(h31,plresp);
+rh32 = resp(h32,plresp);
+rh33 = resp(h33,plresp);
+
+% iplot(rh11,rh12,rh13)
+% iplot(rh21,rh22,rh23)
+% iplot(rh31,rh32,rh33)
+
+% iplot(rh11,rh22,rh33)
+
+%%% Spectra calculation with Kay style method
+
+% Note that the definition of cross-spectral matrix follows that of s M
+% Kay, Modern Spectral Estimation
+fs = 1; % Hz
+% csd11 
+G11 = rh11.y.*conj(rh11.y) + rh12.y.*conj(rh12.y) + rh13.y.*conj(rh13.y);
+G11 = ao(plist('xvals', f, 'yvals', G11, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+G11.setName;
+
+% csd12
+G12 = conj(rh11.y).*rh21.y + conj(rh12.y).*rh22.y + conj(rh13.y).*rh23.y;
+G12 = ao(plist('xvals', f, 'yvals', G12, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+G12.setName;
+
+% csd13
+G13 = conj(rh11.y).*rh31.y + conj(rh12.y).*rh32.y + conj(rh13.y).*rh33.y;
+G13 = ao(plist('xvals', f, 'yvals', G13, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+G13.setName;
+
+% csd21
+G21 = conj(G12);
+G21.setName;
+
+% csd22
+G22 = rh21.y.*conj(rh21.y) + rh22.y.*conj(rh22.y) + rh23.y.*conj(rh23.y);
+G22 = ao(plist('xvals', f, 'yvals', G22, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+G22.setName;
+
+% csd23
+G23 = conj(rh21.y).*rh31.y + conj(rh22.y).*rh32.y + conj(rh23.y).*rh33.y;
+G23 = ao(plist('xvals', f, 'yvals', G23, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+G23.setName;
+
+% csd31
+G31 = conj(G13);
+G31.setName;
+
+% csd32
+G32 = conj(G23);
+G32.setName;
+
+% csd33
+G33 = conj(rh31.y).*rh31.y + conj(rh32.y).*rh32.y + conj(rh33.y).*rh33.y;
+G33 = ao(plist('xvals', f, 'yvals', G33, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+G33.setName;
+
+% iplot(G11,G12,G13,G22,G23,G33)
+
+% iplot(G11,G22,G33)
+
+%%% Build CSD matrix
+
+CSD = zeros(3,3,length(f));
+for hh = 1:length(f)
+  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)];
+end
+
+% for hh = 1:length(f)
+%   CSD(:,:,hh) = [G11.y(hh) 0 0;0 G22.y(hh) 0;0 0 G33.y(hh)];
+% end
+
+%%% 3-dim Eigcsd with kay method stule
+h = utils.math.ndeigcsd(CSD,'MTD','KAY');
+
+% %%% Eigenshuffle
+% 
+% [l,m,kk] = size(CSD);
+% suppr = ones(l,l);
+% for ii = 2:l
+%   k = ones(l,1);
+%   
+%   for jj = ii-1:-1:1
+%     k(jj) = min(sqrt(CSD(jj,jj,:)./CSD(ii,ii,:)));
+%     if k(jj)>=1
+%       suppr(jj,ii) = floor(k(jj));
+%     else
+%       n=0;
+%       while k(jj)<1
+%         k(jj)=k(jj)*10;
+%         n=n+1;
+%       end
+%       k(jj) = floor(k(jj));
+%       suppr(jj,ii) = k(jj)*10^(-n);
+%       %       suppr(ii) = suppr(ii)*suppr(ii-1);
+%     end
+%   end
+%   %     csuppr(ii) = prod(suppr(:,ii));
+% end
+% csuppr = prod(suppr,2);
+% supmat = diag(csuppr);
+% % ssup = rot90(rot90(supmat));
+% ssup = supmat;
+% % ssup = eye(3);
+% % ssup = supmat*supmat.';
+% issup = inv(ssup);
+% 
+% for jj = 1:kk
+%   nCSD(:,:,jj) = ssup*CSD(:,:,jj)*ssup;
+% end
+% 
+% [Vseq,Dseq] = eigenshuffle(nCSD);
+% 
+% % get h = V*sqrt(D);
+% [nn,mm,kk] = size(Vseq);
+% h = zeros(nn,mm,kk);
+% for dd = 1:kk
+% %   h(:,:,dd) = rot90(issup*Vseq(:,:,dd)*sqrt(diag(Dseq(:,dd))));
+%   h(:,:,dd) = issup*Vseq(:,:,dd)*sqrt(diag(Dseq(:,dd)));
+% end
+
+%% Build TF AOs
+
+nh11 = ao(plist('xvals', f, 'yvals', squeeze(h(1,1,:)), 'fs', fs, 'dtype', 'fsdata'));
+nh11.setName;
+
+nh12 = ao(plist('xvals', f, 'yvals', squeeze(h(1,2,:)), 'fs', fs, 'dtype', 'fsdata'));
+nh12.setName;
+
+nh13 = ao(plist('xvals', f, 'yvals', squeeze(h(1,3,:)), 'fs', fs, 'dtype', 'fsdata'));
+nh13.setName;
+
+nh21 = ao(plist('xvals', f, 'yvals', squeeze(h(2,1,:)), 'fs', fs, 'dtype', 'fsdata'));
+nh21.setName;
+
+nh22 = ao(plist('xvals', f, 'yvals', squeeze(h(2,2,:)), 'fs', fs, 'dtype', 'fsdata'));
+nh22.setName;
+
+nh23 = ao(plist('xvals', f, 'yvals', squeeze(h(2,3,:)), 'fs', fs, 'dtype', 'fsdata'));
+nh23.setName;
+
+nh31 = ao(plist('xvals', f, 'yvals', squeeze(h(3,1,:)), 'fs', fs, 'dtype', 'fsdata'));
+nh31.setName;
+
+nh32 = ao(plist('xvals', f, 'yvals', squeeze(h(3,2,:)), 'fs', fs, 'dtype', 'fsdata'));
+nh32.setName;
+
+nh33 = ao(plist('xvals', f, 'yvals', squeeze(h(3,3,:)), 'fs', fs, 'dtype', 'fsdata'));
+nh33.setName;
+
+iplot(rh11,nh11)
+iplot(rh12,nh12)
+iplot(rh13,nh13)
+iplot(rh21,nh21)
+iplot(rh22,nh22)
+iplot(rh23,nh23)
+iplot(rh31,nh31)
+iplot(rh32,nh32)
+iplot(rh33,nh33)
+
+%% Build CSD AOs
+
+% Note that the definition of cross-spectral matrix follows that of s M
+% Kay, Modern Spectral Estimation
+
+% csd11 
+H11 = nh11.y.*conj(nh11.y) + nh12.y.*conj(nh12.y) + nh13.y.*conj(nh13.y);
+H11 = ao(plist('xvals', f, 'yvals', H11, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+H11.setName;
+
+% csd12
+H12 = conj(nh11.y).*nh21.y + conj(nh12.y).*nh22.y + conj(nh13.y).*nh23.y;
+H12 = ao(plist('xvals', f, 'yvals', H12, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+H12.setName;
+
+% csd13
+H13 = conj(nh11.y).*nh31.y + conj(nh12.y).*nh32.y + conj(nh13.y).*nh33.y;
+H13 = ao(plist('xvals', f, 'yvals', H13, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+H13.setName;
+
+% csd21
+H21 = conj(H12);
+H21.setName;
+
+% csd22
+H22 = nh21.y.*conj(nh21.y) + nh22.y.*conj(nh22.y) + nh23.y.*conj(nh23.y);
+H22 = ao(plist('xvals', f, 'yvals', H22, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+H22.setName;
+
+% csd23
+H23 = conj(nh21.y).*nh31.y + conj(nh22.y).*nh32.y + conj(nh23.y).*nh33.y;
+H23 = ao(plist('xvals', f, 'yvals', H23, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+H23.setName;
+
+% csd31
+H31 = conj(H13);
+H31.setName;
+
+% csd32
+H32 = conj(H23);
+H32.setName;
+
+% csd33
+H33 = conj(nh31.y).*nh31.y + conj(nh32.y).*nh32.y + conj(nh33.y).*nh33.y;
+H33 = ao(plist('xvals', f, 'yvals', H33, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
+H33.setName;
+
+iplot(G11,H11)
+iplot(G12,H12)
+iplot(G13,H13)
+iplot(G22,H22)
+iplot(G23,H23)
+iplot(G33,H33)
+
+
+
+
+
+
+
+