1 ## Copyright (C) 2009, 2010 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 ## Common code for the time response functions step, impulse and initial.
21 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
22 ## Created: October 2009
25 function [y, t, x_arr] = __time_response__ (sys, resptype, plotflag, tfinal, dt, x0, sysname)
27 if (! isa (sys, "ss"))
28 sys = ss (sys); # sys must be proper
31 if (is_real_vector (tfinal) && length (tfinal) > 1) # time vector t passed
32 dt = tfinal(2) - tfinal(1); # assume that t is regularly spaced
36 [A, B, C, D, tsam] = ssdata (sys);
38 discrete = ! isct (sys); # static gains are treated as analog systems
39 tsam = abs (tsam); # use 1 second if tsam is unspecified (-1)
43 warning ("time_response: argument dt has no effect on sampling time of discrete system");
49 [tfinal, dt] = __sim_horizon__ (A, discrete, tfinal, dt);
52 sys = c2d (sys, dt, "zoh");
55 [F, G] = ssdata (sys); # matrices C and D don't change
57 n = rows (F); # number of states
58 m = columns (G); # number of inputs
59 p = rows (C); # number of outputs
62 t = reshape (0 : dt : tfinal, [], 1);
67 str = ["Response of ", sysname, " to Initial Conditions"];
68 yfinal = zeros (p, 1);
72 x_arr = zeros (l_t, n);
75 x = reshape (x0, [], 1); # make sure that x is a column vector
77 if (n != length (x0) || ! is_real_vector (x0))
78 error ("initial: x0 must be a real vector with %d elements", n);
89 str = ["Step Response of ", sysname];
90 yfinal = dcgain (sys);
93 y = zeros (l_t, p, m);
94 x_arr = zeros (l_t, n, m);
96 for j = 1 : m # for every input channel
104 y(k, :, j) = C * x + D * u;
111 str = ["Impulse Response of ", sysname];
112 yfinal = zeros (p, m);
114 ## preallocate memory
115 y = zeros (l_t, p, m);
116 x_arr = zeros (l_t, n, m);
118 for j = 1 : m # for every input channel
119 ## initial conditions
124 x = zeros (n, 1); # zero by definition
125 y(1, :, j) = D * u / dt;
129 x = B * u; # B, not G!
149 error ("time_response: invalid response type");
154 if (plotflag) # display plot
156 ## TODO: Set correct titles, especially for multi-input systems
158 stable = isstable (sys);
159 outname = get (sys, "outname");
160 outname = __labels__ (outname, "y_");
162 if (strcmp (resptype, "initial"))
168 if (discrete) # discrete system
172 subplot (p, cols, (k-1)*cols+j);
175 stairs (t, [y(:, k, j), yfinal(k, j) * ones(l_t, 1)]);
177 stairs (t, y(:, k, j));
182 if (k == 1 && j == 1)
187 ylabel (sprintf ("Amplitude %s", outname{k}));
195 else # continuous system
199 subplot (p, cols, (k-1)*cols+j);
202 plot (t, [y(:, k, j), yfinal(k, j) * ones(l_t, 1)]);
204 plot (t, y(:, k, j));
209 if (k == 1 && j == 1)
214 ylabel (sprintf ("Amplitude %s", outname{k}));
228 function [tfinal, dt] = __sim_horizon__ (A, discrete, tfinal, Ts)
230 ## code based on __stepimp__.m of Kai P. Mueller and A. Scottedward Hodel
232 TOL = 1.0e-10; # values below TOL are assumed to be zero
233 N_MIN = 50; # min number of points
234 N_MAX = 2000; # max number of points
235 N_DEF = 1000; # default number of points
236 T_DEF = 10; # default simulation time
242 ## perform bilinear transformation on poles in z
245 if (abs (pol + 1) < TOL)
248 eigw(k) = 2 / Ts * (pol - 1) / (pol + 1);
253 ## remove poles near zero from eigenvalue array eigw
256 if (abs (real (eigw(k))) < TOL)
263 if (isempty (tfinal))
271 eigw = eigw(find (eigw));
272 eigw_max = max (abs (eigw));
275 dt = 0.2 * pi / eigw_max;
278 if (isempty (tfinal))
279 eigw_min = min (abs (real (eigw)));
280 tfinal = 5.0 / eigw_min;
283 yy = 10^(ceil (log10 (tfinal)) - 1);
284 tfinal = yy * ceil (tfinal / yy);
300 if (! isempty (Ts)) # catch case cont. system with dt specified