關(guān)注微信公眾號：腦機(jī)接口研習(xí)社
了解腦機(jī)接口最近進(jìn)展

系列文章目錄

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

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

Python專欄 | MNE數(shù)據(jù)預(yù)處理方法——獨(dú)立成分分析

Python專欄 | 獨(dú)立成分分析(ICA)的實(shí)例應(yīng)用：消除偽影信號

持續(xù)更新中……

---

文章目錄

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

---

識別偽影信號的方法

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

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

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

一旦設(shè)置了排除項(xiàng)，即使未傳遞exclude參數(shù)，比如plot_overlay的ICA方法，也會排除指定的Components。

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

代碼示例：

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

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

即現(xiàn)在進(jìn)入了下圖中的ica.apply環(huán)節(jié)。

▲點(diǎn)擊圖片可查看大圖

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

代碼示例：

#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)#原始數(shù)據(jù) raw.plot(order=artifact_picks, n_channels=len(artifact_picks),show_scrollbars=False) #重塑數(shù)據(jù)，和原始數(shù)據(jù)進(jìn)行對比，來表明心跳和眨眼偽影已得到修復(fù)。 reconst_raw.plot(order=artifact_picks, n_channels=len(artifact_picks),show_scrollbars=False) del reconst_raw #刪除變量reconst_raw

輸出結(jié)果：

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的方法似乎比較簡單，但是當(dāng)需要處理成百上千個被試的數(shù)據(jù)時，手動排除的方式就略顯繁瑣。

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

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

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

代碼示例：

#自動選擇與ECG和EOG信號對應(yīng)的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是計(jì)算每個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)

輸出結(jié)果：

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實(shí)例以及提取了EOG偽影的Evoked實(shí)例上都使用了plot_sources。這可以是另一種確認(rèn)find_bads_eog已標(biāo)識正確component的方法。

總結(jié)

今天用到的代碼總結(jié)：

#手動選擇要排除的ICs ica.exclude = [0, 1] # indices chosen based on various plots above#現(xiàn)在我們已經(jīng)設(shè)置好了排除項(xiàng)（ica.exclude） #接著我們可以使用ica.apply的方法重塑已經(jīng)去除偽影的傳感器信號。 #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)#原始數(shù)據(jù) raw.plot(order=artifact_picks, n_channels=len(artifact_picks),show_scrollbars=False) #重塑數(shù)據(jù)，和原始數(shù)據(jù)進(jìn)行對比，來表明心跳和眨眼偽影已得到修復(fù)。 reconst_raw.plot(order=artifact_picks, n_channels=len(artifact_picks),show_scrollbars=False) del reconst_raw #刪除變量reconst_raw#自動選擇與ECG和EOG信號對應(yīng)的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是計(jì)算每個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

未完待續(xù)……

欲知后事如何，請關(guān)注我們的公眾號：腦機(jī)接口研習(xí)社

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。