% function test_ao_firwhiten
% A test script for ao/firwhiten
% 
% M Hewitson 21-04-08
% 
% $Id: test_ao_firwhiten.m,v 1.5 2009/08/26 18:27:06 luigi Exp $
% 

mc

%% Make an AO

fs    = 10;
nsecs = 1e3;

% pole zero model
pzm = pzmodel(2, [pz(0.1, 2) pz(0.5, 1)] , [pz(0.8) pz(1, 2)]);

% parameter list for ltpda_noisegen
pl = plist();
pl = append(pl, param('nsecs', nsecs));
pl = append(pl, param('fs', fs));

% calling the noisegenerator
a = fngen(pzm, pl);
% a = set(a, 'yunits', 'V');

% Make series of sine waves

% s = ao(plist('waveform', 'sine wave', ...
%              'A', [0.1 1 0.1], ...
%              'f', [3.3  0.42 1.23], ...
%              'phi', [0 0 0], ...
%              'fs', fs, 'nsecs', nsecs));

s = ao(plist('tsfcn', '0.05*sin(2*pi*1.23*t)', ...
             'fs', fs, 'nsecs', nsecs));
% s = ao(plist('waveform', 'sine wave', ...
%              'A', [0.05], ...
%              'f', [1.23], ...
%              'phi', [0], ...
%              'fs', fs, 'nsecs', nsecs));

% Add signal and noise
a = a + s;

%% Whiten
Nfft = fs*100;
[aw, ffs, nfs] = firwhiten(a, plist('Ntaps', 1024, 'Nfft', Nfft));

% plot filter response
iplot(resp(ffs))


%% Remove end effects

spl = plist('split_type', 'times', 'times', [10 -20]);
a   = timeshift(split(a, spl));
aw  = timeshift(split(aw, spl));

%% Compare time-series
iplot(a,aw, plist('Arrangement', 'subplots', 'XRanges', {'All', [0 10]}))

%% Compare spectra

pxx = pwelch(a, aw, plist('Nfft', fs*100));
iplot(pxx)

%% Help test

fs = 1;
a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e5, 'yunits','m'));

% filter
pl = plist('type','bandpass', ...
  'order',3, ...
  'gain',1, ...
  'fc',[0.03 0.1], ...
  'fs',fs);
ft = miir(pl);

% coloring white noise
af = filter(a, ft);

% firwhiten
[aw, ffs, nfs] = firwhiten(af, plist('Ntaps', 5000, 'Nfft', 1e5, 'BW', 5));

awxx = aw.psd;
afxx = af.psd;

iplot(afxx,awxx)

% END