covar.m

数字通信第四版原书的例程

function [p,q] = covar(a,b,c,d,w)
%COVAR  Covariance response of continuous system to white noise.
%	[P,Q] = COVAR(A,B,C,D,W) computes the covariance response of the 
%	continuous state-space system (A,B,C,D) to Gaussian white noise
%	inputs with intensity W,
%
%		E[w(t)w(tau)'] = W delta(t-tau) 
%
%	where delta(t) is the dirac delta.  P and Q are the output and 
%	state covariance responses:
%
%		P = E[yy'];  Q = E[xx'];
%
%	P = COVAR(NUM,DEN,W) computes the output covariance response of 
%	the polynomial transfer function system.  W is the intensity of 
%	the input noise.
%
%	Warning: Unstable systems or systems with a non-zero D matrix have
%	infinite covariance response.
%
%	See also: DCOVAR,LYAP and DLYAP.

%	Clay M. Thompson  7-3-90
%       Revised by Wes Wang 8-5-92
%	Copyright (c) 1986-93 by the MathWorks, Inc.

error(nargchk(3,5,nargin));
if ~((nargin==3)|(nargin==5)), error('Wrong number of input arguments.'); end

% --- Determine which syntax is being used ---
if (nargin == 3) % T.F. syntax
  [num,den] = tfchk(a,b);
  w = c;
  [a,b,c,d] = tf2ss(num,den);
end

if (nargin==5), error(abcdchk(a,b,c,d)); end

q = lyap(a,b*w*b');
p = c*q*c';

% Systems with non-zero D matrix have infinite output covariance.
if any(abs(d(:))>eps),
  disp('Warning: Systems with non-zero D matrix have infinite output covariance.')
  pinf = inf*ones(length(p));
  if ~isempty(p), p((abs(d*w*d')>eps)) = pinf(abs(d*w*d')>eps); end
end

% A valid covariance must be positive semi-definite.
if min(real(eig(q))) < -eps,
  if max(real(eig(a))) > 0
    disp('Warning: Unstable system. Returning Infinity.')
  else
    disp('Warning: Invalid covariance - not positive semi-definite.  Returning infinity.')
  end
  q = inf*ones(length(q)); p = inf*ones(length(p));
end