---

前言

       共空間模式演算法（CSP）廣泛應用於腦機介面（BCI）中腦電訊號（EEG）的特徵提取，因此，在介紹CSP演算法前，有必要對腦電介面，腦電訊號等相關概念做個簡述。

一、BCI與EEG的基礎概念

• BCI： 腦機介面（brain-computer interface， BCI），指在人或動物大腦與外部裝置之間建立的直接連線，實現腦與裝置的資訊交換；
  
• EEG： 腦電訊號（EEG）是一種 5-100μv和低頻的生物電訊號，需放大後才能顯示和處理。在腦電訊號處理與模式識別系統中，為了正確的識別 EEG 訊號，訊號的處理應該包括以下三個部分：預處理、特徵選擇與提取以及 特徵分類。訊號預處理主要為了去除低頻噪聲干擾，如利用空間濾波器（CAR）濾除眼電、肌電等低頻噪聲干擾。特徵提取與選擇主要是為了降低腦電資料的維數和提取出與分類相關的特徵。
         目前，EEG資料的特徵主要有三種：時域特徵、頻域特徵與空域特徵，不同的特徵需要採取不同的特徵提取方法，如空域特徵一般採用空域濾波器(共同空間模式,CSP)進行提取，頻域特徵一般採用傅立葉變換、小波變換或自迴歸(Auto-Regressive, AR)模型獲取。特徵分類主要是利用分類演算法對提取到的特徵進行分類，主要分為兩個步驟：首先，利用訓練樣本特徵進行模型的訓練，獲取分類的引數，然後，用訓練好的分類器來獲取測試樣本特徵的類別。目前，較常用的分類器有 Fisher、支援向量機(SVM)、神經網路分類器(Neural network classifier)和貝葉斯分類器(Bayesian classifier)等。
  
• EEG-BCI： 系統組成示意圖如下：
  

二、CSP演算法簡介

2.1 儲備知識

1. 空間濾波器（spatial filters）

       濾波是指通過操作接受或拒絕一定的頻率。具體到空間濾波器，指的是通過空間濾波器對矩陣的元素進行修改，達到預期的效果。
空間濾波器由一下兩部分組成：
1）一個鄰域（通常為規模較小矩陣）
2）對該鄰域的覆蓋的矩陣元素執行的預定義操作

下圖為使用3 × 3的鄰域的線性空間濾波舉例，過程就是對原矩陣的元素逐個按照預選定義的濾波器的操作進行運算：

2. 協方差與協方差矩陣

3. 方陣特徵值與特徵向量

4. 矩陣對角化及同時對角化

5. 白化

2.2 基本概念

共空間模式（CSP）： 一種對兩分類任務下的空域濾波特徵提取演算法，能夠從多通道的腦機介面資料裡面提取出每一類的空間分佈成分。
共空間模式演算法的基本原理： 利用矩陣的對角化，找到一組最優空間濾波器進行投影，使得兩類訊號的方差值差異最大化，從而得到具有較高區分度的特徵向量。
共空間模式演算法旨在： 設計空間濾波器使得兩組腦電時空訊號矩陣濾波後，方差值差異最大化，從而得到具有較高區分度的特徵向量。用於下一步將特徵向量送入分類器進行分類。

       上圖中，腦電時空訊號矩陣的維數N × T ，N代表腦電通道數（導聯數目），T代表時間層面的取樣點個數。一個矩陣的單位為一個trail，表示一次測試。濾波器的維數是2 m × N ，N同樣代表導聯數目，2m為生成該濾波器時，人為設定的特徵選取個數。

       上面說明了CSP演算法的結果的用途，下面說一下其生成的過程。其利用的基本原理是協方差矩陣對角化。輸入為訓練集的全部trail矩陣，每個矩陣對應的任務標籤（比如左手或右手運動想象）規定的特徵數。下面詳細介紹一下CSP演算法的步驟

• CSP演算法步驟：
  
  1） 混合空間協方差矩陣R：
  
  2）白化特徵值矩陣P：
  
  3）空間濾波器W：
  
  4）特徵提取f：
  

2.3 CSP演算法python例項

原始碼參考來源：Python中MNE庫利用CSP分析運動想象資料

import numpy as np
import matplotlib.pyplot as plt
from sklearn.pipeline import Pipeline
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.model_selection import ShuffleSplit, cross_val_score
from mne import Epochs, pick_types, events_from_annotations
from mne.channels import make_standard_montage
from mne.io import concatenate_raws, read_raw_edf
from mne.datasets import eegbci
from mne.decoding import CSP

###############################################################################
tmin, tmax = -1., 4.  # 設定引數，記錄點的前1秒後4秒用於生成epoch資料
event_id = dict(hands=2, feet=3) # 設定事件的對映關係
subject = 1
runs = [6, 10, 14]
# 獲取想要讀取的檔名稱，這個應該是沒有會預設下載的資料
raw_fnames = eegbci.load_data(subject, runs)
# 將3個檔案的資料進行拼接
raw = concatenate_raws([read_raw_edf(f, preload=True) for f in raw_fnames])
# 去掉通道名稱後面的（.），不知道為什麼預設情況下raw.info['ch_names']中的通道名後面有的點
eegbci.standardize(raw)
montage = make_standard_montage('standard_1005')
raw.set_montage(montage)

raw.rename_channels(lambda x: x.strip('.'))
# 對原始資料進行FIR帶通濾波
raw.filter(7., 30., fir_design='firwin', skip_by_annotation='edge')
# 從annotation中獲取事件資訊
events, _ = events_from_annotations(raw, event_id=dict(T1=2, T2=3))
# 剔除壞道，提取其中有效的EEG資料
picks = pick_types(raw.info, meg=False, eeg=True, stim=False, eog=False, exclude='bads')
# 根據事件生成對應的Epochs資料
epochs = Epochs(raw, events, event_id, tmin, tmax, proj=True, picks=picks, baseline=None, preload=True)
# 擷取其中的1秒到2秒之間的資料，也就是提示音後1秒到2秒之間的資料（這個在後面滑動視窗驗證的時候有用）
epochs_train = epochs.copy().crop(tmin=1., tmax=2.)
# 將events轉換為labels,event為2,3經過計算後也就是0,1
labels = epochs.events[:, -1] - 2
###########################第二部分，特徵提取和分類####################################################
scores = []
# 獲取epochs的所有資料，主要用於後面的滑動視窗驗證
epochs_data = epochs.get_data()
# 獲取訓練資料
epochs_data_train = epochs_train.get_data()
# 設定交叉驗證模型的引數
cv = ShuffleSplit(10, test_size=0.2, random_state=42)
# 根據設計的交叉驗證引數,分配相關的訓練集和測試集資料
cv_split = cv.split(epochs_data_train)
# 建立線性分類器
lda = LinearDiscriminantAnalysis()
# 建立CSP提取特徵，這裡使用4個分量的CSP
csp = CSP(n_components=4, reg=None, log=False, norm_trace=False)
# 建立機器學習的Pipeline,也就是分類模型，使用這種方式可以把特徵提取和分類統一整合到了clf中
clf = Pipeline([('CSP', csp), ('LDA', lda)])
# 獲取交叉驗證模型的得分
scores = cross_val_score(clf, epochs_data_train, labels, cv=cv, n_jobs=1)
# 輸出結果，準確率和不同樣本的佔比
class_balance = np.mean(labels == labels[0])
class_balance = max(class_balance, 1. - class_balance)
print("Classification accuracy: %f / Chance level: %f" % (np.mean(scores), class_balance))
# csp提取特徵,用於繪製CSP不同分量的模式圖（地形圖）
# 如果沒有這一步csp.plot_patterns將不會執行
csp.fit_transform(epochs_data, labels)
# lay檔案的存放路徑，這個檔案不是計算生成的，是mne庫提供的點選分佈描述檔案在安裝路徑下（根據個人安裝路徑查詢）：
# D:\ProgramData\Anaconda3\Lib\site-packages\mne\channels\data\layouts\EEG1005.lay
# layout = read_layout('EEG1005')
# csp.plot_patterns(epochs.info, layout=layout, ch_type='eeg', units='Patterns (AU)', size=1.5)
csp.plot_patterns(epochs.info, ch_type='eeg', units='Patterns (AU)', size=1.5)

#########################驗證演算法的效能###########################################
# 獲取資料的取樣頻率
sfreq = raw.info['sfreq']
# 設定滑動視窗的長度，也就是資料視窗的長度
w_length = int(sfreq * 0.5)
# 設定滑動步長，每次滑動的資料間隔
w_step = int(sfreq * 0.1)
# 每次滑動視窗的起始點
w_start = np.arange(0, epochs_data.shape[2] - w_length, w_step)
# 得分列表用於儲存模型得分
scores_windows = []
# 交叉驗證計算模型的效能
for train_idx, test_idx in cv_split:
    # 獲取測試集和訓練集資料
    y_train, y_test = labels[train_idx], labels[test_idx]
    # 設定csp模型的引數，提取相關特徵，用於後面的lda分類
    X_train = csp.fit_transform(epochs_data_train[train_idx], y_train)
    # 擬合lda模型
    lda.fit(X_train, y_train)
    # 用於記錄本次交叉驗證的得分
    score_this_window = []
    for n in w_start:
        # csp提取測試資料相關特徵
        X_test = csp.transform(epochs_data[test_idx][:, :, n:(n + w_length)])
        # 獲取測試資料得分
        score_this_window.append(lda.score(X_test, y_test))
    # 新增到總得分列表
    scores_windows.append(score_this_window)
 
# 設定繪圖的時間軸，時間軸上的標誌點為視窗的中間位置
w_times = (w_start + w_length / 2.) / sfreq + epochs.tmin
# 繪製模型分類結果的效能圖（得分的均值）
plt.figure()
plt.plot(w_times, np.mean(scores_windows, 0), label='Score')
plt.axvline(0, linestyle='--', color='k', label='Onset')
plt.axhline(0.5, linestyle='-', color='k', label='Chance')
plt.xlabel('time (s)')
plt.ylabel('classification accuracy')
plt.title('Classification score over time (w_length = {0}s)'.format(w_length/sfreq))
plt.legend(loc='lower right')
plt.show()

總結

• 共空間模式（CSP）是一種對兩分類任務下的空域濾波特徵提取演算法，能夠從多通道的腦機介面資料裡面提取出每一類的空間分佈成分；
• 共空間模式演算法（CSP）是利用矩陣的對角化，找到一組最優空間濾波器進行投影，使得兩類訊號的方差值差異最大化，從而得到具有較高區分度的特徵向量；
• 共空間模式演算法旨在：設計空間濾波器使得兩組腦電時空訊號矩陣濾波後，方差值差異最大化，從而得到具有較高區分度的特徵向量。用於下一步將特徵向量送入分類器進行分類。

本文參考來源：
1.運動想象| EEG訊號、共空間模式演算法（CSP）
2.CSP(Common spatial patterns)共空間模式演算法簡介
3.運動想象中共空間模式演算法（CSP）的實現