Python中MNE库利⽤CSP分析运动想象数据
最近在研究EEG运动想象相关的内容，都说CSP在运动想象领域很好⽤，刚好最近接触了MNE库，就尝试了⼀下在这个库中的实现。⽜⽪就是⽜⽪，有现成的⽅法可以调⽤，既然有⽅法那当然也有⽰例了。
官⽅给的⽰例，后⾯的内容就是对⽂档中的⽰例程序进⾏分析。
第⼀部分：导⼊数据
使⽤官⽅提供的⽰例数据，将其转换为scikit中可以使⽤的数据格式。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.pipeline import Pipeline
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
del_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']中的通道名后⾯有的点
# 对原始数据进⾏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 = py().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 = _data()
# 获取训练数据
epochs_data_train = _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 = ansform(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 + in
# 绘制模型分类结果的性能图（得分的均值）
plt.figure()
plt.plot(w_times, np.mean(scores_windows, 0), label='Score')
python默认安装路径
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提取特征的结果
fi为提取的特征，是由zi计算得来，⽽zi⼜是由矩阵W和Xi计算得来，看到这⼤概明⽩了，W是我们的投影矩阵，Xi是我们的输⼊数据，根据矩阵乘法的运算规则：前⼀个的列数等于后⼀个的⾏数，不难看出即使Xi的列数发⽣改变也不会影响计算。这⾥Xi的列数对应的就是时间窗⼝，我想这也就是为什么测试集的时间窗⼝和训练集的时间窗⼝可以不⼀致的原因吧。下⾯分别设置测试集时间窗⼝为1秒（和训练集相同）和1.5秒，和上⾯0.5秒的时候相⽐，在这3种时间窗⼝中1秒的更加平稳，我想这也是其和训练集时间窗⼝⼀致的原因，毕竟测试集和训练集在数据维度上应该保持⼀致。不过这个特性倒是⼗分适合EEG这种具有时间序列特点的数据，这种使⽤⽅法或许在其他⽅⾯有意想不到的效果，呵呵呵。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。

Linux+Django+Python+Wsgi配置过程
« 上一篇
anacondapycharm_Anaconda及Pycharm的安装及配置教程
下一篇 »

发表评论