关于G-P算法计算关联维数

6958192460

以下是我改写的gp算法计算关联维数的程序：
function  [ln_r,ln_C]=G_P(data,tau,m,ss)n=length(data);     % 时间序列长度
    N=n-(m-1)*tau;     % 相空间中点的个数   
for j=1:N      
for i=1:m         
Y(i,j)=data((i-1)*tau+j);     % 相空间重构      
end   
end      
for a=1:N-1      
for b=a+1:N         
D(a,b)=max(abs(Y(:,a)-Y(:,b)));     % 计算相空间中每两点之间的距离      
end   
end
max_D=max(max(D));     % 得到所有点间的最大距离   
D(1,1)=max_D;   min_D=min(min(D));     % 得到所有点间的最短距离   
delt=(max_D-min_D)/ss;     % 得到r的步长
k=1:1:ss      
r=min_D+k*delt;     % r的范围   
sum=0;  
for h=1:N     
for g=h+1:N         
d=norm((Y(:,h)-Y(:,g)),inf);     % 以范数计算重构空间两点间的距离
          site=heaviside(r-d);     % 计算Heaviside函数的值
          sum=sum+site;     
end   
end  
C=2*sum/(N^2);     % 关联函数的值  
CC(k)=log(C);  
rr(k)=log(r);   
CC=log(C);   
rr=log(r);  
plot(rr,CC);     % 画出lnr—lnC(r)关系图   
hold on;
% 以下是最小斜率法拟合曲线求斜率的程序：
function xielv(CC(k),rr(k))
y=polyfit(rr,CC,4);
y1=polyder(y);
y2=polyder(y1);
for k=1:ss   
K(k)=(abs(polyval(y2,rr(k))))/(1+(polyval(y1,rr(k)))^2)^(3/2);
end
min_K=min(K);
a=0;
for aa=1:ss   
if K(aa)==min_K        
a=aa;      
break   
end
end
a2=0;
for i=a:ss   
bb=rr(i)-rr(a);   
cc=abs(bb*polyval(y2,rr(i))/(1+(polyval(y1,rr(i)))^2)^(1/2));   
if cc>=(pi/90)      
a2=i;      
break   
end
end
a1=0;
for p=a:-1:1   
BB=rr(a)-rr(p);   
DD=abs(BB*polyval(y2,rr(i))/(1+(polyval(y1,rr(i)))^2)^(1/2));   
if DD>=(pi/90)        
a1=p;        
break   
end  
end
h=0;
for hh=a1:a2   
h=h+1;   
rrr(h)=rr(hh);   
CCC(h)=CC(hh);
end
KK=polyfit(rrr,CCC,1);
利用Lorenz系统生成的时间序列计算m=3时k=1.6782，而我看其他人算的都是2.0左右，请高手解答一下，程序有问题吗？特别是前半部分另外，这个程序运行大数据很慢，求高手改进。

[ 本帖最后由 6958192460 于 2010-7-18 20:00 编辑 ]