function ms = fms(w,f);
%FMS Measure of skewness of MIMO frequency response.
%
%       MS = FMS(W,F) calculates a measure of skewness of an MVFR matrix, F. 
%       W is the associated frequency response matrix.
%
%       Any square matrix, F, has the Schur decomposition:
%
%                       *
%          F = S(D + T)S
%
%       where S is unitary, D is diagonal, and T is strictly triangular.
%
%          MS(F) = ||T||/||F||   where ||.|| denotes the Frobenius norm.
%
%       The result is usually plotted against frequency: SEMILOGX(W,MS).

%       J-M. Boyle, 8 Sep 1987. Revised by J.M.Maciejowski, 3 Jan 1988.
% Copyright (c) 1987 by GEC Engineering Research Centre & Cambridge Control Ltd.

ms = zeros(1,length(w));           % Set-up a vector to hold the result

for i = 1:length(w);
  ftemp = fgetf(w,f,i);           % Get each component matrix
  [s,tu] = schur(ftemp);          % Compute the Schur form of the
                                  % component matrix.

  if ~any(any(imag(ftemp)))       % If the Schur form is real at any freq.
    [s,tu] = rsf2csf(s,tu);       % use RSF2CSF to convert to a complex
                                  % Schur form.
  end

  t = tu - diag(diag(tu));        %  T is the strictly upper triangular part
                                  %  of the Schur matrix Tu.

  ms(i) = norm(t,'fro')/norm(ftemp,'fro');  % MS is the Frobenius norm of T
end                               % divided by the Frobenius norm of F.