Python-EEG工具庫MNE中文教程(11)-信號空間投影SSP 應用

@[toc]

本分享爲腦機學習者Rose整理髮表於公衆號：腦機接口社區(微信號：Brain_Computer).QQ交流羣：903290195

信號空間投影（SSP）

在前面一篇分享信號空間投影SSP數學原理中提到，投影矩陣將根據您試圖投射出的噪聲種類而變化。信號空間投影（SSP）是一種經過比較有無感興趣信號的測量值來估算投影矩陣應該是什麼的方法。例如，您能夠進行其餘「空房間」測量，以記錄沒有對象存在時傳感器上的活動。經過查看空房間測量中各MEG傳感器的活動空間模式，能夠建立一個或多個N維向量，以給出傳感器空間中環境噪聲的「方向」（相似於上面示例中「觸發器的影響」的向量）。SSP一般也用於消除心跳和眼睛運動僞影，在用於消除心跳和眼睛運動僞影的案例中，就不是經過空房間錄製，而是經過檢測僞影，提取僞影周圍的時間段(epochs)並求平均值來估計噪聲的方向。有關示例，請參見使用SSP修復工件。

一旦知道了噪聲向量，就能夠建立一個與其正交的超平面，並構造一個投影矩陣，將實驗記錄投影到該超平面上。這樣，測量中與環境噪聲相關的部分就能夠被移除。一樣，應該清楚的是，投影下降了數據的維數-你仍然會有相同數量的傳感器信號，但它們不會都是線性獨立的-但一般有數十或數百個傳感器，而你要消除的噪聲子空間只有3-5維，所以自由度的損失一般是沒有問題的。

MNE Python中的投影(projector)

在示例數據中，已經使用空房間記錄執行了SSP，可是投影與原始數據一塊兒存儲，而且還沒有應用（或者說，投影還沒有激活）。
在這裏，我將加載示例數據並將其裁剪爲60秒；

能夠在如下read_raw_fif（）的輸出中看到投影：

1.導入工具庫

import os
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D  # noqa
from scipy.linalg import svd
import mne

2.加載數據

sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = os.path.join(sample_data_folder, 'MEG', 'sample',
                                    'sample_audvis_raw.fif')
raw = mne.io.read_raw_fif(sample_data_raw_file)
raw.crop(tmax=60).load_data()

在MNE-Python中，環境噪聲矢量是經過主成分分析（一般縮寫爲PCA）來計算的，這就是爲何SSP投影儀一般有「PCA-v1」這樣的名稱。（順便說一句，因爲執行主成分分析的過程在幕後使用了奇異值分解，所以在已發表的論文中也常常看到相似「投影儀是使用SVD計算的」這樣的短語。）投影儀存儲在raw.info的projs字段中：

在MNE-Python中，使用主成分分析(一般縮寫爲"PCA")來計算環境噪聲向量，這就是SSP投影一般使用"PCA-v1"之類的名稱的緣由。 (順便說一下，因爲執行PCA的過程在後臺使用了奇異值分解，因此在已發表的論文中常常會看到相似"使用SVD計算投影"之類的短語。） 投影(projector)存儲在raw.info的projs字段中：

print(raw.info['projs'])

[<Projection | PCA-v1, active : False, n_channels : 102>, <Projection | PCA-v2, active : False, n_channels : 102>, <Projection | PCA-v3, active : False, n_channels : 102>]

raw.info['projs']是投影對象的普通Python列表，能夠經過索引來訪問各個投影。投影對象(Projection object)自己相似於Python dict，所以可使用其.keys()方法查看它包含哪些字段(一般不須要直接訪問其屬性，但若有必要，能夠這樣作)：

first_projector = raw.info['projs'][0]
print(first_projector)
print(first_projector.keys())

<Projection | PCA-v1, active : False, n_channels : 102> dict_keys(['kind', 'active', 'desc', 'data', 'explained_var'])

Raw，Epoch和Evoked對象都有一個布爾類型的 proj屬性，該屬性指示對象中是否存儲有任何未應用/不活動的投影。換句話說，若是至少有一個投影而且全部投影都處於活動狀態，則proj屬性爲True。此外，每一個投影還具備一個布爾活動字段：

print(raw.proj)
print(first_projector['active'])

False False

3.計算投影

在MNE Python中，SSP向量可使用如下通用函數計算： mne.compute_proj_raw； mne.compute_proj_epochs； mne.compute_proj_evoked。 這些函數所作的通常假設是，傳遞的數據包含要經過投影修復的工件的原始數據、時間段或平均值。 在實踐中，這一般涉及空房間記錄或平均ECG或EOG僞影的連續原始數據。

"""
經過比較使用和不使用投影的曲線圖，能夠看到投影儀對測量信號的影響。
默認狀況下，`raw.plot（）`將在繪圖前在後臺應用投影儀(不修改：class：`mne.io.Raw`對象);
能夠經過以下所示的布爾``proj``參數來控制它，
也能夠經過繪圖窗口右下角的：kbd：`Proj`按鈕訪問投影界面，
以交互方式打開和關閉它們.
這裏咱們只看磁力計，還有一個文件開頭的2秒樣本。
"""
mags = raw.copy().crop(tmax=2).pick_types(meg='mag')
for proj in (False, True):
    fig = mags.plot(butterfly=True, proj=proj)
    fig.subplots_adjust(top=0.9)
    fig.suptitle('proj={}'.format(proj), size='xx-large', weight='bold')

如上圖，未進行投影的數據proj=Flase的效果展現，點擊"proj"紅框彈出SSP projectors vectors能夠發現都未選中激活；進行了投影的數據圖爲proj=True點擊"proj"紅框能夠發現投影的信息。

4.加載和保存投影

SSP除了能夠減小環境噪聲外，還能夠用於其餘類型的信號清洗。能夠發如今上一個圖中的磁力計信號中有兩個較大的偏移，這些偏移沒有被空房間的投影消除，這是受試者心跳的僞影。SSP也能夠用於移除這些工件。樣本數據包括用於下降心跳噪聲的投影，這些投影與原始數據保存在單獨的文件中，可使用mne.read_proj()函數加載該文件：

ecg_proj_file = os.path.join(sample_data_folder, 'MEG', 'sample',
                             'sample_audvis_ecg-proj.fif')
ecg_projs = mne.read_proj(ecg_proj_file)
print(ecg_projs)

"""
利用mne.write_proj()函數，可用於將投影數據以.fif格式保存到磁盤：

MNE Python推薦使用以-proj.fif(或-proj.fif.gz)來保存投影數據
"""
mne.write_proj('heartbeat-proj.fif', ecg_projs)

5.添加和移除投影

上面，當咱們打印從文件加載的ecg_projs列表時，它顯示了兩臺用於梯度計的投影（前兩臺，標爲"planar"），兩臺用於磁力計的投影（中兩臺，標爲"axial"），兩臺用於EEG 傳感器（最後兩個，標記爲"eeg"）。咱們可使用add_proj()方法將它們添加到Raw對象：

raw.add_proj(ecg_projs)

要刪除投影，能夠利用del_proj()方法，它是根據raw.info['projs']列表中的索引刪除投影。 若是想要用新的投影替換現有投影，可使用raw.add_proj(ecg_projs,remove_existing = True)來實現。

想要了解心電圖(ECG)投影儀如何影響測量的信號，咱們能夠再次對使用投影和不使用投影的數據進行繪圖（注：plot()方法只能臨時應用投影進行可視化，而不會永久更改基礎數據）。咱們將上面建立的mags變量（只有空房間SSP投影）與空房間和ECG投影儀的數據進行比較：

mags_ecg = raw.copy().crop(tmax=2).pick_types(meg='mag')
for data, title in zip([mags, mags_ecg], ['Without', 'With']):
    fig = data.plot(butterfly=True, proj=True)
    fig.subplots_adjust(top=0.9)
    fig.suptitle('{} ECG projector'.format(title), size='xx-large',
                 weight='bold')

在without ECG projector中，meg數據中的ECG部分沒有進行projector，而在with ECG projector中，meg數據中ECG部分進行了projector，結果要平滑一些。

參考 Uusitalo MA and Ilmoniemi RJ. (1997). Signal-space projection method for separating MEG or EEG into components. Med Biol Eng Comput 35(2), 135–140. doi:10.1007/BF02534144

Python-EEG工具庫MNE中文教程(11)-信號空間投影SSP 應用

腦機學習者Rose筆記分享,QQ交流羣：903290195 更多分享，請關注公衆號