+
−% A test script for ao/whiten1D
+
−% 
+
−% M Hewitson 10-11-08
+
−% 
+
−% $Id: test_ao_whiten1D.m,v 1.14 2010/05/05 04:19:42 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)
+
−
+
−%% Whiten it with no model
+
−
+
−pl = plist(...
+
−    'model', [], ...
+
−    'MaxIter', 30, ...
+
−    'PoleType', 1, ...
+
−    'MinOrder', 2, ...
+
−    'MaxOrder', 9, ...
+
−    'Weights', 2, ...
+
−    'Plot', false,...
+
−    'Disp', false,...
+
−    'MSEVARTOL', 1e-2,...
+
−    'FITTOL', 1e-1); % tolerance on MSE Value
+
−  
+
−aw = whiten1D(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, ...
+
−    'model', mdl, ...
+
−    'MaxIter', 50, ...
+
−    'PoleType', 2, ...
+
−    'MinOrder', 2, ...
+
−    'MaxOrder', 9, ...
+
−    'Weights', 2, ...
+
−    'Plot', false,...
+
−    'Disp', false,...
+
−    'MSEVARTOL', 1e-2,...
+
−    'FITTOL', 1e-2); % tolerance on MSE Value
+
−
+
−awm = whiten1D(af, pl);
+
−
+
−
+
−%%
+
−
+
−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(...
+
−    'model', [], ...
+
−    '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
+
−  
+
−[aw1,aw2] = whiten1D(af1,af2,pl);
+
−
+
−%%
+
−
+
−iplot(aw1,aw2)
+
−
+
−%% 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(...
+
−    'model', [], ...
+
−    '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 = whiten1D(af,pl);
+
−
+
−% Calculate psd of colored and whitened data
+
−afxx = af.psd;
+
−awxx = aw.psd;
+
−
+
−% plotting
+
−iplot(afxx,awxx)
+
−
+
−% END