Mercurial > hg > ltpda
diff m-toolbox/test/test_matrix_MultiChannelNoise.m @ 0:f0afece42f48
Import.
| author | Daniele Nicolodi <nicolodi@science.unitn.it> |
|---|---|
| date | Wed, 23 Nov 2011 19:22:13 +0100 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/m-toolbox/test/test_matrix_MultiChannelNoise.m Wed Nov 23 19:22:13 2011 +0100 @@ -0,0 +1,197 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% TEST matrix/MultiChannelNoise +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% HISTORY: 23-04-2009 L Ferraioli +% Creation +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% VERSION + +'$Id: test_matrix_fromCSD.m,v 1.3 2009/04/23 16:03:55 luigi Exp $'; + +%% Loading spectra + +load ..\m-toolbox\test\mpsd.mat % load mpsd.mat first column is f then psd1, csd and psd2 + +f = mpsd(:,1); +psd = mpsd(:,2); +fs = 10; + +% 1dim model +mod1D = ao(plist('xvals', f, 'yvals', psd, 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH1')); +mod1D.setName; + +% 2dim model +csd11 = ao(plist('xvals', f, 'yvals', mpsd(:,2), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH1 PSD')); +csd11.setName; +csd12 = ao(plist('xvals', f, 'yvals', mpsd(:,3), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH12 CSD')); +csd12.setName; +csd21 = ao(plist('xvals', f, 'yvals', conj(mpsd(:,3)), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH21 CSD')); +csd21.setName; +csd22 = ao(plist('xvals', f, 'yvals', mpsd(:,4), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH2 PSD')); +csd22.setName; + +mod2D = [csd11 csd12;csd21 csd22]; + +%% +%-------------------------------------------------------------------------- +% 1 Dimensional Noisegen Filter +%-------------------------------------------------------------------------- + +%% 1 dim noisegen filter in s domain - symbolic precision All pass + +% plist for noise generation +plns = plist(... + 'csd', mod1D, ... + 'targetobj', 'parfrac', ... + 'Nsecs', 1e4, ... + 'fs', fs, ... + 'MaxIter', 70, ... + 'PoleType', 3, ... + 'MinOrder', 15, ... + 'MaxOrder', 30, ... + 'Weights', 2, ... + 'Plot', true,... + 'MSEVARTOL', 1e-2,... + 'FITTOL', 1e-3,... + 'UseSym', 'on'); + +na = matrix(plns); + + + +%% 1 dim noisegen filter in z domain - symbolic precision All pass + +% plist for noise generation +plns = plist(... + 'csd', mod1D, ... + 'targetobj', 'miir', ... + 'Nsecs', 1e4, ... + 'fs', fs, ... + 'MaxIter', 70, ... + 'PoleType', 3, ... + 'MinOrder', 10, ... + 'MaxOrder', 30, ... + 'Weights', 2, ... + 'Plot', true,... + 'MSEVARTOL', 1e-2,... + 'FITTOL', 1e-3,... + 'UseSym', 'on'); + +na = matrix(plns); + +%% +%-------------------------------------------------------------------------- +% 2 Dimensional Noisegen Filter +%-------------------------------------------------------------------------- + + +%% 2 dim noisegen filter in s domain - symbolic precision All pass + +% plist for noise generation +plns = plist(... + 'csd', mod2D, ... + 'targetobj', 'parfrac', ... + 'MaxIter', 70, ... + 'PoleType', 3, ... + 'MinOrder', 25, ... + 'MaxOrder', 40, ... + 'Weights', 2, ... + 'Plot', false,... + 'MSEVARTOL', 1e-2,... + 'FITTOL', 1e-4,... + 'UseSym', 'on'); + +na = matrix(plns); + +%% Check response + +na11 = resp(na.objs(1,1),plist('f',f,'bank','parallel')); +na12 = resp(na.objs(1,2),plist('f',f,'bank','parallel')); +na21 = resp(na.objs(2,1),plist('f',f,'bank','parallel')); +na22 = resp(na.objs(2,2),plist('f',f,'bank','parallel')); + +mna = matrix([na11 na12;na21 na22]); + +ECSD = mna*conj(transpose(mna)); + +iplot(csd11,ECSD.objs(1,1)) +iplot(csd12,ECSD.objs(1,2)) +iplot(csd22,ECSD.objs(2,2)) + + +%% 2 dim noisegen filter in z domain - symbolic precision All pass + +% plist for noise generation +plns = plist(... + 'csd', mod2D, ... + 'targetobj', 'miir', ... + 'fs', fs, ... + 'MaxIter', 70, ... + 'PoleType', 3, ... + 'MinOrder', 25, ... + 'MaxOrder', 40, ... + 'Weights', 2, ... + 'Plot', false,... + 'MSEVARTOL', 1e-1,... + 'FITTOL', 1e-4,... + 'UseSym', 'on'); + +na = matrix(plns); + +%% Check response + +na11 = resp(na.objs(1,1).filters,plist('f',f,'bank','parallel')); +na12 = resp(na.objs(1,2).filters,plist('f',f,'bank','parallel')); +na21 = resp(na.objs(2,1).filters,plist('f',f,'bank','parallel')); +na22 = resp(na.objs(2,2).filters,plist('f',f,'bank','parallel')); + +mna = matrix([na11 na12;na21 na22]); + +ECSD = mna*conj(transpose(mna)); + +iplot(csd11,ECSD.objs(1,1)) +iplot(csd12,ECSD.objs(1,2)) +iplot(csd22,ECSD.objs(2,2)) + +%% +%-------------------------------------------------------------------------- +% 2 Dimensional Noise Generation +%-------------------------------------------------------------------------- + +plng = plist('model',na,... + 'Nsecs',2e5,... + 'fs',fs); + +no = matrix(plng); + +%% Check psd + +no1 = no.objs(1); +no2 = no.objs(2); + +no1xx = no1.psd; +no2xx = no2.psd; + +iplot(no1xx,csd11,no2xx,csd22) + +%% Check cross-coherence + +coh = csd12./sqrt(csd11.*csd22); +rcoh = real(coh); +icoh = imag(coh); + +no12xx = cpsd(no1,no2,plist('nfft',1e4)); +no1xx = no1.psd(plist('nfft',1e4)); +no2xx = no2.psd(plist('nfft',1e4)); + +ecoh = no12xx./(sqrt(no1xx.*no2xx)); +recoh = real(ecoh); +iecoh = imag(ecoh); + +iplot(recoh,rcoh,plist('Yscales',{'All','lin'})) +iplot(iecoh,icoh,plist('Yscales',{'All','lin'})) + + + +