脑电信号特征提取——近似熵

一、近似熵

$\Rightarrow$点击此处——Python代码实现
  \:   \: 若干个相似向量从 k 维空间增加至 k+1 维向量时仍然保持着较高的相似性的概率

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 是嵌入维数，是参与比较序列的长度，即窗口长度。
• 计算 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+1。
• 对每个求得的$C_i^m(r)$，计算其对数值，并求对数的平均值，如公式所示：
  $${\varphi }^m(i)=(\sum _{i=1}^{N-M+1}ln{C}_i^m(r))/(N-m+1)$$
• 将嵌入维数变为 m+1，重复步骤（1）到（4），得到$C_i^m(r)$ 和${\varphi }^{m+1}(i)$。所以该序列的近似熵可以表示为如下公式：
  $$ApEn(m,r)={\varphi }^m(i)-{\varphi }^{m+1}(i)$$

注：选取嵌入维数 m=2，有效阈值 r=0.15*std
  \:

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 _normalization(data):
    """将数据标准化，所有值减去平均值除以标准差"""
    _mean = np.mean(data)
    data_std = _std(data)
    return np.array([(x - _mean) / data_std for x in data])

def phi(x, m, r):
    """计算熵值"""
    length = len(x)-m+1
    # 重构m维向量
    y = [x[i:i+m] for i in range(length)]
     
    # 获取所有的比值列表
    C = [ len([1 for y_j in y if maxdist(y_i, y_j) < r]) / length for y_i in y ]
    
    for i in range(length):
        C[i] = np.log(C[i])
    return np.mean(np.sum(C))

def Approximate_Entropy(data, m):
    """计算近似熵"""
    x = data
    data_std = _std(x)
    r = data_std * 0.15
    return phi(x,m,r) - phi(x,m+1,r)

• 读取数据集

def solve_1():
    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 = Approximate_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 = Approximate_Entropy(data, 2)
            stress.append(ep)
            print(ep)
            
    return calm,stress

calm, stress = solve_1()

• 画图

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

y = calm
z = stress
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()

  \:
  \:

3、参考文献

[1] 苏建新. 基于脑电信号的情绪识别研究[D].
[2] 金子敦. 基于脑电信号的情绪特征提取与分类研究[D].
[3] python算法之近似熵、互近似熵算法