function eigen = computefloquet(data)
% COMPUTEFLOQUET Calculate Floquet multipliers from experimental data.
%   EIGEN = COMPUTEFLOQUET(DATA) returns the estimated eigenvalues of a linear 
%   map in descending magnitude. DATA is a (m x n) matrix of class double where 
%   <m> is the number of events, ordered sequentially, and <n> is the number of 
%   states.
% 
%   Input
%       DATA  | m x n | double |
%
%   Output
%       EIGEN | n x 1 | double | Eigenvalues
%   
%   Example        
%       computefloquetex; % Run the example script
%
%   Reference
%     [1] Hurmuzlu, Y. and Basdogan, C., 1994, "On the measurement of dynamic
%         stability of human locomotion," J Biomech.Eng, Vol. 116, No. 1, 
%         pp. 30-36.
%
%   See also COMPUTEFLOQUETEX.

% Created by: Jesse Norris
%             Kevin P. Granata Musculoskeletal Biomechanics Laboratory
%             http://jesse.musculoskeletaldynamics.com
% Last Revision: 25 May 2007
  
% Establish states at event i and at i+1
  si   = data(1:end-1,:); % States at i
  sip1 = data(2:end,:);   % States at i+1    

  [nEvents,nStates] = size(data); % Number of events and number of states 

% sip1(:,j) = [si(:,:) ones]*[J(j,1:nStates) C(j)]';    
% 
% [J(j,1:nStates) C(j)]' = [si(:,:) ones]\sip1(:,j); % Least-squares fit
%
% Thus, the coefficients need to be transposed before they are the Jacobian.
%
  coeff = [si ones(nEvents-1,1)]\sip1;
  coeff = transpose(coeff);    
  J     = coeff(:,1:nStates);  % Jacobian (remove the column of constants)
  eigen = eig(J);              % Eigenvalues of Jacobian (Floquet multipliers)    

  [junk,i] = sort(abs(eigen),1,'descend'); % Sort by magnitude
  eigen    = eigen(i);             

% Check the linear fit by examining the assumed equilibrium point
% mean(sip1,1)'-J*mean(si,1)'-fit(:,end)

end % END computefloquet