home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts / control / stepimp.m < prev    next >
Encoding:
Text File  |  1999-12-15  |  7.2 KB  |  281 lines
  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. 
  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
  32.  
  33. ## See also: step, impulse
  34.  
  35. function [y, t] = stepimp (sitype, sys, inp, tstop, n)
  36.  
  37.   ## Written by Kai P. Mueller October 2, 1997
  38.   ## based on lsim.m of Scottedward Hodel
  39.  
  40.   if (sitype == 1)         IMPULSE = 0;
  41.   elseif (sitype == 2)     IMPULSE = 1;
  42.   else                  error("stepimp: illegal sitype argument.")
  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
  50.  
  51.   ## collect useful information about the system
  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")
  57.   endif
  58.  
  59.   DIGITAL = is_digital(sys);
  60.   if (DIGITAL)
  61.     NSTATES = ndstates;
  62.     if (TSAMPLE < eps)
  63.       error("stepimp: sampling time of discrete system too small.")
  64.     endif
  65.   else        NSTATES = ncstates;       endif
  66.   if (NSTATES < 1)
  67.     error("step: pure gain block (n_states < 1), step response is trivial");
  68.   endif
  69.   if (nargin < 5)
  70.     ## we have to compute the time when the system reaches steady state
  71.     ## and the step size
  72.     ev = eig(sys2ss(sys));
  73.     if (DIGITAL)
  74.       ## perform bilinear transformation on poles in z
  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
  84.     ## remove poles near zero from eigenvalue array ev
  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;
  94.       ## printf("##STEPIMP-DEBUG: using defaults.\n");
  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;
  101.       ## round up
  102.       yy = 10^(ceil(log10(t_sim)) - 1);
  103.       t_sim = yy * ceil(t_sim / yy);
  104.       ## printf("##STEPIMP-DEBUG: nk=%d   t_step=%f  t_sim=%f\n",
  105.       ##   nk, t_step, t_sim);  
  106.     endif
  107.   endif
  108.  
  109.   if (DIGITAL)
  110.     ## ---- sampled system
  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)
  118.         ## n is unknown
  119.       elseif (nargin >= 1)
  120.         ## tstop and n are unknown
  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
  139.     ## ---- continuous system
  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)
  148.         ## only n in unknown
  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
  156.         ## tstop and n are unknown
  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;
  177.     [jnk,B] = sys2ss(sys);
  178.     B = B(:,inp);
  179.     sys = c2d(sys, t_step);
  180.   endif
  181.   ## printf("##STEPIMP-DEBUG: t_step=%f n=%d  tstop=%f\n", t_step, n, tstop);
  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)
  195.       y(:,1) = D / t_step;
  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)
  215.     ## Plot the information
  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);
  235.     title(sprintf("%s: | %s -> %s", tt,sysgetsignals(sys,"in",inp,1), ...
  236.       sysgetsignals(sys,"out",i,1)));
  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
  261.       ## leave gnuplot in multiplot mode is bad style
  262.       oneplot();
  263.     else
  264.       ## plot everything in one diagram
  265.       title([tt, " response | ", sysgetsignals(sys,"in",inp,1), ...
  266.     " -> all outputs"]);
  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
  279.   ## printf("##STEPIMP-DEBUG: gratulations, successfull completion.\n");
  280. endfunction
  281.