function [dfa,tau,alpha] = computedetrended(data)
% COMPUTEDETRENDED Perform detrended fluctuation analysis.
%   COMPUTEDETRENDED(DATA) returns the root-mean-square fluctuation for
%   a given box size.
%
%   [DFA, TAU] = COMPUTEREDETRENDED(DATA) additionally returns the box sizes.
%
%   [DFA, TAU,ALPHA] = COMPUTEREDETRENDED(DATA) returns detrended fluctuation 
%   coefficient, ALPHA.
%
%   Input
%       DATA | m x 1 | double | m is # samples
%
%   Output
%       DFA   | floor(m/2)-9 x 1 | double | root-mean-square of fluctuation 
%       TAU   | floor(m/2)-9 x 1 | double | interval (or box) length
%       ALPHA | 1 x 1            | double | Hurst rescaled range exponent
%   
%   Example 
%      data = randn(500,1); 
%      [dfa,tau,alpha] = computedetrended(data); 
%
%   Reference
%     [1] Peng, C-K, "Quantification of scaling exponents and crossover 
%         phenomena in nonstationary heartbeat time series," Chaos, Vol. 5, 
%         No. 1, 1995. 
%
%   See also COMPUTERESCALEDRANGE.  

%Created by: Jesse Norris, 
%            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;  
  
  avg = mean(data);       %Store the average of the entire data set

% Create the cumulative sum of the difference from the average            
  X = zeros(size(data)); 
  for ind = 1:length(data)
      X(ind) = sum(data(1:ind)-avg);
  end
  X = X - mean(X);  %Mean zero data (not sure if this is necessary)

% Compute the interval lengths & pre-allocate for dfa
  tau = 10:floor(m/2);         % Interval lengths
  dfa = zeros(length(tau),1);  % Pre-allocate

  for i = tau 
      
    % Inteval start points 
    %  start = 1:i:i*floor(m/i); 
      nIntervals = floor(m/i); % Number of intervals of box size i         
      
    % The squared difference between the linear fit and data points summed for
    % all boxes of length tau   
      sqdiff = 0;
    
    % For a given interval length, tau, loop through all the boxes, points at 
    % the end that do not fit a box are discarded    
      for j = 1:nIntervals
          
         % Fit the data in the box of length tau beginning at start(j)          
           xfit  = [(1:i)' ones(i,1)];    
           coeff = xfit\X((j-1)*i+(1:i));
            
          % Compute the sum of squared difference between the fit and data
            diff   = X((j-1)*i+(1:i))-xfit*coeff;
            sqdiff = sqdiff + diff'*diff;
      end
      
    % Calculate the root-mean-square fluctation by dividing by the number of 
    % data points included in the sum (Note: because points are discarded for 
    % some box lengths, the number of data points included is not always the 
    % same as the total number of data points, e.g. length of x)            
      dfa(i-tau(1)+1) = sqrt(1/(nIntervals*i)*sqdiff);      
  end

% Fit the loglog series
  tau   = tau(:);
  fit   = polyfit(log10(tau),log10(dfa),1);
  alpha = fit(1);    
  
end % END computedetrended