Python实现60Hz陷波滤波器无效问题排查及优化咨询

核心问题定位 #

你的代码出现滤波无效/效果差的核心原因集中在采样频率准确性和滤波器参数适配上：

• 采样频率（fs）错误：初始假设的1000Hz和后来计算的~192.3Hz（1/0.0052）都可能和ECG数据的真实采样率不匹配。signal.iirnotch完全依赖准确的fs来计算陷波的中心频率，fs错了，滤波器实际抑制的频率根本不是真实的60Hz干扰。
• Q值设置过窄：Q=30会让陷波带宽极窄，实际环境中的60Hz干扰往往带有谐波或频率漂移，窄带宽滤波器无法有效覆盖这些成分。
• 缺少预处理：ECG信号的基线漂移可能掩盖滤波效果，先去除基线再做陷波会更明显。

修复步骤与完整代码 #

1. 先确认真实采样率

如果CSV文件包含时间戳列（比如第一列是时间），直接从数据计算准确的采样间隔：

# 假设第一列是时间戳，提取并计算采样间隔
timestamps = [float(row[0]) for row in data]
T = np.mean(np.diff(timestamps))  # 计算平均采样间隔
fs = 1 / T
print(f"真实采样频率：{fs:.2f} Hz")

如果没有时间戳，必须从采集设备说明书或实验记录中获取准确fs，不能凭假设。

2. 优化滤波器参数与预处理

调整Q值到5-15，同时加入高通滤波器去除基线漂移，再应用陷波：

import csv
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# 加载ECG数据
with open('experiment_data60(copy).csv', 'r') as f:
    reader = csv.reader(f)
    data = list(reader)

# 截取数据（保持你的原有范围）
data = data[900:1501] 
# 提取信号值和时间戳（如果有）
values = np.array([float(row[1]) for row in data])
timestamps = np.array([float(row[0]) for row in data])

# 计算真实采样频率
T = np.mean(np.diff(timestamps))
fs = 1 / T
print(f"采样频率：{fs:.2f} Hz")
N = len(values)
t = timestamps  # 直接用原始时间戳更准确

# --- 第一步：去除基线漂移（高通滤波）---
# 设计0.5Hz高通滤波器，去除低频基线波动
high_pass_b, high_pass_a = signal.butter(2, 0.5, 'high', fs=fs)
values_no_baseline = signal.filtfilt(high_pass_b, high_pass_a, values)

# --- 第二步：设计60Hz陷波滤波器 ---
f0 = 60.0
Q = 10.0  # 调整Q值，拓宽陷波带宽
b, a = signal.iirnotch(f0, Q, fs)
filtered_values = signal.filtfilt(b, a, values_no_baseline)

# --- 绘图对比（时域+频谱）---
plt.figure(figsize=(12, 8))

# 原始信号 vs 滤波后信号
plt.subplot(2,1,1)
plt.plot(t, values_no_baseline, 'b', alpha=0.5, label='去基线后信号')
plt.plot(t, filtered_values, 'r', label='陷波滤波后信号')
plt.xlabel('时间 (s)')
plt.ylabel('幅值')
plt.title('ECG信号时域对比')
plt.legend()
plt.grid(True)

# 频谱对比
plt.subplot(2,1,2)
freq, psd_original = signal.welch(values_no_baseline, fs=fs, nperseg=256)
freq, psd_filtered = signal.welch(filtered_values, fs=fs, nperseg=256)
plt.plot(freq, psd_original, 'b', alpha=0.5, label='去基线后频谱')
plt.plot(freq, psd_filtered, 'r', label='陷波滤波后频谱')
plt.xlabel('频率 (Hz)')
plt.ylabel('功率谱密度')
plt.title('ECG信号频谱对比（重点看60Hz处）')
plt.xlim(0, 100)  # 聚焦低频段
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

验证与调试技巧 #

• 频谱分析优先：时域图可能看不出细微变化，用signal.welch绘制功率谱，直接观察60Hz峰值是否被抑制，这是最直观的验证方法。
• 调整Q值：如果60Hz干扰有谐波（比如120Hz），可以再添加一个120Hz的陷波滤波器，或者降低Q值拓宽带宽。
• 检查数据完整性：确认CSV数据没有缺失、异常值，比如有没有非数值的行导致信号中断。

内容的提问来源于stack exchange，提问作者Lily