home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / COVAR.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.8 KB  |  61 lines
  1. function [p,q] = covar(a,b,c,d,w)
  2. %COVAR  Covariance response of continuous system to white noise.
  3. %    [P,Q] = COVAR(A,B,C,D,W) computes the covariance response of the 
  4. %    continuous state-space system (A,B,C,D) to Gaussian white noise
  5. %    inputs with intensity W,
  6. %
  7. %        E[w(t)w(tau)'] = W delta(t-tau) 
  8. %
  9. %    where delta(t) is the dirac delta.  P and Q are the output and 
  10. %    state covariance responses:
  11. %
  12. %        P = E[yy'];  Q = E[xx'];
  13. %
  14. %    P = COVAR(NUM,DEN,W) computes the output covariance response of 
  15. %    the polynomial transfer function system.  W is the intensity of 
  16. %    the input noise.
  17. %
  18. %    Warning: Unstable systems or systems with a non-zero D matrix have
  19. %    infinite covariance response.
  20. %
  21. %    See also: DCOVAR,LYAP and DLYAP.
  22.  
  23. %    Clay M. Thompson  7-3-90
  24. %       Revised by Wes Wang 8-5-92
  25. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  26.  
  27. error(nargchk(3,5,nargin));
  28. if ~((nargin==3)|(nargin==5)), error('Wrong number of input arguments.'); end
  29.  
  30. % --- Determine which syntax is being used ---
  31. if (nargin == 3) % T.F. syntax
  32.   [num,den] = tfchk(a,b);
  33.   w = c;
  34.   [a,b,c,d] = tf2ss(num,den);
  35. end
  36.  
  37. if (nargin==5), error(abcdchk(a,b,c,d)); end
  38.  
  39. q = lyap(a,b*w*b');
  40. p = c*q*c';
  41.  
  42. % Systems with non-zero D matrix have infinite output covariance.
  43. if any(abs(d(:))>eps),
  44.   disp('Warning: Systems with non-zero D matrix have infinite output covariance.')
  45.   pinf = inf*ones(length(p));
  46.   if ~isempty(p), p((abs(d*w*d')>eps)) = pinf(abs(d*w*d')>eps); end
  47. end
  48.  
  49. % A valid covariance must be positive semi-definite.
  50. if min(real(eig(q))) < -eps,
  51.   if max(real(eig(a))) > 0
  52.     disp('Warning: Unstable system. Returning Infinity.')
  53.   else
  54.     disp('Warning: Invalid covariance - not positive semi-definite.  Returning infinity.')
  55.   end
  56.   q = inf*ones(length(q)); p = inf*ones(length(p));
  57. end
  58.  
  59.  
  60.  
  61.