简介这是一份围绕心电信号处理的中文技术文档面向数字信号处理学习者、生物医学工程初学者及需要完成相关实验报告的高校学生系统讲解从基础理论到频谱分析、能量谱分析、滤波去噪的完整链路内容覆盖傅里叶变换、小波变换、时域与频域特征提取并结合MIT-BIH数据库案例展示工程实现细节。资源为1个pdf文件压缩包整体约874KB方便直接阅读和检索适合作为课程作业、实验设计或论文写作的参考资料。文档还包含幅度谱、相位谱、功率谱、自相关与互相关分析等关键知识点并对滤除基线漂移和噪声干扰给出具体处理思路能帮助读者理解随机信号处理的核心方法。该资源已有4078人学习浏览是一份实用性强、性价比高的入门与进阶参考文档。1. 一套能跑的 ECG 信号处理基线从 MIT-BIH 的 2048 点切片说起做心电信号处理的人手上大概率都下载过 MIT-BIH 心律失常数据库的公开数据但真到自己动手做特征提取时很多人第一步就卡在“数据怎么读、谱怎么画、滤波器参数怎么定”上。这篇博文要拆的项目就是基于 MIT-BIH 数据库四组病例数据100-2-3、105-2-3、109-2-3、111-2-3的完整 ECG 信号分析流程每组数据长度 N2048采样率 fs360Hz即约 5.69 秒的心电片段。这个采样率不是随便定的360Hz 是 MIT-BIH 原始记录设备的硬件采集频率意味着设计 FIR、IIR 滤波器时归一化频率必须按 360 来算否则滤波器通带会整体偏移。项目覆盖了频谱分析、相关性分析、线性滤波器设计、维纳滤波四个模块对需要快速搭建 ECG 信号处理基线、做课程设计或算法预研的工程师来说这套流程能直接当脚手架用——替换数据路径改几个参数就能复用到你自己的心电采集设备数据上。信号极其微弱典型幅值约 1mV、频率范围低0.135Hz能量集中在 520Hz、易受工频 50Hz 干扰和基线漂移影响这三条是 ECG 处理里的老生常谈但真正在代码里同时处理它们的人并不多。下面按项目原流程逐步展开。2. 频谱分析幅度谱、相位谱和功率谱的计算与解读2.1 为什么 ECG 分析先看频谱而不直接看波形原始 ECG 波形肉眼可见“毛刺”和基线起伏但仅靠时域波形判断噪声来源非常困难。频谱分析的价值在于它把信号分解成不同频率分量的叠加让你能直接看到能量集中在哪些频段、哪些频段是被噪声污染的。这个项目里明确要求分析幅度谱和相位谱。幅度谱是傅里叶变换结果取模后关于频率的函数反映各频率分量的强度相位谱是辐角关于频率的函数反映各分量的相位信息。能量谱和功率谱则是从能量角度观察信号。这里要区分两件事信号在 t 时刻的瞬时功率是 |f(t)|²而信号的能量是瞬时功率在整个时间轴上的积分。对能量无限、平均功率有限的信号比如平稳随机过程只能通过功率谱密度来看它的频域分布。心电信号本质上是随机过程所以功率谱分析比能量谱更常用。2.2 用 numpy 计算幅度谱、相位谱和功率谱项目附录的代码 1 给了完整实现。核心计算逻辑如下import numpy as np from scipy.io import loadmat # 载入数据数据结构为 mat 文件中的 y 字段 data_100 loadmat(./100-2-3.mat) Y_100 data_100[y] # shape 为 (2048, 1) Fs 360 # 采样率来自硬件采集配置 N 2048 # 单组数据长度对应约 5.69s frq np.arange(N) # 频率轴原始索引 # 傅里叶变换注意这里用的是 fft2对一维信号与 fft 等价 mf_100 np.fft.fft2(Y_100) mag_100 np.abs(mf_100) # 双边幅度谱未归一化 gui_mag_100 mag_100 / N # 归一化幅度谱双边 gui_half_mag_100 gui_mag_100[:N//2] # 取单边利用对称性 # 相位谱单位换算为角度 angle_mag_100 180 * np.angle(mf_100) / np.pi # 功率谱P |X|^2 / N ps_100 (mag_100[:N//2] ** 2) / N ps_db_100 20 * np.log10(ps_100) # 转 dB 方便观察动态范围参数说明和逻辑解释np.fft.fft2对一维数组和np.fft.fft结果相同项目中用它可能是为了统一代码写法实际一维信号用fft更快且更直观。归一化的目的是让变换结果不随 N 变化而失去可比性。mag_100 / N对应的是周期图法中的幅值归一化功率谱要除以 N 而不是 N²这是由帕塞瓦尔定理确定的sum(|x|²) (1/N) * sum(|X[k]|²)。20dB 换算只用于观察不改变数据本身的物理意义。比如幅度谱中-40 到 60 这个范围在 dB 视角下对应从极微弱分量到主导频率分量的巨大跨度。相位谱计算用np.angle得到的弧度值乘以 180/π 转角度便于和项目中图例一致。相位谱在值域上是 (-180,180] 的区间直接绘图时会出现跳变这不是 bug而是反正切函数的周期截断。2.3 频域结果如何对应到 ECG 病理特征项目中明确观察到 100-2-3 与 105-2-3 两组信号的幅度谱有明显差异幅度范围不同相位谱也表现出不同的跳变模式。一个值得注意的细节是除了看幅度谱的峰值位置和高度相位谱的模式往往携带更多关于波形形态的信息。比如正常的 QRS 波群具有陡峭的上升沿和下降沿对应相位谱在特定频段快速变化而基线漂移导致的缓变信号则在低频段表现为相位平缓。实际分析我的做法是先打印几个关键频点的具体数值# 找出前 5 个幅度最大的频率分量 top_k_idx np.argsort(gui_half_mag_100.flatten())[-5:] top_k_freqs top_k_idx * Fs / N # 换算为实际物理频率 print(Top-5 频率分量(Hz):, top_k_freqs)这里把索引换算成实际频率用的是索引 * Fs / N的公式。2048 点做完 FFT 后频率分辨率为Fs/N 360/2048 ≈ 0.176Hz意味着你能区分的最小频率差约为 0.176Hz。如果两组信号的频谱差异在 0.176Hz 以内当前数据长度无法分辨需要更长的观察窗或零填充。这是一个经常被忽略的分辨率陷阱。3. 时域统计与相关性分析109 与 111 组数据的相似度度量3.1 均值、方差在 ECG 分析里具体看什么均值反映了信号的直流偏置水平。ECG 信号本身有约 1mV 的典型幅值但硬件采集时可能叠加一个直流偏置这个偏置在均值上直接体现。方差则描述信号围绕均值的波动程度对应 ECG 信号的动态范围。如果两组数据的均值接近但方差差异很大说明它们的基线水平一致但信号振幅或噪声能量不同。项目的代码 2 直接调用了np.mean和np.var。需要注意np.var默认计算总体方差除以 N而我们在统计推断中常用的无偏样本方差除以 N-1需要设置ddof1。对 N2048 的数据两者差距小于 0.1%但如果你后续要和其他研究的统计量对比最好统一用无偏估计。3.2 自相关与互相关用 numpy 实现并从结果读信息自相关函数描述信号自身在不同时刻取值之间的相关程度。对心电信号它的最大价值在于周期性的心跳会在自相关函数中表现为等间距的峰值且这些峰值不会像原始波形那样被噪声淹没。这是因为自相关运算相当于对信号做了能量累积随机噪声在自相关中趋于零而周期分量被保留。互相关函数则度量两个信号在不同相对位移下的相似程度。项目中用np.correlate(Y_109, Y_111, modefull)得到完整互相关序列峰值位置对应两个信号最大相似时的延迟量。原理上互相关和卷积只差一个翻转卷积是f(t)与g(-t)的滑动内积互相关是f(t)与g(t)的直接滑动内积。所以互相关不需要翻转任何一个序列。Y_109 np.array(data_109[y]).flatten() Y_111 np.array(data_111[y]).flatten() # 归一化互相关标准化相关系数消除幅值差异的影响 def normalized_cross_correlation(x, y): x (x - np.mean(x)) / np.std(x) y (y - np.mean(y)) / np.std(y) # 有效长度范围内的互相关 return np.correlate(x, y, modefull) / len(x) corr normalized_cross_correlation(Y_109, Y_111) # 定位峰值和对应延迟 delay np.argmax(corr) - (len(Y_109) - 1) print(f最大相关系数: {np.max(corr):.4f}, 延迟点数: {delay}, 延迟时间: {delay / 360:.3f}s)这里做了两个关键改进减均值、除标准差把相关系数归一化到 [-1, 1] 区间不同幅值水平的信号才能横向比较。np.argmax(corr) - (len(Y_109) - 1)是因为modefull的输出长度是2N-1索引 N-1 对应零延迟偏移量才能换算为正负延迟。如果互相关峰值接近 1说明 109 和 111 两组信号在波形形态上高度相似如果峰值位置偏离中心点说明两个记录之间存在时间错位可能是硬件采集时触发点不同造成的。项目正文只说了“信号 109,111 的互相关函数”这个观察项但具体量化用上述方法实现最直接。3.3 项目中遗漏的周期估计技巧自相关函数还有一个高频用途估计瞬时心率。已知 fs360N2048如果自相关函数在零延迟附近出现第一个明显峰值排除 τ0 的极大值按该峰值对应的延迟索引可计算心搏周期from scipy.signal import find_peaks ry np.correlate(Y_109, Y_109, modefull) ry_mid ry[N-1:] # 取零延迟及之后的半段 # 找第一个非零延迟的局部极大值高度设为最大值的 30% peaks, props find_peaks(ry_mid, height0.3 * np.max(ry_mid), distance30) if len(peaks) 1: rr_interval peaks[1] # 第二个峰第一个是零延迟 heart_rate 60 * Fs / rr_interval print(f估计心率: {heart_rate:.1f} bpm)distance30对应最小 30 个采样点约 83ms的峰间距这相当于 720bpm 的上限能防住误检。这个方法在信号信噪比低时比直接检测 R 波更稳定因为自相关做了噪声平均。4. 数字滤波器设计与去噪FIR、IIR 和零相移滤波选型实践4.1 线性相位在 ECG 里的实际意义项目正文比较了 FIR 和 IIR 的优缺点核心分歧在相位特性。IIR 滤波器如 Butterworth、Chebyshev能用低阶达到陡峭的过渡带计算量小但非线性的相位响应会改变 ECG 各频率分量的相对时间关系。这对 ECG 来说不是小事QRS 波群是陡峭的瞬态波形相位失真会把波的起点和终点模糊化影响 ST 段的测量。FIR 滤波器可以设计成严格线性相位也就是群延迟恒定各频率分量经过滤波后延迟相同的时间波形形态保持不变。另一个关键差异是稳定性。IIR 因为有反馈递归结构极点必须在单位圆内才能稳定FIR 没有反馈是无条件稳定的。在 ECG 这个场景里信号本身是非平稳的IIR 在极端输入下可能出现溢出而 FIR 不存在这个问题。具体设计时我倾向于用scipy.signal.butter的filtfilt做零相移滤波但要清楚它和普通lfilter的区别。下面给出一套处理 100-2-3 信号的完整示例import numpy as np from scipy import signal def design_ecg_filter(fs, f_low0.5, f_high35.0, order4, ftypebutter): 设计带通滤波器滤除基线漂移(0.5Hz)和高频噪声(35Hz) 返回 SOS 格式系数数值稳定性更好 nyq fs / 2.0 # 归一化频率截止频率 / 奈奎斯特频率 low f_low / nyq high f_high / nyq sos signal.butter(order, [low, high], btypebandpass, outputsos) b, a signal.butter(order, [low, high], btypebandpass, outputba) return sos, (b, a) sos, _ design_ecg_filter(fs360) # 零相移滤波前向反向各跑一遍相位失真抵消 Y_100_filtered signal.sosfiltfilt(sos, Y_100.flatten()) # 对比普通滤波会引入群延迟 Y_100_filtered_lfilter signal.sosfilt(sos, Y_100.flatten())参数逻辑说明截止频率 0.5Hz 和 35Hz 根据 ECG 的能量分布确定0.5Hz 以下主要是呼吸引起的基线漂移35Hz 以上多为肌电噪声和机器噪声。50Hz 工频干扰在 35Hz 之上所以这个带通滤波器本身就把它削弱了不需要单独的陷波滤波器。用sos二阶节级联而不是直接的[b, a]是因为高阶滤波器直接展开成多项式时数值精度急剧下降sos把系统拆成多个二阶节的级联每个节都有独立的极点避免因极点靠太近引起的数值误差。sosfiltfilt是零相移滤波的实现等价于先正向滤波再反向滤波零相位失真的代价是引入了非因果性用到了“未来”的样本。这在离线处理比如分析已采集的 MIT-BIH 数据时完全可行但不能用于实时监护设备。4.2 工频干扰 50Hz 的定向抑制策略项目最后在维纳滤波部分让 111-2-3 信号混入 50Hz 噪声这是刻意模拟工频干扰的场景。直接用上面 0.535Hz 的带通滤波器可以压掉它但如果你的通带需要扩展到 50Hz 以上就得用陷波滤波器。def notch_filter(fs, f050.0, Q30.0): 设计 50Hz 陷波滤波器Q 值控制陷波带宽 b, a signal.iirnotch(f0, Q, fs) return b, a b_notch, a_notch notch_filter(fs360) Y_111_notched signal.lfilter(b_notch, a_notch, Y_111.flatten()) # 也可用 filtfilt 零相移版本 Y_111_notched_2 signal.filtfilt(b_notch, a_notch, Y_111.flatten())Q 值在这里的含义Qf0/BWQ30 对应约 1.67Hz 的陷波带宽50±0.83Hz。Q 值太高容易把 50Hz 邻近的有效心电分量也削掉Q 值太低则抑制效果不明显。项目中在维纳滤波对比时混合进了 50Hz 噪声我的验证做法是滤波后对输出做 FFT看 50Hz 处的幅度下降了多少而不是只凭肉眼观察时域波形。4.3 FIR 设计用 window 法做等波纹之外的快速方案如果选型方向是 FIRscipy.signal.firwin是最方便的起点。它用 window 法默认 Hamming 窗设计线性相位 FIR 滤波器指定截止频率、过渡带宽和滤波器阶数即可。# 40Hz 低通 FIR阶数 101足够陡的过渡带窗函数选 Hamming N_fir 101 fir_coeff signal.firwin(N_fir, 35.0 / (360 / 2), windowhamming) # 频域响应验证 w, h signal.freqz(fir_coeff, worN8192) freqs_hz w * 360 / (2 * np.pi) # 找出 -3dB 点的实际频率 idx_3db np.where(np.abs(h) 0.707)[0] if len(idx_3db) 0: f_3db freqs_hz[idx_3db[0]] print(f实际 -3dB 频点约 {f_3db:.2f} Hz)阶数 101 对应群延迟 50 个采样点约 139ms。如果后续要对滤波后的信号做 R 波检测或者和另一路信号做同步对比这个固定延迟要补偿掉否则时间轴对不上。FIR 滤波后的波形没有 IIR 那种相位畸变但代价是同样的过渡带需要更高阶数——在这个 2048 点、fs360Hz 的数据集上101 阶 FIR 的运算时间远小于毫秒级完全不构成计算瓶颈。4.4 项目里“线性滤波去掉基线漂移”的完整复现基线漂移是 ECG 处理里最常见的干扰。项目提到“线性滤波去掉基线漂移频谱”我的标准做法是两步先用高通或带通把漂移滤掉再看频谱确认低频段被压制。# 高通滤波器截止频率 0.5Hz压掉基线漂移 b_hp, a_hp signal.butter(3, 0.5 / (360 / 2), btypehigh) Y_109_hp signal.filtfilt(b_hp, a_hp, Y_109.flatten()) # 滤波前后的低频段功率对比 freq np.fft.rfftfreq(N, d1/360) before np.fft.rfft(Y_109.flatten()) after np.fft.rfft(Y_109_hp) low_band_before np.sum(np.abs(before[(freq 0.1) (freq 0.5)]) ** 2) low_band_after np.sum(np.abs(after[(freq 0.1) (freq 0.5)]) ** 2) print(f0.1~0.5Hz 频段功率衰减: {10 * np.log10(low_band_after / low_band_before):.1f} dB)rfft对实信号只计算正频率部分输出长度 N/21比fft节省近一半计算量。功率衰减用 dB 表示负值越大说明滤得越干净。我一般要求至少衰减 20dB否则基线漂移依然会对后续 ST 段分析产生干扰。5. 维纳滤波基于最小均方误差的 ECG 去噪5.1 维纳滤波的适用边界为什么项目把 111 组数据混入 50Hz 噪声再滤波维纳滤波是线性最小均方误差估计器也就是说它在“信号与噪声均为平稳随机过程且二者不相关”的假设下是最优的。这个场景与 ECG 高度匹配心电信号本身是准平稳的且与工频干扰、高斯白噪声在统计上不相关。其核心是维纳-霍夫方程R_sx(m) Σ h(i) R_xx(m-i)其中 R_sx 是期望信号与观测信号的互相关R_xx 是观测信号的自相关。求解这个方程得到最优冲激响应 h(i)滤波输出 y(n) Σ h(i) x(n-i)。项目附录中把 111-2-3 混入 50Hz 噪声后用维纳滤波这比单纯用陷波器更接近“利用信号统计特性”的思路陷波器是固定频带切除不区分信号和噪声维纳滤波器会根据观测信号的自相关结构自动推断可能在哪些频点有信号从而在去除噪声的同时尽量保留信号分量。from scipy.signal import wiener # 给 111-2-3 加入 50Hz 正弦噪声和少量白噪声 t np.arange(N) / Fs noise_50hz 0.5 * np.sin(2 * np.pi * 50 * t) Y_111_clean Y_111.flatten() Y_111_noisy Y_111_clean noise_50hz 0.05 * np.random.randn(N) # scipy.signal.wiener 基于局部方差做维纳滤波 Y_111_wiener wiener(Y_111_noisy, mysize11, noiseNone)这里需要解释mysize和noise的含义。mysize11表示用 11 个点的滑动窗口估计局部均值和方差窗口越大滤波器越平滑但信号细节丢失越多。noiseNone时算法用整段信号的中位绝对偏差自动估计噪声方差。如果想手动指定噪声级别可以传入一个浮点数。缺点是 scipy 的wiener是局部实现和理论上的全局维纳滤波有差距在离线场景效果尚可但不推荐在实时系统中直接使用。5.2 自己实现全局维纳滤波更贴近教材公式如果想让结果更精确可控我一般自己算用 Toeplitz 矩阵解维纳-霍夫方程。在 2048 点数据上估计 32 阶滤波器系数就够用from scipy.linalg import toeplitz, solve_toeplitz def wiener_optimal_filter(noisy, clean, taps32): 根据已知参考信号计算最优维纳滤波器系数 noisy noisy - np.mean(noisy) clean clean - np.mean(clean) # 估计自相关和互相关 rxx np.correlate(noisy, noisy, modefull) / len(noisy) rxx rxx[len(noisy) - 1: len(noisy) taps] # 取前 taps 个延迟的自相关 rxs np.correlate(noisy, clean, modefull) / len(noisy) rxs rxs[len(noisy) - 1: len(noisy) taps] # 构造 Toeplitz 矩阵并求解 R toeplitz(rxx, rxx[:taps]) h_opt np.linalg.solve(R, rxs) return h_opt h_opt wiener_optimal_filter(Y_111_noisy, Y_111_clean, taps32) Y_111_wiener_custom np.convolve(Y_111_noisy, h_opt, modesame)核心逻辑rxx是对称的 Toeplitz 矩阵的第一行它编码了观测信号的自相关结构。rxs是观测与期望的互相关。最小二乘求解R h rxs得到最优系数。和 scipy 的局部维纳相比这种方法利用的是全局统计量滤波器是因果的可复用到流式数据的块处理。用np.convolve(..., modesame)输出长度与输入一致但首尾 taps/2 个点的滤波结果不完全准确因为卷积在这里用了零填充。5.3 维纳滤波效果验证信噪比改善度怎么算滤波做完了不能只看图说“看起来干净了”要在频域和时域同时量化。def snr_db(signal, noise): signal: 有用信号, noise: 噪声 ps np.sum(signal ** 2) pn np.sum(noise ** 2) return 10 * np.log10(ps / pn) noise_before Y_111_noisy - Y_111_clean noise_after Y_111_wiener_custom - Y_111_clean snr_before snr_db(Y_111_clean, noise_before) snr_after snr_db(Y_111_clean, noise_after) print(f滤波前 SNR: {snr_before:.2f} dB) print(f滤波后 SNR: {snr_after:.2f} dB) print(fSNR 改善: {snr_after - snr_before:.2f} dB)如果 SNR 改善是负值说明滤波器引入了额外失真通常是 taps 数选取不当或信号不满足平稳假设。这时我会把 taps 从 32 降 16 重新计算或者先做一次 0.535Hz 带通预处理把明显不在信号频带内的噪声先减掉再跑维纳滤波效果一般会好一截。项目中的“功率谱密度和频谱图像”对比也是同样的目的——滤波后的功率谱应在 50Hz 处明显凹陷但在 520Hz 的心电主能量区间保持平坦。5.4 验证滤波结果的干净程度检查残差自相关一个快速判定滤波是否过度的技巧计算残差滤波输出减原始干净信号的自相关函数。如果残差在零延迟处有明显尖峰且其他位置迅速衰减到零说明残差主要是白噪声滤波有效且没“滤过头”。如果残差自相关出现周期性振荡说明信号分量也被削掉了。residual Y_111_wiener_custom - Y_111_clean residual_ac np.correlate(residual, residual, modefull) residual_ac residual_ac[len(residual) - 1:] # 只看后半段 residual_ac / np.max(residual_ac) # 归一化 # 找零延迟之后 100ms 内是否有超过 0.1 的旁瓣 side_lobe_region residual_ac[10:int(0.1 * Fs)] if np.max(side_lobe_region) 0.1: print(警告: 残差中可能存在周期性信号泄漏) else: print(残差近似白噪声滤波效果良好)这个验证方法适用于任何线性去噪手段IIR、FIR、维纳滤波比单纯看输出波形更客观。它利用的是白噪声自相关函数在非零延迟处为零的数学性质——如果残差里还有结构化的信号成分自相关会表现出来。本文还有配套的精品资源点击获取