mv2fr.m

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

function cout=mv2fr(a,b,c,d,w,iu)
%MV2FR  Frequency response of MIMO system
%       MV2FR(A,B,C,D,W)  calculates the MVFR matrix of the system:
%               .
%               x = Ax + Bu                          -1
%               y = Cx + Du             G(s) = C(sI-A) B + D
%
%       Vector W contains the frequencies, in radians, at which the 
%       frequency response is to be evaluated. Each component matrix 
%       of the resulting MVFR MVFR matrix has as many columns as 
%       inputs and as many rows as outputs.
%
%       MV2FR(NUM,COMDEN,W) calculates the MVFR matrix from the
%       transfer function description  G(s) = NUM(s)/COMDEN(s)
%       where NUM is multivariable matrix of polynomial coefficients and
%       COMDEN is the vector containing Common Denominator polynomial
%       coefficients in descending powers of s.
%
%       MV2FR(A,B,C,D,W,iu) and MV2FR(NUM,COMDEN,W,iu) are provided for
%       compatibility 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.
%
%       The CONTROL TOOLBOX's BODE function is provided by
%       [MAG,PH]=R2P(MV2FR(A,B,C,D,W,iu))

%       Dr M.P. Ford 4th August 1987
% Copyright (c) 1987 by GEC Engineering Research Centre & Cambridge Control Ltd
% MRN0019

nargs = nargin;
error(nargchk(3,6,nargs));
if nargs <= 4
  % 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 = sqrt(-1)*c;
  k=1:lb;
  if nargs == 4    % just one input iu=d
    if max(size(d))~=1
       error('Input specified is not a scalar');
    end
    if (d<=0)|(d>ni)
       error('Input specified does not exist in this system');
    end
    a=a(:,(d-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~=4
     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 = w * sqrt(-1);

  % 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==6  % 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~=6
     cout=shpf(cout,md,nd);
  end
end  % else