line source
% 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