❤️

心電図のピーク検出方法　紹介のみ

に公開
Python
ECG
tech

はじめに

情報整理ため,一般的な心電図（ECG）のRピーク取得方法を記載する。
Rピークとは、心電図におけるQRS波のR波部分を指し、心拍のタイミングを示す重要な指標である。

代表的アルゴリズム

  • Pan–Tompkins: 微分・二乗・移動平均・適応閾値を組み合わせた古典的手法。
  • Pan-Tompkins++: より高精度なRピーク検出を実現するための改良版。
  • NeuroKit2（default）: 生理信号処理のためのPythonライブラリ,論文ではないが、NeuroKit2 (JOSS paper) 精度について報告されている.
  • Wavelet Transform: Weblet変換
  • vg-beat-detectors:Emrichら (2023) によるFastNVGアルゴリズム.
  • 他にも様々なアルゴリズムがあるが、一旦保留

実装例（Python：シンプルな流れ）

使用するだなら、neurokit2が便利,method引数でアルゴリズムを選択できるので、他のアルゴリズムも試せます。
neurokit2のデフォルト性能が良いのがわかります。

pip install neurokit2 wfdb



#!/usr/bin/env python3
"""
Simple script to read an ECG record with wfdb, detect R-peaks with neurokit2,
and optionally plot the signal with detected peaks.

Usage:
    python ecg_peaks.py [--record RECORD] [--plot-samples N]

Defaults:
    RECORD: '100' (MIT-BIH arrhythmia database, requires wfdb access)
    N: 2000 (number of initial samples to plot and save; 0 = no plot)
"""
import argparse
import sys
import os
import numpy as np

def main():
    parser = argparse.ArgumentParser()
    parser.add_argument('--record', default='100', help='WFDB record name (default: 100)') # MIT-BIH arrhythmia databaseのID番号
    parser.add_argument('--plot-samples', type=int, default=2000, help='Number of initial samples to plot/save (0 = no plot)')
    args = parser.parse_args()

    try:
        import wfdb
    except Exception as e:
        print('Missing wfdb:', e)
        sys.exit(2)

    try:
        import neurokit2 as nk
    except Exception as e:
        print('Missing neurokit2:', e)
        sys.exit(2)

    # Force a non-GUI backend to avoid Qt plugin errors in headless environments
    try:
        os.environ.setdefault('QT_QPA_PLATFORM', 'offscreen')
        import matplotlib
        matplotlib.use('Agg')
        import matplotlib.pyplot as plt
    except Exception as e:
        print('Warning: matplotlib not available or failed to initialize:', e)
        plt = None

    # Read record from mitdb (requires internet or local WFDB path)
    try:
        record = wfdb.rdrecord(args.record, pn_dir='mitdb')
    except Exception as e:
        print(f'Error reading WFDB record "{args.record}":', e)
        sys.exit(1)

    if record.p_signal is None or record.p_signal.size == 0:
        print('No signal found in record')
        sys.exit(1)

    ecg_signal = record.p_signal[:, 0]
    sampling_rate = getattr(record, 'fs', None)
    if sampling_rate is None:
        print('Warning: record has no fs attribute; using 360 Hz')
        sampling_rate = 360

    # Detect R-peaks
    signals, info = nk.ecg_peaks(ecg_signal, sampling_rate=sampling_rate,method="pantompkins1985") # default method is "neurokit"
    rpeaks = info.get('ECG_R_Peaks')
    if rpeaks is None:
        print('No R-peak information returned')
        sys.exit(1)

    # Normalize rpeaks into an array of peak indices
    try:
        if isinstance(rpeaks, dict):
            # neurokit sometimes wraps results in a dict-like structure
            # try to extract any array-like contained value
            # pick the first array-like value
            vals = None
            for v in rpeaks.values():
                try:
                    arr = np.asarray(v)
                    if arr.size > 0:
                        vals = arr
                        break
                except Exception:
                    continue
            peak_indices = np.asarray(vals, dtype=int) if vals is not None else np.array([], dtype=int)
        else:
            arr = np.asarray(rpeaks)
            if arr.dtype == bool:
                peak_indices = np.where(arr)[0]
            else:
                peak_indices = arr.astype(int)
    except Exception as e:
        print('Error normalizing R-peaks:', e)
        peak_indices = np.array([], dtype=int)

    print(f'Read record: {args.record}')
    print(f'Signal length: {len(ecg_signal)} samples')
    print(f'Sampling rate: {sampling_rate} Hz')
    print(f'Detected R-peaks (count): {len(peak_indices)}')
    if len(peak_indices) > 0:
        print('First 10 R-peak indices:', peak_indices[:10])

    plot_samples = int(args.plot_samples)
    if plot_samples > 0:
        if plt is None:
            print('Matplotlib not available; cannot plot or save image')
        else:
            N = min(plot_samples, len(ecg_signal))
            plt.figure(figsize=(10, 4))
            plt.plot(ecg_signal[:N], label='ECG')
            # use normalized peak_indices to get indices within range
            peak_idx_in_range = peak_indices[peak_indices < N]
            plt.vlines(peak_idx_in_range, ymin=ecg_signal[:N].min(), ymax=ecg_signal[:N].max(), color='r', label='R-peaks')
            plt.legend()
            plt.title(f'ECG record {args.record} (first {N} samples)')
            plt.xlabel('Sample')
            plt.ylabel('Amplitude')
            outname = f'ecg_{args.record}_peaks.png'
            plt.savefig(outname)
            print('Saved plot to', outname)


if __name__ == '__main__':
    main()

python ecg_peaks.py --record 100 --plot-samples 2000

注意点 / Tips

  • サンプリング周波数や電極配置で最適パラメータが異なるため実データでチューニングを行う。
  • 高ノイズ環境ではpan-tompkins法などの古典的な手法では不十分な場合がある。(過検出/見逃し)
  • 複数のアルゴリズムを組み合わせて誤検出を減らす工夫（ポストフィルタリング）を検討する。

参考 (peak検出ライブラあるライブラリ)

  • NeuroKit2: GitHub
  • SciPy にもピーク検出関数がある.

Discussion