line source
+ − % 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.3 2009/02/19 17:46:41 luigi Exp $
+ − %
+ −
+ − %% General use variables and vectors
+ −
+ − f = logspace(-6,log10(5),300);
+ − 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','f1',1e-6,'f2',5,'nf',300));
+ − 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', 35, ...
+ − 'MaxOrder', 40, ...
+ − 'Weights', 2, ...
+ − 'Plot', false,...
+ − 'FitTolerance', 2,...
+ − 'RMSEVar', 7,...
+ − 'UseSym', 0,...
+ − 'Disp', false);
+ −
+ − ac = noisegen2D(a, pl);
+ −
+ − %% Checking results and starting data
+ −
+ − % iplot(a)
+ − iplot(ac)
+ −
+ − %% Making cross-spectrum
+ −
+ − acxx = ac.cpsd(pls);
+ − acch = ac.cohere(pls);
+ −
+ − %% Plotting spectra
+ −
+ − % iplot(acxx);
+ −
+ − iplot(abs(acxx(1,1)),abs(CSD(1,1))) % model data need to be multiplied by 2 because acxx is the onesided cpsd
+ − iplot(abs(acxx(1,2)),abs(CSD(1,2)))
+ − iplot(abs(acxx(2,2)),abs(CSD(2,2)))
+ −
+ − iplot(abs(acxx(1,1)),abs(acxx(1,2)),abs(acxx(2,2)))
+ −
+ − %% Plotting coherence
+ −
+ − iplot(acch(2,1),(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
+ − etf = tfe(a,ac,pls);
+ −
+ − %% Comparing Filters Responses with estimated TFs (e-TFs)
+ −
+ − % Comparing filters responses and calculated TFs
+ − pl = plist('Legends', {'Filter Response','e-TF'});
+ − iplot(rFilt11,etf(1,3),pl)
+ − iplot(rFilt12,etf(2,3),pl)
+ − iplot(rFilt21,etf(1,4),pl)
+ − iplot(rFilt22,etf(2,4),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(1,2),pl)
+ − iplot(rFilt12,etf12(1,2),pl)
+ − iplot(rFilt21,etf21(1,2),pl)
+ − iplot(rFilt22,etf22(1,2),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