Mercurial > hg > ltpda

comparison m-toolbox/test/LTPDA_training/TrainingSessionAll.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 % Training session
2 %
3 % 1) Topic 1 - The basics of LTPDA
4 % 2) Topic 2 - Pre-processing of data
5 % 3) Topic 3 - Spectral Analysis
6 % 4) Topic 4 - Transfer function models and digital filtering
7 % 5) Topic 5 - Model fitting
8 %
9 % HISTORY:
10 %
11 % $Id: TrainingSessionAll.m,v 1.10 2010/02/24 17:55:39 ingo Exp $
12 %
13 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14
15 %%
16 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
17 %%% Topic 1 - The basics of LTPDA %%%
18 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
19
20 % Clear all variables and figures
21 mc()
22 %%
23 % Reading the interferometer data
24 pl_file = plist('FILENAME', 'ifo_temp_example/ifo_training.dat', ...
25 'TYPE', 'tsdata', ...
26 'COLUMNS', [1 2], ...
27 'XUNITS', 's', ...
28 'YUNITS', '', ...
29 'ROBUST', 'no', ...
30 'DESCRIPTION', 'Interferometer data');
31 ifo = ao(pl_file);
32 ifo.setName(); % Set the object name to the variable name (here: 'ifo')
33
34 % Calibrating the interferometer data
35 lambda = 1064e-9;
36 pl_scale = plist('factor', lambda/(2*pi), 'yunits', 'm');
37 ifo.scale(pl_scale);
38 ifo.setName(); % Set the object name to the variable name (here: 'ifo')
39
40 % Plot ifo
41 ifo.iplot(plist('XUNITS', 'h'));
42
43 % Save ifo to 'ifo_temp_example/ifo_disp.xml'
44 ifo.save('ifo_temp_example/ifo_disp.xml');
45
46
47 % Reading the interferometer data
48 pl_fileT = plist('FILENAME', 'ifo_temp_example/temp_training.dat', ...
49 'TYPE', 'tsdata', ...
50 'COLUMNS', [1 2], ...
51 'XUNITS', 's', ...
52 'YUNITS', 'degC', ...
53 'ROBUST', 'no', ...
54 'DESCRIPTION', 'Temperature data');
55 temp = ao(pl_fileT);
56
57 % Add offset
58 temp.offset(plist('offset', 273.15));
59 temp.setYunits('K');
60 temp.setName(); % Set the object name to the variable name (here: 'temp')
61
62 % Plot Tcel
63 temp.iplot(plist('XUNITS', 'h'));
64
65 % Save Tcel
66 temp.save(plist('filename', 'ifo_temp_example/temp_kelvin.xml'));
67
68
69 %%
70 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
71 %%% Topic 2 - Pre-processing of data %%%
72 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
73
74 % Clear all variables and figures
75 mc()
76
77 % Load data from topic 1
78 ifo = ao('ifo_temp_example/ifo_disp.xml');
79 temp = ao('ifo_temp_example/temp_kelvin.xml');
80
81 % plot the data
82 pl_plot1 = plist('arrangement', 'subplots');
83 iplot(ifo, temp, pl_plot1)
84
85 pl_plot2 = plist('ARRANGEMENT', 'subplots', ...
86 'LINESTYLES', {'none','none'}, ...
87 'MARKERS', {'+','+'}, ...
88 'XRANGES', {'all', [200 210]}, ...
89 'YRANGES', {[2e-7 3e-7], [200 350]});
90
91 iplot(ifo, temp, pl_plot2)
92
93 % The temperature data is unevenly sampled.
94 dt = diff(temp.x);
95 min(dt)
96 max(dt)
97
98 % Run 'data fixer' method ao/consolidate
99 [temp_fixed ifo_fixed] = consolidate(temp, ifo, plist('fs',1));
100
101 % Plot fixed data
102 iplot(ifo_fixed, temp_fixed, pl_plot1);
103 iplot(ifo_fixed, temp_fixed, pl_plot2);
104
105 % Save fixed data
106 save(temp_fixed,'ifo_temp_example/temp_fixed.xml');
107 save(ifo_fixed,'ifo_temp_example/ifo_fixed.xml');
108
109 %%
110 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
111 %%% Topic 3 - Spectral Analysis %%%
112 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
113
114 % Clear all variables and figures
115 mc()
116
117 % Get the consolidated data
118 % Using the xml format
119 T_filename = 'ifo_temp_example/temp_fixed.xml';
120 x_filename = 'ifo_temp_example/ifo_fixed.xml';
121
122 pl_load_T = plist('filename', T_filename);
123 pl_load_x = plist('filename', x_filename);
124
125 % Build the data aos
126 T = ao(pl_load_T);
127 x = ao(pl_load_x);
128
129 % PSD
130 x_psd = lpsd(x);
131 x_psd.setName('Interferometer');
132
133 T_psd = lpsd(T);
134 T_psd.setName('Temperature');
135
136 % Plot estimated PSD
137 pl_plot = plist('Arrangement', 'subplots', 'LineStyles', {'-','-'},'Linecolors', {'b', 'r'});
138 iplot(sqrt(x_psd), sqrt(T_psd), pl_plot);
139
140 % Skip some IFO glitch from the consolidation
141 pl_split = plist('split_type', 'interval', ...
142 'start_time', x.t0 + 40800, ...
143 'end_time', x.t0 + 193500);
144 x_red = split(x, pl_split);
145 T_red = split(T, pl_split);
146
147 % PSD
148 x_red_psd = lpsd(x_red);
149 x_red_psd.setName('Interferometer');
150
151 T_red_psd = lpsd(T_red);
152 T_red_psd.setName('Temperature');
153
154 % Plot estimated PSD
155 pl_plot = plist('Arrangement', 'stacked', 'LineStyles', {'-','-'},'Linecolors', {'b', 'r'});
156 iplot(sqrt(x_psd), sqrt(x_red_psd), pl_plot);
157 iplot(sqrt(T_psd), sqrt(T_red_psd), pl_plot);
158
159 % CPSD estimate
160 CTx = lcpsd(T_red, x_red);
161 CxT = lcpsd(x_red, T_red);
162
163 % Plot estimated CPSD
164 iplot(CTx);
165 iplot(CxT);
166
167 % Coherence estimate
168 coh = lcohere(T_red, x_red);
169
170 % Plot estimated cross-coherence
171 iplot(coh, plist('YScales', 'lin'))
172
173 % transfer function estimate
174 tf = ltfe(T_red, x_red);
175
176 % Plot estimated TF
177 iplot(tf);
178
179 % Noise projection in frequency domain
180 proj = T_red_psd.*(abs(tf)).^2;
181 proj.simplifyYunits;
182 proj.setName('temp. contrib. projection')
183
184 % Plotting the noise projection in frequency domain
185 iplot(x_red_psd, proj);
186
187 % Save the PSD data
188 % Plists for the xml format
189 pl_save_x_PSD = plist('filename', 'ifo_temp_example/ifo_psd.xml');
190 pl_save_T_PSD = plist('filename', 'ifo_temp_example/T_psd.xml');
191 pl_save_xT_CPSD = plist('filename', 'ifo_temp_example/ifo_T_cpsd.xml');
192 pl_save_xT_cohere = plist('filename', 'ifo_temp_example/ifo_T_cohere.xml');
193 pl_save_xT_TFE = plist('filename', 'ifo_temp_example/T_ifo_tf.xml');
194
195 x_red_psd.save(pl_save_x_PSD);
196 T_red_psd.save(pl_save_T_PSD);
197 CxT.save(pl_save_xT_CPSD);
198 coh.save(pl_save_xT_cohere);
199 tf.save(pl_save_xT_TFE);
200
201 %%
202 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
203 %%% Topic 4 - Transfer function models and digital filtering %%%
204 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
205
206 mc()
207
208 % Temperature noise PZMODEL
209 TMP = pzmodel(10,1e-5,[]);
210 TMP.setName();
211 TMP.setOunits('K');
212
213 % Interferometer noise PZMODEL
214 IFO = pzmodel(1e-3, {0.4}, []);
215 IFO.setName();
216 IFO.setOunits('rad');
217
218 % Temperature to interferometer coupling PZMODEL
219 K2RAD = pzmodel(1e-1, {5e-4}, []);
220 K2RAD.setName();
221 K2RAD.setOunits('rad');
222 K2RAD.setIunits('K');
223
224 % Plot the response
225 pl = plist('f1',1e-5,'f2',0.01);
226 resp(K2RAD*TMP,IFO,pl);
227
228 % Discretize the three transfer (TMP,IFO,K2RAD) with the MIIR constructor
229 pl_miir = plist('fs', 1);
230 TMPd = miir(TMP, pl_miir);
231 IFOd = miir(IFO, pl_miir);
232 K2RADd = miir(K2RAD, pl_miir);
233
234 % Generate white noise with the AO constructor
235 pl_ao = plist('tsfcn', 'randn(size(t))', ...
236 'fs', 1, ...
237 'nsecs', 250000);
238 WN1 = ao(pl_ao);
239 WN2 = ao(pl_ao);
240
241 % Filter white noise WN1 with the TMP filter
242 T = filter(WN1,TMPd);
243
244 % Filter white noise WN2 with the IFO filter
245 T2 = filter(WN2, IFOd);
246
247 % Filter white noise WN2 with the TMP and the K2RAD filter
248 T3 = filter(WN1, TMPd, K2RADd, plist('bank','serial'));
249
250 % Add Noise
251 IFO = T2 + T3;
252
253 % Split data stream
254 pl_split = plist('times', [1e5 2e5]);
255 IFO = IFO.split(pl_split);
256 IFO.setName('Interferometer');
257 T = T.split(pl_split);
258 T.setName('Temperature');
259
260 % Plot
261 pl_plot1 = plist('ARRANGEMENT', 'subplots');
262 IFO.iplot(pl_plot1, T);
263
264 % Compute power spectral estimates for the temperature and interferometric data
265 pl_lpsd = plist('order', 1, 'scale', 'ASD');
266 lpsd_T = lpsd(T, pl_lpsd);
267 lpsd_ifo = lpsd(IFO, pl_lpsd);
268 iplot(lpsd_T, lpsd_ifo, plist('Arrangement', 'subplots'))
269
270 % Compute transfer function estimate for the temperature and interferometric data
271 tfe_T = ltfe(T, IFO);
272 tfe_T.setName('Transfer function');
273
274 pl = plist('f1',1e-5,'f2',1);
275 iplot(tfe_T, resp(K2RAD,pl));
276
277 % Compute projection
278 Projection = abs(tfe_T).*lpsd_T;
279 Projection.simplifyYunits;
280 Projection.setName;
281 % Plot against interferometer noise
282 iplot(lpsd_ifo,Projection);
283
284 %%
285 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
286 %%% Topic 5 - Model fitting %%%
287 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
288
289 % Clear all variables and figures
290 mc()
291
292 % Load test data from topic 1
293 ifo = ao(plist('filename', 'ifo_temp_example/ifo_fixed.xml'));
294 ifo.setName;
295 T = ao(plist('filename', 'ifo_temp_example/temp_fixed.xml'));
296 T.setName;
297
298 % Split out the good part of the data
299 pl_split = plist('split_type', 'interval', ...
300 'start_time', ifo.t0 + 40800, ...
301 'end_time', ifo.t0 + 193500);
302
303 ifo_red = split(ifo, pl_split);
304 T_red = split(T, pl_split);
305
306 % Plot
307 iplot(ifo_red,T_red,plist('arrangement', 'subplots'))
308
309 % Load transfer function from topic 3
310 tf = ao('ifo_temp_example/T_ifo_tf.xml');
311
312 % split the transfer function to extract only meaningful data
313 tfsp = split(tf,plist('frequencies', [2e-5 1e-3]));
314 iplot(tf,tfsp)
315
316 % force zDomainFit to fit a stable model
317 plfit = plist('FS',1, ...
318 'AutoSearch','off', ...
319 'StartPolesOpt','clin',...
320 'maxiter',20, ...
321 'minorder',3, ...
322 'maxorder',3, ...
323 'weightparam','abs', ...
324 'Plot','on', ...
325 'ForceStability','on',...
326 'CheckProgress','off');
327
328 fobj = zDomainFit(tfsp,plfit);
329
330 fobj.filters.setIunits('K');
331 fobj.filters.setOunits('m');
332
333 % Detrend after the filtering
334 ifoT = filter(T_red,fobj,plist('bank','parallel'));
335 ifoT.detrend(plist('order',0));
336 ifoT.simplifyYunits;
337 ifoT.setName;
338
339 % Substract temperature contribution from measured interferometer data
340 ifonT = ifo_red - ifoT;
341 ifonT.setName;
342
343 % Plot data
344 iplot(ifo_red,ifoT,ifonT)
345
346 % LPSD
347 ifoxx = ifo_red.lpsd;
348 ifonTxx = ifonT.lpsd;
349 iplot(ifoxx,ifonTxx)
350
351
352
353
354
355
356
357
358
359
360
361