相空间重构

##更新：最近重新看了这部分的内容，发现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

这个函数是用互信息法求时延的，好像是找到第一个极小值就是最佳时延

function [mi] = VectorMI(data,delta_t1,N)
% 该函数用来计算一个向量之间的互信息随时延的变化
% delta_t1是计算互信息的最大时延
% N是参与计算的点数

for j=1:delta_t1
    u1=data(1:N);
    u2=data(1+j:N+j);

    wind_size=floor(power(length(u1),1/3)+0.5);%floor向负无穷取整，power幂运算

    x = [u1, u2];
    n = wind_size;
    [xrow, xcol] = size(x);
    bin = zeros(xrow,xcol);
    pmf = zeros(n, 2);

    for i = 1:2
        minx = min(x(:,i));
        maxx = max(x(:,i));
        binwidth = (maxx - minx) / n;
        edges = minx + binwidth*(0:n);
        histcEdges = [-Inf edges(2:end-1) Inf];
        [occur,bin(:,i)] = histc(x(:,i),histcEdges,1);
        pmf(:,i) = occur(1:n)./xrow;
    end

    % Calculate the joint probability density of u1 and u2
    jointOccur = accumarray(bin,1,[n,n]);

    %(xi?yi) The number of two data falling into the n * n aliquot is the joint probability density
    jointPmf = jointOccur./xrow;
    Hx = -(pmf(:,1))'*log2(pmf(:,1)+eps);
    Hy = -(pmf(:,2))'*log2(pmf(:,2)+eps);
    Hxy = -(jointPmf(:))'*log2(jointPmf(:)+eps);
    MI = Hx+Hy-Hxy;
    mi(j) = MI/sqrt(Hx*Hy);
end

这个函数用cao法来求最佳嵌入维，网上找来自己修改的，反正是能用，原理还挺复杂。

function cao_m(data,min_m,max_m,tau)
% 该程序计算时间序列的嵌入维数
% data为原始数据，是列向量
% min_m,max_m分别为最小和最大嵌入维数
% tau为时间延迟
% 计算值会稳定在某个值附近，开始稳定的位置就是嵌入维数
 
N=length(data);
k=1;
 
for m=min_m:max_m
    Y=phasespace(data,m,tau);%相重构
    Y=Y';
    [m,M]=size(Y);
    for i=1:N-m*tau
        %-----计算第i个向量与每个向量的距离存于d(j)中--------
        for j=1:N-m*tau
            d(j)=norm(Y(:,i)-Y(:,j),inf);
        end
        %-----求距离最短向量的下标-----------
        temp=sort(d);
        D(i,1)=i;      %D的第一列为向量序号
        temp1=find(temp>0);
        temp2=find(d==temp(temp1(1)));
        D(i,2)=temp2(1);  %第二列为与之对应的最短距离向量的序号
        D(i,3)=temp(temp1(1));%第三列为与之对应的最短距离
        %-----------计算a(i,m)-----------------------
        Y1=[Y(:,i);data(m*tau+i)];
        Y2=[Y(:,D(i,2));data(D(i,2)+m*tau)];
        ad(i)=norm(Y1-Y2,inf)/D(i,3);
        clear d Y1 Y2 temp temp1 temp2
    end
     %---------------求E(d)-----------------
    E(k,1)=m;
    E(k,2)=sum(ad)/(N-m*tau);
    %---------------求E*(d)-----------------
    En(k,1)=m;
    En1(k,1)=m;
    for kk=1:N-m*tau
         dd(kk)=abs(data(D(kk,1)+m*tau)-data(D(kk,2)+m*tau));
         %dd1(kk)=abs(data(D(kk,1)+1)-data(D(kk,2)+1));
    end
    En(k,2)=sum(dd)/(N-m*tau);
    k=k+1;
 
    clear D
end
 
%-----------求E1(d)-------------------
for i=1:(max_m-min_m)
    E1(i,1)=E(i,1);
    E1(i,2)=E(i+1,2)/E(i,2);
end
%-----------求E2(d)-------------------
for i=1:(max_m-min_m)
    E2(i,1)=En(i,1);
    E2(i,2)=En(i+1,2)/En(i,2);
end
%-----------画图------------------------
figure(1)
plot(E1(:,1),E1(:,2),'-bs',E2(:,1),E2(:,2),'-r*');xlabel('维数');ylabel('E1(d)&E2(d)');
grid on