function [RSV,dtV] = RSanalysisFB(xV)
% [RSV,dtV] = RSanalysisFB(xV)
% RSanalysisFB computes the Hurst exponent for a given time series 'xV'
% using the method of rescaled range analysis (R/S analysis). 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).
% OUTPUTS
% - RSV     : The vector of average fluctuations (ratio of range over SD).
% - dtV     : the array of increasing time segment lengths.

n = length(xV);
dtV = 2.^(2:log2(n))';
ndt = length(dtV);
if dtV(ndt)==n
    RSV = NaN*ones(ndt,1);
else
    RSV = NaN*ones(2*ndt,1);
    zV = flipud(xV);
end
for i=1:ndt
    w = dtV(i);
    nseg=floor(n/w);  % number of segments of the given length
    xM=reshape(xV(1:w*nseg),w,nseg); % a segment at each column
    isegV=find(std(xM,1)~=0); % check for segments with all values equal (no variance)
    if ~isempty(isegV)
        xM=xM(:,isegV); % make calculations only on those segments with non-zero variance
        mxV = mean(xM);
        yM=cumsum(xM-kron(mxV,ones(w,1)),1);
        RSV(i)=mean(range(yM)./std(xM,1));
    end
end
if dtV(ndt)<n
    for i=1:ndt
        w = dtV(i);
        nseg=floor(n/w);  % number of segments of the given length
        zM=reshape(zV(1:w*nseg),w,nseg); % a segment at each column
        isegV=find(std(zM,1)~=0); % check for segments with all values equal (no variance)
        if ~isempty(isegV)
            zM=zM(:,isegV); % make calculations only on those segments with non-zero variance
            mzV = mean(zM);
            yM=cumsum(zM-kron(mzV,ones(w,1)),1);
            RSV(ndt+i)=mean(range(yM)./std(zM,1));
        end
    end
    RSM = reshape(RSV,ndt,2);
    RSV = (mean(RSM'))';
end

