function [d,k,rSq] = ornsteinuhlenbeck(msd,dt)
% ORNSTEINUHLENBECK Fit mean square displacement to Ornstein-Uhlenbeck model.
%   [D,K] = ORNSTEINUHLENBECK(MSD,DT) where MSD is the mean square displacement 
%   data and DT is the sampling increment, returns the diffusion, D, and drift, 
%   K coefficients for the following fit:
%
%                           /                \
%                      D   |      (-2 K tau)  |
%         MSD(tau) = ----- | 1 - e            |
%                     2 K   \                /
%
%         where: 
%           tau is the time between samples
%
%   [D,K,RSQ] = ORNSTEINUHLENBECK(MSD,DT) returns RSQ value for the fit.
%
%   Input
%       MSD | m x 1 | double | mean square displacements for m time lags
%       DT  | 1 x 1 | double | time-increment
%
%   Output
%       D   | 1 x 1 | double | diffusion coefficient
%       K   | 1 x 1 | double | drift coefficient
%       RSQ | 1 x 1 | double | R^2: coefficient of determination
%   
%   Example 
%      dt  = 0.01;                                 % Time increment in seconds
%      tau = (0:dt:10)';                           % Lag times in seconds
%      b   = 0.015*[1 1]; a   = [1 -0.97];         % Filter coefficients
%      noise = tau.*filtfilt(b,a,randn(size(tau)));% Appropriate noise
%      msd = 10/(2*0.5)*(1-exp(-2*0.5*tau))+noise; % Simulated mean-square-disp
%      [d,k]  = ornsteinuhlenbeck(msd,dt);         % Perform fit
%      msdfit = d/(2*k)*(1-exp(-2*k*tau));         % Fit
%      plot(tau,msd,tau,msdfit);                   % Plot fit
%
%   Reference
%     [1] Newell, K.M., et al. "Stochastic processes in postural center-of-
%         pressure profiles," Experimental Brain Research, Vol. 113, 
%         pp. 158-164, 1997.
%  
%     [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: 15 July 2007

% Make sure input MSD is a column vector  
  [m n] = size(msd);
  if (n ~= 1)
      error('Input should be a column vector (i.e. dimension m x 1)');     
  end;  

% Ornstein-Uhlenbeck fit
% x(1) is the diffusion coefficient, d
% x(2) is the drift coefficient, k
  fun = @(x,tau) x(1)/(2*x(2))*(1-exp(-2*x(2)*tau));
  xo  = [10 0.5];
  
  
% Compute the Ornstein-Uhlenbeck fit 
  tau = dt*(0:length(msd)-1)'; 
  x   = lsqcurvefit(fun,xo,tau,msd);
  d   = x(1);
  k   = x(2); 
  
% Compute Coefficient of Determination, R^2
  sse = msd-fun(x,tau); % Sum of squared errors
  sst = msd-mean(msd); % Total sum of squares
  rSq = 1 - sse'*sse/(sst'*sst);
  
end % END ornsteinuhlenbeck