EEG信号分析：关于时频分析方面的个人理解

2023/1/12-2023/1/脑机接口学习内容一览：

        这一篇博客里，主要在写博客的过程中总结和思考自己在前一段时间所进行的频域分析工作有何意义，以及探明时频分析几种主要方法的基本工作原理，最后做一下代码方面的总结。

1.频谱估计

        在脑电信号分析的频谱估计方面，以傅里叶变换为基础实现多个方法。

1.1 傅里叶变换

        在傅里叶变换方面，之前一直对该变换到底是如何实现的持有怀疑态度：很多教程都只讲解了傅里叶变换的原理，但是这只让我理解了傅里叶变换这个概念，傅里叶变换到底如何在计算机上实现？如何用傅里叶变换来处理一段信号？它是怎么将时域信号转化为频域信号的？这其中的具体实现过程让人很迷惑。但是参考文献【1】很好的解答了我这方面的疑问。

        在Steven W. Smith, Ph.D.的解释中，实数傅里叶变换可以将一段长度为N的信号最多分解为N/2+1个正余弦信号，分解出的信号可以理解为是通过下标 0 - N/2 总共 N/2+1 个 Re[X] + Im[X] * i（这个i为复数）式子这种形式进行表示，其中前一个未知数表示cos频率幅度值，后一个表示sin频率幅度值，至于为什么把正弦塞到复数那里，这只是为了让式子表达和计算更加方便，而从N/2这个数量，我们可以联想到频谱中的最大频率为 freqs / 2。

        其中频率幅度值与我之前一直疑惑的问题有关：傅里叶变换是如何从原始信号中判断出各个频率正余弦波的振幅的。从我的理解上来说，傅里叶变换就是把一个函数通过各频率的正余弦之和来等价表示。

1.1.2 DFT合成运算

        从这里可以看出每一个时间点的振幅都由全部频率的某点振幅之和来综合表示。这同样也体现了傅里叶变换的缺陷之一，所有频率都对时段上的某点信号存在影响，难以分析非周期信号。

        同样，从这个式子出发，我们可以看出如果k取0~N/2来表示Re和Im数组，则两个数组内下标为k的数代表的频率可计算为k / N hz，这也稍微解答了我之前的一些疑问。

         

        如果我们手头已经有了通过信号分解得到的各个正余弦波的信息，那么可以用这个公式来合成源信号。证明不再赘述，此公式形式与傅里叶级数相似。

        DFT中的相关转化同样参考【1】的频谱密度理解方式，我觉得未尝不可将频谱密度的离散点理解为以各点的带宽作为柱宽的柱状图，频谱密度其实是将每一个带宽范围内的频率幅度相加近似地以一个点的纵坐标来表示，因此DFT的合成公式可以看作是一种逼近的合成，这种合成的精细程度很明显取决于带宽的长度，带宽越小越精确。

      

 1.1.2 DFT分解运算

        1 通过联立方程进行求解, 但这是这种方法计算量非常的大且极其复杂，很少被采用

        2 利用信号的相关性（correlation）进行计算    

        3 快速傅立叶变换（FFT）

        相关性的计算其实比较好理解，就是把一个待检测信号波乘以另一个信号波，得到一个新的信号波，再把这个新的信号波所有的点进行相加，从相加的结果就可以判断出这两个信号的相似程度。可以通过把输入信号和每一种频率的正余弦信号进行相乘（关联操作），从而得到原始信号与每种频率的关联程度（即总和大小），这个结果便是我们所要的傅立叶变换结果，下面两个等式便是我们所要的计算方法：

        为了便于理解，这里写一个特殊的例子：根据正交的概念，如果某次计算得到的相加和为0，则可以说明这个频率对原始信号没有贡献。

        相关性计算虽然比较容易理解，但是根据我这些天写下的代码，相对于相关性计算，还是FFT即快速傅里叶变换使用的频率较高。

1.1.3快速傅里叶变换（FFT）

        快速傅里叶变换是离散傅里叶变换（DFT）的快速算法，这部分只看【1】不太懂，所以参考了文献【2】。

        按照参考资料的讲解FFT是一种特殊的复数形式DFT，因此上一小节似乎不应该把它和联立方程和相关性计算方法列在一起，容易造成误解，在【2】中使用的是复杂度简化为nlgn后的分治快速算法，与【1】中的相关性算法不同，但是更加便于理解。

      

        具体不再赘述，注意事项可查看【3】。

        另外，在选取fft点数方面，N越大频率的分辨率越好，末尾补零能使数据更加平滑。

1.2 周期图法

        周期图主要通过FFT变换计算得到功率值来展现功率谱。加窗周期图只取有限的一段数据进行傅里叶变换，加窗这个操作十分重要，能缓解傅里叶变换难以解决的跳变问题，具体参考【4】。平均周期图将总长为N的数据分为M个长为L段的数据段分别运算取平均，使方差显著减小，但是分辨率同样减小。

1.3 welch法和multitaper（多窗口）法

        这两个方法在我们使用mne库的时候能经常看见，是函数中某一项参数的两个可选择项，应用范围较广，通过平均加窗数据来降低频率估计的方差。在使用的过程中，会自动计算我们所需要的最佳窗口长度，具体可参考【4】与【5】。

1.4 自回归ar模型

        这种方法的原理看起来与神经网络比较相似，使用参数法通过特定模型来描述信号，估计模型参数，最后计算频谱。但是这个方法在mne中似乎并没有对应函数，非参数法引用比它更加广泛。

---

2.时频分析

        频谱分析是建立在假设脑电信号是频谱固定且不随时间变化的稳定信号，然而我们清楚脑电信号不可能是稳定的，且频谱随时间变化。

2.1 短时傅里叶变换（STFT）

        STFT建立在滑动窗口分析的基础上，即将数据分为多个短数据段，运用滑动窗口进行傅里叶变换来分析，从这方面看来前面的加窗估计似乎与之相似，其区别是在时域方面操作。这个方法综合前面的加窗法来看很容易理解。

2.2 小波变换

        小波变换的自适应窗口变换使其更适合对于脑电信号的分析，具体mne函数可查看往期博客【6】。

学习总结：

        在这一模块的学习中，对傅里叶变换等时频分析方式具有了更深层次的原理理解，但是仍然需要在代码方面有所加强，接下来应该聚焦于机器学习结合脑电分析方面的学习，使自己的水平有更大的提升。

参考文献：

【1】理解傅里叶变换

【2】十分简明易懂的FFT（快速傅里叶变换）

【3】脑电的频谱分析和时频分析-EEG Processing and Feature 5

【4】理解FFT, STFT, 加窗的含义

【5】脑电分析mne库函数compute_psd()记录

【6】mne库脑电时频信号分析函数解读：小波变换及方法比较