相空间重构

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

Python实现相空间重构求关联维数——GP算法、自相关法求时间延迟tau、最近邻算法求嵌入维数m

GP算法：

若有一维时间序列为{x1,x2,…,xn},对其进行相空间重构得到高维相空间的一系列向量：

$x_i(\tau ,m)=\left(x_i,x_{i1},\cdots ,x_{i+{(m-1)}_{\tau }}\right)$

式中：$\tau$为时间延迟，$\tau$=k$\mathrm{\Delta }t$，其中k为整数，为采样时间间隔；m为嵌入维数；i=1，2，⋯，N；N为重构后向量的个数，$N=n-(m-1)\tau$。
重构相空间关联维数为：

$D_2=\underset{r\to 0}{\lim }\frac{\ln c_r}{\ln r}$

$c_r=\frac{1}{N^2}$$\sum \sum H$$\left(r-||x_j-x_k||\right)$

式中：j≠k；r为m维超球半径；H为Heaviside函数。

def GP(imf,tau):            #GP算法求关联维数
    N=2000
    if (len(imf) != N):
        print('请输入指定的数据长度！')   # N为指定数据长度
        return
    elif (isinstance(imf, np.ndarray) != True):
        print('数据格式错误！')
        return
    else:
        m_max=10                  #最大嵌入维数
        ss=50                     #r的步长
        fig=plt.figure()
        for m in range(1,m_max+1):
            i_num = N - (m - 1) * tau
            kj_m = np.zeros((i_num, m))  # m维重构相空间
            for i in range(i_num):
                for j in range(m):
                    kj_m[i][j] = imf[i + j * tau]
            dist_min, dist_max = np.linalg.norm(kj_m[0] - kj_m[1]), np.linalg.norm(kj_m[0] - kj_m[1])
            Dist_m = np.zeros((i_num, i_num))  # 两向量之间的距离
            for i in range(i_num):
                for k in range(i_num):
                    D= np.linalg.norm(kj_m[i] - kj_m[k])
                    if(D>dist_max):
                        dist_max=D
                    elif(D>0 and D<dist_min):
                        dist_min=D
                    Dist_m[i][k] = D
            dr=(dist_max-dist_min)/(ss-1)           #r的间距
            r_m=[]
            Cr_m=[]
            for r_index in range(ss):
                r=dist_min+r_index*dr
                r_m.append(r)
                Temp=np.heaviside(r-Dist_m,1)
                for i in range(i_num):
                    Temp[i][i]=0
                Cr_m.append(np.sum(Temp))
            r_m=np.log(np.array((r_m)))
            Cr_m=np.log(np.array((Cr_m))/(i_num*(i_num-1)))
            plt.plot(r_m,Cr_m)
        plt.show()

自相关法确定$\tau$：

计算时间序列{x1,x2,…,xn}的自相关函数：

$R(j\tau )=\frac{1}{N}\sum$$x(i)x(i+j\tau )$

当自相关函数值下降到初始函数值的1-${\mathrm{e}}^{-1}$时。所对应的$\tau$即为时间延迟参数。

# 计算GP算法的时间延迟参数（自相关法）
def get_tau(imf):
    N=2000
    if (len(imf) != N):
        print('请输入指定的数据长度！')  # N为指定数据长度
        return 0
    elif (isinstance(imf, np.ndarray) != True):
        print('数据格式错误！')
        return 0
    else:
        j = 1  # j为固定值
        tau_max = 20
        Rall = np.zeros(tau_max)
        for tau in range(tau_max):
            R = 0
            for i in range(N - j * tau):
                R += imf[i] * imf[i + j * tau]
            Rall[tau] = R / (N - j * tau)
        for tau in range(tau_max):
            if Rall[tau] < (Rall[0] * 0.6321):
                break
        return tau

假近邻算法确定m

对m维相空间每一个向量$X_{i(m)}=\left\{x_i,x_{i+\tau },\cdots ,x_{i+(m-1)\tau }\right\}$,i=1,2,…,N,N为向量总数，找出它的最近向量$X_{j(m)}$，计算两者欧氏距离$R_m(i)=||X_{i(m)}-X_{j(m)}||$，它们在m+1维空间的距离为：

$R_{m+1}(i)=||X_{i(m+1)}-X_{j(m+1)}||$

如果$R_{m+1}(i)$>>$R_m(i)$,则为虚假近邻点，定义比值：

$R(i)=$$\sqrt{\frac{{\left[R_{m+1}(i)\right]}^2-{\left[R_m(i)\right]}^2}{{\left[R_m(i)\right]}^2}}$

若$R(i)>R_0$,则称$X_j$为$X_i$的假近邻点，$R_0$为阈值通常取大于10.计算该m下虚假近邻点占点比例，直到虚假近邻点百分比很小或不随m增大而减少时，此时的m即为所需嵌入维数。

#计算GP算法的嵌入维数(假近邻算法)
def get_m(imf, tau):
    N=2000
    if (len(imf) != N):
        print('请输入指定的数据长度！')  # N为指定数据长度
        return 0, 0
    elif (isinstance(imf, np.ndarray) != True):
        print('数据格式错误！')
        return 0, 0
    else:
        m_max = 10
        P_m_all = []  # m_max-1个假近邻点百分率
        for m in range(1, m_max + 1):
            i_num = N - (m - 1) * tau
            kj_m = np.zeros((i_num, m))  # m维重构相空间
            for i in range(i_num):
                for j in range(m):
                    kj_m[i][j] = imf[i + j * tau]
            if (m > 1):
                index = np.argsort(Dist_m)
                a_m = 0  # 最近邻点数
                for i in range(i_num):
                    temp = 0
                    for h in range(i_num):
                        temp = index[i][h]
                        if (Dist_m[i][temp] > 0):
                            break
                    D = np.linalg.norm(kj_m[i] - kj_m[temp])
                    D = np.sqrt((D * D) / (Dist_m[i][temp] * Dist_m[i][temp]) - 1)
                    if (D > 10):
                        a_m += 1
                P_m_all.append(a_m / i_num)
            i_num_m = i_num - tau
            Dist_m = np.zeros((i_num_m, i_num_m))  # 两向量之间的距离
            for i in range(i_num_m):
                for k in range(i_num_m):
                    Dist_m[i][k] = np.linalg.norm(kj_m[i] - kj_m[k])
        P_m_all = np.array(P_m_all)
        m_all = np.arange(1, m_max)
        return m_all, P_m_all

三连、三连、三连

完整测试代码如下：

import numpy as np
from scipy.fftpack import fft
from scipy import fftpack
import matplotlib.pyplot as plt

N_ft=2000         #时频域的点数

# 计算GP算法的时间延迟参数（自相关法）
def get_tau(imf):
    if (len(imf) != N_ft):
        print('请输入指定的数据长度！')  # 需要更改，比如弹出对话框
        return 0,0,0
    elif (isinstance(imf, np.ndarray) != True):
        print('数据格式错误！')
        return 0,0,0
    else:
        j = 1  # j为固定值
        tau_max = 20
        Rall = np.zeros(tau_max)
        for tau in range(tau_max):
            R = 0
            for i in range(N_ft - j * tau):
                R += imf[i] * imf[i + j * tau]
            Rall[tau] = R / (N_ft - j * tau)
        for tau in range(tau_max):
            if Rall[tau] < (Rall[0] * 0.6321):
                break
            tauall=np.arange(tau_max)
        return tauall,Rall,tau


# 计算GP算法的嵌入维数(假近邻算法)
def get_m(imf, tau):
    if (len(imf) != N_ft):
        print('请输入指定的数据长度！')  # 需要更改，比如弹出对话框
        return 0, 0
    elif (isinstance(imf, np.ndarray) != True):
        print('数据格式错误！')
        return 0, 0
    else:
        m_max = 10
        P_m_all = []  # m_max-1个假近邻点百分率
        for m in range(1, m_max + 1):
            i_num = N_ft - (m - 1) * tau
            kj_m = np.zeros((i_num, m))  # m维重构相空间
            for i in range(i_num):
                for j in range(m):
                    kj_m[i][j] = imf[i + j * tau]
            if (m > 1):
                index = np.argsort(Dist_m)
                a_m = 0  # 最近邻点数
                for i in range(i_num):
                    temp = 0
                    for h in range(i_num):
                        temp = index[i][h]
                        if (Dist_m[i][temp] > 0):
                            break
                    D = np.linalg.norm(kj_m[i] - kj_m[temp])
                    D = np.sqrt((D * D) / (Dist_m[i][temp] * Dist_m[i][temp]) - 1)
                    if (D > 10):
                        a_m += 1
                P_m_all.append(a_m / i_num)
            i_num_m = i_num - tau
            Dist_m = np.zeros((i_num_m, i_num_m))  # 两向量之间的距离
            for i in range(i_num_m):
                for k in range(i_num_m):
                    Dist_m[i][k] = np.linalg.norm(kj_m[i] - kj_m[k])
        P_m_all = np.array(P_m_all)
        m_all = np.arange(1, m_max)
        return m_all, P_m_all


# GP算法求关联维数(时频域特征)
def GP(imf, tau):
    if (len(imf) != N_ft):
        print('请输入指定的数据长度！')  # 需要更改，比如弹出对话框
        return
    elif (isinstance(imf, np.ndarray) != True):
        print('数据格式错误！')
        return
    else:
        m_max = 10  # 最大嵌入维数
        ss = 50  # r的步长
        fig = plt.figure(1)
        for m in range(1, m_max + 1):
            i_num = N_ft - (m - 1) * tau
            kj_m = np.zeros((i_num, m))  # m维重构相空间
            for i in range(i_num):
                for j in range(m):
                    kj_m[i][j] = imf[i + j * tau]
            dist_min, dist_max = np.linalg.norm(kj_m[0] - kj_m[1]), np.linalg.norm(kj_m[0] - kj_m[1])
            Dist_m = np.zeros((i_num, i_num))  # 两向量之间的距离
            for i in range(i_num):
                for k in range(i_num):
                    D = np.linalg.norm(kj_m[i] - kj_m[k])
                    if (D > dist_max):
                        dist_max = D
                    elif (D > 0 and D < dist_min):
                        dist_min = D
                    Dist_m[i][k] = D
            dr = (dist_max - dist_min) / (ss - 1)  # r的间距
            r_m = []
            Cr_m = []
            for r_index in range(ss):
                r = dist_min + r_index * dr
                r_m.append(r)
                Temp = np.heaviside(r - Dist_m, 1)
                for i in range(i_num):
                    Temp[i][i] = 0
                Cr_m.append(np.sum(Temp))
            r_m = np.log(np.array((r_m)))
            Cr_m = np.log(np.array((Cr_m)) / (i_num * (i_num - 1)))
            plt.plot(r_m, Cr_m)
        plt.xlabel('ln(r)')
        plt.ylabel('ln(C)')
        plt.show()

if __name__=='__main__':
    # 检验关联维数程序
    t = []
    f1 = 25
    f2 = 30
    for i in range(N_ft):
        t.append(i * 0.001)
    t = np.array(t)
    # yu = np.ones(M * N)
    AEall = np.sin(t * 2 * np.pi * f1) + np.sin(t * 2 * np.pi * f2)  #在这里直接改信号


    tauall, Rall, tau = get_tau(AEall)
    m, P = get_m(AEall, tau)
    GP(AEall, 1)
    print(tau)
    fig2 = plt.figure(2)
    yu = np.ones(len(tauall)) * Rall[0] * 0.6321
    plt.plot(tauall, Rall)
    plt.plot(tauall, yu)
    plt.xlabel('tau')
    plt.ylabel('R')
    plt.show()
    fig3 = plt.figure(3)
    plt.plot(m, P)
    plt.xlabel('m')
    plt.ylabel('P')
    plt.show()

闲话几句

相空间重构听起来高大上，但是一旦确定了嵌入维数和延迟时间，之后的操作就十分简单了。

下边给个代码示例一下“之后的操作”，至于嵌入维数和延迟时间的确定，还是需要查查论文的，后期看情况再写。

python程序

def embed_vectors_1d(self, lag, embed):
    """Embeds vectors from a one dimensional time series in m-dimensional
    space.

    Parameters
    ----------
    X : 1d array
        Training or testing set.
    lag : int
        Lag value as calculated from the first minimum of the mutual info.
    embed : int
        Embedding dimension. How many lag values to take.
    predict : int
        Distance to forecast (see example).

    Returns
    -------
    features : 2d array
        Contains all of the embedded vectors. Shape (num_vectors,embed).

    Example
    -------
    >>> X = [0,1,2,3,4,5,6,7,8,9,10]
    em = 3
    lag = 2
    predict=3

    >>> embed_vectors_1d
    features = [[0,2,4], [1,3,5], [2,4,6], [3,5,7]]
    """

    tsize = self.X.shape[0]
    t_iter = tsize-(lag*(embed-1))

    features = np.zeros((t_iter,embed))

    for ii in range(t_iter):

        end_val = ii+lag*(embed-1)+1

        part = self.X[ii : end_val]

        features[ii,:] = part[::lag]

    return features

这个程序是从skccm里边摘抄出来的，程序中用了for循环，这个在效率上不太高。

朴素的解释

相空间重构是要解决什么问题呢？最近看的两篇网文^{1 2}很好的解释了这个问题，先把图片粘贴到这里，有人想看的话我再加上文字解释。

see also

1. skccm：https://github.com/NickC1/skccm（GitHub可能无法打开）