三、排列熵

⇒ \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∑k​Pj​lnPj​)
  其中，当 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