c-cmethod 相空间重构

School

% 此程序用来测试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);

2.子函数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)

break;

end

end

% 寻找时间窗口tw：即Scor最小值对应的t

for i=1:length(Scor)

if Scor(i)==min(Scor)

tw=i;

break;

end

end

3.子函数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

4.子函数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

5.子函数4

function C=correlation_integral_inf(Y,M,r)

% 此函数用于计算关联积分,取无穷范数