line source
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');