octave_packages/control-2.3.52/fitfrd.m

   1 ## Copyright (C) 2011   Lukas F. Reichlin
   2 ##
   3 ## This file is part of LTI Syncope.
   4 ##
   5 ## LTI Syncope is free software: you can redistribute it and/or modify
   6 ## it under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation, either version 3 of the License, or
   8 ## (at your option) any later version.
   9 ##
  10 ## LTI Syncope is distributed in the hope that it will be useful,
  11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 ## GNU General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with LTI Syncope.  If not, see <http://www.gnu.org/licenses/>.
  17
  18 ## -*- texinfo -*-
  19 ## @deftypefn{Function File} {[@var{sys}, @var{n}] =} fitfrd (@var{dat}, @var{n})
  20 ## @deftypefnx{Function File} {[@var{sys}, @var{n}] =} fitfrd (@var{dat}, @var{n}, @var{flag})
  21 ## Fit frequency response data with a state-space system.
  22 ## If requested, the returned system is stable and minimum-phase.
  23 ##
  24 ## @strong{Inputs}
  25 ## @table @var
  26 ## @item dat
  27 ## LTI model containing frequency response data of a SISO system.
  28 ## @item n
  29 ## The desired order of the system to be fitted.  @code{n <= length(dat.w)}.
  30 ## @item flag
  31 ## The flag controls whether the returned system is stable and minimum-phase.
  32 ## @table @var
  33 ## @item 0
  34 ## The system zeros and poles are not constrained.  Default value.
  35 ## @item 1
  36 ## The system zeros and poles will have negative real parts in the
  37 ## continuous-time case, or moduli less than 1 in the discrete-time case.
  38 ## @end table
  39 ## @end table
  40 ##
  41 ## @strong{Outputs}
  42 ## @table @var
  43 ## @item sys
  44 ## State-space model of order @var{n}, fitted to frequency response data @var{dat}.
  45 ## @item n
  46 ## The order of the obtained system.  The value of @var{n}
  47 ## could only be modified if inputs @code{n > 0} and @code{flag = 1}.
  48 ## @end table
  49 ##
  50 ## @strong{Algorithm}@*
  51 ## Uses SLICOT SB10YD by courtesy of
  52 ## @uref{http://www.slicot.org, NICONET e.V.}
  53 ## @end deftypefn
  54
  55 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  56 ## Created: October 2011
  57 ## Version: 0.1
  58
  59 function [sys, n] = fitfrd (dat, n, flag = 0)
  60
  61   if (nargin == 0 || nargin > 3)
  62     print_usage ();
  63   endif
  64
  65   if (! isa (dat, "frd"))
  66     dat = frd (dat);
  67   endif
  68
  69   if (! issiso (dat))
  70     error ("fitfrd: require SISO system");
  71   endif
  72
  73   if (! issample (n, 0) || n != round (n))
  74     error ("fitfrd: second argument must be an integer >= 0");
  75   endif
  76
  77   [H, w, tsam] = frdata (dat, "vector");
  78   dt = isdt (dat);
  79
  80   if (n > length (w))
  81     error ("fitfrd: require n <= length (dat.w)");
  82   endif
  83
  84   [a, b, c, d, n] = slsb10yd (real (H), imag (H), w, n, dt, logical (flag));
  85
  86   sys = ss (a, b, c, d, tsam);
  87
  88 endfunction
  89
  90
  91 %!shared Yo, Ye
  92 %! SYS = ss (-1, 1, 1, 0);
  93 %! T = 0:0.1:50;
  94 %! Ye = step (SYS, T);
  95 %! W = logspace (-2, 2, 100);
  96 %! FR = frd (SYS, W);
  97 %! N = 1;
  98 %! SYSID = fitfrd (FR, N, 1);
  99 %! Yo = step (SYSID, T);
 100 %!assert (Yo, Ye, 1e-2);