--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/test_ao_whiten2D_2.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,224 @@
+% A test script for ao/whiten2D
+% 
+% L. Ferraioli 18-11-08
+% 
+% $Id: test_ao_whiten2D_2.m,v 1.3 2009/08/26 18:27:06 luigi Exp $
+% 
+
+%%
+
+clear all
+
+%% CSD Noise models and TF
+
+[CSD,TF] = get_2D_test_obj_tf_psd();
+
+% get noise coherence
+coh = CSD(1,2)./sqrt(abs(CSD(1,1).*CSD(2,2)));
+rcoh = real(coh);
+rcoh.setName;
+icoh = imag(coh);
+icoh.setName;
+
+%% Colored Noise generation
+
+% Make white noise
+
+a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 5e5));
+a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 5e5));
+
+% Generation from white noise
+
+pl = plist(...
+    'csd11', CSD(1,1),...
+    'csd12', CSD(1,2),...
+    'csd21', CSD(2,1),...
+    'csd22', CSD(2,2),...
+    'MaxIter', 70, ...
+    'PoleType', 2, ...
+    'MinOrder', 25, ...
+    'MaxOrder', 35, ...
+    'Weights', 2, ...
+    'Plot', false,...
+    'MSEVARTOL', 1e-2,...
+    'FITTOL', 1e-4,...
+    'UseSym', 0,...
+    'Disp', false);
+
+[ac1,ac2] = noisegen2D(a1,a2, pl);
+
+ac1.setYunits('m');
+ac2.setYunits('m');
+
+%% Checking generated noise
+
+acxx1 = ac1.psd(plist('Nfft',1e6));
+acxx2 = ac2.psd(plist('Nfft',1e6));
+
+acxy = cohere(ac1,ac2,plist('Nfft',5e4));
+racxy = real(acxy);
+racxy.setName;
+iacxy = imag(acxy);
+iacxy.setName;
+
+iplot(abs(acxx1),abs(CSD(1,1)))
+iplot(abs(acxx2),abs(CSD(2,2)))
+
+iplot(racxy,rcoh,plist('Yscales',{'All','lin'}))
+iplot(iacxy,icoh,plist('Yscales',{'All','lin'}))
+
+%%
+iplot(abs(acxx1),abs(acxx2))
+
+%% Get Coloring filters
+
+cf11 = ac1.procinfo.find('Filt11');
+cf12 = ac1.procinfo.find('Filt12');
+cf21 = ac2.procinfo.find('FILT21');
+cf22 = ac2.procinfo.find('FILT22');
+
+%% Get Theoretical CSD from coloring filters
+
+rc11 = resp(cf11,plist('f',CSD(1,1).x,'bank','parallel'));
+rc12 = resp(cf12,plist('f',CSD(1,1).x,'bank','parallel'));
+rc21 = resp(cf21,plist('f',CSD(1,1).x,'bank','parallel'));
+rc22 = resp(cf22,plist('f',CSD(1,1).x,'bank','parallel'));
+
+
+psd11 = rc11.*conj(rc11)+rc12.*conj(rc12);
+psd22 = rc21.*conj(rc21)+rc22.*conj(rc22);
+psd12 = rc11.*conj(rc21)+rc12.*conj(rc22);
+psd21 = conj(psd12);
+coher = psd12./sqrt(abs(psd11).*abs(psd22));
+rcoher = real(coher);
+rcoher.setName;
+icoher = imag(coher);
+icoher.setName;
+
+% compare with data
+iplot(racxy,rcoher,plist('Yscales',{'All','lin'}))
+iplot(iacxy,icoher,plist('Yscales',{'All','lin'}))
+
+%% Noise whitening
+
+pl = plist(...
+    'csd11', CSD(1,1),...
+    'csd12', CSD(1,2),...
+    'csd21', CSD(2,1),...
+    'csd22', CSD(2,2),...
+    'MaxIter', 75, ...
+    'PoleType', 3, ...
+    'MinOrder', 25, ...
+    'MaxOrder', 40, ...
+    'Weights', 2, ...
+    'Plot', false,...
+    'MSEVARTOL', 1e-2,...
+    'FITTOL', 1e-3,...
+    'UseSym', 0,...
+    'Disp', false);
+
+[aw1,aw2] = whiten2D(ac1,ac2, pl);
+
+aw1.setYunits('m');
+aw2.setYunits('m');
+
+%% Checking whitened noise
+
+awxx1 = aw1.psd(plist('Nfft',1e6));
+awxx2 = aw2.psd(plist('Nfft',1e6));
+%%
+iplot(awxx1,awxx2)
+
+%% Get coherence
+
+awxy = cohere(aw1,aw2,plist('Nfft',6e4));
+rawxy = real(awxy);
+rawxy.setName;
+iawxy = imag(awxy);
+iawxy.setName;
+
+%% Check whitened data coherence
+racxy.setName('a_xy');
+iacxy.setName('a_xy');
+rawxy.setName('wa_xy');
+iawxy.setName('wa_xy');
+iplot(racxy.split(plist('split_type','frequencies','frequencies',[4e-4 5])),rcoh.split(plist('split_type','frequencies','frequencies',[1e-4 5])),rawxy.split(plist('split_type','frequencies','frequencies',[4e-4 5])),plist('Yscales',{'All','lin'},'Ylabels',{'All','Re(coh)'}))
+iplot(iacxy.split(plist('split_type','frequencies','frequencies',[4e-4 5])),icoh.split(plist('split_type','frequencies','frequencies',[1e-4 5])),iawxy.split(plist('split_type','frequencies','frequencies',[4e-4 5])),plist('Yscales',{'All','lin'},'Ylabels',{'All','Im(coh)'}))
+
+%% Get whitening filters
+
+wf11 = aw1.procinfo.find('Filt11');
+wf12 = aw1.procinfo.find('Filt12');
+wf21 = aw2.procinfo.find('FILT21');
+wf22 = aw2.procinfo.find('FILT22');
+
+% filters respose
+rw11 = resp(wf11,plist('f',CSD(1,1).x,'bank','parallel'));
+rw12 = resp(wf12,plist('f',CSD(1,1).x,'bank','parallel'));
+rw21 = resp(wf21,plist('f',CSD(1,1).x,'bank','parallel'));
+rw22 = resp(wf22,plist('f',CSD(1,1).x,'bank','parallel'));
+
+%% Build expected csd and coherence of whitened data
+
+wp11 = conj(rw11).*psd11+psd12.*conj(rw12);
+wp12 = psd11.*conj(rw21)+psd12.*conj(rw22);
+wp21 = psd21.*conj(rw11)+psd22.*conj(rw12);
+wp22 = psd21.*conj(rw21)+psd22.*conj(rw22);
+
+u11 = rw11.*wp11+rw12.*wp21;
+u11.setName;
+u12 = rw11.*wp12+rw12.*wp22;
+u12.setName;
+u21 = rw21.*wp11+rw22.*wp21;
+u21.setName;
+u22 = rw21.*wp12+rw22.*wp22;
+u22.setName;
+wcoher = u12./sqrt(abs(u11).*abs(u22));
+rwcoher = real(wcoher);
+rwcoher.setName;
+iwcoher = imag(wcoher);
+iwcoher.setName;
+
+% compare with data
+iplot(rawxy,rwcoher,plist('Yscales',{'All','lin'}))
+iplot(iawxy,iwcoher,plist('Yscales',{'All','lin'}))
+
+%% Coherence of starting data
+
+axy = cohere(a1,a2,plist('Nfft',1e4));
+raxy = real(axy);
+raxy.setName;
+iaxy = imag(axy);
+iaxy.setName;
+
+% compare with whitened
+iplot(rawxy,raxy,plist('Yscales',{'All','lin'}))
+iplot(iawxy,iaxy,plist('Yscales',{'All','lin'}))
+
+%% Comparing Theoretical data
+
+iplot(rcoher,rwcoher,plist('Yscales',{'All','lin'}))
+iplot(icoher,iwcoher,plist('Yscales',{'All','lin'}))
+
+iplot(rcoh,rwcoher,plist('Yscales',{'All','lin'}))
+iplot(icoh,iwcoher,plist('Yscales',{'All','lin'}))
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+% END
\ No newline at end of file