function [data,principal,sigma] = principalcomponents(data,mse)
% PRINCIPALCOMPONENTS Perform principal components analysis.
%   PRINCIPALCOMPONENTS(DATA) performs principal components analysis to 
%   determine if the data can be represented in a lower-dimensional orthogonal
%   frame that still accounts for 95% of the mean-square-error. Returns the 
%   data in the lower-dimensional frame. 
%   
%   PRINCIPALCOMPONENTS(DATA,MSE) choose components to meet MSE, expressed as a 
%   fraction, i.e. 0.90, instead of the default 95%. If MSE is 1.0, all 
%   components are returned and used in reducing data set.
%
%   [DATA,PRINCIPAL] = PRINCIPALCOMPONENTS(...) returns the principal component 
%   vectors.
%
%   [DATA,PRINCIPAL,SIGMA] = PRINCIPALCOMPONENTS(...) returns the percent
%   variance accounted for by inclusion of corresponding PRINCIPAL vector.
%
%   This function mean-zeros the input, DATA, but it does not scale DATA by 
%   standard deviation. 
%
%   Input
%       DATA  | m x n | double | m is # of samples and n is # of channels
%       MSE   | 1 x 1 | double | mean-square-error to satisfy
%
%   Output
%       DATA      | m x n-reduced | double | data in the reduced frame
%       PRINCIPAL | n x n-reduced | double | eigenvectors rotate to this frame
%       SIGMA     | n-reduced x 1 | double | MSE accounted for 
%   
%   Example        
%       t     = (0:0.01:30)';        
%       data  = [sin(2*pi*t) sin(2*pi*t) zeros(size(t))];     
%       data  = data + 0.1*randn(size(data)); % Add noise
%       [reduced,principal,sigma] = principalcomponents(data,1);
%       figure('Name','Three-dimensional');
%       plot3(data(:,1),data(:,2),data(:,3));
%       figure('Name','Raw & Reduced-data');
%       subplot(3,2,1); plot(t,data(:,1));
%       subplot(3,2,3); plot(t,data(:,2));
%       subplot(3,2,5); plot(t,data(:,3));
%       subplot(3,2,2); plot(t,reduced(:,1));
%       subplot(3,2,4); plot(t,reduced(:,2));
%       subplot(3,2,6); plot(t,reduced(:,3)); 
%
%   Reference
%     [1] Daffertshofer, A, Lamoth, CJ, Meijer, OG, and Beek, PJ, "PCA in
%         studying coordination and variability: a tutorial," Clin.Biomech., 
%         Vol. 19, No. 4, pp. 415-428, 2004.
%
%   See also COV,PRINCOMP.

% Created by: Kevin P. Granata Musculoskeletal Biomechanics Laboratory
%             http://jesse.musculoskeletaldynamics.com
% Last Revision: 23 July 2007

% I wrote this before I realized Matlab statistics include PRINCOMP - oops.

% Handle input arguments
  if (nargin == 1)
      mse = 0.95; % Default of 95% variance explained
  end;
 
% Number of samples & channels 
  [nSamples nChannels] = size(data);  
  
% Remove mean from each column and normalize by standard deviation
  data = data-repmat(sum(data,1)/nSamples,nSamples,1);

% Compute covariance matrix (nChannels x nChannels)
  R = (data'*data)/(nChannels-1);

% Eigenvectors: Principal components
% Eigenvalues:  Mean square error (MSE) reduction of each component
  [vector,eigen] = eig(R);

  [eigen,I] = sort(diag(eigen),'descend'); % Sort eigenvalues, greatest to least
  vector    = vector(:,I);                 % Re-organize eigenvectors
  
% Normalize eigen by trace(R) == sum(eigen)
  eigen = eigen/sum(eigen);        

% Find index for components which result in mean-square-error of >= MSE
  if (mse ~= 1)
      I = find(cumsum(eigen) >= mse,1);
  else
      I = nChannels;
  end;
  
% Create matrix of components to keep (principal)
  principal = vector(:,1:I);
  sigma     = eigen(1:I);

% Data is transformed to the new (possibly lower-dimensional) reference frame
  data = data*principal;
  
end % END principalcomponents