脑信号分析系列(1)-听觉P300实验 / 四六文摘

更多技术，第一时间送达

听觉P300实验与视觉P300相似，但使用听觉刺激来产生oddball

刺激时间为200ms，时间间隔400ms,随机抖动±100ms, 任务是计算玩奇数球刺激的次数，记录单个参与者进行的6次2分钟的实验。

非目标刺激的音调约为523Hz(C5)，目标刺激的音调约为1174 Hz(D6)

下面利用muse库来进行数据分析

第一步：导入相应的包

import osimport sysfrom collections import OrderedDict

import pandas as pdimport numpy as npimport seaborn as snsfrom matplotlib import pyplot as plt%matplotlib inline

from eeg_utils.utils import utils as utils

import warningswarnings.filterwarnings("ignore", category=Warning)

第二步：加载数据

subject = 1session = 1raw = utils.load_data('auditory/P300', sfreq=256., subject_nb=subject, session_nb=session)

查看功率谱(Power Spectrum)

raw.plot_psd(tmax=np.inf);

第三步：过滤

过滤掉1到30Hz之间的数据

raw.filter(1, 30, method='iir')

Epoching

我们在刺激之后将数据从-100ms缩短到800ms，无需进行基线校正(信号经过带通滤波)，并且在信号超过75uV时，我们拒绝每个epochs，主要是拒绝眨眼。

from mne import Epochs, find_events

events = find_events(raw)event_id = {'Non-Target': 1, 'Target': 2}

epochs = Epochs(raw, events=events, event_id=event_id, tmin=-0.1, tmax=0.8, baseline=None, reject={'eeg': 75e-6}, preload=True, verbose=False, picks=[0,1,2,3]

Epoch average

我们绘制两种情况下的平均ERP

conditions = OrderedDict()conditions['Non-target'] = [1]conditions['Target'] = [2]

fig, ax = utils.plot_conditions(epochs, conditions=conditions, ci=97.5, n_boot=1000, title='', diff_waveform=(1, 2))

在两种情况下都可以看到清晰的P300

解码

通过上述的平均epochs，可以很清楚识别ERP。但如何了解有关P300的SNR的任何信息，可以通过分类管道(classification pipline)了解P300响应的强度。

下面我们将使用4个不同的管道。

Vect+RL: 实验矢量化+Logistic回归。可以将其视为MEG/EEG中的标准解码管道。

Vect+RegLDA:实验矢量化+正则化LDA.P300 BCI中经常利用这一点，但维数过多，则可能无法使用。

ERPCov+MDM: ERP协方差+MDM.黎曼几何分类器，一种简单有效的方法(用于较少的通道数)

ERPCov+TS: ERP协方差+切线空间映射。这是一种基于黎曼几何的管道之一。

以AUC作为度量，以交叉验证的方式进行评估(AUC可能是针对二进制和非平衡分类问题的最佳度量标准)

from sklearn.pipeline import make_pipeline

from mne.decoding import Vectorizer

from sklearn.linear_model import LogisticRegressionfrom sklearn.preprocessing import StandardScalerfrom sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA

from sklearn.model_selection import cross_val_score, StratifiedShuffleSplit

from pyriemann.estimation import ERPCovariances, XdawnCovariancesfrom pyriemann.tangentspace import TangentSpacefrom pyriemann.classification import MDM

clfs = OrderedDict()

clfs['Vect + LR'] = make_pipeline(Vectorizer(), StandardScaler(), LogisticRegression())clfs['Vect + RegLDA'] = make_pipeline(Vectorizer(), LDA(shrinkage='auto', solver='eigen'))clfs['ERPCov + TS'] = make_pipeline(ERPCovariances(estimator='oas'), TangentSpace(), LogisticRegression())clfs['ERPCov + MDM'] = make_pipeline(ERPCovariances(estimator='oas'), MDM())clfs['XdawnCov + TS'] = make_pipeline(XdawnCovariances(estimator='oas'), TangentSpace(), LogisticRegression())clfs['XdawnCov + MDM'] = make_pipeline(XdawnCovariances(estimator='oas'), MDM())

# format dataepochs.pick_types(eeg=True)X = epochs.get_data() * 1e6times = epochs.timesy = epochs.events[:, -1]

# 定义交叉验证cv = StratifiedShuffleSplit(n_splits=20, test_size=0.25, random_state=42)

# 对每个管道(pipeline)运行交叉验证auc = []methods = []for m in clfs: try: res = cross_val_score(clfs[m], X, y==2, scoring='roc_auc', cv=cv, n_jobs=-1) auc.extend(res) methods.extend([m]*len(res)) except: pass results = pd.DataFrame(data=auc, columns=['AUC'])results['Method'] = methods

fig = plt.figure(figsize=[8,4])sns.barplot(data=results, x='AUC', y='Method')plt.xlim(0.4, 0.75)sns.despine()

脑机接口BCI爱好者交流群：QQ群：903290195

脑信号分析系列(1)-听觉P300实验

Epoching

Epoch average

解码

相关推荐