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

2020-06-29 17:52:22 浏览数 (1)

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

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

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

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

第一步:导入相应的包

代码语言:javascript复制
import os
import sys
from collections import OrderedDict

import pandas as pd
import numpy as np
import seaborn as sns
from matplotlib import pyplot as plt
%matplotlib inline

from eeg_utils.utils import utils as utils

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

第二步:加载数据

代码语言:javascript复制
subject = 1
session = 1
raw = utils.load_data('auditory/P300', sfreq=256., 
                      subject_nb=subject, session_nb=session)

查看功率谱(Power Spectrum)

代码语言:javascript复制
raw.plot_psd(tmax=np.inf);

第三步:过滤

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

代码语言:javascript复制
raw.filter(1, 30, method='iir')

Epoching

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

代码语言:javascript复制
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

代码语言:javascript复制
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可能是针对二进制和非平衡分类问题的最佳度量标准)

代码语言:javascript复制
from sklearn.pipeline import make_pipeline

from mne.decoding import Vectorizer

from sklearn.linear_model import LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA

from sklearn.model_selection import cross_val_score, StratifiedShuffleSplit

from pyriemann.estimation import ERPCovariances, XdawnCovariances
from pyriemann.tangentspace import TangentSpace
from 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 data
epochs.pick_types(eeg=True)
X = epochs.get_data() * 1e6
times = epochs.times
y = 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
代码语言:javascript复制
fig = plt.figure(figsize=[8,4])
sns.barplot(data=results, x='AUC', y='Method')
plt.xlim(0.4, 0.75)
sns.despine()

0 人点赞