MNE入门

2024-04-29  本文已影响0人  韧心222

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)

其中:

该函数的返回值类型为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)

其中

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)

%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]

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

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]] 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

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

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)

其中

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)

其中：

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)
              )

下面的代码将创建这些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)

其中：

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 包

上一篇 下一篇

猜你喜欢

热点阅读