function [rra,tau,Hrs] = computerescaledrange(data)
% COMPUTERESCALEDRANGE Performs the Hurst rescaled range analysis.
%   COMPUTERESCALEDRANGE(DATA) returns the root-mean-square fluctuation for
%   a given box size.
%
%   [RRA, TAU] = COMPUTERESCALEDRANGE(DATA) additionally returns the box sizes.
%
%   [RRA, TAU, HRS] = COMPUTERESCALEDRANGE(DATA) returns rescaling range
%   coefficient HRS, which is also know as alpha.
%
%   Input
%       DATA | m x 1 | double | m is # samples
%
%   Output
%       RRA | floor(m/2)-9 x 1 | double | root-mean-square of fluctuation 
%       TAU | floor(m/2)-9 x 1 | double | interval (or box) length
%       HRS | 1 x 1            | double | Hurst rescaled range exponent
%   
%   Example 
%      data = randn(500,1); 
%      [rra,tau,hrs] = computerescaledrange(data); 
%
%   Reference
%     [1] Eke, A., "Fractal characterization of complexity in temporal 
%         physiological signal," Physiological Measurement, Vol. 23, R1-R38, 
%         2002 (Topical Review). 
%   
%     [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 COMPUTEREDETRENDED.  

% Created by: Kevin P. Granata Musculoskeletal Biomechanics Laboratory
%             http://www.biomechanics.esm.vt.edu
% Last Revision: 6 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;  

% Inteval lengths & pre-allocate for the rescaled range
  tau = 10:floor(m/2);        % Interval lengths
  rra = zeros(length(tau),1); % Pre-allocate
  
  for i = tau          
      nIntervals = floor(m/i); % Number of intervals of box size i         
      
    % For a given interval length, i, loop through the data computing the mean 
    % for each interval and cumulative departure from the mean, X  
      R = zeros(nIntervals,1); % Pre-allocate
      S = zeros(nIntervals,1); % Pre-allocate    
      
      for j = 1:nIntervals                  
          indicies  = (j-1)*i+(1:i);    % Desired indicies                  
          X = cumsum(data(indicies)-mean(data(indicies))); % Departure from mean                       
          R(j) = max(X)-min(X);       % Local "range"                  
          S(j) = std(data(indicies)); % Local stardard deviation          
      end  
    
    % Remove the boxes that have a standard deviation of zero 
    % Note:  This should only happen for small tau        
      R = R(S~=0);
      S = S(S~=0);    
      
      zerosRemoved = j-length(S);
      if zerosRemoved
          warning('rescaledrange:zerosrm',...
                  '%d zeros were removed when tau = %3d\n',zerosRemoved,i);
      end
      
    % Scale the ranges by the standard deviation
      rra(i-tau(1)+1) = mean(R./S);               
  end

% Fit the loglog series
% rra (is proportional to) tau^(Hrs) 
  tau   = tau(:);
  coeff = polyfit(log10(tau),log10(rra),1);
  Hrs = coeff(1); 
  
end % END computerescaledrange