数学与非线性科学 - 声振论坛 - 振动,动力学,声学,信号处理,故障诊断

cailiang 发表于 2009-7-4 09:50

请教时间延迟

这是从yuling给出的计算时间延迟的程序，我搬到机子上运行，不出结果，不知程序错哪了？还望高手指出。
%主程序
clear all
clear all
t0=0;
tf=130;
=ode45(@Lorenz,,[-1,0,1]);
Lorenz_data=x(10002:end,1);
max_d=200;
% 调用C_CMethod_inf，求tau和tw
tic
=C_CMethod_inf(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 dy = Lorenz(t,y)%方程的定义
a=16.0;
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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
function =C_CMethod_inf(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,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 C_I=correlation_integral(X,M,r)
%该函数用来计算计算关联积分
%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));%关联积分的值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 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

cailiang 发表于 2009-7-4 14:08

出结果了，yuling原来的程序里有几个小错误，修改了，可以的。呵呵

心灯发表于 2009-7-4 15:00

回复沙发 cailiang 的帖子

建议把修改后的程序在论坛上给出来，从而方便后来者，他们会感谢你的，呵呵。

ChaChing 发表于 2009-7-4 21:13

回复沙发 cailiang 的帖子

个人认同楼上的意见!
除此之外, 若能给出yuling原有程序的连接, 并在那亦回覆错误处更正, 那才真正方便后来者!

cailiang 发表于 2009-7-5 12:43

好的，只是中间有几个大小写的错误。
程序如下，只是转载，呵呵呵。
function =C_CMethod_inf(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,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 C_I=correlation_integral(X,M,r)
%该函数用来计算计算关联积分
%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));%关联积分的值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 dy = Lorenz(t,y)%方程的定义
a=16.0;
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;
=ode45(@Lorenz,,[-1,0,1]);
Lorenz_data=x(10002:end,1);
max_d=200;
% 调用C_CMethod_inf，求tau和tw
tic
=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 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

页: [1]

声振论坛's Archiver

请教时间延迟

回复 沙发 cailiang 的帖子

回复 沙发 cailiang 的帖子

回复沙发 cailiang 的帖子

回复沙发 cailiang 的帖子