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