非官方中文翻译
本页译自 MNE-Python 官方文档：
查看英文原文 ↗。
文档版权归 MNE Developers 所有（BSD-3-Clause）。
说明
转到末尾下载完整示例代码。
使用主成分分析—最优基组（PCA-OBS）去除心脏伪迹#
本脚本展示如何使用 PCA-OBS 的改编版本 [1]。PCA-OBS 最初用于去除同步 EEG-fMRI 中的心冲击伪迹。这里做了调整，以消除检测到的 R 峰与心冲击伪迹之间的延迟，从而可把该算法用于去除 EEG（脑电）与 ESG（脊髓电）数据中的心脏伪迹。我们将把算法应用于 ESG 数据来说明其效果，因为在该数据中去除效果最明显。
数据集与 ESG 应用的更多细节见：https://www.biorxiv.org/content/10.1101/2024.09.05.611423v1。
# Authors: Emma Bailey <bailey@cbs.mpg.de>,
# Steinn Hauser Magnusson <hausersteinn@gmail.com>
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
import glob
import numpy as np
若尚未下载，请从 OpenNeuro 下载示例受试者数据。这会下载单个受试者单次 run 的同步 EEG 与 ESG 数据（左腕正中神经刺激）。
import openneuro
from matplotlib import pyplot as plt
import mne
from mne import Epochs, events_from_annotations
from mne.io import read_raw_eeglab
from mne.preprocessing import find_ecg_events, fix_stim_artifact
# add the path where you want the OpenNeuro data downloaded. Each run is ~2GB of data
ds = "ds004388"
target_dir = mne.datasets.default_path() / ds
run_name = "sub-001/eeg/*median_run-03_eeg*.set"
if not glob.glob(str(target_dir / run_name)):
target_dir.mkdir(exist_ok=True)
openneuro.download(dataset=ds, target_dir=target_dir, include=run_name[:-4])
block_files = glob.glob(str(target_dir / run_name))
assert len(block_files) == 1
Using default location ~/mne_data...
定义 ESG 通道（分两片贴于颈部与下背部）。
esg_chans = [
"S35",
"S24",
"S36",
"Iz",
"S17",
"S15",
"S32",
"S22",
"S19",
"S26",
"S28",
"S9",
"S13",
"S11",
"S7",
"SC1",
"S4",
"S18",
"S8",
"S31",
"SC6",
"S12",
"S16",
"S5",
"S30",
"S20",
"S34",
"S21",
"S25",
"L1",
"S29",
"S14",
"S33",
"S3",
"L4",
"S6",
"S23",
]
# Interpolation window in seconds for ESG data to remove stimulation artefact
tstart_esg = -7e-3
tmax_esg = 7e-3
# Define timing of heartbeat epochs in seconds relative to R-peaks
iv_baseline = [-400e-3, -300e-3]
iv_epoch = [-400e-3, 600e-3]
接着做最少的预处理，包括去除刺激伪迹、降采样与滤波。
raw = read_raw_eeglab(block_files[0], verbose="error")
raw.set_channel_types(dict(ECG="ecg"))
# Isolate the ESG channels (include the ECG channel for R-peak detection)
raw.pick(esg_chans + ["ECG"])
# Trim duration and downsample (from 10kHz) to improve example speed
raw.crop(0, 60).load_data().resample(2000)
# Find trigger timings to remove the stimulation artefact
events, event_dict = events_from_annotations(raw)
trigger_name = "Median - Stimulation"
fix_stim_artifact(
raw,
events=events,
event_id=event_dict[trigger_name],
tmin=tstart_esg,
tmax=tmax_esg,
mode="linear",
stim_channel=None,
)
Reading 0 ... 600000 = 0.000 ... 60.000 secs...
Used Annotations descriptions: ['Median - Stimulation']
查找 ECG 事件，并作为事件注释加入原始数据结构。
ecg_events, ch_ecg, average_pulse = find_ecg_events(raw, ch_name="ECG")
ecg_event_samples = np.asarray(
[[ecg_event[0] for ecg_event in ecg_events]]
) # Samples only
qrs_event_time = [
x / raw.info["sfreq"] for x in ecg_event_samples.reshape(-1)
] # Divide by sampling rate to make times
duration = np.repeat(0.0, len(ecg_event_samples))
description = ["qrs"] * len(ecg_event_samples)
raw.annotations.append(
qrs_event_time, duration, description, ch_names=[esg_chans] * len(qrs_event_time)
)
Using channel ECG to identify heart beats.
Setting up band-pass filter from 5 - 35 Hz
FIR 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: 5.00
- Lower transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 4.75 Hz)
- Upper passband edge: 35.00 Hz
- Upper transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 35.25 Hz)
- Filter length: 20000 samples (10.000 s)
Number of ECG events detected : 45 (average pulse 45.0 / min.)
在心脏伪迹校正前后，围绕检测到的 R 峰创建诱发响应。
events, event_ids = events_from_annotations(raw)
event_id_dict = {key: value for key, value in event_ids.items() if key == "qrs"}
epochs = Epochs(
raw,
events,
event_id=event_id_dict,
tmin=iv_epoch[0],
tmax=iv_epoch[1],
baseline=tuple(iv_baseline),
)
evoked_before = epochs.average()
# Apply function - modifies the data in place. Optionally high-pass filter
# the data before applying PCA-OBS to remove low frequency drifts
raw = mne.preprocessing.apply_pca_obs(
raw, picks=esg_chans, n_jobs=5, qrs_times=raw.times[ecg_event_samples.reshape(-1)]
)
epochs = Epochs(
raw,
events,
event_id=event_id_dict,
tmin=iv_epoch[0],
tmax=iv_epoch[1],
baseline=tuple(iv_baseline),
)
evoked_after = epochs.average()
Used Annotations descriptions: ['Median - Stimulation', 'qrs']
Not setting metadata
45 matching events found
Applying baseline correction (mode: mean)
0 projection items activated
[Parallel(n_jobs=5)]: Using backend ThreadingBackend with 5 concurrent workers.
[Parallel(n_jobs=5)]: Done 13 tasks | elapsed: 0.2s
[Parallel(n_jobs=5)]: Done 37 out of 37 | elapsed: 0.5s finished
Not setting metadata
45 matching events found
Applying baseline correction (mode: mean)
0 projection items activated
比较诱发响应，评估伪迹去除的完整性。
fig, axes = plt.subplots(1, 1, layout="constrained")
data_before = evoked_before.get_data(units=dict(eeg="uV")).T
data_after = evoked_after.get_data(units=dict(eeg="uV")).T
hs = list()
hs.append(axes.plot(epochs.times, data_before, color="k")[0])
hs.append(axes.plot(epochs.times, data_after, color="green", label="after")[0])
axes.set(ylim=[-500, 1000], ylabel="Amplitude (µV)", xlabel="Time (s)")
axes.set(title="ECG artefact removal using PCA-OBS")
axes.legend(hs, ["before", "after"])
plt.show()
参考文献#
脚本总运行时间：（0 分 3.162 秒）