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

GP算法：

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

xi(τ,m)=(xi,xi1,⋯,xi+(m−1)τ){x_i}(\tau ,m) = \left( {{x_i},{x_{i1}}, \cdots ,{x_{i + {{(m - 1)}_\tau }}}} \right)xi(τ,m)=(xi,xi1,⋯,xi+(m−1)τ)

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

D2=lim⁡r→0ln⁡crln⁡r{D_2} = \mathop {\lim }\limits_{r \to 0} \frac{{\ln {c_r}}}{{\ln r}}D2=r→0limlnrlncr

cr=1N2{c_r} = \frac{1}{{{N^2}}}cr=N21∑∑H\sum\sum H∑∑H(r−∣∣xj−xk∣∣)\left( {r - ||{x_j} - {x_k}||} \right)(r−∣∣xj−xk∣∣)

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

def GP(imf,tau):            #GP算法求关联维数N=2000if (len(imf) != N):print('请输入指定的数据长度！')   # N为指定数据长度returnelif (isinstance(imf, np.ndarray) != True):print('数据格式错误！')returnelse:m_max=10                  #最大嵌入维数ss=50                     #r的步长fig=plt.figure()for m in range(1,m_max+1):i_num = N - (m - 1) * taukj_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=Delif(D>0 and D<dist_min):dist_min=DDist_m[i][k] = Ddr=(dist_max-dist_min)/(ss-1)           #r的间距r_m=[]Cr_m=[]for r_index in range(ss):r=dist_min+r_index*drr_m.append(r)Temp=np.heaviside(r-Dist_m,1)for i in range(i_num):Temp[i][i]=0Cr_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τ)=1N∑R(j\tau )= \frac{1}{{{N}}}\sumR(jτ)=N1∑x(i)x(i+jτ)x(i)x(i + j\tau )x(i)x(i+jτ)

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

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

假近邻算法确定m

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

Rm+1(i)=∣∣Xi(m+1)−Xj(m+1)∣∣{R_{m + 1}}(i) = ||{X_{i(m + 1)}} - {X_{j(m + 1)}}||Rm+1(i)=∣∣Xi(m+1)−Xj(m+1)∣∣

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

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

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

#计算GP算法的嵌入维数(假近邻算法)
def get_m(imf, tau):N=2000if (len(imf) != N):print('请输入指定的数据长度！')  # N为指定数据长度return 0, 0elif (isinstance(imf, np.ndarray) != True):print('数据格式错误！')return 0, 0else:m_max = 10P_m_all = []  # m_max-1个假近邻点百分率for m in range(1, m_max + 1):i_num = N - (m - 1) * taukj_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 = 0for h in range(i_num):temp = index[i][h]if (Dist_m[i][temp] > 0):breakD = 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 += 1P_m_all.append(a_m / i_num)i_num_m = i_num - tauDist_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 pltN_ft=2000         #时频域的点数# 计算GP算法的时间延迟参数（自相关法）
def get_tau(imf):if (len(imf) != N_ft):print('请输入指定的数据长度！')  # 需要更改，比如弹出对话框return 0,0,0elif (isinstance(imf, np.ndarray) != True):print('数据格式错误！')return 0,0,0else:j = 1  # j为固定值tau_max = 20Rall = np.zeros(tau_max)for tau in range(tau_max):R = 0for 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):breaktauall=np.arange(tau_max)return tauall,Rall,tau# 计算GP算法的嵌入维数(假近邻算法)
def get_m(imf, tau):if (len(imf) != N_ft):print('请输入指定的数据长度！')  # 需要更改，比如弹出对话框return 0, 0elif (isinstance(imf, np.ndarray) != True):print('数据格式错误！')return 0, 0else:m_max = 10P_m_all = []  # m_max-1个假近邻点百分率for m in range(1, m_max + 1):i_num = N_ft - (m - 1) * taukj_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 = 0for h in range(i_num):temp = index[i][h]if (Dist_m[i][temp] > 0):breakD = 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 += 1P_m_all.append(a_m / i_num)i_num_m = i_num - tauDist_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('请输入指定的数据长度！')  # 需要更改，比如弹出对话框returnelif (isinstance(imf, np.ndarray) != True):print('数据格式错误！')returnelse:m_max = 10  # 最大嵌入维数ss = 50  # r的步长fig = plt.figure(1)for m in range(1, m_max + 1):i_num = N_ft - (m - 1) * taukj_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 = Delif (D > 0 and D < dist_min):dist_min = DDist_m[i][k] = Ddr = (dist_max - dist_min) / (ss - 1)  # r的间距r_m = []Cr_m = []for r_index in range(ss):r = dist_min + r_index * drr_m.append(r)Temp = np.heaviside(r - Dist_m, 1)for i in range(i_num):Temp[i][i] = 0Cr_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 = 25f2 = 30for 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.6321plt.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实现相空间重构求关联维数——GP算法、自相关法求时间延迟tau、最近邻算法求嵌入维数m相关推荐

相空间重构中延迟时间tau的选择：自相关法（matlab实现）
相空间重构技术中有两个关键参数:延迟时间tau和嵌入维数m 使用自相关法确定延迟时间tau,主要思想在于通过计算原信号的自相关函数,找到自相关函数值下降到初始值R(0)时的1-1/e倍时的延迟时间ta ...
相空间重构与几个常用非线性模型实现
本文主要介绍了Lorenz模型.Rossler模型.Logistic映射以及其代码实现 1.Lorenz模型 %sy.m t = -10*pi:pi/250:10*pi; comet3((cos(2* ...
相空间重构 matlab 程序源,matlab求相空间重构延迟时间和嵌入维数
关联积分计算 function C_I=correlation_integral(X,M,r) %该函数用来计算关联积分 %C_I:关联积分的返回值 %X:重构的相空间矢量,是一个m*M的矩阵 %M: ...
混沌相空间重构中求时延和嵌入维的matlab代码
这个函数是用互信息法求时延的,好像是找到第一个极小值就是最佳时延 function [mi] = VectorMI(data,delta_t1,N) % 该函数用来计算一个向量之间的互信息随时延的变化 ...
相空间重构matlab实现
相空间重构 ##更新:最近重新看了这部分的内容,发现dts=max(ds)-min(ds);有误,应该是用r最大的时候的s减去最小时候的s.改为dts = ds(end)-ds(1).以及C2应该有个 ...
时间序列模型之相空间重构
一般的时间序列主要是在时间域中进行模型的研究,而对于混沌时间序列,无论是混沌不变量的计算,混沌模型的建立和预测都是在所谓的相空间中进行,因此相空间重构就是混沌时间序列处理中非常重要的一个步骤.所谓混沌 ...
混沌性时间序列的分析方法：EEMD+相空间重构
一.引言上一篇文章中,我们理解了混沌理论的发展.定义以及特点. 接下来,要结合我的研究方向,在机械振动时间序列信号的基础上,做出故障的诊断和预判. 时间序列也是结构化的数据,每一个时间戳下就有一个值 ...
基于相空间重构的混沌背景下微弱信号检测方法仿真
1.1算法参数取值对系统性能的影响在研究算法性能之前,首先需要分析各个参数对算法整体性能的影响,本文将重点考虑相空间重构参数和m,SVM支持向量机参数C和.这里分别对四个参数进行性能影响测试,首先对 ...
python机器学习案例系列教程——关联分析（Apriori、FP-growth）
全栈工程师开发手册 (作者:栾鹏) python数据挖掘系列教程关联分析的基本概念关联分析(Association Analysis):在大规模数据集中寻找有趣的关系. 频繁项集(Frequent ...

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

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

GP算法：

自相关法确定τ\tauτ：

假近邻算法确定m

三连、三连、三连

完整测试代码如下：

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

最新文章

热门文章