ERP实验设计编程和数据采集

  • ERP实验设计和数据采集与行为实验存在较大的差别，具有一些特别的要求，如实验刺激的试次和类型。如果这两个步骤出现问题，则后期的数据处理就失去了大楼的基石。
  • 在哔哩哔哩上看到一个up主的一句话，甚为赞同。安全的食品是生产出来的，不是检测出来的。科学实验数据也一样。干净的实验是做出来的，而不是处理出来的。这一点在ERP实验和fMRI实验上体现最深。前期的数据采集工作结束的时候，就已经决定了这批数据的命运。实验设计的不严谨，无关实验因素控制不严格，都可能导致得到的是一堆垃圾数据。面对垃圾数据，就算是再好的数据处理方法也不可能变废为宝！因此，在设计ERP实验和fMRI实验时，从材料选取，呈现流程，Mark标记到数据采集，要事无巨细，方方面面都要考虑清楚。任何一个方面的疏忽，都会导致百分百的实验失败。

Experimental Design

  1. 实验设计的第一步是根据实验目的选取可靠的实验范式。实验范式选取的第一选择是依据前人文献选择经典的，被反复验证的实验范式。如果前人范式无法满足实验目的，再在前人范式的基础上进行创新改进。另外，需要注意行为范式不能直接应用到ERP或fmri实验中，后者都有其各自的技术特点。

  2. Ask yourself the following questions.

    1. Which cognitive processes do you want to study? Why do you want to use ERP approach to study it?

    2. Which ERP components you want to use to represent the cognitive processes of interest?

    3. Is your design in line with the ERP princinples (see Luck book for more detail)?

    4. How many trials? How many participants?

Data collection

  1. 准备哪些材料？ 脑电膏，棉棒，胶布，去角质膏，脑电帽和配套识别文件
  2. 脑电实验特别注意事项：头动（超级影响数据），眼睛和屏幕间距，刺激呈现在屏幕中的位置和大小，房间电器的控制（电磁干扰）
  3. 打电极膏注意事项，保证电阻达到理想较低水平
  4. 数据采集过程中实时观察数据质量并记录

ERP Components in the field of emotion regulation

  1. 每个成分需要注意的事项：
    1. 基线的标准（-200ms）。
    2. 其参考点有没有特别要求，平均参考还是双侧乳突？
    3. 其空间分布的经典位置在哪些电极点上？
    4. 其时间分布的经典时间窗是哪些时间段？
    5. 其时间分段的通用标准是什么？
    6. 其常用的指标又哪些，振幅还是潜伏期？每个指标的心理学含义是什么？
    7. 其操作化定义，具体计算的方法是什么？
  2. 建议：不要仅仅靠语句计算出ERP指标然后就去做统计。再这之前，肉眼观察每个测量到的 ERP成分。画出每个被试的波形图，观察比较其成分的波幅和潜伏期。这是因为被试间的波形 往往差异很大。

Early components (<300ms)

  1. C1, P100, P200, N100, N170（面孔特异性）, N200,

  2. Early Posterior Negativity (EPN), highly responsive to phylogenetically threatening stimuli

  3. MisMatched Negativity (MMN)

  4. Error-related negativity (ERN) typically peaks from 80-150 milliseconds (ms) after the erroneous response begins

  5. 关注哪一个早期ERP成分取决于实验刺激呈现的感官通道。不同的感官通道，如视觉和听觉，所诱发的 ERP成分（如N1），其到达顶峰的时间可能有所重合，但是头皮分布往往差别很大。 例如，与任务相关的听觉和视觉刺激都会诱发出N2b成分，但听觉刺激的反应在中央区最大， 而视觉刺激的反应在头后部最大。

Late Components (>300ms)

  1. P300 LPP N400（语言加工成分）
  2. LPP： early (300-1000) and late (>1000) LPP
  3. 必须以双侧乳突为参考点
  4. 建议以图片呈现时间来进行分段，例如图片呈现4s，就以4s来分段。
  5. LPP可以分为早期LPP（<1000ms），和晚期LPP
  6. 计算方法：建议从300ms开始，以200ms为最小时间单位分段。

ERP成分经典指标的计算

  1. ERP成分的分析是基于头皮电位的分析。所以首先需要确定电极点（electrodes）-根据参考文献和综述。

  2. 确定peak点。local peak amplitude method should be used, because it reflect the peak better than the maximum voltage method (See Luck 2014). 可以使用matlab中的findpeaks函数实现。

[pks,locs] = findpeaks(data) additionally returns the indices at which the peaks occur.

  1. 平均振幅。选用某些电极点peak前后某段时间窗内的电位平均值。 注意，这一方法不宜选用过窄的时间窗（<40ms）。这是因为虽然 更窄的时间窗能够更好地避免重叠成分的影响，但其受条件间、电极点之间以及不同被试间 变异的影响更大。

  2. 潜伏期（latency），即经典成分peak出现的时间。

ERP-EEG Data 整体流程

EEG数据处理的两大核心思想就是，降噪（排除干扰因素）和平均。

EEG数据处理百科

  1. CIMeC Wiki on EEG Data pre-processing
  2. EEGLAB Wiki

第一步，预处理流程

  • 预处理之前，需要对EEG的数据结构（被试试次通道时间空间*频域）有一个清晰的认识，搞清楚EEG数据有几个维度。
  • 预处理方法的选择和步骤取决于实验设计的细节，收集数据的设备，您计划执行的分析，以及您或您的研究小组开发的特殊协议和偏好。所以预处理由大概的步骤，但没有绝对的答案，需要结合自己的研究目的来决定。 例如，ERP分析和时频分析的预处理步骤，在重采样、滤波（ERP主要分析的是低频，0.01-1Hz，而TFA是1-50HZ）以及分段（TFA分段要长一点）上存在一些差异。
  • EEG数据的每一个维度上，都可以进行预处理以实现降噪。
  • 不同软件设计的预处理流程会存在一定的差异。现在总结常用的几种预处理和后期处理流程：
    1. EEGLAB代码进行预处理，包含分段（epoch）；使用EEGLAB或fieldtrip代码进行组分析；
    2. EEGLAB代码进行预处理，不包含分段（epoch）；使用ERPLAB进行分段（erplab的术语叫建立erplist，assiganbins）和ERP组分析
  • 第一个流程基本上是全部代码流程，最后的绘图和数据导出都需要代码实现；流程2ERPLAB自带绘图和数据导出功能，而且目前还有新的数据质量检测功能加入。

准备工作

  1. Prepare data1. 检查被试的命名是否统一（e.g., sub01_name.*）, 检查所有被试的marker命名是否正确；

  2. Prepare data2. 如果存在不一致的问题，例如由于编程时问题，marker名称不统一或重复， 使用语句对marker进行重命名。确保数据中，实验条件与marker一一对应，没有错误。

  3. 确定预处理参数：离线参考点。主要离线参考点如何选择（全脑平均还是双侧乳突，或者鼻尖参考）。分析目的和指标不一一，参考也不一样。
    • the target ERP component: linked mastoids are good for components with vertically oriented dipoles (e.g., P300), whereas average reference is better for lateralized components (e.g., P1/N1);
    • the planned analysis: if you plan to run topographical analysis or estimate neural sources, you have to use average reference. For more information, see: Issues in the application of the average reference: Review, critiques, and recommendations (Dien, 1998, BRMIC) For a short, informal discussion on the topic, see https://sccn.ucsd.edu/pipermail/eeglablist/2008/002510.html
  4. 确定预处理参数：高通滤波。
    • High-pass filtering above 1 Hz had dramatic effects on ERP shapes. 0.1Hz is common for ERP study Does filtering preclude us from studying ERP time-courses?
    • There is a published study showing that high-pass filter up to 2Hz improves ICA’s performance. See this paper by Winkler et al.(2015).
  5. 确定预处理参数：重采样 * Downsampling to something like 100 Hz is generally recommended for connectivity analyses.
    • TFA分析最好不要重采样

第一次预处理

  1. 第一次预处理目的：生成基于1HZ高通滤波数据的ICA矩阵信息。基于1-2Hz做ICA的论文
  2. EEG的数据的信息可以分为两大类别：时间信息和通道（空间）信息。注意，时间和空间信息的处理不是串行进行的 而是你中有我，我中有你。
  3. 想哥思路，想哥之前代码的思路是使用cleanline，clean_rawdata等对时间上的信息进行去噪； 空间上，使用set_channel_types函数把通道分为EOG,bad,eeg等几个大类： 第一大类非脑通道，包括VEOG,HEOH，ECG等； 第二大类为坏道，例如前额叶的点由于阻抗很大，一般信噪比很差； 第三大类就是正常的大脑通道，如Cz等。 想哥是在ICA之前，就把三大类的通道区分开来。ICA的时候只输入正常大脑的通道。
  4. 需要学习的是，想哥代码在处理时会把去除的坏段，坏导等，都记录下来（存到bk中）。这有利于数据的QA。例如，在预处理后对所有数据中的坏段数据进行汇总分析，可以计算出每个被试数据的留存率，以及组水平上有效数据的留存率。
  5. SCCN Makoto’s preprocessing pipeline Makoto在2020年的处理流程中，跟想哥的处理流程有几点改善的地方：
    1. Makoto强烈建议跑ICA的时候，如果EOG跟eeg通道用的同一参考，就使用这些通道。EEGLAB论坛里也有不少人支持 不要在ICA之前删除EOG通道，而是把EOG数据输入到ICA中，这样可以提高ICA区分眼动成分的能力。
    2. 去线性噪声的插件由cleanline换为Prep
    3. Makoto建议（可选）在ICA之前对数据进行分段，以便ICA的结果在时间维度上更契合ERP分段数据；否则ICA就是在连续数据上所得出。我想这里可以只分段，但不需要去除坏段。 因为ICA前的数据包括了EOG等绝对属于non-brain噪声，而去坏段的时候是把所有通道平均再去掉，所以EOG这些non-brain通道的存在会导致epoch的信噪比很差，很容易全部去掉。

  6. cleanRaw参数的设置需要非常注意。设置过于激进会删除大部分的数据。而且clean_rawdata()的帮助给的参数并不很好。 cleanRaw参数具体可以参考此教程： EEG = clean_rawdata(EEG, 1, -1, 0.8, -1, 8, 0.25);

    %The first box provides information on any flat channels %The second box is disabled (-1) because FIR high-pass filter is already applied. %The third box means that if a channel is correlated to the surrounding channels less than 0.8, the channels is rejected. If it rejects too many channels, use less value.

    %The fourth box is disabled (-1) because line noise is supposed to be already taken care of (with Cleanline)

    %The fifth box is the ASR （该选项不会删除数据，而是会用插值计算的方法修补不良数据）. This is the point. Finds calibration data (i.e. the cleanest part of data) %Apply 1-s sliding window principal component analysis (PCA) on the continuous data. Classify PCs into high variance (in this case, 8 standard deviation from the % calibration data) or normal variance. %Reconstruct high variance subspace (i.e. group of channels associated %by a PC) with normal variance subspace. An important assumption is that %neighboring channels are heavily correlated due to volume conductance and scalp mixing %because of this, the central limit theorem works and the scalp channel signals %become closer to Gaussian than source signals; this guarantees ICA. See also Arnaud Delorme’s %explanation of ICA as non-Gaussianity maximization). This is how you %can estimate a signal of one channel from signals of adjacent channels. % Makoto advises a value that is conservative (from 10-20), for the %purpose of correcting only absolutely bad glitches that critically influences ICA’s performance. %Nima’s result showed that ASR with SD=20 produced the best result. % The understanding is to let ICA explain the data as much % as possible so the correction using PCA-based method should be minimized, and the % empirical balancing point is somewhere around 10-20.

    %The sixth box is the window rejection criterion（default=0.25）. 如果第四个选项的ASR无法修补的数据，该选项设置后就会将其删除掉。

  7. 第一次预处理的时候是否要用clean_rawdata先去除坏导，再补回坏导，甚至说要不要还有clean_rawdata进行ASR的数据修补？ 我想至少不能去除坏导再补回，因为这样EOG等很可能被直接去掉了，补回来的数据就不再能反映眼动。失去了在ICA中加入EOG的意义。ASR感觉还是得做， 因为ICA是一种平稳信号降维（stationary signal decomposition）技术（ICA也能应用到非稳态数据上，只是效果好不好的问题）。不使用ASR修补数据的话，数据的不稳态性质很高。ASR应该可以提高数据的稳态性。
  8. ICA参数的计算。主要需要计算的输入参数就是数据的data rank，其取决于通道数目。务必保持两次预处理通道数目的一致。
  9. ICA方法的选择：

第二次预处理

  1. 应用上一次生成的ICA矩阵信息，并对数据进行分段，剔除坏段。
  2. 注意不能在ICA之前把EOG等通道删除，否则会因为ICA矩阵维度不一致报错。
  3. 两次预处理的重参考方法需要保持一致。

眼电剔除

  1. 手动（半自动）： 利用EEGLAB Adjust插件半自动化地去除眼电的ICA成分。
  2. 全自动：ICALabel

第二步，检查数据，保证质量

有些被试由于头动等莫名其妙的原因，其电位可能比正常被试的数据大几十倍。In this case, 这一个极端被试的混入就能左右结果的显著性。为了避免这种极端值，需要对数据进行质量控制， 剔除极端被试，避免污染结果。

  1. 对每个被试的感兴趣成分进行绘图，观察电位强度范围，以及ERP波的走向。 例如LPP的波幅强度（就情绪唤起而言），一般在10-15左右。如果一个被试的电位到了100， 这个被试很可能就会影响结果的显著性。而且这种被试的原始数据一般也会有问题，怎么剔除坏段 都救不回来。Rubbish in, rubbish out.
  2. 求出所有被试某个成分的值，直接观察原始值。然后重点对极端值绘图。
  3. 注意，第一步和第二步并不是严格的时间先后顺序，两者应该是交互进行的。

组分析一：时域-ERP成分分析

1. 时域分析简介

时域分析只考虑信号电压强度随时间变化的规律，不考虑频域信息。时域分析大体可以分为以下几个类别。

  1. 常规ERP分析。比较不同实验条件下的特定ERP成分的强度。
    • 锁时（time-locked）指锁时一般是按照刺激的onset叠加平均。ERP成分的锁时特性是试次平均之后能够得到ERP的根本前提，因为如果大脑对刺激不出现在同一个时间段内，平均就没有意义了。量化锁时可以通过分析时域的ERP 波形或者时频域内的ERSP(event-related spectral perturbation)来获得。ERSP能够衡量刺激相关的能量相比基线的变化。具体而言，实验能量（谱分布）的叠加之后得到ERSP图，如果ERSP相对于基线有所提高或降低，就说明能量对于时间锁定得很好；而如果ERSP没有随着时间变化则说明能量对于时间没有锁时关系。如：ERP和ERS相对于基线能量都有所提高，而ERD相对于基线能量有所降低，所以ERP ERD/ERS都是time-locked的。
    • 锁相(phase-locked)。相位是反映波形的走势，是跟能否获得好的ERP相关的。如果得到好的ERP波形，那么可以说多次的实验中，大脑对刺激的反应的相位也是很一致的，亦即，波形的走势是很一致的，也就是说相位很好地锁定了，同时可以说不同的试验数据相关性很强。量化锁相可以通过分析时频域内的ITC(inter trial coherence)来获得。一般来说，ITC越接近于1，越说明不同的实验，大脑对刺激的反应的波形走势越一致，相位就锁定得越好，那么就很可能通过平均获得很好的时域内的ERP波；相反，如果接近于0，说明相位没有锁定关系not phase-locked，则很难获得ERP了。
    • 总之，叠加平均能出得来的东西，必须既time-locked又phase-locked，如ERP。因此锁时锁相是能够得到良好ERP波形的前提。但是，锁时锁相并不是和ERP以及频谱分析相对应的，不锁时锁相也能得到结果，不过可能是很烂的结果。
    • 参考文献： David et al., 2006 NeuroImage 31 1580 – 1591；Basar et al. (2001)和Yener et al. (2007)讲述了一些关于锁相的问题。
  2. 时间信号相关分析，认知神经科学一般称之为功能连接（functional connectivity）。计算功能连接有两大基本要素，第一选择ROI并提取该ROI内的时间序列信号；第二，计算ROI之间的相关系数。
    • ROI的选取。如何通过不独立的EEG信号得到独立的大脑皮层ROI信号？
    • 相关方法的选择。直接使用相关方法存在许多问题。同步似然指数，信号能量的包络（Hipp et al., 2012 Nature Neuroscience: Large-scale cortical correlation structure of spontaneous oscillatory activity）。参考雷旭老师博客。
    • 推荐HERMES工具包。HERMES是由西班牙马德里技术大学（Technical University of Madrid）的Centre for Biomedical Technology团队研发的基于Matlab的开源EEG工具包，其主要的功能和特点是计算基于各种方法的功能连接，HERMES官方网址：http://hermes.ctb.upm.es/

2. 时域分析操作：使用ERPLAB进行组分析（时域）

  • 使用ERPLAB进行组分析的优点在于ERPLAB自带较好的后期绘图和数据导出代码
  • 使用ERPLAB进行组分析，需要注意的是前面预处理的时候不可以进行分段。ERPLAB只识别连续数据，分段后无法识别。
  • ERPLAB处理基本流程：
    1. 单个被试：ERPLAB——CreatEventlist——Assignbins——Extract bin-bassed epochs（相当于EEGLAB的分段）-Artifact detection-average。单个被试数据分析完之后，就可以输出代码，添加循环改成批处理即可。
    2. 组水平平均（平均被试维度）：载入ERPsets或者选择本地ERPsets文件-总平均（Average across ERPsets (Grand Average)）-生成新的Grand Average ERPset。后面就是在总平均数据上进行绘图和指标导出工作。
  • erplab的bin相当于事件名称（如负性和中性图片）加数字化mark（如101，102）的组合。
  • 中文知乎教程
  • 英文ERPLAB官网manual
  • 练习数据和代码集ERP CORE

组分析二：时频分析

1. 时频分析简介

  1. 时频分析(Time-Frequency Analysis, TFA)同时考虑时域和频域信息，对特定时间段和特定频率范围内的信号能量强度进行分析。
  2. 常见术语：
    • Time-Frequency Analysis, TFA
    • Event-related spectral perturbations, ERSP
  3. In signal processing, Time–Frequency Analysis (TFA) comprises those techniques that study a signal in both the Time and Frequency Domains (TFD) simultaneously.
  4. 脑电波中α频段（8~12Hz）和β频段（13~30Hz）的振荡(降低和增强)可以显示出事件相关的中止反应或事件相关的增强, 其中前者被称为事件相关同步化(event-related synchronization, ERS), 后者称之为事件相关去同步化(event-related desynchronization, ERD).
  5. ERS和ERD具有以下特性(参考自刘刚&孙伟, 2011):
    1. 频段特异性: alpha和beta频段的ERD可作为大脑激活的指标， 而ERS的出现则表明大脑处于失活状态，提示皮质区域恢复到静息或惰性状态(Pineda 2005) 。 与之相反，gamma频段的ERS 与皮质激活相关，可能参与多区域和多种模式的信息整合加工过程(Cheyne, et al., 2008) 。 不同频段在正常和病理的运动加工过程中表现出不同的振荡模式，表明在大脑皮质和皮质下水平可能有着相互作用的独立功能调节机制来参与运动神经环路的信息加工。
  6. 锁时性: 前人EEG和MEG的研究发现感觉运动皮层的神经反应振荡活动在运动开始前几百毫秒出现并持续到运动结束后几秒钟。 例如，在运动开始前几百毫秒，alpha和beta频段的能量会迅速降低（ERD），在运动结束之后１ ～２ ｓ beta频段能量会迅速增高（ERS）(Pineda 2005)
  7. 情绪相关研究中, ERS和ERD往往用于确定感兴趣的时间-频率窗(TFD of interest)
  8. 常见因变量指标：
    1. 能量/功率谱密度（power spectral density, PSD）. 确定好感兴趣的TFD之后, 就需要从感兴趣的TFD中提取每个实验条件的节律能量.
    2. 相位一致性(Phase Locking Index, PLI)
    3. 耦合(Coherence). 不同电极同频率信号的相位同步性(Phase Synchrony)

2 时域到频域的常见方法

  1. 短时傅立叶(short-time Fourier transform,STFT)变换，可以对感兴趣的较低的频率进行分析。STFT有诸多优点, 但是也有两大缺点:
    1. First, it treats the power within the window as if it was the power at the center of the window, even though the entire window contributes equally to the power measurement.
    2. STFT依赖于一个固定的时间窗. 而时间窗的长短直接影响时间分辨率和频率分辨率. 窄窗口时间分辨率高但频率分辨率低, 宽窗口时间分辨率低但频率分辨率高. 对于时变的非稳态信号(EEG信号就是典型), 高频适合用小窗口, 低频适合用大窗口. 然而STFT的时间窗是固定的,在一次TFA分析中是不会变化的. 因此, STFT无法满足非稳态信号变化频率的要求.
  2. 小波变换。小波有很多种，脑电TFA主要使用的是Morlet wavelets。小波变换并没有采用窗的思想，更没有做傅里叶变换。小波直接把傅里叶变换的基给换了——将无限长的三角函数基换成了有限长的会衰减的小波基。这样不仅能够获取频率，还可以定位到时间了。 用高斯函数(Gaussian function, bell curve) function)乘以正弦波就成了Morlet小波。Morlet小波也有很多种。但是，小波分析中使用到的小波函数多种多样。小波分析在工程应用中,一个十分重要的问题就是最优小波基的选择问题在使用小波变换时，小波基函数的选择非常关键（丛丰裕老师讲座）。There are several advantages of Morlet wavelets for time-frequency analysis:
    • One is that the Morlet wavelet is Gaussian-shaped in the frequency domain (Figure 2, bottom row). The absence of sharp edges minimizes ripple effects that can be misinterpreted as oscillations (this is a potential danger associated with plateau-shaped filters).
    • Second, the results of Morlet wavelet convolution retain the temporal resolution of the original signal.
    • Third, wavelet convolution is more computationally efficient and requires less code compared to other methods, because it involves the smallest number of computations, most of which are implemented using the fast Fourier transform (Cohen, 2019).
  3. 小波变换有两个变量：
    1. 尺度a（scale），尺度a控制小波函数的伸缩，尺度就对应于频率（反比）
  4. 平移量 τ（translation）。平移量 τ控制小波函数的平移，平移量 τ就对应于时间。
  5. 那么怎么选择合适的小波基函数？光靠想是不行的，需要结合实际数据进行绘图检查：

3. 分析操作

  • 时频分析预处理的注意事项：
    • 使用小波变换的话，分段要长一点，以避免边际效应。例如你对0-1000ms感兴趣，则TFA的分段可以在-1500-2500ms，前后各加1500ms。所加时长的长短，取决于研究者感兴趣的频段。感兴趣的频段越低，则加的时长越长。 例如对于0.5Hz，则可能要加好几秒；而对于100Hz，则十几ms就足够了。如果不确定的话，可以用一个被试的数据反复试几次找到最优解（Cohen, 2014，pp76-77）。一般而言，低频时，边际效应会持续三个周期（cycles）。以2Hz来计算，一个周期时长0.5s，三个周期就是1.5s，也就是1500ms。但这也要注意试次的时间间隔一定要大于分段时长，否则会造成数据重叠。虽然这些重叠的数据在绘图和提取指标时会删掉，对于数据分析没有影响，但会影响ICA分析。
    • 注意尽量保证不同条件之间的试次数目差异不显著。对于时频分析，原始power值总是正值。因此，额外的噪声会对power信号产生积极偏差。
    • 滤波，高通 0.1 or 0.5 Hz。注意高通必须在连续EEG（没有分段前）上进行。因为0.5Hz的低频信号会持续6s，而一般EEG实验的epoch达不到6s。TFA分析的低频滤波一般没有必要。
    • 重参考
  • 时频分析牢记的原则：: 每个时间点的活动是对瞬时活动的一个估计值，且其受到左右相邻时间点的影响。
  • 组分析有多种选择，常见的如: 1. 使用EEGLAB自带的功能进行时频分析. 但是EEGLAB的组分析study功能特别难用. 优点是只要study建好了, 后续也可以直接点点点. 2. letswave教程：脑电数据的时频分析/组平均与统计分析 letswave的优点在于不用编写代码, GUI全程点点点就可以得到结果. 3. Fieldtrip 4. 目前Alexander等开发了“Time Frequency Analysis Plugin for EEGLAB”（以下简称：TFA）。TFA既支持EEGLAB数据，也支持ERPLAB数据

操作方案1：使用EEGLAB TFA插件

第一步 构建小波函数
    1. 原理：先根据实验设计，构建合适的小波函数参数（小波基函数），以在时间和频域分辨率之间取得一个最佳平滑。
      • The crucial parameter of Morlet wavelets is the width of the Gaussian that tapers the sine wave （Cohen, 2019）. This width parameter controls the trade-off between temporal precision and frequency precision.
      • 小波循的次数（number of cycles）指从左侧0起到右侧归于0。循环次数这一指标本身没有对错之分。其问题在于不能清晰地与时频分辨率对应起来。
      • FWHM基本可以与循环次数对等，区别在于FWHM这个表达更加清晰，能够更直观地反映时频分辨率。This parameter is typically defined as the “number of cycles,” but the purpose of this paper is to argue that it would be better to define the Gaussian width as the full-width at half-maximum (FWHM), which is the distance in time (or frequency) between 50% gain before the peak to 50% gain after the peak （Cohen, 2019）.
      • 无论是时间FWHMt，还是频率FWHMf，都应该是值越小，对应的分辨率越高。这就类似于CPU的制作工艺，即纳米级越小，加工精度越高。
    1. 操作。TFA界面中 MWD Wavelet Design tool对应pop_WavPlot()函数。pop_WavPlot()来确定小波函数的核心参数。其理念与Cohen (2019)提出的用FWHM来表达基本是一样的。
      • 输入：（1）正弦函数的频率（如10或20Hz）；（2）number of cycles（如5或9）。
      • 输出：不同cycle下的FWHMt和FWHMf，以及图形。注意，这里的FWHMt和FWHMf后面文章中应该进行报告（Cohen, 2019）
      • 决定用哪个
第二步 分解EEG数据（时频变换）
    1. 原理：Morlet wavelet decomposition (MWD)，即基于小波变换的信号分解（时域到频域）。相当于把每个被试的EEG四维数据（条件-epoch-时间点-通道），变换为五维数据（条件-epoch-时间点-通道-频率）。
    1. 操作：函数：pop_MWD() 。
      • 输入：EEG/ERP数据，通道、频率、以及第一步计算出的小波函数的参数，以及输出类型等
      • 输出：OUTSET = 输入EEG/ERP + EEG.etc.MWD field that contains information abou tht edecomposition as well as the results of the decomposition
    1. 去基线

第三步 使用ERPLAB合并ERD/ERS（去掉trial/epoch维度）, 并进行绘图

  1. 使用ERPLAB对数据进行合并
  2. 根据提取的数据进行绘图
  3. 个体分析完成，后续用ERPLAB合并所有被试后，再在合并数据集上进行绘图操作、数据提取

操作方案2：使用fieldtrip进行组分析（时频）

  • EEGLAB的组分析代码，或者ERPLAB插件可以进行较好的时域分析，但时频分析不友好。所以时频分析使用fieldtrip来做。
  • 注意fieldtrip的预处理流程包含分段。
  • 一般流程：
    1. EEGLAB代码进行预处理
    2. 将EEGLAB预处理后的数据转换为fieldtrip可以识别的格式
    3. 使用fieldtrip进行时频分析

操作方案3：使用EEGLAB代码做时频分析

  1. 总体步骤:
    • eeglab里面前面的preprocessing与ERP分析大体相似，但细节不一。例如分段。
    • TFA组分析
  2. 函数. 利用每个被试的预处理后的文件，可以进行时间-频率分析。EEGLAB主要采用自带的newtimef()函数进行时频分析. newtimef()函数中进行TFA分析的核心函数有如下选择:
    • 第一, 短时傅立叶(short-time Fourier transform,STFT)变换
    • 第二, 采用小波变换,主要使用的是Morlet wavelets。

操作方案4：使用Evoked ERP-ERO Toolbox汇总数据和时频统计分析

数据汇总

  • Evoked ERP-ERO Toolbox该工具包只进行时域，频域或时频的组分析，并不提供预处理功能。因而可以很好地替代EEGLAB的Study功能。

  • 首先需要将预处理后的数据中的试次进行平均，然后按照条件拆分为每个被试每个条件一个数据。
  • 然后，使用Evoked ERP-ERO 提供的汇总工具Evoked ERP ERO_v1.1_Forming_fourth_order_tensor_demo data来合并数据，得到四维的汇总数据。就是一个代码，无需修改直接运行出现GUI界面，填写所需参数，选择所需数据，就可以傻瓜导出后续分析所需的四维数据。
  • 最终得到的用于ERP分析的数据是一个四维数据，四个维度分别为通道（chanells），时间点（time points），被试（subjects）和条件（conditions）。
  • Evoked ERP-ERO 工具箱提供的合并数据代码会生成两个数据：
    1. 一个是同时包含实验设计信息和数据信息的‘…_structure.mat’文件，可以直接在>File>Import the processed data’ 中导入Evoked ERP-ERO 工具箱，无需二次输入实验设计信息；
    2. 二是只包含被试数据信息的四维数据，‘test_fourdiemension_fourth_order_tensor’。需要在‘>File> Import fourth-order tensor & parameters’导入。导入时，需要再次填写实验设计信息。
  • 汇总原理：平均试次维度，计算出所有被试，所有条件，所有时间点，所有通道的平均值。 最后的数据就是一个四维数据，在Matlab中为mat格式文件。时频分析的数据为五维数据，相比ERP汇总数据多了一个频率维度。但是注意Evoked ERP-ERO Toolbox不需要先进行傅立叶变换得到五维数据再进行分析，其可以直接在四维数据的基础上进行时频变换和时频分析。
Time domain analysis (ERPs分析)
  1. 输入ERP_ERO打开数据，导入之前计算好的四维数据，默认为Dataset 1, 可以添加备注便于识别,如Dataset Raw.
  2. 可以先在>Time domain analysis> Conventional time domain analysis > Plot waveform & topo绘制波形图和地形图查看下结果。
  3. 然后，在Matlab命令行窗口按任意键，出现set the parameters for plotting top窗口，输入感兴趣的时间段以及指标类型（均值还是peak）。确认后，绘制出地形图和被试间相似性矩阵。之后，会在最初数据上生成新的数据集，如Dataset_Cz。
  4. 最后，可以选择在ERP_ERO内>Statistical analysis> Within-subject rm-ANOVA直接进行统计分析；也可以选择Outputdata导出数据到excel，再使用SPSS进行分析
  5. 如果想要换电极点重新分析，需要在Dataset里回到最初导入的数据集开始分析。
  6. 此外，还可以先对Dataset Raw在时域上做一个PCA分析，即主成份分析。Time domain analysis>t-PCA>Run t-PCA. 此处成分数目的选取根据The explanation of varience界面中第四个成分贡献百分比图来确定。例如图中说22个成分贡献了99%，那么就可以确定22个成分即可。后面利用被试间的相似性以及该成分的贡献度来确定分析成分。确定成分后，在t-PCA选择该成分，绘图。最后可以在该成分上进行统计分析。
Time-Frequency Analysis (TFA)分析
  1. 使用ERP_ERO载入数据后，>Time-frequency analysis > Morlet CWT> Run Morlet CWT，即可以使用连续小波变换（Morlet CWT）对数据进行时频变换，得到每个被试每个条件下频域的信息。
  2. 数据结构解读： 与ERP不同，TFA数据分析还需要在预处理后再进行时频变换，最后得到汇总的一个五维数据。与ERP汇总数据相比，多了一个频率维度。因此，TFA分析中ERP_ERO的数据结构EEG_ERO.data有五个维度，而ERP分析中该数据只有四个维度。至此，个人也可以完全使用代码提取出五维数据，根据自己感兴趣的时间-频域进行后续分析
  3. 在对数据进行时频变换之后，就需要确定自己感兴趣的时间（起止时间）-频域（起止频域），类似于fMRI分析中的确定ROI。ERP_ERO提供了两种方法来确定ROI：
    1. 传统方法：结合前人文献和当前数据结果，确定一个长方形来圈定ROI。
    2. 新方法，利用边界检测方法确定感兴趣震荡区域。
  4. 最后，可以对不同实验条件在ROI内的能量值（power）进行统计比S较，或者导出每个条件下每个被试ROI内的能量值进行后续分析。

第五步，绘图

  • 绘图首先要了解每种分析需要哪种图来表达，每一类图形的含义
  • 绘图需要两类软件：
    • 绘图软件
    • 排版软件，即把多个单独条件的图片整合到一起，排版并添加注释。

5.1 ERP分析绘图

  1. 了解常见的图形类型，含义以及绘图要求：
    1. 地形图
    2. 波形图
  2. 方案一：ERPLAB绘图，然后用canvas或者PPT排版。
    • ERPLAB绘图之前，先用ERP Operations求出自己绘图所需的数据集，例如A和B条件的差值，CP1和CP2两个通道的均值等等
    • 然后再进行绘图
  3. 方案二：Matlab输出到R，然后用R（ggplot2）绘图和排版

5.2 TFA绘图

  1. TFA分析主要涉及两个图:
    • 第一个是总平均时间-频率图. 该图是指计算每个被试每种条件下的时频图，然后进行平均得出的。基于总平均时间-频率图, 可以通过肉眼得到感兴趣的TFD. 但是注意，这一方法只能作为先验方法的辅助验证，不能作为选择TFD的主要依据。否则就会有循环论证的问题。
    • 第二个是感兴趣的频率的时间能量波形图。该图是根据自己感兴趣的TFD，只计算和分析自己感兴趣的频带的能量分布图。
  2. 数据. 要画出上述两个图，首先必须建立总平均时间频率文件, 即被试条件channeltimefrequency的数据. 这个文件是进行后续数据处理的基础. 基于这一数据, 后续分析如下:
    • 第一, 平均掉总数据集中的被试,通道和条件维度, 获得time*frequency的二维数据. 基于这一二维数据, 即可绘制总平均时间-频率图.
    • 第二, 根据总平均时间-频率图得到感兴趣的TFD之后, 即确定了所需要分析的具体波段和具体时间窗口.
    • 第三, 根据所确定的TFD(已确定的时间和频率), 从总数据集合中提取出每个被试在每个条件下该TFD的ERSP的值，然后根据该数据进行进一步的统计分析(如ANOVA)。
  3. 基线矫正(Baseline Correction). 主要基线矫正方法有如下三种:
    • some studies subtract the average baseline power at a given frequency from the power at each poststimulus time point for that frequency.
    • Other studies use division rather than subtraction (i.e., for each frequency, the power at each poststimulus time point is divided by the average prestimulus power).
    • In other studies, the change between the prestimulus baseline and the poststimulus period is represented on a log scale (decibels) to take into account the fact that power typically falls off as the frequency increases.

Notes

  1. what are ERPs and what are they good for？
  2. ERP研究中如何设置刺激呈现时间，ITI、SOA？
  3. EEG来自于哪里？有什么意义？
  4. [技术帖 如何在任何一项ERP实验中发现显著的效应](https://mp.weixin.qq.com/s/ifyxia0BkC4niGPCjaOw7w)
  5. ERP实验设计逻辑与实验范式 6.

References

  1. Uusberg, A., Thiruchselvam, R., & Gross, J. J. (2014). Using distraction to regulate emotion: Insights from EEG theta dynamics. International Journal of Psychophysiology, 91(3), 254-260.
  2. https://www.cnblogs.com/minks/p/5946442.html
  3. 知乎: 如何通俗地讲解傅立叶分析和小波分析间的关系?
  4. 武侠, 钟楚鹏, 丁玉珑, & 曲折. (2018). 利用时频分析研究非相位锁定脑电活动. 心理科学进展, v.26；No.216(08), 23-38.
  5. https://neuroimage.usc.edu/brainstorm/Tutorials/TimeFrequency
  6. 丛丰裕老师主页： http://www.escience.cn/people/cong/index.html
  7. 悦影科技 https://zhuanlan.zhihu.com/p/97394678
  8. 脑电信号功能连接 http://blog.sina.com.cn/s/blog_60a7516201015m7q.html
  9. EEGLAB batch1 EEGLAB batch2
  10. Cohen, M. X. (2019). A better way to define and describe Morlet wavelets for time-frequency analysis. NeuroImage, 199, 81-86.