% Test file for ltpda_timedomain fit. This does a noise subtraction from
% engineering model data.
% 
% M Hewitson 18-03-07
% 
% $Id: test_timeDomain.m,v 1.3 2007/12/07 08:54:07 hewitson Exp $
% 

clear all;

%% Load data into analysis objects
disp('*** Loading data....')
x12   = ao('x12.txt'); %x12
phi1  = ao('phi1.txt'); %phi1
eta1  = ao('eta1.txt'); %eta1
phi12 = ao('phi12.txt'); %phi12
eta12 = ao('eta12.txt'); %eta12
disp('*** Done.')

% get the sample frequency for designing bandpass filter
d1 = x12.data;
fs = d1.fs;

%% Load bandpass
pl = plist();
pl = append(pl, param('type', 'bandpass'));
% pl = append(pl, param('fs', fs));
pl = append(pl, param('fc', [0.003 0.03]));
pl = append(pl, param('order', 2));
bp = miir(pl);

%% Filter data and remove startup transients

% parameter list for split() 
% to remove startup transient of filters
spl    = plist(param('times', [1000 2000]));

% parameter list for filter()
fpl    = plist(param('filter', bp));

% filter and split data
x12f   = split(filter(x12, fpl), spl);
phi1f  = split(filter(phi1, fpl), spl);
eta1f  = split(filter(eta1, fpl), spl);
phi12f = split(filter(phi12, fpl), spl);
eta12f = split(filter(eta12, fpl), spl);


%% Make Timedomainfit

coeffsAO = ltpda_timedomainfit(x12f, phi1f, eta1f, phi12f, eta12f);

%% Make linear combination 

x12ns = ltpda_lincom(x12, phi1, eta1, phi12, eta12, coeffsAO);

%% Make LPSD of each

% Window function
w = specwin('Kaiser', 10, 150);

% parameter list for lpsd
pl = plist();
pl = append(pl, param('Kdes', 100));
pl = append(pl, param('Kmin', 2));
pl = append(pl, param('Jdes', 500));
pl = append(pl, param('Win', w));

% LPSD of filtered x12
x12fxx = ltpda_lpsd(x12f, pl);
% LPSD of x12
x12dt = ltpda_polydetrend(x12, plist(param('N', 1)));
x12xx = ltpda_lpsd(x12dt, pl);
% LPSD of noise-subtracted x12
x12nsdt = ltpda_polydetrend(x12ns, plist(param('N', 3)));
x12nsxx = ltpda_lpsd(x12ns, pl);

%% Plot results

figure
plot([x12fxx x12xx x12nsxx])
ltpda_allxaxis(0.001, 1);
subplot(3,1,1:2)
legend('Bandpassed X12', 'X12', 'X12 with noise subtracted');
%%
figure
plot(x12ns.hist)

%%

save(x12fxx, 'x12fxx.xml');
save(x12xx, 'x12xx.xml');
save(x12nsxx, 'x12nsxx.xml');

return

%% Reproduce from history

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

%% now run it
clear;
test

figure
plot(a1)
figure
plot(a1.hist)