關注微信公眾號：腦機接口研習社
了解腦機接口最近進展

系列文章目錄

Python專欄 | 腦電圖和腦磁圖(EEG/MEG)的數據分析方法之載入數據

Python專欄 | MNE腦電數據(EEG/MEG)可視化

Python專欄 | MNE數據預處理方法——獨立成分分析

Python專欄 | 獨立成分分析(ICA)的實例應用：消除偽影信號

持續更新中……

---

文章目錄

• 系列文章目錄
• 識別偽影信號的方法
• • 1. Selecting ICA components manually （手動識別）
  • 2. Using an EOG channel to select ICA components
• 總結

---

識別偽影信號的方法

1. Selecting ICA components manually （手動識別）

確定要排除的Components后，可以通過設置ica.exclude屬性手動指定。

僅設置ica.exclude不會立即執行任何操作，只會將排除的Independent Components (IC)添加到列表中，以便之后在需要時使用。

一旦設置了排除項，即使未傳遞exclude參數，比如plot_overlay的ICA方法，也會排除指定的Components。

使用mne.preprocessing.ICA.save和mne.preprocessing.read_ica時可以保留被排除的Components列表。

代碼示例：

ica.exclude = [0, 1] #indices chosen based on various plots above

現在我們已經設置好了排除項（ica.exclude），接著我們可以使用ica.apply的方法重塑已經去除偽影的傳感器信號。

即現在進入了下圖中的ica.apply環節。

▲點擊圖片可查看大圖

我們將把原始數據與重塑的數據一起繪制，來表明心跳和眨眼偽影已得到修復。

代碼示例：

#ica.apply() changes the Raw object in-place, so let's make a copy first: raw.load_data() reconst_raw = raw.copy() ica.apply(reconst_raw)#原始數據 raw.plot(order=artifact_picks, n_channels=len(artifact_picks),show_scrollbars=False) #重塑數據，和原始數據進行對比，來表明心跳和眨眼偽影已得到修復。 reconst_raw.plot(order=artifact_picks, n_channels=len(artifact_picks),show_scrollbars=False) del reconst_raw #刪除變量reconst_raw

輸出結果：

Applying ICA to Raw instanceTransforming to ICA space (15 components)Zeroing out 2 ICA componentsProjecting back using 364 PCA components

2. Using an EOG channel to select ICA components

通過上面手動選擇要排除的IC的方法似乎比較簡單，但是當需要處理成百上千個被試的數據時，手動排除的方式就略顯繁瑣。

在這種情況下，比較合理的做法是：首先使用專用的EOG或ECG傳感器進行實驗，把由此測得的信號作為一種“模型（pattern）”，然后將“模型”與每一種IC相對應，利用算法自動標記出和“EOG / ECG模型信號”匹配的所有ICs。

在mne-python包里，我們將使用find_bads_eog自動查找與EOG信號最匹配的ICs，然后再使用plot_scores以及其他一些繪圖功能來查看選擇了哪些IC。

由于在手動選擇步驟里，我們使用了ica.exclude = [0, 1]，手動選擇了ICA001和ICA000，因此在使用自動選擇功能函數前，我們首先需要將ica.exclude重置為空列表。

代碼示例：

#自動選擇與ECG和EOG信號對應的ICs ica.exclude = []#首先將ica.exclude變成空列表#find which ICs match the EOG pattern eog_indices, eog_scores = ica.find_bads_eog(raw) ica.exclude = eog_indices# barplot of ICA component "EOG match" scores ica.plot_scores(eog_scores) #eog_scores是計算每個ICs（本例子里共有15個ICs）與EOG信號的match scores。絕對值越接近1代表匹配度越高，越接近0代表匹配度越低。 #由畫出的條形圖可以看出，ICA000的值是0.94763376，匹配度最高。# plot diagnostics ica.plot_properties(raw, picks=eog_indices)# plot ICs applied to raw data, with EOG matches highlighted ica.plot_sources(raw, show_scrollbars=False)# plot ICs applied to the averaged EOG epochs, with EOG matches highlighted ica.plot_sources(eog_evoked)

輸出結果：

filtering target Setting up band-pass filter from 1 - 10 HzFIR filter parameters --------------------- Designing a two-pass forward and reverse, zero-phase, non-causal bandpass filter: - Windowed frequency-domain design (firwin2) method - Hann window - Lower passband edge: 1.00 - Lower transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 0.75 Hz) - Upper passband edge: 10.00 Hz - Upper transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 10.25 Hz) - Filter length: 6007 samples (10.001 sec)C:\ProgramData\Anaconda3\lib\site-packages\mkl_fft\_numpy_fft.py:331: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.output = mkl_fft.rfft_numpy(a, n=n, axis=axis)Using multitaper spectrum estimation with 7 DPSS windows Not setting metadata Not setting metadata 138 matching events found No baseline correction applied 0 projection items activated 0 bad epochs dropped C:\ProgramData\Anaconda3\lib\site-packages\mkl_fft\_numpy_fft.py:331: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.output = mkl_fft.rfft_numpy(a, n=n, axis=axis)Creating RawArray with float64 data, n_channels=16, n_times=166800Range : 25800 ... 192599 = 42.956 ... 320.670 secs Ready.

以上我們在原始Raw實例以及提取了EOG偽影的Evoked實例上都使用了plot_sources。這可以是另一種確認find_bads_eog已標識正確component的方法。

總結

今天用到的代碼總結：

#手動選擇要排除的ICs ica.exclude = [0, 1] # indices chosen based on various plots above#現在我們已經設置好了排除項（ica.exclude） #接著我們可以使用ica.apply的方法重塑已經去除偽影的傳感器信號。 #ica.apply() changes the Raw object in-place, so let's make a copy first: raw.load_data() reconst_raw = raw.copy() ica.apply(reconst_raw)#原始數據 raw.plot(order=artifact_picks, n_channels=len(artifact_picks),show_scrollbars=False) #重塑數據，和原始數據進行對比，來表明心跳和眨眼偽影已得到修復。 reconst_raw.plot(order=artifact_picks, n_channels=len(artifact_picks),show_scrollbars=False) del reconst_raw #刪除變量reconst_raw#自動選擇與ECG和EOG信號對應的ICs ica.exclude = []#首先將ica.exclude變成空列表#find which ICs match the EOG pattern eog_indices, eog_scores = ica.find_bads_eog(raw) ica.exclude = eog_indices# barplot of ICA component "EOG match" scores ica.plot_scores(eog_scores) #eog_scores是計算每個ICs（本例子里共有15個ICs）與EOG信號的match scores。絕對值越接近1代表匹配度越高，越接近0代表匹配度越低。 #由畫出的條形圖可以看出，ICA000的值是0.94763376，匹配度最高。# plot diagnostics ica.plot_properties(raw, picks=eog_indices)# plot ICs applied to raw data, with EOG matches highlighted ica.plot_sources(raw, show_scrollbars=False)# plot ICs applied to the averaged EOG epochs, with EOG matches highlighted ica.plot_sources(eog_evoked)

參考鏈接

https://mne.tools/stable/auto_tutorials/preprocessing/plot_40_artifact_correction_ica.html#fitting-and-plotting-the-ica-solution

未完待續……

欲知后事如何，請關注我們的公眾號：腦機接口研習社

總結

以上是生活随笔為你收集整理的Python专栏 | ICA应用：如何识别伪影信号？（一）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。