function test_ltpda_arma_freq()
% A test script for the ltpda_arma_freq function
% 
% miquel 25-02-08
% 
% $Id: test_ltpda_arma_freq.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-4*an1;

pn2 = plist('waveform', 'noise','type','Uniform','fs', fs, 'nsecs', nsecs);
an2 = ao(pn2);

aout = aout + 1e-4*an2;


%% Fourier transform input varialbles

% Initialise variables
in=ain.data.y;
out=aout.data.y;

% Fourier Transform
L=length(in);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
NFFT=length(in);

x_ = fft(in,NFFT)/L;
y_ = fft(out,NFFT)/L;

x = 2*abs(x_(1:NFFT/2)); % single-sided amplitude spectrum.
y = 2*abs(y_(1:NFFT/2)); % single-sided amplitude spectrum.

f = fs/2*linspace(0,1,NFFT/2)';

Px=ao(fsdata(f,x,fs));
Py=ao(fsdata(f,y,fs));

%% Find ARMA parameters for the transfer function

% Parameter list for ltpda_arma_freq
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_freq
ao_arma = ltpda_arma_freq(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');

% Compute time domain fit
flt = 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',flt)));

% Compute frequency domain fit
% [Pfit, Pf] = filter(Px,plist(param('filter',f)));  % problems with ao.filter
yfit= freqfun([ao_arma.data.y(1:p)],[1 ao_arma.data.y(p+1:p+q)],x,f,fs); % using adhoc function
Pfit=ao(fsdata(f,yfit,fs));

%  Plot time domain
iplot(ain,aout,afit,plist('Legends', {'Input','Output','Fit'}))

%  Plot frequency domain
iplot(Px,Py,Pfit,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');


% Compute filter
function F= freqfun(ps,qs,xdata,f,fs)
Hnum=0;
Hden=0;
for j = 1:length(ps)
Hnum = Hnum+ps(j)*exp(-(j-1)*i*2*pi*f/fs);
end
for j = 1:length(qs)
Hden = Hden+qs(j)*exp(-(j-1)*i*2*pi*f/fs);
end
F=abs(Hnum./Hden).*xdata;