Mercurial > hg > ltpda

comparison m-toolbox/html_help/help/ug/ltpda_training_topic_5_2.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>Generation of noise with given PSD (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 "ltpda_training_topic_5_1.html"><img src="b_prev.gif" border="0" align=
30 "bottom" alt="System identification in z-domain"></a>&nbsp;&nbsp;&nbsp;<a href=
31 "ltpda_training_topic_5_3.html"><img src="b_next.gif" border="0" align=
32 "bottom" alt="Fitting time series with polynomials"></a></td>
33 </tr>
34 </table>
35
36 <h1 class="title"><a name="f3-12899" id="f3-12899"></a>Generation of noise with given PSD</h1>
37 <hr>
38
39 <p>
40
41
42
43 <p>
44 Generation of model noise is performed with the function <tt>ao/noisegen1D</tt>.
45 Details on the algorithm can be found in <a href="ng1D.html">noisegen1D help page</a>.
46 </p>
47
48 <h2> Generation of noise with given PSD</h2>
49
50 <p>
51 During this exercise we will:
52 <ol>
53 <li> Load from file an fsdata object with the model (obtained with a fit to the the PSD of test data with zDomainFit)
54 <li> Genarate noise from this model
55 <li> Compare PSD of the generated noise with original PSD
56 </ol>
57 </p>
58
59 <p>
60 Let's open a new editor window and load the test data.
61 </p>
62
63 <div class="fragment"><pre>
64 tn = ao(plist(<span class="string">'filename'</span>, <span class="string">'topic5/T5_Ex03_TestNoise.xml'</span>));
65 tn.setName;
66 </pre></div>
67
68 <p>
69 This command will load an Analysis Object containing a test time series
70 10000 seconds long, sampled at 1 Hz. The command <tt>setName</tt> sets the name
71 of the AO to be the same as the variable name, in this case <tt>tn</tt>.
72 </p>
73 <p>
74 Now let's calculate the PSD of our data. We apply some averaging, in order to decrease the
75 fluctuations in the data.
76 </p>
77
78 <div class="fragment"><pre>
79 tnxx = tn.psd(plist(<span class="string">'Nfft'</span>,2000));
80 </pre></div>
81
82 <p>
83 Additionally, we load a smooth model that represents well our data. It was obtained,
84 as described <a href="ltpda_training_topic_5_2.html#brbme84">here</a>,
85 by fitting the target PSD with z-domain fitting.
86 We load the data from disk, and plot them against the target PSD. Please note that
87 the colouring filters (whose response represents our model) have no units, so we force them to
88 be the same as the PSD we compare with:
89 <div class="fragment"><pre>
90 fmod = ao(plist(<span class="string">'filename'</span>, <span class="string">'topic5/T5_Ex03_ModelNoise.xml'</span>));
91 iplot(tnxx, fmod.setYunits(tnxx.yunits))
92 </pre></div>
93 </p>
94
95 <p>
96 The comparison beyween the target PSD and the model should look like:
97 <div align="center">
98 <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_2_4.png" align="center" border="0">
99 </div>
100 </p>
101
102 We can now start the noise generation process. The first step is to
103 generate a white time series Analysis Object, with the desired duration and units:
104 </p>
105 <div class="fragment"><pre>
106 a = ao(plist(<span class="string">'tsfcn'</span>,<span class="string">'randn(size(t))'</span>,<span class="string">'fs'</span>,1,<span class="string">'nsecs'</span>,10000,<span class="string">'yunits'</span>,<span class="string">'m'</span>));
107 </pre></div>
108 <p>
109 Then we run the noise coloring process calling <tt>noisegen1D</tt>
110 </p>
111
112 <div class="fragment"><pre>
113 plng = plist(...
114 <span class="string">'model'</span>, fmod, ... <span class="comment">% model for colored noise psd</span>
115 <span class="string">'MaxIter'</span>, 50, ... <span class="comment">% maximum number of fit iteration per model order</span>
116 <span class="string">'PoleType'</span>, 2, ... <span class="comment">% generates complex poles distributed in the unitary circle</span>
117 <span class="string">'MinOrder'</span>, 20, ... <span class="comment">% minimum model order</span>
118 <span class="string">'MaxOrder'</span>, 50, ... <span class="comment">% maximum model order</span>
119 <span class="string">'Weights'</span>, 2, ... <span class="comment">% weight with 1/abs(model)</span>
120 <span class="string">'Plot'</span>, false,... <span class="comment">% on to show the plot</span>
121 <span class="string">'Disp'</span>, false,... <span class="comment">% on to display fit progress on the command window</span>
122 <span class="string">'RMSEVar'</span>, 7,... <span class="comment">% Root Mean Squared Error Variation</span>
123 <span class="string">'FitTolerance'</span>, 2); <span class="comment">% Residuals log difference</span>
124
125 ac = noisegen1D(a, plng);
126 </pre></div>
127
128 <p>
129 Let's check the result. We calculate the PSD of the generated noise and compare it
130 with the PSD of the target data.
131 </p>
132
133 <div class="fragment"><pre>
134 acxx = ac.psd(plist(<span class="string">'Nfft'</span>,2000));
135 iplot(tnxx,acxx)
136 </pre></div>
137
138 <p>
139 As can be seen, the result is in quite satisfactory agreement with the original data
140 <div align="center">
141 <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_2_3.png" align="center" border="0">
142 </div>
143 </p>
144
145 <a name="brbme84"></a><h3 class="title" id="brbme84">Appendix: evaluation of the model for the noise PSD</h3>
146 <p>
147 The smooth model for the data, that we used to reproduce the synthesized noise, was actually obtained
148 by applying the procedure of z-domain fitting that we discussed in <a href="ltpda_training_topic_5_1.html">the previous section</a>.
149 If you want to practise more with this fitting technique, we repost here the steps.
150 </p>
151 <p>
152 In order to extract a reliable model from PSD data we need to discard the first frequency bins;
153 we do that by means of the <tt>split</tt> method.
154 </p>
155
156 <div class="fragment"><pre>
157 tnxxr = split(tnxx,plist(<span class="string">'frequencies'</span>, [2e-3 +inf]));
158 iplot(tnxx,tnxxr)
159 </pre></div>
160
161 <p>
162 The result should look like:
163 <div align="center">
164 <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_2_1.png"align="center" border="0">
165 </div>
166 </p>
167 <p>
168 Now it's the moment to fit our PSD to extract a smooth model to pass to the noise generator.
169
170 First of all we should define a set of proper weights for our fit process.
171 We smooth our PSD data and then define the weights as the inverse of the absolute value
172 of the smoothed PSD. This should help the fit function to do a good job with noisy data.
173 It is worth noting here that weights are not univocally defined and there could be better
174 ways to define them.
175 </p>
176 <div class="fragment"><pre>
177 stnxx = smoother(tnxxr);
178 iplot(tnxxr, stnxx)
179 wgh = 1./abs(stnxx);
180 </pre></div>
181 <p>
182 The result of the <tt>smoother</tt> method is shown in the plot below:
183 </p>
184 <div align="center">
185 <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_2_smoother.png" align="center" border="0">
186 </div>
187 <p>
188 Now let's run an automatic search for the proper model and pass the set
189 of externally defined weights. The first output of <tt>zDomainFit</tt> is a <tt>miir</tt> filter model;
190 the second output is the model response. Note that we are setting
191 <span class="string">'ResFlat'</span> parameter to define the exit condition.
192 <span class="string">'ResFlat'</span> check the spectral flatness of the
193 absolute value of the fit residuals.
194 </p>
195 <div class="fragment"><pre>
196 plfit = plist(<span class="string">'fs'</span>,1,...
197 <span class="string">'AutoSearch'</span>,<span class="string">'on'</span>,...
198 <span class="string">'StartPolesOpt'</span>,<span class="string">'clog'</span>,...
199 <span class="string">'maxiter'</span>,50,...
200 <span class="string">'minorder'</span>,30,...
201 <span class="string">'maxorder'</span>,45,...
202 <span class="string">'weights'</span>,wgh,... <span class="comment">% assign externally calculated weights</span>
203 <span class="string">'rmse'</span>,5,...
204 <span class="string">'condtype'</span>,<span class="string">'MSE'</span>,...
205 <span class="string">'msevartol'</span>,0.1,...
206 <span class="string">'fittol'</span>,0.01,...
207 <span class="string">'Plot'</span>,<span class="string">'on'</span>,...
208 <span class="string">'ForceStability'</span>,<span class="string">'off'</span>,...
209 <span class="string">'CheckProgress'</span>,<span class="string">'off'</span>);
210
211 <span class="comment">% Do the fit</span>
212 fit_results = zDomainFit(tnxxr,plfit);
213 </pre></div>
214
215 <p>
216 Fit result should look like:
217 <div align="center">
218 <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_2_2.png" align="center" border="0">
219 </div>
220 </p>
221
222 <p>
223 The fit results consist in a <tt>filterbank</tt> object; we can evaluate the absolute values of
224 the response of these filters at the frequencies defined by the <tt>x</tt> field of the PSD we want to match.
225
226 <div class="fragment"><pre>
227 <span class="comment">% Evaluate the absolute value of the response of the colouring filter</span>
228 b = resp(fit_results,plist(<span class="string">'f'</span>,tnxxr.x));
229 b.abs;
230
231 <span class="comment">% Save the model on disk</span>
232 b.save(plist(<span class="string">'filename'</span>, <span class="string">'topic5/T5_Ex03_ModelNoise.xml'</span>));
233 </pre></div>
234
235
236
237
238
239
240
241
242
243
244
245
246 </p>
247
248 <br>
249 <br>
250 <table class="nav" summary="Navigation aid" border="0" width=
251 "100%" cellpadding="0" cellspacing="0">
252 <tr valign="top">
253 <td align="left" width="20"><a href="ltpda_training_topic_5_1.html"><img src=
254 "b_prev.gif" border="0" align="bottom" alt=
255 "System identification in z-domain"></a>&nbsp;</td>
256
257 <td align="left">System identification in z-domain</td>
258
259 <td>&nbsp;</td>
260
261 <td align="right">Fitting time series with polynomials</td>
262
263 <td align="right" width="20"><a href=
264 "ltpda_training_topic_5_3.html"><img src="b_next.gif" border="0" align=
265 "bottom" alt="Fitting time series with polynomials"></a></td>
266 </tr>
267 </table><br>
268
269 <p class="copy">&copy;LTP Team</p>
270 </body>
271 </html>