脑机接口：尝试通过脑电预处理提高特征识别准确率

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的识别坏成分并处理的手段。寒假的学习可能就告一段落了，希望下次进行代码实践时能有更大的提高。

---

参考文献：

【1】MNE-Python | 使用 ICLabel 模型进行 ICA 并自动修复样本

【2】手把手教你EEG脑电数据预处理-操作篇