用Python和NumPy分析心电图：手把手教你找出QRS波的核心频率（附完整代码）

用Python和NumPy分析心电图手把手教你找出QRS波的核心频率附完整代码在生物医学信号处理领域心电图ECG分析一直是研究热点。QRS波作为ECG信号中最显著的特征之一其频率分布直接反映了心脏电活动的关键信息。本文将带你用Python实现一套完整的ECG频谱分析流程从原始信号处理到QRS波核心频段定位最后计算信号质量指标pSQI。不同于教科书式的理论讲解我们更关注实际编码中的坑点规避和结果解读。1. 环境准备与数据加载开始前确保已安装必要的Python库pip install numpy matplotlib scipy典型的ECG数据通常以CSV或MAT格式存储。我们使用来自MIT-BIH心律失常数据库的样本数据可通过PhysioNet获取import numpy as np import matplotlib.pyplot as plt # 模拟生成ECG信号实际应用应替换为真实数据 fs 250 # 采样频率250Hz t np.arange(0, 10, 1/fs) # 10秒时长 ecg_signal 0.5 * np.sin(2 * np.pi * 1 * t) # 基础节律 ecg_signal 1.2 * np.sin(2 * np.pi * 12 * t) # QRS主频 ecg_signal 0.3 * np.random.randn(len(t)) # 添加噪声注意实际项目中建议使用wfdb库直接读取标准ECG数据库import wfdb record wfdb.rdrecord(mitdb/100, sampto3000) ecg_signal record.p_signal[:,0]2. 频谱分析核心步骤2.1 快速傅里叶变换实现NumPy的FFT模块让频域分析变得简单但有几个关键参数需要注意n len(ecg_signal) fft_result np.fft.fft(ecg_signal) # 执行FFT freqs np.fft.fftfreq(n, d1/fs) # 生成频率轴 # 取单边频谱利用对称性 half_n n // 2 magnitude np.abs(fft_result[:half_n]) * 2 / n freqs freqs[:half_n]常见问题排查频谱出现异常峰值 → 检查信号是否包含NaN值幅度量级过小 → 确认是否漏乘2/n的缩放系数频率轴不正确 → 验证采样频率参数2.2 QRS波频段定位通过频谱能量分布确定QRS波核心频段def find_peak_band(freqs, magnitude, min_width5): # 寻找超过平均能量3倍的连续频段 threshold np.mean(magnitude) * 3 mask magnitude threshold changes np.where(np.diff(mask))[0] if len(changes) 2: return None start_idx changes[0] 1 end_idx changes[1] while (end_idx - start_idx) min_width and len(changes) 2: start_idx changes[0] 1 end_idx changes[2] changes changes[2:] return freqs[start_idx], freqs[end_idx] qrs_low, qrs_high find_peak_band(freqs, magnitude) print(fQRS核心频段: {qrs_low:.1f}-{qrs_high:.1f} Hz)3. 信号质量评估pSQI功率谱质量指数(pSQI)是评估ECG信号质量的重要指标def calculate_psqi(freqs, magnitude, qrs_band(6,30), full_band(0,125)): # 计算频段能量 qrs_mask (freqs qrs_band[0]) (freqs qrs_band[1]) full_mask (freqs full_band[0]) (freqs full_band[1]) qrs_energy np.sum(magnitude[qrs_mask]**2) total_energy np.sum(magnitude[full_mask]**2) return qrs_energy / total_energy psqi_score calculate_psqi(freqs, magnitude) print(fpSQI指数: {psqi_score:.3f} (0.5-0.7为优质信号))4. 完整可视化分析流程将上述步骤整合为可交互的分析工具def analyze_ecg(ecg_signal, fs): plt.figure(figsize(12, 8)) # 时域信号 plt.subplot(3, 1, 1) plt.plot(np.arange(len(ecg_signal))/fs, ecg_signal) plt.title(原始ECG信号) plt.xlabel(时间(s)) # 频域分析 n len(ecg_signal) fft_result np.fft.fft(ecg_signal) freqs np.fft.fftfreq(n, d1/fs)[:n//2] magnitude np.abs(fft_result[:n//2]) * 2 / n plt.subplot(3, 1, 2) plt.plot(freqs, magnitude) plt.xlim(0, 50) # 聚焦0-50Hz plt.title(频谱分析) plt.xlabel(频率(Hz)) # QRS频段标记 qrs_low, qrs_high 6, 30 # 或使用find_peak_band结果 plt.axvspan(qrs_low, qrs_high, colorred, alpha0.2) # pSQI计算 psqi_score calculate_psqi(freqs, magnitude) plt.subplot(3, 1, 3) plt.bar([QRS能量,其他能量], [psqi_score, 1-psqi_score], color[red,grey]) plt.title(f能量分布 (pSQI{psqi_score:.3f})) plt.tight_layout() plt.show() analyze_ecg(ecg_signal, fs)5. 实战技巧与性能优化处理长时程ECG信号时可采用分段分析方法def segment_analysis(signal, fs, window_sec5): window_size window_sec * fs results [] for i in range(0, len(signal), window_size): segment signal[i:iwindow_size] if len(segment) window_size: continue n len(segment) fft_result np.fft.fft(segment) freqs np.fft.fftfreq(n, d1/fs)[:n//2] magnitude np.abs(fft_result[:n//2]) * 2 / n qrs_band find_peak_band(freqs, magnitude) psqi calculate_psqi(freqs, magnitude, qrs_band) results.append({ start_time: i/fs, qrs_band: qrs_band, psqi: psqi }) return pd.DataFrame(results)内存优化技巧对超长信号使用np.lib.stride_tricks.sliding_window_view实时处理场景考虑scipy.signal.stft短时傅里叶变换使用numba加速关键计算循环