二、样本熵
⇒ \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
)
C_i^m(r)
Cim(r),即统计 D{y(i),y(j)}<r 的个数
N
m
(
i
)
N^m(i)
Nm(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 ) N^m(i) Nm(i)也叫做模板匹配数,其中 j<=N-m。并且对所有的 C i m ( r ) C_i^m(r) Cim(r)求平均值:
B m ( r ) = 1 N − m ( ∑ i = 1 N − m C i m ( r ) ) B^m(r) = {1 \over {N-m}} (\displaystyle \sum^{N-m}_{i=1} C_i^m(r)) Bm(r)=N−m1(i=1∑N−mCim(r)) - 将嵌入维数设为 m+1,重复步骤(1)到(3),同理得到
B
m
+
1
(
r
)
B^{m+1}(r)
Bm+1(r)。所以可以得到该序列的样本熵为:
S a m p E n ( m , r ) = lim N → ∞ [ − l n B m + 1 ( r ) B m ( r ) ] SampEn(m,r) = \displaystyle \lim_{N \rightarrow \infty}[-ln{{B^{m+1}(r)} \over {B^m(r)}]} SampEn(m,r)=N→∞lim[−lnBm(r)]Bm+1(r)
若 N 为有限值时,样本熵可为:
S a m p E n ( m , r ) = − l n B m + 1 ( r ) B m ( r ) SampEn(m,r) = -ln{{B^{m+1}(r)} \over {B^m(r)}} SampEn(m,r)=−lnBm(r)Bm+1(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)