目录
(二)数据集下载以及OpenBMI-master.zip工具箱介绍下载,matlab下载:
1、数据集介绍:
(一)数据集背景:
项目邀请54名健康受试者(年龄24-35岁;25名女性)参加了实验。其中38名被试者为初始脑机接口使用者,其他都有脑机接口实验的经验。所有被试都没有神经学、精神病学和其他相关疾病的病史。受试者坐在60(±5)cm的椅子上面对21英寸液晶显示器(刷新率: 60 Hz;分辨率: 1600 × 1200)。水平和垂直的近似视角分别为 37.7度和28.1度。在本项目只涉及到MI范例,故只说明MI范例的过程:对于每一个被试,每个实验的前3秒在监视器中心出现一个黑色的固定十字,之后出现左右箭头视觉提示,被试进行用合适的手抓握得想象任务4秒。每个任务结束后,屏幕空白时间为6s(±1.5 s)(见图1)。实验分为train和test两个阶段(注意这里并不是训练集和测试集);每个阶段有100个实验,同时平衡左手和右手的图像任务。本实验分为两天,分别为session01和session02,每个被试者在两天分别做100个实验。
Motor-imagery 任务过程
(二)数据集下载以及OpenBMI-master.zip工具箱介绍下载,matlab下载:
数据集下载
数据集下载之前先了解一下数据集细节:
基于脑电图(EEG)的脑机接口(BCI)系统主要分为三大范式:运动想象(MI)、事件相关电位(ERP)和稳态视觉诱发电位(SSVEP)。我们记录了脑电图信号
(1)大量受试者(54名受试者)
(2)多个会议(两个不同天数[session1,session2])
(3)使用多种范式(MI、ERP和SSVEP)
本次实验只需用到MI范式,下载的时候只需要下载MI范式的数据集即可。这个数据的文件名是sess01_subj1_EEG_MI.mat的格式 数据集及工具箱下载链接:https://ftp.cngb.org/pub/gigadb/pub/10.5524/100001_101000/100542/
两个session里面分别是54个被试的实验数据,看要求你们需要多少个被试的数据,自己下载。然后每个被试都有三种范式的数据,只需要下载MI范式即可。
OpenBMI工具箱
OpenBMI工具箱提供了多种功能和工具,包括:
- 脑电图(EEG)信号处理:包括数据采集、预处理、特征提取等。
- 实时脑信号处理:支持实时监测和处理脑部信号,用于实时控制应用程序或设备。
- 机器学习算法:提供多种机器学习算法,用于分析和解释脑部信号,实现脑控制接口。
- 开放数据集:包括多个开放的脑电数据集,用于研究和开发BMI算法。
下载链接在上面。
matlab下载
参考博客:MATLAB下载+安装教程
参考
参考论文:如果想详细了解数据集可以仔细研读一下这篇论文,会有很大收获。
2、数据预处理:
数据集结果以及详细介绍:
每一个文件,也就是单个session中的某一个subject,例如sess01_subj1_EEG_MI.mat,这个文件中其实包含了两个文件,通过Matlab读取的时候也能发现,一个是train,一个是test,但其实并不是真正意义上的训练集和测试集。如下图所示:
在Matlab打开之后可以看到数据的结构:
读取数据后,得到这么一个struct数据(不管是train还是test),因为脑电信号是随时间变化的图像,以下为部分理解
X:输入,可知维度
t:取样时间,这里1X100,表示有100个取样点
fs:取样频率
y_dec:y标签,1表示左,2是右。(这里也知道是一个二分类问题了)
y_logic:y标签的逻辑值
y_class:y标签的类
class:类别
chan:通道数,1X62表示是62个通道
time_interval:采样持续时间,这里是4000ms
这个时候,我们能知道特征和标签,之后会经常提到trial(这个概念对后续数据划分设计很重要),实际上对应就是y标签的数量,这里实际上是进行了100次间歇性实验,所以有100个结果,那么,一个train或者test对应就是100个trial
特征提取
看到x的维度非常奇怪,脑电信号有很多通道,自然是非常奇怪。我们要使用一些算法,把某些频率的信号进行增强,然后提取出其中重要的特征,然后作为处理后的输入。
常见的算法有CSP,FCSP等等,这些东西都已经集成在工具箱内,你需要看懂那些函数的输入和输出是什么,然后就可以对x进行滤波操作
代码(matlab)
clear;
sum_loss = double(0);
for i = 1:54
if i < 10
idx = strcat('0', string(i));
else
idx = string(i);
end
% 读取数据集
str = sprintf('E:/pre_process_data/session1/s%d/sess01_subj%s_EEG_MI.mat', i, idx);
[CNT_tr,CNT_te] = Load_MAT(str);
% 训练集预处理
CNT_tr = prep_selectChannels(CNT_tr,{'Name', {'FC5', 'FC3', 'FC1', 'FC2', 'FC4', 'FC6', 'C5', 'C3', 'C1', 'Cz', 'C2', 'C4', 'C6', 'CP5', 'CP3', 'CP1', 'CPz', 'CP2', 'CP4', 'CP6'}});
CNT_tr = prep_filter(CNT_tr,{'frequency',[8,30]});划分频率8到30赫兹
SMT_tr = prep_segmentation(CNT_tr,{'interval',[1000,3500]});%数据分段,10ms分一段,为250X20X100(trials)
%得到SMT下一步数据
%测试集预处理
CNT_te = prep_selectChannels(CNT_te,{'Index',1:20});
CNT_te = prep_filter(CNT_te,{'frequency',[8,30]});
SMT_te = prep_segmentation(CNT_te,{'interval',[1000,3500]});
%CSP训练
[CSP_tr,CSP_W] = func_csp(SMT_tr,{'nPatterns',[2]});%csp空间滤波,得到滤波后的训练集、滤波参数,滤波评分
FT_tr=func_featureExtraction(CSP_tr,{'feature','logvar'});%对滤波后的训练集进行特征提取
CF_PARAW=func_train(FT_tr,{'classifier','LDA'});%训练LDA模型,得到分类器参数
%CSP测试
CSP_te = func_projection(SMT_te,CSP_W);
FT_te = func_featureExtraction(CSP_te,{'feature','logvar'});
cf_out = func_predict(FT_te,CF_PARAW);
loss = eval_calLoss(FT_te.y_dec,cf_out);
%计算平均loss
sum_loss = sum_loss + double(loss);
loss
FT_tr_x = FT_tr.x;
FT_tr_y = FT_tr.y_dec;
FT_te_x = FT_te.x;
FT_te_y = FT_te.y_dec;
strx_tr = sprintf('D:/rzkx_lab/data1/train_sess1/sess01_subj%s_FT_tr_x.mat', idx);
stry_tr = sprintf('D:/rzkx_lab/data1/train_sess1/sess01_subj%s_FT_tr_y.mat', idx);
strx_te = sprintf('D:/rzkx_lab/data1/test_sess1/sess01_subj%s_FT_te_x.mat', idx);
stry_te = sprintf('D:/rzkx_lab/data1/test_sess1/sess01_subj%s_FT_te_y.mat', idx);
%保存处理好的结果
save(strx_tr,"FT_tr_x");
save(stry_tr,"FT_tr_y");
save(strx_te,"FT_te_x");
save(stry_te,"FT_te_y");
sprintf('第%d被试的loss是', i)
end
avg_loss = sum_loss / 54
然后你会发现在Matlab左侧出现了一堆变量
Matlab会存储你的每一行代码运行结果,这点我觉得比较可以的,然后,FT变量就是提取出来的特征变量(就是feature的缩写),比如FT_te
进入这个变量,你会发现某些变量变化了,比如x,现在变成了4X100大小,也就是说,CSP滤波提取出了4个特征。
你需要做的就是把这里的x和y_dec作为数据的特征和标签导出来,也就是代码部分里的save。
3、跨被试和不跨被试介绍:
(一)划分方式
假如我们把所有的数据集都已经提取过,然后我们需要对数据进行拼接和划分,确定哪些是模型的训练集和测试集.
这里介绍一下概念,不跨被试的意思就是,仅仅对单个被试者的数据进行建模和预测,不把其它的被试者数据拼接。跨被试的意思就是,将你手上被试者的数据拼接在一起,进行建模和预测
第一种不跨被试
把每个subj的train和test拼接在一起(也就是200个trial),作为不跨被试的数据集,然后对数据集进行交叉验证,至于是几折,则依据情况确定(这个概念不懂的自己去查)
第二种不跨被试
把每个subj的train和test拼接在一起(200个trial),然后以session01作为训练集,session02作为测试集
跨被试
1个分类器,把每个subj的train和test拼接在一起(200个trial),然后9个subj拼接在一起作为训练集,剩余一个作为测试集(其实是留一验证法)
所有算法针对跨被试和不跨被试的一般处理流程:
1、跨被试:针对某一个具体的session,每次采用留一法选取其中一个subject的数据作为测试集,将其余所有subject信息合并作为训练集,放入所设计的分类器中,最后以准确率的平均值作为最终结果呈现。
2、不跨被试(不跨session):针对每一个subject,放入分类算法中利用交差验证训练subject,以交叉验证得到的准确率的平均值作为最终呈现结果。
3、不跨被试(跨session):使用session01的数据作为训练集,选取session02中对应subject的数据作为测试集,放入分类算法中进行预测得出准确率,以平均值作为最终呈现结果。
4、分类器的选择:
我这里用到的有SVM、逻辑回归、随机森林、朴素贝叶斯、KNN,你们也可以用其他更多的分类器模型。
这个是跨被试的参考代码。
import statistics
from scipy.io import loadmat
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score
from sklearn import svm
from sklearn.neighbors import KNeighborsClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.ensemble import RandomForestClassifier
import pyecharts.options as opts
from pyecharts.charts import Bar
import pandas as pd
import matplotlib.pyplot as plt
plt.rcParams["font.sans-serif"] = ["SimHei"]
plt.rcParams["axes.unicode_minus"] = False
# 将mat文件转化为dataframe,其中包含转置操作,返回值是dataframe数据集
def preprocess(data):
for i in data.keys():
if i == "FT_tr_x":
data = data["FT_tr_x"]
elif i == "FT_tr_y":
data = data["FT_tr_y"]
elif i == "FT_te_x":
data = data["FT_te_x"]
elif i == "FT_te_y":
data = data["FT_te_y"]
else:
continue
data = pd.DataFrame(data)
data = data.T
return data
# 拼接每个subj的train和test成为200个trial,返回值是200个trial
def one_mat(filename1, filename2):
x1_train = loadmat(filename1)
x1_train = preprocess(x1_train)
x2_train = loadmat(filename2)
x2_train = preprocess(x2_train)
x_train = pd.concat([x1_train, x2_train])
return x_train
# 每次选择一个i当作test,其余当作train数据,把中间数据放入一个list中返回
def left_one(all_data: list):
all_train_test_data = []
for i in range(len(all_data)):
test_data = all_data[i]
train_data = pd.DataFrame([])
for j in range(len(all_data)):
if i != j:
train_data = pd.concat([train_data, all_data[j]])
all_train_test_data.append((train_data, test_data))
return all_train_test_data
# 通过svm计算 accuracy
def accuracy_model_svm(x_train, x_test, y_train, y_test):
model = svm.SVC(C=10, kernel='linear', probability=True)
model.fit(x_train, pd.DataFrame(y_train).values.ravel())
y_pre = model.predict(x_test)
accuracy = accuracy_score(y_test, y_pre)
return accuracy
# 通过logic计算 accuracy
def accuracy_model_logic(x_train, x_test, y_train, y_test):
log = LogisticRegression(max_iter=1000)
log.fit(x_train, pd.DataFrame(y_train).values.ravel())
y_pre = log.predict(x_test)
accuracy = accuracy_score(y_test, y_pre)
return accuracy
# 通过knn计算 accuracy
def accuracy_model_KNN(x_train, x_test, y_train, y_test, n_neibor):
knn = KNeighborsClassifier(n_neighbors=n_neibor)
knn.fit(x_train, pd.DataFrame(y_train).values.ravel())
y_pre = knn.predict(x_test)
accuracy = accuracy_score(y_test, y_pre)
return accuracy
# 通过 朴素贝叶斯 计算 accuracy
def accuracy_model_GaussianNB(x_train, x_test, y_train, y_test):
gaussian = GaussianNB()
gaussian.fit(x_train, pd.DataFrame(y_train).values.ravel())
y_pre = gaussian.predict(x_test)
accuracy = accuracy_score(y_test, y_pre)
return accuracy
# 通过 随机森林 计算 accuracy
def accuracy_model_RandomForest(x_train, x_test, y_train, y_test):
rf = RandomForestClassifier()
rf.fit(x_train, pd.DataFrame(y_train).values.ravel())
y_pre = rf.predict(x_test)
accuracy = accuracy_score(y_test, y_pre)
return accuracy
sess = 1
total = 54
folderpath_train = f'D:\\rzkx_lab\\data\\train_sess{sess}'
folderpath_test = f'D:\\rzkx_lab\\data\\test_sess{sess}'
x_train = []
y_train = []
# 获取所有的 x_train 数据和 y_train 数据
for i in range(1, total + 1):
i = f'0{i}' if i < 10 else str(i)
tr_data_x = f'{folderpath_train}\sess0{sess}_subj{i}_FT_tr_x.mat'
te_data_x = f'{folderpath_test}\sess0{sess}_subj{i}_FT_te_x.mat'
tr_data_y = f'{folderpath_train}\sess0{sess}_subj{i}_FT_tr_y.mat'
te_data_y = f'{folderpath_test}\sess0{sess}_subj{i}_FT_te_y.mat'
x_train.append(one_mat(tr_data_x, te_data_x))
y_train.append(one_mat(tr_data_y, te_data_y))
train_test_x_data = left_one(x_train)
train_test_y_data = left_one(y_train)
# 获取所有的 accuracy
accuracy_svm = []
accuracy_logic = []
accuracy_knn = []
accuracy_GaussianNB = []
accuracy_randomforest = []
for i in range(len(train_test_x_data)):
train_x_data, test_x_data = train_test_x_data[i]
train_y_data, test_y_data = train_test_y_data[i]
accuracy_svm.append(accuracy_model_svm(train_x_data, test_x_data, train_y_data, test_y_data))
accuracy_logic.append(accuracy_model_logic(train_x_data, test_x_data, train_y_data, test_y_data))
accuracy_knn.append(accuracy_model_KNN(train_x_data, test_x_data, train_y_data, test_y_data, 6))
accuracy_GaussianNB.append(accuracy_model_GaussianNB(train_x_data, test_x_data, train_y_data, test_y_data))
accuracy_randomforest.append(accuracy_model_RandomForest(train_x_data, test_x_data, train_y_data, test_y_data))
print(f'总被试人数: {total}')
print('svm: ', statistics.mean(accuracy_svm), statistics.stdev(accuracy_svm))
print('逻辑回归: ', statistics.mean(accuracy_logic), statistics.stdev(accuracy_logic))
print('knn: ', statistics.mean(accuracy_knn), statistics.stdev(accuracy_knn))
print('朴素贝叶斯: ', statistics.mean(accuracy_GaussianNB), statistics.stdev(accuracy_GaussianNB))
print('随机森林: ', statistics.mean(accuracy_randomforest), statistics.stdev(accuracy_randomforest))
x_labels = []
for i in range(total):
x_labels.append(f'sub-{i + 1}')
line = (
Bar(init_opts=opts.InitOpts(width="1600px"))
.add_xaxis(xaxis_data=x_labels)
.add_yaxis(series_name='svm', y_axis=accuracy_svm)
.add_yaxis(series_name='逻辑回归', y_axis=accuracy_logic)
.add_yaxis(series_name='knn', y_axis=accuracy_knn)
.add_yaxis(series_name='朴素贝叶斯', y_axis=accuracy_GaussianNB)
.add_yaxis(series_name='随机森林', y_axis=accuracy_randomforest)
.set_global_opts(title_opts=opts.TitleOpts(title=f'跨被试结果---session{sess}'), xaxis_opts=opts.AxisOpts(name='被试者编号'), yaxis_opts=opts.AxisOpts(name='accuracy(%)'))
.set_series_opts(label_opts=opts.LabelOpts(is_show=False))
)
line.render()
5、模型结果评估指标:
准确率:准确率是指模型正确预测的数量所占全部预测的数量的比例,其公式为:
其中TP是正确做出positive的数量,TN是正确做出Negative的数量,FP是指错误地将Negative判定为Positive的数量,FN是指错误地将Positive判定为Negative的数量。
6、结果(准确度+标准差)做表格:
这个的话你们就按照自己的结果和审美来做。注意加减标准差就行。
举例:
完整代码:
链接:https://pan.baidu.com/s/17bDPJYiSWP15-Prx_FTcZA?pwd=2268
提取码:2268