Mercurial > hg > ltpda

comparison m-toolbox/test/diagnostics/test_ltpda_arma_freq.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 function test_ltpda_arma_freq()
2 % A test script for the ltpda_arma_freq function
3 %
4 % miquel 25-02-08
5 %
6 % $Id: test_ltpda_arma_freq.m,v 1.1 2008/03/04 12:55:05 miquel Exp $
7 %
8
9 %% Make test AOs
10
11 fs=1;
12 nsecs=10000;
13
14 pl = plist();
15 pl = append(pl, param('nsecs',nsecs));
16 pl = append(pl, param('fs',fs));
17 pl = append(pl,param('tsfcn','(heaviside(t-1000.5)).*(1-(exp(-(t-1000.5)/300))+(exp(-(t-1000.5)/500)-1))'));
18
19 ain = ao(pl);
20
21 % Filter input time-series
22 pl = append(pl, param('type', 'highpass'));
23 pl = append(pl, param('ripple', 5e-1));
24 pl = append(pl, param('order', 1));
25 pl = append(pl, param('gain', 1e-1));
26 pl = append(pl, param('fc', 1e-3));
27
28 f = miir(pl);
29 [aout, fout] = filter(ain,plist(param('filter',f)));
30
31 % Add independent noise to both channels
32 pn1 = plist('waveform', 'noise','type','Uniform','fs', fs, 'nsecs', nsecs);
33 an1 = ao(pn1);
34
35 ain = ain + 1e-4*an1;
36
37 pn2 = plist('waveform', 'noise','type','Uniform','fs', fs, 'nsecs', nsecs);
38 an2 = ao(pn2);
39
40 aout = aout + 1e-4*an2;
41
42
43 %% Fourier transform input varialbles
44
45 % Initialise variables
46 in=ain.data.y;
47 out=aout.data.y;
48
49 % Fourier Transform
50 L=length(in);
51 NFFT = 2^nextpow2(L); % Next power of 2 from length of y
52 NFFT=length(in);
53
54 x_ = fft(in,NFFT)/L;
55 y_ = fft(out,NFFT)/L;
56
57 x = 2*abs(x_(1:NFFT/2)); % single-sided amplitude spectrum.
58 y = 2*abs(y_(1:NFFT/2)); % single-sided amplitude spectrum.
59
60 f = fs/2*linspace(0,1,NFFT/2)';
61
62 Px=ao(fsdata(f,x,fs));
63 Py=ao(fsdata(f,y,fs));
64
65 %% Find ARMA parameters for the transfer function
66
67 % Parameter list for ltpda_arma_freq
68 pl = append(pl, param('MaxIter',1e3));
69 pl = append(pl, param('MaxFunEvals',1e3));
70 pl = append(pl, param('TolX',1e-7));
71 pl = append(pl, param('TolFun',1e-7));
72 pl = append(pl, param('ARMANum',2));
73 pl = append(pl, param('ARMADen',1));
74
75 % Use ltpda_arma_freq
76 ao_arma = ltpda_arma_freq(ain, aout, pl);
77 %
78 disp(sprintf('\r! Filter parameters:\r'))
79 disp(sprintf('%d\t', fout.a,fout.b(2:length(fout.b))))
80 disp(sprintf('\r! Estimated parameters:\r'))
81 disp(sprintf('%d\t',ao_arma.data.y(1:length(ao_arma.data.y)-1)))
82
83 %% Compute fitted output
84 p = find(pl, 'ARMANum');
85 q = find(pl, 'ARMADen');
86
87 % Compute time domain fit
88 flt = miir([ao_arma.data.y(1:p)],[1 ao_arma.data.y(p+1:p+q)],find(pl,'fs'));
89 [afit, ffit] = filter(ain,plist(param('filter',flt)));
90
91 % Compute frequency domain fit
92 % [Pfit, Pf] = filter(Px,plist(param('filter',f))); % problems with ao.filter
93 yfit= freqfun([ao_arma.data.y(1:p)],[1 ao_arma.data.y(p+1:p+q)],x,f,fs); % using adhoc function
94 Pfit=ao(fsdata(f,yfit,fs));
95
96 % Plot time domain
97 iplot(ain,aout,afit,plist('Legends', {'Input','Output','Fit'}))
98
99 % Plot frequency domain
100 iplot(Px,Py,Pfit,plist('Legends', {'Input','Output','Fit'}))
101
102
103 % Plot history
104 figure
105 plot(ao_arma.hist)
106
107 %% Reproduce from history
108
109 % Write an m-file from AO
110 ao2m(ao_arma, 'ltpda_test.m');
111
112 % now run it
113 clear all;
114 a_out = ltpda_test;
115
116 % iplot(a_out)
117
118 figure
119 plot(a_out.hist)
120
121 %% Delete test files
122 delete('ltpda_test.m');
123
124
125 % Compute filter
126 function F= freqfun(ps,qs,xdata,f,fs)
127 Hnum=0;
128 Hden=0;
129 for j = 1:length(ps)
130 Hnum = Hnum+ps(j)*exp(-(j-1)*i*2*pi*f/fs);
131 end
132 for j = 1:length(qs)
133 Hden = Hden+qs(j)*exp(-(j-1)*i*2*pi*f/fs);
134 end
135 F=abs(Hnum./Hden).*xdata;
136
137
138