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