--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/diagnostics/test_ltpda_arma_time.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,92 @@
+function test_ltpda_arma_time()
+% A test script for the ltpda_arma_time function
+% 
+% miquel 20-02-08
+% 
+% $Id: test_ltpda_arma_time.m,v 1.1 2008/03/04 12:55:05 miquel Exp $
+% 
+
+%% Make test AOs
+
+fs=1;
+nsecs=10000;
+
+pl = plist();
+pl = append(pl, param('nsecs',nsecs));
+pl = append(pl, param('fs',fs));
+pl = append(pl,param('tsfcn','(heaviside(t-1000.5)).*(1-(exp(-(t-1000.5)/300))+(exp(-(t-1000.5)/500)-1))'));
+
+ain = ao(pl);
+
+% Filter input time-series
+pl = append(pl, param('type', 'highpass'));
+pl = append(pl, param('ripple', 5e-1));
+pl = append(pl, param('order', 1));
+pl = append(pl, param('gain', 1e-1));
+pl = append(pl, param('fc', 1e-3));
+
+f = miir(pl);
+[aout, fout] = filter(ain,plist(param('filter',f)));
+
+% Add independent noise to both channels
+pn1 = plist('waveform', 'noise','type','Uniform','fs', fs, 'nsecs', nsecs);
+an1 = ao(pn1);
+
+ain = ain + 1e-3*an1;
+
+pn2 = plist('waveform', 'noise','type','Uniform','fs', fs, 'nsecs', nsecs);
+an2 = ao(pn2);
+
+aout = aout + 1e-3*an2;
+
+
+%% Find ARMA parameters for the transfer function
+
+% Parameter list for ltpda_arma_time
+pl = append(pl, param('MaxIter',1e3));
+pl = append(pl, param('MaxFunEvals',1e3));
+pl = append(pl, param('TolX',1e-7));
+pl = append(pl, param('TolFun',1e-7));
+pl = append(pl, param('ARMANum',2));
+pl = append(pl, param('ARMADen',1));
+
+% Use ltpda_arma_time
+ao_arma = ltpda_arma_time(ain, aout, pl);
+%  
+disp(sprintf('\r! Filter parameters:\r')) 
+disp(sprintf('%d\t', fout.a,fout.b(2:length(fout.b))))
+disp(sprintf('\r! Estimated parameters:\r')) 
+disp(sprintf('%d\t',ao_arma.data.y(1:length(ao_arma.data.y)-1)))
+
+
+%  Compute fitted output
+p = find(pl, 'ARMANum'); 
+q = find(pl, 'ARMADen'); 
+
+f = miir([ao_arma.data.y(1:p)],[1 ao_arma.data.y(p+1:p+q)],find(pl,'fs'));
+[afit, ffit] = filter(ain,plist(param('filter',f)));
+
+%  Plot
+iplot(ain,aout,afit,plist('Legends', {'Input','Output','Fit'}))
+
+% %% Plot history 
+% figure
+% plot(ao_arma.hist)
+% 
+% %% Reproduce from history
+% 
+% % Write an m-file from AO
+% ao2m(ao_arma, 'ltpda_test.m');
+% 
+% % now run it
+% clear all;
+% a_out = ltpda_test;
+% 
+% % iplot(a_out)
+% 
+% figure
+% plot(a_out.hist)
+
+%% Delete test files
+% delete('ltpda_test.m');
+