Mercurial > hg > ltpda

comparison m-toolbox/test/test_matrix_fromCSD.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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 % TEST matrix/fromCSD
3 %
4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 % HISTORY: 23-04-2009 L Ferraioli
6 % Creation
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8 %% VERSION '$Id: test_matrix_fromCSD.m,v 1.3 2009/11/06 17:00:30 luigi Exp $';
9
10 %% Loading spectra
11
12 load ..\m-toolbox\test\mpsd.mat % load mpsd.mat first column is f then psd1, csd and psd2
13
14 f = mpsd(:,1);
15 psd = mpsd(:,2);
16 fs = 10;
17
18 % 1dim model
19 mod1D = ao(plist('xvals', f, 'yvals', psd, 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH1'));
20 mod1D.setName;
21
22 % 2dim model
23 csd11 = ao(plist('xvals', f, 'yvals', mpsd(:,2), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH1 PSD'));
24 csd11.setName;
25 csd12 = ao(plist('xvals', f, 'yvals', mpsd(:,3), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH12 CSD'));
26 csd12.setName;
27 csd21 = ao(plist('xvals', f, 'yvals', conj(mpsd(:,3)), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH21 CSD'));
28 csd21.setName;
29 csd22 = ao(plist('xvals', f, 'yvals', mpsd(:,4), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH2 PSD'));
30 csd22.setName;
31
32 mod2D = [csd11 csd12;csd21 csd22];
33
34 %%
35 %--------------------------------------------------------------------------
36 % 1 Dimensional Noisegen Filter
37 %--------------------------------------------------------------------------
38
39 %% 1 dim noisegen filter in s domain - double precision All pass
40
41 % plist for noise generation
42 plns = plist(...
43 'csd', mod1D, ...
44 'targetobj', 'parfrac', ...
45 'Nsecs', 1e4, ...
46 'fs', fs, ...
47 'MaxIter', 70, ...
48 'PoleType', 3, ...
49 'MinOrder', 10, ...
50 'MaxOrder', 30, ...
51 'Weights', 2, ...
52 'Plot', true,...
53 'MSEVARTOL', 1e-2,...
54 'FITTOL', 1e-3,...
55 'UseSym', 'off');
56
57 na = matrix(plns);
58
59 %% 1 dim noisegen filter in s domain - symbolic precision All pass
60
61 % plist for noise generation
62 plns = plist(...
63 'csd', mod1D, ...
64 'targetobj', 'parfrac', ...
65 'Nsecs', 1e4, ...
66 'fs', fs, ...
67 'MaxIter', 70, ...
68 'PoleType', 3, ...
69 'MinOrder', 15, ...
70 'MaxOrder', 30, ...
71 'Weights', 2, ...
72 'Plot', true,...
73 'MSEVARTOL', 1e-2,...
74 'FITTOL', 1e-3,...
75 'UseSym', 'on');
76
77 na = matrix(plns);
78
79 %% 1 dim noisegen filter in z domain - double precision All pass
80
81 % plist for noise generation
82 plns = plist(...
83 'csd', mod1D, ...
84 'targetobj', 'miir', ...
85 'Nsecs', 1e4, ...
86 'fs', fs, ...
87 'MaxIter', 70, ...
88 'PoleType', 3, ...
89 'MinOrder', 10, ...
90 'MaxOrder', 30, ...
91 'Weights', 2, ...
92 'Plot', true,...
93 'MSEVARTOL', 1e-2,...
94 'FITTOL', 1e-3,...
95 'UseSym', 'off');
96
97 na = matrix(plns);
98
99 %% 1 dim noisegen filter in z domain - symbolic precision All pass
100
101 % plist for noise generation
102 plns = plist(...
103 'csd', mod1D, ...
104 'targetobj', 'miir', ...
105 'Nsecs', 1e4, ...
106 'fs', fs, ...
107 'MaxIter', 70, ...
108 'PoleType', 3, ...
109 'MinOrder', 10, ...
110 'MaxOrder', 30, ...
111 'Weights', 2, ...
112 'Plot', true,...
113 'MSEVARTOL', 1e-2,...
114 'FITTOL', 1e-3,...
115 'UseSym', 'on');
116
117 na = matrix(plns);
118
119 %%
120 %--------------------------------------------------------------------------
121 % 2 Dimensional Noisegen Filter
122 %--------------------------------------------------------------------------
123
124
125 %% 2 dim noisegen filter in s domain - double precision All pass
126
127 % plist for noise generation
128 plns = plist(...
129 'csd', mod2D, ...
130 'targetobj', 'parfrac', ...
131 'MaxIter', 70, ...
132 'PoleType', 3, ...
133 'MinOrder', 20, ...
134 'MaxOrder', 40, ...
135 'Weights', 2, ...
136 'Plot', true,...
137 'MSEVARTOL', 1e-2,...
138 'FITTOL', 1e-3,...
139 'UseSym', 'off');
140
141 na = matrix(plns);
142
143 %% 2 dim noisegen filter in s domain - symbolic precision All pass
144
145 % plist for noise generation
146 plns = plist(...
147 'csd', mod2D, ...
148 'targetobj', 'parfrac', ...
149 'MaxIter', 70, ...
150 'PoleType', 3, ...
151 'MinOrder', 20, ...
152 'MaxOrder', 30, ...
153 'Weights', 2, ...
154 'Plot', true,...
155 'MSEVARTOL', 1e-2,...
156 'FITTOL', 1e-3,...
157 'UseSym', 'on');
158
159 na = matrix(plns);
160
161 %% Check response
162
163 na11 = resp(na.objs(1,1),plist('f',f));
164 na12 = resp(na.objs(1,2),plist('f',f));
165 na21 = resp(na.objs(2,1),plist('f',f));
166 na22 = resp(na.objs(2,2),plist('f',f));
167
168 est_csd11 = na11.*conj(na11) + na12.*conj(na12);
169 est_csd11.setName;
170 est_csd12 = na11.*conj(na21) + na12.*conj(na22);
171 est_csd12.setName;
172 est_csd22 = na12.*conj(na12) + na22.*conj(na22);
173 est_csd22.setName;
174
175 iplot(csd11,est_csd11)
176 iplot(csd12,est_csd12)
177 iplot(csd22,est_csd22)
178
179
180 %% 2 dim noisegen filter in z domain - double precision All pass
181
182 % plist for noise generation
183 plns = plist(...
184 'csd', mod2D, ...
185 'targetobj', 'miir', ...
186 'fs', fs, ...
187 'MaxIter', 70, ...
188 'PoleType', 3, ...
189 'MinOrder', 25, ...
190 'MaxOrder', 30, ...
191 'Weights', 2, ...
192 'Plot', true,...
193 'MSEVARTOL', 1e-2,...
194 'FITTOL', 1e-3,...
195 'UseSym', 'off');
196
197 na = matrix(plns);
198
199 %% 2 dim noisegen filter in z domain - symbolic precision All pass
200
201 % plist for noise generation
202 plns = plist(...
203 'csd', mod2D, ...
204 'targetobj', 'miir', ...
205 'fs', fs, ...
206 'MaxIter', 70, ...
207 'PoleType', 3, ...
208 'MinOrder', 25, ...
209 'MaxOrder', 35, ...
210 'Weights', 2, ...
211 'Plot', false,...
212 'MSEVARTOL', 1e-1,...
213 'FITTOL', 1e-3,...
214 'UseSym', 'on');
215
216 na = matrix(plns);
217
218 %% Check response
219
220 na11 = resp(na.objs(1,1).filters,plist('f',f,'bank','parallel'));
221 na12 = resp(na.objs(1,2).filters,plist('f',f,'bank','parallel'));
222 na21 = resp(na.objs(2,1).filters,plist('f',f,'bank','parallel'));
223 na22 = resp(na.objs(2,2).filters,plist('f',f,'bank','parallel'));
224
225 mna = matrix([na11 na12;na21 na22]);
226
227 ECSD = mna*conj(transpose(mna));
228
229 iplot(csd11,(2/fs).*ECSD.objs(1,1))
230 iplot(csd12,(2/fs).*ECSD.objs(1,2))
231 iplot(csd22,(2/fs).*ECSD.objs(2,2))
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287