function [upper,lower]=fadj(w,f)
%FADJ   Estimates the order of a pseudo-decoupling pre-compensator
%       [UPPER,LOWER]=FADJ(W,F) estimates the difference in the
%       orders of the pre-compensator elements required at the
%       upper and lower frequencies to diagonalize the system.
%       The system needs at least as many columns as rows.
%
%       The matrix UPPER contains the relative orders of
%       the compensator elements for the upper frequencies while
%       LOWER contains the relative orders for the lower frequencies.
%
%       W needs to be monotonic and for accurate estimates needs to
%       extend above and below the system's non-zero and finite
%       poles and zeros.
%
%       If the numbers in UPPER and LOWER are not approximately integers
%       then the frequency range does not extend high enough or low enough
%       respectively.
%
%       See FPSEUDO

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


[m,n]=fsize(w,f);
if m>n
   error('System has more row then columns')
end
lw=length(w);
if lw<=1
   error('Error not enough frequency points.')
end
k=1:m;         % vector of rows of each matrix in F

% Calculate upper
if w(lw)<=w(lw-1)
   error('Frequency points not monotonically increasing.')
end
fm=f((lw-1)*m+k,:);
u2=log10(abs(pinv(fm)));   % abs inverse of last freq. matrix
fm=f((lw-2)*m+k,:);
u1=log10(abs(pinv(fm)));   % abs inverse of next to last
upper=(u2-u1)./(log10(w(lw)/w(lw-1)));  % divide by log freq
di=diag(upper);    % normalize the columns so that
		   % the diagonal elements are zero
for i=1:m
    upper(:,i)=upper(:,i)-di(i);
end

% Calculate lower
if w(2)<=w(1)
   error('Frequency points not monotonically increasing.')
end
fm=f(m+k,:);
u2=log10(abs(pinv(fm)));   % abs inverse of second freq. matrix
fm=f(k,:);
u1=log10(abs(pinv(fm)));   % abs inverse of first
lower=(u2-u1)./(log10(w(1)/w(2)));  % divide by log freq
di=diag(lower);    % normalize the columns so that
		   % the diagonal elements are zero
for i=1:m
    lower(:,i)=lower(:,i)-di(i);
end