<!-- $Id: zdomainfit_content.html,v 1.6 2009/08/27 11:38:58 luigi Exp $ -->

  <!-- ================================================== -->
  <!--                 BEGIN CONTENT FILE                 -->
  <!-- ================================================== -->
  <!-- ===== link box: Begin ===== -->
  <p>
    <table border="1"  width="80%">
      <tr>
        <td>
          <table border="0" cellpadding="5" class="categorylist" width="100%">
            <colgroup>
              <col width="37%"/>
              <col width="63%"/>
            </colgroup>
            <tbody>
              <tr valign="top">
                <td>
                  <a href="#description">Description</a>
                </td>
                <td>Z-domain system identification in LTPDA.</td>
              </tr>
              <tr valign="top">
                <td>
                  <a href="#algorithm">Algorithm</a>
                </td>
                <td>Fit Algorithm.</td>
              </tr>
              <tr valign="top">
                <td>
                  <a href="#examples">Examples</a>
                </td>
                <td>Usage example of z-domain system identification tool.</td>
              </tr>
              <tr valign="top">
                <td>
                  <a href="#references">References</a>
                </td>
                <td>Bibliographic references.</td>
              </tr>
            </tbody>
          </table>
        </td>
      </tr>
    </table>
  </p>
  <!-- ===== link box: End ====== -->



<h2><a name="description">Z-domain system identification in LTPDA</a></h2>
<p>
  System identification in z-domain is performed with the function 
  <a href="matlab:doc('ao/zDomainFit')">zDomainFit</a>. 
  It is based on a modeified version of the vector fitting algorithm that was 
  adapted to fit in z-domain. Details on the core agorithm can be found in [1 - 3].
</p>
<p>
  If you provide more than one AO as input, they will be fitted
  together with a common set of poles.
  Only frequency domain (<a href="matlab:doc('fsdata')">fsdata</a>) data can be
  fitted. Each non fsdata object is ignored. Input
  objects must have the same number of elements.
</p>


<h2><a name="algorithm">Fit algorithm</a></h2>

<p>
  The function performs a fitting loop to automatically identify model 
  order and parameters in z-domain. Output is a z-domain model expanded 
  in partial fractions:
</p>
<p>
  <div>
    <IMG src="images/zdomainfit_1.gif" border="0">
  </div>
</p>
<p>
  Each element of the partial fraction expansion can be seen as a
  <a href="sigproc_iir.html">miir</a> filter. Therefore the complete expansion
  is simply a parallel <a href="sigproc_filterbanks.html">filterbank</a> of
  <a href="sigproc_iir.html">miir</a> filters.
  Since the function can fit more than one input analysis object at a time
  with a common set of poles, output filterbank are embedded in a 
  <a href="class_desc_matrix.html">matrix</a> (note that this characteristic
  will be probably changed becausse of the introduction of the
  <a href="class_desc_collection.html">collection</a> class).
</p>
<p>
  Identification loop stops when the stop condition is reached. 
  Stop criterion is based on three different approaches:
  <ol>
    <li> Mean Squared Error and variation <br>
      Check if the normalized mean squared error is lower than the value specified in
      <tt>FITTOL</tt> and if the relative variation of the mean squared error is lower
      than the value specified in <tt>MSEVARTOL</tt>.
      E.g. <tt>FITTOL = 1e-3</tt>, <tt>MSEVARTOL = 1e-2</tt> search for a fit with
      normalized meam square error lower than <tt>1e-3</tt> and <tt>MSE</tt> relative
      variation lower than <tt>1e-2</tt>. 
    </li>
    <li> Log residuals difference and root mean squared error
      <ul>
        <li> Log Residuals difference <br>
        Check if the minimum of the logarithmic difference between data and 
        residuals is larger than a specified value. ie. if the conditioning 
        value is <tt>2</tt>, the function ensures that the difference between data and 
        residuals is at lest two order of magnitude lower than data itsleves. 
        <li> Root Mean Squared Error <br>
        Check that the variation of the root mean squared error is lower than 
        <tt>10^(-1*value)</tt>.
      </ul>
    </li>
    <li> Residuals spectral flatness and root mean squared error
      <ul>
        <li> Residuals Spectral Flatness <br>
        In case of a fit on noisy data, the residuals from a good fit are
        expected to be as much as possible similar to a white noise. This
        property can be used to test the accuracy of a fit procedure. In
        particular it can be tested that the spectral flatness coefficient of
        the residuals is larger than a certain qiantity sf such that <tt>0 < sf < 1</tt>.
        <li> Root Mean Squared Error <br>
        Check that the variation of the root mean squared error is lower than 
        <tt>10^(-1*value)</tt>.
      </ul>
    </li>
  </ol>
  Fitting loop stops when the two stopping conditions are satisfied, in both cases.
</p>
<p>
  The function can also perform a single loop without taking care of
  the stop conditions. This happens when <span class="string">'AUTOSEARCH'</span> parameter is
  set to <span class="string">'OFF'</span>.
</p>



<h2><a name="examples">Usage example of z-domain system identification tool</a></h2>
<p>
  In this example we fit a given frequency response to get a stable <tt>miir</tt> filter.
  For the meaning of any parameter please refer to
  <a href="matlab:doc('ao')">ao</a> and
  <a href="matlab:doc('ao/zDomainFit')">zDomainFit</a>
  documentation pages.
</p>

<div class="fragment"><pre>
    pl = plist(...
      <span class="string">'fsfcn'</span>, <span class="string">'(1e-3./(2.*pi.*1i.*f).^2 + 1e3./(0.001+2.*pi.*1i.*f) + 1e5.*(2.*pi.*1i.*f).^2).*1e-10'</span>,...
      <span class="string">'f1'</span>, 1e-6,...
      <span class="string">'f2'</span>, 5,...
      <span class="string">'nf'</span>, 100);
    
    a = ao(pl);
    a.setName;

    <span class="comment">% Fit parameter list</span>
    pl_fit = plist(<span class="string">'FS'</span>,10,...
      <span class="string">'AutoSearch'</span>,<span class="string">'on'</span>,...
      <span class="string">'StartPolesOpt'</span>,<span class="string">'clog'</span>,...
      <span class="string">'maxiter'</span>,50,...
      <span class="string">'minorder'</span>,15,...
      <span class="string">'maxorder'</span>,30,...
      <span class="string">'weightparam'</span>,<span class="string">'abs'</span>,...
      <span class="string">'CONDTYPE'</span>,<span class="string">'MSE'</span>,...
      <span class="string">'FITTOL'</span>,1e-2,...
      <span class="string">'MSEVARTOL'</span>,1e-1,...
      <span class="string">'Plot'</span>,<span class="string">'on'</span>,...
      <span class="string">'ForceStability'</span>,<span class="string">'on'</span>);

    <span class="comment">% Do fit</span>
    mod = zDomainFit(a, pl_fit);
</pre></div>

<p>
  <tt>mod</tt> is a <tt>matrix</tt> object containing a <tt>filterbank</tt> object.
</p>

<div class="fragment"><pre>
    >> mod
    ---- matrix 1 ----
    name: fit(a)
    size: 1x1
    01: filterbank | filterbank(fit(a)(fs=10.00, ntaps=2.00, a=[-1.19e+005 0], b=[1 0.0223]), fit(a)(fs=10.00, ntaps=2.00, a=[1.67e+005 0], b=[1 0.137]), fit(a)(fs=10.00, ntaps=2.00, a=[-5.41e+004 0], b=[1 0.348]), fit(a)(fs=10.00, ntaps=2.00, a=[1.15e+004 0], b=[1 0.603]), fit(a)(fs=10.00, ntaps=2.00, a=[-1.69e+005 0], b=[1 0.639]), fit(a)(fs=10.00, ntaps=2.00, a=[1.6e+005 0], b=[1 0.64]), fit(a)(fs=10.00, ntaps=2.00, a=[9.99e-009 0], b=[1 -1]), fit(a)(fs=10.00, ntaps=2.00, a=[-4.95e-010 0], b=[1 1]), fit(a)(fs=10.00, ntaps=2.00, a=[9.4e+003-i*3.7e+003 0], b=[1 -0.0528-i*0.0424]), fit(a)(fs=10.00, ntaps=2.00, a=[9.4e+003+i*3.7e+003 0], b=[1 -0.0528+i*0.0424]), fit(a)(fs=10.00, ntaps=2.00, a=[1.66e+003-i*1.45e+004 0], b=[1 0.0233-i*0.112]), fit(a)(fs=10.00, ntaps=2.00, a=[1.66e+003+i*1.45e+004 0], b=[1 0.0233+i*0.112]), fit(a)(fs=10.00, ntaps=2.00, a=[-1.67e+004+i*432 0], b=[1 0.171-i*0.14]), fit(a)(fs=10.00, ntaps=2.00, a=[-1.67e+004-i*432 0], b=[1 0.171+i*0.14]), fit(a)(fs=10.00, ntaps=2.00, a=[7.61e+003+i*7.36e+003 0], b=[1 0.378-i*0.112]), fit(a)(fs=10.00, ntaps=2.00, a=[7.61e+003-i*7.36e+003 0], b=[1 0.378+i*0.112]), fit(a)(fs=10.00, ntaps=2.00, a=[3.67e-015-i*4.61e-006 0], b=[1 -1-i*1.08e-010]), fit(a)(fs=10.00, ntaps=2.00, a=[3.67e-015+i*4.61e-006 0], b=[1 -1+i*1.08e-010]))
    description: 
    UUID: 9274455a-68e8-4bf1-b1ad-db81551f3cd6
    ------------------
</pre></div>

<p>
  The <tt>filterbank</tt> object contains a parallel bank of 18 filters.
</p>

<div class="fragment"><pre>
    >> mod.objs
    ---- filterbank 1 ----
    name: fit(a)
    type: parallel
    01: fit(a)(fs=10.00, ntaps=2.00, a=[-1.19e+005 0], b=[1 0.0223])
    02: fit(a)(fs=10.00, ntaps=2.00, a=[1.67e+005 0], b=[1 0.137])
    03: fit(a)(fs=10.00, ntaps=2.00, a=[-5.41e+004 0], b=[1 0.348])
    04: fit(a)(fs=10.00, ntaps=2.00, a=[1.15e+004 0], b=[1 0.603])
    05: fit(a)(fs=10.00, ntaps=2.00, a=[-1.69e+005 0], b=[1 0.639])
    06: fit(a)(fs=10.00, ntaps=2.00, a=[1.6e+005 0], b=[1 0.64])
    07: fit(a)(fs=10.00, ntaps=2.00, a=[9.99e-009 0], b=[1 -1])
    08: fit(a)(fs=10.00, ntaps=2.00, a=[-4.95e-010 0], b=[1 1])
    09: fit(a)(fs=10.00, ntaps=2.00, a=[9.4e+003-i*3.7e+003 0], b=[1 -0.0528-i*0.0424])
    10: fit(a)(fs=10.00, ntaps=2.00, a=[9.4e+003+i*3.7e+003 0], b=[1 -0.0528+i*0.0424])
    11: fit(a)(fs=10.00, ntaps=2.00, a=[1.66e+003-i*1.45e+004 0], b=[1 0.0233-i*0.112])
    12: fit(a)(fs=10.00, ntaps=2.00, a=[1.66e+003+i*1.45e+004 0], b=[1 0.0233+i*0.112])
    13: fit(a)(fs=10.00, ntaps=2.00, a=[-1.67e+004+i*432 0], b=[1 0.171-i*0.14])
    14: fit(a)(fs=10.00, ntaps=2.00, a=[-1.67e+004-i*432 0], b=[1 0.171+i*0.14])
    15: fit(a)(fs=10.00, ntaps=2.00, a=[7.61e+003+i*7.36e+003 0], b=[1 0.378-i*0.112])
    16: fit(a)(fs=10.00, ntaps=2.00, a=[7.61e+003-i*7.36e+003 0], b=[1 0.378+i*0.112])
    17: fit(a)(fs=10.00, ntaps=2.00, a=[3.67e-015-i*4.61e-006 0], b=[1 -1-i*1.08e-010])
    18: fit(a)(fs=10.00, ntaps=2.00, a=[3.67e-015+i*4.61e-006 0], b=[1 -1+i*1.08e-010])
    description: 
    UUID: 21af6960-61a8-4351-b504-e6f2b5e55b06
    ----------------------
</pre></div>

<p>
  Each object of the <tt>filterbank</tt> is a <tt>miir</tt> filter.
</p>

<div class="fragment"><pre>
    filt = mod.objs.filters.index(3)
    ------ miir/1 -------
    b: [1 0.348484501572296]
    histin: 0
    version: $Id: zdomainfit_content.html,v 1.6 2009/08/27 11:38:58 luigi Exp $
    ntaps: 2
    fs: 10
    infile: 
    a: [-54055.7700068032 0]
    histout: 0
    iunits: [] [1x1 unit]
    ounits: [] [1x1 unit]
    hist: miir.hist [1x1 history]
    procinfo: (empty-plist) [1x1 plist]
    plotinfo: (empty-plist) [1x1 plist]
    name: (fit(a)(3,1))(3)
    description: 
    mdlfile: 
    UUID: 6e2a1cd8-f17d-4c9d-aea9-4d9a96e41e68
    ---------------------
</pre></div>


<h2><a name="references">References</a></h2>
<p>
  <ol>
    <li> B. Gustavsen and A. Semlyen, "Rational approximation of frequency
    domain responses by Vector Fitting", IEEE Trans. Power Delivery
    vol. 14, no. 3, pp. 1052-1061, July 1999.
    <li> B. Gustavsen, "Improving the Pole Relocating Properties of Vector
    Fitting", IEEE Trans. Power Delivery vol. 21, no. 3, pp.
    1587-1592, July 2006.
    <li> Y. S. Mekonnen and J. E. Schutt-Aine, "Fast broadband
    macromodeling technique of sampled time/frequency data using
    z-domain vector-fitting method", Electronic Components and
    Technology Conference, 2008. ECTC 2008. 58th 27-30 May 2008 pp.
    1231 - 1235.
  </ol>
</p>