相空间重构matlab实现

相空间重构

##更新：最近重新看了这部分的内容，发现dts=max(ds)-min(ds);有误，应该是用r最大的时候的s减去最小时候的s。改为dts = ds(end)-ds(1)。以及C2应该有个m次幂。还有每次计算r应该是子序列对应的领域半径。

对于相空间重构需要确定的嵌入维数（m）和时间延迟($\tau$)
使用C-C法可同时确定m和$\tau$,也可使用求自相关系数当ACF的值低于$1-e^{-1}$，从而确定时间延迟。然后通过G-P算法确定m。由于使用G-P算法所需时间较长。故本文直接使用C-C法确定两个参数。
时间窗口$t_m=(m-1)\ast \tau$
本文所给的函数是自己编写的，故大家可用于参考，对于相空间理论部分，懒得编写了，直接百度即可。

确定参数$t_m$和$\tau$的主函数

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