--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/MDC2/test_ao_noisegen2D_mdc2.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,247 @@
+% A test script for ao/noisegen2D
+% 
+% DESCRIPTION: Run noisegen2D with mdc2 models and test procedure accuracy
+% 
+% L. Ferraioli 04-02-2009
+% 
+% $Id: test_ao_noisegen2D_mdc2.m,v 1.1 2009/04/20 14:31:43 luigi Exp $
+% 
+
+%% General use variables and vectors
+
+f = logspace(-6,log10(5),300);
+f = f.';
+fs = 10;
+Nsecs = 1e5; % number of seconds
+Nfft = 1e5; % number of samples for the fft
+pls  = plist('Nfft', Nfft,'Order',0); % plist for spectra
+
+%% MDC2 Models
+
+b = ao(plist('built-in','mdc2r2_fd_ltpnoise','f',f,'fs',fs));
+CSD = [b(1) b(2);conj(b(2)) b(3)];
+
+%% Make white noise
+
+a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs));
+a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs));
+a1.setYunits('m');
+a2.setYunits('m');
+a = [a1 a2];
+
+% axx = a.cpsd(pls);
+
+%% 
+iplot(a)
+
+%% some ploting
+
+iplot(axx)
+
+iplot(CSD(1,1),CSD(1,2),CSD(2,1),CSD(2,2))
+
+%% Noise generation
+
+pl = plist(...
+    'csd11', CSD(1,1), ...
+    'csd12', CSD(1,2), ...
+    'csd21', CSD(2,1), ...
+    'csd22', CSD(2,2), ...
+    'MaxIter', 80, ...
+    'PoleType', 2, ...
+    'MinOrder', 30, ...
+    'MaxOrder', 40, ...
+    'Weights', 2, ...
+    'Plot', false,...
+    'FITTOL', 1e-5,...
+    'MSEVARTOL', 1e-2,...
+    'UseSym', 0,...
+    'Disp', false);
+
+ac = noisegen2D(a, pl);
+
+%% Checking results and starting data
+
+% iplot(a)
+iplot(ac)
+
+%% Making cross-spectrum
+
+acxx1 = ac(1).psd(pls);
+acxx2 = ac(2).psd(pls);
+acxx12 = cpsd(ac(1),ac(2),pls);
+acch = ac.cohere(pls);
+
+%% Plotting spectra
+
+% iplot(acxx);
+
+iplot(abs(acxx1),abs(CSD(1,1))) % model data need to be multiplied by 2 because acxx is the onesided cpsd
+iplot(abs(acxx12),abs(CSD(1,2)))
+iplot(abs(acxx2),abs(CSD(2,2)))
+
+iplot(abs(acxx1),abs(acxx12),abs(acxx2))
+
+%% Plotting coherence
+
+iplot(acch,(abs(CSD(1,2)).^2)./(CSD(1,1).*CSD(2,2)))
+
+%%
+
+m1=mean(acxx(2,2).data.y(end-10,end)) % calculate average on the tail of channel 2
+m2=mean(CSD(2,2).data.y(end-5,end))% calculate average on the tail of channel 2
+m1/m2 % verify that the ratio is near 1
+
+%% 
+% ************************************************************************
+% Some more analysis for testing the accuracy of noise generation procedure
+% ************************************************************************
+
+%% Extracting filters from data
+
+Filt11 = find(ac(1).procinfo,'Filt11');
+Filt12 = find(ac(1).procinfo,'Filt12');
+Filt21 = find(ac(2).procinfo,'Filt21');
+Filt22 = find(ac(2).procinfo,'Filt22');
+
+%% Calculating filters responses
+
+tr11 = resp(Filt11,plist('f',f));
+rFilt11 = tr11(1);
+for ii = 2:numel(tr11)
+  rFilt11 = rFilt11 + tr11(ii);
+end
+rFilt11.setName('rFilt11', 'internal');
+
+tr12 = resp(Filt12,plist('f',f));
+rFilt12 = tr12(1);
+for ii = 2:numel(tr12)
+  rFilt12 = rFilt12 + tr12(ii);
+end
+rFilt12.setName('rFilt12', 'internal');
+
+tr21 = resp(Filt21,plist('f',f));
+rFilt21 = tr21(1);
+for ii = 2:numel(tr21)
+  rFilt21 = rFilt21 + tr21(ii);
+end
+rFilt21.setName('rFilt21', 'internal');
+
+tr22 = resp(Filt22,plist('f',f));
+rFilt22 = tr22(1);
+for ii = 2:numel(tr22)
+  rFilt22 = rFilt22 + tr22(ii);
+end
+rFilt22.setName('rFilt22', 'internal');
+
+%% Obtaining transfer functions
+
+% calculating transfer functions from data
+pl  = plist('Nfft', Nt);
+etf11 = tfe(a(1),ac(1),pl);
+etf12 = tfe(a(2),ac(1),pl);
+etf21 = tfe(a(1),ac(2),pl);
+etf22 = tfe(a(2),ac(2),pl);
+
+%% Comparing Filters Responses with estimated TFs (e-TFs)
+
+% Comparing filters responses and calculated TFs
+pl = plist('Legends', {'Filter Response','e-TF'});
+iplot(rFilt11,etf11,pl)
+iplot(rFilt12,etf12,pl)
+iplot(rFilt21,etf21,pl)
+iplot(rFilt22,etf22,pl)
+
+%% Filtering data separately
+
+% This operation is performed internally to the noisegen2D. Output data are
+% then obtained by b1 = b11 + b12 and b2 = b21 + b22
+b11 = filter(a1,plist('filter',Filt11,'bank','parallel'));
+b12 = filter(a2,plist('filter',Filt12,'bank','parallel'));
+b21 = filter(a1,plist('filter',Filt21,'bank','parallel'));
+b22 = filter(a2,plist('filter',Filt22,'bank','parallel'));
+
+%% Extracting transfer functions from separately filtered data se-TFs
+
+etf11 = tfe(a1,b11,pls);
+etf12 = tfe(a2,b12,pls);
+etf21 = tfe(a1,b21,pls);
+etf22 = tfe(a2,b22,pls);
+
+%% Comparing separately-estimated TFs (se-TFs) with filter responses
+
+pl = plist('Legends', {'Filter Response','se-TF'});
+iplot(rFilt11,etf11,pl)
+iplot(rFilt12,etf12,pl)
+iplot(rFilt21,etf21,pl)
+iplot(rFilt22,etf22,pl)
+
+%% Comparing filters with TFs obtained by eigendecomposition
+
+% This function output transfer functions as they are obtained by the
+% eigendecomposition process. i.e. before the fitting process
+
+icsd11 = CSD(1,1).data.y*fs/2;
+icsd12 = CSD(1,2).data.y*fs/2;
+icsd21 = CSD(2,1).data.y*fs/2;
+icsd22 = CSD(2,2).data.y*fs/2;
+
+[tf11,tf12,tf21,tf22] = utils.math.eigcsd(icsd11,icsd12,icsd21,icsd22,'USESYM',0,'DIG',50,'OTP','TF');
+
+% Making AOs
+eigtf11 = ao(fsdata(f,tf11,fs));
+eigtf12 = ao(fsdata(f,tf12,fs));
+eigtf21 = ao(fsdata(f,tf21,fs));
+eigtf22 = ao(fsdata(f,tf22,fs));
+
+%% Comparing eig-TFs with output filters
+
+% Compare TFs before and after the fitting process
+
+pl = plist('Legends', {'eig-TF','Filter Response'});
+iplot(eigtf11,rFilt11,pl)
+iplot(eigtf12,rFilt12,pl)
+iplot(eigtf21,rFilt21,pl)
+iplot(eigtf22,rFilt22,pl)
+
+%% Phase difference between eig-TFs and output filters
+
+% checking that phase differences between TFs are preserved by the fitting
+% process
+eigph1 = unwrap(angle(eigtf11)-angle(eigtf21));
+eigph1.setYunits('rad')
+filtph1 = unwrap(angle(rFilt11)-angle(rFilt21));
+filtph1.setYunits('rad')
+eigph2 = unwrap(angle(eigtf22)-angle(eigtf12));
+eigph2.setYunits('rad')
+filtph2 = unwrap(angle(rFilt22)-angle(rFilt12));
+filtph2.setYunits('rad')
+
+pl = plist('Legends',{'eig-TF','Filter'},'YScales',{'All','lin'});
+iplot(eigph1+2*pi,filtph1,pl)
+iplot(eigph2,filtph2,pl)
+
+
+%% Comparing eig-TFs with se-TFs
+
+% Compare eigendecomposition results with separately estimated TFs (se-TFs)
+pl = plist('Legends', {'eig-TF','se-TF'});
+iplot(eigtf11,etf11(1,2),pl)
+iplot(eigtf12,etf12(1,2),pl)
+iplot(eigtf21,etf21(1,2),pl)
+iplot(eigtf22,etf22(1,2),pl)
+
+%% Phase difference between eig-TFs and se-TFs
+
+% checking that phase differences between TFs are preserved by the fitting
+% process also for the filtered data (se-TFs)
+eigph1 = unwrap(angle(eigtf11)-angle(eigtf21));
+filtph1 = unwrap(angle(etf11(1,2))-angle(etf21(1,2)));
+eigph2 = unwrap(angle(eigtf22)-angle(eigtf12));
+filtph2 = unwrap(angle(etf22(1,2))-angle(etf12(1,2)));
+
+pl = plist('Legends',{'eig-TF \Delta\phi','se-TF \Delta\phi'},'YScales',{'All','lin'});
+iplot(eigph1,filtph1,pl)
+iplot(eigph2,filtph2,pl)
+
+% END
\ No newline at end of file