Mercurial > hg > ltpda

comparison m-toolbox/test/test_ao_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 ao/fromCSD
3 %
4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 % HISTORY: 23-04-2009 L Ferraioli
6 % Creation
7 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
8 %% VERSION
9
10 '$Id: test_ao_fromCSD.m,v 1.3 2009/04/23 16:03:55 luigi Exp $';
11
12 %% Loading spectra
13
14 load ..\m-toolbox\test\mpsd.mat % load mpsd.mat first column is f then psd1, csd and psd2
15
16 f = mpsd(:,1);
17 psd = mpsd(:,2);
18 fs = 10;
19
20 % 1dim model
21 mod1D = ao(plist('xvals', f, 'yvals', psd, 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH1'));
22 mod1D.setName;
23
24 % 2dim model
25 csd11 = ao(plist('xvals', f, 'yvals', mpsd(:,2), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH1 PSD'));
26 csd11.setName;
27 csd12 = ao(plist('xvals', f, 'yvals', mpsd(:,3), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH12 CSD'));
28 csd12.setName;
29 csd21 = ao(plist('xvals', f, 'yvals', conj(mpsd(:,3)), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH21 CSD'));
30 csd21.setName;
31 csd22 = ao(plist('xvals', f, 'yvals', mpsd(:,4), 'fs', fs, 'dtype', 'fsdata','description','MDC1 IFO CH2 PSD'));
32 csd22.setName;
33
34 mod2D = [csd11 csd12;csd21 csd22];
35
36
37 %% Generate 1D Noise
38
39 % plist for noise generation
40
41 plns = plist(...
42 'csd', mod1D, ...
43 'Nsecs', 1e4, ...
44 'fs', fs, ...
45 'xunits', 's', ...
46 'yunits', '', ...
47 'MaxIter', 50, ...
48 'PoleType', 3, ...
49 'MinOrder', 7, ...
50 'MaxOrder', 35, ...
51 'Weights', 2, ...
52 'Plot', false,...
53 'MSEVARTOL', 1e-2,...
54 'FITTOL', 1e-2);
55
56 na = ao(plns);
57
58 %% Compare with Noisegen1D
59
60 aw = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e4));
61
62 pl = plist(...
63 'model', mod1D, ...
64 'MaxIter', 50, ...
65 'PoleType', 3, ...
66 'MinOrder', 7, ...
67 'MaxOrder', 35, ...
68 'Weights', 2, ...
69 'Plot', false,...
70 'Disp', false,...
71 'MSEVARTOL', 1e-2,...
72 'FITTOL', 1e-2);
73
74 nat = noisegen1D(aw, pl);
75
76 %% Make psd
77
78 naxx = na.lpsd;
79 natxx = nat.lpsd;
80
81
82 iplot(naxx,natxx,mod1D)
83
84 %% 2 dim noise generation
85
86 % plist for noise generation
87
88 plns = plist(...
89 'csd', mod2D, ...
90 'Nsecs', 1e4, ...
91 'fs', fs, ...
92 'xunits', 's', ...
93 'yunits', '', ...
94 'MaxIter', 70, ...
95 'PoleType', 3, ...
96 'MinOrder', 20, ...
97 'MaxOrder', 40, ...
98 'Weights', 2, ...
99 'Plot', false,...
100 'MSEVARTOL', 1e-2,...
101 'FITTOL', 1e-2);
102
103 na = ao(plns);
104
105 %% Compare with Noisegen2D
106
107 aw1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e4));
108 aw2 = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', 1e4));
109
110 pl = plist(...
111 'csd11', csd11, ...
112 'csd12', csd12, ...
113 'csd21', csd21, ...
114 'csd22', csd22, ...
115 'MaxIter', 70, ...
116 'PoleType', 3, ...
117 'MinOrder', 20, ...
118 'MaxOrder', 40, ...
119 'Weights', 2, ...
120 'Plot', false,...
121 'Disp', false,...
122 'MSEVARTOL', 1e-2,...
123 'FITTOL', 1e-2);
124
125 [nat1,nat2] = noisegen2D(aw1,aw2, pl);
126
127 %% Make psd
128
129 naxx1 = na(1).lpsd;
130 natxx1 = nat1.lpsd;
131 naxx2 = na(2).lpsd;
132 natxx2 = nat2.lpsd;
133
134
135 iplot(naxx1,natxx1,csd11)
136 iplot(naxx2,natxx2,csd22)
137
138
139 %% 3 Dim noise generation
140
141 % sampling frequency
142 fs = 1; % Hz
143 % frequency vector
144 f = logspace(-6,log10(0.5),300);
145 f = f.';
146
147 %%% Set the 3 channel model
148 h11 = pzmodel(plist('GAIN', [0.01], 'POLES', [pz(9.9999999999999995e-007,NaN) pz(0.01,NaN) pz(0.02,NaN)], 'ZEROS', [pz(0.0001,NaN) pz(0.002,NaN)]));
149 h12 = pzmodel(plist('GAIN', [0.0001], 'POLES', [pz(0.10000000000000001,NaN) pz(0.0030000000000000001,NaN) pz(0.0050000000000000001,NaN) pz(1.0000000000000001e-005,NaN)], 'ZEROS', [pz(0.00050000000000000001,NaN) pz(0.001,NaN)]));
150 h13 = pzmodel(plist('GAIN', [0.00001], 'POLES', [pz(5.0000000000000004e-006,NaN) pz(0.02,NaN)], 'ZEROS', pz(0.00040000000000000002,NaN)));
151
152 h21 = pzmodel(plist('GAIN', [0.00001], 'POLES', [pz(5.0000000000000004e-006,NaN) pz(0.10000000000000001,NaN) pz(0.050000000000000003,NaN)], 'ZEROS', [pz(0.002,NaN) pz(0.00020000000000000001,NaN)]));
153 h22 = pzmodel(plist('GAIN', [0.001], 'POLES', [pz(0.01,NaN) pz(9.9999999999999995e-008,NaN) pz(0.002,NaN) pz(0.001,NaN)], 'ZEROS', [pz(1.0000000000000001e-005,NaN) pz(2.0000000000000002e-005,NaN) pz(5.0000000000000002e-005,NaN) pz(0.20000000000000001,NaN)]));
154 h23 = pzmodel(plist('GAIN', [0.0001], 'POLES', [pz(0.01,NaN) pz(0.029999999999999999,NaN) pz(0.10000000000000001,NaN) pz(5.0000000000000002e-005,NaN)], 'ZEROS', [pz(0.0001,NaN) pz(0.0050000000000000001,NaN)]));
155
156 h31 = pzmodel(plist('GAIN', [0.001], 'POLES', [pz(0.01,NaN) pz(0.029999999999999999,NaN) pz(0.10000000000000001,NaN) pz(0.00050000000000000001,NaN)], 'ZEROS', [pz(0.0001,NaN) pz(0.0050000000000000001,NaN)]));
157 h32 = pzmodel(plist('GAIN', [0.00001], 'POLES', [pz(0.01,NaN) pz(0.029999999999999999,NaN) pz(0.10000000000000001,NaN) pz(0.00050000000000000001,NaN) pz(0.00029999999999999997,NaN)], 'ZEROS', [pz(0.002,NaN) pz(0.059999999999999998,NaN)]));
158 h33 = pzmodel(plist('GAIN', [0.05], 'POLES', [pz(0.0001,NaN) pz(9.9999999999999995e-008,NaN) pz(0.002,NaN)], 'ZEROS', [pz(3.0000000000000001e-005,NaN) pz(5.0000000000000002e-005,NaN) pz(0.10000000000000001,NaN)]));
159
160 %%% TFs resps
161 % filters response
162 plresp = plist('f',f);
163
164 rh11 = resp(h11,plresp);
165 rh12 = resp(h12,plresp);
166 rh13 = resp(h13,plresp);
167
168 rh21 = resp(h21,plresp);
169 rh22 = resp(h22,plresp);
170 rh23 = resp(h23,plresp);
171
172 rh31 = resp(h31,plresp);
173 rh32 = resp(h32,plresp);
174 rh33 = resp(h33,plresp);
175
176 %%% Spectra calculation with Kay style method
177 % Note that the definition of cross-spectral matrix follows that of s M
178 % Kay, Modern Spectral Estimation
179
180 % csd11
181 G11 = rh11.y.*conj(rh11.y) + rh12.y.*conj(rh12.y) + rh13.y.*conj(rh13.y);
182 G11 = ao(plist('xvals', f, 'yvals', G11, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
183 G11.setName;
184
185 % csd12
186 G12 = conj(rh11.y).*rh21.y + conj(rh12.y).*rh22.y + conj(rh13.y).*rh23.y;
187 G12 = ao(plist('xvals', f, 'yvals', G12, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
188 G12.setName;
189
190 % csd13
191 G13 = conj(rh11.y).*rh31.y + conj(rh12.y).*rh32.y + conj(rh13.y).*rh33.y;
192 G13 = ao(plist('xvals', f, 'yvals', G13, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
193 G13.setName;
194
195 % csd21
196 G21 = conj(G12);
197 G21.setName;
198
199 % csd22
200 G22 = rh21.y.*conj(rh21.y) + rh22.y.*conj(rh22.y) + rh23.y.*conj(rh23.y);
201 G22 = ao(plist('xvals', f, 'yvals', G22, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
202 G22.setName;
203
204 % csd23
205 G23 = conj(rh21.y).*rh31.y + conj(rh22.y).*rh32.y + conj(rh23.y).*rh33.y;
206 G23 = ao(plist('xvals', f, 'yvals', G23, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
207 G23.setName;
208
209 % csd31
210 G31 = conj(G13);
211 G31.setName;
212
213 % csd32
214 G32 = conj(G23);
215 G32.setName;
216
217 % csd33
218 G33 = conj(rh31.y).*rh31.y + conj(rh32.y).*rh32.y + conj(rh33.y).*rh33.y;
219 G33 = ao(plist('xvals', f, 'yvals', G33, 'fs', fs, 'dtype', 'fsdata', 'yunits', unit('m^2').*unit('Hz^-1')));
220 G33.setName;
221
222 %%% Build CSD matrix
223 mod3D = [G11 G12 G13;G21 G22 G23;G31 G32 G33];
224
225 %% noise generation
226
227 % plist for noise generation
228
229 plns = plist(...
230 'csd', mod3D, ...
231 'Nsecs', 1e4, ...
232 'fs', fs, ...
233 'xunits', 's', ...
234 'yunits', '', ...
235 'MaxIter', 90, ...
236 'PoleType', 3, ...
237 'MinOrder', 20, ...
238 'MaxOrder', 50, ...
239 'Weights', 2, ...
240 'Plot', false,...
241 'MSEVARTOL', 1e-2,...
242 'FITTOL', 1e-2);
243
244 na = ao(plns);
245
246 %% Make psd
247
248 naxx1 = na(1).psd;
249
250 naxx2 = na(2).psd;
251
252 naxx3 = na(3).psd;
253
254
255
256 iplot(naxx1,G11)
257 iplot(naxx2,G22)
258 iplot(naxx3,G33)
259
260
261
262
263
264
265