function test_ltpda_lcpsd()

% Test ltpda_lcpsd functionality.
% 
% M Hewitson 13-02-07
% 
% $Id: test_ltpda_lcpsd.m,v 1.5 2008/06/13 13:51:46 hewitson Exp $
% 

%% Make test AOs

nsecs = 10000;
fs    = 10;

pl = plist();
pl = append(pl, param('nsecs', nsecs));
pl = append(pl, param('fs', fs));
pl = append(pl, param('tsfcn', 'randn(size(t))'));

a1 = ao(pl);
a2 = ao(pl);

%% Filter one time-series


% Make a filter
pl = plist(param('type', 'bandpass'));
pl = append(pl, param('fs', fs));
pl = append(pl, param('order', 4));
pl = append(pl, param('fc', [0.05 0.2]));
f2 = miir(pl);

% filter the input data
[a3, f2] = filter(a1,plist(param('filter', f2)));

%% make some cross-power

a4 = a3+a2;


%% Make CPSD from a6 to a7

pl = plist();
pl = append(pl, param('Kdes', 100));
pl = append(pl, param('Kmin', 2));
pl = append(pl, param('Jdes', 1000));
pl = append(pl, param('Order', 1));
pl = append(pl, param('Win', specwin('Hanning', 10)));

a8 = lcpsd(a1, a4, pl);
a9 = lcpsd(a1, a4);


%% Plot results and history

iplot(a9(1,2), a8(1,2));
plot(a8(1,2).hist);


%% Reproduce from history

% Write an m-file from AO
ao2m(a8(1,2), 'test.m');

% now run it
clear all;
a_out = test

iplot(a_out(1,2))
plot(a_out(1,2).hist)