home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / impulse2 < prev    next >
Encoding:
Text File  |  1995-11-15  |  2.9 KB  |  132 lines
  1. //-----------------------------------------------------------------------------
  2. //
  3. // impulse
  4. //
  5. // Syntax: impulse(a,b,c,d,ntim,iu)
  6. //
  7. // This routine finds the impulse response of a continuous-time
  8. // linear system given by,
  9. //         .
  10. //         x = Ax + Bu
  11. //         y = Cx + Du
  12. //
  13. // Calling the routine as,
  14. //
  15. //     impulse(A,B,C,D,NTIM)
  16. //
  17. // plots the time response of the linear system to an impulse
  18. // response applied to each input. The variable NTIM is the number
  19. // of time steps to be taken in the response.
  20. //
  21. // Calling the routine as,
  22. //
  23. //     impulse(A,B,C,D,NTIM,IU)
  24. //
  25. // plots the time response of the linear system to an impulse apllied
  26. // to the single input IU.
  27. //
  28. // The routine can also be used to find the impulse response of a polynomial
  29. // transfer function as,
  30. //
  31. //     impulse(NUM,DEN,NTIM)
  32. //
  33. // Version JBL 940519
  34. //-----------------------------------------------------------------------------
  35.  
  36. rfile tfchk
  37. rfile tf2ss
  38. rfile abcdchk
  39. rfile c2d
  40. rfile simdata
  41. rfile timeplot
  42. rfile ltitr
  43.  
  44. impulse2 = function(a,b,c,d,ntim,iu)
  45. {
  46.    local(nargs,Adum,num,den,nu,ny,msg,estr,aa,bb,cc,dd,Y,i,j,f,iiu,...
  47.          x,y,t,A,B,C,D)
  48.  
  49. // Count number of input arguments
  50.    nargs=0;
  51.    if (exist(a)) {nargs=nargs+1;}
  52.    if (exist(b)) {nargs=nargs+1;}
  53.    if (exist(c)) {nargs=nargs+1;}
  54.    if (exist(d)) {nargs=nargs+1;}
  55.    if (exist(ntim)) {nargs=nargs+1;}
  56.    if (exist(iu)) {nargs=nargs+1;}
  57.  
  58.    if ( (nargs <= 2) || (nargs > 6) ) {
  59.        error("impulse: Error in number of input arguments.");
  60.    }
  61.  
  62. // Convert to S-S if in TF form
  63.    if (nargs < 5) {
  64.        Adum=tfchk(a,b);
  65.        num=Adum.numc;
  66.        den=Adum.denc;
  67.        Adum=tf2ss(num,den);
  68.        A=Adum.a;
  69.        B=Adum.b;
  70.        C=Adum.c;
  71.        D=Adum.d;
  72.        nu=1;
  73.        ny=1;
  74.        iiu=1;
  75.    else
  76. // See if a,b,c, and d are compatible.
  77.        A=a;
  78.        B=b;
  79.        C=c;
  80.        D=d;
  81.        msg="";
  82.        msg=abcdchk(A,B,C,D);
  83.        if (msg != "") {
  84.            estr="IMPULSE: "+msg;
  85.            error(estr);
  86.        }
  87.        ny=c.nr;
  88.        nu=b.nc;
  89.        if (!exist(iu)) {
  90.            iiu=-1;
  91.        else
  92.            iiu=iu;
  93.        }
  94.    }
  95.  
  96. // Determine sample frequency
  97.    f=50*max(abs(eig(a).val)) / (2*pi);
  98. // Convert continuous model to discrete with 1/f as the sample time
  99.    Adum=c2d(A,B,1/f);
  100.    aa=Adum.phi;
  101.    bb=Adum.gamma;
  102.    cc=C;
  103.    dd=zeros(ny,nu);
  104.  
  105. // Set-up the time axis vector
  106.    t=[0:(ntim-1)/f:1/f]';
  107.  
  108.  
  109. // Plot the simulation output
  110.    pstart();
  111.  
  112.    if (iiu < 0) {
  113.        for (i in 1:ny) {
  114.             for (j in 1:nu) {
  115. // Simulate the response by using ltitr
  116.                  x=ltitr(aa,bb[;j],zeros(ntim,1),b[;j]);
  117.                  y=x*C[i;].';
  118.                  timeplot(i,j,f,ntim,y,"Impulse Response");
  119.                  pause;
  120.             }
  121.        }
  122.    else
  123.        for (i in 1:ny) {
  124.             x=ltitr(aa,bb[;iiu],zeros(ntim,1),b[;iiu]);
  125.             y=x*C[i;].';
  126.             timeplot(i,iiu,f,ntim,y,"Impulse Response");
  127.             pause;
  128.        }
  129.    }
  130.  
  131. };
  132.