求救：哪位大侠有改进的CC算法，感激不尽

zhl20062159

最近用吕金虎那本书里介绍的C-C算法求嵌入维和时间延迟，然后代入G-P算法中求关联维，并用Chaos Toolbox Ver.2.0工具箱里的C-C、G-P程序仿真运行了下，结果始终不尽如意。我用logistics混沌序列进行验证，具体过程如下：首先当u=4时，我产生40000个混沌序列点，然后代入C-C算法，画出了三条曲线图如下：
file:///C:/DOCUME~1/ME/LOCALS~1/Temp/msohtml1/01/clip_image002.gif
蓝线----smean(t)

绿点---delt_s
红*-----s_cor     由图求的delt_s第一个极小值点：4
s_cor最小值点在9

即：Tw=9,tau=4.
然后我代入G_P算法算出m=2:10下的ln_C~ln_r，如图：
file:///C:/DOCUME~1/ME/LOCALS~1/Temp/msohtml1/01/clip_image002.gif

但是画出的不同嵌入维下ln_C~ln_r的斜率始终不趋于稳定。
并且下面几条曲线用最小二乘法拟合的时候，输出均为Inf，不知道为什么。


起初以为数据格式不够，我也试着改了下，还是不好使。后来看到一篇论文里用改进的C-C算法，就是在关联积分计算过程中引入了权衡计算精度与速度的可调参数。效果还不错，不知哪位大侠有没有相关程序，发给小妹一份，不胜感激。



还有在论坛上看到用fft求最佳延迟时间，如果哪位大侠有程序也发给小妹看一下。谢谢，谢谢啊~~~


我的邮zhl20062159@126.com


logistic 序列产生程序：
clear all;
N=40000;
% initial conditions:
x(1)=0.2;
% parameter:
% r=3.8;  % chaos
% r=2.0;  % steady state
% r=3.2;   % period 2
r=4;   % period 4?
for n=1:N
  x(n+1)=r*x(n)*(1-x(n));
end
plot(x,'k')
title('logistic map')
xlabel('n')
ylabel('x(n)')
figure
plot(x(2:N),x(1:N-1),'k.','markersize',5)
title('logistic map')
xlabel('x(n+1)')
ylabel('x(n)')

fid=fopen('logistic.txt','w');
fprintf(fid,'%f\n',x);
fclose(fid);


C-C算法：
function [s,delt_s,s_cor]=C_CMethod(data)
%this function calculate time delay and embedding demension with C-C
% 该函数用C-C算法来计算延迟时间和嵌入维数
%Method,which proved by H.S.Kim
%skyhawk&flyinghawk
% %****************调试程序段****************************
clear all;
data1=load('E:\chaos\4.15\Chaos Toolbox Ver.2.0_zhu\Chaos Toolbox Ver.2.0\Main\logistic.txt');
data=data1(1:20:length(data1));         %%数据采样
% max(data)
% min(data)
% for i=1:length(data)                    %%数据的归一化
%     data(i)=(data(i)-min(data))/(max(data)-min(data));
% end
% %************************************************

N=length(data)                                      % 数据向量长度
max_d=20;    %the maximum value of the time delay     时间延迟的最大值

sigma=std(data) %calcute standard deviation s_d       计算标准偏差s_d

  for t=1:max_d
    t
    s_t=0;
    delt_s_s=0;
    for m=2:5
        m
        s_t1=0;
        for j=1:4  %zhu:对r的循环
            j
            r=sigma*j/2;
            data_d=disjoint(data,N,t);    %将时间序列分解成t个不相交的时间序列
            [ll,N_d]=size(data_d)
            s_t3=0;
            for i=1:t                     %zhu:对t个序列分别计算关联积分C
                i
                Y=data_d(i,:);            %zhu:取data_d中的第i行，即第i个子序列
                C_1(i)=correlation_integral(Y,N_d,r);            %计算C(1,N_d,r,t)  zhu:即长度为N_d的一维序列
                X=reconstitution(Y,N_d,m,t);                     %相空间重构，zhu:新的序列为总长度为N_d产生的m维序列
                N_r=N_d-(m-1)*t;                                 %zhu:t=tao      新的序列个数为N_r
                C_I(i)=correlation_integral(X,N_r,r);            %计算C(m,N_r,r,t)
                s_t3=s_t3+(C_I(i)-C_1(i)^m);                     %对t个不相关的时间序列求和
            end
            s_t2(j)=s_t3/t;                                      %计算S(m,N,r,t)
            s_t1=s_t1+s_t2(j);                                   %对rj求和
         end
        delt_s_m(m)=max(s_t2)-min(s_t2);                         %求delt S(m,t)
        delt_s_s=delt_s_s+delt_s_m(m);                           %delt S(m,t)对m求和
        s_t0(m)=s_t1;
        s_t=s_t+s_t0(m);                                         %S对m求和
    end
    s(t)=s_t/16;
    delt_s(t)=delt_s_s/4;
    s_cor(t)=delt_s(t)+abs(s(t));
  end

t=1:max_d;
plot(t,s,t,delt_s,'.',t,s_cor,'*')

fid=fopen('s(t).txt','w');
fprintf(fid,'%f\n',s(t));
fclose(fid);
fid=fopen('delt_s(t).txt','w');
fprintf(fid,'%f\n',delt_s(t));
fclose(fid);
fid=fopen('s_cor(t).txt','w');
fprintf(fid,'%f\n',s_cor(t));
fclose(fid);

% 寻找时间延迟tau：即Sdeltmean第一个极小值点对应的t
for i=2:length(delt_s)-1
    if delt_s(i)<delt_s(i-1)&delt_s(i)<delt_s(i+1)
        tau=i;
        break;
    end
end
% 寻找时间延迟tau：即Sdeltmean第一个极小值点对应的t



% 寻找时间窗口tw：即Scor最小值对应的t
for i=1:length(s_cor)
    if s_cor(i)==min(s_cor)
        tw=i;
        break;
    end
end
% 寻找时间窗口tw：即Scor最小值对应的t


G-P算法：
function [ln_r,ln_C]=G_P(data,N,tau,min_m,max_m,ss)
% the function is used to calculate correlation dimention with G-P algorithm
%    计算关联维数的G－P算法
% data:the time series                       时间序列
% N: the length of the time series           时间序列长度
% tau: the time delay                        时间延迟
% min_m:the least embedded dimention m       最小的嵌入维数
% max_m:the largest embedded dimention m     最大的嵌入维数
% ss:the stepsize of r                       r的步长        !!!我的理解—ss—r的循环数目!!!
%skyhawk
%%%%%%%%%%%zhu    test %%%%%%%%%%%%%%%%%%%%%
clear all;
data1=load('E:\chaos\4.15\Chaos Toolbox Ver.2.0_zhu\Chaos Toolbox Ver.2.0\Main\logistic.txt');
data=data1(1:20:length(data1));
% max(data)
% min(data)
% for i=1:length(data)                    %%数据的归一化
%     data(i)=(data(i)-min(data))/(max(data)-min(data));
% end

min_m=2;
max_m=10;
N=length(data)
tau=4;
ss=20;
%%%%%%%%%zhu    test %%%%%%%%%%%%%%%%%%%%%

for m=min_m:max_m
    m
    Y=reconstitution(data,N,m,tau);        %reconstitute state space
    M=N-(m-1)*tau;                         %the number of points in state space
    for i=1:M-1
%         i
        for j=i+1:M
            d(i,j)=max(abs(Y(:,i)-Y(:,j)));%calculate the distance of each two           
        end                 %points in state space  计算状态空间中每两点之间的距离
    end
    max_d=max(max(d))    %the max distance of all points   得到所有点之间的最大距离
    d(1,1)=max_d;
    min_d=min(min(d))     %the min distance of all points   得到所有点间的最短距离
    delt=(1.2*max_d-min_d)/ss;     %the stepsize of r            得到r的步长
    for k=1:ss      
        k
        r=min_d+k*delt;
        C(k)=correlation_integral(Y,M,r);               %calculate the correlation integral
        ln_C(m,k)=log(C(k));                            %lnC(r)
        ln_r(m,k)=log(r);                               %lnr
        fprintf('%d/%d/%d/%d\n',k,ss,m,max_m);
        fprintf('%d/%d/%d\n',k,ss,m);
    end
    plot(ln_r(m,:),ln_C(m,:));
    xlabel('ln_r');
    ylabel('ln_C(r)');   
    hold on;
end   

%%%%%%%%%%%%%%%%%%% 画 Correlation dimensional versus embedding dimension.
figure
S=[];ii=1;
for mm=min_m:max_m
    x=1:floor(length(ln_C(mm,:)));
    pp=polyfit(ln_r(mm,x),ln_C(mm,x),1) %%曲线拟合，最小二乘多项式拟合，n=1，为直线拟合。
    S(ii)=pp(1)  %%% 返回斜率。
    ii=ii+1;
end
MM=min_m:max_m
plot(MM,S,'*');
xlabel('嵌入维m');
ylabel('关联维D');
%%%%%%%%%%%%%%%%%%% 画 Correlation dimensional versus embedding dimension.


fid=fopen('lnr.txt','w');
fprintf(fid,'%f %f\n',ln_r);
fclose(fid);
fid = fopen('lnC.txt','w');
fprintf(fid,'%f %f\n',ln_C);
fclose(fid);