Mercurial > hg > ltpda

comparison m-toolbox/test/MDC2/test_mdc2_filtimpresp.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 % A test script for to check the filters impulse response
2 %
3 %
4 % L. Ferraioli 04-02-2009
5 %
6 % $Id: test_mdc2_filtimpresp.m,v 1.3 2009/04/20 14:34:12 luigi Exp $
7 %
8
9 %% General use variables and vectors
10
11 fs = 1;
12 f = logspace(-7,log10(fs/2),500);
13 f = f.';
14 Nsecs = 1e6; % number of seconds
15
16 %% MDC2 Models
17
18 b = ao(plist('built-in','mdc2r2_fd_ltpnoise','f',f));
19
20 %%
21
22 iplot(b)
23
24 %%
25
26 CSD = [b(1) b(2);conj(b(2)) b(3)];
27 MSC = b(2).*conj(b(2))./(b(1).*b(3));
28
29 %% some ploting
30
31 iplot(CSD(1,1),CSD(2,2),MSC)
32 iplot(CSD(1,1),CSD(2,2),CSD(1,2),CSD(2,1))
33
34 %% Make white noise
35
36 a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
37 a2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
38
39 %% Noise generation
40
41 plng = plist(...
42 'csd11', CSD(1,1), ...
43 'csd12', CSD(1,2), ...
44 'csd21', CSD(2,1), ...
45 'csd22', CSD(2,2), ...
46 'MaxIter', 90, ...
47 'PoleType', 2, ...
48 'MinOrder', 35, ...
49 'MaxOrder', 50, ...
50 'Weights', 2, ...
51 'Plot', false,...
52 'MSEVARTOL', 1e-2,...
53 'FITTOL', 1e-5,...
54 'UseSym', 0,...
55 'Disp', false);
56
57 [ac1,ac2] = noisegen2D(a1,a2, plng);
58
59 % extracting filters
60 h11 = ac1.procinfo.find('Filt11');
61 h12 = ac1.procinfo.find('Filt12');
62 h21 = ac2.procinfo.find('Filt21');
63 h22 = ac2.procinfo.find('Filt22');
64
65 %% Calculate Impulse response with LTPDA filter
66
67 % LTPDA filter initialize the filters before the filtering process
68
69 imp = zeros(Nsecs,1);
70 imp(1,1) = 1;
71 imp = tsdata(imp,fs);
72 imp = ao(imp);
73 imp.setYunits('m');
74 imp.setXunits('s');
75
76 irh11 = filter(imp,h11);
77 irh12 = filter(imp,h12);
78 irh21 = filter(imp,h21);
79 irh22 = filter(imp,h22);
80
81 irh11.setName;
82 irh12.setName;
83 irh21.setName;
84 irh22.setName;
85
86 % iplot(irh11,irh12,irh21,irh22)
87 %% Plot impulse response
88
89 iplot(irh11,plist('Yscales','linear'))
90 iplot(irh12,plist('Yscales','linear'))
91 iplot(irh21,plist('Yscales','linear'))
92 iplot(irh22,plist('Yscales','linear'))
93
94 %% Impluse response with Matlab filter
95
96 % No filter initialization
97
98 idat = [1; zeros(999999,1)];
99
100 Nfilt11 = numel(h11);
101 num11 = h11(1).a;
102 den11 = h11(1).b;
103 irh11 = filter(num11,den11,idat);
104
105 for ii = 2:Nfilt11
106 num11 = h11(ii).a;
107 den11 = h11(ii).b;
108 todat = filter(num11,den11,idat);
109 irh11 = irh11 + todat;
110 end
111
112 Nfilt12 = numel(h12);
113 num12 = h12(1).a;
114 den12 = h12(1).b;
115 irh12 = filter(num12,den12,idat);
116
117 for ii = 2:Nfilt12
118 num12 = h12(ii).a;
119 den12 = h12(ii).b;
120 todat = filter(num12,den12,idat);
121 irh12 = irh12 + todat;
122 end
123
124 Nfilt21 = numel(h21);
125 num21 = h21(1).a;
126 den21 = h21(1).b;
127 irh21 = filter(num21,den21,idat);
128
129 for ii = 2:Nfilt21
130 num21 = h21(ii).a;
131 den21 = h21(ii).b;
132 todat = filter(num21,den21,idat);
133 irh21 = irh21 + todat;
134 end
135
136 Nfilt22 = numel(h22);
137 num22 = h22(1).a;
138 den22 = h22(1).b;
139 irh22 = filter(num22,den22,idat);
140
141 for ii = 2:Nfilt22
142 num22 = h22(ii).a;
143 den22 = h22(ii).b;
144 todat = filter(num22,den22,idat);
145 irh22 = irh22 + todat;
146 end
147
148 %%
149
150 figure
151 plot(irh11)
152
153 figure
154 plot(irh12)
155
156 figure
157 plot(irh21)
158
159 figure
160 plot(irh22)
161
162 %% Filter state
163
164 ac11 = filter(a1,h11);
165 h11 = ac11.procinfo.find('Filters'); % prooagate histout
166 ac12 = filter(a2,h12);
167 h12 = ac12.procinfo.find('Filters');
168 ac21 = filter(a1,h21);
169 h21 = ac21.procinfo.find('Filters');
170 ac22 = filter(a2,h22);
171 h22 = ac22.procinfo.find('Filters');
172
173 iplot(ac11)
174 iplot(ac12)
175 iplot(ac21)
176 iplot(ac22)
177
178 %% Filter state evolution
179
180 Nrun = 20;
181
182 for jj = 1:Nrun % filter state evolution
183 disp(num2str(jj))
184 ac11 = filter(a1,h11);
185 h11 = ac11.procinfo.find('Filters'); % prooagate histout
186 ac12 = filter(a2,h12);
187 h12 = ac12.procinfo.find('Filters');
188 ac21 = filter(a1,h21);
189 h21 = ac21.procinfo.find('Filters');
190 ac22 = filter(a2,h22);
191 h22 = ac22.procinfo.find('Filters');
192 end
193 %%
194 iplot(ac11)
195 iplot(ac12)
196 iplot(ac21)
197 iplot(ac22)
198
199 %%
200
201 ac1 = ac11 + ac12;
202 ac2 = ac21 + ac22;
203
204 iplot(ac1)
205 iplot(ac2)
206
207 %%
208
209 ac1xx = ac1.lpsd(plist('Nfft',1e5,'order',2));
210 ac2xx = ac2.lpsd(plist('Nfft',1e5,'order',2));
211
212 %%
213 iplot(CSD(1,1),ac1xx,CSD(2,2),ac2xx)
214
215