MIT-BIH-AF房颤数据集实战：用Python库wfdb和PyWavelets搞定ECG信号处理（附完整代码）

Lang Run

542人浏览 · 2026-05-22 12:52:41

Lang Run · 2026-05-22 12:52:41 发布

MIT-BIH-AF房颤数据集实战：Python信号处理全流程解析

1. 环境准备与数据加载

处理MIT-BIH-AF数据集的第一步是搭建合适的Python环境。推荐使用Anaconda创建独立环境，避免库版本冲突：

conda create -n ecg python=3.8
conda activate ecg
pip install wfdb PyWavelets scipy matplotlib numpy

数据集文件通常包含三种关键格式：

.dat ：存储原始ECG信号
.atr ：包含医生标注的心律事件
.qrs ：记录R峰位置信息

使用wfdb库加载数据的正确方式：

import wfdb

# 示例：加载患者04015的数据
record = wfdb.rdrecord('04015', pb_dir='afdb') 
annotations = wfdb.rdann('04015', 'atr', pb_dir='afdb')

常见问题排查表：

错误类型	可能原因	解决方案
FileNotFoundError	文件路径错误	确认PhysioNet文件结构
ValueError	文件格式不匹配	检查文件扩展名完整性
ImportError	依赖库缺失	重新安装wfdb库

提示：首次使用PhysioNet数据集需注册账号并签署数据使用协议

2. 信号可视化与基础分析

原始ECG信号通常包含多种噪声成分，需要先进行初步观察：

import matplotlib.pyplot as plt

plt.figure(figsize=(12,4))
plt.plot(record.p_signal[:2000,0])
plt.title('Raw ECG Signal (Lead I)')
plt.xlabel('Samples')
plt.ylabel('Amplitude (mV)')
plt.grid()
plt.show()

典型信号特征分析流程：

识别采样率 ： record.fs (通常为250Hz)
检查导联数 ： record.n_sig (MIT-BIH-AF为双导联)
提取时间轴 ： t = np.arange(len(ecg))/record.fs

关键统计指标计算：

import numpy as np

ecg = record.p_signal[:,0]
print(f"信号长度: {len(ecg)} samples")
print(f"持续时间: {len(ecg)/record.fs/60:.1f} minutes")
print(f"幅值范围: {np.min(ecg):.2f} ~ {np.max(ecg):.2f} mV")

3. 高级信号处理技术

3.1 小波去噪实战

PyWavelets提供多级小波分解能力，特别适合ECG信号处理：

import pywt

def wavelet_denoise(signal, wavelet='db4', level=3):
    coeffs = pywt.wavedec(signal, wavelet, level=level)
    sigma = np.median(np.abs(coeffs[-1]))/0.6745
    uthresh = sigma * np.sqrt(2*np.log(len(signal)))
    coeffs[1:] = [pywt.threshold(c, uthresh, mode='soft') for c in coeffs[1:]]
    return pywt.waverec(coeffs, wavelet)

去噪效果对比参数：

小波基	计算效率	保真度	适用场景
db4	高	优	常规ECG
sym5	中	极佳	高精度分析
haar	极高	良	实时处理

3.2 R峰检测优化

改进的Pan-Tompkins算法实现：

def detect_r_peaks(ecg, fs=250):
    # 带通滤波
    b = signal.firwin(31, [5, 15], fs=fs, pass_zero=False)
    filtered = signal.lfilter(b, 1, ecg)
    
    # 微分处理
    diff = np.diff(filtered, prepend=0)
    
    # 平方增强
    squared = diff ** 2
    
    # 移动平均
    window = int(0.15 * fs)
    integrated = np.convolve(squared, np.ones(window)/window, mode='same')
    
    # 峰值检测
    peaks, _ = signal.find_peaks(integrated, distance=0.6*fs)
    return peaks

4. 特征工程与模型输入构建

4.1 时域特征提取

典型RR间期特征计算方法：

def extract_rr_features(r_peaks, fs):
    rr_intervals = np.diff(r_peaks)/fs*1000  # 转换为毫秒
    features = {
        'mean_rr': np.mean(rr_intervals),
        'sdnn': np.std(rr_intervals),
        'rmssd': np.sqrt(np.mean(np.square(np.diff(rr_intervals)))),
        'nn50': sum(np.abs(np.diff(rr_intervals))>50),
        'pnn50': sum(np.abs(np.diff(rr_intervals))>50)/len(rr_intervals)*100
    }
    return features

4.2 频域分析实现

使用Welch方法计算功率谱密度：

from scipy import signal

def compute_psd(rr_intervals):
    f, pxx = signal.welch(rr_intervals, fs=4, nperseg=256)
    lf = np.trapz(pxx[(f>=0.04)&(f<0.15)])
    hf = np.trapz(pxx[(f>=0.15)&(f<0.4)])
    return {'lf': lf, 'hf': hf, 'lf_hf_ratio': lf/hf}

5. 完整处理流程示例

整合所有步骤的端到端处理管道：

# 1. 数据加载
record = wfdb.rdrecord('04015', pb_dir='afdb')
ecg = record.p_signal[:,0]

# 2. 信号预处理
denoised = wavelet_denoise(ecg)
detrended = signal.detrend(denoised)

# 3. 特征提取
r_peaks = detect_r_peaks(detrended, fs=record.fs)
rr_features = extract_rr_features(r_peaks, record.fs)
psd_features = compute_psd(np.diff(r_peaks)/record.fs)

# 4. 结果可视化
plt.figure(figsize=(10,6))
plt.plot(detrended[5000:7000])
plt.plot(r_peaks[(r_peaks>5000)&(r_peaks<7000)], 
         detrended[r_peaks[(r_peaks>5000)&(r_peaks<7000)]], 'rx')
plt.title('Processed ECG with R-peaks')
plt.show()