关于CC_method的改进算法

附：1.主程序%此程序用来测试CC_method%2008-12-01%zhangliclearallclearall%利用方程获得%产生Lorenz时间序列%dx/dt=sigma*(y-x)%dy/dt=r*x-y-x*z%dz/dt=-b*z+x*ysigma=16;%Lorenz方程参数b=4;r=45.92;y=[-1,0,1];%起始点(1x3的行向量)h=0.01;%积分时间步长k1=10000;%前面的迭代点数k2=3000;%后面的迭代点数Z=LorenzData(y,h,k1+k2,sigma,r,b);X=Z(k1+1:end,1);max_d=200;%最大延迟时间%调用C_CMethod_inf，求tautic[Smean_inf,Sdeltmean_inf,Scor_inf,tau_inf,tw_inf]=C_CMethod_inf(X,max_d);toctau_inftw_inf%相关作图figure('name','CC法求时间延迟');plot(1:max_d,Smean_inf,'-b');holdon;plot(1:max_d,Sdeltmean_inf,'-*c');holdon;plot(1:max_d,Scor_inf,'-m');holdon;plot(1:max_d,zeros(1,max_d),'r');title('C_CMethod_inf');xlabel('Lag');legend('S(t)平均值','ΔS(t)平均值','Scor_inf');%将数据保持下来fid=fopen('Smean_inf.txt','w');fprintf(fid,'%f\n',Smean_inf);fclose(fid);fid=fopen('Sdeltmean_inf.txt','w');fprintf(fid,'%f\n',Sdeltmean_inf);fclose(fid);fid=fopen('Scor_inf.txt','w');fprintf(fid,'%f\n',Scor_inf);fclose(fid);2.子函数1function[Smean,Sdeltmean,Scor,tau,tw]=C_CMethod_inf(X,max_d)%用于求延迟时间tau%X为输入时间序列%max_d为最大时间延迟%Smean，Sdeltmean,Scor为返回值%tau为计算得到的延迟时间%tw为时间窗口%zhangli%2008-11-30N=length(X);Smean=zeros(1,max_d);Scmean=zeros(1,max_d);Scor=zeros(1,max_d);delt=std(X);%计算Smean,Sdeltmean,Scorfort=1:max_dS=zeros(4,4);Sdelt=zeros(1,4);form=2:5forj=1:4r=delt*j/2;Xdt=disjoint(X,N,t);%将时间序列X分解成t个不相交的时间序列Xdt=Xdt';s=0;fortau=1:tN_t=floor(N/t);%分成的子序列长度Y=Xdt(:,tau);%每个子序列Cs1(tau)=correlation_integral_inf(Y,N_t,r);%计算C(1,N/t,r,t)Z=reconstitution(Y,N_t,m,1);%相空间重构Z=Z';M=N_t-(m-1);Cs(tau)=correlation_integral_inf(Z,M,r);%计算C(m,N/t,r,t)s=s+(Cs(tau)-Cs1(tau)^m);%对t个不相关的时间序列求和endS(m-1,j)=s/tau;endSdelt(m-1)=max(S(m-1,:))-min(S(m-1,:));%差量计算endSmean(t)=mean(mean(S));%计算平均值Sdeltmean(t)=mean(Sdelt);%计算平均值Scor(t)=abs(Smean(t))+Sdeltmean(t);end%寻找时间延迟tau：即Sdeltmean第一个极小值点对应的tfori=2:length(Sdeltmean)-1ifSdeltmean(i)Sdeltmean(i-1)&Sdeltmean(i)Sdeltmean(i+1)tau=i;break;endend%寻找时间窗口tw：即Scor最小值对应的tfori=1:length(Scor)ifScor(i)==min(Scor)tw=i;break;endend3.子函数2functiondata_d=disjoint(data,N,t)%此函数用于将时间序列分解成t个不相交的时间序列%data:输入时间序列%N:data的长度%t:theindexlag%data_d：返回分解后的t个不相交的时间序列%2008-11-28%zhanglifori=1:tforj=1:(N/t)data_d(i,j)=data(i+(j-1)*t);endend4.子函数3functionData=reconstitution(data,N,m,tau)%该函数用来重构相空间%data为输入时间序列%N为时间序列长度%m为嵌入空间维数%tau为时间延迟%Y为输出,是M*m维矩阵%2008-11-26%zhangliM=N-(m-1)*tau;Data=zeros(m,M);fori=1:mData(i,:)=data([((i-1)*tau+1):1:((i-1)*tau+M)]);end5.子函数4functionC=correlation_integral_inf(Y,M,r)%此函数用于计算关联积分,取无穷范数%Y为重构的相空间%M为相空间中点的个数%r为搜索半径%Y为输出,是M*m维矩阵%2008-11-30%zhangliC=0;fori=1:M-1forj=i+1:Md1=norm((Y(i,:)-Y(j,:)),inf);%计算状态空间中每两点之间的距离,取无穷范数ifr-d1=0C=C+1;endendendC=2*C/(M*(M-1));6.运行结果图1原始CC_method算法结果7.关于参考文献《非线性时间序列的相空间重构技术研究》(秦奕青,蔡卫东,杨炳儒,2008)文章中提出了一种CC的改进算法。由于不支持pdf的文档，所以我没有办法传上来，大家可以自己下载。根据文章的思路，我对我的程序进行了改变，得到结果如图2：却跟参考文献《非线性时间序列的相空间重构技术研究》中不一样。请教各位高手。谢谢！图2改进后的CC算法