Python中MNE库利用CSP分析运动想象数据

最新推荐文章于 2024-07-31 19:34:02 发布

zhoudapeng01

最新推荐文章于 2024-07-31 19:34:02 发布

阅读量1.2w

点赞数 17

分类专栏： MNE 算法 MI

本文链接：https://blog.csdn.net/zhoudapeng01/article/details/103979381

版权

算法同时被 3 个专栏收录

19 篇文章 7 订阅

订阅专栏

MNE

15 篇文章 52 订阅

订阅专栏

1 篇文章 0 订阅

订阅专栏

最近在研究EEG运动想象相关的内容，都说CSP在运动想象领域很好用，刚好最近接触了MNE库，就尝试了一下在这个库中的实现。牛皮就是牛皮，有现成的方法可以调用，既然有方法那当然也有示例了。

官方给的示例，后面的内容就是对文档中的示例程序进行分析。

https://mne.tools/stable/auto_examples/decoding/plot_decoding_csp_eeg.html#sphx-glr-auto-examples-decoding-plot-decoding-csp-eeg-py

第一部分：导入数据

使用官方提供的示例数据，将其转换为scikit中可以使用的数据格式。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.pipeline import Pipeline
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.model_selection import ShuffleSplit, cross_val_score
from mne import Epochs, pick_types, events_from_annotations
from mne.channels import read_layout
from mne.io import concatenate_raws, read_raw_edf
from mne.datasets import eegbci
from mne.decoding import CSP

###############################################################################
tmin, tmax = -1., 4.  #设置参数，记录点的前1秒后4秒用于生成epoch数据
event_id = dict(hands=2, feet=3) #设置事件的映射关系
subject = 1
runs = [6, 10, 14]
# 获取想要读取的文件名称，这个应该是没有会默认下载的数据
raw_fnames = eegbci.load_data(subject, runs)
# 将3个文件的数据进行拼接
raw = concatenate_raws([read_raw_edf(f, preload=True) for f in raw_fnames])
# 去掉通道名称后面的（.），不知道为什么默认情况下raw.info['ch_names']中的通道名后面有的点
raw.rename_channels(lambda x: x.strip('.'))
# 对原始数据进行FIR带通滤波
raw.filter(7., 30., fir_design='firwin', skip_by_annotation='edge')#
# 从annotation中获取事件信息
events, _ = events_from_annotations(raw, event_id=dict(T1=2, T2=3))
# 剔除坏道，提取其中有效的EEG数据
picks = pick_types(raw.info, meg=False, eeg=True, stim=False, eog=False, exclude='bads')
# 根据事件生成对应的Epochs数据
epochs = Epochs(raw, events, event_id, tmin, tmax, proj=True, picks=picks, baseline=None, preload=True)
# 截取其中的1秒到2秒之间的数据，也就是提示音后1秒到2秒之间的数据（这个在后面滑动窗口验证的时候有用）
epochs_train = epochs.copy().crop(tmin=1., tmax=2.)
# 将events转换为labels,event为2,3经过计算后也就是0,1
labels = epochs.events[:, -1] - 2

第二部分：利用CSP提取特征，建立线性分类器

利用CSP算法提取特征，将特征作为输入，结合线性分类器进行运动想象的分类，这里使用了pipeline的方式可以简化实现过程。还有一个点比较有意思，官方给出csp的fit_transform函数第二个参数可以为空，但实际上如果你运行csp.fit_transform(epochs_data)是会报错的，个人认为文档描述有问题，csp应该是有两个参数的，因为公共空间模式（CSP）算法采用监督的方法创建一个最优的公共空间滤波器，在最大化一类方差的同时最小化另一类方差，采用同时对角化两类任务协方差矩阵的方式，得到可区分程度最大的特征向量。适用于二分类任务的特征提取。既然是监督的那labels就一定要传递csp.fit_transform(epochs_data, labels)

###########################第二部分，特征提取和分类####################################################
scores = []
# 获取epochs的所有数据，主要用于后面的滑动窗口验证
epochs_data = epochs.get_data()
# 获取训练数据
epochs_data_train = epochs_train.get_data()
# 设置交叉验证模型的参数
cv = ShuffleSplit(10, test_size=0.2, random_state=42)
# 根据设计的交叉验证参数,分配相关的训练集和测试集数据
cv_split = cv.split(epochs_data_train)
# 创建线性分类器
lda = LinearDiscriminantAnalysis()
# 创建CSP提取特征，这里使用4个分量的CSP
csp = CSP(n_components=4, reg=None, log=False, norm_trace=False)
# 创建机器学习的Pipeline,也就是分类模型，使用这种方式可以把特征提取和分类统一整合到了clf中
clf = Pipeline([('CSP', csp), ('LDA', lda)])
# 获取交叉验证模型的得分
scores = cross_val_score(clf, epochs_data_train, labels, cv=cv, n_jobs=1)
# 输出结果，准确率和不同样本的占比
class_balance = np.mean(labels == labels[0])
class_balance = max(class_balance, 1. - class_balance)
print("Classification accuracy: %f / Chance level: %f" % (np.mean(scores), class_balance))
# csp提取特征,用于绘制CSP不同分量的模式图（地形图）
# 如果没有这一步csp.plot_patterns将不会执行
csp.fit_transform(epochs_data, labels)
# lay文件的存放路径，这个文件不是计算生成的，是mne库提供的点击分布描述文件在安装路径下（根据个人安装路径查找）：
# D:\ProgramData\Anaconda3\Lib\site-packages\mne\channels\data\layouts\EEG1005.lay
layout = read_layout('EEG1005')
csp.plot_patterns(epochs.info, layout=layout, ch_type='eeg', units='Patterns (AU)', size=1.5)

第三部分：验证算法的性能

还记得前面在提取数据得时候，epochs_data是事件发生前1秒到发生后4秒之间的数据，而我们的训练数据epochs_data_train使用的是事件发生后1秒到发生后2秒之间的数据。有趣的是测试数据这部分，这里以0.5秒为窗长，0.1秒为步长在epochs_data上滑动生成。不知道你发现没有这里的测试数据时间窗口0.5秒，而训练数据的时间窗口为1秒，纳尼，这意味这训练数据和测试数据的维度不一致啊！！！训练集数据是64*161，然鹅我们的测试数据是64*80，竟然不一样，竟然还能不一样。详细的分析放到后面，这一部分先看一下验证的结果。结果符合预期，测试集窗口滑动到1秒以后才能和训练集对应，识别的准确率才会提高。

#########################验证算法的性能###########################################
# 获取数据的采样频率
sfreq = raw.info['sfreq']
# 设置滑动窗口的长度，也就是数据窗口的长度
w_length = int(sfreq * 0.5)
# 设置滑动步长，每次滑动的数据间隔
w_step = int(sfreq * 0.1)
# 每次滑动窗口的起始点
w_start = np.arange(0, epochs_data.shape[2] - w_length, w_step)
# 得分列表用于保存模型得分
scores_windows = []
# 交叉验证计算模型的性能
for train_idx, test_idx in cv_split:
    # 获取测试集和训练集数据
    y_train, y_test = labels[train_idx], labels[test_idx]
    # 设置csp模型的参数，提取相关特征，用于后面的lda分类
    X_train = csp.fit_transform(epochs_data_train[train_idx], y_train)
    # 拟合lda模型
    lda.fit(X_train, y_train)
    # 用于记录本次交叉验证的得分
    score_this_window = []
    for n in w_start:
        # csp提取测试数据相关特征
        X_test = csp.transform(epochs_data[test_idx][:, :, n:(n + w_length)])
        # 获取测试数据得分
        score_this_window.append(lda.score(X_test, y_test))
    # 添加到总得分列表
    scores_windows.append(score_this_window)

# 设置绘图的时间轴，时间轴上的标志点为窗口的中间位置
w_times = (w_start + w_length / 2.) / sfreq + epochs.tmin
# 绘制模型分类结果的性能图（得分的均值）
plt.figure()
plt.plot(w_times, np.mean(scores_windows, 0), label='Score')
plt.axvline(0, linestyle='--', color='k', label='Onset')
plt.axhline(0.5, linestyle='-', color='k', label='Chance')
plt.xlabel('time (s)')
plt.ylabel('classification accuracy')
plt.title('Classification score over time (w_length = {0}s)'.format(w_length/sfreq))
plt.legend(loc='lower right')
plt.show()

下面来看一下为什么测试数据集的时间窗口和训练数据集的时间窗口不一致也可以运行呢？lda的输入是经过csp提取后的特征，所以问题的重点是为什么csp的输入时间窗口可以不一致呢?参考下面的博文中csp提取特征的结果

https://blog.csdn.net/oh__NO/article/details/84310982

fi为提取的特征，是由zi计算得来，而zi又是由矩阵W和Xi计算得来，看到这大概明白了，W是我们的投影矩阵，Xi是我们的输入数据，根据矩阵乘法的运算规则：前一个的列数等于后一个的行数，不难看出即使Xi的列数发生改变也不会影响计算。这里Xi的列数对应的就是时间窗口，我想这也就是为什么测试集的时间窗口和训练集的时间窗口可以不一致的原因吧。下面分别设置测试集时间窗口为1秒（和训练集相同）和1.5秒，和上面0.5秒的时候相比，在这3种时间窗口中1秒的更加平稳，我想这也是其和训练集时间窗口一致的原因，毕竟测试集和训练集在数据维度上应该保持一致。不过这个特性倒是十分适合EEG这种具有时间序列特点的数据，这种使用方法或许在其他方面有意想不到的效果，呵呵呵。