导读
近年来，人类睡眠脑电图(EEG)研究激增，采用了越来越复杂的分析策略将电生理活动与认知和疾病联系起来。然而，正确计算和解释当代睡眠EEG中使用的指标需要注意许多理论和实际的信号处理细节。本研究回顾了与频谱分析、蒙太奇选择、相位和振幅信息提取、替代结构相关的几个方法学问题，说明了方法选择对结果的影响，以及通过可视化和简化示例检查处理步骤的重要性。本文以非数学的形式，结合睡眠特定示例，以及代码实现来展示这些问题，旨在使新手和非技术性研究人员能够对方法学有更深层次的理解，从而帮助提高未来睡眠EEG研究的质量。
前言
睡眠参与认知和神经精神疾病的证据越来越多，将这种独特的大脑状态推向了基础和临床神经科学的前沿。虽然可以使用多种技术测量睡眠大脑，但EEG仍然是人类睡眠研究的首选方法。随着EEG硬件(如高密度放大器和脑电帽)和分析程序(如时频、功能连接、跨频率耦合分析)的进步及其在神经科学领域中的广泛使用，这些硬件和分析方法也开始渗透到睡眠领域。事实上，通过同时捕获电生理活动的时间、频谱和空间方面，现代EEG技术(结合其他神经科学工具)在理解非快速眼动(NREM)和快速眼动(REM)睡眠的复杂组织和功能方面具有重要作用。例如，NREM活动的某些方面(慢波振荡(SOs)，睡眠纺锤波及其耦合)与衰老和认知衰退、精神分裂症、创伤后应激障碍、自闭症、以及健康个体的认知、记忆和情绪功能方面有关。
尽管人们对睡眠EEG的兴趣日益浓厚，但仍需要大量的理论和实用的信号处理专业知识来正确分析睡眠EEG并评估已发表结果的有效性。然而，信号处理培训在最可能进行或评估睡眠脑电图研究的领域(如心理学、神经科学、医学)中并不常见，有时会对已发表结果的质量产生负面影响。虽然犯错对于提高分析技能的过程中是不可避免的，但对方法问题的认识不足可能会完全阻止错误的发现。再加上统计不严谨，分析工具的不当使用将产生错误和不可复制的结果。本文通过相关示例和代码实现回顾了一些方法论问题，希望为非技术受众提供一个睡眠EEG分析的切入点。
本文有两个目的。首先，强调不起眼的分析选择如何对结果测量产生重大影响，进而对基础和临床研究的解释也产生重大影响的。虽然分析选择会产生一些下游影响是微不足道的，但我们经常对结果的变异性程度感到惊讶。同时，在检查大量潜在有趣的结果变量时，以不同的方式重复分析将增加假阳性的风险。其次，强调检查处理步骤是否产生所需输出的重要性。分析脚本可能运行平稳，产生合理的结果值，并生成引人注目的图表，但根本没有执行研究人员设想的计算，这可能会影响研究结果。这些问题很难追踪，但通过可视化中间结果并将分析应用于剥离或模拟数据，可以极大地促进这些问题的解决。为了实现这两个目的，本研究讨论了睡眠EEG研究中常见的一组数据分析程序，提请注意经常被忽视的问题，并指出定义中的歧义。附带的Matlab脚本进一步说明了这些问题，并可作为自定义分析的起点，解决基本和转化睡眠EEG研究中的问题。
就范围而言，首先，本文不提供所有相关方法学问题的详尽概述，也不假设涵盖了最重要的问题。虽然所选的问题当然希望其具有实际用途，但它们的主要目的是说明上述目标。其次，许多涉及到的问题早已在其他信号处理和EEG领域中广为人知，甚至也被睡眠专家所熟知，并且没有声称有新的见解。相反，本研究希望以非数学形式，从专门的睡眠EEG角度和代码实现来呈现这些问题，使非技术受众更容易理解这些问题。第三，虽然所包含的主题和示例源于我们在人类无创头皮水平睡眠EEG方面的经验，但其中许多问题同样适用于源重建和侵入性EEG记录、脑磁图(MEG)、其他物种和觉醒状态。第四，由于本研究的重点是数据分析，因而将忽略电极阻抗问题和通常在生成干净数据之前的许多预处理步骤。但是，有关滤波、伪影校正、通道插值等等的选择已经构成了分析链中的重要决策点。
代码和数据
本文附带的所有代码都是在Matlab 2018a(Mathworks，Natick)中编写的，需要用到Matlab的信号处理、统计和机器学习、曲线拟合工具箱，以及EEGlab功能。纳入的数据(58通道NREM数据分为N2和N3阶段，采样率为400Hz，带通滤波范围为0.3-35Hz)来自三名健康被试，他们接受了一整晚的多导睡眠脑电图记录(EEG、眼电图和肌电图)，根据美国睡眠医学学会指南对睡眠进行评分。合并的代码和数据(EEGlab格式)可以从http://doi.org/10.5281/zenodo.3712074下载。
频谱分析
频谱分析是描述具有节律成分信号的最著名方法之一，在睡眠EEG中有着悠久的历史。它提供了信号特性的广泛而直接的概述，包括频谱成分、它们的变异性和整体信号质量。虽然信息量很大，但分析方法的变化阻碍了研究之间的比较。事实上，“功率谱”这个术语本身是相当模糊的，因为它通常用于指代功率或功率谱密度(PSD)，这是两个密切相关但具有不同单位的截然不同的概念(PSD：μV2/Hz；power：μV2)。虽然PSD被认为是更合理的表示形式，但离散信号的PSD和功率谱(如EEG)通过一个简单的换算系数相关联，并且在大多数实际用途中(组间比较，将PSD/power与认知测量相关联)，它们可以互换使用。
原始和相对/标准化PSD
睡眠EEG的PSD/功率估计通常使用Welch方法(或Welch周期图)。这些估计值总是正的，在最慢和最快频率之间存在许多数量级的差异(图1A：i)。对于EEG信号，频率f和PSD/功率之间的关系近似于1/fa分布(其中指数a表示频谱下降的速度)，也称为幂律缩放。节律信号成分的存在表示为与该背景“1/f”活动的正偏差。请注意，虽然每个频率(频段)始终具有非零功率，但这并不意味着存在振荡活动。当不是正弦形状时，节律活动也可能由多个频谱峰值(谐波)反映。此外，频谱功率受各种伪影(例如，肌肉活动、眼动、汗液伪影)的影响。因此，虽然频谱形状提供了振荡节律存在与否的粗略指示，但其本身并不是决定性的。
图1.功率谱密度、功率和归一化。
由于1/f缩放，原始PSD/功率通常在对数y刻度上绘制(图1A：ii)，或对数变换后在线性刻度上绘制(图1A：iii)。虽然ii和iii的形状相同，但对数变换的PSD/功率在较高频率下为负值，而未变换的PSD/功率＜1，这可能导致不直观的行为(例如，较窄的0.5-20Hz范围的PSD/功率总和大于较宽的0.5-30Hz范围的PSD/功率)。因此，理想情况下，应该使用未转换的PSD值进行进一步的分析。当覆盖更宽的频率范围时，将x轴按对数间隔也会很有帮助(图1A：iv和v)。
在确定特定频段内的“功率”(例如，0.5-4Hz慢波活动或SWA)时，会产生其他复杂性。对于连续信号，功率定义为PSD函数相关部分的积分。对于离散采样信号，如EEG，可以通过对相关频率箱的PSD值求和并乘以适当的箱大小(或等效地，对箱之间的功率求和)来近似。但由于PSD和功率只是彼此的缩放版本，因此直接相加PSD值通常没有什么影响。许多研究人员对各个箱的值进行平均，而不是求和。虽然当频宽不相等时，用求和或平均法得到的PSD或功率的不同模式没有什么实际差别。此外，一些研究对对数变换值(而不是对未变换值)进行求和或平均。由于这些细节并不总是被记录下来，所以不清楚有限频宽活动是否涉及PSD或功率、原始或对数变换值、平均值、总和或积分值。
原始PSD/功率与信号幅值的关系相对直接，较大信号幅值的通道通常具有较大的功率，如图1B和C中所示的原始PSD。因此，当信号幅值的绝对差异被认为是有意义的(例如，地形分析)时，原始频谱是有用的。另一方面，与参考点的距离影响信号幅值时，个体通常会表现出相当大的信号振幅差异，至少部分是由于不感兴趣的因素(脑回折叠，脑电帽定位，颅骨形状和厚度等)造成的。因此，频谱通常按总功率进行归一化，从而指示每个频谱分量对信号的相对贡献。
无论相对功率的确切定义如何，功率归一化的后果都可能很大。比较图1C和D(左侧)的PSD频谱形状，很明显，与原始PSD相比，相对PSD的使用使频谱向上或向下移动。这不仅减少了SO频段中的通道差异，而且还影响了这些示例通道中显示最大的快速纺锤功率(原始：Pz；相对：Oz)。考虑到一些频率的PSD地形图(图1C和D，右侧)，相对PSD消除了0.8Hz SO活动的正面优势，并将13.3Hz快速纺锤体活性的热图向后移动。令人惊讶的是，原始和相对PSD在6Hz θ和25Hz β频率下产生完全不同的地形，分别在前部或后部头皮区域显示出最大功率。图1的代码可在powerDefined.m和powerNormalization.m中找到。
当使用原始或相对PSD比较N2和N3之间的频谱分布时，情况也类似。因此，在哪个睡眠阶段，人们认为功率最大取决于所选择的归一化。因此，应根据分析目标定制归一化。一般而言，归一化的基准(分母)不应随感兴趣的因素而变化。由于1/f分布，振荡分量有时可能表现为肩峰而不是清晰的频谱峰。增强峰值的一种简单方法是采用时域信号梯度(作为时间导数的粗略近似值)，并根据该信号估计PSD。尽管与原始PSD相比，这不会影响通道光谱和地形的相对顺序(图1E)，但它可能有助于频谱分量的分离。抵消1/f活动的更复杂的方法包括拟合该成分并从原始PSD中减去它。然而，到目前为止，还没有关于如何处理不同频道、不同个体和不同睡眠阶段的指南，但几乎可以肯定的是，这些选择也会影响结果。
蒙太奇选择
蒙太奇或参考选择的争议问题由来已久。尽管参考影响信号特性这一基本概念已广为人知，但人们并不总是能认识到其后果有多深远。重要的是，所有蒙太奇都是不考虑时间动态的瞬时变换(即，在每个样本上独立执行计算)。尽管如此，蒙太奇的选择对信号的时间属性有相当大的影响。本研究考虑了三种常用的参考方案如何影响睡眠EEG的各个方面。
蒙太奇影响信号振幅和极性
图2显示了N3信号在四个中线通道上的30s记录，采用乳突(LM；A)、共同平均(CA；B)和表面拉普拉斯(SL；C)蒙太奇。简而言之，LM可以被认为是睡眠EEG的“默认”蒙太奇，参考位置位于相对不活跃的部位。对于CA，每个电极都参考所有电极的平均值，假设所有电位之和为零。最后，SL蒙太奇突出了每个记录位点的局部活动，同时减少了共同参考和容积传导的影响。
图2.不同蒙太奇的EEG比较。
如图2A所示，LM参考显示了大振幅SOs的典型额部表达，其向后通道的表达变弱。同时，跨通道的信号形状和极性相对相似(第一条垂直线)。相反，虽然CA参考产生的额叶SO看起来与LM参考相似，但Cz和Pz之间会发生极性反转，导致前后SO(第二和第三垂直线)之间存在反相关系。此外，与LM参考相比，后部SOs的振幅要大得多，接近额叶信号的幅值，而Cz SO则小得多。这是由CA计算得出的：因为根据定义，在每个样本上，所有电极的活动总和必须为零，所以前额区域的大负振幅必须与其他地方的大正振幅相匹配。与其对局部活动的关注一致，SL蒙太奇在通道之间产生了更多的变化，但也产生了相对于LM蒙太奇的轻微极性反转。重现图2的代码可以在referencePolarity.m中找到。值得注意的是，这些与蒙太奇相关的相对信号振幅变化对PSD地形图有影响，与功率归一化类似，详见referencePSDTopography.m。
蒙太奇影响功能连接
参考选择对功能连接也有很大的影响。这里将锁相值(PLV)作为相位同步的测量方法。相同频率脑区之间振荡相位的一致关系被认为能够使底层神经元群之间进行有效传播。图3A的地形图显示了SO(左)、慢纺锤(中)和快速纺锤活动(右)的种子通道和头皮其余部分之间LM蒙太奇的PLV，种子通道的选择与图4A的PSD地形图相对应。在每个频率上，与附近通道的连通性最大，并随着到种子的距离增加而降低。种子和头皮其余部分之间的平均相位差地形图表明几乎不存在相位差(0°)，与图2A中通道之间的相似性相匹配。与此形成鲜明对比的是，CA参考产生了从种子到附近和远处通道的强连通性(图3B)，其模式大致遵循图4B的PSD地形图。此外，与远距离区域的连通性在每个频率上都显示出明显的反相位关系，与前面描述的极性反转一致。最后，SL图显示连接的范围要小得多，相位差图噪声更大。重现图3的代码可以在referencePhaseSynchrony.m中找到。
图3.蒙太奇对相位同步的影响。
图4.蒙太奇和功率归一化对PSD的影响。
蒙太奇影响跨频率耦合
最后，参考的选择影响跨频相位-振幅耦合(PAC)，即较快振荡的幅值(或功率)取决于较慢振荡的相位。这被认为可以实现跨多个时空尺度的大脑交流。图5显示了使用去偏相位-振幅耦合(dPAC)方法计算得到的慢纺锤(左)和快速纺锤(右)最大表达的SO相位地形图。LM蒙太奇地形图具有非常均匀的外观(图5A)，慢纺锤在SO波谷约180°时达到峰值，快速纺锤在SO波峰0°附近达到峰值。该结果再次与容积传导一致，确保SO和纺锤波在不同电极上具有最大振幅和功率(图3A和图4A)。相比之下，CA蒙太奇产生大致的双峰图，使得前部和后部区域的纺锤体活动与相反的SO相关联，反映了该参考选择固有的极性反转(图5B)。最后，SL相位图出现了噪声(图5C)，因为拉普拉斯会产生SO或纺锤体活动的相对集中的热点，因此很难从所有通道中提取有意义的相位耦合估计(图4C)。实际上，SO-纺锤耦合的程度(而不是相位)也取决于蒙太奇。通常，应将耦合相位的解释限制在实际显示这种耦合的通道上。然而，原始耦合强度通常难以解释，这将在替代检验中进行讨论。重现图5的代码可以在referencePhaseAmplitudeCoupling.m中找到。
图5.蒙太奇对慢波振荡(SO)-纺锤耦合相位的影响。
瞬时相位和振幅
准确提取瞬时相位和振幅信息对于回答与功能连接和跨频相互作用相关的基础和临床睡眠研究问题至关重要。例如，先前的睡眠研究发现精神分裂症患者的纺锤体相干性减弱，创伤后应激障碍患者的α带相位同步改变，以及与记忆力下降相关的SO-纺锤体PAC精度较低。根据所使用的指标，需要提取特定频段或一系列频段的带限相位或振幅信息。对于许多指标，这些计算可以在时域或频域中等效地执行。频域计算依赖于信号的全PSD和交叉频谱密度，而时域方法需要提取逐采样(或瞬时)相位和振幅信息。获得这些瞬时估计的两种流行方法是：①在感兴趣的频段内对EEG进行带通滤波并应用希尔伯特变换，②将EEG与一组不同频率的复值小波进行卷积。这两种技术都提供复值输出，从中获得对相位和振幅的估计。但是，在许多方面，返回的估计可能是不准确的。例如，在评估SO-纺锤耦合相位时，或者当振幅估计值被提供给下游检测算法时，精确的相位/幅度估计是至关重要的。接下来将讨论与相位/振幅估计相关的几个问题。
提取相位
相位表示波在周期中的相对位置，用弧度或角度表示。相位角通常是相对于余弦或正弦波来定义的。当从(正确应用的)希尔伯特滤波器或小波卷积技术的复值结果中提取相位时，它是根据余弦隐式定义的，峰值对应于0弧度(或 0°)，正零交叉对应于−0.5π/1.5π弧度(或−90/270°)(图6A：i和ii)。如果需要正弦相对相位角，则必须手动添加0.5π弧度(图6A：iii)。然而，从文献中并不总是能看出采用了哪种约定。当转置运算符(')应用于复值时，Matlab中会出现另一个问题：当该运算符正式执行复共轭转置时，虚值被符号翻转，相位沿相反方向前进，导致相位估计错误(图：6A：iv)。相反，应该使用非共轭转置(.')。最后，重要的是要使用专用的圆弧工具计算角度变量的统计属性(例如，均值，标准差，相关性)。重现图6A的代码可以在extractPhase.m中找到。
图6.提取相位和振幅。
提取希尔伯特振幅
图6B显示了模拟的SO-纺锤相位-振幅耦合，纺锤振幅在SO峰值处最大(红色迹线)。一旦在SO(图6B：i)和纺锤(图6B：ii)范围内对数据进行滤波，并应用希尔伯特变换，就可以得到滤波后的信号及其振幅包络分别作为希尔伯特变换后数据的实部(黑色)和振幅(绿色)。然而，滤波必然会引入频谱和时间的不精确性，这从实际包含这种活动的窗口外的明显的纺锤体活动(橙色)和相对于其指定振幅1的纺锤体包络线有所减少就可以看出。重现图6B的代码可以在extractHilbertAmplitude.m中找到。
提取小波振幅和相位
虽然从小波中提取振幅/包络信息在概念上与希尔伯特方法相似，但小波需要在频域中重新调节以获得与原始数据一致的值(图6C：i和ii)。此外，当在频域中执行小波卷积时，如果时域小波通常以最大值为中心，则相位估计将不正确(相移)(图6C：iii)。准确的偏移量取决于小波特性，包括频率和小波长度。相反，小波应该被包裹起来，使其两个时域的一半基本上是交换的(图6C：i；插图)。详细解释这些问题的代码可以在extractWaveletPhaseAmplitude.m中找到。这些问题是否会影响分析取决于在下游分析中如何准确地使用相位/振幅估计。
创建小波
尽管存在多种类型的小波，但这里只考虑(复)Morlet小波，它是由高斯窗口逐渐变细的(复)正弦波。除了频率之外，它们还由许多参数(长度、数量、频率间隔、时间分辨率、频率分辨率、时频)定义，这些参数通常需要仔细考虑，而睡眠EEG则需要考虑不同的因素。特别是，应该选择与感兴趣的振荡成分相匹配的小波频率，这可能需要非常慢的小波来捕获SO，并需要适当间隔的小波具有相对较高的频谱精度来识别慢速和快速纺锤。此外，相对较长(例如，20s)的小波可确保即使在最慢的频率下也能在频域中保持高斯形状。因此，各种工具箱中使用的默认值通常不是针对睡眠EEG进行优化的。灵活构建和可视化Morlet小波的代码可以在createWavelets.m中找到。本研究还包括用于创建小波的代码，这些代码具有适合睡眠头皮EEG的属性(createWavelets_sleepScalpEEG.m)。
替代检验
在计算了功能连接或跨频率耦合的指标后，通常会在通道、频率、组之间比较原始值。但是，原始值并不总是易于解释，因为在没有连接/耦合的情况下不知道预期值是什么。虽然对于许多指标来说，理论值为0表示没有连接/耦合，但几乎没有观察到确切的0值，而且不清楚值应该距离0(0.1，0.01，0.001？)多远才能表示有意义的连接/耦合。此外，在没有真正影响的情况下，期望值通常取决于功率、频率、通道和个体等因素，如果要进行公平的比较，就应该考虑到这些因素。这些问题通常可以通过替代检验技术来解决，这种技术特别适用于基础和转化睡眠EEG研究，因为它们通常有大量的数据。
替代结构
测试连接/耦合是否存在的常用方法是构建大量替代项：原始数据的修改版本，保持某些属性不变。对每个替代项重新计算感兴趣的指标，并将实际数据中的观测值与替代项值的分布进行比较。当观测值在替代分布(也称为零分布，因为它反映了无效的零假设)的情况下完全不可能时，可以推断连接/耦合。然而，替代检验需要考虑许多因素，其中一些影响很大。
有不同的方法可以将观测值与替代分布进行对比。首先，可以简单地评估替代值大于或等于观测值的比例，提供一个非参数的单侧p值。在图7A中，在SO相位和慢速纺锤振幅之间的额通道AFz处观察到的PAC大于1000个替代值中的每一个(P＜0.001)。其次，如果替代值近似正态分布，则可以用分布的均值和标准差对观测值进行z评分(图7A，右)，此处为23.0。由于每个z分数都有一个关联的p值，因此该方法提供了另一种确定显著性的方法(参数)。此外，由于z评分表示某个值与均值的距离(就标准差而言)，因此当替代分布在这些方面发生变化时，它提供了更容易比较的标准化分数。当考虑后部通道Oz时，这一点很明显(图7B)。虽然与AFz相比，原始PAC大幅降低(0.05)，但原始替代分布也发生了变化，仍然产生了1.9的中等正z分数。
图7.替代选择对不同通道SO-慢纺锤耦合的影响。
替代影响耦合动力学
图8A显示了SO-纺锤相位-振幅耦合强度(dPAC指标)的地形图。虽然原始和归一化耦合强度在慢速纺锤体上显示出相似的额叶地形，但快速纺锤体在顶枕(原始)或额叶(归一化)区域上似乎受到SO相位的调制最为强烈。原始耦合强度的通道差异可能部分反映了SO/纺锤密度/振幅之间的差异。相反，替代归一化dPAC表示在纺锤包络分布和SO相位存在情况下的耦合程度。虽然两者都是相位振幅耦合的有效表征，但它们可能导致耦合最强出现在何处的观点存在本质上的不同。
图8.基于替代的归一化对耦合动力学的影响。
图8B显示了考虑原始或归一化加权相位滞后指数值(wPLI)时的远距离AFz-Oz跨频率相位同步。原始wPLI(蓝线)在SO波段显示了最高的连接值(~0.55)，而在纺锤范围内的峰值为(~0.4)，这可能导致SO比纺锤更同步的结论。相比之下，时移替代归一化wPLIZ(红线)表示纺锤范围内的连接性最强(~1.75)，而SO则小得多(~0.75)。在许多情况下，比较原始指标是合适的(例如，比较具有相似频谱功率的被试内条件之间的连通性)，但也需要仔细考虑，以确保所采用的指标与特定的研究目标相匹配。
结论
本文回顾了人类睡眠EEG分析中的各种方法学问题，特别是与频谱分析、蒙太奇选择、相位/振幅信息的提取、替代结构有关的问题。首先，通常被认为在概念上不感兴趣的因素可能会对结果和结论产生巨大影响。因此，报告方法需要足够详细，以便其他研究者可以进行再分析。其次，小的计算细节可能导致结果指标提供对感兴趣因素的无意义的估计。所以可视化中间处理步骤对于确保代码按预期运行是至关重要的。最后，将睡眠EEG与认知和疾病联系起来很重要，但也很困难。虽然已经取得了一定的进展，但如果想要真正了解睡眠的作用，则可以而且应该提高方法学质量。总之，本研究希望所提出的问题有助于未来研究人员更好地进行睡眠脑电研究。
原文：Analyzing human sleep EEG: A methodological primer with code Implementation.
https://doi.org/10.1016/j.smrv.2020.101353