Mercurial > hg > ltpda
view m-toolbox/test/test_ao_noisegen2D_mdc2.m @ 10:75007001cbfe database-connection-manager
Check for binary only objects
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
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