Mercurial > hg > ltpda
view m-toolbox/test/diagnostics/test_ltpda_arma_freq.m @ 36:5eb86f6881ef database-connection-manager
Remove commented-out code
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Mon, 05 Dec 2011 16:20:06 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
line source
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;