1 ## Copyright (C) 2011 Lukas F. Reichlin
3 ## This file is part of LTI Syncope.
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.
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.
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/>.
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.
27 ## LTI model containing frequency response data of a SISO system.
29 ## The desired order of the system to be fitted. @code{n <= length(dat.w)}.
31 ## The flag controls whether the returned system is stable and minimum-phase.
34 ## The system zeros and poles are not constrained. Default value.
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.
44 ## State-space model of order @var{n}, fitted to frequency response data @var{dat}.
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}.
50 ## @strong{Algorithm}@*
51 ## Uses SLICOT SB10YD by courtesy of
52 ## @uref{http://www.slicot.org, NICONET e.V.}
55 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
56 ## Created: October 2011
59 function [sys, n] = fitfrd (dat, n, flag = 0)
61 if (nargin == 0 || nargin > 3)
65 if (! isa (dat, "frd"))
70 error ("fitfrd: require SISO system");
73 if (! issample (n, 0) || n != round (n))
74 error ("fitfrd: second argument must be an integer >= 0");
77 [H, w, tsam] = frdata (dat, "vector");
81 error ("fitfrd: require n <= length (dat.w)");
84 [a, b, c, d, n] = slsb10yd (real (H), imag (H), w, n, dt, logical (flag));
86 sys = ss (a, b, c, d, tsam);
92 %! SYS = ss (-1, 1, 1, 0);
94 %! Ye = step (SYS, T);
95 %! W = logspace (-2, 2, 100);
98 %! SYSID = fitfrd (FR, N, 1);
99 %! Yo = step (SYSID, T);
100 %!assert (Yo, Ye, 1e-2);