本文是关于小波熵(Wavewlet Entropy, WE)的初步理解,由于笔者所能找到最早提出小波熵的论文是 Wavelet entropy: a new tool for analysis of short duration brain electrical signals,因此,本文便是以此论文为基础进行的说明。
1 预备知识
小波熵是在小波变换和信息熵的基础上形成的一套理论。
1.1 信息熵
信息熵(Information Entropy)最早由香农在他的论文A Mathematical Theory of Communication中提出,在Cover的《信息论基础》这本书中,对熵的定义如下:
设 X X X是一个离散型随机变量,其字母表(即概率论中的取值空间)为 X \mathcal{X} X,概率密度函数 p ( x ) = Pr ( X = x ) , x ∈ X p(x)=\Pr({X}=x),x\in\mathcal{X} p(x)=Pr(X=x),x∈X;
定义一个离散型随机变量 X {X} X的熵 H ( X ) H({X}) H(X)为
H ( X ) = − ∑ x ∈ X p ( x ) log p ( x ) H(X)=-\sum\limits_{x\in\mathcal{X}}{p(x)\log p(x)} H(X)=−x∈X∑p(x)logp(x)
同时书中定义离散型随机变量
X
{X}
X的取值
x
x
x自信息
I
(
x
)
=
−
log
p
(
x
)
I(x)=-\log p(x)
I(x)=−logp(x)
我们可以理解为离散型随机变量
X
{X}
X的熵
H
(
X
)
H({X})
H(X)为其各取值对应的互信息根据概率分布情况的加权平均值。
不难推断,当
p
(
x
)
→
1
p(x) \rightarrow 1
p(x)→1时,
I
(
x
)
→
0
I(x) \rightarrow 0
I(x)→0,当
p
(
x
)
→
0
p(x) \rightarrow 0
p(x)→0时,
I
(
x
)
→
∞
I(x) \rightarrow \infty
I(x)→∞,即当一个事件发生的概率越大时,其蕴含的信息量越小,反之,蕴含的信息量越大。
比如说:“明天太阳会从东边升起”和“明天太阳会从西边升起”这两句话的信息量是完全不同的,后者明显比前者更有价值(在听者相信这些话的前提下)。
同时,对于信息熵,我们可以推断出:当
x
x
x等概分布时,信息熵达到最大值,即
H
(
X
)
=
log
m
H(X)=\log m
H(X)=logm,其中
m
m
m为字母表
X
\mathcal{X}
X的取值个数。
以上是信息熵的定义,如涉及到其他拓展知识,在后续用到时会作说明。
1.2 小波变换
我们知道利用傅里叶变换可以分析一段信号的频域情况,但是无法知道各频率分量在时域上的分布情况,于是有人提出了小波变换(Wavelet Transform)。
小波变换要求引入合适的基,利用基在不同时间段上与信号进行相关性分析,根据相关性的强弱,分析得到信号在这一频段上时域的分布情况。
如果小波变换的基为正交基,则可以实现对任意信号的唯一分解,并且是可逆的。
一般要求小波函数是一个平滑且快速收敛的振荡函数,在时域和频域上都具有很好的局域性。定义小波族(Wavelet family)
ψ
a
,
b
\psi_{a,b}
ψa,b是由母小波
ψ
\psi
ψ伸缩和平移变换得到的函数集合:
ψ
a
,
b
(
t
)
=
∣
a
∣
−
1
/
2
ψ
(
t
−
b
a
)
\psi_{a,b}(t)=\begin{vmatrix}a\end{vmatrix}^{-1/2}\psi\biggl(\frac{t-b}{a}\biggr)
ψa,b(t)=
a
−1/2ψ(at−b)
其中
a
,
b
∈
R
,
a
≠
0
a,b\in\mathbb{R}, a\neq0
a,b∈R,a=0,
a
,
b
a,b
a,b分别为尺度参数和平移参数。
信号
S
(
t
)
S(t)
S(t)的连续小波变换(CWT)定义为函数
S
(
t
)
S(t)
S(t)与
ψ
a
,
b
\psi_{a,b}
ψa,b的相关性:
(
W
ψ
S
)
(
a
,
b
)
=
∣
a
∣
−
1
/
2
∫
−
∞
∞
S
(
t
)
ψ
∗
(
t
−
b
a
)
d
t
=
⟨
S
,
ψ
a
,
b
⟩
(W_{\psi}S)(a,b)=\begin{vmatrix}a\end{vmatrix}^{-1/2}\int_{-\infty}^{\infty}S(t)\psi^{*}\biggl(\frac{t-b}{a}\biggr)\mathrm{d}t=\langle S,\psi_{a,b}\rangle
(WψS)(a,b)=
a
−1/2∫−∞∞S(t)ψ∗(at−b)dt=⟨S,ψa,b⟩
因为在计算机运算中,2的幂次运算效率高,因此一般选取尺度参数
a
j
=
2
−
j
a_j=2^{-j}
aj=2−j,
b
j
,
k
=
2
−
j
k
b_{j,k}=2^{-j}k
bj,k=2−jk,其中
j
,
k
∈
Z
j,k \in \mathbb{Z}
j,k∈Z,因此小波函数可以写为:
ψ
j
,
k
(
t
)
=
2
j
/
2
ψ
(
2
j
t
−
k
)
j
,
k
∈
Z
\psi_{j,k}(t)=2^{j/2}\psi(2^{j}t-k)\quad j,k\in \mathbb{Z}
ψj,k(t)=2j/2ψ(2jt−k)j,k∈Z
与之相关地,还有离散小波变换(DWT)。
关于小波变换的进一步了解,可以参考B站UP主“有趣的理工男”的视频
一次讲透小波变换原理,全新角度切入,20分钟时长警告!
和
版本T0的离散小波变换解说,父小波母小波是什么?高低通滤波是怎么回事?时频图怎么画?具体计算原理又是什么?你的疑问在这都有答案!
2 新的定义
在说明新的定义之前,需要说明:本文参考文献使用的母小波是三次样条函数,如果对此不了解,可以参考此文章:样条函数。
假设信号
S
=
{
s
0
(
n
)
,
n
=
1
,
.
.
.
,
M
}
S=\{s_{0}(n), n=1, ..., M\}
S={s0(n),n=1,...,M}对应采样频率为
t
s
t_s
ts(一般取
M
M
M为
2
2
2的幂次),为简单起见,取
t
s
=
1
t_s=1
ts=1,若在所有分辨率上对信号进行分解,则分解次数
N
=
log
2
(
M
)
N=\log_2(M)
N=log2(M),且信号
S
(
t
)
S(t)
S(t)可展开为:
S
(
t
)
=
∑
j
=
−
N
−
1
∑
k
C
j
(
k
)
ψ
j
,
k
(
t
)
=
∑
j
=
−
N
−
1
r
j
(
t
)
S(t)=\sum\limits_{j=-N}^{-1}\sum\limits_{k}C_{j}(k)\psi_{j,k}(t)=\sum\limits_{j=-N}^{-1}r_{j}(t)
S(t)=j=−N∑−1k∑Cj(k)ψj,k(t)=j=−N∑−1rj(t)
其中小波系数
C
j
(
k
)
C_j(k)
Cj(k)可理解为在尺度
j
j
j和
j
+
1
j+1
j+1上信号逼近的局部残差,
r
j
(
t
)
r_j(t)
rj(t)是尺度
j
j
j上的残差信号,包含了频率范围
2
j
−
1
ω
s
≤
∣
ω
∣
≤
2
j
ω
s
2^{j-1}\omega_{s}\leq\left|\omega\right|\leq2^{j}\omega_{s}
2j−1ωs≤∣ω∣≤2jωs对应的
S
(
t
)
S(t)
S(t)的信息。
2.1 相对小波能量
在介绍小波熵之前,我们需要先引入一个定义:相对小波能量(Relative Wavelet Energy):
由信号与系统的知识可以知道,对于每个分辨率
j
j
j的残差信号的能量
E
j
E_j
Ej为:
E
j
=
∥
r
j
∥
2
=
∑
k
∣
C
j
(
k
)
∣
2
E_{j}=\|r_{j}\|^{2}=\sum\limits_{k}|C_{j}(k)|^{2}
Ej=∥rj∥2=k∑∣Cj(k)∣2;
并且每个采样时刻
k
k
k的能量为:
E
(
k
)
=
∑
j
=
−
N
−
1
∣
C
j
(
k
)
∣
2
E(k)=\sum\limits_{j=-N}^{-1}\lvert C_j(k)\rvert^2
E(k)=j=−N∑−1∣Cj(k)∣2;
因此,总能量可以表示为:
E
t
o
t
=
∥
S
∥
2
=
∑
j
<
0
∑
k
∣
C
j
(
k
)
∣
2
=
∑
j
<
0
E
j
E_{\mathrm{tot}}=\begin{Vmatrix}S\end{Vmatrix}^2=\sum\limits_{j<0}\sum\limits_{k}\bigl|C_{j}(k)\bigr|^2=\sum\limits_{j<0}E_{j}
Etot=
S
2=j<0∑k∑
Cj(k)
2=j<0∑Ej;
表示相对小波能量的归一化值
p
j
=
E
j
E
t
o
t
p_j=\frac{E_j}{E_{\mathrm{tot}}}
pj=EtotEj,其中
j
=
−
1
,
−
2
,
.
.
.
,
−
N
j=-1, -2, ..., -N
j=−1,−2,...,−N,表示在某一特定分辨率
j
j
j上能量占总能量的比例,显然
∑
j
p
j
=
1
\sum_{j}p_j=1
∑jpj=1,我们可以理解它为能量在尺度层面的概率分布,也就是说这个分布反映了信号在各个频带(分辨率,或者说尺度)的能量的分布情况。
2.2 小波熵
2.2.1 总小波熵
定义总小波熵(Total WE)为:
S
W
T
≡
S
W
T
(
p
)
=
−
∑
j
<
0
p
j
⋅
log
(
p
j
)
S_{\mathrm{WT}}\equiv S_{\mathrm{WT}}(p)=-\sum_{j<0}p_j\cdot\log(p_j)
SWT≡SWT(p)=−j<0∑pj⋅log(pj)
像信息熵一样,小波熵可以量化信号的有序(或无序)程度,因此它可以提供与信号相关的潜在动态过程的有用信息。
举个例子,假如被分析信号为一段单频率的周期信号,则它在小波分析中将得到很高的分辨率,即除了表征该信号频率的分辨率,相对小波能量基本为0,对于该分辨率,相对小波能量基本为1,因此总小波熵接近于0;与之相反,若被分析信号为一段全频带下的均匀噪声,则可以推断出,各分辨率的相对小波能量的贡献基本相同,此时总小波熵为最大值。
2.2.2 相对小波熵
在Cover的书中,对相对熵是如下定义的:
相对熵是两个随机分布之间距离的度量。相对熵 D ( p ∥ q ) D(p\|q) D(p∥q)度量当真实分布为 p p p而假定分布为 q q q时的无效性。
两个概率密度函数为 p ( x ) p(x) p(x)和 q ( x ) q(x) q(x)之间的相对熵或Kullback-Leibler距离定义为
D ( p ∥ q ) = ∑ x ∈ X p ( x ) log p ( x ) q ( x ) = E p log p ( x ) q ( x ) D(p\|q)=\sum\limits_{x\in\mathcal{X}}{p(x)\log\frac{p(x)}{q(x)}}=E_p\log\frac{p(x)}{q(x)} D(p∥q)=x∈X∑p(x)logq(x)p(x)=Eplogq(x)p(x)
与之相似地,定义相对小波熵(RWE):
假设有2个不同的概率分布
{
p
j
}
\{p_j\}
{pj}和
{
q
j
}
\{q_j\}
{qj},且
∑
j
p
j
=
1
\sum_{j}p_j=1
∑jpj=1,
∑
j
q
j
=
1
\sum_{j}q_j=1
∑jqj=1,它们可以表示为用尺度表示一个信号的两段或两个不同信号的小波能量的概率分布。
S
W
T
(
p
∣
q
)
=
∑
j
<
0
p
j
⋅
log
[
p
j
q
j
]
S_{\mathrm{WT}}(p|q)=\sum_{j<0}p_{j}\cdot\log\biggl[\frac{p_{j}}{q_{j}}\biggr]
SWT(p∣q)=j<0∑pj⋅log[qjpj]
相对小波熵可以用于度量两个概率分布之间相似程度,当且仅当
p
j
=
q
j
p_j=q_j
pj=qj时,RWE为零。
3 简单理解
就如同Haar最早提出小波变换时,解决了傅里叶变换无法直观显示不同频率的分量在时域中的分布情况的问题,小波熵则解决了谱熵(spectral entropy)由于受到傅里叶变换的局限,不适用非平稳信号的问题。
小波熵继承了小波变换的优点,对时域和频域均很敏感,在处理一些特殊信号时,比傅里叶变换要有用得多。
举个实例,对高斯脉冲信号和高斯白噪声信号分别求它们的小波熵:
# Wavelet Entropy
import numpy as np
import pywt
from scipy.stats import entropy
# 设置参数
N = 100 # 总点数
sigma = 5
# 生成信号
t = np.arange(N)
A_signal = np.exp(-((t - N // 2) ** 2) / (2 * sigma ** 2)) # 信号A为高斯脉冲信号
B_signal = np.random.randn(N) # 信号B为白噪声信号
# 小波变换
wavelet = 'db4'
levels = 3
A_wavelet = pywt.wavedec(A_signal, wavelet, level=levels)
B_wavelet = pywt.wavedec(B_signal, wavelet, level=levels)
# 计算每个尺度的能量
def calculate_energy(wave):
energy = [np.sum(np.abs(w) ** 2) for w in wave]
return energy
A_energy = calculate_energy(A_wavelet)
B_energy = calculate_energy(B_wavelet)
# 总能量
A_total_energy = np.sum(A_energy)
B_total_energy = np.sum(B_energy)
# 相对小波能量
B_RWE = [e / B_total_energy for e in B_energy]
A_RWE = [e / A_total_energy for e in A_energy]
# A_RWE = [0.8,0.1,0.1]
# B_RWE = [0.1,0.2,0.7]
# 小波熵
A_WE = entropy(A_RWE, base=2)
B_WE = entropy(B_RWE, base=2)
# 输出结果
print(f"信号A的相对小波能量:{A_RWE},小波熵:{A_WE:.4f}")
print(f"信号B的相对小波能量:{B_RWE},小波熵:{B_WE:.4f}")
# 确保相对小波能量之和为1
assert np.isclose(np.sum(A_RWE), 1.0)
assert np.isclose(np.sum(B_RWE), 1.0)
运行结果如图所示:
由此可见,高斯脉冲信号的小波熵明显低于白噪声信号,且相对小波能量更集中。
其他理解暂未想出,更新随缘。