三、排列熵
⇒ \Rightarrow ⇒点击此处——Python代码实现
\quad 优点:抗噪能力强、计算量小、所需数据长度较短
1、算法步骤
- 设信号时间序列为{X(i),i=1,2,…,n},进行相空间重构后,得到如下矩阵:
[ x ( 1 ) x ( 1 + τ ) ⋯ x ( 1 + ( m − 1 ) τ ) x ( 2 ) x ( 2 + τ ) ⋯ x ( 2 + ( m − 1 ) τ ) ⋮ ⋮ ⋱ ⋮ x ( k ) x ( k + τ ) ⋯ x ( k + ( m − 1 ) τ ) ] \left[ \begin{matrix} x(1) & x(1+\tau) & \cdots & x(1+(m-1)\tau)\\ x(2) & x(2+\tau) & \cdots & x(2+(m-1)\tau) \\ \vdots & \vdots & \ddots & \vdots \\ x(k) & x(k+\tau) & \cdots & x(k+(m-1)\tau) \end{matrix} \right] ⎣⎢⎢⎢⎡x(1)x(2)⋮x(k)x(1+τ)x(2+τ)⋮x(k+τ)⋯⋯⋱⋯x(1+(m−1)τ)x(2+(m−1)τ)⋮x(k+(m−1)τ)⎦⎥⎥⎥⎤
其中,m是嵌入维数, τ \tau τ是延时因子, k = n − ( m − 1 ) τ k=n-(m-1)\tau k=n−(m−1)τ,j=1,2,…,k。矩阵共有 k 个重构分量,每个重构分量有 m 维嵌入元素。 - 将矩阵中的第j个分量 x ( j ) , x ( j + τ ) , ⋯ , x ( j + ( m − 1 ) τ ) x(j),x(j+\tau),\cdots,x(j+(m-1)\tau) x(j),x(j+τ),⋯,x(j+(m−1)τ)按照数值的大小进行升序排列,得到:
x ( i + ( j 1 − 1 ) τ ) ≤ x ( i + ( j 2 − 1 ) τ ) ≤ ⋯ x ( i + ( j m − 1 ) τ ) x(i+(j_1-1)\tau) \leq x(i+(j_2-1)\tau) \leq \cdots x(i+(j_m-1)\tau) x(i+(j1−1)τ)≤x(i+(j2−1)τ)≤⋯x(i+(jm−1)τ)
其中 j 1 , j 2 , ⋯ , j m j_1,j_2,\cdots,j_m j1,j2,⋯,jm代表各元素在重构分量中的下标索引值。如果在重构分量中有两个或多个相等的值,比如 x ( i + ( j 1 − 1 ) τ ) = x ( i + ( j 2 − 1 ) τ ) x(i+(j_1-1)\tau)=x(i+(j_2-1)\tau) x(i+(j1−1)τ)=x(i+(j2−1)τ)时,则需要根据 j 1 、 j 2 j_1、j_2 j1、j2的大小来进行排序。满足 j 1 < j 2 j_1<j_2 j1<j2即可,此时有: x ( i + ( j 1 − 1 ) τ ) ≤ x ( i + ( j 2 − 1 ) τ ) x(i+(j_1-1)\tau) \leq x(i+(j_2-1)\tau) x(i+(j1−1)τ)≤x(i+(j2−1)τ) - 每个重构分量都可以得到一个重构符号序列:
S ( l ) = ( j 1 , j 2 , ⋯ . j m ) S(l)=(j_1,j_2,\cdots.j_m) S(l)=(j1,j2,⋯.jm)
其中,l=1.2. ⋯ \cdots ⋯,k,满足 k ≤ m ! k \leq m! k≤m! 。每个重构分量是 m 维空间,映射到 m 维符号序列,共有 m!种排列方式。 - 计算每一种m维符号序列的概率 P 1 , P 2 , ⋯ , P K P_1,P_2,\cdots,P_K P1,P2,⋯,PK,那么,根据香浓熵的定义,信号时间序列 X(i)的排列熵(PE)定义为:
H p ( m ) = − ( ∑ j = 1 k P j l n P j ) H_p(m) = -(\displaystyle \sum^k_{j=1}P_jlnP_j) Hp(m)=−(j=1∑kPjlnPj)
其中,当 P j = 1 / m ! P_j=1/m! Pj=1/m!时,排列熵达到最大值ln(m!)。实际的处理过程中,通常将 H p ( m ) H_p(m) Hp(m)进行归一化的处理:
0 ≤ H p = H p / l n ( m ! ) ≤ 1 0 \leq H_p=H_p/ln(m!)\leq1 0≤Hp=Hp/ln(m!)≤1
\quad\quad
\quad\quad 排列熵 H p H_p Hp 的大小衡量信号时间序列的随机变化程度, H p H_p Hp 的值越大,表示信号时间序列越随机,信号越复杂;反之,则说明信号序列越规则,复杂度较小。
\quad\quad
2、核心代码
- 排列熵算法
def func(n):"""求阶乘"""if n == 0 or n == 1:return 1else:return (n * func(n - 1))def compute_p(S):"""计算每一种 m 维符号序列的概率"""_map = {}for item in S:a = str(item)if a in _map.keys():_map[a] = _map[a] + 1else:_map[a] = 1freq_list = []for freq in _map.values():freq_list.append(freq / len(S))return freq_listdef Permutation_Entropy(x, m, t):"""计算排列熵值"""length = len(x) - (m-1) * t# 重构 k*m 矩阵y = [x[i:i+m*t:t] for i in range(length)]# 将各个分量升序排序S = [np.argsort(y[i]) for i in range(length)]# 计算每一种 m 维符号序列的概率freq_list = compute_p(S)# 计算排列熵pe = 0for freq in freq_list:pe += (- freq * np.log(freq))return pe / np.log(func(m))
- 读取数据集
# 自定义函数计算排列熵
def solve_3():calm = []stress = []print("calm:")for i in range(7):with open('F:\EEG\paper1\calm\s01\_' + channel_names[i] + '.npy', 'rb') as file:sub = np.load(file, encoding='latin1', allow_pickle=True)data = sub[1][0][0:1000]ep = Permutation_Entropy(data, 5, 10)calm.append(ep)print(ep)print("stress:")for i in range(7):with open('F:\EEG\paper1\stress\s01\_' + channel_names[i] + '.npy', 'rb') as file:sub = np.load(file, encoding='latin1', allow_pickle=True)data = sub[2][0][0:1000]ep = Permutation_Entropy(data, 5, 10)stress.append(ep)print(ep)return calm,stress
calm_5, stress_5 = solve_3()
- 画图
plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号y = calm_5 # sub[1]
z = stress_5 # sub[2]
fig = plt.figure(figsize=(8, 6))labels = ['FP1', 'F3', 'F7', 'FP2', 'FZ', 'F4', 'F8']plt.xlabel("导联", fontsize=14)
plt.ylabel("排列熵", fontsize=14)plt.plot(labels, y, c='b', label="平静")
plt.plot(labels, z, c='r', label="压力")
plt.tick_params(axis='both', labelsize=14)
plt.legend()
plt.show()
\quad\quad
3、参考文献
[1] 苏建新. 基于脑电信号的情绪识别研究[D].
本文链接:https://my.lmcjl.com/post/18723.html
展开阅读全文
4 评论