步骤 11:坏导检测(第二次)

import mne
import numpy as np

# ====== 与步骤 07 复用 LOF 函数 ======
# find_bad_channels_lof 在 MNE 1.7+ 可用
from mne.preprocessing import find_bad_channels_lof

# 1) LOF 检测
bad_lof = find_bad_channels_lof(
    raw,
    n_neighbors=20,
    threshold=1.5,
    picks=['eeg'],
)

# 2) RANSAC 检测 — 注意：find_bad_channels_ransac 不在 MNE 核心 API
# 推荐替代：autoreject.Ransac 或 pyprep.NoisyChannels.find_bad_by_ransac

# ====== 基于滤波后 PSD 的精细检测 ======
psd = raw.compute_psd(picks='eeg', fmax=80)
log_power = 10 * np.log10(psd.get_data(power=True))

# 以全局中位数为基准
median_power = np.median(log_power.mean(axis=1))
mad = np.median(np.abs(log_power.mean(axis=1) - median_power))

# 滤波后更稳定的检测（sigma 调小）
threshold_high = median_power + 3 * mad
threshold_low = median_power - 3 * mad

log_pow_mean = log_power.mean(axis=1)
flat = [raw.ch_names[i] for i, p in enumerate(log_pow_mean) if p < threshold_low]
noisy = [raw.ch_names[i] for i, p in enumerate(log_pow_mean) if p > threshold_high]
print(f"PSD 平导: {flat}, 噪导: {noisy}")

# ====== 合并到 raw.info['bads'] ======
raw.info['bads'] = list(set(raw.info['bads']) | set(bad_lof) | set(flat) | set(noisy))
print(f"累计坏导: {raw.info['bads']}")

# 验证 1：坏导总数
n_bad = len(raw.info['bads'])
n_eeg = sum(1 for ch in raw.info['chs'] if ch['kind'] == 2)
print(f"坏导率: {n_bad}/{n_eeg} = {n_bad/n_eeg*100:.1f}%")
assert n_bad / n_eeg <= 0.20, "坏导过多，请人工复核"

# 验证 2：可视化（坏导标红）
raw.plot(duration=5, n_channels=30)
# 人工核对：坏导确实呈现明显异常

# 验证 3：分布合理性
print(raw.info['bads'])
# 期望：无聚集性（不在同一区域大量出现）