说明
转到末尾下载完整示例代码。
使用 SSP 修复伪迹#
本教程介绍信号空间投影（SSP）的基础知识，说明如何将 SSP 用于伪迹修复；并通过较完整的示例演示如何用 SSP 降低环境噪声，以及修复眼动和心跳伪迹。
和往常一样，首先导入所需的 Python 模块。为了避免反复书写 mne.preprocessing，我们直接从该子模块导入几个函数：
# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
import os
import matplotlib.pyplot as plt
import numpy as np
import mne
from mne.preprocessing import (
compute_proj_ecg,
compute_proj_eog,
create_ecg_epochs,
create_eog_epochs,
)
说明
在应用 SSP（或任何伪迹修复策略）之前，请务必先观察数据中的伪迹，以确保选择了合适的修复工具。有时最好的办法就是不做处理——如果伪迹足够小，你甚至可能无需修复它们就能获得良好的分析结果。关于检测和可视化各类伪迹的指导，参见 伪迹检测概览。
什么是 SSP？#
信号空间投影（SSP）[1] 是一种通过将信号投影到低维子空间来去除 EEG 和 MEG 信号中噪声的技术。该子空间的选择方式是：计算噪声存在时各传感器上的平均模式，将该模式视为传感器空间中的一个“方向”，并构造与该噪声方向正交的子空间（关于投影的详细讲解参见 投影器与投影背景知识）。
SSP 最常见的用法是在 MEG 信号中去除来自环境源（被试身体和 MEG 系统之外的源，例如附近电气设备产生的电磁场）的噪声，且要求该噪声是平稳的（在记录期间变化不大）。不过，SSP 也可以用来去除生物伪迹，例如心跳（ECG）和眼动（EOG）伪迹。下面分别给出这些用法的示例。
示例：利用空房间记录降低环境噪声#
示例数据由 Neuromag 系统记录，该系统会在系统配置中保存用于去除环境噪声的 SSP 投影器（以便在采集过程中实时查看较为干净的原始数据）。因此，示例数据集中的所有 Raw 数据都已包含 SSP 投影器，加载数据时输出中会注明这一点：
sample_data_folder = mne.datasets.sample.data_path()
sample_data_raw_file = os.path.join(
sample_data_folder, "MEG", "sample", "sample_audvis_raw.fif"
)
# here we crop and resample just for speed
raw = mne.io.read_raw_fif(sample_data_raw_file).crop(0, 60)
raw.load_data().resample(100)
Opening raw data file /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw.fif...
Read a total of 3 projection items:
PCA-v1 (1 x 102) idle
PCA-v2 (1 x 102) idle
PCA-v3 (1 x 102) idle
Range : 25800 ... 192599 = 42.956 ... 320.670 secs
Ready.
Reading 0 ... 36037 = 0.000 ... 60.000 secs...
Finding events on: STI 014
86 events found on stim channel STI 014
Event IDs: [ 1 2 3 4 5 32]
Finding events on: STI 014
86 events found on stim channel STI 014
Event IDs: [ 1 2 3 4 5 32]
示例数据还包含一份与记录被试数据同一天采集的“空房间”（empty room）记录。与系统自带的投影器（通常在年度维护和调谐时生成）相比，这份记录能提供更准确的环境噪声估计。既然我们有这份针对具体被试的空房间记录，就用自己的记录来创建投影器，并丢弃系统提供的 SSP 投影器（先将其保存下来，以便之后与自定义投影器比较）：
system_projs = raw.info["projs"]
raw.del_proj()
empty_room_file = os.path.join(sample_data_folder, "MEG", "sample", "ernoise_raw.fif")
# cropped to 60 s just for speed
empty_room_raw = mne.io.read_raw_fif(empty_room_file).crop(0, 30)
Opening raw data file /home/circleci/mne_data/MNE-sample-data/MEG/sample/ernoise_raw.fif...
Isotrak not found
Read a total of 3 projection items:
PCA-v1 (1 x 102) idle
PCA-v2 (1 x 102) idle
PCA-v3 (1 x 102) idle
Range : 19800 ... 85867 = 32.966 ... 142.965 secs
Ready.
请注意，空房间记录本身也含有系统提供的 SSP 投影器——我们也会把它们从空房间文件中移除。
可视化空房间噪声#
下面来看看空房间噪声的频谱。可以查看每个传感器单独的频谱，也可以查看跨传感器的平均频谱（带置信带）：
raw.info["bads"] = ["MEG 2443"]
spectrum = empty_room_raw.compute_psd()
for average in (False, True):
spectrum.plot(
average=average,
dB=False,
amplitude=True,
xscale="log",
picks="data",
exclude="bads",
)
Effective window size : 3.410 (s)
Plotting amplitude spectral density (dB=False).
Plotting amplitude spectral density (dB=False).
创建空房间投影器#
我们使用 compute_proj_raw 创建 SSP 向量，并通过参数 n_grad 和 n_mag 控制投影器的数量。创建完成后，可用 plot_projs_topomap 方便地可视化投影器的场分布模式。我们传入参数 vlim='joint'，使同一通道类型的所有投影器共用同一色标范围；这样更便于比较它们之间的相对平滑程度。注意，要让函数知道投影器中包含哪些通道类型，还必须提供相应的 Info 对象：
empty_room_projs = mne.compute_proj_raw(empty_room_raw, n_grad=3, n_mag=3)
mne.viz.plot_projs_topomap(
empty_room_projs, colorbar=True, vlim="joint", info=empty_room_raw.info
)
Not setting metadata
29 matching events found
No baseline correction applied
Dropped 0/29 epochs
Adding projection: planar-Raw-0.000-30.001-PCA-01 (exp var=5.1%)
Adding projection: planar-Raw-0.000-30.001-PCA-02 (exp var=4.4%)
Adding projection: planar-Raw-0.000-30.001-PCA-03 (exp var=3.0%)
Adding projection: axial-Raw-0.000-30.001-PCA-01 (exp var=80.6%)
Adding projection: axial-Raw-0.000-30.001-PCA-02 (exp var=14.6%)
Adding projection: axial-Raw-0.000-30.001-PCA-03 (exp var=3.9%)
No channels 'eeg' found. Skipping.
请注意，基于梯度计的投影器似乎反映的是个别传感器单元的问题，而不是全局噪声源（事实上，平面梯度计对远处源的敏感度要低得多）。这正是系统提供的噪声投影器只针对磁强计计算的原因。把系统提供的投影器与针对具体被试的投影器相比较，可以看出两者相当相似（尽管顺序不同），并且左右成分的极性似乎发生了改变。
fig, axs = plt.subplots(2, 3)
for idx, _projs in enumerate([system_projs, empty_room_projs[3:]]):
mne.viz.plot_projs_topomap(
_projs, axes=axs[idx], colorbar=True, vlim="joint", info=empty_room_raw.info
)
可视化投影器对信号的影响#
我们可以用 apply_proj 把不同的投影器分别应用于 Raw 对象的不同副本，从而可视化它们对数据产生的不同影响。不过，plot 方法有一个 proj 参数，允许在绘图时临时应用投影器，因此可以借此可视化差异而无需复制数据。由于这些投影器非常相似，需要把数据放大很多才能看出差异：
mags = mne.pick_types(raw.info, meg="mag")
for title, projs in [
("system", system_projs),
("subject-specific", empty_room_projs[3:]),
]:
raw.add_proj(projs, remove_existing=True)
with mne.viz.use_browser_backend("matplotlib"):
fig = raw.plot(proj=True, order=mags, duration=1, n_channels=2)
fig.subplots_adjust(top=0.9) # make room for title
fig.suptitle(f"{title} projectors", size="xx-large", weight="bold")
3 projection items deactivated
Using matplotlib as 2D backend.
Using qt as 2D backend.
3 projection items deactivated
Using matplotlib as 2D backend.
Using qt as 2D backend.
有时在平均后的数据上更容易看出效果。这里我们使用 mne.Evoked.plot_topomap 的交互功能来开启和关闭投影器，以观察其对数据的影响。当然，这种交互功能在教程网站上无法使用，但你可以下载教程并在本地尝试：
events = mne.find_events(raw, stim_channel="STI 014")
event_id = {"auditory/left": 1}
# NOTE: appropriate rejection criteria are highly data-dependent
reject = dict(
mag=4000e-15, # 4000 fT
grad=4000e-13, # 4000 fT/cm
eeg=150e-6, # 150 µV
eog=250e-6,
) # 250 µV
# time range where we expect to see the auditory N100: 50-150 ms post-stimulus
times = np.linspace(0.05, 0.15, 5)
epochs = mne.Epochs(raw, events, event_id, proj="delayed", reject=reject)
fig = epochs.average().plot_topomap(times, proj="interactive")
Finding events on: STI 014
86 events found on stim channel STI 014
Event IDs: [ 1 2 3 4 5 32]
Not setting metadata
20 matching events found
Setting baseline interval to [-0.2, 0.0] s
Applying baseline correction (mode: mean)
Entering delayed SSP mode.
Created an SSP operator (subspace dimension = 3)
Rejecting epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 006', 'EEG 007']
此外，使用 evoked.plot() 或 evoked.plot_joint() 分别绘制应用和不应用投影器时的 ERP/F 也可能提供有用信息；使用 proj='reconstruct' 绘图同样如此，它可以减轻投影引入的信号偏差（参见下文“通过信号重建可视化 SSP 在传感器空间的偏差”）。
示例：修复 EOG 与 ECG 伪迹#
可视化伪迹#
正如 ICA 教程中所说，重要的第一步是可视化你想要修复的伪迹。下面是原始数据中的伪迹：
# pick some channels that clearly show heartbeats and blinks
regexp = r"(MEG [12][45][123]1|EEG 00.)"
artifact_picks = mne.pick_channels_regexp(raw.ch_names, regexp=regexp)
raw.plot(order=artifact_picks, n_channels=len(artifact_picks))
使用 SSP 修复 ECG 伪迹#
MNE-Python 提供了多个函数用于检测和去除 EEG、MEG 数据中的心跳。正如 伪迹检测概览 中所见，create_ecg_epochs 既可以检测心跳伪迹，又能将其提取到 Epochs 对象中，从而可以可视化心跳伪迹在各传感器上的表现：
ecg_evoked = create_ecg_epochs(raw).average()
ecg_evoked.plot_joint()
Reconstructing ECG signal from Magnetometers
Setting up band-pass filter from 8 - 16 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: 8.00
- Lower transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 7.75 Hz)
- Upper passband edge: 16.00 Hz
- Upper transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 16.25 Hz)
- Filter length: 1000 samples (10.000 s)
Number of ECG events detected : 59 (average pulse 59.0 / min.)
Not setting metadata
59 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 3)
Using data from preloaded Raw for 59 events and 101 original time points ...
0 bad epochs dropped
Created an SSP operator (subspace dimension = 3)
3 projection items activated
SSP projectors applied...
看起来 EEG 通道的分布相当分散；我们先做基线校正再重新绘图：
ecg_evoked.apply_baseline((None, None))
ecg_evoked.plot_joint()
Applying baseline correction (mode: mean)
Created an SSP operator (subspace dimension = 3)
3 projection items activated
SSP projectors applied...
要为心跳伪迹计算 SSP 投影器，可以使用 compute_proj_ecg。它以 Raw 对象为输入，返回指定数量的磁强计、梯度计和 EEG 通道投影器（默认每种通道类型两个投影器）。compute_proj_ecg 还会返回一个 events 数组，其中包含每个检测到的心跳 R 波峰值对应的采样点编号。
projs, events = compute_proj_ecg(raw, n_grad=1, n_mag=1, n_eeg=1, reject=None)
Including 3 SSP projectors from raw file
Running ECG SSP computation
Reconstructing ECG signal from Magnetometers
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: 1000 samples (10.000 s)
Number of ECG events detected : 58 (average pulse 58.0 / min.)
Computing projector
Filtering a subset of channels. The highpass and lowpass values in the measurement info will not be updated.
Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 1 - 35 Hz
FIR filter parameters
---------------------
Designing a two-pass forward and reverse, zero-phase, non-causal bandpass filter:
- Windowed frequency-domain design (firwin2) method
- Hamming window
- Lower passband edge: 1.00
- Lower transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 0.75 Hz)
- Upper passband edge: 35.00 Hz
- Upper transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 35.25 Hz)
- Filter length: 1000 samples (10.000 s)
Not setting metadata
58 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 3)
3 projection items activated
Using data from preloaded Raw for 58 events and 61 original time points ...
0 bad epochs dropped
Adding projection: planar--0.200-0.400-PCA-01 (exp var=74.5%)
Adding projection: axial--0.200-0.400-PCA-01 (exp var=95.0%)
Adding projection: eeg--0.200-0.400-PCA-01 (exp var=77.4%)
Done.
输出的第一行告诉我们，compute_proj_ecg 发现 Raw 对象中已有三个投影器，并将它们包含在返回的投影器列表中（把新的 ECG 投影器追加到列表末尾）。如果你不希望如此，可以通过布尔参数 no_proj 改变这一行为。由于我们已经运行了计算，也可以直接通过对投影器列表建立索引来分离出 ECG 投影器：
[<Projection | ECG-planar--0.200-0.400-PCA-01, active : False, n_channels : 203, exp. var : 74.55%>, <Projection | ECG-axial--0.200-0.400-PCA-01, active : False, n_channels : 102, exp. var : 94.96%>, <Projection | ECG-eeg--0.200-0.400-PCA-01, active : False, n_channels : 60, exp. var : 77.44%>]
与空房间投影器一样，我们可以可视化其头皮分布：
此外，由于这些投影器是使用专门挑选的、包含时间锁定伪迹的 epoch 创建的，我们可以把投影器及其对时间平均 epoch 的影响画成联合图。该图有三列：
左侧显示投影前（黑色）和投影后（绿色）的数据波形。可以看到，每种通道类型各一个投影器就能很好地抑制 ECG 伪迹。
中间显示与投影器相关的头皮地形图；在本例中，由于每种通道类型只有一个投影器，因此只有一幅地形图。
右侧再次显示数据波形（黑色），但这次还把这些波形投影到每种通道类型的第一个投影器上（红色），并叠加了一条作为 ECG 通道替代真值的信号（
MEG 0111）。
# ideally here we would just do `picks_trace='ecg'`, but this dataset did not
# have a dedicated ECG channel recorded, so we just pick a channel that was
# very sensitive to the artifact
fig = mne.viz.plot_projs_joint(ecg_projs, ecg_evoked, picks_trace="MEG 0111")
fig.suptitle("ECG projectors")
1 projection items deactivated
1 projection items deactivated
1 projection items deactivated
由于在 Raw 对象中未检测到专门的 ECG 传感器通道，compute_proj_ecg 默认使用磁强计来估计 ECG 信号（如上面输出第三行所述）。你也可以提供 ch_name 参数来限定用于 ECG 伪迹检测的通道；当你实际有 ECG 传感器但它在 Raw 文件中未被如此标记时，这一点最为有用。
接下来的几行输出描述了用于分离 ECG 事件的滤波器。默认设置通常已经足够，但可以通过参数 ecg_l_freq、ecg_h_freq 和 filter_length 自定义滤波器（详见 compute_proj_ecg 的文档）。
识别出 ECG 事件后，compute_proj_ecg 还会在提取每次心跳周围的 epoch 之前对数据通道进行滤波，所用的参数值由 l_freq、h_freq、filter_length、filter_method 和 iir_params 给定。同样，默认参数值通常已经足够。
默认情况下，在计算投影之前会先对滤波后的 epoch 取平均；这可以通过布尔参数 average 控制。一般来说，这会提高信噪比（这里的“信号”指的就是我们的伪迹！），因为伪迹的时间波形在各个 epoch 之间相当相似，并且与检测到的事件在时间上高度锁定。
要了解心跳如何影响每个传感器上的信号，可以分别绘制应用和不应用 ECG 投影器时的数据：
raw.del_proj()
for title, proj in [("Without", empty_room_projs), ("With", ecg_projs)]:
raw.add_proj(proj, remove_existing=False)
with mne.viz.use_browser_backend("matplotlib"):
fig = raw.plot(order=artifact_picks, n_channels=len(artifact_picks))
fig.subplots_adjust(top=0.9) # make room for title
fig.suptitle(f"{title} ECG projectors", size="xx-large", weight="bold")
6 projection items deactivated
Using matplotlib as 2D backend.
Using qt as 2D backend.
3 projection items deactivated
Using matplotlib as 2D backend.
Using qt as 2D backend.
最后请注意，上面我们向 compute_proj_ecg 函数传入了 reject=None，这意味着计算投影器时会使用所有检测到的 ECG epoch（无论这些 epoch 期间数据传感器中的信号质量如何）。默认行为是基于信号幅值剔除 epoch：EEG 通道峰峰值超过 50 µV、EOG 通道超过 250 µV、梯度计通道超过 2000 fT/cm、磁强计通道超过 3000 fT 的 epoch 会被剔除。你可以通过传入一个以 eeg、eog、mag 和 grad 为键的字典来更改这些阈值（注意阈值需以伏特、特斯拉或特斯拉/米为单位）。通常，在计算 ECG 投影器时剔除这类 epoch 是个好主意（因为通道中这些高幅波动大概属于噪声，而非反映脑活动）；上面传入 reject=None 只是为了避免数十行额外输出（逐一列出每个被剔除 epoch 的责任传感器）干扰教程阅读。
说明
compute_proj_ecg 还有一个类似的参数 flat，用于指定每种通道类型可接受的最小峰峰值幅值。
虽然 compute_proj_ecg 方便地把多项操作合并到一个函数中，MNE-Python 也提供了执行其中每个步骤的函数。具体来说：
mne.preprocessing.find_ecg_events：在Raw对象中检测心跳并返回相应的 events 数组mne.preprocessing.create_ecg_epochs：在Raw对象中检测心跳并返回Epochs对象mne.compute_proj_epochs：从任意Epochs对象创建投影器
更多细节请参阅各函数的文档。
说明
在可用于分析的电极十分有限的情况下，使用依赖空间信息的技术（如 SSP）来去除心电伪迹可能无法实现。此时，可以考虑仅需时域中心跳信息的算法，例如 mne.preprocessing.apply_pca_obs()。
使用 SSP 修复 EOG 伪迹#
再次强调，在尝试修复之前，先可视化我们的伪迹。上文我们已经看到原始数据中额区 EEG 通道的大幅偏转；下面是眼动伪迹在所有传感器上的表现：
eog_evoked = create_eog_epochs(raw).average(picks="all")
eog_evoked.apply_baseline((None, None))
eog_evoked.plot_joint()
Using EOG channel: EOG 061
EOG channel index for this subject is: [375]
Filtering the data to remove DC offset to help distinguish blinks from saccades
Selecting channel EOG 061 for blink detection
Setting up band-pass filter from 1 - 10 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: 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: 1000 samples (10.000 s)
Now detecting blinks and generating corresponding events
Found 10 significant peaks
Number of EOG events detected: 10
Not setting metadata
10 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 9)
Using data from preloaded Raw for 10 events and 101 original time points ...
0 bad epochs dropped
Applying baseline correction (mode: mean)
Created an SSP operator (subspace dimension = 9)
9 projection items activated
SSP projectors applied...
与处理心跳伪迹时一样，我们可以使用 compute_proj_eog 为眼动伪迹计算 SSP 投影器。它同样以 Raw 对象为输入，返回指定数量的磁强计、梯度计和 EEG 通道投影器（默认每种通道类型两个投影器）。这次我们传入 no_proj 参数（这样只返回新的 EOG 投影器，而不包含 Raw 对象中已有的投影器），并通过将其赋给 _ 来忽略返回的 events 数组（这是在 Python 中处理不需要的返回值元素的惯用做法）。
eog_projs, _ = compute_proj_eog(
raw, n_grad=1, n_mag=1, n_eeg=1, reject=None, no_proj=True
)
Running EOG SSP computation
Using EOG channel: EOG 061
EOG channel index for this subject is: [375]
Filtering the data to remove DC offset to help distinguish blinks from saccades
Selecting channel EOG 061 for blink detection
Setting up band-pass filter from 1 - 10 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: 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: 1000 samples (10.000 s)
Now detecting blinks and generating corresponding events
Found 10 significant peaks
Number of EOG events detected: 10
Computing projector
Filtering a subset of channels. The highpass and lowpass values in the measurement info will not be updated.
Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 1 - 35 Hz
FIR filter parameters
---------------------
Designing a two-pass forward and reverse, zero-phase, non-causal bandpass filter:
- Windowed frequency-domain design (firwin2) method
- Hamming window
- Lower passband edge: 1.00
- Lower transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 0.75 Hz)
- Upper passband edge: 35.00 Hz
- Upper transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 35.25 Hz)
- Filter length: 1000 samples (10.000 s)
Not setting metadata
10 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 9)
9 projection items activated
Using data from preloaded Raw for 10 events and 41 original time points ...
0 bad epochs dropped
Adding projection: planar--0.200-0.200-PCA-01 (exp var=85.8%)
Adding projection: axial--0.200-0.200-PCA-01 (exp var=86.3%)
Adding projection: eeg--0.200-0.200-PCA-01 (exp var=98.8%)
Done.
与空房间投影器和 ECG 投影器一样，我们可以可视化其头皮分布：
还可以做联合图：
fig = mne.viz.plot_projs_joint(eog_projs, eog_evoked, "eog")
fig.suptitle("EOG projectors")
1 projection items deactivated
1 projection items deactivated
1 projection items deactivated
最后，我们可以与 EOG 的 Evoked 一起做联合可视化。这里我们还会做一个糟糕的选择：为 EEG 和磁强计各选两个 EOG 投影器，然后会看到它们在图里表现为噪声。尽管投影后的时间进程（左列）看起来也许还行，但问题会出现在中间（地形图）和右侧（把通道数据投影到投影向量上）的图中：
第二幅磁强计地形图呈现出双侧听觉场模式。
统一缩放后的投影时间进程（实线）显示，虽然第一个投影器对应的波形（红色）具有较大的类 EOG 幅值，但第二个投影器的波形（蓝绿色）要小得多。
重新归一化后的投影时间进程显示，相对于 EOG 通道数据（黄色），第二个 PCA 波形噪声很大。
eog_projs_bad, _ = compute_proj_eog(
raw, n_grad=1, n_mag=2, n_eeg=2, reject=None, no_proj=True
)
fig = mne.viz.plot_projs_joint(eog_projs_bad, eog_evoked, picks_trace="eog")
fig.suptitle("Too many EOG projectors")
Running EOG SSP computation
Using EOG channel: EOG 061
EOG channel index for this subject is: [375]
Filtering the data to remove DC offset to help distinguish blinks from saccades
Selecting channel EOG 061 for blink detection
Setting up band-pass filter from 1 - 10 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: 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: 1000 samples (10.000 s)
Now detecting blinks and generating corresponding events
Found 10 significant peaks
Number of EOG events detected: 10
Computing projector
Filtering a subset of channels. The highpass and lowpass values in the measurement info will not be updated.
Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 1 - 35 Hz
FIR filter parameters
---------------------
Designing a two-pass forward and reverse, zero-phase, non-causal bandpass filter:
- Windowed frequency-domain design (firwin2) method
- Hamming window
- Lower passband edge: 1.00
- Lower transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 0.75 Hz)
- Upper passband edge: 35.00 Hz
- Upper transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 35.25 Hz)
- Filter length: 1000 samples (10.000 s)
Not setting metadata
10 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 9)
9 projection items activated
Using data from preloaded Raw for 10 events and 41 original time points ...
0 bad epochs dropped
Adding projection: planar--0.200-0.200-PCA-01 (exp var=85.8%)
Adding projection: axial--0.200-0.200-PCA-01 (exp var=86.3%)
Adding projection: axial--0.200-0.200-PCA-02 (exp var=3.3%)
Adding projection: eeg--0.200-0.200-PCA-01 (exp var=98.8%)
Adding projection: eeg--0.200-0.200-PCA-02 (exp var=0.6%)
Done.
1 projection items deactivated
2 projection items deactivated
2 projection items deactivated
现在我们重复上面的绘图（含空房间和 ECG 投影器），并将其与同时包含空房间、ECG 和 EOG 投影器的图比较，看看眼动伪迹修复得如何：
for title in ("Without", "With"):
if title == "With":
raw.add_proj(eog_projs)
with mne.viz.use_browser_backend("matplotlib"):
fig = raw.plot(order=artifact_picks, n_channels=len(artifact_picks))
fig.subplots_adjust(top=0.9) # make room for title
fig.suptitle(f"{title} EOG projectors", size="xx-large", weight="bold")
Using matplotlib as 2D backend.
Using qt as 2D backend.
3 projection items deactivated
Using matplotlib as 2D backend.
Using qt as 2D backend.
请注意，前两个磁强计通道（MEG 1411 和 MEG 1421）中与 EEG 大幅偏转同时出现的小峰也已被去除。
选择投影器的数量#
在上面的示例中，我们用了 3 个投影器（全部针对磁强计）来捕捉空房间噪声，并看到为梯度计计算的投影器未能捕捉全局模式（因此我们丢弃了梯度计投影器）。接着我们计算了 3 个投影器（每种通道类型 1 个）来捕捉心跳伪迹，又用 3 个来捕捉眼动伪迹。我们是怎样确定这些数量的？简短的回答是“凭经验”——了解心跳伪迹在传感器阵列上通常如何表现，能让我们在看到它时识别出来，也能识别出额外的投影器何时捕捉的是心跳伪迹之外的其他东西（因而可能正在去除脑信号，应当丢弃）。
通过信号重建可视化 SSP 在传感器空间的偏差#
由于 SSP 执行的是正交投影，数据中任何与 SSP 空间方向不完全正交的空间成分，其整体幅值都会被投影操作降低。换句话说，SSP 通常会在传感器空间数据中引入一定程度的幅值降低偏差。
在对 M/EEG 数据进行源定位时，这些投影会被正确地考虑在内——不仅应用于 M/EEG 数据，也应用于正向解，因此 SSP 不应使估计的源幅值产生偏差。然而，对传感器空间分析而言，可视化 SSP 投影使数据产生偏差的程度可能很有用。可以在 evoked 绘图函数中使用 proj='reconstruct' 来考察这一点，例如通过 evoked.plot()，这里为加快速度只限于 EEG 通道：
evoked_eeg = epochs.average().pick("eeg")
evoked_eeg.del_proj().add_proj(ecg_projs).add_proj(eog_projs)
fig, axes = plt.subplots(
1, 3, figsize=(8, 3), sharex=True, sharey=True, layout="constrained"
)
for pi, proj in enumerate((False, True, "reconstruct")):
ax = axes[pi]
evoked_eeg.plot(proj=proj, axes=ax, spatial_colors=True)
parts = ax.get_title().split("(")
ylabel = (
f"{parts[0]} ({ax.get_ylabel()})\n{parts[1].replace(')', '')}"
if pi == 0
else ""
)
ax.set(ylabel=ylabel, title=f"proj={proj}")
ax.yaxis.set_tick_params(labelbottom=True)
for text in list(ax.texts):
text.remove()
Rejecting epoch based on EEG : ['EEG 001', 'EEG 002', 'EEG 003', 'EEG 006', 'EEG 007']
3 projection items deactivated
3 projection items deactivated
Created an SSP operator (subspace dimension = 2)
6 projection items activated
SSP projectors applied...
Automatic origin fit: head of radius 91.2 mm
Computing dot products for 60 EEG channels...
Computing cross products for 60 → 60 EEG channels...
Preparing the mapping matrix...
Truncating at 58/60 components and regularizing with α=1.0e-01
当重建仅用于可视化时，使用 proj="reconstruct" 很方便。若要获得重建后的传感器空间数据以用于进一步处理，应直接使用 reconstruct_proj()。它还可以把操作限制在选定的投影器上，这在比较不同伪迹校正的效果时很有用。这里 evoked_eeg 同时包含 ECG 和 EOG 投影器，但我们只应用并重建 EOG 投影器：
evoked_eog_reconstructed = evoked_eeg.copy().reconstruct_proj(projs=eog_projs)
evoked_eog_reconstructed.plot(spatial_colors=True)
Created an SSP operator (subspace dimension = 1)
1 projection items activated
SSP projectors applied...
Automatic origin fit: head of radius 91.2 mm
Computing dot products for 60 EEG channels...
Computing cross products for 60 → 60 EEG channels...
Preparing the mapping matrix...
Truncating at 59/60 components and regularizing with α=1.0e-01
注意，这里 EEG 和磁强计通道的偏差因重建而减小。这表明应用 SSP 确实略微降低了信号在传感器空间中的幅值，但不应使源空间中的幅值产生偏差。
参考文献#
脚本总运行时间：（0 分钟 18.242 秒）