进行分类"/>
利用python对脑电数据进行分类
本篇文章利用到的数据是mne包中的多导睡眠图(PSG)数据进行睡眠阶段分类
作者处于大二阶段,文中解读字数较少,可能一些东西讲的不清楚,大家直接看代码就行,这也是本人第一次关于脑电数据的实战,参考MNE文档官方的例子,链接如下:
Sleep stage classification from polysomnography (PSG) data — MNE 1.5.1 documentation --- 根据多导睡眠图 (PSG) 数据对睡眠阶段进行分类 — MNE 1.5.1 文档e/stable/auto_tutorials/clinical/60_sleep.html
希望通过本文能帮助你了解脑电数据分类的大致流程
一、数据集介绍
给定Sleep Physionet数据集2 3中的两个受试者,我们默认命名为Alice和Bob(实际上可能并不是这两个人),从Alice的数据中预测出Bob的睡眠有多少。这是一个监督式的多类分类任务来解决,从30秒的数据中5个可能阶段预测睡眠阶段
利用dataset包中的sleep_physionet.age.fetch_data()加载数据。
-PSG.edf
包含多导睡眠图。来自脑电图头盔的原始数据
-Hypnogram.edf
包含专家记录的注释
将这两者结合在一个 mne.io.Raw
对象中,然后我们可以根据注释的描述提取事件以获得纪元
Sleep Physionet 数据集使用 8 个标签进行注释:Wake (W)、Stage 1、Stage 2、Stage 3、Stage 4,对应于从浅睡眠到深度睡眠的范围,REM 睡眠 (R),其中 REM 是快速眼动睡眠、运动 (M) 和任何未评分片段的阶段 (?) 的缩写
但是这里我们只处理五个阶段,分别是:Wake (W)、stage1、stage2、stage3/4 和REM sleep (R)
二、数据处理与特征工程
1.导包
from mne.datasets.sleep_physionet.age import fetch_data
import mne
import numpy as np
import matplotlib.pyplot as plt
from mne.time_frequency import psd_array_welch
from sklearn.preprocessing import LabelEncoder
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 sklearn.ensemble import RandomForestClassifier
from xgboost import XGBClassifier
from sklearn.svm import SVC
from sklearn.model_selection import cross_val_score
2.读取PSG数据并创建对象
# 对实验对象进行ALICE BOB进行临时命名
ALICE,BOB = 0,1
# 加载两位受试者的实验数据 recording表示夜间记录的编号 有效值为[1] [2] [1,2]
[alice_files,bob_files] = fetch_data(subjects=[ALICE,BOB],recording=[1])
# 通道名称
mapping = {'EOG horizontal': 'eog','Resp oro-nasal': 'misc','EMG submental': 'misc','Temp rectal': 'misc','Event marker': 'misc'}# 读取ALICE的edf文件
raw_train = mne.io.read_raw_edf(alice_files[0])
annot_train = mne.read_annotations(alice_files[1])
print(raw_train)
print(annot_train)
print(alice_files[0])
raw_train.set_annotations(annot_train,emit_warning=False)
raw_train.set_channel_types(mapping)# 绘制空0s开始 时间窗口为40s的连续通道数据波形图
raw_train.plot(duration=40,scalings='auto')
plt.show()
3.绘制睡眠与事件的图像
利用绘图直观地得出事件的分布
#进行睡眠事件与事件映射annotation_desc_2_event_id = {'Sleep stage W': 1,'Sleep stage 1': 2,'Sleep stage 2': 3,'Sleep stage 3': 4,'Sleep stage 4': 4,'Sleep stage R': 5}events_train,_ = mne.events_from_annotations(raw_train,event_id=annotation_desc_2_event_id,chunk_duration=30
)
# 自定义事件标签
event_id = {'Sleep stage W': 1,'Sleep stage 1': 2,'Sleep stage 2': 3,'Sleep stage 3/4': 4,'Sleep stage R': 5}# 绘制事件数据
mne.viz.plot_events(events_train,event_id=event_id,sfreq=raw_train.info['sfreq'])# 保留颜色代码以便进一步绘制
stage_colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
4.创建Epoch对象
tmax = 30 - 1 / raw_train.info['sfreq']# 所创建的是时间从tmin=0开始,到tmax为止epochsepochs_train = mne.Epochs(raw=raw_train,events=events_train,event_id=event_id,tmin=0,tmax=tmax,baseline=None)print(epochs_train)
raw_test = mne.io.read_raw_edf(bob_files[0])
annot_test = mne.read_annotations(bob_files[1])
raw_test.set_annotations(annot_test,emit_warning=False)
raw_test.set_channel_types(mapping)
events_test,_ = mne.events_from_annotations(raw_test,event_id=annotation_desc_2_event_id,chunk_duration=30)
epochs_test = mne.Epochs(raw=raw_test,events=events_test,event_id=event_id,tmin=0,tmax=tmax,baseline=None)
print(epochs_test)
print(epochs_test.get_data)
进行功率谱密度(PSD)图像的绘制
fig,(ax1,ax2) = plt.subplots(ncols=2)stages = sorted(event_id.keys())
for ax,title,epochs in zip([ax1,ax2],['Alice','Bob'],[epochs_train,epochs_test]):for stage,color in zip(stages,stage_colors):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)')
plt.tight_layout()
plt.show()
三、从python函数设计scikit-learn转换器
1.设计一个函数用于Epochs的处理
def eeg_power_band(epochs):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 = epochspute_psd(picks='eeg',fmin=0.5,fmax=30.0)psds,freqs = spectrum.get_data(return_freqs=True) # 归一化 PSDspsds /= 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))return np.concatenate(X,axis=1)
2.构造pipe
将函数和随机森林分类器注入到pipe当中构建出一个管道,并利用Bob的数据进行模型评估得出准确率,该过程运行的事件可能比较长。
利用循环遍历不同的param得到树木的棵树对精确度的影响,并将其存储在列表当中
params = [10,20,30,40,50,60,70,80,90,100,110,120,130,140,150]
acc_list = []
for param in params:pipe = make_pipeline(FunctionTransformer(eeg_power_band,validate=False),RandomForestClassifier(n_estimators=param,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)acc_list.append(acc)
绘制出树木棵数与分类器的准确率之间的关系
plt.plot(params,acc_list)
plt.show()
四、总结与思考
是否利用SVM、CSP、ICA去伪迹得到的效果会更好呢?这个是后面继续做的改进与测试,争取提高获取较高准确率。
本文只是简单地对脑电数据进行处理与分类预测,欢迎大家发表看法。最后,致敬脑机领域的所有人。
更多推荐
利用python对脑电数据进行分类
发布评论