Mercurial > hg > ltpda
diff m-toolbox/test/test_matrix_fromCSD.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_fromCSD.m Wed Nov 23 19:22:13 2011 +0100 @@ -0,0 +1,287 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% TEST matrix/fromCSD +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% HISTORY: 23-04-2009 L Ferraioli +% Creation +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% VERSION '$Id: test_matrix_fromCSD.m,v 1.3 2009/11/06 17:00:30 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 - double precision All pass + +% plist for noise generation +plns = plist(... + 'csd', mod1D, ... + 'targetobj', 'parfrac', ... + 'Nsecs', 1e4, ... + 'fs', fs, ... + 'MaxIter', 70, ... + 'PoleType', 3, ... + 'MinOrder', 10, ... + 'MaxOrder', 30, ... + 'Weights', 2, ... + 'Plot', true,... + 'MSEVARTOL', 1e-2,... + 'FITTOL', 1e-3,... + 'UseSym', 'off'); + +na = matrix(plns); + +%% 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 - double 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', 'off'); + +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 - double precision All pass + +% plist for noise generation +plns = plist(... + 'csd', mod2D, ... + 'targetobj', 'parfrac', ... + 'MaxIter', 70, ... + 'PoleType', 3, ... + 'MinOrder', 20, ... + 'MaxOrder', 40, ... + 'Weights', 2, ... + 'Plot', true,... + 'MSEVARTOL', 1e-2,... + 'FITTOL', 1e-3,... + 'UseSym', 'off'); + +na = matrix(plns); + +%% 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', 20, ... + 'MaxOrder', 30, ... + 'Weights', 2, ... + 'Plot', true,... + 'MSEVARTOL', 1e-2,... + 'FITTOL', 1e-3,... + 'UseSym', 'on'); + +na = matrix(plns); + +%% Check response + +na11 = resp(na.objs(1,1),plist('f',f)); +na12 = resp(na.objs(1,2),plist('f',f)); +na21 = resp(na.objs(2,1),plist('f',f)); +na22 = resp(na.objs(2,2),plist('f',f)); + +est_csd11 = na11.*conj(na11) + na12.*conj(na12); +est_csd11.setName; +est_csd12 = na11.*conj(na21) + na12.*conj(na22); +est_csd12.setName; +est_csd22 = na12.*conj(na12) + na22.*conj(na22); +est_csd22.setName; + +iplot(csd11,est_csd11) +iplot(csd12,est_csd12) +iplot(csd22,est_csd22) + + +%% 2 dim noisegen filter in z domain - double precision All pass + +% plist for noise generation +plns = plist(... + 'csd', mod2D, ... + 'targetobj', 'miir', ... + 'fs', fs, ... + 'MaxIter', 70, ... + 'PoleType', 3, ... + 'MinOrder', 25, ... + 'MaxOrder', 30, ... + 'Weights', 2, ... + 'Plot', true,... + 'MSEVARTOL', 1e-2,... + 'FITTOL', 1e-3,... + 'UseSym', 'off'); + +na = matrix(plns); + +%% 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', 35, ... + 'Weights', 2, ... + 'Plot', false,... + 'MSEVARTOL', 1e-1,... + 'FITTOL', 1e-3,... + '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,(2/fs).*ECSD.objs(1,1)) +iplot(csd12,(2/fs).*ECSD.objs(1,2)) +iplot(csd22,(2/fs).*ECSD.objs(2,2)) + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +