function phase = fpp(w,f);
%FPP    Calculates principal phases.
%       FPP(W,F) returns the principal phases (in degrees) of
%       F, an MVFR matrix, as a set of column vectors.
%       W is the associated frequency vector.
%
%       The principal phases are only available for square systems.
%       Algorithm:   Let G have SVD  = Y Sigma U'
%                                    = YU' U Sigma U'
%                where  Thi := YU' := Phase matrix (of G)
%
%           The eigenvalue decomposition of Thi is
%                              YU' = PQP'
%                  where        Q  = diag(exp(sqrt(-1)*p(i))
%        The angles p(i) are  principal phases of G
%
%       (The principal phases may need sorting using UNWRAP and CSORT.)

%       J-M Boyle 7th September 1987
% Copyright (c) 1987,1993 by GEC Engineering Research Centre and 
% Cambridge Control Ltd. All Rights Reserved.
% History:
%   Reference to `fixphase' replaced by one to `unwrap', 29.5.93, JMM.

[m,n]=fsize(w,f);
if m~=n
   error('System is not square');
end
[yf, sf, uf] = fsvd(w, f);
pm = fmulf(w, yf, ftrn(w,uf));
[mag,phase] = r2p(feig(w, pm));