--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/test_ao_whiten1D.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,149 @@
+% 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