function [FV,dtV] = DFAmFB(xV,m)
% [FV,dtV] = DFAmFB(xV,m)
% DFAmFB computes the Hurst exponent for a given time series 'xV' using
% the method of detrended fluctuation analysis (DFA). The segments are
% powers of 2 and if the time series length is not exactly a power of 2
% then segments are formed both in forward and backward direction and the
% statistics are computed on all these segments.
% INPUTS:
% - xV      : The given scalar time series (vector of size n x 1).
% - m       : Order of the polynomial to fit the segments and remove the
%             trend. If not specified, the default m=1 is used. 
% OUTPUT
% - FV      : The vector of the square root of the average of mean square
%             fluctuations.
% - dtV     : the array of increasing time segment lengths.

if nargin==1
    m=1;
end
n = length(xV);
dtV = 2.^(2:log2(n))';
ndt = length(dtV);
yV=cumsum(xV-mean(xV));
FV = NaN*ones(ndt,1);
if dtV(ndt)==n
    FV = NaN*ones(ndt,1);
else    
    FV = NaN*ones(2*ndt,1);
    zV = flipud(yV);
end
if m==0
    for i=1:ndt
        w = dtV(i);
        nseg=floor(n/w);  % number of segments of the given length
        FdtV = (yV(1:w:(nseg-1)*w+1)-yV(w:w:nseg*w)).^2; % square of segment edge differences
        FV(i)=sqrt(mean(FdtV)); % square root of the mean of segment fluctuations
    end
    if dtV(ndt)<n        
        for i=1:ndt
            w = dtV(i);
            nseg=floor(n/w);  % number of segments of the given length
            FdtV = (zV(1:w:(nseg-1)*w+1)-zV(w:w:nseg*w)).^2; % square of segment edge differences
            FV(ndt+i)=sqrt(mean(FdtV)); % square root of the mean of segment fluctuations
        end
    end    
else
    for i=1:ndt
        w = dtV(i);
        nseg=floor(n/w);  % number of segments of the given length
        fityV=NaN*ones(nseg*w,1); 
        for j=1:nseg
            bV=polyfit([1:w]',yV((j-1)*w+1:j*w),m); % parameters of polynomial fit of degree m at each segment
            fityV((j-1)*w+1:j*w)=polyval(bV,[1:w]'); % the fitted values at each segment
        end
        resV=yV(1:nseg*w)-fityV; % residual of fit over all consecutive segments
        FV(i)=sqrt(resV'*resV/(nseg*w)); % square root of the mean of segment fluctuations
    end
    if dtV(ndt)<n
        for i=1:ndt
            w = dtV(i);
            nseg=floor(n/w);  % number of segments of the given length
            fitzV=NaN*ones(nseg*w,1); 
            for j=1:nseg
                bV=polyfit([1:w]',zV((j-1)*w+1:j*w),m); % parameters of polynomial fit of degree m at each segment
                fitzV((j-1)*w+1:j*w)=polyval(bV,[1:w]'); % the fitted values at each segment
            end
            resV=zV(1:nseg*w)-fitzV; % residual of fit over all consecutive segments
            FV(ndt+i)=sqrt(resV'*resV/(nseg*w)); % square root of the mean of segment fluctuations
        end
    end
end
if dtV(ndt)<n
    FM = reshape(FV,ndt,2);
    FV = sqrt(mean(FM.^2,2));
    % FV = (mean(FM'))'; changed after comment of Christos Matsoukas matsoukas@aegean.gr
end
    