声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 8157|回复: 25

[分形与混沌] 用CC_method计算时延和时延窗

[复制链接]
发表于 2009-3-4 17:48 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
今天使用CC_method方法计算了一下时延和时延窗

方程和参数如下:
function dy = Lorenz(t,y)
a=16;
b=4.0;
c=45.92;
dy=zeros(3,1);
dy(1)=-a*(y(1)-y(2));
dy(2)=-y(1)*y(3)+c*y(1)-y(2);
dy(3)=y(1)*y(2)-b*y(3);

计算使用的程序如下:
clear all
clear all
t0=0;
tf=130;
[t,x]=ode45(@Lorenz,[t0:0.01:tf],[-1,0,1]);
Lorenz_data=x(10002:end,1);
max_d=200;
% 调用C_CMethod_inf,求tau
tic
[Smean_inf,Sdeltmean_inf,Scor_inf,tau_inf,tw_inf]=C_CMethod(Lorenz_data,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('CC法求Lorenz时间延迟');xlabel('Lag');
legend('S(t)平均值','ΔS(t)平均值','Scor');

调用的子函数如下:
function [Smean,Sdeltmean,Scor,tau,tw]=C_CMethod(data,max_d)
% 本函数用于求延迟时间tau和时间窗口tw
% data:输入时间序列
% max_d:最大时间延迟
% Smean,Sdeltmean,Scor为返回值
% tau:计算得到的延迟时间
% tw:时间窗口
N=length(data);     %时间序列的长度
Smean=zeros(1,max_d);    %初始化矩阵
Scmean=zeros(1,max_d);
Scor=zeros(1,max_d);
sigma=std(data);      %计算序列的标准差
% 计算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=sigma*j/2;
            Xdt=disjoint(data,N,t);          % 将时间序列data分解成t个不相交的时间序列
            s=0;
           for tau=1:t
                N_t=floor(N/t);                          % 分成的子序列长度
                Y=Xdt(:,tau);                            % 每个子序列           
               
                %计算C(1,N/t,r,t),相当于调用Cs1(tau)=correlation_integral1(Y,N_t,r)            
                Cs1(tau)=0;
                for ii=1:N_t-1
                    for jj=ii+1:N_t
                        d1=abs(Y(ii)-Y(jj));     % 计算状态空间中每两点之间的距离,取无穷范数
                        if r>d1
                            Cs1(tau)=Cs1(tau)+1;            
                        end
                    end
                end
                Cs1(tau)=2*Cs1(tau)/(N_t*(N_t-1));
              
                Z=reconstitution(Y,N_t,m,1);             % 相空间重构
                M=N_t-(m-1);
                Cs(tau)=correlation_integral(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

function Data=reconstitution(data,N,m,tau)
%该函数用来重构相空间
% m为嵌入空间维数
% tau为时间延迟
% data为输入时间序列
% N为时间序列长度
% Data为输出,是m*n维矩阵
M=N-(m-1)*tau;%相空间中点的个数
Data=zeros(m,M);
for j=1:M
  for i=1:m           %相空间重构
    %Data(i,:)=data(((i-1)*tau+1):1:((i-1)*tau+M));
    Data(i,j)=data((i-1)*tau+j);
  end
end

%该函数用来计算计算关联积分
%C_I:关联积分的返回值
%X:重构的相空间矢量,是一个m*M的矩阵
%M::M是重构的m维相空间中的总点数
%r:Heaviside 函数中的搜索半径
sum_H=0;
for i=1:M-1
    for j=i+1:M
        d=norm((X(:,i)-X(:,j)),inf);%计算相空间中每两点之间的距离,其中NORM(V,inf) = max(abs(V)).
        if r>d   
        %sita=heaviside(r,d);%计算Heaviside 函数之值n
           sum_H=sum_H+1;
        end
    end
end
C_I=2*sum_H/(M*(M-1));%关联积分的值

运行结果所得与Kim论文(tau=10,tw=100)中的一致,但当把方程参数改为a=10,b=28,c=8/3是结果时延为5,时延窗为21,与文中结果(tau=18,tw=123)相差很大,请教各位高手

结果图.doc

34 KB, 下载次数: 98

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2009-3-5 10:02 | 显示全部楼层

回复 楼主 yuling 的帖子

参数变化会使结果变化很大的!
 楼主| 发表于 2009-3-5 15:10 | 显示全部楼层

回复 沙发 xiaocheng_2007 的帖子

首先感谢你的回复。
       我先是取参数a=16;b=4.0;c=45.92,其结果与论文中的一致,同时我也用文中Rossler方程的所给参数进行了验证,结果也和文中结果一致,这说明程序的编制是没有问题的。但是当用文中所给的另一组参数a=10,b=28,c=8/3计算时,结果大相径庭。鉴于Kim论文的权威性(他是在这篇论文中提出的CC_method),我才有此疑惑啊!
发表于 2009-3-7 16:46 | 显示全部楼层

回复 板凳 yuling 的帖子

cc上面程序是不是没贴全啊 ? 比如heaviside(r,d)函数 和计算关联积分函数?
还有你的cc程序运行需要多才时间呀?
 楼主| 发表于 2009-3-8 11:39 | 显示全部楼层

回复 地板 sandman 的帖子

谢谢你的提醒,我重新整理了一下程序和问题。(运行时间和结果在最后)
function dy = Lorenz(t,y)  %方程的定义
a=16;
b=4.0;
c=45.92;
dy=zeros(3,1);
dy(1)=-a*(y(1)-y(2));
dy(2)=-y(1)*y(3)+c*y(1)-y(2);
dy(3)=y(1)*y(2)-b*y(3);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%主程序
clear all
clear all
t0=0;
tf=130;
[t,x]=ode45(@Lorenz,[t0:0.01:tf],[-1,0,1]);
Lorenz_data=x(10002:end,1);
max_d=200;

% 调用C_CMethod_inf,求tau和tw
tic
[Smean_inf,Sdeltmean_inf,Scor_inf,tau_inf,tw_inf]=C_CMethod(Lorenz_data,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('CC法求Lorenz时间延迟');xlabel('Lag');
legend('S(t)平均值','ΔS(t)平均值','Scor');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
function [Smean,Sdeltmean,Scor,tau,tw]=C_CMethod(data,max_d)
% 本函数用于求延迟时间tau和时间窗口tw
% data:输入时间序列
% max_d:最大时间延迟
% Smean,Sdeltmean,Scor为返回值
% tau:计算得到的延迟时间
% tw:时间窗口
%yuling
N=length(data);     %时间序列的长度
Smean=zeros(1,max_d);    %初始化矩阵
Scmean=zeros(1,max_d);
Scor=zeros(1,max_d);
sigma=std(data);      %计算序列的标准差

% 计算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=sigma*j/2;
            Xdt=disjoint(data,t);        % 将时间序列data分解成t个不相交的时间序列
            s=0;
           for tau=1:t
                N_t=floor(N/t);                          % 分成的子序列长度
                Y=Xdt(:,tau);                            % 每个子序列           
               
                %计算C(1,N/t,r,t)         
                Cs1(tau)=0;
                for ii=1:N_t-1
                    for jj=ii+1:N_t
                        d1=abs(Y(ii)-Y(jj));     % 计算状态空间中每两点之间的距离,取无穷范数
                        if r>d1
                            Cs1(tau)=Cs1(tau)+1;            
                        end
                    end
                end
                Cs1(tau)=2*Cs1(tau)/(N_t*(N_t-1));
              
                Z=reconstitution(Y,m,1);             % 相空间重构
                M=N_t-(m-1);
                Cs(tau)=correlation_integral(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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Data=disjoint(data,t)
% 此函数用于将时间序列分解成t个不相交的时间序列
% data:输入时间序列
% t:延迟,也是不相交时间序列的个数
% Data:返回分解后的t个不相交的时间序列
%yuling
N=length(data);   %data的长度
for i=1:t
    for j=1:(N/t)
        Data(j,i)=data(i+(j-1)*t);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function C_I=correlation_integral(X,M,r)
%该函数用来计算计算关联积分
%C_I:关联积分的返回值
%X:重构的相空间矢量,是一个m*M的矩阵
%M::M是重构的m维相空间中的总点数
%r:Heaviside 函数中的搜索半径
%yuling
sum_H=0;
for i=1:M-1
    for j=i+1:M
        d=norm((X(:,i)-X(:,j)),inf);%计算相空间中每两点之间的距离,其中NORM(V,inf) = max(abs(V)).
        if r>d   
           sum_H=sum_H+1;
        end
    end
end
C_I=2*sum_H/(M*(M-1));%关联积分的值

运行结果如下,(见图Lorenz1.JPG)
Elapsed time is 662.015000 seconds.
tau_inf = 10   ,tw_inf =191,  这和Kim论文中的结果基本一致

当取参数a=10 ,b=28 ,c=8/3时。运行结果如下(见图Lorenz2.JPG)
Elapsed time is 659.718000 seconds.
tau_inf =  5   tw_inf =21,而这和Kim论文中的结果相差很大

同时对Rossler方程进行了验证,
function ff=Rossler(t,y)
d=0.2;e=0.2;f=5.0;
ff=[-y(2)-y(3);y(1)+d*y(2);e+y(3)*(y(1)-f)];
结果tau=16,tw=117,这和Kim论文中的结果完全一致,(见图Rossler.JPG)

Kim论文的结果见Kim_lunwen.JPG

[ 本帖最后由 yuling 于 2009-3-8 11:45 编辑 ]
Lorenz1.jpg
Lorenz2.jpg
Rossler.jpg
Kim_lunwen.JPG
发表于 2009-3-9 09:36 | 显示全部楼层

回复 5楼 yuling 的帖子

我看了一下你的程序。我认为你的function C_I=correlation_integral(X,M,r)
程序内容可能不正确!如果d=0的话,你的sum_H也相应地加了1。你可以使用heaviside函数
发表于 2009-3-9 09:45 | 显示全部楼层
另外function Data=disjoint(data,t)
你这个函数名是不是少了data的长度呀?
 楼主| 发表于 2009-3-9 12:26 | 显示全部楼层

回复 6楼 xiaocheng_2007 的帖子

d=norm((X(:,i)-X(:,j)),inf);%计算相空间中每两点之间的距离,其中NORM(V,inf) = max(abs(V)),因此对时间序列来说,d几乎不可能为零;退一步说,就算使用heaviside函数,也不能避免你所说的情况,其实这两种用法是等价的。程序中之所以避免使用heaviside函数,是为了加快程序运行速度,当使用heaviside函数时,因为大量的函数调用,运行时间会大大增加,你可以试一下。
 楼主| 发表于 2009-3-9 12:30 | 显示全部楼层
请仔细看一下程序,程序中有“N=length(data);   %data的长度”就是求时间序列长度。这样做是为了节省函数调用时参数传递的时间。
发表于 2009-3-9 15:21 | 显示全部楼层
我的QQ50785952,交流下
 楼主| 发表于 2009-3-9 16:27 | 显示全部楼层

回复 10楼 xiaocheng_2007 的帖子

你好,我已经加了(QQ1015014470)
发表于 2009-3-12 21:07 | 显示全部楼层

回复 yuling帖子

Cs1(tau)=0;
                for ii=1:N_t-1
                    for jj=ii+1:N_t
                        d1=abs(Y(ii)-Y(jj));     % 计算状态空间中每两点之间的距离,取无穷范数
                        if r>d1
                            Cs1(tau)=Cs1(tau)+1;            
                        end
                    end
                end
                Cs1(tau)=2*Cs1(tau)/(N_t*(N_t-1));
***************************************************************************************
难道 这一步计算是省了不少时间。。。。。
我看过的程序在这是用 Cs1(tau)=correlation_integral_inf(Y,N_t,r)来实现的  采用上面的方法是少一些函数之间的调用,


ps: if Sdeltmean(i)<SDELTMEAN(I-1)&SDELTMEAN(I)<SDELTMEAN(I+1)
     
改为 :if Sdeltmean(i)<Sdeltmean(i-1)<Sdeltmean(i)<Sdeltmean(i+1)

最后:Elapsed time is 1317.287872 seconds.

tau_inf =

    10


tw_inf =

   191
似乎是比我以前的要快点

[ 本帖最后由 sandman 于 2009-3-12 21:18 编辑 ]
 楼主| 发表于 2009-3-12 21:21 | 显示全部楼层

回复 12楼 sandman 的帖子

我主要是减少子函数的调用,尤其是H函数
发表于 2009-3-12 21:57 | 显示全部楼层

回复 13楼 yuling 的帖子

哦 ,这个程序想直接在matlab里面缩短时间的确是有限,但是如果用vc写成dll格式的话,可能情况就大不一样了。  我是不懂vc不知道你有没兴趣?
 楼主| 发表于 2009-3-13 11:10 | 显示全部楼层

回复 14楼 sandman 的帖子

这的确能大幅度提高运行速度,因为我曾用简单的多重循环验证过这一点。不过我对VC不是多熟悉,只知道点入门知识。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-10 17:12 , Processed in 0.112816 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表