function [F alpha pred] = detrended_fluctuation_analysis(X,nsamps,overlap)

% [F alpha pred] = detrended_fluctuation_analysis(X,nsamps,overlap) 
% 
% Inputs:
% X       - vector of power over time
% nsamps  - vector of sample numbers defining the time segments used for the
%           analysis; segment lengthes should be logarithmic
% overlap - expressed as proportion, e.g., 0 - no overlap, 0.5 - 50% overlap
%
% Outputs:
% F     - dfa fluctuation parameter
% alpha - slope of linear fit to log10(dfa) as a function log10(nsamps);
%         alpha ~0.5 --> uncorrelated signal (white noise)
%         alpha >0.5 --> temporal correlations over the range of tau
% pred  - predicted dfa values from linear fit
%
% Steps:
% (1) mean subtraction, (2) cummulative sum, (3) detrend segments, (4) rms
%  of detrended sgement, (5) fitting a linear function to dfa
%
% References:
% Linkenkaer-Hansen K et al. (2001) Long-Range Temporal Correlations and 
%    Scaling Behavior in Human Brain Oscillations. J Neurosci 21:1370–1377.
% Smit DJA et al. (2013) Long-Range Temporal Correlations in Resting-State Alpha
%    Oscillations Predict Human Timing-Error Dynamics. J Neurosci 33:11212–11220.
%
% Description:
% The script calculates long-range temporal correlations of oscillatory
% activity. It relates to the slope of the 1/f scale-free properties of
% power spectra. It is an autocorrelation analysis for non-stationary data.
% ---------
%
%    Copyright (C) 2015, B. Herrmann
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
% ----------------------------------------------------------------------
% B. Herrmann, Email: herrmann.b@gmail.com, 2015-10-23

% subtract mean from amplitude envelope
X = X(:);
X = X - mean(X);
nsamps = nsamps(:);

% calculate cummulative sum
CX = cumsum(X);

% get shifting vector
shift = round(nsamps*(1-overlap));

F = zeros([length(nsamps) 1]);
for ii = 1 : length(nsamps)
	nn = 1:nsamps(ii);
	counter = 0;
	while nn(end) < length(CX)
		% detrend segment
		tmpX = detrend(CX(nn));
		
		% calculate RMS
		F(ii)   = F(ii) + sqrt(sum(tmpX.^2)/length(tmpX));
		nn      = nn + shift(ii);
		counter = counter + 1;
	end
	F(ii) = F(ii)/counter;
end

coef  = polyfit(log10(nsamps),log10(F),1);   % linear fit
pred  = polyval(coef,log10(nsamps));         % predicted dfa values
alpha = coef(1);                             % slope
