Mercurial > octave-nkf
annotate scripts/control/base/__stepimp__.m @ 11658:db22340e1f24 release-3-0-x
__stepimp__: don't call subplot for single plot
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Fri, 22 Feb 2008 04:02:37 -0500 |
parents | 4a375de63f66 |
children | df9519e9990c 72830070a17b |
rev | line source |
---|---|
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 | |
11658
db22340e1f24
__stepimp__: don't call subplot for single plot
John W. Eaton <jwe@octave.org>
parents:
7126
diff
changeset
|
240 if (nrows > 1 || ncols > 1) |
db22340e1f24
__stepimp__: don't call subplot for single plot
John W. Eaton <jwe@octave.org>
parents:
7126
diff
changeset
|
241 subplot (nrows, ncols, i); |
db22340e1f24
__stepimp__: don't call subplot for single plot
John W. Eaton <jwe@octave.org>
parents:
7126
diff
changeset
|
242 endif |
6448 | 243 if (DIGITAL) |
244 [ts, ys] = stairs (t, y(i,:)); | |
245 ts = ts(1:2*n-2)'; | |
246 ys = ys(1:2*n-2)'; | |
247 if (length (gm) > 0) | |
248 yy = [ys; gm(i)*ones(size(ts))]; | |
249 else | |
250 yy = ys; | |
251 endif | |
252 plot (ts, yy); | |
253 grid ("on"); | |
254 xlabel ("time [s]"); | |
255 ylabel ("y(t)"); | |
4422 | 256 else |
6448 | 257 if (length (gm) > 0) |
258 yy = [y(i,:); gm(i)*ones(size(t))]; | |
259 else | |
260 yy = y(i,:); | |
261 endif | |
262 plot (t, yy); | |
263 grid ("on"); | |
264 xlabel ("time [s]"); | |
265 ylabel ("y(t)"); | |
4789 | 266 endif |
6448 | 267 title (sprintf ("%s: | %s -> %s", tt, |
268 sysgetsignals (sys, "in", inp, 1), | |
269 sysgetsignals (sys, "out", i, 1))); | |
270 endfor | |
271 y = []; | |
272 t = []; | |
273 endif | |
7126 | 274 |
3679 | 275 endfunction |