【脑电信号处理与特征提取】P6-张治国：频谱分析和时频分析

频谱分析和时频分析

背景

脑电（尤其是静息态脑电）一般在频域进行分析，以刻画脑电信号的周期性特征，需要使用频谱分析来描述脑电信号功率沿频率的分布特征。
任务态脑电实验中，任务可增强或减弱脑电在特定频段的节律幅度。事件相关的频谱变化被称为事件相关同步话/去同步化（ERS/ERD），通常表示为在时间-频率域中随时间变化的频谱功率，可以通过时频分析方法进行估计。

频谱估计基本概念

时间序列信号：例如在某通道连续记录的脑电信号，可以在时间域中表征为信号幅度（或其它量值）相对于时间的变化，也可以在频率域中表征为信号功率（或其它量值）沿频率变化的分布。
频率：描述震荡波形在单位时间内周期活动的基本参数。单位事赫兹（Hz），即每秒一个周期。
频谱：时序信号的功率、幅度或相位等在频率域沿着频率等分布曲线。
频谱估计：将时域信号变换到频域频谱的估计方法，目的是观察对应周期的频率峰值，以检测信号周期性。

下面是一个模拟信号的时域波形和频域变换的例子：
在这里插入图片描述
下面是一个真实脑电信号的时域波形和频域变换：

傅立叶变换

频谱估计的基础是傅立叶变换
简单的说，傅立叶变化从“波形”的角度看世界，将时间序列信号表示为一系列波形（余弦函数）的总和。

比如下面的红色波形，可以看做六个蓝色波形的组合。
在这里插入图片描述

这六个蓝色波形都有自己的频率和幅值，从侧面看，便可以得到信号在频率域上的变化。

因此，时间序列信号就可以表示为一系列余弦函数的总和。

频谱分析：细节与注意事项

在这里插入图片描述

频谱估计方法：周期图

脑电分析中常用的几种方法：
- 周期图
- Welch法
- 多窗口法
- 基于自回归模型的估计
周期图法是一个最简单且最常用的频谱估计方法

其中x[n]是时序信号，N是信号点数，Fs是采样频率
周期图可以简单地计算为信号FFT幅值的平方，再乘以一个常数。
周期图可以显示在两种尺度上
- 线性尺度可突出显示主要的频谱峰值（如下图10Hz附近的峰值），但是会使其他频谱分量难以辨识。
- 对数尺度可以使不同频段上的分量更具视觉可比性，但是峰值较难突出。对数尺度显示脑电频谱更常见。

在这里插入图片描述

去趋势（移除超低频）

注意：脑电频谱功率随着频率增加而迅速减小，所以脑电的超低频成分（<1Hz）经常主导整个频谱
因此，在脑电频谱估计之前须执行一个去趋势，以移除超低频（通常是一个线性趋势）。
如下图所示，上方是去除超低频之前的频谱图，可以看到信号有一个逐渐上升的趋势，超低频成分非常强的幅值，将其他成分都淹没了，去除超低频之后，这个逐渐上升的趋势也消失了。

频谱估计方法：Welch法

周期法有一个缺点是方差过大，使用Welch法可以解决这个问题，步骤如下：

将信号的N个采样分成k个数据段（可重叠），每个数据段包含M个采样，相邻两个数据段有D个采样点重叠；
将数据段加窗（即将数据段点乘一个对称的钟形窗口）；
计算每个加窗数据段的FFT并获得其周期图；
所有加窗数据段周期图的平均值作为该信号最终频谱估计。

下面是一段睡眠脑电的Welch频谱估计，将一段睡眠脑电划分为多个小段，每一小段进行加窗后分别计算FFT获得周期图，再将所有小段周期图求平均值得到该信号最终频谱估计。

下面是一个闭眼脑电信号的分别用周期图和Welch法得到的功率谱对比，以及使用Welch法设置不同分M和D得到的功率谱对比。

由上面的左图可以看到Welch法得到的功率谱波形更平滑，峰值也更突出，而周期法有很大的方差。

频谱特征提取

估计频谱后，仍需从频谱中提取少量的、有意义的特征
最常见的脑电频谱特征是特定频段内的频谱功率。
- 频段范围可以遵循常规定义（即delta、theta、alpha、beta、gamma等，或更细分如alpha-1和alpha-2）
- 频段范围也可由统计分析检测出最重要或具有最强区分度的频率点（例如在该频率点两组被试或者两个实验条件下的脑电功率统计意义上显著不同）。
某一感兴趣的频段内的带限功率可以计算为该频段内所有频率点上功率的平均值或总和。
为消除个体之间的差异，经常使用相对功率（即某频段功率与整个频率范围内总功率之间的比率）作为特征。
另一类重要的脑电频谱特征是频谱峰值的频率、幅度和带宽。
- 例如，在基于稳态视觉诱发电位的脑机接口的应用中，需要定位频谱峰的频率和幅度以确定被试的行为意图
此外，还可将脑电信号频谱看作一个随机过程，从中计算其统计值或每个节律的统计值，如平均值、方差和熵等，作为脑电频谱特征。

那上面进行Welch频谱估计等睡眠脑电举例进行特征提取，可以将得到的频谱图划分为多个频段，再将每一段频谱计算一个平均值。

时频分析

频谱分析的假设：脑电信号的频谱是固定的，不随时间而改变。

显然，这种假设过分简化了脑电图的非平稳和动态特性。事实上，脑电是信号是高度非稳态的，其频谱随时间变化。

比如下面的两个信号，蓝色波形在0～0.5秒频率较低，在0.5秒之后频率较高，而绿色相反，但是二者的频谱图是一样的。所以频谱分析没有考虑时间信息。
在这里插入图片描述
另外一个例子，在两种状态下采集脑电信号（0～2s：睁眼；2～4s：闭眼）。脑电在Oz电极记录，采样率为1000Hz。

下面是频谱分析，整个脑电信号（0～4s）的频谱（红色虚线）以及脑电信号在睁眼（0～2s）和闭眼（2～4s）两个状态下的频谱（黑色实线）。频谱估计方法是Welch法，可看出，两个状态下的频谱估计大不相同，尤其是在alpha波段。
在这里插入图片描述
下图是使用短时傅立叶变换STFT估计的脑电信号时频分布，在2s左右可以清楚地观察到alpha节律功率到升高，这表明闭眼状态的开始时间。

时频分析技术

为了揭示非稳态信号的时变频谱，需要时频分析（time- frequency analysis，TFA）技术
时频分析涵盖一系列联合研究信号时域特征的方法。一般而言，有两种类型的时频分析技术：
- 时频功率分析：估计意味着估计信号在时频域上每一个特定时间-频率点的功率量，包括短时傅立叶变换（short-time Fourier transform，STFT）和连续小波变换（continuous wavelet transform，CWT）
- 时频信号分解：将信号分解为一组具有特定时频特性的叠加成分，其中典型的方法是离散小波变换（discrete wavelet transform，DWT）

短时傅立叶变换

短时傅立叶变换（STFT）是一种简单实用的时域分析法，它建立在一个通用的滑动窗口分析的思路上。
滑动窗口法假设非平稳信号可以被分成一系列短数据段，每个短数据的频谱是固定的。因此，可以在每一个短数据段上执行常规的频谱估计方法。最后，所有短数据段的频谱估计值堆叠在一起，形成在联合时频域上的一个频率功率分布图。

在这里插入图片描述
下图是一个STFT的例子：

时频分辨率与窗口选择

时频分析方法最重要的特性是它们的时间分辨率和频率分辨率（分别指的是时频分布在区分时域或频域中紧密相邻的两个信号分量的能力）。
在STFT（以及其它滑动窗口方法）中，时间分辨率和频率分辨率之间存在折中。如果选择大的分析窗口（即数据段较长），则频率分辨率较高，但时间分辨率较低；如果选择小的分析窗口（即数据段较短），则时间分辨率较高，但频率分辨率较低。

连续小波变换

与STFT不同的是，连续小波变换（CWT）的窗口大在时频域内自适应的变化，窗口选择策略是：在低频范围使用长窗口，在高频范围使用短窗口。
在这里插入图片描述

其它常用的时频分析方法

在这里插入图片描述

事件相关同步化/去同步化

什么是事件相关同步化/去同步化
ERS（Event-Related Synchronization）和ERD（Event-Related Desynchronization）是指在特定事件发生后，脑电信号在特定频率带内出现的同步或去同步现象。是对ERPs进行进一步分析的方法，用于研究不同时间段或频带中的脑电活动的同步和去同步状态。
ERS（Event-Related Synchronization）是指在特定事件发生后，脑电信号在特定频率带内出现同步增强的情况。这意味着在该时间段内，大脑的神经元活动呈现出更加同步的状态，导致脑电信号在特定频带内的幅值增加。
ERD（Event-Related Desynchronization）是指在特定事件发生后，脑电信号在特定频率带内出现去同步的情况。这意味着在该时间段内，大脑的神经元活动呈现出更加去同步的状态，导致脑电信号在特定频带内的幅值减小。