基于时频变换的脑波信号（EEG）处理方法

离散傅里叶变换（DFT）

在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作其周期延拓的变换。在实际应用中通常采用快速傅里叶变换计算DFT 。
对于N点序列 x[n]0≤n≤Nx[n]_{0 \leq n \leq N}x[n]0≤n≤N, 它的离散傅立叶变换为 (DFT)(\mathrm{DFT})(DFT) 为:

X[k]=∑n=0N−1e−j2πNnkx[n]X[k]=\sum_{n=0}^{N-1} e^{-j \frac{2 \pi}{N} n k} x[n] X[k]=n=0∑N−1e−jN2πnkx[n]

其中 k=0,1,….,N−1\mathrm{k}=0,1, \ldots ., \mathrm{N}-1k=0,1,….,N−1, 上面的式子展开一下：

X[k]=∑n=0N−1x[n]⋅[cos⁡(2πNkn)−jsin⁡(2πNkn)]X[k]=\sum_{n=0}^{N-1} x[n] \cdot\left[\cos \left(\frac{2 \pi}{N} k n\right)-j \sin \left(\frac{2 \pi}{N} k n\right)\right] X[k]=n=0∑N−1x[n]⋅[cos(N2πkn)−jsin(N2πkn)]

DFT的时间复杂度位O(N2)O(N^2)O(N2)

快速傅里叶变换(FFT)

假设多项式f(x)=a0+a1∗x+⋯+an−2∗xn−2+an−1∗xn−1f(x)=a_0+a_1*x+\dots+a_{n-2}*x^{n-2}+a_{n-1}*x^{n-1}f(x)=a0+a1∗x+⋯+an−2∗xn−2+an−1∗xn−1，那么我们可以按aaa下标的奇偶性将f(x)f(x)f(x)分成两部分

f(x)=(a0+a2x2+a4x4+⋯+an−2xn−2)+(a1x+a3x3+a5x5+⋯+an−1xn−1)\begin{aligned} f(x)&=(a_0+a_2x^2+a_4x^4+\dots+a_{n-2}x^{n-2})\\ &+(a_1x+a_3x^3+a_5x^5+\dots+a_{n-1}x^{n-1}) \end{aligned} f(x)=(a0+a2x2+a4x4+⋯+an−2xn−2)+(a1x+a3x3+a5x5+⋯+an−1xn−1)

令

f1(x)=(a0+a2x1+a4x2+⋯+an−2xn2−1)f2(x)=x(a1+a3x+a5x2+⋯+an−1xn2−1)\begin{aligned} f_1(x)&=(a_0+a_2x^1+a_4x^2+\dots+a_{n-2}x^{\dfrac{n}{2}-1})\\ f_2(x)&=x(a_1+a_3x+a_5x^2+\dots+a_{n-1}x^{\dfrac{n}{2}-1}) \end{aligned} f1(x)f2(x)=(a0+a2x1+a4x2+⋯+an−2x2n−1)=x(a1+a3x+a5x2+⋯+an−1x2n−1)

则有f(x)=f1(x2)+xf2(x2)f(x)=f_1(x^2)+xf_2(x^2)f(x)=f1(x2)+xf2(x2)

带入 ωnk(k<n2)\omega_{n}^{k}\left(k<\frac{n}{2}\right)ωnk(k<2n) 可得

f(ωnk)=f1(ωn2k)+ωnkf2(ωn2k)=f1(ωn2k)+ωnkf2(ωn2k)\begin{aligned} f\left(\omega_{n}^{k}\right) &=f_{1}\left(\omega_{n}^{2 k}\right)+\omega_{n}^{k} f_{2}\left(\omega_{n}^{2 k}\right) \\ &=f_{1}\left(\omega_{\frac{n}{2}}^{k}\right)+\omega_{n}^{k} f_{2}\left(\omega_{\frac{n}{2}}^{k}\right) \end{aligned} f(ωnk)=f1(ωn2k)+ωnkf2(ωn2k)=f1(ω2nk)+ωnkf2(ω2nk)

带入 ωnk+n2(k<n2)\omega_{n}^{k+\frac{n}{2}}\left(k<\frac{n}{2}\right)ωnk+2n(k<2n) 可得

f(ωnk+n2)=f1(ωn2k+n)+ωnk+n2f2(ωn2k+n)=f1(ωn2k∗ωnn)−ωnkf2(ωn2k∗ωnn)=f1(ωn2k)−ωnkf2(ωn2k)\begin{aligned} f\left(\omega_{n}^{k+\frac{n}{2}}\right) &=f_{1}\left(\omega_{n}^{2 k+n}\right)+\omega_{n}^{k+\frac{n}{2}} f_{2}\left(\omega_{n}^{2 k+n}\right) \\ &=f_{1}\left(\omega_{n}^{2 k} * \omega_{n}^{n}\right)-\omega_{n}^{k} f_{2}\left(\omega_{n}^{2 k} * \omega_{n}^{n}\right) \\ &=f_{1}\left(\omega_{n}^{2 k}\right)-\omega_{n}^{k} f_{2}\left(\omega_{n}^{2 k}\right) \end{aligned} f(ωnk+2n)=f1(ωn2k+n)+ωnk+2nf2(ωn2k+n)=f1(ωn2k∗ωnn)−ωnkf2(ωn2k∗ωnn)=f1(ωn2k)−ωnkf2(ωn2k)

上面的两个式子只有一个常数项不同，当求出第一个式子的值后可以以 O(1)\mathrm{O}(1)O(1) 的时间复杂度求出第二个式子的值。所以原来的问题缩小到了之前的一半，FFT的时间复杂度即为 O(nlog⁡n)\mathrm{O}(n \log n)O(nlogn) 。

对于以上的变换，存在一个缺陷，即仅分析了频率成分，忽略了每个频率出现的时间，这就表示了它仅适用于平稳信号，而无法应用于非平稳信号。

示例

对于图上三种不同的信号，得到的频谱图是几乎相同的，这表示傅里叶变换忽略了信号的时间信息。为了分析非平稳信号，短时傅里叶变换诞生了。

短时傅里叶变换（STFT）

离散样本的短时傅里叶变换的定义：

X(n,w)=∣X(n,w)∣2X(n, w)=|X(n, w)|^2 X(n,w)=∣X(n,w)∣2

式中， X(m)X(m)X(m)是输入信号；w(m)w(m)w(m)是窗函数；它在时间上翻转并且有n个样本的偏移量； X(n,w)X(n, w)X(n,w)是定义在样本（时间）和频率上的二维函数；输入信号的时频图（Spectrogram)定义为S(n,w)=∣X(n,w)∣2S(n, w)=|X(n, w)|^2S(n,w)=∣X(n,w)∣2

窗函数

汉宁窗

wHann (n−m)=12[1+cos⁡(πn−mm)]w_{\text {Hann }}(n-m)=\frac{1}{2}\left[1+\cos \left(\pi \frac{n-m}{m}\right)\right] wHann (n−m)=21[1+cos(πmn−m)]

海明窗

wHamming (n−m)=12[1.08+0.92cos⁡(πn−mm)]w_{\text {Hamming }}(n-m)=\frac{1}{2}\left[1.08+0.92 \cos \left(\pi \frac{n-m}{m}\right)\right] wHamming (n−m)=21[1.08+0.92cos(πmn−m)]

高斯窗

wGauss (n−m)=e−12(n−mσ)2w_{\text {Gauss }}(n-m)=e^{-\frac{1}{2}\left(\frac{n-m}{\sigma}\right)^{2}} wGauss (n−m)=e−21(σn−m)2

其中汉宁窗和汉明窗定义 for ∣n−m∣<m|n-m|<m∣n−m∣<m ，高斯窗 for ∣n−m∣<3σ|n-m|<3 \sigma∣n−m∣<3σ

对于帧长固定的短时傅里叶变换，在全局范围内的时间分辨率和频率分辨率是固定的。如果我们想要在低频区域具有高频率分辨率，在高频区域具有高时间分辨率，显然STFT是不能满足要求的。我们继续引入另一种时频分析方法——小波变换。

小波变换（WT）

对于任意能量有限信号f(t)f(t)f(t)，其连续小波变换（CWT)定义为

Wf(a,b)=1af−∞+∞ψ∗(t−ba)dtW_f(a, b)=\dfrac{1}{\sqrt a}f_{-\infty}^{+\infty}\psi*(\dfrac{t-b}{a})dt Wf(a,b)=a1f−∞+∞ψ∗(at−b)dt

其中ψ(t)\psi(t)ψ(t)是母小波或者基本小波，满足ψ(±∞)=0,ψ(0)=0,f−∞+∞ψ(t)dt=0\psi(\pm \infty)=0, \psi(0)=0, f_{-\infty}^{+\infty}\psi(t)dt=0ψ(±∞)=0,ψ(0)=0,f−∞+∞ψ(t)dt=0

小波变换将无限长的三角函数基转换成有限长会衰减的小波基，解决了傅里叶变换的吉布斯效应。

吉布斯效应

对于变换突然且剧烈的信号，即使只有一小段变换，傅里叶也不得不用大量的三角波去拟合。

而对于衰减的小波，当突变时只有小波函数和信号突变处重叠时，系数不为0。

信号的频率特性

频谱

频谱是频率的分布曲线，复杂振荡分解为振幅不同和频率不同的谐振荡，这些谐振荡按照频率排列的图形就称为频谱，它描述了信号在各个频率上的分布大小。信号的频谱通常可以由傅里叶变换来得到。

功率谱密度（PSD）

PSD——Power Spectral Density 是表征信号的功率能量与频率的关系的物理量，是指用密度的概念表示信号功率在各频率点的分布情况，是对随机变量均方值的量度，是单位频率的平均功率量纲；也就是说，对功率谱在频域上积分就可以得到信号的平均功率，而不是能量。PSD经常用来研究随机振动信号或根据频率分辨率做归一化。

功率谱密度的单位通常用每赫兹的瓦特数（W/Hz）表示，另一种单位 dB，当单位为dB时是因为对数据做了对数处理(10logX)，为了拉高低振幅成分，便于观察噪声中的周期信号

能量谱密度

单位频率的幅值平方和量纲，能量谱密度曲线下面的面积才是这个信号的总能量。

能量谱是功率谱密度函数在相位上的卷积，也是幅值谱密度函数的平方在频率上的积分；功率谱是信号自相关函数的傅里叶变换，能量谱是信号本身傅立叶变换幅度的平方。

熵理论

近似熵

**近似熵（ApEn）**是一种用于量化时间序列波动的规律性和不可预测性的非线性动力学参数，它用一个非负数来表示一个时间序列的复杂性，反映了时间序列中新信息发生的可能性，越复杂的时间序列对应的近似熵越大.

算法表述如下:

设存在一个以等时间间隔采样获得的 NNN 维的时间序列 u(1),u(2),…,u(N)u(1), u(2), \ldots, u(N)u(1),u(2),…,u(N).
定义算法相关参数 m,rm, rm,r ，其中, mmm 为整数，表示比较向量的长度, rrr 为实数, 表示"相似度的度量值.
重构 mmm 维向量 X(1),X(2),…,X(N−m+1)X(1), X(2), \ldots, X(N-m+1)X(1),X(2),…,X(N−m+1),其中 X(i)=[u(i),u(i+1),…,u(i+m−1)]X(i)=[u(i), u(i+1), \ldots, u(i+m-1)]X(i)=[u(i),u(i+1),…,u(i+m−1)].
对于 1≤i≤N−m+11 \leq i \leq N-m+11≤i≤N−m+1 ，统计满足以下条件的向量个数
Cim(r)=(C_{i}^{m}(r)=(Cim(r)=( number of X(j)X(j)X(j) such that d[X(i),X(j)]≤r)/(N−m+1)d[X(i), X(j)] \leq r) /(N-m+1)d[X(i),X(j)]≤r)/(N−m+1)
其中, d[X,X∗]d\left[X, X^{*}\right]d[X,X∗] 定义为 d[X,X∗]=max⁡a∣u(a)−u∗(a)∣d\left[X, X^{*}\right]=\max _{a}\left|u(a)-u^{*}(a)\right|d[X,X∗]=maxa∣u(a)−u∗(a)∣
u(a)u(a)u(a) 为向量 XXX 的元素， ddd 表示向亘 X(i)X(i)X(i) 与 X(j)X(j)X(j) 的距离，由对应元素的最大差值决定， jjj 的取值范围为 [1,N−m+1][1, N-m+1][1,N−m+1], 包括 j=j=j= iii.
定义 Φm(r)=(N−m+1)−1∑i=1N−m+1log⁡(Cim(r))\Phi^{m}(r)=(N-m+1)^{-1} \sum_{i=1}^{N-m+1} \log \left(C_{i}^{m}(r)\right)Φm(r)=(N−m+1)−1∑i=1N−m+1log(Cim(r))
则近似樀 (ApEn)(\mathrm{ApEn})(ApEn) 定义为

ApEn=Φm(r)−Φm+1(r)A p E n=\Phi^{m}(r)-\Phi^{m+1}(r) ApEn=Φm(r)−Φm+1(r)

log⁡\loglog 表示自然对数， mmm 和 rrr 由第2步定义 ...
参数选择：通常选择参数 m=2m=2m=2 或 m=3;rm=3 ; rm=3;r 的选择在很大程度上取决于实际应用场景，通常选择 r=0.2∗stdr=0.2 * s t dr=0.2∗std, 其中std表示原时间序列的标准差.
根据相关文献 [1][1][1], 在实际应用中选择 d[X(i),X(j)]<rd[X(i), X(j)]<rd[X(i),X(j)]<r 或 d[X(i),X(j)]≤rd[X(i), X(j)] \leq rd[X(i),X(j)]≤r 都是可行的. 如果一个时间序列的规律性比较强，则其近似樀值(ApEn)比较小，对应地，一个比较复杂的时间序列则对应一个较大的熵值.

近似熵的优势:
a. 对数据长度的依赖性较小. ApEn可以用于小数据样本(n < 50 n<50n<50)，并可实现实时计算.
b. 抗噪声能力较强. 如果数据含有噪声，则可以将ApEn与噪声水平进行比较，以确定原始数据中真实信息的表达程度.

样本熵

样本熵(SampEn)是基于近似熵(ApEn)的一种用于度量时间序列复杂性的改进方法，在评估生理时间序列的复杂性和诊断病理状态等方面均有应用.由于样本熵是近似熵的一种改进方法，因此可以将其与近似熵联系起来理解.

设存在一个以等时间间隔采样获得的 NNN 维的时间序列 u(1),u(2),…,u(N)u(1), u(2), \ldots, u(N)u(1),u(2),…,u(N).
定义算法相关参数 m,rm, rm,r, 其中， mmm 为整数，表示比较向量的长度， rrr 为实数，表示"相似度”的度量值.
重构 mmm 维向量 X(1),X(2),…,X(N−m+1)X(1), X(2), \ldots, X(N-m+1)X(1),X(2),…,X(N−m+1),其中 X(i)=[u(i),u(i+1),…,u(i+m−1)]X(i)=[u(i), u(i+1), \ldots, u(i+m-1)]X(i)=[u(i),u(i+1),…,u(i+m−1)].
对于 1≤i≤N−m+11 \leq i \leq N-m+11≤i≤N−m+1 ，统计满足以下条件的向量个数

Bim(r)=(number of X(j)such that d[X(i),X(j)]≤r)/(N−m),i≠jB_{i}^{m}(r)=(\text { number of } X(j) \text { such that } d[X(i), X(j)] \leq r) /(N-m), i \neq j Bim(r)=( number of X(j) such that d[X(i),X(j)]≤r)/(N−m),i=j

其中, d[X,X∗]d\left[X, X^{*}\right]d[X,X∗] 定义为 d[X,X∗]=max⁡a∣u(a)−u∗(a)∣,X≠X∗d\left[X, X^{*}\right]=\max _{a}\left|u(a)-u^{*}(a)\right|, X \neq X^{*}d[X,X∗]=maxa∣u(a)−u∗(a)∣,X=X∗
u(a)u(a)u(a) 为向量 XXX 的元素， ddd 表示向量 X(i)X(i)X(i) 与 X(j)X(j)X(j) 的距离，由对应元素的最大差值决定， jjj 的取值范围为 [1,N−m+1][1, N-m+1][1,N−m+1], 但是 j≠j \neqj= iii.
5. 求 Bim(r)B_{i}^{m}(r)Bim(r) 对所有 iii 值的平均值，记为 Bm(r)B^{m}(r)Bm(r), 即

Bm(r)=(N−m+1)−1∑i=1N−m+1Bim(r)B^{m}(r)=(N-m+1)^{-1} \sum_{i=1}^{N-m+1} B_{i}^{m}(r) Bm(r)=(N−m+1)−1i=1∑N−m+1Bim(r)

Aik(r)=(number of X(j)such that d[X(i),X(j)]≤r)/(N−k),i≠jA_{i}^{k}(r)=(\text { number of } X(j) \text { such that } d[X(i), X(j)] \leq r) /(N-k), i \neq j Aik(r)=( number of X(j) such that d[X(i),X(j)]≤r)/(N−k),i=j

则样本熵(SampEn)定义为

SampEn=lim⁡N→∞{−ln⁡[Ak(r)/Bm(r)]}S a m p E n=\lim _{N \rightarrow \infty}\left\{-\ln \left[A^{k}(r) / B^{m}(r)\right]\right\} SampEn=N→∞lim{−ln[Ak(r)/Bm(r)]}

由于在实际计算应用过程中， NNN 不可能为 ∞\infty∞ ，因此当 NNN 取有限值时，样本熵估计为

SampEn⁡=−ln⁡[Ak(r)/Bm(r)]}\left.\operatorname{SampEn}=-\ln \left[A^{k}(r) / B^{m}(r)\right]\right\} SampEn=−ln[Ak(r)/Bm(r)]}

其中， lnl nln 表示自然对数， mmm 和 rrr 由第2肯定义 stds t dstd 表示原时间序列的标准差.

近似熵与样本熵理论上的对比
设B为维数为m 时，时间序列的自相似概率；A 为维数为k = m + 1时，时间序列的自相似概率，得出C P = A / B。近似熵的计算是以− l n ( C P )为模型，并且计算出了所有模型的平均值。为了防止出现计算l n ( 0 ) ln(0)ln(0)的情况，近似熵在算法的第4步中包含了与自身向量的比较，这种方式与新信息观点是不相容的，所以一定会存在偏差。不同的是，样本熵计算的是和的对数，且不包含与自身向量的比较，所以其优势在于包含更大的A、B，以及更加准确的C P估计.

与近似熵相比，样本熵具有两个优势：样本熵的计算不依赖数据长度；样本熵具有更好的一致性，即参数m mm和r rr的变化对样本熵的影响程度是相同的.

后续

喜欢的话可以关注一下我的公众号技术开发小圈，尤其是对深度学习以及计算机视觉有兴趣的朋友，我会把相关的源码以及更多资料发在上面，希望可以帮助到新入门的大家！