function FV = DFAm(xV,dtV,m)
% FV = DFAm(xV,dtV,m)
% DFAm computes the Hurst exponent for a given time series 'xV' using
% the method of detrended fluctuation analysis (DFA).
% INPUTS:
% - xV      : The given scalar time series (vector of size n x 1).
% - dtV     : the array of increasing time segment lengths.
% - 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.

if nargin==2
    m=1;
end
n = length(xV);
ndt = length(dtV);
FV = NaN*ones(ndt,1);
yV=cumsum(xV-mean(xV));
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
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
end
    