步骤 18:剔心电成分

import mne

# ====== 方式 A：自动识别（推荐新手，默认 ctps） ======
# 注意：method 仅接受 'ctps'（默认）/'correlation'，
#       'variance' 不是合法值（从未存在过）
ecg_idx, ecg_scores = ica.find_bads_ecg(
    raw,
    ch_name=None,             # None = 自动选 ecg 通道
    method='ctps',            # cross-trial phase statistics（默认）
    threshold='auto',         # 'auto'（默认）或 float
    # 注意：threshold 对两种 method 都生效：
    # - 'ctps': 表示 Kuiper 检验的显著性水平；'auto' 会基于采样频率自动算
    # - 'correlation': 表示 z-score 阈值
    l_freq=8,                 # 真实默认值（注意是 8 不是 8.0，Hz）
    h_freq=16,                # 真实默认值（不是 30）
)

print(f"自动识别的心电成分: {ecg_idx}")
print(f"各成分分数: {ecg_scores}")

# ====== 方式 B：基于简单相关（最快） ======
ecg_idx, ecg_scores = ica.find_bads_ecg(
    raw,
    ch_name='ECG001',         # 显式指定 ECG 通道
    method='correlation',
    threshold=2.0,
    l_freq=8,
    h_freq=16,
)

# ====== 人工核对（必做） ======
ica.plot_scores(ecg_scores, exclude=ecg_idx)
ica.plot_components(picks=ecg_idx, inst=raw)
ica.plot_sources(inst=raw, picks=ecg_idx, start=0, stop=30)
# 期望：成分时间序列呈现明显的 QRS 周期节律

# ====== 合并到 exclude 列表 ======
ica.exclude.extend(ecg_idx)
print(f"当前排除列表（含 EOG + ECG）: {ica.exclude}")

# 验证 1：候选成分数量
print(f"自动检出 ECG 成分数: {len(ecg_idx)}")
assert 0 <= len(ecg_idx) <= 3, "ECG 成分数异常"

# 验证 2：topomap 是否合理
# ECG 成分常见模式：单极强（V-lead 类似）或扩散型
for idx in ecg_idx:
    ica.plot_components(picks=[idx], inst=raw)

# 验证 3：时间序列显示 QRS 节律
ica.plot_sources(inst=raw, picks=ecg_idx, start=0, stop=60)
# 期望：清晰可见周期性心跳尖峰

# 验证 4：与 ECG 通道相关性
import numpy as np
ecg_data = raw.get_data(picks='ecg')[0]
ic_data = ica.get_sources(raw).get_data()[ecg_idx[0]]
corr = np.corrcoef(ecg_data, ic_data)[0, 1]
print(f"IC{ecg_idx[0]} vs ECG 相关系数: {corr:.2f}")
# 期望：|r| > 0.2（ECG 信号较弱时）