% 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