← Back to branch summary

~ubuntu-branches/debian/sid/octave3.0/sid

« back to all changes in this revision

Viewing changes to scripts/control/base/__stepimp__.m

  • Committer: Bazaar Package Importer
  • Author(s): Rafael Laboissiere
  • Date: 2007-12-23 16:04:15 UTC
  • Revision ID: james.westby@ubuntu.com-20071223160415-n4gk468dihy22e9v
Tags: upstream-3.0.0
Import upstream version 3.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
scripts/control/base/__stepimp__.m
 
1
## Copyright (C) 1996, 1998, 2000, 2003, 2004, 2005, 2006, 2007
 
2
##               Auburn University.  All rights reserved.
 
3
##
 
4
## This file is part of Octave.
 
5
##
 
6
## Octave is free software; you can redistribute it and/or modify it
 
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.
 
10
##
 
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.
 
15
##
 
16
## You should have received a copy of the GNU General Public License
 
17
## along with Octave; see the file COPYING.  If not, see
 
18
## <http://www.gnu.org/licenses/>.
 
19
 
 
20
## Undocumented internal function.
 
21
 
 
22
## -*- texinfo -*-
 
23
## @deftypefn {Function File} {[@var{y}, @var{t}] =} __stepimp__ (@var{sitype}, @var{sys} [, @var{inp}, @var{tstop}, @var{n}])
 
24
## Impulse or step response for a linear system.
 
25
## The system can be discrete or multivariable (or both).
 
26
## This m-file contains the ``common code'' of step and impulse.
 
27
##
 
28
## Produces a plot or the response data for system @var{sys}.
 
29
##
 
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.
 
33
## @seealso{step, impulse}
 
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
 
 
40
function [y, t] = __stepimp__ (sitype, sys, inp, tstop, n)
 
41
 
 
42
  if (sitype == 1)
 
43
    IMPULSE = 0;
 
44
  elseif (sitype == 2)
 
45
    IMPULSE = 1;
 
46
  else
 
47
    error ("__stepimp__: invalid sitype argument");
 
48
  endif
 
49
  sys = sysupdate (sys, "ss");
 
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
 
57
  [ncstates, ndstates, NIN, NOUT] = sysdimensions (sys);
 
58
  TSAMPLE = sysgettsam (sys);
 
59
 
 
60
  if (nargin < 3)
 
61
    inp = 1;
 
62
  elseif (inp < 1 || inp > NIN)
 
63
    error ("__stepimp__: argument inp out of range");
 
64
  endif
 
65
 
 
66
  DIGITAL = is_digital (sys);
 
67
  if (DIGITAL)
 
68
    NSTATES = ndstates;
 
69
    if (TSAMPLE < eps)
 
70
      error ("__stepimp__: sampling time of discrete system too small")
 
71
    endif
 
72
  else
 
73
    NSTATES = ncstates;
 
74
  endif
 
75
  if (NSTATES < 1)
 
76
    error ("__stepimp__: pure gain block (n_states < 1), step response is trivial");
 
77
  endif
 
78
  if (nargin < 5)
 
79
    ## we have to compute the time when the system reaches steady state
 
80
    ## and the step size
 
81
    ev = eig (sys2ss (sys));
 
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
 
96
      if (abs (real (ev(i))) < 1.0e-10)
 
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
 
105
      ev = ev(find (ev));
 
106
      x = max (abs (ev));
 
107
      t_step = 0.2 * pi / x;
 
108
      x = min (abs (real (ev)));
 
109
      t_sim = 5.0 / x;
 
110
      ## round up
 
111
      yy = 10^(ceil (log10 (t_sim)) - 1);
 
112
      t_sim = yy * ceil (t_sim / yy);
 
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)
 
121
      n = round (n);
 
122
      if (n < 2)
 
123
        error ("__stepimp__: n must not be less than 2.")
 
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
 
136
      n = floor (tstop / TSAMPLE) + 1;
 
137
      if (n < 2)
 
138
        n = 2;
 
139
      endif
 
140
      if (n > N_MAX)
 
141
        n = N_MAX;
 
142
        printf ("Hint: number of samples limited to %d by default.\n", \
 
143
                N_MAX);
 
144
        printf ("  ==> increase \"n\" parameter for longer simulations.\n");
 
145
      endif
 
146
    endif
 
147
    tstop = (n - 1) * TSAMPLE;
 
148
    t_step = TSAMPLE;
 
149
  else
 
150
    ## ---- continuous system
 
151
    if (nargin == 5)
 
152
      n = round (n);
 
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
 
164
          n = floor (tstop / t_step) + 1;
 
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;
 
174
          n = floor (tstop / t_step) + 1;
 
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;
 
188
    [jnk,B] = sys2ss (sys);
 
189
    B = B(:,inp);
 
190
    sys = c2d (sys, t_step);
 
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);
 
198
  y = zeros (NOUT, n);
 
199
  t = linspace (0, tstop, n);
 
200
 
 
201
  if (IMPULSE)
 
202
    if (! DIGITAL && D'*D > 0)
 
203
      error ("impulse: D matrix is nonzero, impulse response infinite.")
 
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
 
217
    if (DIGITAL)
 
218
      y *= t_step; 
 
219
    endif 
 
220
  else
 
221
    x = zeros (NSTATES, 1);
 
222
    for i = 1:n
 
223
      y(:,i) = C * x + D;
 
224
      x = F * x + G;
 
225
    endfor
 
226
  endif
 
227
  
 
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)");
 
254
      else
 
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)");
 
264
      endif
 
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
 
272
 
 
273
endfunction