浅谈一下小波变换和EEG信号处理

小波变换(Wavelet Transform)

谈及小波变换，说不开的就是时域关系、频域关系和时频域关系(一些Papers上可能会提及Wavelet Domain，不过从我的理解来看，都是F-T-Power图谱，都大相径庭了)了。

首先简单说一说时域信号(Time Domain)

顾名思义，时域信号就是以时间作为横轴，记录信号的幅值的一种记录方法，其实许多信号在被记录时，就已经是时域信号了。这里以EEG(脑电)信号做简单描述。
EEG 原始信号就是一种时域上的数据。波形幅值随着时间变换而变换，这种分析方法具有较高的时间精度。

这段波形是我读取Bonn大学公开脑电数据集里所形成的图像，读取的是A_Z中的Z001.txt数据，共有4097个sample点。代码比较简单，我也贴出来，方便大家对照：

import numpy as np
import matplotlib.pyplot as plt

#读取信号
eeg_signal = np.loadtxt("./A_Z/Z001.txt")

#绘制图像
plt.figure(figsize=(10,5))
plt.plot(eeg_signal)
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.title('EEG signal')

plt.show()

在计算机中，数据一定会是以离散的形式存储起来的，根据香农采样定理，当
$$f\ge 2f_{max}$$
时，原始的连续数据转存为离散(进行A-D转换)后，再转换回来就不会出现失真的情况。
而对EEG数据进行时域分析时，就会离不开时间相关电位(ERP)方法，以下是一个简单的ERP处理流程，

1. 滤波，滤去噪声干扰；
2. 分段 & 加标签；
3. 去基线；
4. 叠加平均；
   关于ERP，读者可以参考其它博主的博客，题主在这不作过多赘叙。比如这篇事件相关电位（ERP）的简单处理流程（含MATLAB代码）。

频域分析

频域分析使用的是傅里叶变换(Fourier Transform)。根据傅里叶定理，任何连续测量的时序或信号，都可以表示为不同频率的正弦波。
信号的无限叠加。脑电信号可以看成就是由不同的正弦信号混合而成的混合信号，通过傅里叶变换，就能够将这个混合信号重新分解成具有不同频率的正弦波，从而获得频域上的信息。
频域分析不仅可以用于分析任务态的数据，还常用于分析静息态的数据。
而关于傅里叶级数和傅里叶变换，其实推导并不难，这里我只简单叙述一下，有兴趣详细推导的，可以看一下网上的推导过程，都很详细。

傅里叶级数

傅里叶级数其本质上是因为三角函数族的正交性所求得的，对于一般的傅里叶级数来说
$$\begin{gathered} f(x)=\sum _{n=0}^{\infty }{a}_n\cos nx+\sum _{n=0}^{\infty }{b}_n\sin nx \\ =a_0+\sum _{n=1}^{\infty }{a}_n\cos nx+\sum _{n=1}^{\infty }{b}_n\sin nx. \end{gathered}$$
其中，
$$\begin{gathered} a_0=\frac{1}{\pi }{\int }_{-\pi }^{\pi }f(x)dx, \\ {a}_n=\frac{1}{\pi }{\int }_{-\pi }^{\pi }f(x)\cos nxdx,(n=1,2,\dots ) \\ {b}_n=\frac{1}{\pi }{\int }_{-\pi }^{\pi }f(x)\sin nxdx,(n=1,2,\dots ) \end{gathered}$$
,上述的傅里叶级数的周期仅为2Π，而在整个函数域上，不可能只是周期为2Π的，甚至有可能没有周期，所以这做了一个相关的变换。
将周期2Π拓宽到2L,得相应的级数为
$$f(x)=\frac{a_0}{2}+\sum _{n=1}^{\infty }{a}_n\cos \frac{n\pi x}{L}+\sum _{n=1}^{\infty }{b}_n\sin \frac{n\pi x}{L}.$$
其中，
$$\begin{gathered} a_0=\frac{1}{L}{\int }_{-L}^{L}f(x)dx, \\ {a}_n=\frac{1}{L}{\int }_{-L}^{L}f(x)\cos \frac{n\pi x}{L}dx,(n=1,2,\dots ) \\ {b}_n=\frac{1}{L}{\int }_{-L}^{L}f(x)\sin \frac{n\pi x}{L}dx,(n=1,2,\dots ) \end{gathered}$$
即在函数域上所有的函数都能够被级数表示出来。再由欧拉公式，
$$e^{ix}=\cos x+i\sin x$$
就能得出
$$\begin{gathered} f(x)=\sum _{-\infty }^{\infty }{C}_n{e}^{\frac{n\pi x}{L}i} \\ {C}_n=\frac{1}{2L}{\int }_{-L}^{L}f(x)\cdot e^{-i\cdot \frac{n\pi x}{L}}dx. \end{gathered}$$
再由角频率公式，即可得
$$\begin{gathered} f(x)=\sum _{-\infty }^{\infty }{C}_n{e}^{i\cdot n{\omega }_{0}x}, \\ {C}_n=\frac{1}{2L}{\int }_{-L}^{L}f(x)\cdot e^{-i\cdot n{\omega }_{0}x}dx \end{gathered}$$
而当L趋向于无穷大时，整个级数就能表示所有函数域上的函数了，即傅里叶变换公式为
$$F(\omega )=F[f(x)]={\int }_{-\infty }^{\infty }f(x)\cdot e^{-i\cdot \omega x}dx$$
对应的逆变换为
$$f(x)=F^{-1}[F(\omega )]=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }F(\omega )\cdot e^{i\cdot \omega x}d\omega$$

而傅里叶变换即是在整个频域分析上最为重要的一环。但由于傅里叶变换仅能将时间信号和频率信号相互转换，且即使用到了短时傅里叶变换将时间加窗，也不能完全将高频和低频的信息完全地呈现出来。而小波变换就存在这样的优点。

时频域分析——小波变换(Wavelet Tranform)

首先来看一张关于FT，STFT，WT之间的对比图(这张图来自台湾阳明大学卢家峰老师的Notes)，方便大家来看。

一般地讲，信号的高频特征往往需要较少的时间来表现；而对应的，低频特征则需要更多的时间来表现；一般的傅里叶变换或者说快速傅里叶变换(Fast Fourier Transform)仅仅是将时域信号和频域信号之间进行相互转换，并没有具体的去划分来表现特征。且傅里叶变换仅仅能够处理平稳信号，对于非平稳信号并不能有较好地体现。
而短时傅里叶变换(Short-Time Fourier Transform)也只是在时域上加窗(Window)，并不能很好地体现信号的特征(高频低时，低频高时)。
那么小波变换(Wavelet Transform)就相对来说解决了这么个问题。

连续小波变换(Continous Wavelet Transform)

在这里，就直接给出公式，但不作过多推导和讲述，因为，计算机能够处理的往往只是离散信号，对于一般的连续信号也需要对这些信号进行采样来获得离散信号来进行处理。

一般连续小波变换公式：
$$\begin{gathered} W_{a,b}={\int }_{-\infty }^{+\infty }f(t){\psi }_{a,b}(t)dt, \\ \psi a,b=\frac{1}{\sqrt{|a|}}\psi (\frac{t-b}{a}) \end{gathered}$$
从小波图形上来说，连续小波变换较离散小波变换所得到的波形其实更平滑，但由于连续小波变换具有如下缺点：

1. CWT 具有冗余性；
2. CWT 的小波数量无穷；
3. CWT 不具备解析解。
   所以，一般会将连续小波变换进行离散化，一般是将伸缩因子(a)和平移因子(b)进行离散化；
   即
   $$\begin{gathered} aj=2^j, \\ {b}_{j,k}=k\cdot 2^j \end{gathered}$$
   得离散小波函数为
   $${\psi }_{j,k}=\frac{1}{\sqrt{2^j}}\psi (\frac{t}{2^j}-k)$$

离散小波变换(Discrete Wavelet Transform)

从上面的公式上去理解小波变换可能会非常抽象，下面简单讲一个离散小波变换的例子。
离散小波变换其本质就是一个相关(correlation)运算,只是使用的滤波器不相同。

这是一张Level1的小波变换图，使用一个4系数的高(低)通滤波器对数据进行处理，其实就是一个步长为2的相关运算。x是一个16位的数据，而滤波器是4位的，滤波器在数据上滑动，每次与其对应的x上的4位数据做点击，就获得w中的系数，以步长为2，就可以将数据分成等长的两部分了。因为只做了一次小波变换，所以它是level1的。

EEG信号和DWT

对于小波变换来说，使用不同的小波(其实就是滤波器不同)，意味着所用的小波变换就不相同，对于EEG信号来说，小波中与之相关性最好的就是Daubechies-4小波或者db_4小波。
db_4小波的低通滤波器为
$$h=[\frac{1+\sqrt{3}}{4\sqrt{2}},\frac{3+\sqrt{3}}{4\sqrt{2}},\frac{3-\sqrt{3}}{4\sqrt{2}},\frac{1-\sqrt{3}}{4\sqrt{2}}]$$
利用小波函数(滤波器)之间的相互正交关系，就可以推出对应的高通滤波器g了。这里给出公式，但不做计算。
$$g_k=(-1{)}^k{h}_{n-k-1},k\in 0,1,\dots ,n-1$$
然后，信号再经过小波变换后会有两个系数A(approximation)和D(Detail)，A系数是经过低通滤波器获得的，而D系数是通过高通滤波器所得到的。以下是我是用db_4小波对一段EEG原始信号进行level4变换得到一个时域图(基于matlab实现，很简单调包就行)。

从小波变换的角度来理解，就是说，A3信号是A4信号和D4信号的叠加，A2信号是A3信号和D3信号的叠加，依此类推，即能得到原始信号。
可以看出，EEG信号在经过小波变换后，被分解成了高频(可能存在噪声)和低频成分，所以可说，小波变换其实是对信号进行降频降噪的处理。

而一般来说，不同研究所需的EEG信号频带是不同的，所以用小波变换做EEG信号处理，其实也需要控制小波类型和使用小波变换的次数。

对于癫痫检测来说，使用小波变换将 EEG 信号分解到不同频带可以提取不同频率范围内的信号特征，有助于分析和诊断癫痫。通常情况下，将 EEG 信号分解到 theta（4-8 Hz）、alpha（8-13 Hz）和 beta（13-30 Hz）频带可以提供有用的信息。

一般使用小波变换的方法如下：先是选择基本小波，然后将基本小波平移，缩放，生成一个函数族，构成一个框架，将信号向该框架上投影得到分解信号。选择合适的小波函数，对脑电信号变换，突出需要的特征波形。通过阈值判断，确定是否为待检测的特征波形。