Mercurial > hg > ltpda

comparison m-toolbox/html_help/help/ug/whitening.html @ 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 <!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
2 "http://www.w3.org/TR/1999/REC-html401-19991224/loose.dtd">
3
4 <html lang="en">
5 <head>
6 <meta name="generator" content=
7 "HTML Tidy for Mac OS X (vers 1st December 2004), see www.w3.org">
8 <meta http-equiv="Content-Type" content=
9 "text/html; charset=us-ascii">
10
11 <title>Noise whitening (LTPDA Toolbox)</title>
12 <link rel="stylesheet" href="docstyle.css" type="text/css">
13 <meta name="generator" content="DocBook XSL Stylesheets V1.52.2">
14 <meta name="description" content=
15 "Presents an overview of the features, system requirements, and starting the toolbox.">
16 </head>
17
18 <body>
19 <a name="top_of_page" id="top_of_page"></a>
20
21 <p style="font-size:1px;">&nbsp;</p>
22
23 <table class="nav" summary="Navigation aid" border="0" width=
24 "100%" cellpadding="0" cellspacing="0">
25 <tr>
26 <td valign="baseline"><b>LTPDA Toolbox</b></td><td><a href="../helptoc.html">contents</a></td>
27
28 <td valign="baseline" align="right"><a href=
29 "gapfill.html"><img src="b_prev.gif" border="0" align=
30 "bottom" alt="Data gap filling"></a>&nbsp;&nbsp;&nbsp;<a href=
31 "sigproc.html"><img src="b_next.gif" border="0" align=
32 "bottom" alt="Signal Processing in LTPDA"></a></td>
33 </tr>
34 </table>
35
36 <h1 class="title"><a name="f3-12899" id="f3-12899"></a>Noise whitening</h1>
37 <hr>
38
39 <p>
40
41 <!-- ================================================== -->
42 <!-- BEGIN CONTENT FILE -->
43 <!-- ================================================== -->
44 <!-- ===== link box: Begin ===== -->
45 <p>
46 <table border="1" width="80%">
47 <tr>
48 <td>
49 <table border="0" cellpadding="5" class="categorylist" width="100%">
50 <colgroup>
51 <col width="37%"/>
52 <col width="63%"/>
53 </colgroup>
54 <tbody>
55 <tr valign="top">
56 <td>
57 <a href="#WhitenIntro">Introduction</a>
58 </td>
59 <td>Noise whitening in LTPDA.</td>
60 </tr>
61 <tr valign="top">
62 <td>
63 <a href="#WhitenAlgo">Algorithm</a>
64 </td>
65 <td>Whitening Algorithms.</td>
66 </tr>
67 <tr valign="top">
68 <td>
69 <a href="#Whiten1D">1D data</a>
70 </td>
71 <td>Whitening noise in one-dimensional data.</td>
72 </tr>
73 <tr valign="top">
74 <td>
75 <a href="#Whiten2D">2D data</a>
76 </td>
77 <td>Whitening noise in two-dimensional data.</td>
78 </tr>
79 </tbody>
80 </table>
81 </td>
82 </tr>
83 </table>
84 </p>
85 <!-- ===== link box: End ====== -->
86
87
88 <!-- ===== Intro ====== -->
89 <h2><a name="WhitenIntro">Noise whitening in LTPDA</a></h2>
90 <p>
91 A random process <i>w(t)</i> is considered white if it is zero mean
92 and uncorrelated:
93 </p>
94 <p>
95 <IMG src="images/whitening01.gif" align="center" border="0">
96 </p>
97 <p>
98 As a consequence, the power spectral density of a white process is a
99 constant at every frequency:
100 </p>
101 <p>
102 <IMG src="images/whitening02.gif" align="center" border="0">
103 </p>
104 <p>
105 In other words, The power per unit of frequency associated to a white noise
106 process is uniformly distributed on the whole available frequency range.
107 An example is reported in figure 1.
108 </p>
109
110 <div align="center">
111 <table border="0">
112 <caption align="bottom">
113 <b> Figure 1:</b> Power spectral density (estimated with the welch method)
114 of a gaussian unitary variance zero mean random process.
115 The process <i>w(t)</i> is assumed to have
116 physical units of <tt>m</tt> therefore its power spectral density has
117 physical units of <tt>m^2/Hz</tt>. Note that the power spectral density
118 average value is 2 instead of the expected 1 (unitary variance process)
119 since we calculated one-sided power spectral density.
120 </caption>
121 <tr>
122 <td>
123 <IMG src="images/whitening03.png" align="center" border="0">
124 </td>
125 </tr>
126 </table>
127 </div>
128
129 <p>
130 </p>
131
132 <p>
133 A non-white (colored) noise process is instead characterized by a given
134 distribution of the power per unit of frequency along the available frequency
135 bandwidth. <br>
136 Whitening operation on a given non-white process corresponds to force
137 such a process to satisfy the conditions described above for a white process.
138 </p>
139 <p>
140 In LTPDA there are different methods for noise whitening:
141 <ul>
142 <li> <a href="matlab:doc('ao/buildWhitener1D')"> buildWhitener1D.m</a>
143 <li> <a href="matlab:doc('ao/whiten1D')"> whiten1D.m</a>
144 <li> <a href="matlab:doc('ao/firwhiten')">firwhiten.m</a>
145 <li> <a href="matlab:doc('ao/whiten2D')">whiten2D.m</a>
146 </ul>
147 They accept time series analysis objects as an input and they output noise
148 whitening filters or whitened time series analysis objects.
149 </p>
150
151 <!-- ===== Algorithm ====== -->
152 <h2><a name="WhitenAlgo">Whitening Algorithms</a></h2>
153
154 <h3>buildWhitener1D</h3>
155 <p>
156 <tt>buildWhitener1D</tt> performs a frequency domain identification of the system
157 in order to extract the proper whitening filter. The function needs a model
158 for the one-sided power spectral density of the given process. If no model
159 is provided, the power spectral density of the process is calculated with
160 the <a href="matlab:doc('ao/psd')">psd</a> and <a href="matlab:doc('ao/bin_data')">bin_data</a> algorithm. <br>
161 <ol>
162 <li> The inverse of the square root of the model for the power spectral
163 density is fit in z-domain in order to determine a whitening
164 filter.
165 <li> Unstable poles are removed by an all-pass stabilization procedure.
166 <li> Whitening filter is provided at the output.
167 </ol>
168 </p>
169
170 <h3>Whiten1D</h3>
171 <p>
172 <tt>whiten1D</tt> implements the same functionality of <tt>buildWhitener1D</tt>
173 but it adds the filtering step so input data are filtered with the identified filter
174 internally to the method.
175 </p>
176
177 <h3>Firwhiten</h3>
178 <p>
179 <tt>firwhiten</tt> whitens the input time-series by building an FIR
180 whitening filter. <br>
181 <ol>
182 <li> Make ASD of time-series.
183 <li> Perform running median to get noise-floor estimate <a href="matlab:doc('ao/smoother')">ao/smoother</a>.
184 <li> Invert noise-floor estimate.
185 <li> Call <a href="matlab:doc('mfir')">mfir()</a> on noise-floor estimate to produce whitening filter.
186 <li> Filter data.
187 </ol>
188 </p>
189
190 <h3>Whiten2D</h3>
191 <p>
192 <tt>whiten2D</tt> whitens cross-correlated time-series. Whitening
193 filters are constructed by a fitting procedure to the models
194 for the corss-spectral matrix provided.
195 In order to work with <tt>whiten2D</tt> you must provide
196 a model (frequency series analysis objects) for the cross-spectral density
197 matrix of the process.
198 <ol>
199 <li> Whitening filters frequency response is calculated by the
200 eigendecomposition of the cross-spectral matrix.
201 <li> Calculated responses are fit in z-domain in order to identify
202 corresponding autoregressive moving average filters.
203 <li> Input time-series is filtered. The filtering process corresponds to:<br>
204 w(1) = Filt11(a(1)) + Filt12(a(2))<br>
205 w(2) = Filt21(a(1)) + Filt22(a(2))
206 </ol>
207 </p>
208
209
210 <!-- ===== 1D Examples ====== -->
211 <h2><a name="buildWhitener1D">Whitening noise in one-dimensional data</a></h2>
212 <p>
213 We can now test an example of the one-dimensinal whitening filters capabilities.
214 With the following commands we can generate a colored noise data series
215 for parameters description please refer to the
216 <a href="matlab:doc('ao')">ao</a>,
217 <a href="matlab:doc('miir')">miir</a> and
218 <a href="matlab:doc('ao/filter')">filter</a>
219 documentation pages.
220 </p>
221 <div class="fragment"><pre>
222
223 fs = 1; <span class="comment">% sampling frequency</span>
224
225 <span class="comment">% Generate gaussian white noise</span>
226 pl = plist(<span class="string">'tsfcn'</span>, <span class="string">'randn(size(t))'</span>, ...
227 <span class="string">'fs'</span>, fs, ...
228 <span class="string">'nsecs'</span>, 1e5, ...
229 <span class="string">'yunits'</span>, <span class="string">'m'</span>);
230 a = ao(pl);
231
232 <span class="comment">% Get a coloring filter</span>
233 pl = plist(<span class="string">'type'</span>, <span class="string">'bandpass'</span>, ...
234 <span class="string">'fs'</span>, fs, ...
235 <span class="string">'order'</span>, 3, ...
236 <span class="string">'gain'</span>, 1, ...
237 <span class="string">'fc'</span>, [0.03 0.1]);
238 ft = miir(pl);
239
240 <span class="comment">% Coloring noise</span>
241 af = filter(a, ft);
242
243 </pre></div>
244
245 <p>
246 Now we can try to white colored noise.
247 </p>
248
249
250 <h3>buildWhitener1D</h3>
251 <p>
252 If you want to try <tt>buildWhitener1D</tt> to get a whitening filter for
253 the present colored noise, you can try the following code. Please refer to the
254 <a href="matlab:doc('ao/buildWhitener1D')">buildWhitener1D</a> documentation page
255 for the meaning of any parameter. The result of the whitening procedure
256 is reported in figure 2.
257 </p>
258
259 <div class="fragment"><pre>
260
261 pl = plist(...
262 <span class="string">'MaxIter'</span>, 30, ...
263 <span class="string">'MinOrder'</span>, 9, ...
264 <span class="string">'MaxOrder'</span>, 15, ...
265 <span class="string">'FITTOL'</span>, 5e-2);
266
267 wfil = buildWhitener1D(af,pl);
268
269 aw = filter(af,wfil);
270
271 </pre></div>
272
273 <div align="center">
274 <table border="0">
275 <caption align="bottom">
276 <b> Figure 2:</b> Power spectral density (estimated with the welch method)
277 of colored and whitened processes.
278 </caption>
279 <tr>
280 <td>
281 <IMG src="images/whitening04.png" align="center" border="0">
282 </td>
283 </tr>
284 </table>
285 </div>
286
287
288 <h3>Firwhiten</h3>
289 <p>
290 As an alternative you can try <tt>firwhiten</tt> to whiten the present
291 colored noise. Please refer to the
292 <a href="matlab:doc('ao/firwhiten')">firwhiten</a> documentation page
293 for the meaning of any parameter. The result of the whitening procedure
294 is reported in figure 3.
295 </p>
296
297 <div class="fragment"><pre>
298
299 pl = plist(...
300 <span class="string">'Ntaps'</span>, 5000, ...
301 <span class="string">'Nfft'</span>, 1e5, ...
302 <span class="string">'BW'</span>, 5);
303
304 aw = firwhiten(af, pl);
305
306 </pre></div>
307
308 <div align="center">
309 <table border="0">
310 <caption align="bottom">
311 <b> Figure 3:</b> Power spectral density (estimated with the welch method)
312 of colored and whitened processes.
313 </caption>
314 <tr>
315 <td>
316 <IMG src="images/whitening05.png" align="center" border="0">
317 </td>
318 </tr>
319 </table>
320 </div>
321
322
323 <!-- ===== 2D Examples ====== -->
324 <h2><a name="Whiten2D">Whitening noise in two-dimensional data</a></h2>
325
326 <p>
327 We consider now the problem of whitening cross correlated data series.
328 As a example we consider a typical couple of x-dynamics LTP data series.
329 <tt>a1</tt> and <tt>a2</tt> are interferometer output noise data series.
330 In oreder to whiten data we must input a frequency response model of the
331 cross spectral matrix of the cross-correlated process.
332 </p>
333 <p>
334 <IMG src="images/whitening10.gif" align="center" border="0">
335 </p>
336 <p>
337 Refer to <a href="matlab:doc('ao/firwhiten')">whiten2D</a> documentation page
338 for the meaning of any parameter.
339 </p>
340
341 <div class="fragment"><pre>
342
343 pl = plist(...
344 <span class="string">'csd11'</span>, mod11, ...
345 <span class="string">'csd12'</span>, mod12, ...
346 <span class="string">'csd21'</span>, mod21, ...
347 <span class="string">'csd22'</span>, mod22, ...
348 <span class="string">'MaxIter'</span>, 75, ...
349 <span class="string">'PoleType'</span>, 3, ...
350 <span class="string">'MinOrder'</span>, 20, ...
351 <span class="string">'MaxOrder'</span>, 40, ...
352 <span class="string">'Weights'</span>, 2, ...
353 <span class="string">'Plot'</span>, false,...
354 <span class="string">'Disp'</span>, false,...
355 <span class="string">'MSEVARTOL'</span>, 1e-2,...
356 <span class="string">'FITTOL'</span>, 1e-3);
357
358 [aw1,aw2] = whiten2D(a1,a2,pl);
359
360 </pre></div>
361
362 <div align="center">
363 <table border="0">
364 <caption align="bottom">
365 <b> Figure 4:</b> Power spectral density of the noisy data series
366 before (left) and after (right) the whitening.
367 </caption>
368 <tr>
369 <td>
370 <IMG src="images/whitening06.png" align="center" border="0">
371 </td>
372 <td>
373 <IMG src="images/whitening07.png" align="center" border="0">
374 </td>
375 </tr>
376 </table>
377 </div>
378
379 <div align="center">
380 <table border="0">
381 <caption align="bottom">
382 <b> Figure 5:</b> Real (left) and Imaginary (right) part of the
383 <a href="matlab:doc('ao/cohere')">coherence</a> function.
384 Blue line refers to theoretical expectation for colored noise data.
385 Red line refers to calculated values for colored noise data.
386 Green line refers to calculated values for whitened noise data.
387
388 </caption>
389 <tr>
390 <td>
391 <IMG src="images/whitening08.png" align="center" border="0">
392 </td>
393 <td>
394 <IMG src="images/whitening09.png" align="center" border="0">
395 </td>
396 </tr>
397 </table>
398 </div>
399
400
401
402
403
404
405
406
407
408 </p>
409
410 <br>
411 <br>
412 <table class="nav" summary="Navigation aid" border="0" width=
413 "100%" cellpadding="0" cellspacing="0">
414 <tr valign="top">
415 <td align="left" width="20"><a href="gapfill.html"><img src=
416 "b_prev.gif" border="0" align="bottom" alt=
417 "Data gap filling"></a>&nbsp;</td>
418
419 <td align="left">Data gap filling</td>
420
421 <td>&nbsp;</td>
422
423 <td align="right">Signal Processing in LTPDA</td>
424
425 <td align="right" width="20"><a href=
426 "sigproc.html"><img src="b_next.gif" border="0" align=
427 "bottom" alt="Signal Processing in LTPDA"></a></td>
428 </tr>
429 </table><br>
430
431 <p class="copy">&copy;LTP Team</p>
432 </body>
433 </html>