时频分析方法及其在EEG脑电中的应用

EEG提供了一种测量丰富的大脑活动即神经元振荡的方法。然而，目前大多数的脑电研究工作都集中在分析脑电数据的事件相关电位（ERPs）或基于傅立叶变换的功率分析，但是它们没有利用EEG信号中包含的所有信息——ERP分析忽略了非锁相信号，基于傅里叶的功率分析忽略了时间信息。而时频分析（TF）通过分离不同频率上功率和相位信息，可以更好地表征脑电数据中包含的振荡，TF提供了对神经生理机制更接近的解释，促进神经生理学学科之间的连接，并能够捕获ERP或基于傅里叶分析未观察到的过程（如连通性）。但是，本文献综述表明，脑电时频分析尚未被发展认知神经科学领域所广泛应用。因此，本文从概念上介绍时频分析，为了让研究人员便于使用时频分析，还提供了一个可访问脚本教程，用于计算时频功率（信号强度）、试次间相位同步（信号一致性）和两种基于相位的连接类型（通道间相位同步和加权相位滞后指数）。

简介

       EEG（EEG）是一种成熟且实用的工具，用于研究大脑功能，心理学和精神病学的发展。迄今为止，大多数发育性EEG研究主要集中在事件相关电位（ERP）分析或基于傅立叶变换的功率分析。虽然这些方法已经证明是有效的，但它们没有利用EEG信号中包含的所有信息。即ERP分析忽略了非锁相信号，基于傅里叶的功率分析忽略了时间信息。时频（TF）分析可以更好地表征EEG数据中包含的三个振荡特征的时间动力学:频率、功率和相位。本文的目标是为发展性EEG研究人员提供一个全面、易懂的教程和自动化但易于调试的Matlab脚本，来计算TF功率、试次间相位同步（ITPS）和两种类型的相位同步（通道间相位同步（ICPS）和加权相位滞后指数（wPLI））。

1 为什么要做时频分析

       在从婴儿到成年的整个生命周期中，EEG可以在实验室和移动环境以相对低的成本收集，同时保持良好的时间分辨率，对运动和噪声具有相当强的鲁棒性。在过去的几十年里，EEG被用于研究认知、社会情绪能力和精神病理学的发展。历史上，由于计算能力的限制，大多数EEG研究集中在ERP或基于傅立叶的功率分析。

虽然ERP和基于傅立叶的功率提供了丰富的信息，并为我们理解心理现象的发展提供了重要的见解，但ERP和基于傅立叶的功率并不能充分利用EEG信号中的所有信息。ERP分析假设感兴趣成分在不同trials中是时间同步的，只关注与感兴趣事件时间锁定的神经活动，而忽略在不同trials中不同步的信号。例如，如图1，假设有两个ERPs，一个在0 ms，另一个在500 ms。如果这些组成部分在不同trials的延迟中略有不同（即，在不同trials中不是时间一致的），在我们平均trials时在0 ms时的ERP活动将作为噪声丢失。相比之下，在500 ms时，第二部分在各个trials中是完全同步的，并且当我们在各个trial中平均时，可以清楚地保留下来。虽然我们的例子是假设的，但在不同trials中大脑反应并不总是符合时间同步的假设，在研究发育过程时尤其值得怀疑。因此，在发育人群中研究ERP时，尤其是在比较不同年龄的ERP时，考虑到这种差异是特别重要的。

图1 ERP和TF分析比较 

2 用振荡表征大脑活动

有一种假设是将大脑活动表征为振荡，而不是特定时间的ERPs。将EEG表征为振荡的优势在于，大脑活动可以被表征为几个独立的测量指标，如频率、振幅和相位。在图2A中，我们展示了一个每秒完成两个完整周期的正弦波，所以它的频率是2 Hz。相反，图2D所示的正弦波振荡更快，每秒完成4个周期，频率为4 Hz。这些正弦波的另一个特点是它们在0附近从1到 -1循环。振荡的高度代表振幅，振幅是平衡点到最高点和最低点的距离，所以图2A和图2D所示的波的振幅都为1。另一方面，图2B所示的振荡频率与图2A相同，但振幅较小，为0.5。最后，我们可以根据相位来描述振荡。相位是振荡在特定时间（通常为0）上的位置。这样，我们可以使用相位来估计振荡相对于特定时间、事件或其他振荡的对齐情况。例如，图2A和图2B中的波具有相同的频率和相位（即它们是完美排列的），因此它们都从振荡的峰值开始。这些通常被称为锁相信号。另一方面，图2C所示的波与图2A所示的波相同，但它从平衡态或零点开始发生了位移。相位以360度的角度正式测量，代表一个完整的振荡周期，用极坐标图说明。如图2所示，与图2A和图2B相比，图2C所示的正弦波为90度（或360度）。这些相移的信号彼此相比较，通常被称为非锁相信号。

图2 振荡特性

用振荡来表征大脑活动最常见的方法是用傅里叶变换测量给定频率下的活动量。该方法通过将不同频率的正弦波与EEG数据卷积来测量其在给定频率下的振幅。在给定频率下，信号的能量通常用功率来测量，功率是振幅的平方。然而，它假设信号是平稳的，失去了脑电信号中丰富的时间信息。

TF分析可以测量不同频率神经振荡的振幅和相位的动态变化。通过区分振幅和相位信息，锁相和非锁相信号都可以对感兴趣的事件进行研究。例如，回到图1，我们在0 ms有一个6 Hz的非相位锁定振荡，ERP没有捕捉到它，但在TF表示中可以清楚地观察到，在0到200 ms之间，6hz的信号强度或功率增加了。类似地，TF分析捕捉到500 ms时功率和频率增加时的10 Hz锁相振荡。作为一个单独的测量，TF分析还提供了给定频率的相位随时间的估计。可以检查该阶段信息在各个trials中的一致性或同步性,即试次阶段间同步（ITPS）。在图1的例子中，非锁相分量在6Hz时产生的ITPS值为约从0到200毫秒。相比之下，在10Hz的锁相分量的ITPS为1.0，表明在整个trials中完美的相位一致性。除了检查各试次间相位振荡的一致性（ITPS）， TF分析还可以检查不同电极间相位振荡在各试次或时间的一致性（ICPS）。

3 TF分析的优势

首先，TF分析相对于其他EEG方法的一个主要优势是其可解释性。因为神经元振荡是大脑的基本属性，TF分析提供了关于EEG数据捕获的过程下的神经生理机制更直接的信息。这也促进了神经生理学多学科的连接。其次，因为TF分析提供了信号强度（功率）、各试次阶段一致性（ITPS）和连通性（ICPS）的单独测量，它们可能为神经认知过程的发展提供了独特的理论。最后，有一些处理过程只能通过TF来捕获。比如捕捉神经元系统活动的EEG测量，称为mu去同步。Mu节律反映了运动皮层（即C3和C4）发生在α范围内的EEG活动（成人为约8-13 Hz，婴儿和儿童为约6-9 Hz），当个体执行或观察一个动作时，其振幅下降。在一个婴儿样本中，科学家通过TF分析发现在执行和观察动作时mu不同步。

4 TF分析的计算

4.1卷积

TF分析主要操作是卷积运算。卷积是一个将两个信号结合起来产生另一个信号的过程，即两个信号之间相似性的度量。在数学上，卷积先是做点积运算——两个信号中元素与元素的乘积，然后它们的乘积相加：

其中a是第一个信号，如EEG信号，b是第二个信号，是一个具有预定义或已知函数的信号，称为核函数。对于基于傅里叶的分析，使用的核函数是一个正弦波。对于TF分析，使用的是时变正弦波——小波，一种振幅从零开始，增加，然后减少的振荡。最常用的小波是Morlet小波或Gabor小波，Morlet小波是正弦波和高斯窗的组合，Morlet小波及其与正弦波和高斯窗的关系如图3所示。

图3 Morlet小波及其与正弦波的关系

卷积就是在一段时间内重复地做点积。如图4所示，通过在EEG信号上滑动或移动小波来实现的。例如，在图4中，我们将图1的第一次trial和一个6 Hz的小波进行卷积，捕捉所有6 Hz的振荡，这些振荡在0 ms附近更大。对于TF分析，卷积实际上是用复Morlet小波来进行的。将EEG信号与复数Morlet小波进行卷积，得到一个复数的时间序列。在数学上，卷积过程可以用以下公式来表征：

其中M对应振幅，其中ϕ表示相位的角度，e是自然对数的底，i是虚单位。此外,所有的项都有下标t,表明这是将整个EEG信号应用小波,建立一个频率f的时间序列。虽然我们用6 Hz的小波来说明卷积过程，但这个过程可以用不同特征的小波（例如，不同频率）来重复，从而允许研究人员在时间和频率上捕捉EEG信号的振幅和相位信息。因此，仔细选择小波的特征是很重要的。

图4 使用复数Morlet小波进行卷积，以测量每个频率随时间变化的振荡幅度和相位

4.2 时频功率

       分析TF功率是研究人员使用TF分析最常用的方法之一。为了测量不同频段内一段时间内TF功率的时间动态，先利用时频分解来隔离信号在每个时间和频率的幅值。然后将每次trials的TF分解平均起来，以显示总功率。研究人员通常使用分贝转换对原始总功率值进行基线归一化：

其中，用时间t和频率f的trials平均活动除以频率f的基线周期的平均活动，log10转换后再乘以10，产生分贝（dB）。在基线标准化后，所有活动将在相同的尺度上，任何与平均基线时段共同的活动将被移除，从而把TF功率解释为相对于基线时段的活动水平。因此，研究人员应该仔细选择基线时间段。

与ERP相比，TF分析需要考虑TF分析的基线窗口大小，窗口通常应该在相关事件开始之前。TF分析通常使用几百毫秒，并在感兴趣的事件之前结束。理想情况下，基线条件应该与感兴趣的时期具有相似的长度。最后，除了分贝转换，还有百分比变化或Zscore等方法做基线标准化。

4.3 相位同步

TF分析的另一个主要结果是EEG信号的相位信息，相位信息可以在不同频率的每个时间点上被捕获。这个相位测量需要在各个trials中进行检查，以捕捉振荡的相位在各个trials中是如何一致或同步的。相位同步可以通过平均各个trials的相位角值来估计，因为相位角是圆的，而不是线性的，相位角在一个具有实虚分量的复平面中表示，相位角差的平均值通过矢量数学计算。结果向量是一个介于0和1之间的值，较低的值表示较低的一致性，较高的值表示较高的一致性。

4.3.1 试次间相位同步（ITPS）

ITPS是在整个trials中检查特定时间和频率上相位一致性的度量。ITP捕获的一致性EEG活动感兴趣的事件在一个特定时间和频率较高的ITP值反映更多的一致性（1 =完美之间的一致性的）,或者减少ITP代表接近随机相位对齐锁相事件感兴趣的阶段（0 =随机频率的观点试次）。在数学上，ITPS的计算公式如下：

其表示时间点t和频率f上试次x上n个trials的相位角向量的平均值。发展研究人员使用ITPS检查大脑对特定事件反应的一致性，以检查错误处理、奖励处理和听觉辨别。

需要注意的是，在条件或参与者之间的trails数量很大不同的情况下，ITPS的估计可能是有偏差的。对于trail较少的情况，可以人为地夸大该条件/参与者的ITPS估计。因此，在不同的条件下采用相同的试次次数是很有用的。此外，为了进一步减轻试次数量不平衡的影响，并隔离感兴趣的事件，也可以做ERP类似的基线修正。

4.3.2 相位同步作为连接性的衡量标准

量化通道间相位同步（ICPS）与ITPS类似，但ICPS计算不同通道间的相位同步，而不是检查试次间相位的一致性。如果两个通道的相位角随时间同步波动，它们的差异将保持不变，产生接近1的ICPS值。相反，如果两个通道的相位角随时间随机波动，它们的差异将会很大，产生接近于0的ICPS值。ICPS测量的形式如下：

它表示在时间点t和频率f时，试次x上n次trails的通道j和通道k之间的相位角（其中n）差异的平均值。

在电极间进行脑电分析的一个局限性是脑电的空间分辨率较低。由于电容传导，由一个信号源可以传播到头皮的不同区域，在头皮上一个位置记录的读数有多个活动。为了计算具有更大空间特异性的ICPS，必须在计算ICPS之前对数据进行这些转换，比如加权相位滞后指数（wPLI），wPLI对滞后的幅度进行加权，而不是去除所有0或180度相位差，减少了可能由于体积传导而对噪声不敏感的估计，并与真实相一致性显示出更可靠的关系。

这些基于阶段的连接度量（ICPS和wPLI）可以在头部的每个电极之间进行计算，称为all-to-all连接。相反，如果研究人员有一个关于两个特定区域之间连接的先验假设，ICPS和wPLI可以从一个seed电极计算到其他感兴趣的电极，称为基于seed的连接。最后，除了通过trials对ICPS和wPLI进行量化外，这些基于阶段的测量还可以对每个试次进行随时间的平均估计，称为随时间的连接性，而不是通过trails。

5 TF分析的挑战和局限

首先，在实际操作TF中，神经信号的这些不同特征的计算是计算密集型和耗时的，需要使用高性能计算机或高性能计算集群来运行这些分析。

其次，TF分析产生了大量特征，需要对多个特征比较。可以使用现有的文献来先验地定义不同的感兴趣区域，或者使用非参数聚类方法来考虑样本点之间跨地形、时间和频率的依赖性，以及使用降维技术，如主成分分析，以捕获有意义的活动。

TF分析的另一个局限性是频率和时间之间的权衡。为了改善这一问题，可以使用不同类型的时频分解，如Cohen类减少干扰分布，它在时间和频率上产生一致的高分辨率。

在进行TF分析时要考虑的另一个重要问题是振荡的性质。在本文中，我们假设神经元振荡具有正弦形状，但大脑振荡可能并不完全像正弦波，大脑活动，特别是在更高频率如贝塔和伽玛（15 Hz）下，可能会以burst形式发生。一个解决方案是检查trial-by-trail的试次活动，以确保支持对平均功率的解释。

因为TF分析在发展认知神经科学研究中相对较新，存在较少的标准（例如，最小trials数量）。一项对4- 9岁儿童对错误反应的研究的初步发现发现TF分析和ERP需要相同数量的试次来观察错误与正确对比的影响。根据测量方法、兴趣条件和参与者的年龄，范围从8到超过32个trails，以达到可接受的内部一致性可靠性（0.60）。

最后，要强调，研究人员必须检查在他们自己的数据中计算出的TF心理测量特性，并认识到TF心理测量在不同任务、不同年龄和不同人群中的潜在细微差别。

6 使用Matlab计算TF

TF度量可以使用各种软件进行计算，包括基于脚本和基于GUI的程序，如FieldTrip和EEGLAB。我们选择在Matlab中实现我们的TF分析脚本,以便用户可以轻松地根据特定需求编辑代码。我们的脚本使用EEGLAB数据格式，主要基于Cohen（2014）。

输入的数据是清洗过的epoched数据。初始TF分析脚本（mainscript_calculate_tf_itps_icp.m）包括各种设置和参数供研究人员选择。如图5所示,这些脚本开始与用户输入的路径表明：

1） data_location路径文件夹包含数据的分析,

2） save_location路径文件夹,他们想要保存的数据,

3） scripts_location路径脚本的位置,

4） EEGLAB位置的路径（关于如何安装这个工具箱，请参阅EEGLAB的文档）。

接下来，用户输入一些关于数据集和分析过程的信息（图5）。研究人员设置通道的数量（nbchan）以及它是静止的还是事件相关的数据（RestorEvent）。如果范式包含多个条件，则研究人员还需要指定每个条件的名称（Conds）。为了方便起见，我们还提供了一个脚本（Edit_events.m），该脚本将在运行这些脚本所需的命名约定中重新标记感兴趣的事件。研究人员还需要指出参与者在要分析的条件中必须进行的最小试次次数。如果参与者的条件之一不满足此阈值，则将为该参与者保存一个名为“notenoughdata.mat”的空白文件，并且脚本将跳到列表中的下一个参与者。在这里，研究人员还需要输入他们是否希望数据被基线校正（BaselineCorrect）以及基线周期的时间段。

图5 “MainScript_Calculate_TF_ITPS_ICPS.m”脚本的输入设置和数据集信息

研究人员还需要决定是否要对数据进行降采样。在脚本开始时，数据将自动向下采样到250 Hz，以减少文件大小和计算时间;然而，如果研究人员想在TF计算后进一步降采样到125 Hz，以减少文件大小和存储空间，而不损失大量的分辨率，他们可以选择再次降采样。我们建议在TF分解后下采样到125 Hz，特别是对于具有许多条件和/或大样本容量的数据集。最后，研究人员可以选择将研究标识符的扩展名附加到最终的文件名中

接下来，研究人员必须确定TF分析的参数（图6）。首先，用户需要确定要分析的最低频率。这一决定受到几个因素的影响，包括先验的兴趣事件长度和epoch的长度。其次，用户也要选择频率最高的进行分析。

图6 MainScript_Calculate_TF_ITPS_ICPS.m脚本中输入计算TF功率、ITPS 和 ICPS/wPLI 的设置和参数 

第三，研究人员需要确定频率之间的分辨率，或者在最低频率和最高频率之间应该有多少个频率bins。最后，研究人员需要确定小波的宽度。小波的宽度对时间和频率精度都有影响。更长的宽度增加频率精度，而更短的宽度导致增加时间精度。

最后，研究人员需要设置基于相位分析的参数（图6）。首先，他们将决定是否计算或跳过ITPS。研究人员还将决定，如果他们的条件有不同的试次次数，他们是否愿意使用分抽样程序来等同试次次数。对于次级抽样，他们还需要决定

1） 对每个次级抽样进行多少次随机抽样，

2） 次级抽样的数量。

在选择这些参数时，我们建议研究人员确保他们选择的数字有可能使所有数据都被使用。例如，如果在一种情况下有300个试次，但在10个子样本中，每个样本只进行了5个试次，那么很多数据就被遗漏了。子采样程序设置好后，用户将决定是否计算或跳过ICPS/wPLI计算。如果选择了ICPS，我们建议对数据进行CSD或拉普拉斯变换。研究人员还需要确定ICPS/wPLI是否应随着时间或试次进行计算。研究人员还可以选择是否在所有电极之间（all-to-all）或在seeds电极和其他特定电极之间计算ICPS/wPLI。最后，研究人员需要使用正则表达式（例如，搜索索引将使用的字符序列）作为脚本，以查找数据位置路径中感兴趣的文件，并创建要循环浏览的列表。选择这些设置和参数后，MainScript_Calculate_TF_ITPS_ICPS.m脚本将使用timefreq.m、ITPS.m、ICPS.m 和 wPLI.m 脚本计算感兴趣量，并保存 Matlab （.mat） 文件以及每个条件和目标的最终度量值。

除了用于计算各种 TF 度量的脚本之外，我们还提供了另外两个脚本。一个脚本（Edit_events.m）将确保以易于读入脚本的方式标记EEG文件中不同条件下的事件/标志。在通过主脚本运行此脚本之前，应先对所有epochs文件运行此脚本。其次，在 TF 脚本计算出适当的度量值后，名为 Compiling_and_plotting.m 的最终脚本将创建一个矩阵，该矩阵按时间按通道编译所有主体和条件的频率数据。然后，脚本使用此矩阵为所有主体的平均值绘制 TF 表面和地形图。

Compiling_and_plotting.m脚本要求研究人员再次输入一些信息（图 7）。首先，研究人员将输入三个路径：

1）data_location路径，该路径应指向包含每个目标和主脚本计算的条件的.mat文件的文件夹，

2）save_location路径，应指向研究人员希望保存编译数据的文件夹，

3）将路径添加到EEGLAB的位置。

接下来，研究人员指定数据是静止的还是与事件相关的，Edit_events.m脚本中指定的条件名称以及条件的通用名称。研究人员还将根据主脚本执行的内容指出有关 TF 计算的详细信息。首先，他们将选择其TF数据是否经过基线校正（TF_bln）。他们还将指示他们是否选择计算 ITPS、ICPS 或 wPLI（ITPS_calc、ICPS_calc 或 wPLI_calc），以及他们是否选择基线校正这些度量值（ITPS_bln、ICPS_bln 和wPLI_bln）。最后，他们将决定是否选择对数据进行降采样。

图7 由研究人员在Compiling_and_plotting.m脚本中输入设置、数据集信息和绘图参

最后，研究人员将做出一些关于绘图参数（图 7） 画surface 图（即按频率计的时间）和头皮地形图（基于通道导联的 TF 测量头皮位置值）。他们需要选择在曲面图的 x 轴上绘制的时间窗口（time_window_surface）。通常，此时间窗口应包括感兴趣事件之前和感兴趣事件之后的时间。他们还需要选择在曲面图的 y 轴上绘制的频率窗口（freq_window_surface）。时间窗口和频率窗口的选择允许研究人员自定义其图以说明感兴趣的效果。此外，研究人员需要决定他们想要绘制表面的电极（或电极簇）。他们需要分别为TF（chans2plot_TF），ITPS（chans2plot_ITPS）和ICPS / wPLI（chans2plot_ICPSwpli）执行此操作。为了以类似的方式绘制all-to-all的ICPS或wPLI，研究人员必须选择Seed电极。在绘制ICPS或wPLI的表面图时，研究人员应仅绘制非Seed电极。

除了表面图的参数外，研究人员还将为头皮地形图设置一些参数。同样，研究人员需要输入要在头皮地形中绘制的时间窗口（time_window_for_topo）和频率窗口（freq_window_for_topo）。虽然表面图通常包括较大的时间和频率窗口（例如，整个epoch和1-30 Hz），但头皮地形图通常会放大特定感兴趣事件的时间段和频带（例如，0-100 ms和4-8 Hz）。

7 例子

接下来，我们将通过TF surface plots和地形图来说明我们的TF处理脚本的输出，这些图来自于一项使用奖励提示侧翼者范式研究奖励对青少年男性认知控制的影响的研究。为了说明这一点，我们将重点研究25名男性儿童的不一致刺激的刺激锁定数据（Mean age = 10.55, SD =1.10, Range = 9.08 − 13.92）。所有参与者均获得父母同意和儿童同意。所有数据均使用马里兰州发育EEG（MADE）流水线分析进行预处理。

在固定刺激的ERP图中（图8A），在前额中心区域有一个大约300 ms的明显的负偏转，表明N2的存在。在不一致的情况下，TF功率（图8B）和ITPS（图8C）表面的theta （4-8 Hz）从大约0-400 ms升高。通过检查（4-8 Hz）功率和ITPS在整个头皮的分布情况，我们发现额中央和枕部区域的信号强度（功率）和一致性（ITPS）增加。作为连接的一个例子，我们使用了基于seed的wPLI，它带有一个额中心seed（电极E6）。结果显示，从0到300 ms的θ波的连通性增加。值得注意的是，头皮地形图看起来有些不同，因为这些地形图与感兴趣的seed（E6）有关。因此，观察到的模式可能代表了在不一致条件下增加的双边连接额叶区域。

8 总结

    本文主要做了TF分析的介绍，强调了这种方法与传统的ERP和基于傅里叶功率相比，对发展认知神经科学研究的独特贡献。ERP分析忽略了非锁相信号，基于傅里叶的功率方法忽略了时间信息。TF分析分离了特定频率下EEG数据中包含的功率和相位信息，为EEG信号提供了更全面的表征。通过回顾利用TF分析发展性EEG数据的新兴研究，我们强调了这些分析在增进我们对发育过程的理解方面的潜力。我们还讨论了TF方法面临的挑战和未来的发展方向。最后，为了促进发展认知神经科学社区使用TF分析，我们提供了一个易于使用的Matlab脚本教程，基于Cohen（2014），以青少年样本计算最广泛使用的TF测量值。

图8 TF分析的例子

 参考文献：Time-frequency analysis methods and their application in developmental EEG data；