home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts / control / stepimp.m < prev    next >
Text File  |  1999-03-05  |  7KB  |  276 lines
  1. # Copyright (C) 1996,1998 A. Scottedward Hodel
  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, 675 Mass Ave, Cambridge, MA 02139, USA. 
  18.  
  19. function [y, t] = stepimp(sitype, sys, inp, tstop, n)
  20. # step: Impulse or step response for a linear system.
  21. #       The system can be discrete or multivariable (or both).
  22. #       This m-file contains the "common code" of step and impulse.
  23. #
  24. # [y, t] = stepimp(sitype, sys[, inp, tstop, n])
  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. # Use step or impulse instead.
  29. #
  30. # See also:  step, impulse
  31.  
  32. # Written by Kai P. Mueller October 2, 1997
  33. # based on lsim.m of Scottedward Hodel
  34.  
  35.   if (sitype == 1)         IMPULSE = 0;
  36.   elseif (sitype == 2)     IMPULSE = 1;
  37.   else                  error("stepimp: illegal sitype argument.")
  38.   endif
  39.   sys = sysupdate(sys,"ss");
  40.  
  41.   USE_DEF = 0;   # default tstop and n if we have to give up
  42.   N_MIN = 50;    # minimum number of points
  43.   N_MAX = 2000;  # maximum number of points
  44.   T_DEF = 10.0;  # default simulation time
  45.  
  46.   # collect useful information about the system
  47.   [ncstates,ndstates,NIN,NOUT] = sysdimensions(sys);
  48.   TSAMPLE = sysgettsam(sys);
  49.  
  50.   if (nargin < 3)                      inp = 1;
  51.   elseif (inp < 1 | inp > NIN)         error("Argument inp out of range")
  52.   endif
  53.  
  54.   DIGITAL = is_digital(sys);
  55.   if (DIGITAL)
  56.     NSTATES = ndstates;
  57.     if (TSAMPLE < eps)
  58.       error("stepimp: sampling time of discrete system too small.")
  59.     endif
  60.   else        NSTATES = ncstates;       endif
  61.   if (NSTATES < 1)
  62.     error("step: pure gain block (n_states < 1), step response is trivial");
  63.   endif
  64.   if (nargin < 5)
  65.     # we have to compute the time when the system reaches steady state
  66.     # and the step size
  67.     ev = eig(sys2ss(sys));
  68.     if (DIGITAL)
  69.       # perform bilinear transformation on poles in z
  70.       for i = 1:NSTATES
  71.         pole = ev(i);
  72.     if (abs(pole + 1) < 1.0e-10)
  73.       ev(i) = 0;
  74.     else
  75.       ev(i) = 2 / TSAMPLE * (pole - 1) / (pole + 1);
  76.     endif
  77.       endfor
  78.     endif
  79.     # remove poles near zero from eigenvalue array ev
  80.     nk = NSTATES;
  81.     for i = 1:NSTATES
  82.       if (abs(ev(i)) < 1.0e-10)
  83.         ev(i) = 0;
  84.         nk = nk - 1;
  85.       endif
  86.     endfor
  87.     if (nk == 0)
  88.       USE_DEF = 1;
  89.       #printf("##STEPIMP-DEBUG: using defaults.\n");
  90.     else
  91.       ev = ev(find(ev));
  92.       x = max(abs(ev));
  93.       t_step = 0.2 * pi / x;
  94.       x = min(abs(real(ev)));
  95.       t_sim = 5.0 / x;
  96.       # round up
  97.       yy = 10^(ceil(log10(t_sim)) - 1);
  98.       t_sim = yy * ceil(t_sim / yy);
  99.       #printf("##STEPIMP-DEBUG: nk=%d   t_step=%f  t_sim=%f\n",
  100.       #       nk, t_step, t_sim);  
  101.     endif
  102.   endif
  103.  
  104.   if (DIGITAL)
  105.     # ---- sampled system
  106.     if (nargin == 5)
  107.       n = round(n);
  108.       if (n < 2)
  109.         error("stepimp: n must not be less than 2.")
  110.       endif
  111.     else
  112.       if (nargin == 4)
  113.         # n is unknown
  114.       elseif (nargin >= 1)
  115.         # tstop and n are unknown
  116.         if (USE_DEF)
  117.           tstop = (N_MIN - 1) * TSAMPLE;
  118.         else
  119.           tstop = t_sim;
  120.         endif
  121.       endif
  122.       n = floor(tstop / TSAMPLE) + 1;
  123.       if (n < 2)  n = 2;  endif
  124.       if (n > N_MAX)
  125.     n = N_MAX;
  126.     printf("Hint: number of samples limited to %d by default.\n", \
  127.            N_MAX);
  128.     printf("  ==> increase \"n\" parameter for longer simulations.\n");
  129.       endif
  130.     endif
  131.     tstop = (n - 1) * TSAMPLE;
  132.     t_step = TSAMPLE;
  133.   else
  134.     # ---- continuous system
  135.     if (nargin == 5)
  136.       n = round(n);
  137.       if (n < 2)
  138.         error("step: n must not be less than 2.")
  139.       endif
  140.       t_step = tstop / (n - 1);
  141.     else
  142.       if (nargin == 4)
  143.         # only n in unknown
  144.         if (USE_DEF)
  145.           n = N_MIN;
  146.       t_step = tstop / (n - 1);
  147.         else
  148.           n = floor(tstop / t_step) + 1;
  149.         endif
  150.       else
  151.         # tstop and n are unknown
  152.         if (USE_DEF)
  153.           tstop = T_DEF;
  154.       n = N_MIN;
  155.       t_step = tstop / (n - 1);
  156.         else
  157.           tstop = t_sim;
  158.           n = floor(tstop / t_step) + 1;
  159.         endif
  160.       endif
  161.       if (n < N_MIN)
  162.     n = N_MIN;
  163.         t_step = tstop / (n - 1);
  164.       endif
  165.       if (n > N_MAX)
  166.         tstop = (n - 1) * t_step;
  167.     t_step = tstop / (N_MAX - 1);
  168.     n = N_MAX;
  169.       endif
  170.     endif
  171.     tstop = (n - 1) * t_step;
  172.     [jnk,B] = sys2ss(sys);
  173.     B = B(:,inp);
  174.     sys = c2d(sys, t_step);
  175.   endif
  176.   #printf("##STEPIMP-DEBUG: t_step=%f n=%d  tstop=%f\n", t_step, n, tstop);
  177.  
  178.   F = sys.a;
  179.   G = sys.b(:,inp);
  180.   C = sys.c;
  181.   D = sys.d(:,inp);
  182.   y = zeros(NOUT, n);
  183.   t = linspace(0, tstop, n);
  184.  
  185.   if (IMPULSE)
  186.     if (!DIGITAL && (D'*D > 0))
  187.       error("impulse: D matrix is nonzero, impulse response infinite.")
  188.     endif
  189.     if (DIGITAL)
  190.       y(:,1) = D;
  191.       x = G / t_step;
  192.     else
  193.       x = B;
  194.       y(:,1) = C * x;
  195.       x = F * x;
  196.     endif
  197.     for i = 2:n
  198.       y(:,i) = C * x;
  199.       x = F * x;
  200.     endfor
  201.   else
  202.     x = zeros(NSTATES, 1);
  203.     for i = 1:n
  204.       y(:,i) = C * x + D;
  205.       x = F * x + G;
  206.     endfor
  207.   endif
  208.  
  209.   if(nargout == 0)
  210.     # Plot the information
  211.     oneplot();
  212.     gset nogrid
  213.     gset nologscale
  214.     gset autoscale
  215.     gset nokey
  216.     clearplot();
  217.     if (gnuplot_has_multiplot)
  218.       if (IMPULSE)
  219.         gm = zeros(NOUT, 1);
  220.     tt = "impulse";
  221.       else
  222.         ssys = ss2sys(F, G, C, D, t_step);
  223.         gm = dcgain(ssys);
  224.     tt = "step";
  225.       endif
  226.       ncols = floor(sqrt(NOUT));
  227.       nrows = ceil(NOUT / ncols);
  228.       for i = 1:NOUT
  229.         subplot(nrows, ncols, i);
  230.     title(sprintf("%s: | %s -> %s", tt,sysgetsignals(sys,"in",inp,1), ...
  231.       sysgetsignals(sys,"out",i,1)));
  232.     if (DIGITAL)
  233.       [ts, ys] = stairs(t, y(i,:));
  234.       ts = ts(1:2*n-2)';  ys = ys(1:2*n-2)';
  235.       if (length(gm) > 0)
  236.         yy = [ys; gm(i)*ones(size(ts))];
  237.       else
  238.         yy = ys;
  239.       endif
  240.       grid("on");
  241.       xlabel("time [s]");
  242.       ylabel("y(t)");
  243.           plot(ts, yy);
  244.     else
  245.       if (length(gm) > 0)
  246.         yy = [y(i,:); gm(i)*ones(size(t))];
  247.       else
  248.         yy = y(i,:);
  249.       endif
  250.       grid("on");
  251.       xlabel("time [s]");
  252.       ylabel("y(t)");
  253.       plot(t, yy);
  254.     endif
  255.       endfor
  256.       # leave gnuplot in multiplot mode is bad style
  257.       oneplot();
  258.     else
  259.       # plot everything in one diagram
  260.       title([tt, " response | ", sysgetsignals(sys,"in",inp,1), ...
  261.     " -> all outputs"]);
  262.       if (DIGITAL)
  263.         stairs(t, y(i,:));
  264.       else
  265.     grid("on");
  266.     xlabel("time [s]");
  267.     ylabel("y(t)");
  268.     plot(t, y(i,:));
  269.       endif
  270.     endif
  271.     y=[];
  272.     t=[];
  273.   endif
  274.   #printf("##STEPIMP-DEBUG: gratulations, successfull completion.\n");
  275. endfunction
  276.