function [msd,lags] = meansquaredisplacement(data,lagmax)
% MEANSQUAREDISPLACEMENT Calculate the mean square displacement.
%   MEANSQUAREDISPLACEMENT(DATA,LAGMAX) returns the mean square displacement as 
%   a function of time increments from 0 to LAGMAX-1.
%
%                           /                            \
%                       1   | m-j                        |
%         MSD(j+1)  = ----- | sum (data(i+j) - data(i))^2|
%                      m-j  |i = 1                       |
%                           \                            /
%
%         where: 
%           j is [0, lagmax-1] and corresponds to delta time increments of dt*j
%           m is the number of sample points
%
%   [MSD,LAGS] = MEANSQUAREDISPLACEMENT(DATA,LAGMAX) returns lags which are from
%   0 to LAGMAX-1.
%
%   Input
%       DATA | m x 1 | double | m is # samples
%
%   Output
%       MSD   | LAGMAX x 1 | double | mean square displacements
%       LAGS  | LAGMAX x 1 | double | lag indicies
%   
%   Example 
%      data = randn(500,1); 
%      [msd,lags] = meansquaredisplacement(data,100); 
%      plot(lags,msd);
%
%   Reference
%     [1] Collins, JJ, and De Luca, CJ, "Open-loop and closed-loop control of
%         posture: a random-walk analysis of center-of-pressure trajectories,"
%         Experimental Brain Research, Vol. 95, pp. 308-318, 1993.
%  
%     [2] Norris, J., et al. "Ability of static and statistical mechanics
%         posturographic measures to distinguish between age and fall-risk," 
%         Journal of Biomechanics, Vol. 38, No. 6, pp. 1263-1272, 2005.
%
%   See also .  

%Created by: Jesse Norris, 
%            Kevin P. Granata Musculoskeletal Biomechanics Laboratory
%            http://www.biomechanics.esm.vt.edu
% Last Revision: 9 July 2007

% Make sure input DATA is a column vector  
  [m n] = size(data);
  if (n ~= 1)
      error('Input should be a column vector (i.e. dimension m x 1)');     
  end;  
  
% Handle error where LAGMAX is greater than number of samples
  if (lagmax >= m)
      error('Input LAGMAX needs to be less than the number of samples');
  end;
      
% Pre-allocate output
  lags = 0:lagmax-1;
  msd  = zeros(lagmax,1);  
  
  for j = lags
      disp     = data(j+1:m)-data(1:m-j);
      msd(j+1) = disp'*disp/(m-j);
  end
  
% Return as column vector
  lags = lags(:);
  
end % END meansquaredisplacement