%____________________________________________________________________________________________
%Written by Paolo B. DePetrillo (2003)
%Revised to include slope standard deviation 2007
%If this routine is used in your work, please reference the following:
%DePetrillo PB, Speers D'A, Ruttiman UE. Determining the hurst exponent
%of a fractal time series and its application to EKG signal analsyis.
%Comput Biol Med 29:6, 393-406. 1999.
%_____________________________________________________________________________________________


%input vector T is a column containing the magnitudes of the interbeat interval, such as in ms
%input scalar E is the embedding dimension for phase space
function [ H ] = hurst(T, E)

T = T';											%B assumes column structure for input vector
B=T;											%transposed to row vector

%matrix representing phase space
for i = 1:E-1									%Generate array from input vector T
T(:,[end])=[];									%Remove last point of T
B(:,[1])=[];									%remove first point of B
T=[T;B];										%concatenate T and B
end												%E is embedding dimension
												

Max=round(sqrt(size(T,2)));						%figure out maximum hop
												%which is the sqrt of the number of columns of final array T

%curve length																							
for i = 1:Max									%i is the hop size			
	D=diff(T(:,[1:i:size(T,2)]),1,2);			%take the successive differences based on hop size i
	S=D.^2;										%and square them
	L(i)=(sum(sqrt(sum(S))))/i;					%now L(i) is the length of the curve
end												%at progressive hop size i
												%divide by hop size to normalize
											

%plotting
Win=log([1:Max])'							%usual linear regression
Y=log(L)'										%matrix transpose so work with columns
X=[ones(size(Win)) Win];
a=X\Y;

SX=sum(Win);
SY=sum(Y);
SXX=sum(Win .* Win)
SYY=sum(Y.*Y);
SXY=sum(Win.*Y);
LX=length(Win);
LY=length(Y);
XM=SX/LX;
XY=SY/LY;
QX=SXX-SX*SX/LX;
QY=SYY-SY*SY/LY;
QXY=SXY-SX*SY/LX;
fractaldim=QXY/QX
slope=fractaldim;
QYdotX=QY-slope*QXY;
SYdotX=sqrt(QYdotX/(LX-2));
stdevfractaldimension=SYdotX/sqrt(QX)


%output
D = -a(2);										%D is fractal dimension
H = 2-D											%H is Hurst exponent

plot(Win,Y)
XLABEL('log(r)')
YLABEL('log((L(r)/r))')