浑沌分析 关于CC_method的改进算法

附：
主程序
% 此程序用来测试CC_method
% 2008-12-01
% zhangli
clear all
clear all
%利用方程获得
% 产生Lorenz 时间序列
% dx/dt = sigma*(y-x)
% dy/dt = r*x - y - x*z
% dz/dt = -b*z + x*y
sigma=16; % Lorenz 方程参数
b=4;
r=45.92;
y=[-1,0,1]; % 起始点(1 x 3 的行向量) 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，求tau
tic
[Smean_inf,Sdeltmean_inf,Scor_inf,tau_inf,tw_inf]=C_CMethod _inf(X,max_d);
toc
tau_inf
tw_inf
% 相关作图
figure('name','CC法求时间延迟');
plot(1:max_d,Smean_inf,'-b');hold on;
plot(1:max_d,Sdeltmean_inf,'-*c');hold on;
plot(1:max_d,Scor_inf,'-m');hold on;
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);
子函数1
function
[Smean,Sdeltmean,Scor,tau,tw]=C_CMethod_inf(X,max_d) % 用于求延迟时间tau
% X为输入时间序列
% max_d为最大时间延迟
% Smean，Sdeltmean,Scor为返回值
% tau为计算得到的延迟时间
% tw为时间窗口
% zhangli
% 2008-11-30
N=length(X);
Smean=zeros(1,max_d);
Scmean=zeros(1,max_d);
Scor=zeros(1,max_d);
delt=std(X);
% 计算Smean,Sdeltmean,Scor
for t=1:max_d
S=zeros(4,4);
Sdelt=zeros(1,4);
for m=2:5
for j=1:4
r=delt*j/2;
Xdt=disjoint(X,N,t); % 将时间序列X分解成t个不相交的时间序列
Xdt=Xdt';
s=0;
for tau=1:t
N_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个不相关的时间序列求和
end
S(m-1,j)=s/tau;
end
Sdelt(m-1)=max(S(m-1,:))-min(S(m-1,:)); % 差量计算
end
Smean(t)=mean(mean(S)); % 计算平均值
Sdeltmean(t)=mean(Sdelt); % 计算平均值
Scor(t)=abs(Smean(t))+Sdeltmean(t);
end
% 寻找时间延迟tau：即Sdeltmean第一个极小值点对应的t for i=2:length(Sdeltmean)-1
if
Sdeltmean(i)<Sdeltmean(i-1)&Sdeltmean(i)<Sdeltmean(i+1) tau=i;
break;
end
end
% 寻找时间窗口tw：即Scor最小值对应的t
for i=1:length(Scor)
if Scor(i)==min(Scor)
tw=i;
break;
end
end
子函数2
function data_d=disjoint(data,N,t)
% 此函数用于将时间序列分解成t个不相交的时间序列
% data:输入时间序列
% N:data的长度
% t:the index lag
% data_d：返回分解后的t个不相交的时间序列% 2008-11-28
% zhangli
for i=1:t
for j=1:(N/t)
data_d(i,j)=data(i+(j-1)*t);
end
end
子函数3 function Data=reconstitution(data,N,m,tau) % 该函数用来重构相空间
% data为输入时间序列
% N为时间序列长度
% m为嵌入空间维数
% tau为时间延迟
% Y为输出,是M*m维矩阵
% 2008-11-26
% zhangli
M=N-(m-1)*tau;
Data=zeros(m,M);
for i=1:m
Data(i,:)=data([((i-1)*tau+1):1:((i-1)*tau+M)]);
end
子函数4
function C=correlation_integral_inf(Y,M,r)
% 此函数用于计算关联积分,取无穷范数
% Y为重构的相空间
% M为相空间中点的个数
% r为搜索半径
% Y为输出,是M*m维矩阵
% 2008-11-30
% zhangli
C=0;
for i=1:M-1
for j=i+1:M
d1=norm((Y(i,:)-Y(j,:)),inf); % 计算状态空间中每两点之间的距离,取无穷范数
if r-d1>=0
C=C+1;
end
end
end
C=2*C/(M*(M-1));
6. 运行结果
图1 原始CC_method算法结果
7. 关于参考文献《非线性时间序列的相空间重构技术研究》(秦奕青,蔡卫东,杨炳儒,2008)
文章中提出了一种CC的改进算法。

由于不支持pdf的文档，所以我没有办法传上来，大家可以自己下载。

根据文章的思路，我对我的程序进行了改变，得到结果如图2：
却跟参考文献《非线性时间序列的相空间重构技术研究》中不一样。

请教各位高手。

谢谢！
图2 改进后的CC算法。