2023/2/1 -2/4脑机接口学习内容一览:
本文主要记录 在改进上个星期所发博客对脑电信号进行特征提取并分类(二分类)中的预处理流程来提高特征识别并分类的准确率 的过程中,对预处理方式的一些自己的理解与展示。
原代码:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import mne
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 pre(raw):
return raw
def transform(raw, event_id):
picks = mne.pick_types(raw.info, eeg=True)
# channel_types = {'EOG1': 'eog', 'EOG2': 'eog'}
# raw.set_channel_types(channel_types)
reject_criteria = dict(eeg=100e-6) # 100 µV
flat_criteria = dict(eeg=1e-6) # 1 µV
# 根据print得到线索将annotation转化为events
events, _ = mne.events_from_annotations(raw, event_id=event_id, chunk_duration=None)
# print(events)
# 绘制事件数据
# mne.viz.plot_events(events, event_id=event_id,
# sfreq=raw.info['sfreq'])
'''
根据事件创建epochs
在这个部分,提取epochs的各项参数需要自己仔细确定
比如tmin和tmax这两个值,可以通过raw的plot图像来大致确定
picks这里选取的是eeg通道
同时使用了基线校正
'''
tmin, tmax = -0.2, 0.3
epochs = mne.Epochs(raw=raw, events=events, event_id=event_id, tmin=tmin,
tmax=tmax, baseline=(None, 0), preload=True, picks=picks)
# epochs.plot_drop_log()
# print(epochs.info)
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))
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'))
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())
# 训练
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))
print(classification_report(y_test, y_pred, target_names=event_id.keys()))
def main(file):
raw = mne.io.read_raw_eeglab(file, preload=True, uint16_codec=None)
print(raw.info)
# raw.plot()
# 切出训练和测试集
event_id = {'rt': 1, 'square': 2}
raw_train = raw.copy()
raw_train = raw_train.crop(0, 100)
raw_test = raw.crop(100, 200)
epochs_train, events_train = transform(raw_train, event_id)
epochs_test, events_tset = transform(raw_test, event_id)
'''
特征工程:
将两个事件的epoch综合展示
'''
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)
if __name__ == "__main__" :
file = "C:\\Users\\86136\\Desktop\\innovation\\python_mne\\eeglab-read-set\\secdata.set"
main(file)
代码中使用的sklearn方法分别为SVC,随机森林与决策树分类。分类效果如下:
Accuracy score: 0.8307692307692308
Accuracy score: 0.7846153846153846
Accuracy score: 0.6307692307692307
可以看到如果我们使用基本没有处理过的数据,SVC分类的准确率最高。
方法1:
处理方式:
根据官方文档Rejecting bad data spans and breaks,我们采用第一种方式来处理数据。
本方案主要通过规定某通道类型的信号最高波幅来筛去不合理的epochs,因为所读取的文件的通道类型被全部标注为eeg,因此只限制eeg信号。
reject_criteria = dict(eeg=100e-6) # 100 µV
flat_criteria = dict(eeg=1e-6) # 1 µV
epochs = mne.Epochs(raw=raw, events=events, event_id=event_id, tmin=tmin,
tmax=tmax, baseline=(None, 0), preload=True, picks=picks, reject_tmax=0,
reject=reject_criteria, flat=flat_criteria, reject_by_annotation=False)
结果分析:
Accuracy score: 0.8163265306122449
Accuracy score: 0.7551020408163265
Accuracy score: 0.673469387755102
在筛除不合理数据后识别正确率如上所示,整体来看有所降低,但是决策树准确率提高。在预处理过程中,两次预处理分别有17,16个epochs被认定为bad。
方法2:
处理方式:
使用ICA去除伪影。
在查询该方法的使用方式时,我发现mne中也同样存在类似于eeglab插件中的自动ICA检测方法,在下载mne_icalabel库之后即可使用。参考资料为使用 ICLabel 模型进行 ICA 自动修复样本,直接下载库会很慢,建议使用清华源下载。
假设有15个成分时,对于训练集和测试集来说,每个成分的icalabel分别如下:
['brain', 'eye blink', 'brain', 'brain', 'brain', 'brain', 'line noise', 'brain', 'line noise', 'brain', 'line noise', 'brain', 'brain', 'line noise', 'brain']
['brain', 'brain', 'eye blink', 'brain', 'brain', 'brain', 'brain', 'brain', 'brain', 'brain', 'line noise', 'brain', 'brain', 'line noise', 'line noise']
可以看到各个主成分类型的数量基本相同,操作所用预处理函数替换如下:
def transform(raw, event_id):
picks = mne.pick_types(raw.info, eeg=True)
channel_types = {'EOG1': 'eog', 'EOG2': 'eog'}
raw.set_channel_types(channel_types)
raw.rename_channels({'FPz': 'Fpz'})
raw.set_montage('standard_1020')
montage = mne.channels.make_standard_montage('standard_1005')
raw.set_montage(montage)
raw.filter(l_freq=1., h_freq=None)
raw = raw.set_eeg_reference("average")
ica = mne.preprocessing.ICA(n_components=15, method='infomax', max_iter='auto',
random_state=97, fit_params=dict(extended=True))
ica.fit(raw)
raw.load_data()
ica.plot_sources(raw, show_scrollbars=False, show=True)
ic_labels = label_components(raw, ica, method="iclabel")
print(ic_labels["labels"])
# ica.plot_properties(raw, picks=[0, 12], verbose=False)
labels = ic_labels["labels"]
exclude_idx = [idx for idx, label in enumerate(labels) if label not in ["brain", "other"]]
print(f"Excluding these ICA components: {exclude_idx}")
reconst_raw = raw.copy()
ica.apply(reconst_raw, exclude=exclude_idx)
raw.plot(show_scrollbars=False)
reconst_raw.plot(show_scrollbars=False)
events, _ = mne.events_from_annotations(reconst_raw, event_id=event_id, chunk_duration=None)
# print(events)
# 绘制事件数据
# mne.viz.plot_events(events, event_id=event_id,
# sfreq=raw.info['sfreq'])
'''
根据事件创建epochs
在这个部分,提取epochs的各项参数需要自己仔细确定
比如tmin和tmax这两个值,可以通过raw的plot图像来大致确定
picks这里选取的是eeg通道
同时使用了基线校正
'''
tmin, tmax = -0.2, 0.3
epochs = mne.Epochs(raw=reconst_raw, events=events, event_id=event_id, tmin=tmin,
tmax=tmax, baseline=(None, 0), preload=True, picks=picks, reject_tmax=0)
epochs.plot_drop_log()
# print(epochs.info)
return epochs, events
在进行ica操作之后三种方法识别准确率如下:
Accuracy score: 0.7692307692307693
Accuracy score: 0.8615384615384616
Accuracy score: 0.6153846153846154
基于ar模型的白化可能会提高ica的效率,但是回头看看ar模型还有好多不清楚的内容。
方法一与方法二一起用时随机森林的识别准确率下降了一些,可能是去除的坏成分中含有部分有用成分。
操作总结:
在本次的代码实践中,虽然识别的准确率并没有明显的上升,但是对预处理流程更加熟悉,相比之前也更加了解了在mne中ica的识别坏成分并处理的手段。寒假的学习可能就告一段落了,希望下次进行代码实践时能有更大的提高。