Pan-Tompkins : A Real-Time QRS Detection Algorithm

ECG beat detection

 

심전도 (ECG) 신호 처리의 기본 구성 요소는 심장 박동을 검출하는 것이다. 비트 검출은 심장박동을 계산하고, 심박수 변동성의 측정치를 도출하고, 신호 품질 지표를 개발하고, 질병을 검출하는 데 사용된다. 다양한 정확도를 가진 ECG 신호에서 심장 박동의 QRS 복합의 R-피크를 검출하기 위한 수천 가지 전략이 있다. 방법은 임계값 기반 피크 검출기에서 웨이블릿 기반 신호 처리, 여러 가지 방법을 확률적으로 결합할 수 있다. ECG 신호의 R-피크를 검출하기 위해 가장 널리 구현된 피크 검출 알고리즘 중 하나는 Pan-Tompkins 알고리즘이다.

 

ECG 파형

 

심전도 분석은 파형 형태와 간격을 이해하는 것으로 시작한다.

 

https://coffeymedical.com/wp-content/uploads/2017/08/DXL_Algorithm_White_Paper.pdf

 

P파는 심방 탈포화를 반영한다. P파의 진폭은 관절통 또는 비정상적인 심장 박동의 일종인 심방 세동과 같은 질병에서 감소한다. 따라서, 일반적으로 AFib 분류를 위해 P파의 진폭과 지속 시간을 정량화하려고 한다.

 

QRS complex의 P파 시작과 시작 사이의 거리는 정상 지속 시간이 120-220ms인 PR 간격이다. QRS complex는 좌심실의 탈편극을 반영한다 (좌심실의 전기 벡터가 우심실의 전기 벡터보다 훨씬 크기 때문). 짧은 QRS 지속 시간은 심실이 제대로 작동하고 있음을 증명하고 광범위한 QRS 지속 시간은 심실 활성화가 느리고 심장의 전기 전도 시스템에 기능 장애가 있을 수 있음을 나타낸다. QRS complex의 R-피크는 후속 R-피크 사이의 간격 (RR 간격)에서 순간 심박수를 계산하는 데 사용된다. 400ms의 RR 간격은 분당 150 비트의 순간 심박수와 동일하다.

 

ST 상승과 우울증은 모두 급성 심근 허혈 또는 ST-고혈압 심근경색 (STEMI)과 같은 심장 기능 장애와 관련이 있기 때문에 ST 세그먼트는 ECG 파형의 또 다른 중요한 형태학적 특징이다. 상승 또는 우울증은 J 지점 (ST 세그먼트가 시작되는 곳)과 PR 세그먼트 사이의 차이 (밀리미터)로 계산된다. 마지막으로, T파는 수축 세포의 재편극을 반영하며 또한 범위 o 심장 상태와 관련이 있다.

 

Pan-Tompkins

 

Pan-Thompkins 알고리즘은 널리 사용되며 실시간 연속 QRS 검출에 사용될 수 있다. 알고리즘은 일련의 필터를 사용하여 심전도의 기울기, 너비 및 진폭 분석을 기반으로 한다. 심전도 신호는 먼저 대역 통과 필터를 통과한 다음 미분기, 제곱 연산, moving window 적분기를 통과하고 마지막으로 적응형 임계값 및 search-back을 통해 R-피크를 찾는다.

 

 

원시 심전도 신호에는 근육 소음 (호흡으로 인한), 운동 인공물, QRS complex 및 P-T파가 포함된다. 대역 통과 필터는 QRS complex 스펙트럼과 일치하도록 설계되었으며, 근육 소음, 60Hz 간섭, 기준선 방랑 및 T파 간섭을 완화한다. 5-15Hz의 통과 대역은 QRS 에너지를 최대화한다.

 

파형에 대한 필터링은 부정적인 영향을 미칠 수 있다. 저역 통과 필터는 성공적으로 ECG 추적의 노이즈를 감소시키는 반면 QRS 진폭도 감소시킨다. 고역 통과 필터 (ex: 0.5Hz에서 낮은 차단)는 기저선 변동을 감소시키지만 ST 왜곡을 도입하기도 한다. 정방향 / 역방향 필터링 (하이패스, 역방향 시간, 하이패스, 역방향 시간)을 사용하면 ST 세그먼트에서 고역 통과 필터에 의해 유입되는 왜곡을 대부분 제거할 수 있다.

 

여기서 대역 통과 필터를 모방하기 위해 저역 통과 필터와 고역 통과 필터를 결합한 캐스케이드 필터를 사용한다. 필터는 P파와 T파 (5Hz 미만에서 피크)를 감쇠시키며, 이는 QRS complex를 검출하는 것이 목표이기 때문에 원하는 특징이다.

 

그런 다음 필터링된 ECG를 미분하고 제곱하여 QRS 복합을 증폭한다. 유도 필터는 P파와 T파의 저주파 성분을 추가로 억제한다. 제곱은 신호를 양수로 만들고 고주파 QRS complex을 증폭하여 도함수를 향상시킨다.

 

다음으로 150ms 윈도우를 통한 이동 평균 필터는 QRS complex의 지속 시간을 캡처하고 통합 신호를 제공한다. 이는 잔류 고주파 성분을 평활화하여 더 작은 진동을 억제한다. 여기서 이동 평균에 대한 최적의 윈도우 길이를 정의해야 한다. 큰 윈도우는 QRS와 T파를 함께 병합하고 작은 윈도우는 R-피크를 찾기 어렵게 하는 QRS complex에서 여러 개의 피크를 생성한다. QRS complex를 검출하는 것 외에도 이동 평균 필터는 QRS complex의 폭을 제공한다.

 

적분 신호에서 피크 검출을 위한 수많은 휴리스틱 (ex: 이동 윈도우 적분의 단순 임계)이 있다. Pan-Thompskins는 신호 및 노이즈 피크의 추정치를 계산하여 ECG의 변화에 적응하여 QRS complex에 해당하는 시간 값의 범위를 선택하기 위해 적응적 임계 및 search-back을 사용할 것을 제안한다. 여기서 대신 신호의 에너지를 얻기 위해 가우스 필터로 통합 신호를 평활한다. 피크는 그러면 첫 번째 차이의 제로 크로싱에 해당한다.

 

import sys
import importlib 
from pathlib import Path
import string

import numpy as np
import pandas as pd
import scipy.signal

# plotting imports
import matplotlib.pyplot as plt
import seaborn as sns
sns.reset_defaults()
plt.rc('font', family='Helvetica')
%matplotlib inline
%config InlineBackend.figure_format = 'retina'

# add parent folder to sys.path so we can import
# config and utils modules
if str(Path.cwd().parent) not in sys.path:
    sys.path.append(str(Path.cwd().parent))

import config
import utils

importlib.reload(config)
importlib.reload(utils)

# load data
fs = 300  # Hz
data = np.load(config.PHYSIONET_2017_ECG, allow_pickle=True).item(0)
record_ids = list(data.keys())

# select a single record
ecg = data[record_ids[0]]
tvec = np.arange(len(ecg)) / fs

max_QRS_duration = 0.150  # sec
low_cutoff = 5
high_cutoff = 15
window_size = int(max_QRS_duration * fs)

# apply a bandpass filter to the ECG signal
lowpass = scipy.signal.butter(1, high_cutoff / (fs / 2.0), "low")
highpass = scipy.signal.butter(1, low_cutoff / (fs / 2.0), "high")
ecg_low = scipy.signal.filtfilt(*lowpass, x=ecg)
ecg_band = scipy.signal.filtfilt(*highpass, x=ecg_low)

diff = np.diff(ecg_band)
squared = np.square(diff)

# moving average filter
# apply padding on both sides of the signal and convolve to get the integrated signal
mwa = np.pad(squared, (window_size - 1, 0), "constant", constant_values=(0, 0))
mwa = np.convolve(mwa, np.ones(window_size), "valid")
for i in range(1, window_size):
    mwa[i - 1] = mwa[i - 1] / i
mwa[window_size - 1 :] = mwa[window_size - 1 :] / window_size
mwa[: int(max_QRS_duration * fs * 2)] = 0

# smooth the moving window integrated signal with a gaussian filter and take the derivative
energy = scipy.ndimage.gaussian_filter1d(mwa, fs / 8.0)
energy_diff = np.diff(energy)

# peaks are the points where the derivative crosses zero, adjust window size
zero_crossings = (energy_diff[:-1] > 0) & (energy_diff[1:] < 0)
zero_crossings = np.flatnonzero(zero_crossings)
zero_crossings -= int(window_size / 2)