基于CSP的运动想象 EEG 特征提取和可视化

news2024/10/7 7:31:57

基于运动想象的公开数据集：Data set IVa (BCI Competition III)¹
数据描述参考前文：https://blog.csdn.net/qq_43811536/article/details/134224005?spm=1001.2014.3001.5501
EEG 信号时频空域分析参考前文：https://blog.csdn.net/qq_43811536/article/details/134273470?spm=1001.2014.3001.5501

本文使用公开数据集 Data set IVa 中的部分被试数据，数据已公开可以从网盘获取：
链接：https://pan.quark.cn/s/5425ee5918f4
提取码：hJFz

1. 实验介绍

本任务的实验数据来自一名健康受试者，代号al。受试者在视觉提示出现后3.5s内完成以下3个运动想象中的一个：（L）左手，（R）右手，（F）右脚。分类任务中的数据只包括了右手和右脚两类，共280个试次。实验过程中使用脑电帽记录了118个通道的EEG信号，电极位置如图1所示。采集到的EEG信号首先经过带通滤波（0.05-200Hz）,再经过数字化和下采样，得到采样率为100Hz的信号。

在这里插入图片描述

图1 电极位置

2. EEG 特征提取和可视化

2.1 比较两类信号的差异

人在运动想象的过程中，大脑皮层会产生两种变化明显的节律信号，分别是 8-15Hz 的 μ 节律信号和 18-24Hz 的 β 节律。在运动想象时，大脑皮层对侧运动感觉区的脑电节律能量会明显降低，而同侧运动感觉区的脑电节律能量增大，这种现象称为事件相关去同步（Event Related Desynchronization，ERD）/事件相关同步（Event Related Synchronization，ERS）。本任务要区分的两类运动想象动作都是右侧，区别在于部位不同。

我们首先计算了两类运动想象的平均功率谱，得到图1所示的结果。可以看出，两类运动想象在C3和C4通道的5-25Hz具有较大的差异，这反映的正是运动感觉区的 μ 节律和 β 节律。我们还计算了5-25Hz的功率谱拓扑图，如图2所示。

在这里插入图片描述

(a) (b) 图1 两类信号的功率谱。（a）C3通道。（b）C4通道。

在这里插入图片描述

图2 功率谱拓扑图。左图为右手运动想象，右图为右脚运动想象。

2.2 CSP特征计算

2.2.1 频域特征

接下来，我们对EEG进行了带通滤波（12-28Hz），之后通过CSP 算法计算运动想象脑电特征，分量数目设置为10。图3展示了CSP模式图。

在这里插入图片描述

图3 CSP模式图

为了验证所提取到的CSP特征对于两类运动想象的区分性，绘制了两类样本的CSP特征。图4（a）展示了CSP特征的前两个维度，图4（b）则对10维CSP特征进行了t-SNE可视化。可以看出提取到的CSP特征具有较高的线性区分性。

在这里插入图片描述

(a) (b) 图4 两类样本的CSP特征可视化。（a）CSP特征的前两个维度。（b）CSP特征经过t-SNE降维。

3. 分析代码

主要使用Python的MNE包进行分析
MNE Tutorials 官网：https://mne.tools/stable/auto_tutorials/index.html
部分变量说明：
- raw：由 mne.io.RawArray() 函数创建，代表原始EEG数据
- epochs：由 mne.Epochs() 函数创建，代表一个事件（event）对应的所有数据，在该数据集中一个事件即 “右手”或者“脚”的想象运动

raw = mne.io.RawArray(eeg_data.T, info)
from scipy.signal import welch


# raw.filter(12, 28, fir_design="firwin", skip_by_annotation="edge")

# Events and epochs
# Events
events = np.vstack((cues_pos, np.zeros(len(cues_pos)), target_label[0, :])).T.astype(int)
picks = mne.pick_types(raw.info, meg=False, eeg=True, stim=False, eog=False, exclude="bads")
# Cue onset
# Epochs
epochs = mne.Epochs(
    raw,
    events,
    events_id,
    0,
    3.5,
    proj=True,
    picks=picks,
    baseline=None,
    preload=True,
)
for i_ch in range(len(ch_name)):
    f, pxx0 = welch(epochs.get_data()[right_id,i_ch,:], fs=eeg_fs)
    f, pxx1 = welch(epochs.get_data()[foot_id,i_ch,:], fs=eeg_fs)
    plt.figure()
    plt.plot(f,10 * np.log10(pxx0.mean(axis=0)),'r',linewidth=2)
    plt.plot(f,10 * np.log10(pxx1.mean(axis=0)),'g',linewidth=2)
    plt.xlim([5,40])
    plt.xlabel('[Hz]',fontsize=16)
    plt.ylabel('[dB]',fontsize=16)
    plt.legend(['Right','Foot'],fontsize=16)
    plt.xticks(fontsize=16)  # 设置 x 轴刻度标签的字体大小为 12
    plt.yticks(fontsize=16)  # 设置 y 轴刻度标签的字体大小为 12
    plt.title(list(ch_name)[i_ch][0], fontsize=16)
    plt.tick_params(axis='both', which='major', width=2)  # 设置主要刻度线的宽度为 2
    # 设置顶部边框的粗细
    plt.gca().spines['top'].set_linewidth(2)

    # 设置底部边框的粗细
    plt.gca().spines['bottom'].set_linewidth(2)

    # 设置左侧边框的粗细
    plt.gca().spines['left'].set_linewidth(2)

    # 设置右侧边框的粗细
    plt.gca().spines['right'].set_linewidth(2)

    plt.savefig(f'Results\\topomap\\{list(ch_name)[i_ch][0]}.png')

plt.figure()
epochs[right_id].compute_psd().plot_topomap(vlim=(155,180),bands={"5-25Hz":(5,25)},ch_type="eeg", agg_fun=np.mean, dB=True)
plt.savefig(f'Results\\topomap_r.png')
plt.show()

plt.figure()
epochs[foot_id].compute_psd().plot_topomap(vlim=(155,180),bands={"5-25Hz":(5,25)},ch_type="eeg", agg_fun=np.mean, dB=True)

plt.savefig(f'Results\\topomap_f.png')

plt.show()

plt.figure()
raw.compute_psd(tmin=0,tmax=3.5).plot_topomap(vlim=(155,180),bands={"5-25Hz":(5,25)},ch_type="eeg", agg_fun=np.mean, dB=True)

plt.savefig(f'Results\\topomap_f.png')

x = epochs.get_data()
y = target_label[0, :] 
csp = CSP(n_components=10, reg=None, log=True, norm_trace=False)
csp.fit_transform(x, y)
csp.plot_patterns(epochs.info, ch_type="eeg", units="Patterns (AU)", size=1.5)
plt.show()