EEG相位检测和预测方法调研
1、需求
调研一下EEG相位检测的不同算法（比如传统的希尔伯特变换）和相位预测的方法（比如用深度学习来预测），要考虑硬件上采样率的需求，采样率太低会导致相位丢失。
2、方法
2.1 脑信号瞬时参数估计MATLAB工具箱
Github链接：
工具箱介绍：
大脑信号瞬时参数估计MATLAB工具箱，该工具箱实现了不同的方法，估计大脑信号（即EEG、EMG等）的瞬时参数，如相位，信号的频率和包络，并在各种应用中计算一些相关的已知量。
其中计算相位的方法有三种：
1、通过FIR滤波和希尔伯特变换使用传统解析表示法进行瞬时相位（IP）的估计。
2、使用Zer极点扰动相位估计法（Z-PPP）通过解析表示、IIR滤波器和前后向滤波进行瞬时相位估计。
3、使用传递函数扰动相位估计法（TFP）通过解析表示、IIR滤波器和前后向滤波进行瞬时相位估计。
注：在调用三种算法时都可以设置采样频率。
2.1.1 传统的瞬时相位估计方法（希尔伯特变换）
对于信号$x(t)$，第一步是进行频率滤波，得到窄带信号$x_f(t)$。因为如果信号不具有比较窄带的频谱，提取的相位信号可能会没有意义。因此，滤波信号的分析形式如下：
其中$H\{x_f(t)\}$是滤波信号的希尔伯特变换。使用解析形式，瞬时振幅和相位的唯一定义如下：
在传统的相位估计过程中，广泛使用线性相位滤波器，即FIR滤波器，以使信号相位的失真最小。
在工具箱中实现此功能的是Phase_Ex_Trad函数（函数具体说明见User_Guide_Ver2.3.pdf 3.1节），通过m文件Test_Phase_Ex_Trad.m可以进行调用。
相关文献：
[1] M.Chavez,M.Besserve,C.Adam,J.Martinerie,“Towards a Proper Estimation of Phase Synchronization From Time Series,”Journal of Neuroscience Methods,vol.154,2006.
[2] B.Picinbono,“On the Instantaneous Amplitude and Phase of Signals, “IEEE Transaction Signal Processesing, vol.45,pp.552-560, 1997.
2.1.2 零极点扰动相位估计法
零极点扰动相位估计法（Z-PPP）是一种提取信号瞬时相位的新方法。该方法基于在零相位滤波过程中使用IIR滤波器，完全没有相位失真，以使信号窄带，然后扰动所用滤波器的频率响应，以减少该过程的副作用。在Z-PPP方法中，对2.1.1中描述的程序进行了修改，以使估计的瞬时相位信号更加可靠和有意义。
在工具箱中实现此功能的是Phase_Ex_ZPPP函数（函数具体说明见User_Guide_Ver2.3.pdf 3.2节），通过m文件Test_Phase_Ex_ZPPP.m可以进行调用。
相关文献：
[1] Esmaeil Seraj and Reza Sameni, \Robust Electroencephalogram Phase Estimation with Applications in Brain-computer Interface Systems,”Physiological Measurement, vol. 38, no. 3, pp.501{523, Feb. 2017. [Published Online] DOI: 10.1088/1361-6579/aa5bba
2.1.3 传递函数扰动相位估计法
利用背景（自发）和前景EEG之间的加性数据模型，推导了瞬时EEG包络、相位和频率的概率密度函数和其他统计特性。分析和数值计算表明，在低解析信号包络中，脑电相位噪声大，易受背景脑电活动的影响。研究表明，尽管EEG相位变化传达了有关EEG的重要信息，但一些瞬时相位跳变是低分析包络中用于EEG相位提取的处理阶段的系统性副作用，与大脑无关。
在此方法中，利用被广泛接受的相位计算方法和统计信号处理的成熟方法，提出了一种窄带前景和背景脑电活动叠加的随机模型。利用该模型，导出了脑电信号瞬时包络（IE）和瞬时相位（IP）的概率密度函数。提出了一种Monte Carlo估计方案，用于相位相关参数的精确计算和平滑。
在工具箱中实现此功能的是Phase_Ex_TFP函数（函数具体说明见User_Guide_Ver2.3.pdf 3.3节），通过m文件Test_Phase_Ex_TFP.m可以进行调用。
相关文献：
[1] Reza Sameni and Esmaeil Seraj, “A Robust Statistical Framework for Instantaneous Electroencephalogram Phase and Frequency Analysis,” Physiological Measurement, vol. 38, no. 12, pp.2141{2163, Nov. 2017. [Published Online] DOI: 10.1088/1361-6579/aa93a1
2.2 通过实小波变换提取相位
文章：实小波变换提取相位信息方法研究
理论上快速离散小波变换(RDWT)的近似系数 、细节系数 中包含被分析信号的所有信息，包括相位信息，因此相位信息可利用RDWF提取。这篇文章提出三种从RDWT的 系数中提取相位信息的方法，而且对每个尺度上的小波变换系数进行正交化处理，便可借助小波域相子提取该尺度上的一种频率分量的真是相位。
该文章对普通实小波变换的实系数进行适当处理，得到一组实部、虚部具有正交性的复系数，再由这组复系数获得相位信息。这样做的好处是，可利用现有简单、快速、信息非冗余的实小波变换，通过简单的运算获得小波变换域的Hilbert变换对。
需要注意的是，某个小波变换系数对应时窗中的被分析波形包含多个频率分量时，该小波系数提供的幅值、相位信息，只能是这些频率分量的幅值、相位值的加权平均值，并且随小波而变，因为不同小波的权值不同。因此，只有用小波变换分析单一频率信号、分段单一频率信号或每个频率分量的频域间隔大于所用小波变换的频域分辨率的信号时，才有可能提供被分析信号的不同频率分量的实际幅值、相位值。
2.3 基于脑电自回归（AR）预测的实时相位估计方法
文章：基于脑电自回归预测的实时相位估计方法
文章的目的是实现闭环锁相刺激系统，其关键是实时获取当前脑电活动的相位信息。然而，由于系统软硬件的延时，采集到的 EEG 信号并不是实时的，实时的 EEG 信号通常采用预测的方式获得。因此，刺激信号的生成过程包括先计算采集到的一段 EEG 数据的平均周期，根据周期与相位的关系预测得到待刺激点的相位，通过延时控制刺激器在预设的相位时间点上输出刺激信号。
该文章提出了一种应用于闭环锁相电刺激系统的相位估计方法。该方法通过自回归模型进行 EEG 信号预测，采用阈值法识别预测信号的相位特征点; 依据预测相位特征点及前一个已知相位特征点之间的时间差和相位差，计算出两个相邻相位特征点时间间隔(TIBP)内任意时刻的相位值。该方法与平均周期法相比，具有更高的瞬时相位预测精度。
基于自回归的相位预测方法（phase prediction based on autoregressive，PP-AＲ) 主要包括EEG信号的窄带滤波、自回归建模、相位计算及相位特征点识
别4部分。
PP-AR算法的实现流程如下图：
文章中使用的采样频率为1kHz。
2.4 FFT相位估计法
文章：A Fast EEG Forecasting Algorithm for Phase-Locked Transcranial Electrical Stimulation of the Human Brain
https://www.readcube.com/articles/10.3389/fnins.2017.00401
该方法是通过首先在一个定义的时间窗口内分析记录的EEG信号，然后预测未来的信号。该预测方法基于这样一个假设：在一个小的时间窗口内，记录的EEG信号感兴趣的频带中的主导振荡的相位和频率将保持大致稳定。一个持续时间为 的EEG信号段用于预测持续时间为 的刺激信号的相位和频率。
如果从适当大小的时间窗口 中提取EEG信号的相位和频率，则可以对EEG信号进行准确预测。
预测的流程如下图所示：
首先，使用十阶椭圆IIR滤波器对EEG信号进行带通滤波，保留需要的频带，使用MATLAB“ellip”函数设计了一个10阶滤波器，通带纹波为0.5dB，阻带衰减为40dB。为了避免采样频率远高于滤波频率时经常遇到的浮点不稳定性，使用MATLAB“ss2sos”函数将滤波器传递函数转换为二阶部分，并利用MATLAB“sosfilt”函数按顺序使用。
其次，计算该信号段的快速傅里叶变换（FFT）；
第三，根据FFT计算信号主分量的频率和相位；FFT计算基于MATLAB的“FFT”函数。对一万个采样点的信号进行零填充，以减少FFT中的bin sizes，通过插值提供更高的频率分辨率。信号的频率由具有最大幅度的FFT单元的频率确定；信号的相位由其复值的角度分量给出。
文章中实验采样频率选择为500 Hz。如果用于锁相的目标EEG信号在从delta到gamma频率（即<60–80 Hz）的较低范围内，则较高的采样率不太可能带来额外的精度，并且可能会将算法的速度降低到阻止实时调谐的程度。
2.5 ETP相位预测算法
文章：Experimental evaluation of methods for real-time EEG phase-specific transcranial magnetic stimulation
https://iopscience.iop.org/article/10.1088/1741-2552/ab9dba
Github：https://github.com/OpitzLab/CL-phase
注：此文章不仅提出了一种新的相位预测算法ETP，还与两种已存在的算法进行了比较，包括FFT和AR算法（上面提到的两种方法），三种算法在Github里都有程序实现
该算法名为Educated Temporal Prediction（ETP），利用实时应用程序之前的一个简短训练节来学习每个个体的目标大脑振荡统计特征。在当前的实现中，ETP算法依赖于相位数据的基本统计特征，即时域中峰间间隔的中心时刻。
算法流程如下：
1、记录三分钟的实时脑电数据；
2、对想要预测的区域使用拉普拉斯蒙太奇（Laplacian montage）；
3、在目标的频带内，使用零相位FIR滤波器对数据进行滤波；该信号被用作后续步骤的真值（ground truth）。
4、识别训练数据（学习阶段）前90秒的信号峰值，并计算每对相邻峰值之间的间隔（Tn）。为了确保发现的峰值是有意义的，并排除由非理想滤波引起的可能的小波动，引入了最小峰值距离的标准（即，α振荡为62.5 ms）；
5、选择这些峰值之间周期的中位数作为峰间间隔（T），作为峰值预测的初始估计值；
6、第二个90秒原始静息脑电图数据用于模拟估计相位预测（验证阶段）的准确性。在验证阶段，调整T的值，以获得该部分数据集的最佳阶段精度。为此，选择了250个长度为500毫秒的重叠窗口；
7、应用目标区域的拉普拉斯蒙太奇，并使用非因果sinc滤波器（砖墙滤波器；频域中的理想锐滤波器）在感兴趣的频带中单独过滤窗口。
8、从窗口末端移除边缘（alpha为40 ms），以避免过滤器造成的边缘伪影；
9、检测到该窗口中具有步骤4中解释的标准的峰值；
10、通过将T添加到最后检测到的峰值来预测下一峰值的时间；
11、根据真值（ground truth）的希尔伯特变换测量预测峰值处信号的实际相位。
12、为了使相位检测中的潜在偏差最小化，T的值以增量的方式增加或减少，并使用新的T重新运行验证阶段（步骤6–11），直到与目标相位的偏差最小化。优化调整值 被用来进行实时相位估计。
文章中使用的脑电数据采样频率为1kHz，分别通过FFT、AR、ETP算法进行处理。
2.7 一种基于机器学习的脑电节律相位估计方法
文章：McIntosh, J. R., & Sajda, P. (2019). Estimation of phase in EEG rhythms for real-time applications. arXiv preprint arXiv:1910.08784
(https://arxiv.org/abs/1910.08784)
Github：GitHub - jrmxn/py_eegepe: Example methods for estimation of phase of brain-rhythms for real-time use
使用机器学习方法模拟非因果信号处理链和纯因果信号处理链。能够从经过非常小的预处理的脑电图信号中估计瞬时相位，其精度高于更标准的信号处理方法。具体来说，我们首先学习最小预处理的EEG信号和非因果处理的瞬时相位之间的关系。然后，利用这种学习到的关系来预测事件发生前瞬间的阶段。这种方法比基于非机器学习的方法提供了更好的相位估计。以因果方式映射原始信号和非因果恢复状态之间的关系，用于实时操作和一般离线分析。
流程如下：
2.7.1 预处理
首先导入数据，然后在50Hz下对所有EEG数据进行非因果低通滤波，以去除任何残留的线路噪声，并重新采样至512Hz（论文中使用的数据集采样速率为2048Hz）。 然后应用了以POz为中心的Hjorth[1]，这意味着从POz中减去PO3、PO4、Pz和Oz处的信号平均值。数据通过除以30-50Hz频段的平均功率进行标准化，以使受试者的噪声水平保持一致。通过在8-12Hz范围内取最大功率下的频率来估计每个受试者的单个α频率（IAF）（论文中讨论的是α节律的相位估计）。
预处理后，处理的信号沿两条路径分开。在一个分支上，数据被分段以生成输入信号（用于训练的特征）。在另一个分支上，我们应用FIR带通滤波器（8-12Hz通带，769阶，约1.5s），并在分段之前进行希尔伯特变换，以生成输出信号（用作我们评估相位估计的信号）。
2.7.2 基于预测的学习
考虑利用频率信息的非线性网络方法，该网络的构造使其能够学习抑制一些滤波器输出，并增强其他滤波器输出，这取决于通过引入多个密集层（能够通过softmax操作相互选通）而呈现的频率混合。尝试使用选通循环单元（GRU）[2]，这是一种循环神经网络（RNN）如下图。
文章中使用的采样频率为2048Hz。
注：关于RNN的结构论文里没有仔细说明，所以对这个方面还有些疑问。另外程序的内容较丰富，暂时还没有搞清楚各个程序的作用和使用方法。
参考文献：
[1] Bo Hjorth. An on-line transformation of EEG scalp potentials into orthogonal source derivations. Bo Hjorth. An on-line transformation of EEG scalp potentials into orthogonal source derivations.
[2] Junyoung Chung, Caglar Gulcehre, KyungHyun Cho, and Yoshua Bengio. Empirical evaluation of gated recurrent neural networks on sequence modeling. arXiv preprint arXiv:1412.3555, 2014.
