function [lyapunov,t] = computelyapunovspectrum(data,nMaxDt,dt,nNearest,radius)
% COMPUTELYAPUNOVSPECTRUM Calculates Lyapunov spectrum for experimental data.
%   COMPUTELYAPUNOVSPECTRUM(DATA,NMAXDT) returns the lyapunov exponent
%   spectrum computed for NMAXDT samples foward in time. 
%
%   COMPUTELYAPUNOVSPECTRUM(DATA,NMAXDT,DT) where DT is the sampling rate. It 
%   defaults to unity if not specified.
%   
%   COMPUTELYAPUNOVSPECTRUM(DATA,NMAXDT,DT,NNEAREST) where NNEAREST is an 
%   optional parameter that may be used to increase the number of nearest
%   neighbors used in the calculation of the state-transition matrix.
%   Default is the minimum number of states.
%
%   COMPUTELYAPUNOVSPECTRUM(DATA,NMAXDT,DT,NNEAREST,RADIUS) where RADIUS
%   is an optional parameter that specifies the minimum and maximum distances 
%   that nearest neighbors are allowed to be from the central point, RADIUS 
%   = [MINRADIUS MAXRADIUS]. It defaults to [0 Inf] and the minimum number of 
%   nearest neighbors to construct the state-transition matrix are found.
%
%   See the note in STATEDISTANCE about the size of the data sets that may be
%   analyzed. If larger data sets are desired, it will be necessary to modify
%   this function to perhaps use FINDNEARESTNEIGHBORS instead.
%
%   Input
%       DATA     | m x n | double | m is # samples and n is # states
%       NMAXDT   | 1 x 1 | double | integer # of samples to follow
%       DT       | 1 x 1 | double | Time between samples, Default: 1
%       NNEAREST | 1 x 1 | double | See [1], Default: n 
%       RADIUS   | 1 x 2 | double | [MINRADIUS MAXRADIUS], Default: [0 Inf]
%
%       [1] NNEAREST is the number of nearest neighbors to search for within the
%           bounds specified by RADIUS. NNEAREST must be greater than the
%           number of states, n. Additionally, if NNEAREST cannot be found
%           within RADIUS, it expands only the upper bound, MAXRADIUS, until
%           NNEAREST are found.
%
%   Output
%       LYAPUNOV | (m-NMAXDT) x n | double | Lyapunov exponent spectrum estimation
%       T        | (m-NMAXDT) x 1 | double | Time
%   
%   Example     
%       lorenzex; % Run example script
%
%   Reference
%     [1] Banbrook, M., Ushaw, G., & McLaughlin, S., 1997, "How to extract 
%         Lyapunov exponents from short and noisy time series," IEEE 
%         Transactions on Signal Processing, Vol. 45, No. 7, pp. 1378-1382.
%
%   See also LORENZEX, COMPUTEDIVERGENCE, ANALYTICALLYAPUNOV.  

% Created by: Jesse Norris 
%             Kevin P. Granata Musculoskeletal Biomechanics Laboratory
%             http://www.biomechanics.esm.vt.edu
%             http://jesse.musculoskeletaldynamics.com
% Last Revision: 4 October 2007

% Check that nMaxDt is an integer
  if (fix(nMaxDt) ~= nMaxDt)
      error('Input argument nMaxDt should be an integer.');
  end;
  
% Handle the optional input arguments radius and nNearest
  if ~exist('radius','var')   || isempty(radius);   radius = [0 inf]; end;
  if ~exist('nNearest','var') || isempty(nNearest); nNearest = size(data,2); end; 
  if ~exist('dt','var'); dt = 1; end; 

% Determine the number of samples to skip to avoid sampling nearest neighbors on
% same cycle. To ensure nSkip is sufficiently large, skip +/- first minimum of 
% autoCorrelation, i.e. +/- nSkip.
  autoCorrelation  = xcorr(data(:,1));
  autoCorrelation  = autoCorrelation(ceil(length(autoCorrelation)/2):end);
  nSkip = find(diff(autoCorrelation)>0,1); % Find first min of autoCorrelation

% Calculate state distance between every sample instant
  stateDistance = statedistance(data,nSkip,nMaxDt);
  
% Estimate the state transition matrix by tracking nearest neighbors.
  nSamples = size(data,1);    % Total # of time samples
  nValid   = nSamples-nMaxDt; % # of samples that may be tracked forward NMAXDT
  verticalIndex = [0 cumsum((nValid-nSkip-2):-1:1)]; % Used to extract appropriate data from the stateDistance vector  
  lyapunov   = zeros(nValid,size(data,2)); % Pre-allocate
  cumulative = zeros(1,size(data,2));
  Q = eye(size(data,2));      % Initialize the Q-matrix
  %distanceMeanStd = zeros(nValid,2);  
  
  for i = 1:nValid
    % Select appropriate row
      k = (i-1)*(nValid-nSkip-i/2)+(1:(nValid-nSkip-i));
      distance = [zeros(1,nValid-length(k)) stateDistance(k)'];
      
    % Append the transposed column from the upper triangular matrix
      l = (i-nSkip-1)+verticalIndex(1:(i-nSkip-1));
      distance = distance+[stateDistance(l)' zeros(1,nValid-length(l))];
      
    % Select all instances where state distance within desired bound
      I = find((distance > radius(1)) & (distance < radius(2)));
      
    % Handle if too many data points were selected
      if (length(I) > nNearest)
          [junk,J] = sort(distance(I));     % Sort increasing distances
          I = I(J(1:nNearest));             % Return the desired number of states to compute state transition matrix          
          
    % Handle if too few data points were selected    
      elseif (length(I) < size(data,2) || length(I) < nNearest)
          warning('computelyapunovspectrum:radius',['Desired/suitable number ',...
                  'of neighbors were not found within desired distances']);
          I = find(distance > radius(1)); % Find nearest neighbors outside of radius(1)
          [junk,J] = sort(distance(I));   % Sort increasing distances
          I = I(J(1:nNearest));           % Return the desired number of states to compute state transition matrix                    
      end;
      clear junk J
      
    % Determine evolution of nearest neighbors after nMaxDt time increments
    % Renormalize using QR decomposition to estimate expansion along orthogonal
    % directions.
    
    % Keep track of the distances used 
    %  distanceMeanStd(i,:) = [mean(distance(I)) std(distance(I))];
    
    % Calculate inital set of vectors, forming a matrix, X, about data(i,:)
      neighbor = data(I,:)';
      X0 = neighbor - repmat(data(i,:)',1,size(neighbor,2));       
      
    % Advance forward nMaxDt time-steps
      neighbor = data(I+nMaxDt,:)';
      X1 = neighbor - repmat(data(i+nMaxDt,:)',1,size(neighbor,2));            
      
    % Estimate state-transition matrix from i to i+nMaxDt
      phi = X1/X0;            
      
    % Peform QR (Gram-Schmidt) orthonormalization Q=orthonorm vector, R=magnitude
    % Note that phi is multiplied by the previous Q to track exapansion of the 
    % rotating frame.
      [Q,R] = qr(phi*Q);
      
    % Vector length growth = abs(diag(R))         
      cumulative = cumulative + log(abs(diag(R)))'/nMaxDt;            
      
    % Estimation of Lyapunov exponents
      lyapunov(i,:) = cumulative/(i*dt);              
  end;  
  
  % Output the time as well
    t = dt*(0:(i-1))';    
    
end % END computelyapunovspectrum