2023/3/1-2023/3/4 脑机接口学习内容一览:
这一篇博客里,将对来自于kaggle的P300数据集(P300-Dataset)开展研究,并且简化了任务流程,从识别P300电位对应字符到识别是否出现P300电位,因此准确率较高,处理难度主要出现在陌生数据集的处理方面。
数据集:
本次代码的数据集具体情况如下:
channelrange 值为1-8,表示有8个通道 fs 采样频率为250hz trial 试验次数为35次,data中包含每一个实验开始的样本点 data.y 在数据y你就有了所有的标记,这些标记表示这个时间点是否与刺激巧合, 而刺激是为了引起P300反应。 如果该值为0,则没有应用刺激。 如果该值为1,则应用了与预期字母不对应的行/列刺激。 如果该值为2,则存在与目标字母对应的行/列。 stims 刺激(具体意义未知) channelnames 通道名称:['FZ', 'CZ', 'P3', 'PZ', 'P4', 'PO7', 'PO8', 'OZ'] samples 8个通道的采样点 word 35次试验中每一次刺激出现的字母['TOKENMIRARJUJUYMANSOCINCOJUEGOQUESO'] data.trial 35个试验中每个试验开始的样本点 data.flash Flash由4个字段组成,[采样点id,持续时间,刺激,命中/nohit]。 第一个是指数时间点,刺激开始的地方。下面是时间点上的持续时间, 一个标记表示激活了什么刺激,如果刺激是目标字母,则命中/未命中。 格式[7486, 31, 11, 1]
处理思路:
对于这个数据集,我们根据flash的最后一个字段区分注释。
注释分为两个类别1和2,1表示命中,2为未命中。
接下来将提取前一段数据的1、2段特征来对后面数据段是否命中做出判断。
以下是处理数据的步骤:
1.将flash整理为data的注释,并将其转化为events对象。
2.根据events构建epoch进行特征提取。
3.根据提取特征进行后半段数据的预测(是否为P300刺激段) 。
因此我们得到了这个提取数据的函数:
def Extractive(file):
original_data = scio.loadmat(file)
samples = original_data['samples'].transpose()
ch_names = ['FZ', 'CZ', 'P3', 'PZ', 'P4', 'PO7', 'PO8', 'OZ']
sfreq = 250
info = mne.create_info(ch_names, sfreq)
raw = mne.io.RawArray(samples, info)
# raw.plot()
# plt.pause(0)
inside_data = original_data['data']
# print(inside_data.dtype)
# np.set_printoptions(threshold=np.inf)
# print(inside_data['flash'][0][0])
flash = inside_data['flash'][0][0]
# 创建注释
onset = flash[:, 0]/250
description = flash[:, 3]
annot = mne.Annotations(onset=onset, duration=0.125, description=description)
# print(annot)
raw.set_annotations(annot)
# print(samples.shape)
raw_train = raw.copy()
raw_test = raw.copy()
raw_train.crop(tmin=200, tmax=500)
raw_test.crop(tmin=500, tmax=800)
'''
raw_train.plot()
raw_test.plot()
plt.pause(0)
'''
return raw_train, raw_test
结果分析:
在经过最麻烦的提取数据工作之后,我们得到了mne支持的raw数据。这给我们的研究带来了很大的便利,因为根据之前的博客代码实践:对脑电信号进行特征提取并分类(二分类)中的函数处理方式,我们可以直接将提取过后的数据套入相同的流程中进行处理(实际上我也这样做了),虽然用了这种看起来有些偷懒的方法,但最后得到的准确率与我想象的有较大差别——只有83%左右,按照我的思路来说,P300这么明显的电位识别率不应该在80%上下。
如下是输出的准确率,分别是SVC,随机森林以及决策树方法:
Accuracy score: 0.8309178743961353
Accuracy score: 0.8309178743961353
Accuracy score: 0.7246376811594203
经过思考,我列出了以下可能造成误差的原因:
1.数据的提取存在不足,有些提供的信息并没有被我用上。
2.在raw.plot()之后,我发现flash的数据标注存在一定的问题——部分1和2的标注是重合的。
3.对P300的了解不足,没有选取合适的数据段。
4.缺少预处理步骤。
准确率低的原因一部分是由于我的处理存在不足,同时也与我对数据集的不了解,数据集制作者没有提供足够详细的信息有关。希望下一次我可以有所进步,尝试一下功率谱之外的提取特征的方法。
完整代码:
import scipy.io as scio
import numpy as np
import matplotlib.pyplot as plt
import mne
from mne_icalabel import label_components
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import FunctionTransformer
from mne.preprocessing import (ICA, create_eog_epochs, create_ecg_epochs, corrmap)
from sklearn.svm import SVC
from sklearn.tree import DecisionTreeClassifier
def Extractive(file):
original_data = scio.loadmat(file)
# print(original_data.keys())
'''
for i in original_data.keys():
print('the', i, 'is', original_data.get(i))
'''
'''
在此数据集中有几个需要关注的关键字:
channelrange 值为1-8,表示有8个通道
fs 采样频率为250hz
trial 试验次数为35次,data中包含每一个实验开始的样本点
data.y 在数据y你就有了所有的标记,这些标记表示这个时间点是否与刺激巧合,而刺激是为了引起P300反应。
如果该值为0,则没有应用刺激。如果该值为1,则应用了与预期字母不对应的行/列刺激。如果该值为2,则存在与目标字母对应的行/列。
stims 刺激(具体意义未知)
channelnames 通道名称:['FZ', 'CZ', 'P3', 'PZ', 'P4', 'PO7', 'PO8', 'OZ']
samples 8个通道的采样点
word 35次试验中每一次刺激出现的字母['TOKENMIRARJUJUYMANSOCINCOJUEGOQUESO']
data.trial 35个试验中每个试验开始的样本点
data.flash Flash由4个字段组成,[采样点id,持续时间,刺激,命中/nohit]。
第一个是指数时间点,刺激开始的地方。下面是时间点上的持续时间,一个标记表示激活了什么刺激,如果刺激是目标字母,则命中/未命中。
格式[7486, 31, 11, 1]
'''
# print(original_data['stims'].shape)
samples = original_data['samples'].transpose()
ch_names = ['FZ', 'CZ', 'P3', 'PZ', 'P4', 'PO7', 'PO8', 'OZ']
sfreq = 250
info = mne.create_info(ch_names, sfreq)
raw = mne.io.RawArray(samples, info)
# raw.plot()
# plt.pause(0)
inside_data = original_data['data']
# print(inside_data.dtype)
# np.set_printoptions(threshold=np.inf)
# print(inside_data['flash'][0][0])
flash = inside_data['flash'][0][0]
"""
根据flash的最后一个字段区分注释
注释分为两个类别1和2,1表示命中,2为未命中
接下来将提取前一段数据的1、2段特征来对后面数据段是否命中做出判断
以下是处理数据的步骤:
1.将flash整理为data的注释,并将其转化为events对象
2.根据events构建epoch进行特征提取
3.根据提取特征进行后半段数据的预测(是否为P300刺激段)
"""
# 创建注释
onset = flash[:, 0]/250
description = flash[:, 3]
annot = mne.Annotations(onset=onset, duration=0.125, description=description)
# print(annot)
raw.set_annotations(annot)
# print(samples.shape)
raw_train = raw.copy()
raw_test = raw.copy()
raw_train.crop(tmin=200, tmax=500)
raw_test.crop(tmin=500, tmax=800)
'''
raw_train.plot()
raw_test.plot()
plt.pause(0)
'''
return raw_train, raw_test
def transform(raw):
event_id = {'1': 1, '2': 2, '1/2': 2}
# 重设通道名称
channel_types = {'FZ': 'eeg', 'CZ': 'eeg', 'P3': 'eeg', 'PZ': 'eeg',
'P4': 'eeg', 'PO7': 'eeg', 'PO8': 'eeg', 'OZ': 'eeg'}
raw.set_channel_types(channel_types)
events, _ = mne.events_from_annotations(raw, event_id=event_id) # 将注释转化为events
# mne.viz.plot_events(events_train, event_id=event_id) # 绘制事件发生时间分布图
# plot一下估计每一个事件持续的长度
# raw.plot()
# plt.pause(0)
epochs = mne.Epochs(raw=raw, events=events, tmin=-0.125, tmax=0.875, event_id=event_id,
preload=True, event_repeated='drop')
return epochs, events
def eeg_power_band(epochs):
"""
该函数根据epochs的特定频段中的相对功率来创建eeg特征
使用welch方法可得到83%左右正确率,使用multi方法只得到70%左右,效果不是很好
"""
# 特定频带
FREQ_BANDS = {"delta": [0.5, 4.5],
"theta": [4.5, 8.5],
# "alpha": [8.5, 11.5],
"sigma": [11.5, 15.5],
"beta": [15.5, 30]}
spectrum = epochs.compute_psd(method='welch', picks='eeg', fmin=0.5, fmax=30., n_fft=64, n_overlap=10)
psds, freqs = spectrum.get_data(return_freqs=True)
# 归一化 PSDs
psds /= np.sum(psds, axis=-1, keepdims=True)
X = []
for fmin, fmax in FREQ_BANDS.values():
psds_band = psds[:, :, (freqs >= fmin) & (freqs < fmax)].mean(axis=-1)
X.append(psds_band.reshape(len(psds), -1))
# print(X)
return np.concatenate(X, axis=1)
def classification(epochs_train, epochs_test, event_id, k):
if k == 0:
pipe = pipe = make_pipeline(FunctionTransformer(eeg_power_band, validate=False),
SVC(C=1.2, kernel='linear', random_state=42))
elif k == 1:
pipe = pipe = make_pipeline(FunctionTransformer(eeg_power_band, validate=False),
RandomForestClassifier(n_estimators=100, random_state=42))
elif k == 2:
pipe = pipe = make_pipeline(FunctionTransformer(eeg_power_band, validate=False),
DecisionTreeClassifier(random_state=42))
# 训练
y_train = epochs_train.events[:, 2]
pipe.fit(epochs_train, y_train)
# 预测
y_pred = pipe.predict(epochs_test)
# 评估准确率
y_test = epochs_test.events[:, 2]
acc = accuracy_score(y_test, y_pred)
print("Accuracy score: {}".format(acc))
def main():
file = 'C:/Users/86136/Desktop/innovation/python_mne/data/archive/P300S01.mat'
# 从数据集中提取数据
raw_train, raw_test = Extractive(file)
# 生成epochs和events对象
epochs_train, events_train = transform(raw_train)
epochs_test, events_test = transform(raw_test)
'''
特征工程:
将两个事件的epoch综合展示
'''
event_id = {'1': 1, '2': 2}
fig, (ax1, ax2) = plt.subplots(ncols=2)
stages = sorted(event_id.keys())
for ax, title, epochs in zip([ax1, ax2], ['train', 'test'], [epochs_train, epochs_test]):
for stage, color in zip(stages, ['red', 'blue']):
epochs[stage].plot_psd(area_mode=None, color=color, ax=ax,
fmin=0.1, fmax=20., show=False,
average=True, spatial_colors=False)
ax.set(title=title, xlabel='Frequency (Hz)')
ax2.set(ylabel='uV^2/hz (dB)')
ax2.legend(ax2.lines[2::3], stages)
plt.tight_layout()
# plt.show()
'''
分类预测
'''
for k in range(3):
classification(epochs_train, epochs_test, event_id, k)
main()