python - 使用 ICA 去除 EEG 信号中的眨眼

Question

我是 scikit-learn 的新手，但我正在尝试消除单个 EEG 通道内的眨眼（噪声峰值）。我在互联网上搜索过，但只看到使用 MNE、PyEEG 或其他 Python 模块的更复杂的读数。我只想要一些简单且仅依赖于 sklearn 的东西。这是我的代码：

#The channel containing some eye-blinks
X = f1ep1_data[:,[4]]

#run ICA on signal
ica = FastICA(n_components=2)
ica.fit(X)

#reconstruct signal with independent components
components = ica.fit_transform(X)
X_restored = ica.inverse_transform(components)

fig1 = plt.figure()
plt.subplot(3,1,1)
plt.title("Original signal")
plt.plot(f1ep1_timescale, X)

plt.subplot(3,1,2)
plt.title("Components")
plt.plot(f1ep1_timescale, components)

plt.subplot(3,1,3)
plt.title("Signal Reconstructed")
plt.plot(f1ep1_timescale, X_restored)
plt.draw()

我所期待的是分离成两个部分，一个更清晰的 EEG 信号和眨眼。我无法弄清楚是什么问题。有人可以伸出援助之手吗？

标签： pythonscikit-learntime-series

Answer 1

您是否注意到您的“组件”正是按比例缩放和倒置的原始信号？那是因为你不能得到比信号更多的组件。

您需要执行以下步骤：

1. 将所有EEG 通道输入 ICA
2. 手动移除包含眨眼或其他伪影的组件
3. 使用逆变换重建

让我们详细看一下第2步：为什么要手动移除组件？ ICA 对眨眼一无所知。它根据统计测量将信号分离成分量。如果幸运的话，其中一些组件看起来像眨眼。

到目前为止还好，但真正的问题是组件的顺序没有定义。运行 ICA，您可能会发现组件 1 包含眨眼。再次运行它，它们在组件 3 中。再次，它们在组件 2 和 5 中......

无法提前知道要删除哪些组件以及要删除多少组件。这就是为什么每次运行时都需要手动告诉算法的原因。

在看起来像这样的代码中：

# Use all channels - they will contain eye blinks to varying degrees
X = f1ep1_data[:, :]

# run ICA on signal
ica = FastICA(n_components=x.shape[1])  # we want *all* the components
ica.fit(X)

# decompose signal into components
components = ica.fit_transform(X)

# plot components and ask user which components to remove
# ...
remove_indices = [0, 1, 3]  # pretend the user selected components 0, 1, and 3

# "remove" unwanted components by setting them to 0 - simplistic but gets the job done
components[:, remove_indices] = 0

#reconstruct signal
X_restored = ica.inverse_transform(components)

您可能对手动步骤不满意。运气不好 :) 2013 年没有强大的自动算法可以标记眨眼成分。我认为这没有改变，但如果有什么东西，你会发现它是特定领域的包之一，如 MNE 或 PyEEG。

---

不过，如果您碰巧有 EOG 记录，那还是有希望的！EEG 记录中存在EOG 伪影的全自动校正方法。该方法基于典型相关或回归（我不记得细节），但您需要将 EOG 信号与 EEG 一起记录。

---

我用模拟的“EEG”数据创建了一个工作示例。数据由三个通道组成：额叶、中央和顶叶。10 Hz 的 alpha 活动在后部最强，一些类似眨眼的尖峰在前部最强。

希望这个例子能更好地说明如何从多通道数据中删除组件。

import numpy as np
import scipy.signal as sps
from sklearn.decomposition import FastICA
import matplotlib.pyplot as plt

np.random.seed(42)

n = 1000
fs = 100
noise = 3

# simulate EEG data with eye blinks
t = np.arange(n)
alpha = np.abs(np.sin(10 * t / fs)) - 0.5
alpha[n//2:] = 0
blink = np.zeros(n)
blink[n//2::200] += -1
blink = sps.lfilter(*sps.butter(2, [1*2/fs, 10*2/fs], 'bandpass'), blink)

frontal = blink * 200 + alpha * 10 + np.random.randn(n) * noise
central = blink * 100 + alpha * 15 + np.random.randn(n) * noise
parietal = blink * 10 + alpha * 25 + np.random.randn(n) * noise

eeg = np.stack([frontal, central, parietal]).T  # shape = (100, 3)

# plot original data
plt.subplot(3, 1, 1)
plt.plot(frontal + 50)
plt.plot(central + 100)
plt.plot(parietal + 150)
plt.yticks([50, 100, 150], ['Fz', 'Cz', 'Pz'])
plt.ylabel('original data')

# decompose EEG and plot components
ica = FastICA(n_components=3)
ica.fit(eeg)
components = ica.transform(eeg)

plt.subplot(3, 1, 2)
plt.plot([[np.nan, np.nan, np.nan]])  # advance the color cycler to give the components a different color :)
plt.plot(components + [0.5, 1.0, 1.5])
plt.yticks([0.5, 1.0, 1.5], ['0', '1', '2'])
plt.ylabel('components')

# looks like component #2 (brown) contains the eye blinks
# let's remove them (hard coded)!
components[:, 2] = 0

# reconstruct EEG without blinks
restored = ica.inverse_transform(components)

plt.subplot(3, 1, 3)
plt.plot(restored + [50, 100, 150])
plt.yticks([50, 100, 150], ['Fz', 'Cz', 'Pz'])
plt.ylabel('cleaned data')