dmv2fr.m

控制系统计算机辅助设计——MATLAB语言与应用(源代码）

function cout=dmv2fr(a,b,c,d,w,dt,iu)
%DMV2FR Frequency response of MIMO discrete time system
%	DMV2FR(A,B,C,D,W,DT)  calculates the MVFR matrix of the system:
%
%		x(t+dt) = Ax(t) + Bu(t) 		    -1
%		   y(t) = Cx(t) + Du(t)        G(z) = C(zI-A) B + D
%
%	at the points z = exp(sqrt(-1)*W*DT). Usually max(W) < pi/DT.
%	Vector W contains the frequencies, in rad/sec, at which the frequency
%	response is to be evaluated. Each component matrix of the resulting 
%	MVFR matrix has as many columns as inputs and as many rows as outputs.
%
%	DMV2FR(NUM,COMDEN,W,DT) calculates the MVFR matrix from the
%	transfer function description  G(z) = NUM(z)/COMDEN(z)
%	where NUM is the multivariable matrix of polynomial coefficients and
%	COMDEN is the vector containing Common Denominator polynomial
%	coefficients in descending powers of z.
%
%	DMV2FR(A,B,C,D,W,DT,iu) and DMV2FR(NUM,COMDEN,W,DT,iu) are provided for
%	compatiblity with the CONTROL TOOLBOX's NYQUIST.
%	They produce an MVFR matrix for a single input, iu, and the results
%	are returned as one row per frequency.

%	P. Phaal, November 1987
% Copyright (c) 1987 by GEC Engineering Research Centre & Cambridge Control Ltd
%       MRN0019

nargs = nargin;
error(nargchk(4,7,nargs));
if nargs <= 5
  % It is in transfer function form.  Do directly, using Horner's method
  % of polynomial evaluation at the frequency points, for each row in
  % the numerator.  Then divide by the denominator.
  [ma,na] = size(a);
  lb = length(b);
  if na > lb
    if rem(na,lb)~=0
      error('Numerator matrix not consistent with Common Denominator.');
    end
    ni=na/lb;	      % number of inputs of system
  else
    ni=1;
    a=[zeros(ma,lb-na),a];  % pad out numerator matrix
  end
  c=c(:);
  c = exp(d*sqrt(-1)*c);
  k=1:lb;
  if nargs == 5    % just one input iu=w
    if max(size(w))~=1
       error('Input specified is not a scalar');
    end
    if (w<=0)|(w>ni)
       error('Input specified does not exist in this system');
    end
    a=a(:,(w-1)*lb+k);
    ni=1;
  end
  lc=length(c);
  cout=zeros(lc,ni*ma);
  for j=1:ni	    % for each column
    for i=1:ma	     % for each row
      cout(:,i+ma*(j-1)) = polyval(a(i,(j-1)*lb+k),c);
    end
  end
  dcout = polyval(b,c);
  for i=1:lc;
    cout(i,:)=cout(i,:)./dcout(i);
  end
  if nargs~=5
     cout=shpf(cout,ma,ni);   % reshape cout into MVFR
  end
else	       %  state space input
  error(abcdchk(a,b,c,d));
  [md,nd]=  size(d);
  [no,ns] = size(c);
  nw = max(size(w));
  index=(0:(nw-1)).*md;  % indices for first input
  w = exp(dt*sqrt(-1)*w);

  % Uncomment these line to Balance A
  %[t,a] = balance(a);
  %b = t \ b;
  %c = c * t;

  [p,a] = hess(a);		% Reduce A to Hessenberg form:
  %	 Apply similarity transformations from Hessenberg
  %	 reduction to B and C:
  b = p' * b;
  c = c * p;
  if nargs==7  % calculate one input only
     if max(size(iu))~=1
	error('Input specified is not a scalar.');
     end
     if (iu<=0)|(iu>nd)
	error('Input specified does not exist in this system.');
     end
     b=b(:,iu);
     d=d(:,iu);
     nd=1;
  end
  cout=zeros(nw,md*nd);
  k=1:md;
  for iu=1:nd  % loop through inputs
  g = ltifr(a,b(:,iu),w);
  g = c * g + diag(d(:,iu)) * ones(no,nw);
  cout(:,k+md*(iu-1))=g.';
  end	% for iu=iufirst:iulast
  if nargs~=7
     cout=shpf(cout,md,nd);
  end
end  % else