1. MNE的安装

下面的方法只能安装2D版本，无法安装3D版本（不过我还没有搞懂是什么意思）

pip install mne

除此以外，还需要安装numpy、pandas和matplotlib，建议安装jupyter notebook和pyqt5

pip install numpy pandas matplotlib jupyter notebook PyQt5

2. 读取数据

2.1 read_raw_eeglab

MNE提供了多种格式的数据读取方法，都放在mne的io中，在此我们先试用针对EEGLab的数据读取方法read_raw_eeglab，其函数定义为：

mne.io.read_raw_eeglab(input_fname, eog=(), preload=False, uint16_codec=None, montage_units='auto', verbose=None)

其中:

• input_fname 要读取的文件名，一般以set为后缀
• eog EOG通道的通道名称或索引。如果设置为'auto'，则使用包含EOG或EYE的通道名称。默认值为空元组。
• preload 如果为真，数据将被预加载到内存中。如果preload是一个字符串，那么preload是存储在硬盘上的内存映射文件的文件名，用于存储数据
• montage_units 通道位置的表示单位。默认值为'mm'（毫米）

该函数的返回值类型为mne.io.Raw

import os
import numpy as np
import pandas as pd
import mne
import matplotlib.pyplot as plt

data_dir = 'D:\\BrainTools\\eeglab2024.0\\sample_data\\'
data_path = data_dir + 'eeglab_data.set'

raw = mne.io.read_raw_eeglab(data_path, preload=True)

读取之后的结果为：

read_eeglab_set.png

2.2 read_custom_montage

用于读取用户定义的脑电图（EEG）电极定位图（montage）。这个函数允许用户导入自己的电极定位信息，以便在MNE-Python中进行数据分析，其定义为：

mne.channels.read_custom_montage(fname, head_size=0.095, coord_frame=None)

其中

• fname 文件名，应该以特定的扩展名之一结尾，这些扩展名对应于不同的电极定位文件格式。以下是每个扩展名对应的数据文件格式的说明：
  • .loc 或 .locs 或 .eloc: 这些扩展名通常用于 EEGLAB 格式的文件，其中包含电极的局部坐标。
  • .sfp: 用于 BESA/EGI 格式的文件，其中包含电极的球面坐标。
  • .csd 或 .elc: 这些扩展名通常用于 BESA 格式的文件，其中包含电极的球面坐标。
  • .txt: 纯文本文件，其中包含电极的名称和它们的位置。
  • .csv 或 .tsv: 这些扩展名用于逗号分隔值（CSV）或制表符分隔值（TSV）文件，其中包含电极的名称和它们的位置。
  • .xyz: 纯文本文件，其中包含电极的名称和它们的三维坐标（x, y, z）。
• head_size 用于指定头部的大小（半径，单位为米）。如果 head_size 设置为 None，则会使用montage文件中设置的值。

locs_info_path = data_dir + "eeglab_chan32.locs"
# 读取电极位置信息
montage = mne.channels.read_custom_montage(locs_info_path)
# 传入数据的电极位置信息
raw.set_montage(montage)

3 数据可视化

如果使用Jupyter的话，需要使用如下命令与MNE进行交互

%matplotlib qt

3.1 绘制脑电原始数据

MNE中最常用的绘制原始脑电数据的命令为plot，其定义为：

plot(events=None, duration=10.0, start=0.0, n_channels=20, bgcolor='w', color=None, bad_color='lightgray', event_color='cyan', scalings=None, remove_dc=True, order=None, show_options=False, title=None, show=True, block=False, highpass=None, lowpass=None, filtorder=4, clipping=1.5, show_first_samp=False, proj=True, group_by='type', butterfly=False, decim='auto', noise_cov=None, event_id=None, show_scrollbars=True, show_scalebars=True, time_format='float', precompute=None, use_opengl=None, *, picks=None, theme=None, overview_mode=None, splash=True, verbose=None)

• duration 要绘制的时窗大小（单位：秒）。这个值和原始文件持续时间中较小的那个将被使用。
• n_channels 同时要绘制的通道数。默认值为20。n_channels 和 len(raw.ch_names) 中较小的那个将被显示。如果 order 是 ‘position’、‘selection’ 或 ‘butterfly’，则此参数将无效果。
• clipping 如果设置为 None，则允许通道值超过其在图中的指定范围。如果设置为"clamp"，则值会被限制在适合显示的范围内。如果设置为"transparent"，则超出范围的值不会显示。
• butterfly 在 MNE-Python 中，butterfly模式是一种特殊的绘图方式，它允许用户在同一图中显示多个通道的波形，以直观地比较它们的波形。在 butterfly模式下，每个通道的波形都被压缩到一个垂直的条形中，这样多个通道的波形就可以并排显示，从而方便地比较它们的形状和幅度。

%matplotlib qt
raw.plot(duration=5, n_channels=16, clipping=None)
plt.show()

3.2 绘制功率谱或振幅谱

使用.compute_psd().plot()函数为每种通道类型绘制单独的图。

compute_psd(method='welch', fmin=0, fmax=inf, tmin=None, tmax=None, picks=None, exclude=(), proj=False, remove_dc=True, reject_by_annotation=True, *, n_jobs=1, verbose=None, **method_kw)[source]

• fmin, fmax 低通和高通滤波器的频率
• tmin, tmax 计算的时域范围，如果为 None，则使用数据中的第一个或最后一个时间点
• picks 要包括的通道，可以使用切片或者整数的列表来表示通道索引。在列表中，可以使用通道类型字符串（例如，[‘meg’, ‘eeg’]）或通道名称字符串（例如，[‘MEG0111’, ‘MEG2623’]）。如果设置为“all”则选择所有通道，设置为“data”则选择数据通道。默认值为None，选择良好的数据通道（不包括参考MEG通道）。请注意，如果在info[‘bads’]出现的通道，在列表中，则这些通道也会被选择。
• exclude 要排除的通道名称。如果设置为'bads'，则 info['bads']中的通道将被排除；传递一个空列表,则包括所有通道（包括'bad'通道）
• remove_dc 如果设置为 True，在计算每个段的频谱之前，将从每个段中减去其平均值。
• reject_by_annotation 是否在进行频谱估计之前省略数据中的坏段。如果设置为True，则以bad开头的描述的注释所标记的段将被省略。

计算之后，可以使用plot()函数进行绘制，其主要参数有:

• average 如果设置为False，将显示所有通道的功率谱密度（PSD），同时参数area_mode和area_alpha将不起作用。此时，可以绘制一个区域（按住左键并拖动）来绘制一个头皮地形图。
  其他绘制参数，可以参考已经过时的plot_psd函数

plot_psd(fmin=0, fmax=inf, tmin=None, tmax=None, picks=None, proj=False, reject_by_annotation=True, *, method='auto', average=False, dB=True, estimate='auto', xscale='linear', area_mode='std', area_alpha=0.33, color='black', line_alpha=None, spatial_colors=True, sphere=None, exclude='bads', ax=None, show=True, n_jobs=1, verbose=None, **method_kw)

# 频率范围设置为10-50Hz，绘制平均功率谱
raw.compute_psd(fmin=10,fmax=50).plot(average=True)

psd_plot.png

# 绘制功率谱
raw.compute_psd().plot()

psd_plot2.png

4. 数据预处理

4.1 伪迹检测（Artifact Detection）

主要参考了https://mne.tools/stable/auto_tutorials/preprocessing/10_preprocessing_overview.html和https://github.com/ZitongLu1996/Python-EEG-Handbook-CN/blob/master/EEG_CN_Python_Handbook_Preprocessing.ipynb

伪迹（也被称为伪影）是脑电记录信号中的一部分，它们来源于除了我们感兴趣的源（即大脑中的神经活动）之外的其他源。例如：

• 环境伪迹
  • 持续的振荡围绕交流电源线频率（通常是50或60 Hz）
  • 由于建筑物振动（例如门砰地关上）导致的信号跳跃
  • 来自附近电梯、手机、地磁场等附近的电磁场噪声
• 仪器伪迹（Instrumentation artifacts）
  • 电磁干扰来自刺激呈现（例如EEG传感器捕捉到未屏蔽耳机产生的场）：如果实验中使用的耳机没有被适当地屏蔽，它们的电磁场可能会被EEG传感器捕捉到，从而在记录的数据中引入额外的噪声。这种类型的干扰可能导致数据中出现异常的信号波动，这些波动与大脑活动无关。
  • 头部位置指示器（HPI）线圈使用的特定频率的连续振荡：HPI线圈用于追踪受试者的头部位置，以便在数据分析时考虑头部移动的影响。这些线圈在特定频率上工作，有时这些频率的信号会作为伪迹出现在EEG数据中。这些伪迹通常表现为在特定频率上的连续振荡，它们可能与大脑活动信号相混淆，因此需要在预处理步骤中识别并去除。
  • 由于传感器故障（例如表面电极、头皮接触不良）导致的单个通道中的随机高幅度波动（或者恒定零信号）：当EEG传感器出现故障时，如电极与头皮接触不良或电极本身损坏，可能会导致数据中出现随机的高幅度波动，或者在某些情况下，信号完全丢失，表现为恒定的零信号。这些伪迹会影响单个通道的数据质量，并且在进行后续分析之前需要被检测和处理。
• 生物学伪迹
  • 周期性QRS样信号模式（特别是在磁强计通道中）由于心脏的电活动：在EEG（脑电图）数据中，由心脏的电活动产生的周期性QRS样信号模式是一种常见的生物学伪迹。这些信号模式通常在磁强计通道中更为明显，因为它们对心脏产生的磁场较为敏感。这些周期性的信号可能与心律有关，且在频率和形态上类似于心电图（ECG）中的QRS复合波。
  • 由于眼动导致的短暂的阶跃式偏转（特别是在前额EEG通道中）：眼动，尤其是快速的眼睛移动或眨眼，会在EEG数据中产生短暂的阶跃式偏转。这种类型的伪迹在前额区域的EEG通道中尤为常见，因为这些通道更接近于眼睛，更容易捕捉到与眼动相关的电活动。
  • 由于眨眼导致的大的瞬时偏转（特别是在前额EEG通道中）：当受试者眨眼时，会在EEG数据中产生大的瞬时偏转。这些偏转通常是突发的，并且幅度较大，它们在前额EEG通道中尤其明显，因为眨眼产生的电活动在这些区域有较强的影响。
  • 跨越多个通道的高频波动的短暂爆发由于吞咽时的肌肉活动：在吞咽过程中，肌肉活动会产生短暂的高频波动，这些波动可能会影响多个EEG通道。这种类型的伪迹在时间上通常是短暂的，但由于涉及到多个肌肉群的协调活动，因此可能会影响到多个通道。

4.1.1 低频漂移（Low-frequency drifts）

在处理EEG（脑电图）数据时，低频漂移是一种常见的伪迹，它可能由多种因素引起，包括环境噪声、受试者的生理活动（如呼吸和心跳），以及设备本身的温度变化等。这些低频变化通常表现为数据中缓慢的、持续的电压变化。为了有效地检测这些低频漂移，研究人员通常会使用plot()函数来进行视觉检查。这种方法允许研究人员观察数据中是否存在任何不期望的、持续的趋势或模式。为了更好地识别这些低频变化，建议在一个较长的时间跨度内绘制数据，例如60秒，这样可以更容易地观察到数据中的长期趋势。

此处我们使用了mne提供的sample数据集，执行下面代码后，在界面上电极proj按钮以呈现更加清晰的漂移：

import os

import numpy as np

import mne

sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = os.path.join(
    sample_data_folder, "MEG", "sample", "sample_audvis_raw.fif"
)
raw = mne.io.read_raw_fif(sample_data_raw_file)
raw.crop(0, 60).load_data()  # just use a fraction of data for speed here

mag_channels = mne.pick_types(raw.info, meg="mag")
raw.plot(duration=60, order=mag_channels, n_channels=len(mag_channels), remove_dc=False)

Pasted image 20240430102122.png

在EEG（脑电图）数据预处理中，低频漂移通常可以通过应用高通滤波器来去除。高通滤波器可以消除低于特定截止频率的信号成分，从而减少或消除低频噪声，包括低频漂移。

在给定的情况中，观察到的低频漂移的波长大约是20秒，这意味着这些漂移的频率成分大约是0.05 Hz（因为频率是周期的倒数，即1/周期时间）。为了去除这些低频漂移，可以选择一个略高于这个频率的截止频率进行高通滤波。在这种情况下，我们选择一个0.1 Hz的截止频率，这个值足够低，可以去除大多数低频漂移，同时又不会过度影响可能与认知过程相关的低频脑电活动。

滤波函数为Raw.filter()，其定义为：

filter(l_freq, h_freq, picks=None, filter_length='auto', l_trans_bandwidth='auto', h_trans_bandwidth='auto', n_jobs=None, method='fir', iir_params=None, phase='zero', fir_window='hamming', fir_design='firwin', skip_by_annotation=('edge', 'bad_acq_skip'), pad='reflect_limited', verbose=None)

其中

• l_freq 对于FIR（有限脉冲响应）滤波器来说，该参数表示通带的下边缘；对于IIR（无限脉冲响应）滤波器来说，该参数是下截止频率。如果为None，则数据只进行低通滤波。
• h_freq 对于FIR（有限脉冲响应）滤波器，是指通带上边缘；对于IIR（无限脉冲响应）滤波器，是指上截止频率。如果为None，则数据只进行高通滤波。
• method 指定滤波器的方法，可选的参数包括fir和iir。

4.1.2 电力线噪声

用前面提到的compute_psd()可以很容易发现电力线的伪迹。如前面绘制过的psd图像中，可以很容易看到在60hz处的伪迹。

psd_plot2.png

同样地，我们也可以看到sample数据集中的伪迹：

fig = raw.compute_psd(tmax=np.inf, fmax=250).plot(
    average=True, amplitude=False, picks="data", exclude="bads"
)
# add some arrows at 60 Hz and its harmonics:
for ax in fig.axes[1:]:
    freqs = ax.lines[-1].get_xdata()
    psds = ax.lines[-1].get_ydata()
    for freq in (60, 120, 180, 240):
        idx = np.searchsorted(freqs, freq)
        ax.arrow(
            x=freqs[idx],
            y=psds[idx] + 18,
            dx=0,
            dy=-12,
            color="red",
            width=0.1,
            head_width=3,
            length_includes_head=True,
        )

在这里，我们看到了在60、120、180和240 Hz处的狭窄频率峰值——这是美国（样本数据记录的地方）的电源线频率及其第二、第三和第四谐波。其他峰值（大约在25到30 Hz，以及这些频率的第二谐波）可能与心跳有关，这在使用专门的心跳检测功能中的时间域中更容易看到。

psd_plot3.png

4.1.3 心跳伪迹（ECG）

MNE-Python 在 mne.preprocessing 子模块中包含了一个专门的函数find_ecg_events()，用于从ECG通道或磁力计（magnetometers，不知道翻译的是否准确）（如果没有 ECG 通道）中检测心跳伪迹。此外，函数create_ecg_epochs()会调用find_ecg_events()，并使用返回的事件数组（[[MNE——Events]]）来提取心跳epoch（[[MNE——Epoch]]），默认以event为中心前后时长为-0.5s~0.5s。

mne.preprocessing.create_ecg_epochs函数的定义为：

mne.preprocessing.create_ecg_epochs(raw, ch_name=None, event_id=999, picks=None, tmin=-0.5, tmax=0.5, l_freq=8, h_freq=16, reject=None, flat=None, baseline=None, preload=True, keep_ecg=False, reject_by_annotation=True, decim=1, verbose=None)

其中：

• ch_name 用于检测ECG的通道名称。当设置为None时，如果存在ECG通道，就会使用ECG通道；如果没有ECG通道，则将从跨通道的平均值创建一个合成ECG通道。这个合成通道只能从MEG通道创建
• picks 设置要包括的通道。
  • 如果该参数为整数列表（或切片），则这些整数将被解释为通道索引
  • 如果该参数被设置为字符串列，需要根据字符串的含义来判断包括的通道：如果字符串表示通道类型（例如['meg', 'eeg']），则将选择这些类型的通道；如果字符串表示通道名称（例如['MEG0111', 'MEG2623']），则将选择给定名称的通道。也可以使用'all'字符串值来选择所有通道，或使用'data'来选择数据通道。
  • 默认情况下该参数为None，表示选择所有通道。请注意，如果在info['bads']中的通道名称或索引在picks参数中被明确指定，那么这些通道也将被包括在内
• l_freq 在检测事件时应用于ECG通道的低通频率
• h_freq 在检测事件时应用于ECG通道的高通频率
• reject 基于最大峰-峰值（peak-to-peak，PTP）信号——即最低和最高信号值之间的绝对差值——幅度拒绝（过滤）epochs。其具体做法为，在每个单独的epoch中，为每个通道计算该epoch的PTP值。如果PTP超过拒绝阈值，则相应的epoch将被丢弃。字典键对应于不同的通道类型；有效键可以是对象中出现的任何通道类型。例如：

reject = dict(grad=4000e-13,  # unit: T / m (gradiometers)
              mag=4e-12,      # unit: T (magnetometers)
              eeg=40e-6,      # unit: V (EEG channels)
              eog=250e-6      # unit: V (EOG channels)
              )

• baseline 在应用基线校正时，视为“基线”的时间间隔。如果为None，则不应用基线校正。如果baseline为一个元组（a, b），则该间隔位于a和b之间（以秒为单位，包括端点）。如果a是None，那么从数据的开始部分开始；如果b是None，则包括数据的结束部分。如果为（None, None），则使用整个时间间隔作为基线。对每个epoch和通道单独应用校正的方法如下：1) 计算基线期间的平均信号。2)从整个epoch中减去这个平均值。

下面的代码将创建这些ECGepoch，然后显示检测到的ECG伪迹的图像，以及所有伪迹的平均 ERF。我们将显示所有三种通道类型，需要注意的是EEG通道受到心跳伪迹的影响较小。

ecg_epochs = mne.preprocessing.create_ecg_epochs(raw)
ecg_epochs.plot_image(combine="mean")

磁力计图像中的水平条纹反映了心跳伪迹叠加在我们在前面章节中看到的低频漂移上的事实；为了避免这种情况，可以在调用 create_ecg_epochs() 时设置 baseline=(-0.5, -0.2)。还可以通过使用 mne.Evoked.plot_topomap() 方法对 ECG epoch 进行平均，快速查看传感器之间与 ECG 相关的场模式。

Pasted image 20240506153340.png

Pasted image 20240506153350.png

Pasted image 20240506153357.png

avg_ecg_epochs = ecg_epochs.average().apply_baseline((-0.5, -0.2))

在得到avg_ecg_epochs后，可以直观地看到与EOG响应峰值相关的不同时间点的场空间分布模式：

avg_ecg_epochs.plot_topomap(times=np.linspace(-0.05, 0.05, 11))

场空间分布模式

也可以使用plot()函数得到一个ERP/F图，或者使用plot_joint()函数得到一个结合的头皮电场图和ERP/F图。代码中已经手动指定了头皮电场图的时间，但如果没有提供，它们将根据信号中的峰值自动选择：

avg_ecg_epochs.plot_joint(times=[-0.25, -0.025, 0, 0.025, 0.25])

Pasted image 20240506155628.png

Pasted image 20240506155634.png

Pasted image 20240506155642.png

4.1.4 眼动伪迹（ Ocular artifacts）

与心跳伪迹检测方法相类似，MNE-Python 也包含了用于检测和提取眼动伪迹的函数：find_eog_events() 和create_eog_epochs()。

eog_epochs = mne.preprocessing.create_eog_epochs(raw, baseline=(-0.5, -0.2))
eog_epochs.plot_image(combine="mean")
eog_epochs.average().plot_joint()

Pasted image 20240514091046.png

Pasted image 20240514091054.png

Pasted image 20240514091100.png

Pasted image 20240514091103.png

Pasted image 20240514091108.png

Pasted image 20240514091112.png

在采集EEG时，有时会因为一些特殊原因，导致个别通道出现故障，提供的数据太嘈杂，无法使用，这些被我们称为故障通道（bad channel）。在MNE中，故障通道的列表通常Raw、Epochs 和 Evoked对象的Info对象的bads字段（该字段为一个普通列表，可以按照常规的python列表进行处理）中。仍然沿用前面的例子：

import os
from copy import deepcopy

import numpy as np

import mne

sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = os.path.join(
    sample_data_folder, "MEG", "sample", "sample_audvis_raw.fif"
)
raw = mne.io.read_raw_fif(sample_data_raw_file, verbose=False, preload=True)

print(raw.info["bads"])

可以很清楚的看到故障通道为：

['MEG 2443', 'EEG 053']

此时，我们可以绘制故障通道的时间序列。在此，我们使用了pick_channels_regexp()函数，利用正则表达式挑选了所有EEG 05开头的通道，并用plot()进行绘制

picks = mne.pick_channels_regexp(raw.ch_names, regexp="EEG 05.")
raw.plot(order=picks, n_channels=len(picks))

Pasted image 20240514092838.png

同样地，我们也可以绘制另一个故障通道：

picks = mne.pick_channels_regexp(raw.ch_names, regexp="MEG 2..3")
raw.plot(order=picks, n_channels=len(picks))

Pasted image 20240514093402.png

可以看到，故障通道通常以灰色显示，从图中可以看到EEG 053 没有拾取头皮电位，而 MEG 2443 则比其邻居有更多的内部噪声。这些都是故障通道的特征。
在调用plot()函数后，可以在通道上点击，使其变为bad通道，这样会更新raw对象，当关闭绘图窗口时会对标记信息进行保留。

MNE的诸多函数中，均包括exclude参数，该参数的默认值一般为'bad'（即exclude='bads'），用于排除故障通道，如果要在计算中涵盖故障通道，则需要将该参数设置为空，即exclude=[]

The bad EEG channel is not so obvious, but the bad gradiometer is easy to see.
故障通道不太明显，但梯度计很容易看出。

插值

对于坏通道的处理，直接排除是一种简单高效的做法，但是更常规的做法是对其进行插值，这能够保证所有被试都有相同的数据维度。

在 Raw 对象中插值故障通道的方法为interpolate_bads() 。其定义为：

interpolate_bads(reset_bads=True, mode='accurate', origin='auto', method=None, exclude=(), verbose=None)

其中：

• reset_bads 如果设置为True，则将从info中将故障通道移除
• mode ‘accurate’ 或 ‘fast’，决定用更快的方法插值，还是更准确的方法插值
• method 对每种通道类型使用的插值方法。可以为dict、str或者为None。
  • 如果为str，则相当于所有通道使用相同的方法
  • 如果设置为None，则相当于

method=dict(meg="MNE", eeg="spline", fnirs="nearest")

示例代码如下所示：

eeg_data = raw.copy().pick(picks="eeg")
eeg_data_interp = eeg_data.copy().interpolate_bads(reset_bads=False)

for title, data in zip(["orig.", "interp."], [eeg_data, eeg_data_interp]):
    with mne.viz.use_browser_backend("matplotlib"):
        fig = data.plot(butterfly=True, color="#00000022", bad_color="r")
    fig.subplots_adjust(top=0.9)
    fig.suptitle(title, size="xx-large", weight="bold")

对于更自动化的故障通道检测和插值方法，可以考虑使用 autoreject 包