function cltfm = ffbf(w,g,h,k);
%FFBF Closed-loop frequency response with arbitrary feedback (MIMO).
%
%       CLTFM = FFBF(W,G,H,K) calculates the closed loop MVFR matrix.
%                CLTFM = inv(I + GKH).GK
%
%       The calculations are performed for the system shown in the FFB section
%       of the MFD Reference manual.
%
%       The input arguments are:
%          G the plant
%          H a non-unity feedback matrix and optionally
%          K a pre-compensator. G, H and K are MVFR matrices.
%
%       CLTFM is returned as an MVFR matrix.

%       J-M. Boyle 8th September 1987
% Copyright (c) 1987 by GEC Engineering Research Centre & Cambridge Control Ltd.

nargs = nargin;
error(nargchk(3,4,nargs));       % Check the number of input arguments

[mg,ng] = fsize(w,g);
lw = length(w);
[mh,nh]=fsize(w,h);

if nargs == 3              % IF 3 input arguments, the compensator K
                            % is an identity matrix.
  if (ng~=mh) | ( nh~=mg)
     error('G*H and H*G must both be square (K = I).')
  end
  k = eyef(w,ng);
  mk = ng;
  nk = ng;
else
  [mk,nk]=fsize(w,k);
  if nk~=mh
     error('Rows of H not equal to columns of K.')
  end
  if mk~=ng
     error('Rows of K not equal to columns of G.')
  end
  if nh~=mg
     error('Rows of G not equal to columns of H.')
  end
end

id = eye(mg);
cltfm =  zeros(mg*lw,nk);
indg = 1:mg;
indk = 1:mk;
indh = 1:mh;

for i=1:lw;
  gm = g(indg+mg*(i-1),:);
  km = k(indk+mk*(i-1),:);
  hm = h(indh+mh*(i-1),:);
  cltfm(indg+mg*(i-1),:) = (id + gm*km*hm)\gm*km;
                % CLTFM = inv(I + GKH) * GK.
end