|
7017
|
1 ## Copyright (C) 1996, 1998, 2000, 2003, 2004, 2005, 2006, 2007 |
|
|
2 ## Auburn University. All rights reserved. |
|
3437
|
3 ## |
|
|
4 ## This file is part of Octave. |
|
|
5 ## |
|
|
6 ## Octave is free software; you can redistribute it and/or modify it |
|
7016
|
7 ## under the terms of the GNU General Public License as published by |
|
|
8 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
9 ## your option) any later version. |
|
3437
|
10 ## |
|
7016
|
11 ## Octave is distributed in the hope that it will be useful, but |
|
|
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
14 ## General Public License for more details. |
|
3437
|
15 ## |
|
|
16 ## You should have received a copy of the GNU General Public License |
|
7016
|
17 ## along with Octave; see the file COPYING. If not, see |
|
|
18 ## <http://www.gnu.org/licenses/>. |
|
3437
|
19 |
|
6945
|
20 ## Undocumented internal function. |
|
|
21 |
|
3437
|
22 ## -*- texinfo -*- |
|
3500
|
23 ## @deftypefn {Function File} {[@var{y}, @var{t}] =} __stepimp__ (@var{sitype}, @var{sys} [, @var{inp}, @var{tstop}, @var{n}]) |
|
3437
|
24 ## Impulse or step response for a linear system. |
|
|
25 ## The system can be discrete or multivariable (or both). |
|
5016
|
26 ## This m-file contains the ``common code'' of step and impulse. |
|
3437
|
27 ## |
|
5016
|
28 ## Produces a plot or the response data for system @var{sys}. |
|
3437
|
29 ## |
|
5016
|
30 ## Limited argument checking; ``do not attempt to do this at home''. |
|
|
31 ## Used internally in @command{impulse}, @command{step}. Use @command{step} |
|
|
32 ## or @command{impulse} instead. |
|
5642
|
33 ## @seealso{step, impulse} |
|
3437
|
34 ## @end deftypefn |
|
|
35 |
|
|
36 ## Author: Kai P. Mueller <mueller@ifr.ing.tu-bs.de> |
|
|
37 ## Created: October 2, 1997 |
|
|
38 ## based on lsim.m of Scottedward Hodel |
|
|
39 |
|
3438
|
40 function [y, t] = __stepimp__ (sitype, sys, inp, tstop, n) |
|
3437
|
41 |
|
6448
|
42 if (sitype == 1) |
|
|
43 IMPULSE = 0; |
|
|
44 elseif (sitype == 2) |
|
|
45 IMPULSE = 1; |
|
|
46 else |
|
|
47 error ("__stepimp__: invalid sitype argument"); |
|
3437
|
48 endif |
|
6448
|
49 sys = sysupdate (sys, "ss"); |
|
3437
|
50 |
|
|
51 USE_DEF = 0; # default tstop and n if we have to give up |
|
|
52 N_MIN = 50; # minimum number of points |
|
|
53 N_MAX = 2000; # maximum number of points |
|
|
54 T_DEF = 10.0; # default simulation time |
|
|
55 |
|
|
56 ## collect useful information about the system |
|
6448
|
57 [ncstates, ndstates, NIN, NOUT] = sysdimensions (sys); |
|
|
58 TSAMPLE = sysgettsam (sys); |
|
3437
|
59 |
|
6448
|
60 if (nargin < 3) |
|
|
61 inp = 1; |
|
|
62 elseif (inp < 1 || inp > NIN) |
|
|
63 error ("__stepimp__: argument inp out of range"); |
|
3437
|
64 endif |
|
|
65 |
|
6448
|
66 DIGITAL = is_digital (sys); |
|
3437
|
67 if (DIGITAL) |
|
|
68 NSTATES = ndstates; |
|
|
69 if (TSAMPLE < eps) |
|
6448
|
70 error ("__stepimp__: sampling time of discrete system too small") |
|
3437
|
71 endif |
|
6448
|
72 else |
|
|
73 NSTATES = ncstates; |
|
|
74 endif |
|
3437
|
75 if (NSTATES < 1) |
|
6448
|
76 error ("__stepimp__: pure gain block (n_states < 1), step response is trivial"); |
|
3437
|
77 endif |
|
|
78 if (nargin < 5) |
|
|
79 ## we have to compute the time when the system reaches steady state |
|
|
80 ## and the step size |
|
6448
|
81 ev = eig (sys2ss (sys)); |
|
3437
|
82 if (DIGITAL) |
|
|
83 ## perform bilinear transformation on poles in z |
|
|
84 for i = 1:NSTATES |
|
|
85 pole = ev(i); |
|
|
86 if (abs(pole + 1) < 1.0e-10) |
|
|
87 ev(i) = 0; |
|
|
88 else |
|
|
89 ev(i) = 2 / TSAMPLE * (pole - 1) / (pole + 1); |
|
|
90 endif |
|
|
91 endfor |
|
|
92 endif |
|
|
93 ## remove poles near zero from eigenvalue array ev |
|
|
94 nk = NSTATES; |
|
|
95 for i = 1:NSTATES |
|
6448
|
96 if (abs (real (ev(i))) < 1.0e-10) |
|
3437
|
97 ev(i) = 0; |
|
|
98 nk = nk - 1; |
|
|
99 endif |
|
|
100 endfor |
|
|
101 if (nk == 0) |
|
|
102 USE_DEF = 1; |
|
|
103 ## printf("##STEPIMP-DEBUG: using defaults.\n"); |
|
|
104 else |
|
6448
|
105 ev = ev(find (ev)); |
|
|
106 x = max (abs (ev)); |
|
3437
|
107 t_step = 0.2 * pi / x; |
|
6448
|
108 x = min (abs (real (ev))); |
|
3437
|
109 t_sim = 5.0 / x; |
|
|
110 ## round up |
|
6448
|
111 yy = 10^(ceil (log10 (t_sim)) - 1); |
|
|
112 t_sim = yy * ceil (t_sim / yy); |
|
3437
|
113 ## printf("##STEPIMP-DEBUG: nk=%d t_step=%f t_sim=%f\n", |
|
|
114 ## nk, t_step, t_sim); |
|
|
115 endif |
|
|
116 endif |
|
|
117 |
|
|
118 if (DIGITAL) |
|
|
119 ## ---- sampled system |
|
|
120 if (nargin == 5) |
|
6448
|
121 n = round (n); |
|
3437
|
122 if (n < 2) |
|
6448
|
123 error ("__stepimp__: n must not be less than 2.") |
|
3437
|
124 endif |
|
|
125 else |
|
|
126 if (nargin == 4) |
|
|
127 ## n is unknown |
|
|
128 elseif (nargin >= 1) |
|
|
129 ## tstop and n are unknown |
|
|
130 if (USE_DEF) |
|
|
131 tstop = (N_MIN - 1) * TSAMPLE; |
|
|
132 else |
|
|
133 tstop = t_sim; |
|
|
134 endif |
|
|
135 endif |
|
6448
|
136 n = floor (tstop / TSAMPLE) + 1; |
|
|
137 if (n < 2) |
|
|
138 n = 2; |
|
|
139 endif |
|
3437
|
140 if (n > N_MAX) |
|
|
141 n = N_MAX; |
|
6448
|
142 printf ("Hint: number of samples limited to %d by default.\n", \ |
|
|
143 N_MAX); |
|
|
144 printf (" ==> increase \"n\" parameter for longer simulations.\n"); |
|
3437
|
145 endif |
|
|
146 endif |
|
|
147 tstop = (n - 1) * TSAMPLE; |
|
|
148 t_step = TSAMPLE; |
|
|
149 else |
|
|
150 ## ---- continuous system |
|
|
151 if (nargin == 5) |
|
6448
|
152 n = round (n); |
|
3437
|
153 if (n < 2) |
|
|
154 error("step: n must not be less than 2.") |
|
|
155 endif |
|
|
156 t_step = tstop / (n - 1); |
|
|
157 else |
|
|
158 if (nargin == 4) |
|
|
159 ## only n in unknown |
|
|
160 if (USE_DEF) |
|
|
161 n = N_MIN; |
|
|
162 t_step = tstop / (n - 1); |
|
|
163 else |
|
6448
|
164 n = floor (tstop / t_step) + 1; |
|
3437
|
165 endif |
|
|
166 else |
|
|
167 ## tstop and n are unknown |
|
|
168 if (USE_DEF) |
|
|
169 tstop = T_DEF; |
|
|
170 n = N_MIN; |
|
|
171 t_step = tstop / (n - 1); |
|
|
172 else |
|
|
173 tstop = t_sim; |
|
6448
|
174 n = floor (tstop / t_step) + 1; |
|
3437
|
175 endif |
|
|
176 endif |
|
|
177 if (n < N_MIN) |
|
|
178 n = N_MIN; |
|
|
179 t_step = tstop / (n - 1); |
|
|
180 endif |
|
|
181 if (n > N_MAX) |
|
|
182 tstop = (n - 1) * t_step; |
|
|
183 t_step = tstop / (N_MAX - 1); |
|
|
184 n = N_MAX; |
|
|
185 endif |
|
|
186 endif |
|
|
187 tstop = (n - 1) * t_step; |
|
6448
|
188 [jnk,B] = sys2ss (sys); |
|
3437
|
189 B = B(:,inp); |
|
6448
|
190 sys = c2d (sys, t_step); |
|
3437
|
191 endif |
|
|
192 ## printf("##STEPIMP-DEBUG: t_step=%f n=%d tstop=%f\n", t_step, n, tstop); |
|
|
193 |
|
|
194 F = sys.a; |
|
|
195 G = sys.b(:,inp); |
|
|
196 C = sys.c; |
|
|
197 D = sys.d(:,inp); |
|
6448
|
198 y = zeros (NOUT, n); |
|
|
199 t = linspace (0, tstop, n); |
|
3437
|
200 |
|
|
201 if (IMPULSE) |
|
6448
|
202 if (! DIGITAL && D'*D > 0) |
|
|
203 error ("impulse: D matrix is nonzero, impulse response infinite.") |
|
3437
|
204 endif |
|
|
205 if (DIGITAL) |
|
|
206 y(:,1) = D / t_step; |
|
|
207 x = G / t_step; |
|
|
208 else |
|
|
209 x = B; |
|
|
210 y(:,1) = C * x; |
|
|
211 x = F * x; |
|
|
212 endif |
|
|
213 for i = 2:n |
|
|
214 y(:,i) = C * x; |
|
|
215 x = F * x; |
|
|
216 endfor |
|
4375
|
217 if (DIGITAL) |
|
|
218 y *= t_step; |
|
|
219 endif |
|
3437
|
220 else |
|
6448
|
221 x = zeros (NSTATES, 1); |
|
3437
|
222 for i = 1:n |
|
|
223 y(:,i) = C * x + D; |
|
|
224 x = F * x + G; |
|
|
225 endfor |
|
|
226 endif |
|
3679
|
227 |
|
6448
|
228 if (nargout == 0) |
|
|
229 if (IMPULSE) |
|
|
230 gm = zeros (NOUT, 1); |
|
|
231 tt = "impulse"; |
|
|
232 else |
|
|
233 ssys = ss (F, G, C, D, t_step); |
|
|
234 gm = dcgain (ssys); |
|
|
235 tt = "step"; |
|
|
236 endif |
|
|
237 ncols = floor (sqrt (NOUT)); |
|
|
238 nrows = ceil (NOUT / ncols); |
|
|
239 for i = 1:NOUT |
|
|
240 subplot (nrows, ncols, i); |
|
|
241 if (DIGITAL) |
|
|
242 [ts, ys] = stairs (t, y(i,:)); |
|
|
243 ts = ts(1:2*n-2)'; |
|
|
244 ys = ys(1:2*n-2)'; |
|
|
245 if (length (gm) > 0) |
|
|
246 yy = [ys; gm(i)*ones(size(ts))]; |
|
|
247 else |
|
|
248 yy = ys; |
|
|
249 endif |
|
|
250 plot (ts, yy); |
|
|
251 grid ("on"); |
|
|
252 xlabel ("time [s]"); |
|
|
253 ylabel ("y(t)"); |
|
4422
|
254 else |
|
6448
|
255 if (length (gm) > 0) |
|
|
256 yy = [y(i,:); gm(i)*ones(size(t))]; |
|
|
257 else |
|
|
258 yy = y(i,:); |
|
|
259 endif |
|
|
260 plot (t, yy); |
|
|
261 grid ("on"); |
|
|
262 xlabel ("time [s]"); |
|
|
263 ylabel ("y(t)"); |
|
4789
|
264 endif |
|
6448
|
265 title (sprintf ("%s: | %s -> %s", tt, |
|
|
266 sysgetsignals (sys, "in", inp, 1), |
|
|
267 sysgetsignals (sys, "out", i, 1))); |
|
|
268 endfor |
|
|
269 y = []; |
|
|
270 t = []; |
|
|
271 endif |
|
7126
|
272 |
|
3679
|
273 endfunction |
