--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/test_ao_noisegen2D.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,301 @@
+% A test script for ao/noisegen2D
+% 
+% DESCRIPTION: Run noisegen2D with test data and test procedure accuracy
+% 
+% L. Ferraioli 10-11-08
+% 
+% $Id: test_ao_noisegen2D.m,v 1.9 2010/05/03 19:04:45 luigi Exp $
+% 
+
+%% General use variables and vectors
+
+userdir = 'C:\Users\Luigi'; % You should set your own dir
+
+f = logspace(-6,log10(5),300);
+fs = 10;
+
+Nt = fs*10000;
+
+%% Make white noise
+
+a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 1e5));
+a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 10, 'nsecs', 1e5));
+a = [a1 a2];
+
+%% CSD Noise models
+
+fundir = fullfile(userdir,'ltp_data_analysis\MDCs\MDC1_UTN');
+cf = cd;
+cd(fundir);
+[TF,CSD] = mdc1_tf_models(plist('f',f,'fs',fs));
+cd(cf);
+
+%% Input data psd
+
+a1xx = a1.psd;
+a2xx = a2.psd;
+
+%% Making csd of input data
+
+pl  = plist('Nfft', Nt);
+iCSD = cpsd(a,pl);
+
+%% plot input data psd
+
+iplot(a1xx,a2xx)
+
+%% Plot Model Tfs
+
+iplot(TF(1,1),TF(1,2),TF(2,1),TF(2,2))
+
+%% Plot Model CSD
+
+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', 15, ...
+    'MaxOrder', 45, ...
+    'Weights', 3, ...
+    'FITTOL', 1e-3,...
+    'MSEVARTOL', 1e-1,...
+    'UseSym', 0,...
+    'Plot', false,...
+    'Disp', false);
+
+ac = noisegen2D(a, pl);
+
+%% Checking results and starting data
+
+% iplot(a)
+iplot(ac)
+
+%% Making cross-spectrum
+
+% acxx = ac.cpsd;
+plpsd = plist('navs',2,'order',1,'olap',50);
+acxx1 = ac(1).psd(plpsd);
+acxx2 = ac(2).psd(plpsd);
+
+iplot(acxx1,CSD(1,1),acxx2,CSD(2,2))
+
+%% Plotting spectra
+
+% iplot(acxx);
+
+iplot(acxx(1,1),CSD(1,1))
+iplot(abs(acxx(1,2)),abs(CSD(1,2)))
+iplot(acxx(2,2),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)
+
+%% Comparing starting TFs (s-TFs) with estimated TFs
+
+pl = plist('Legends', {'s-TF','e-TF'});
+iplot(TF(1,1),etf11/sqrt(fs))
+iplot(TF(1,2),etf12/sqrt(fs))
+iplot(TF(2,1),etf21/sqrt(fs))
+iplot(TF(2,2),etf22/sqrt(fs))
+
+%% Building CSD from estimated TFs (e-TFs)
+
+% Output CSD is obtained as
+% eCSD = [TF11 TF12;TF21 TF22]*iCSD*[TF11 TF12;TF21 TF22]'
+eCSD = [etf(1,3) etf(2,3);etf(1,4) etf(2,4)]*iCSD*[etf(1,3) etf(2,3);etf(1,4) etf(2,4)]';
+
+ecsd11 = eCSD(1,1);
+ecsd12 = eCSD(1,2);
+ecsd22 = eCSD(2,2);
+
+% ecsd11 = etf(1,3).*conj(etf(1,3))+etf(2,3).*conj(etf(2,3));
+% ecsd12 = etf(1,3).*conj(etf(1,4))+etf(2,3).*conj(etf(2,4));
+% ecsd22 = etf(2,4).*conj(etf(2,4))+etf(1,4).*conj(etf(1,4));
+
+%% Comparing original CSD with e-TFs CSD
+
+pl = plist('Legends', {'Original CSD','e-TF CSD'});
+iplot(abs(CSD(1,1)),ecsd11/(fs),pl)
+iplot(abs(CSD(1,2)),ecsd12/(fs),pl)
+iplot(abs(CSD(2,2)),ecsd22/(fs),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
+
+pl  = plist('Nfft', Nt);
+etf11 = tfe(a1,b11,pl);
+etf12 = tfe(a2,b12,pl);
+etf21 = tfe(a1,b21,pl);
+etf22 = tfe(a2,b22,pl);
+
+%% 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)
+
+%% Building CSD from se-TFs
+
+% Output CSD is obtained as
+% eCSD = [TF11 TF12;TF21 TF22]*iCSD*[TF11 TF12;TF21 TF22]'
+seCSD = [etf11 etf12;etf21 etf22]*iCSD*[etf11 etf12;etf21 etf22]';
+
+secsd11 = seCSD(1,1);
+secsd12 = seCSD(1,2);
+secsd22 = seCSD(2,2);
+
+% secsd11 = etf11(1,2).*conj(etf11(1,2))+etf12(1,2).*conj(etf12(1,2));
+% secsd12 = etf11(1,2).*conj(etf21(1,2))+etf12(1,2).*conj(etf22(1,2));
+% secsd22 = etf22(1,2).*conj(etf22(1,2))+etf21(1,2).*conj(etf21(1,2));
+
+%% Comparing original CSD with e-TFs CSD
+
+pl = plist('Legends', {'Original CSD','se-TF CSD'});
+iplot(CSD(1,1),secsd11/(fs),pl)
+iplot(CSD(1,2),secsd12/(fs),pl)
+iplot(CSD(2,2),secsd22/(fs),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));
+filtph1 = unwrap(angle(rFilt11)-angle(rFilt21));
+eigph2 = unwrap(angle(eigtf22)-angle(eigtf12));
+filtph2 = unwrap(angle(rFilt22)-angle(rFilt12));
+
+pl = plist('Legends',{'eig-TF \Delta\phi','Filter \Delta\phi'},'YScales',{'All','lin'});
+iplot(eigph1,filtph1+2*pi,pl)
+iplot(eigph2,filtph2+2*pi,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+2*pi,pl)
+iplot(eigph2,filtph2,pl)
+
+% END
\ No newline at end of file