Mercurial > hg > ltpda

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