相空间重构
##更新:最近重新看了这部分的内容,发现dts=max(ds)-min(ds);有误,应该是用r最大的时候的s减去最小时候的s。改为dts = ds(end)-ds(1)。以及C2应该有个m次幂。还有每次计算r应该是子序列对应的领域半径。
对于相空间重构需要确定的嵌入维数(m)和时间延迟(
τ
τ
τ)
使用C-C法可同时确定m和
τ
τ
τ,也可使用求自相关系数当ACF的值低于
1
−
e
−
1
1-e ^{-1}
1−e−1,从而确定时间延迟。然后通过G-P算法确定m。由于使用G-P算法所需时间较长。故本文直接使用C-C法确定两个参数。
时间窗口
t
m
=
(
m
−
1
)
∗
τ
t_m=(m-1)*τ
tm=(m−1)∗τ
本文所给的函数是自己编写的,故大家可用于参考,对于相空间理论部分,懒得编写了,直接百度即可。
确定参数 t m t_m tm和 τ τ τ的主函数
clc;clear
% C-C方法的主程序
data=xlsread('all_data7.xlsx','all_data','b1:b20000');
[ACF,lags,bounds]=autocorr(data,100);
tau=59;
m=15;
[D,MEAN,STD]=zscore(data);
r1=0.5*STD;
N=length(data);
[a,b]=find(ACF<ACF(1)*(1-exp(-1)));
beifen=data;
%% 求统计量S
SS=[];
DTS=[];
COR=[];
T=[];
NN=1:6000;
for t=1:100
sums=0;
sumdts=0;
% if mod(N,t)~=0
% [a,b]=find(mod(NN,t)==0);
% data=beifen(1:b(end));
% end
for i=1:t
eval(['datat',num2str(i),'=[];'])
for j=1:N/t
eval(['tran=data(',num2str(i),'+',num2str(t),'*(',num2str(j),'-1));'])
eval(['datat',num2str(i),'=[datat',num2str(i),';tran];'])
end
end
for m=2:5
ds=[];
for k=1:4
S=[];
for i=1:t
eval(['S',num2str(i),'=[];'])
eval(['datat=datat',num2str(i),';'])
r=r1*k;
C1=C(m,N/t,r,t,datat);%%计算关联积分
C2=C(1,N/t,r,t,datat);
eval(['S',num2str(i),'=[S',num2str(i),',C1-C2];'])%% 计算检验统计量
eval(['S=[S;S',num2str(i),'];'])
end
S=sum(S)/t;
ds=[ds,S];
sums=sums+S;
end
dts=max(ds)-min(ds);%% 求ΔS
sumdts=sumdts+dts;
end
sums=sums/16;
sumdts=sumdts/4;
scor=sumdts+abs(sums);
T=[T;t];
% S=S/t;
SS=[SS;sums];
DTS=[DTS;sumdts];
COR=[COR;scor];
end
% xieru=["datat1","datat2","datat3","datat4","datat5"];
% for i=1:t
% eval(['tran=datat',num2str(i),';'])
% xlswrite('xkjcg.xlsx',tran,xieru(i));
% end
tt=1:100;
plot(tt,SS,tt,DTS,tt,COR)
legend('SS','DTS','COR')
C子程序
%% 相空间重构确定嵌入维数和时间延迟
function [C]=C(m,N,r,t,data)
%% 将原始数据 重构相空间
M=N-(m-1)*t;%% M为重构相空间中相点的个数
x=[];%x用于保存相空间重构的结果
for i=1:M
for j=1:m
x(i,j)=[data(i+(j-1)*t)];
end
end
sumM=0;
for j=1:M
for i=1:j
dij=sum(abs(x(i,:)-x(j,:)));
if (r-dij)<0
th=0;
else
th=1;
end
sumM=sumM+th;
end
end
C=2*sumM/(M*(M-1));
end