一、近似熵
⇒
\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
)
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+1。 - 对每个求得的
C
i
m
(
r
)
C_i^m(r)
Cim(r),计算其对数值,并求对数的平均值,如公式所示:
ϕ m ( i ) = ( ∑ i = 1 N − M + 1 l n C i m ( r ) ) / ( N − m + 1 ) \phi^m(i) = (\displaystyle \sum^{N-M+1}_{i=1}{lnC_i^m(r)})/(N-m+1) ϕm(i)=(i=1∑N−M+1lnCim(r))/(N−m+1) - 将嵌入维数变为 m+1,重复步骤(1)到(4),得到
C
i
m
(
r
)
C_i^m(r)
Cim(r) 和
ϕ
m
+
1
(
i
)
\phi^{m+1}(i)
ϕm+1(i)。所以该序列的近似熵可以表示为如下公式:
A p E n ( m , r ) = ϕ m ( i ) − ϕ m + 1 ( i ) ApEn(m,r) = \phi^m(i) - \phi^{m+1}(i) ApEn(m,r)=ϕm(i)−ϕ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算法之近似熵、互近似熵算法