0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % Test script for ao/fftfilt
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 % L Ferraioli 09-10-08
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 % $Id: test_ao_cov.m,v 1.2 2008/12/11 08:25:15 hewitson Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 s = ao(plist('tsfcn', '3.*sin(2.*pi.*0.01.*t) + randn(size(t))', 'fs', 10, 'nsecs', 1e4));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10 %% Test with a miir filter
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 % get a miir filter
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 plf = plist('type','lowpass',...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 'order',4,...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 'fs',10,...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 'fc',1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 filt = miir(plf);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 % do classical fitering
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20 fs1 = filter(s,filt);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 iplot(s,fs1)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 % do fftfilt
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 fs2 = fftfilt(s,filt);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 % compare results
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 iplot(fs1,fs2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 iplot(fs1./fs2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29 iplot(fs1-fs2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31 %% Test with a miir filter
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33 % get a miir filter
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 plf = plist('type','lowpass',...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35 'order',4,...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 'fs',10,...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 'fc',1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 filt = miir(plf);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40 % do classical fitering
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 fs1 = filter(s,filt);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 iplot(s,fs1)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 % do fftfilt
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45 plfft = plist('Npad',5.*length(s.data.y));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46 fs2 = fftfilt(s,filt);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48 % compare results
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 iplot(fs1,fs2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50 iplot(fs1./fs2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 iplot(fs1-fs2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53 %%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55 % imp = zeros(10000,1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56 % imp(1) = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 % imp_resp = filter(num,den,imp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61 % fs1 = filter(num,den,s.y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 % fs2 = conv(s.y,imp_resp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63 % X = fft([s.y;zeros(length(imp_resp)-1,1)]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 % Y = fft([imp_resp;zeros(length(s.y)-1,1)]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65 % fs3 = ifft(X.*Y);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 %% Test pzmodel fftfilt
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69 poles = [pz(1e-3,2) pz(1e-2)];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70 zeros = [pz(3e-3,3) pz(5e-2)];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71 pzm = pzmodel(10, poles, zeros);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73 % do fftfilt
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 plfft = plist('Npad',5.*length(s.data.y));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 fs = fftfilt(s,pzm);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77 tf = tfe(s,fs,plist('navs',2));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78 rsp = resp(pzm,plist('f',logspace(-4,log10(5))));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79 iplot(tf,rsp)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81 %% Test smodel
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83 str = 'a.*(2.*i.*pi.*f).^2 + b.*2.*i.*pi.*f + c';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 mod = smodel(str);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85 mod.setParams({'a','b','c'},{1,2,3});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86 mod.setXvar('f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88 freq = logspace(-5,log10(5),100);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 mod.setXvals(freq);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90 mod.setXunits('Hz');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91 mod.setYunits('kg m s^-2');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 % eval smodel
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94 em = eval(mod,plist('output x',freq,'type','fsdata'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96 % make a time series
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97 dt = ao.randn(10000,10);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99 % do fftfilt
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
100 sdt = fftfilt(dt,mod);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
101
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
102 % test output
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
103 tdt = tfe(dt,sdt);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
104 iplot(em,tdt)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
105
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
106 %% Test initial conditions
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
107
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
108 str = 'a.*(2.*i.*pi.*f).^2 + b.*2.*i.*pi.*f + c';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
109 mod = smodel(str);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
110 mod.setParams({'a','b','c'},{1,2,3});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
111 mod.setXvar('f');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
112
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
113 freq = logspace(-5,log10(5),100);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
114 mod.setXvals(freq);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
115 mod.setXunits('Hz');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
116 mod.setYunits('kg m s^-2');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
117
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
118 % eval smodel
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
119 em = eval(mod,plist('output x',freq,'type','fsdata'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
120
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
121 % make a time series
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
122 dt = ao.randn(10000,10);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
123
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
124 % do fftfilt
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
125 sdt = fftfilt(dt,mod);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
126
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
127 % do fftfilt with initial conditions assuming a 2nd order equation
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
128 sdt_inCond = fftfilt(dt,mod,plist('initial conditions',[2000,1000]));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
129
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
130 % test output
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
131 iplot(sdt_inCond,sdt)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
132 iplot(sdt_inCond-sdt,plist('yscales',{'all','log'}))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
133
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
134
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
135
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
136
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
137
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
138
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
139
|