home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 4 / DATAFILE_PDCD4.iso / languages / rlab1_23a / CTB / impulse < prev    next >
Text File  |  1995-11-14  |  3KB  |  119 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.  
  43. impulse = function(a,b,c,d,ntim,iu)
  44. {
  45.    local(nargs,Adum,num,den,nu,ny,msg,estr,aa,bb,cc,dd,Y,i,j,f,iiu)
  46.  
  47. // Count number of input arguments
  48.    nargs=0;
  49.    if (exist(a)) {nargs=nargs+1;}
  50.    if (exist(b)) {nargs=nargs+1;}
  51.    if (exist(c)) {nargs=nargs+1;}
  52.    if (exist(d)) {nargs=nargs+1;}
  53.    if (exist(ntim)) {nargs=nargs+1;}
  54.    if (exist(iu)) {nargs=nargs+1;}
  55.  
  56.    if ( (nargs <= 2) || (nargs > 6) ) {
  57.        error("impulse: Error in number of input arguments.");
  58.    }
  59.  
  60. // Convert to S-S if in TF form
  61.    if (nargs < 5) {
  62.        Adum=tfchk(a,b);
  63.        num=Adum.numc;
  64.        den=Adum.denc;
  65.        Adum=tf2ss(num,den);
  66.        a=Adum.a;
  67.        b=Adum.b;
  68.        c=Adum.c;
  69.        d=Adum.d;
  70.        nu=1;
  71.        ny=1;
  72.        iiu=1;
  73.    else
  74. // See if a,b,c, and d are compatible.
  75.        msg="";
  76.        msg=abcdchk(a,b,c,d);
  77.        if (msg != "") {
  78.            estr="impulse: "+msg;
  79.            error(estr);
  80.        }
  81.        ny=c.nr;
  82.        nu=b.nc;
  83.        if (!exist(iu)) {
  84.            iiu=-1;
  85.        else
  86.            iiu=iu;
  87.        }
  88.    }
  89.  
  90. // Determine sample frequency
  91.    f=50*max(abs(eig(a).val)) / (2*pi);
  92. // Convert continuous model to discrete with 1/f as the sample time
  93.    Adum=c2d(a,b,1/f);
  94.    aa=Adum.phi;
  95.    bb=Adum.gamma;
  96.    cc=c;
  97.    dd=zeros(ny,nu);
  98.  
  99. // Simulate the system by calling simdata
  100.    Y=simdata(aa,bb,cc,dd,ntim);
  101.  
  102. // Plot the impulse response
  103.    pstart();
  104.  
  105.    if (iiu < 0) {
  106.        for (i in 1:ny) {
  107.             for (j in 1:nu) {
  108.                  timeplot(i,j,f,ntim,Y,"Impulse Response");
  109.                  pause;
  110.             }
  111.        }
  112.    else
  113.        for (i in 1:ny) {
  114.             timeplot(i,iiu,f,ntim,Y,"Impulse Response");
  115.        }
  116.    }
  117.  
  118. };
  119.