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

用Python和NumPy分析心电图手把手教你找到QRS波的核心频率附完整代码在生物医学信号处理领域心电图ECG分析一直是临床诊断和健康监测的重要工具。而QRS波作为ECG信号中最显著的特征波群其频率分布直接反映了心脏电活动的关键信息。本文将带您用Python和NumPy实现一套完整的ECG频谱分析流程从原始信号处理到QRS波核心频段提取最后通过pSQI指标评估信号质量。1. 环境准备与数据加载在开始分析前我们需要准备Python环境和示例ECG数据。推荐使用Anaconda创建虚拟环境conda create -n ecg_analysis python3.8 conda activate ecg_analysis pip install numpy matplotlib scipy对于ECG数据我们可以使用scipy内置的示例数据集或者加载自己的数据。以下是加载MIT-BIH心律失常数据库的示例代码import numpy as np from scipy.io import loadmat # 加载MAT格式的ECG数据 ecg_data loadmat(ecg_sample.mat) signal ecg_data[val][0] # 获取第一导联信号 fs 360 # 采样频率360Hz # 绘制原始信号 import matplotlib.pyplot as plt plt.figure(figsize(12,4)) plt.plot(np.arange(len(signal))/fs, signal) plt.xlabel(Time (s)) plt.ylabel(Amplitude (mV)) plt.title(Raw ECG Signal) plt.grid() plt.show()表常见ECG数据库采样频率参考数据库名称采样频率(Hz)信号时长典型用途MIT-BIH36030分钟心律失常分析PTB Diagnostic10002分钟心肌缺血诊断Fantasia250120分钟健康成人静息ECG2. 频谱分析基础与FFT实现理解离散傅里叶变换(DFT)是分析ECG频域特征的关键。NumPy提供的fft模块让我们可以轻松实现快速傅里叶变换(FFT)def compute_spectrum(signal, fs): n len(signal) yf np.fft.fft(signal) xf np.linspace(0, fs, n) magnitude np.abs(yf) return xf, magnitude # 计算频谱 frequencies, spectrum compute_spectrum(signal, fs) # 绘制频谱图 plt.figure(figsize(12,6)) plt.plot(frequencies[:len(frequencies)//2], spectrum[:len(spectrum)//2]) plt.xlabel(Frequency (Hz)) plt.ylabel(Magnitude) plt.title(ECG Frequency Spectrum) plt.grid() plt.show()需要注意的几个关键点频谱对称性实数信号的FFT结果在Nyquist频率(fs/2)处对称频率分辨率Δf fs/N其中N是信号长度频谱泄露可通过加窗函数(如Hamming窗)缓解提示实际分析时通常只使用前半部分频谱(0-fs/2)因为后半部分是前半部分的镜像。3. QRS波核心频段识别通过观察频谱图我们可以识别出QRS波的主要能量分布区域。以下是系统化的分析方法基线去除ECG信号通常包含低频基线漂移(0.5Hz)频带分割将频谱划分为多个关键区间能量计算统计各频段的能量占比def analyze_qrs_band(spectrum, frequencies): # 定义关键频段边界 bands { Baseline: (0, 0.5), P/T Waves: (0.5, 6), QRS Complex: (6, 30), High Frequency: (30, fs//2) } results {} total_energy np.sum(spectrum[(frequencies 0) (frequencies fs//2)]) for name, (low, high) in bands.items(): mask (frequencies low) (frequencies high) band_energy np.sum(spectrum[mask]) results[name] { energy: band_energy, percentage: band_energy / total_energy * 100 } return results band_analysis analyze_qrs_band(spectrum, frequencies) for band, data in band_analysis.items(): print(f{band}: {data[percentage]:.1f}% energy)表典型ECG信号各频段能量分布频段名称频率范围(Hz)能量占比(%)对应生理活动基线漂移0-0.55-15呼吸运动等P/T波0.5-620-30心房活动QRS波6-3050-70心室去极化高频噪声305-15肌电干扰等4. pSQI指标计算与信号质量评估功率谱质量指数(pSQI)是评估ECG信号质量的重要指标它计算QRS波频段(6-30Hz)能量与全频段(0-125Hz)能量的比值def calculate_psqi(signal, fs): n len(signal) yf np.fft.fft(signal) xf np.linspace(0, fs//2, n//2) # 计算QRS频段能量 qrs_mask (xf 6) (xf 30) qrs_energy np.sum(np.abs(yf[:n//2][qrs_mask])**2) # 计算全频段能量 full_mask (xf 0) (xf 125) full_energy np.sum(np.abs(yf[:n//2][full_mask])**2) return qrs_energy / full_energy psqi_score calculate_psqi(signal, fs) print(fpSQI score: {psqi_score:.3f})pSQI值的解读指南0.5-0.7优质信号适合临床分析0.3-0.5中等质量可能需要预处理0.3信号质量差建议重新采集5. 完整代码实现与优化技巧将上述步骤整合为完整的分析流程并添加实用优化import numpy as np from scipy.signal import butter, filtfilt import matplotlib.pyplot as plt def preprocess_ecg(signal, fs): # 高通滤波去除基线漂移 b, a butter(4, 0.5/(fs/2), highpass) filtered filtfilt(b, a, signal) return filtered def full_analysis(signal, fs): # 预处理 signal preprocess_ecg(signal, fs) # 频谱分析 n len(signal) yf np.fft.fft(signal) xf np.linspace(0, fs, n) spectrum np.abs(yf) # 频段分析 bands analyze_qrs_band(spectrum, xf) # 质量评估 psqi calculate_psqi(signal, fs) # 可视化 plt.figure(figsize(15,10)) plt.subplot(211) plt.plot(np.arange(n)/fs, signal) plt.title(Processed ECG Signal) plt.xlabel(Time (s)) plt.grid() plt.subplot(212) plt.plot(xf[:n//2], spectrum[:n//2]) plt.title(Frequency Spectrum) plt.xlabel(Frequency (Hz)) plt.grid() plt.tight_layout() plt.show() return { spectrum: (xf[:n//2], spectrum[:n//2]), band_analysis: bands, psqi: psqi } # 执行完整分析 results full_analysis(signal, fs)实际项目中遇到的几个优化点滤波参数调整根据具体ECG设备调整高通滤波截止频率频谱平滑使用移动平均或Savitzky-Golay滤波器平滑频谱曲线多导联分析当有多个导联数据时可综合各导联结果提高准确性