Mercurial > hg > ltpda
diff m-toolbox/test/test_ao_buildWhitener1D.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_ao_buildWhitener1D.m Wed Nov 23 19:22:13 2011 +0100 @@ -0,0 +1,146 @@ +% A test script for ao/buildWhitener1D +% +% M Hewitson 10-11-08 +% +% $Id: test_ao_buildWhitener1D.m,v 1.1 2010/05/04 07:14:17 mauro Exp $ +% + +%% Make test data +fs = 10; +a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 10000, 'yunits', 'N')); + +% filter +pzm = pzmodel(1e-2, {0.01}, {0.1}); +ft = miir(pzm,plist('fs',fs)); + +af = filter(a, ft); + +%% + +rsp = pzm.resp; +rsmm = resp(ft,plist('f',logspace(-4,log10(5),100))); +iplot(rsp,rsmm) + +%% +pl_psd = plist('scale', 'PSD', 'order', 1, 'navs', 16); +axx = a.psd(pl_psd); +afxx = af.psd(pl_psd); + +iplot(axx,afxx) + +%% Build whitening filter it with no model + +pl = plist(... + 'MaxIter', 30, ... + 'PoleType', 1, ... + 'MinOrder', 2, ... + 'MaxOrder', 9, ... + 'Weights', 2, ... + 'Plot', false,... + 'Disp', false,... + 'MSEVARTOL', 1e-2,... + 'FITTOL', 1e-1); % tolerance on MSE Value + +awf = buildWhitener1D(af, pl); + +%% Whiten with a model + +mdl = (abs(rsp).^2).*(2/fs); % this corresponds to the theoretical psd of the data +mdl.setYunits(axx.yunits); + +pl = plist(... + 'fs', fs, ... + 'MaxIter', 50, ... + 'PoleType', 2, ... + 'MinOrder', 2, ... + 'MaxOrder', 9, ... + 'Weights', 2, ... + 'Plot', false,... + 'Disp', false,... + 'MSEVARTOL', 1e-2,... + 'FITTOL', 1e-2); % tolerancee on MSE Value + +awmf = buildWhitener1D(mdl, pl); + + +%% +aw = filter(af, awf); +awm = filter(af, awmf); +iplot(aw, awm) + +%% + +awxx = aw.psd(pl_psd); +awmxx = awm.psd(pl_psd); + +iplot(axx, afxx, awxx, awmxx); + +%% Testing flatness + +svec = [axx.data.y afxx.data.y awxx.data.y awmxx.data.y]; +fc = utils.math.spflat(svec); + +%% Working with multiple inputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +fs = 10; +a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 10000, 'yunits', 'V')); +a2 = ao(plist('tsfcn', 'randn(size(t)).''', 'fs', fs, 'nsecs', 10000, 'yunits', 'T')); + +% filter +pzm = pzmodel(1e-2, {0.01}, {0.1}); +ft = miir(pzm,plist('fs',fs)); + +af1 = filter(a1, ft); +af2 = filter(a2, ft); + +pl = plist(... + 'MaxIter', 30, ... + 'PoleType', 2, ... + 'MinOrder', 2, ... + 'MaxOrder', 9, ... + 'Weights', 2, ... + 'Plot', false,... + 'Disp', false,... + 'MSEVARTOL', 1e-2,... + 'FITTOL', 1e-2); % tolerance on fit residuals spectral flatness + +[awf1,awf2] = buildWhitener1D(af1,af2,pl); + +%% + +iplot(filter(af1, awf1),filter(af2, awf2)) + +%% help test + +% Generate white noise +fs = 1; +a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e5, 'yunits','m')); + +% filter +ft = miir(plist('type','bandpass','order',3,'gain',1,'fc',[0.03 0.1],'fs',fs)); + +% coloring white noise +af = filter(a, ft); + +% Whitening colored noise +pl = plist(... + 'MaxIter', 30, ... + 'PoleType', 2, ... + 'MinOrder', 9, ... + 'MaxOrder', 15, ... + 'Weights', 2, ... + 'Plot', false,... + 'Disp', false,... + 'MSEVARTOL', 1e-1,... + 'FITTOL', 5e-2); % tolerance on fit residuals spectral flatness + +aw = buildWhitener1D(af,pl); + +% Calculate psd of colored and whitened data +afxx = af.psd; +awxx = af.filter(aw).psd; + +% plotting +iplot(afxx,awxx) + +% END \ No newline at end of file