% 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 vectorsf = logspace(-6,log10(5),300);f = f.';fs = 10;Nsecs = 1e5; % number of secondsNfft = 1e5; % number of samples for the fftpls = plist('Nfft', Nfft,'Order',0); % plist for spectra%% MDC2 Modelsb = ao(plist('built-in','mdc2r2_fd_ltpnoise','f',f,'fs',fs));CSD = [b(1) b(2);conj(b(2)) b(3)];%% Make white noisea1 = 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 plotingiplot(axx)iplot(CSD(1,1),CSD(1,2),CSD(2,1),CSD(2,2))%% Noise generationpl = 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-spectrumacxx1 = 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 cpsdiplot(abs(acxx12),abs(CSD(1,2)))iplot(abs(acxx2),abs(CSD(2,2)))iplot(abs(acxx1),abs(acxx12),abs(acxx2))%% Plotting coherenceiplot(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 2m2=mean(CSD(2,2).data.y(end-5,end))% calculate average on the tail of channel 2m1/m2 % verify that the ratio is near 1%% % ************************************************************************% Some more analysis for testing the accuracy of noise generation procedure% ************************************************************************%% Extracting filters from dataFilt11 = 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 responsestr11 = resp(Filt11,plist('f',f));rFilt11 = tr11(1);for ii = 2:numel(tr11) rFilt11 = rFilt11 + tr11(ii);endrFilt11.setName('rFilt11', 'internal');tr12 = resp(Filt12,plist('f',f));rFilt12 = tr12(1);for ii = 2:numel(tr12) rFilt12 = rFilt12 + tr12(ii);endrFilt12.setName('rFilt12', 'internal');tr21 = resp(Filt21,plist('f',f));rFilt21 = tr21(1);for ii = 2:numel(tr21) rFilt21 = rFilt21 + tr21(ii);endrFilt21.setName('rFilt21', 'internal');tr22 = resp(Filt22,plist('f',f));rFilt22 = tr22(1);for ii = 2:numel(tr22) rFilt22 = rFilt22 + tr22(ii);endrFilt22.setName('rFilt22', 'internal');%% Obtaining transfer functions% calculating transfer functions from datapl = 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 TFspl = 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 + b22b11 = 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-TFsetf11 = 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 responsespl = 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 processicsd11 = 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 AOseigtf11 = 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 processpl = 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% processeigph1 = 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