home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 4 / DATAFILE_PDCD4.iso / languages / rlab1_23a / CTB / dcovar < prev    next >
Encoding:
Text File  |  1995-11-15  |  3.9 KB  |  149 lines
  1. //-----------------------------------------------------------------------
  2. //
  3. // dcovar
  4. //
  5. // Syntax: R=dcovar(A,B,C,D,W)  or  R=dcovar(NUM,DEN,W)
  6. //
  7. // This routine computes the covariance response of a discrete
  8. // system to white noise input. Calling the routine as,
  9. //
  10. //    R=dcovar(A,B,C,D,W)
  11. //
  12. // computes the covariance response of the discrete state-space
  13. // system defined by A,B,C,D to Gaussian white noise input with
  14. // intensity W given by the following,
  15. //  
  16. //              E[w(k)w(tau)'] = W delta(k-tau)
  17. //  
  18. // where delta(k,tau) is the kronecker delta and E is the expectation
  19. // operator.
  20. //
  21. // The routine returns both the state covariance matrix X and the
  22. // output covariance matrix Y given by the following,
  23. //
  24. //     X=E[xx']
  25. //     Y=E[yy']
  26. //
  27. // where x is the state vector and y is the output vector.
  28. //
  29. // This routine can also be used to compute the covariance response
  30. // of a transfer function type system by calling it as,
  31. //
  32. //     R=covar(NUM,DEN,W)
  33. //
  34. // where NUM is a vector containing the coefficients of the numerator
  35. // transfer function polynomial, DEN contains the coefficients of the
  36. // denominator transfer function polynomial, and W is the intensity
  37. // of the white noise input.
  38. //
  39. // The results are returned in a list. For example,  R=covar(A,B,C,D,W)
  40. // produces,
  41. //
  42. //  R.Y = output covariance matrix (Square of Output RMS)
  43. //  R.X = state covariance matrix
  44. //
  45. // Ref: Skelton, R. "Dynamic Systems Control Linear System Analysis and
  46. //      Synthesis," John Wiley and Sons, 1988.
  47. //
  48. // Copyright (C), by Jeffrey B. Layton, 1994
  49. // Version JBL 940918
  50. //-----------------------------------------------------------------------
  51.  
  52. rfile abcdchk
  53. rfile tfchk
  54. rfile tf2ss
  55. rfile dlyap
  56.  
  57. dcovar = function(a,b,c,d,w)
  58. {
  59.    local(nargs,msg,estr,num,den,X,Y,A,B,C,D,W,Dum)
  60.  
  61. // Count number of arguments
  62.    nargs=0;
  63.    if (exist(a)) { nargs=nargs+1; }
  64.    if (exist(b)) { nargs=nargs+1; }
  65.    if (exist(c)) { nargs=nargs+1; }
  66.    if (exist(d)) { nargs=nargs+1; }
  67.    if (exist(w)) { nargs=nargs+1; }
  68.  
  69. // Error trap for wrong number of inputs.
  70.    if ( nargs != 3) {
  71.         if ( nargs != 5) {
  72.             error("DCOVAR: Wrong number of input arguments.");
  73.         }
  74.    }
  75.  
  76. // Error Traps for matrix compatibility
  77.    if (a.nr != b.nr) {
  78.        error("DCOVAR: A and B have different number of rows.");
  79.    }
  80.  
  81.    if (w.nr != w.nc) {
  82.        error("DCOVAR: W needs to be square.");
  83.    }
  84.  
  85.    if (w.nr != b.nr) {
  86.        error("DCOVAR: rows(W) not equal to rows(B).");
  87.    }
  88.  
  89.    if (w.nc != b.nr) {
  90.        error("DCOVAR: cols(W) not equal to rows(B).");
  91.    }
  92.  
  93. // Transfer Function Representation
  94.    if (nargs == 3) {   
  95.        Dum=tfchk(a,b);
  96.        num=Dum.numc;
  97.        den=Dum.denc;
  98.        W=c;
  99. // Convert to state-space form (using MATLAB compatible conversion)
  100.        Dum=tf2ss(num,den,2);
  101.        A=Dum.a;
  102.        B=Dum.b;
  103.        C=Dum.c;
  104.        D=Dum.d;
  105.    else
  106. // State-Space Representation
  107.        if (nargs == 5) {
  108.            A=a;
  109.            B=b;
  110.            C=c;
  111.            D=d;
  112.            W=w;
  113.  
  114.            msg="";
  115.            msg=abcdchk(A,B,C,D);
  116.            if (msg != "") {
  117.                estr="DCOVAR: "+msg;
  118.                error(estr);
  119.            }
  120.        else
  121.            error("DCOVAR: Wrong number of input arguments");
  122.        }
  123.    }
  124.  
  125.  
  126. // Find state covariance by solving discrete lyapunov
  127. // Then find the output covariance matrix Y.
  128.    X=dlyap(a,B*W*B');
  129.    Y=C*X*C'+D*W*D';
  130.  
  131. // A valid covariance must be positive semi-definite.
  132.    if (min(real(eig(X).val)) < -epsilon()) {
  133.    printf("%s","Warning: Invalid covariance - not positive semi-definite. Returning infinity.\n");
  134.        X=inf()*ones(X.nr,X.nc);
  135.        Y=inf()*ones(Y.nr,Y.nc);
  136.        return << X=X; Y=Y >>
  137.    }
  138.  
  139. // The system must be stable for a valid covariance
  140.    if (any(real(eig(a).val) > 0) ) {
  141.        printf("%s","Warning: Unstable system. Returning Infinity.\n");
  142.        X=inf()*ones(length(X));
  143.        Y=inf()*ones(length(Y));
  144.        return << X=X; Y=Y >>
  145.    }
  146.  
  147.    return << X=X; Y=Y >>
  148. };
  149.