脑电信号特征提取——样本熵

二、样本熵

$\Rightarrow$点击此处——Python代码实现

是对近似熵的一种改进算法，是一种不同于近似熵而且不进行自身匹配的统计量方法

• 近似熵有两个缺点：

  • 近似熵在与自身匹配时具有偏差性；
  • 近似熵结果的一致性较差

• 样本熵具有如下特点：

  • 具有比时域统计（均值、方差等）更好的估计效果；
  • 对原始数据处理时无需进行粗粒化提取；
  • 可用于由确定信号与随机信号组成的混合信号。

1、算法步骤

• 设原始信号为{x(i),i=1,2,…,N}，按照下面公式重构出 m 维向量，用 y（i）表示，{ y(i),i=1,2,…,M,M=N-m+1},即：
  y ( i ) = { x ( i ) , x ( i + 1 ) , x ( i + 2 ) , ⋯   , x ( i + m − 1 ) } y(i)=\{x(i),x(i+1),x(i+2),\cdots ,x(i+m-1)\} y(i)={x(i),x(i+1),x(i+2),⋯,x(i+m−1)}
  其中， m 是选定的空间维度矢量，也即窗口长度，表示从第 i 个点开始的连续 m 个值。
• 计算 y(i)与 y(j)任意分量之间的欧式距离 d{y(i),y(j)}，并将各个分量之间最大的距离定义为最大贡献成分距离 D{y(i),y(j)}，如下式所示：
  D y ( i ) , y ( j ) = m a x { ∣ y ( i + k ) − x ( j + k ) ∣ } D{y(i),y(j)}=max\{|y(i+k)-x(j+k)|\} Dy(i),y(j)=max{∣y(i+k)−x(j+k)∣}
  其中，i，j = 1,2, ,N-m+1， k=0,1,…,m-1。
• 给定阈值 r（r>0），给定嵌入维数 m，计算代表序列{y(i)}规律性的概率大小量度$C_i^m(r)$，即统计 D{y(i),y(j)}<r 的个数$N^m(i)$与总数 N-m+1 的比例，也即概率大小，如公式所示：
  KaTeX parse error: Can't use function '$' in math mode at position 9: C_i^m(r)$̲ = $N^m(i) / (N…
  其中， i<=N-m。$N^m(i)$也叫做模板匹配数，其中 j<=N-m。并且对所有的$C_i^m(r)$求平均值：
  $$B^m(r)=\frac{1}{N-m}(\sum _{i=1}^{N-m}{C}_i^m(r))$$
• 将嵌入维数设为 m+1，重复步骤（1）到（3），同理得到$B^{m+1}(r)$。所以可以得到该序列的样本熵为:
  $$SampEn(m,r)=\underset{N\to \infty }{\lim }[-ln\frac{B^{m+1}(r)}{B^m(r)]}$$
  若 N 为有限值时，样本熵可为：
  $$SampEn(m,r)=-ln\frac{B^{m+1}(r)}{B^m(r)}$$

注：一般选取嵌入维数 m=2，有效阈值 r一般选 0.1~0.25 倍的原始数据标准差

2、核心代码

• 样本熵算法

def maxdist(x_i, x_j):
    """计算矢量之间的距离"""
    return np.max([np.abs(np.array(x_i) - np.array(x_j))])

def _std(data):
    """计算标准差"""
    if not isinstance(data, np.ndarray):
            data = np.array(data)
    return np.std(data)

def phi_2(x, m, r):
    """计算熵值"""
    length = len(x) - m + 1
    # 重构m维向量
    y = [x[i:i+m] for i in range(length)]
     
    # 获取所有的比值列表(不与自身匹配。不包括最后一个向量,注意是N-m不是N-m+1)
    C = [ len([1 for y_j in y[:-1] if maxdist(y_i, y_j) < r and (y_i != y_j).any()]) / (length-1) for y_i in y[:-1] ]
    
    ans = 0
    for i in range(length - 1):
        ans += C[i]
    return ans / (length - 1)

def Sample_Entropy(data, m):
    """计算样本熵"""
    x = data
    data_std = _std(x)
    r = data_std * 0.15
    return - np.log(phi_2(x,m+1,r) / phi_2(x,m,r))

• 读取数据集

# 使用自定义函数计算样本熵
def solve_2():
    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 = Sample_Entropy(data, 2)
            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[0][0][0:1000]
            ep = Sample_Entropy(data, 2)
            stress.append(ep)
            print(ep)
            
    return calm,stress

calm_3, stress_3 = solve_2()

• 画图

plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

y = calm_3  # sub[1]
z = stress_3
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

3、参考文献

[1] 苏建新. 基于脑电信号的情绪识别研究[D].
[2] 信号处理算法（2）：样本熵（SampEn）