%本版本為加速版hurst計(jì)算程序
clear
% B=normrnd(0,1,1,1024)
% B=B+5
% y=B(1:1024)
load z1024h5.mat%待分析數(shù)據(jù)
y=z(1:1024)

t=length(y)
tt=[1:t]
for i=1:t %計(jì)算各個(gè)時(shí)間標(biāo)度的平均值
    PJ(i)=sum(y(1:i))/i
end
X=0
for i=2:t %計(jì)算累計(jì)離差
    X=X+(PJ(i-1)-PJ(i))*tt(1:i-1)
    X(i)=sum(y(1:i))-PJ(i)*i
    R(i)=max(X)-min(X)
end
for i=1:t
    S(i)=std(y(1:i))
% YY=y(1:i)-PJ(i)
% YYY=YY.^2
% S(i)=(sum(YYY)/i)^0.5    
    RS(i)=R(i)/S(i)
end
RS_log=log(RS(3:t))
%tt=tt/t
t_log=log(tt(3:t))
P=polyfit(t_log,RS_log,1)%最小二乘法擬合數(shù)據(jù)
xx=[floor(min(t_log)):ceil(max(t_log))]
yy=P(1)*xx+P(2)
plot(xx,yy)
hold on
plot(t_log,RS_log,'.')
grid on