Mercurial > hg > ltpda

comparison m-toolbox/html_help/help/ug/ltpda_training_topic_5_1.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>System identification in z-domain (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.html"><img src="b_prev.gif" border="0" align=
30 "bottom" alt="Topic 5 - Model fitting"></a>&nbsp;&nbsp;&nbsp;<a href=
31 "ltpda_training_topic_5_2.html"><img src="b_next.gif" border="0" align=
32 "bottom" alt="Generation of noise with given PSD"></a></td>
33 </tr>
34 </table>
35
36 <h1 class="title"><a name="f3-12899" id="f3-12899"></a>System identification in z-domain</h1>
37 <hr>
38
39 <p>
40 <p>
41 System identification in Z-domain is performed with the function <tt>ao/zDomainFit</tt>.
42 It is based on a modified version of the vector fitting algorithm that was
43 adapted to fit in the Z-domain. Details of the agorithm can be found in the
44 <a href="zdomainfit.html">Z-domain fit documentation page</a></li>.
45 </p>
46
47 <h2> System identification in Z-domain </h2>
48
49 <p>
50 During this exercise we will:
51 <ol>
52 <li> Generate white noise
53 <li> Filter white noise with a <tt>miir</tt> filter generated by a <tt> pzmodel</tt>
54 <li> Extract the transfer function from data
55 <li> Fit the transfer function with <tt>ao/zDomainFit</tt>
56 <li> Check results
57 </ol>
58 </p>
59
60 <p>
61 Let's start by generating some white noise.
62 </p>
63
64 <div class="fragment"><pre>
65 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>));
66 </pre></div>
67 <p>
68 This command generates a time series of gaussian distributed random noise
69 with a sampling frequency (<span class="string">'fs'</span>) of 1 Hz,
70 10000 seconds long (<span class="string">'nsecs'</span>) and with
71 <span class="string">'yunits'</span> set to meters (<span class="string">'m'</span>).
72 </p>
73 <p>
74 Now we are ready to move on the second step where we will:
75 <ul>
76 <li> Build a pole-zero model (<tt>pzmodel</tt>)
77 <li> Construct a <tt>miir</tt> filter from the <tt>pzmodel</tt>
78 <li> filter white noise data with the <tt>miir</tt> filter in order to obtain a
79 colored noise time series.
80 </ul>
81
82 </p>
83 <div class="fragment"><pre>
84 pzm = pzmodel(1, [0.005 2], [0.05 4]);
85
86 filt = miir(pzm, plist(<span class="string">'fs'</span>, 1, <span class="string">'name'</span>, <span class="string">'None'</span>));
87 filt.setIunits(<span class="string">'m'</span>);
88 filt.setOunits(<span class="string">'V'</span>);
89
90 <span class="comment">% Filter the data</span>
91 ac = filter(a,filt);
92 ac.simplifyYunits;
93 </pre></div>
94 <p>
95 We can calculate the PSD of the data streams in order to check the difference between the coloured noise and the
96 white noise. Let's choose the log-scale estimation method.
97 </p>
98 <div class="fragment"><pre>
99 axx = lpsd(a);
100 acxx = lpsd(ac);
101 iplot(axx,acxx)
102 </pre></div>
103 <p>
104 You should obtain a plot similar to this:
105 <div align="center">
106 <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_1_1.png" align="center" border="0">
107 </div>
108 </p>
109 <p>
110 Let us move to the third step. We will generate the transfer function from the data
111 and split it in order to remove the first 3 bins. The last operation is useful
112 for the fitting process.
113 </p>
114 <div class="fragment"><pre>
115 tf = ltfe(a,ac);
116 tfsp = split(tf,plist(<span class="string">'frequencies'</span>, [5e-4 5e-1]));
117
118 iplot(tf,tfsp)
119 </pre></div>
120 <p>
121 The plot should look like the following:
122 <div align="center">
123 <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_1_2.png" align="center" border="0">
124 </div>
125 </p>
126 <p>
127 It is now the moment to start fitting with <tt>zDomainFit</tt>. As reported in the
128 function help page we can run an automatic search loop to identify proper model order.
129 In such a case we have to define a set of conditions to check fit accuracy and to exit the fitting loop.<br/>
130 We can start checking Mean Squared Error and variation (CONDTYPE = 'MSE').
131 It checks if the normalized mean squared error is lower than the value specified in the parameter FITTOL
132 and if the relative variation of the mean squared error is lower than the value specified in the parameter MSEVARTOL.
133 <p><!-- Table of the set: Default -->
134 <a name="1"/>
135 <table cellspacing="0" class="body" cellpadding="4" summary="" width="100%" border="2">
136 <colgroup>
137 <col width="15%"/>
138 <col width="20%"/>
139 <col width="20%"/>
140 <col width="45%"/>
141 </colgroup>
142 <thead>
143 <tr valign="top">
144 <th bgcolor="#B9C6DD" colspan="4"><h3>Default</h3></th>
145 </tr>
146 <tr valign="top">
147 <th bgcolor="#D7D7D7">Key</th>
148 <th bgcolor="#D7D7D7">Default Value</th>
149 <th bgcolor="#D7D7D7">Options</th>
150 <th bgcolor="#D7D7D7">Description</th>
151 </tr>
152 </thead>
153 <tbody>
154 <tr valign="top">
155 <td bgcolor="#F2F2F2">CONDTYPE</td>
156 <td bgcolor="#F2F2F2">'MSE'</td>
157 <td bgcolor="#F2F2F2">
158 <ul>
159 <li><font color="#1111FF">'MSE'</font></li>
160 <li><font color="#1111FF">'RLD'</font></li>
161 <li><font color="#1111FF">'RSF'</font></li>
162 </ul>
163 </td>
164 <td bgcolor="#F2F2F2">Fit conditioning type. Admitted values are:
165 <ul>
166 <li>'MSE' Mean Squared Error and variation</li>
167 <li>'RLD' Log residuals difference and mean squared error variation</li>
168 <li>'RSF' Residuals spectral flatness and mean squared error variation</li>
169 </ul>
170 </td>
171 </tr>
172 <tr valign="top">
173 <td bgcolor="#F2F2F2">FITTOL</td>
174 <td bgcolor="#F2F2F2">0.001</td>
175 <td bgcolor="#F2F2F2"><i>none</i></td>
176 <td bgcolor="#F2F2F2">Fit tolerance.</td>
177 </tr>
178 <tr valign="top">
179 <td bgcolor="#F2F2F2">MSEVARTOL</td>
180 <td bgcolor="#F2F2F2">0.01</td>
181 <td bgcolor="#F2F2F2"><i>none</i></td>
182 <td bgcolor="#F2F2F2">Mean Squared Error Variation - Check if the<br/>relative variation of the mean squared error is<br/>smaller than the value specified. This<br/>option is useful for finding the minimum of the Chi-squared.</td>
183 </tr>
184 </tbody>
185 </table>
186 </p>
187 <p>You will find a list of all Parameters of zDomainFit here: <a href="matlab:web(ao.getInfo('zDomainFit').tohtml, '-helpbrowser')">Parameter of zDomainFit</a></p>
188
189 Now let's run the fit:
190 </p>
191 <div class="fragment"><pre>
192 <span class="comment">% Set up the parameters</span>
193 plfit = plist(<span class="string">'fs'</span>,1,... <span class="comment">% Sampling frequency for the model filters</span>
194 <span class="string">'AutoSearch'</span>,<span class="string">'on'</span>,... <span class="comment">% Automatically search for a good model</span>
195 <span class="string">'StartPolesOpt'</span>,<span class="string">'clog'</span>,... <span class="comment">% Define the properties of the starting poles - complex distributed in the unitary circle</span>
196 <span class="string">'maxiter'</span>,50,... <span class="comment">% Maximum number of iteration per model order</span>
197 <span class="string">'minorder'</span>,2,... <span class="comment">% Minimum model order</span>
198 <span class="string">'maxorder'</span>,9,... <span class="comment">% Maximum model order</span>
199 <span class="string">'weightparam'</span>,<span class="string">'abs'</span>,... <span class="comment">% Assign weights as 1./abs(data)</span>
200 <span class="string">'condtype'</span>,<span class="string">'MSE'</span>,... <span class="comment">% Mean Squared Error and variation</span>
201 <span class="string">'fittol'</span>,1e-2,... <span class="comment">% Fit tolerance</span>
202 <span class="string">'msevartol'</span>,1e-1,... <span class="comment">% Mean Squared Error Variation tolerance</span>
203 <span class="string">'Plot'</span>,<span class="string">'on'</span>,... <span class="comment">% Set the plot on or off</span>
204 <span class="string">'ForceStability'</span>,<span class="string">'on'</span>,... <span class="comment">% Force to output a stable poles model</span>
205 <span class="string">'CheckProgress'</span>,<span class="string">'off'</span>); <span class="comment">% Display fitting progress on the command window</span>
206
207 <span class="comment">% Do the fit</span>
208 fobj = zDomainFit(tfsp,plfit);
209
210 <span class="comment">% Set the input and output units for fitted model</span>
211 fobj.setIunits(<span class="string">'m'</span>);
212 fobj.setOunits(<span class="string">'V'</span>);
213 </pre></div>
214 <p>
215 When <span class="string">'Plot'</span> parameter is set to <span class="string">'on'</span>
216 the function plots the fit progress.
217 <div align="center">
218 <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_1_3.png" align="center" border="0">
219 </div>
220 </p>
221
222 <p>
223 We can now check the result of our fitting procedures.
224 We calculate the frequency response of the fitted models (filters) and compare them
225 with the starting IIR filter response, then we will plot the percentage error on the filters magnitudes.
226 </p>
227 <p>
228 Note that the result of the fitting procedure is a <tt>matrix</tt> object containing a <tt>filterbank</tt> object,
229 which itself contains a parallel bank of 3 IIR filters.
230 </p>
231 <p>
232 Note that at the moment we need to access the individual filter objects inside the <tt>matrix</tt> object that was the result of the fitting procedure.</p>
233 <div class="fragment"><pre>
234 <span class="comment">% set plist for filter response</span>
235 plrsp = plist(<span class="string">'bank'</span>,<span class="string">'parallel'</span>,<span class="string">'f1'</span>,1e-5,<span class="string">'f2'</span>,0.5,<span class="string">'nf'</span>,100,<span class="string">'scale'</span>,<span class="string">'log'</span>);
236
237 <span class="comment">% compute the response of the original noise-shape filter</span>
238 rfilt = resp(filt,plrsp);
239 rfilt.setName;
240
241 <span class="comment">% compute the response of our fitted filter bank</span>
242 rfobj = resp(fobj.filters,plrsp);
243 rfobj.setName;
244
245 <span class="comment">% compare the responses</span>
246 iplot(rfilt,rfobj)
247
248 <span class="comment">% and the percentage error on the magnitude</span>
249 pdiff = 100.*abs((rfobj-rfilt)./rfilt);
250 pdiff.simplifyYunits;
251 iplot(pdiff,plist(<span class="string">'YRanges'</span>,[1e-2 100]))
252 </pre></div>
253 <p>
254 The first plot shows the response of the original filter and the fitted filter bank,
255 whereas the second plot reports the percentage difference between fitted model and target filter magnitude. <br/>
256 As can be seen, the difference between filters magnitude is at most 10%.
257 <div align="center">
258 <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_1_4.png" align="center" border="0">
259 </div>
260 <div align="center">
261 <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_1_5.png" align="center" border="0">
262 </div>
263 </p>
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282 </p>
283
284 <br>
285 <br>
286 <table class="nav" summary="Navigation aid" border="0" width=
287 "100%" cellpadding="0" cellspacing="0">
288 <tr valign="top">
289 <td align="left" width="20"><a href="ltpda_training_topic_5.html"><img src=
290 "b_prev.gif" border="0" align="bottom" alt=
291 "Topic 5 - Model fitting"></a>&nbsp;</td>
292
293 <td align="left">Topic 5 - Model fitting</td>
294
295 <td>&nbsp;</td>
296
297 <td align="right">Generation of noise with given PSD</td>
298
299 <td align="right" width="20"><a href=
300 "ltpda_training_topic_5_2.html"><img src="b_next.gif" border="0" align=
301 "bottom" alt="Generation of noise with given PSD"></a></td>
302 </tr>
303 </table><br>
304
305 <p class="copy">&copy;LTP Team</p>
306 </body>
307 </html>