Mercurial > hg > ltpda

<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
   "http://www.w3.org/TR/1999/REC-html401-19991224/loose.dtd">

<html lang="en">
<head>
  <meta name="generator" content=
  "HTML Tidy for Mac OS X (vers 1st December 2004), see www.w3.org">
  <meta http-equiv="Content-Type" content=
  "text/html; charset=us-ascii">

  <title>System identification in z-domain (LTPDA Toolbox)</title>
  <link rel="stylesheet" href="docstyle.css" type="text/css">
  <meta name="generator" content="DocBook XSL Stylesheets V1.52.2">
  <meta name="description" content=
  "Presents an overview of the features, system requirements, and starting the toolbox.">
  </head>

<body>
  <a name="top_of_page" id="top_of_page"></a>

  <p style="font-size:1px;">&nbsp;</p>

  <table class="nav" summary="Navigation aid" border="0" width=
  "100%" cellpadding="0" cellspacing="0">
    <tr>
      <td valign="baseline"><b>LTPDA Toolbox</b></td><td><a href="../helptoc.html">contents</a></td>

      <td valign="baseline" align="right"><a href=
      "ltpda_training_topic_5.html"><img src="b_prev.gif" border="0" align=
      "bottom" alt="Topic 5 - Model fitting"></a>&nbsp;&nbsp;&nbsp;<a href=
      "ltpda_training_topic_5_2.html"><img src="b_next.gif" border="0" align=
      "bottom" alt="Generation of noise with given PSD"></a></td>
    </tr>
  </table>

  <h1 class="title"><a name="f3-12899" id="f3-12899"></a>System identification in z-domain</h1>
  <hr>

  <p>
	<p>
  System identification in Z-domain is performed with the function <tt>ao/zDomainFit</tt>.
  It is based on a modified version of the vector fitting algorithm that was
  adapted to fit in the Z-domain. Details of the agorithm can be found in the
  <a href="zdomainfit.html">Z-domain fit documentation page</a></li>.
</p>

<h2> System identification in Z-domain </h2>

<p>
  During this exercise we will:
  <ol>
    <li> Generate white noise
    <li> Filter white noise with a <tt>miir</tt> filter generated by a <tt> pzmodel</tt>
    <li> Extract the transfer function from data
    <li> Fit the transfer function with <tt>ao/zDomainFit</tt>
    <li> Check results
  </ol>
</p>

<p>
  Let's start by generating some white noise.
</p>

<div class="fragment"><pre>
    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>));
</pre></div>
<p>
  This command generates a time series of gaussian distributed random noise
  with a sampling frequency (<span class="string">'fs'</span>) of 1 Hz,
  10000 seconds long (<span class="string">'nsecs'</span>) and with
  <span class="string">'yunits'</span> set to meters (<span class="string">'m'</span>).
</p>
<p>
  Now we are ready to move on the second step where we will:
  <ul>
    <li> Build a pole-zero model (<tt>pzmodel</tt>)
    <li> Construct a <tt>miir</tt> filter from the <tt>pzmodel</tt>
    <li> filter white noise data with the <tt>miir</tt> filter in order to obtain a
    colored noise time series.
  </ul>

</p>
<div class="fragment"><pre>
    pzm = pzmodel(1, [0.005 2], [0.05 4]);

    filt   = miir(pzm, plist(<span class="string">'fs'</span>, 1, <span class="string">'name'</span>, <span class="string">'None'</span>));
    filt.setIunits(<span class="string">'m'</span>);
    filt.setOunits(<span class="string">'V'</span>);

    <span class="comment">% Filter the data</span>
    ac = filter(a,filt);
    ac.simplifyYunits;
</pre></div>
<p>
  We can calculate the PSD of the data streams in order to check the difference between the coloured noise and the
  white noise. Let's choose the log-scale estimation method.
</p>
<div class="fragment"><pre>
    axx  = lpsd(a);
    acxx = lpsd(ac);
    iplot(axx,acxx)
</pre></div>
<p>
  You should obtain a plot similar to this:
  <div align="center">
    <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_1_1.png"  align="center" border="0">
  </div>
</p>
<p>
  Let us move to the third step. We will generate the transfer function from the data
  and split it in order to remove the first 3 bins. The last operation is useful
  for the fitting process.
</p>
<div class="fragment"><pre>
    tf = ltfe(a,ac);
    tfsp = split(tf,plist(<span class="string">'frequencies'</span>, [5e-4 5e-1]));

    iplot(tf,tfsp)
</pre></div>
<p>
  The plot should look like the following:
  <div align="center">
    <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_1_2.png"  align="center" border="0">
  </div>
</p>
<p>
  It is now the moment to start fitting with <tt>zDomainFit</tt>. As reported in the
  function help page we can run an automatic search loop to identify proper model order.
  In such a case we have to define a set of conditions to check fit accuracy and to exit the fitting loop.<br/>
  We can start checking Mean Squared Error and variation (CONDTYPE = 'MSE').
  It checks if the normalized mean squared error is lower than the value specified in the parameter FITTOL
  and if the relative variation of the mean squared error is lower than the value specified in the parameter MSEVARTOL.
  <p><!-- Table of the set: Default -->
    <a name="1"/>
    <table cellspacing="0" class="body" cellpadding="4" summary="" width="100%" border="2">
      <colgroup>
        <col width="15%"/>
        <col width="20%"/>
        <col width="20%"/>
        <col width="45%"/>
      </colgroup>
      <thead>
        <tr valign="top">
          <th bgcolor="#B9C6DD" colspan="4"><h3>Default</h3></th>
        </tr>
        <tr valign="top">
          <th bgcolor="#D7D7D7">Key</th>
          <th bgcolor="#D7D7D7">Default Value</th>
          <th bgcolor="#D7D7D7">Options</th>
          <th bgcolor="#D7D7D7">Description</th>
        </tr>
      </thead>
      <tbody>
        <tr valign="top">
          <td bgcolor="#F2F2F2">CONDTYPE</td>
          <td bgcolor="#F2F2F2">'MSE'</td>
          <td bgcolor="#F2F2F2">
            <ul>
              <li><font color="#1111FF">'MSE'</font></li>
              <li><font color="#1111FF">'RLD'</font></li>
              <li><font color="#1111FF">'RSF'</font></li>
            </ul>
          </td>
          <td bgcolor="#F2F2F2">Fit conditioning type. Admitted values are:
            <ul>
              <li>'MSE' Mean Squared Error and variation</li>
              <li>'RLD' Log residuals difference and mean squared error variation</li>
              <li>'RSF' Residuals spectral flatness and mean squared error variation</li>
            </ul>
          </td>
        </tr>
        <tr valign="top">
          <td bgcolor="#F2F2F2">FITTOL</td>
          <td bgcolor="#F2F2F2">0.001</td>
          <td bgcolor="#F2F2F2"><i>none</i></td>
          <td bgcolor="#F2F2F2">Fit tolerance.</td>
        </tr>
        <tr valign="top">
          <td bgcolor="#F2F2F2">MSEVARTOL</td>
          <td bgcolor="#F2F2F2">0.01</td>
          <td bgcolor="#F2F2F2"><i>none</i></td>
          <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>
        </tr>
      </tbody>
    </table>
  </p>
  <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>

  Now let's run the fit:
</p>
<div class="fragment"><pre>
    <span class="comment">% Set up the parameters</span>
    plfit = plist(<span class="string">'fs'</span>,1,...      <span class="comment">% Sampling frequency for the model filters</span>
      <span class="string">'AutoSearch'</span>,<span class="string">'on'</span>,...       <span class="comment">% Automatically search for a good model</span>
      <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>
      <span class="string">'maxiter'</span>,50,...            <span class="comment">% Maximum number of iteration per model order</span>
      <span class="string">'minorder'</span>,2,...            <span class="comment">% Minimum model order</span>
      <span class="string">'maxorder'</span>,9,...            <span class="comment">% Maximum model order</span>
      <span class="string">'weightparam'</span>,<span class="string">'abs'</span>,...     <span class="comment">% Assign weights as 1./abs(data)</span>
      <span class="string">'condtype'</span>,<span class="string">'MSE'</span>,...        <span class="comment">% Mean Squared Error and variation</span>
      <span class="string">'fittol'</span>,1e-2,...           <span class="comment">% Fit tolerance</span>
      <span class="string">'msevartol'</span>,1e-1,...        <span class="comment">% Mean Squared Error Variation tolerance</span>
      <span class="string">'Plot'</span>,<span class="string">'on'</span>,...             <span class="comment">% Set the plot on or off</span>
      <span class="string">'ForceStability'</span>,<span class="string">'on'</span>,...   <span class="comment">% Force to output a stable poles model</span>
      <span class="string">'CheckProgress'</span>,<span class="string">'off'</span>);     <span class="comment">% Display fitting progress on the command window</span>

    <span class="comment">% Do the fit</span>
    fobj = zDomainFit(tfsp,plfit);

    <span class="comment">% Set the input and output units for fitted model</span>
    fobj.setIunits(<span class="string">'m'</span>);
    fobj.setOunits(<span class="string">'V'</span>);
</pre></div>
<p>
  When <span class="string">'Plot'</span> parameter is set to <span class="string">'on'</span>
  the function plots the fit progress.
  <div align="center">
    <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_1_3.png"  align="center" border="0">
  </div>
</p>

  <p>
    We can now check the result of our fitting procedures.
    We calculate the frequency response of the fitted models (filters) and compare them
    with the starting IIR filter response, then we will plot the percentage error on the filters magnitudes.
  </p>
  <p>
    Note that the result of the fitting procedure is a <tt>matrix</tt> object containing a <tt>filterbank</tt> object,
  which itself contains a parallel bank of 3 IIR filters.
</p>
    <p>
    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>
<div class="fragment"><pre>
    <span class="comment">% set plist for filter response</span>
    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>);

    <span class="comment">% compute the response of the original noise-shape filter</span>
    rfilt = resp(filt,plrsp);
    rfilt.setName;

    <span class="comment">% compute the response of our fitted filter bank</span>
    rfobj = resp(fobj.filters,plrsp);
    rfobj.setName;

    <span class="comment">% compare the responses</span>
    iplot(rfilt,rfobj)

    <span class="comment">% and the percentage error on the magnitude</span>
    pdiff = 100.*abs((rfobj-rfilt)./rfilt);
    pdiff.simplifyYunits;
    iplot(pdiff,plist(<span class="string">'YRanges'</span>,[1e-2 100]))
</pre></div>
<p>
  The first plot shows the response of the original filter and the fitted filter bank,
  whereas the second plot reports the percentage difference between fitted model and target filter magnitude. <br/>
  As can be seen, the difference between filters magnitude is at most 10%.
  <div align="center">
    <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_1_4.png" align="center" border="0">
  </div>
  <div align="center">
    <IMG src="images/ltpda_training_1/topic5/ltpda_training_5_1_5.png" align="center" border="0">
  </div>
</p>


  </p>

  <br>
  <br>
  <table class="nav" summary="Navigation aid" border="0" width=
  "100%" cellpadding="0" cellspacing="0">
    <tr valign="top">
      <td align="left" width="20"><a href="ltpda_training_topic_5.html"><img src=
      "b_prev.gif" border="0" align="bottom" alt=
      "Topic 5 - Model fitting"></a>&nbsp;</td>

      <td align="left">Topic 5 - Model fitting</td>

      <td>&nbsp;</td>

      <td align="right">Generation of noise with given PSD</td>

      <td align="right" width="20"><a href=
      "ltpda_training_topic_5_2.html"><img src="b_next.gif" border="0" align=
      "bottom" alt="Generation of noise with given PSD"></a></td>
    </tr>
  </table><br>

  <p class="copy">&copy;LTP Team</p>
</body>
</html>
author	Daniele Nicolodi <nicolodi@science.unitn.it>
date	Mon, 05 Dec 2011 16:20:06 +0100
parents	f0afece42f48
children