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