A basic component of processing electrocardiogram (ECG) signals is detecting the heart beat. Beat detection is used to calculate the heart rate, to derive measures of heart rate variability, to develop signal quality indicators, and to detect diseases. There are thousands of publications and strategies to detecting the R-peak of the QRS complex of a heart beat from an ECG signal with varying degrees of accuracy (see section on other beat detection algorithms for a survey). Methods can range from threshold based peak detectors, to wavelet-based signal processing, to probabilistically combining multiple methods.

We will use the Pan-Tompkins algorithm, one of the most widely implemented peak-detection algorithms, to detect the R-peak of the ECG signal in this chapter. Data for this exercise is from the 2017 Physionet Challenge which was aimed to classifying atrial fibrillation from single channel ECG signals. The data was sampled at 300 Hz and band pass filtered. First, we start with a short introduction to ECG wave analysis.

ECG waves

ECG analysis starts with understanding the wave morphology and intervals.

The P-wave reflects atrial deploarization. The amplitude of the P-wave is decreases in diseases like atrial fibrillation, which is a type of arrythmia or abnormal heartbeat. Therefore, we typically want to quantify the amplitude and duration of the P-wave for AFib classification. The distance between the P-wave onset and onset of the QRS complex is the PR interval with a normal duration of 120-220 ms.

The QRS complex reflects depolarization of the left ventricle (since the electrical vector of the left ventricle is much larger than that of the right ventricle). A short QRS duration proves the ventricles are functioning properly and a broad QRS duration indicates that ventricular activation is slow and there could be a dysfunction in the electrical conduction system of the heart. The R-peak of the QRS complex is used to calculate the instantaneous heart rate from the interval between subsequent R-peaks (RR-interval). An RR-interval of 400 ms is equivalent to an instantaneous heart rate of 150 beats per minute (\(60 s / 400 ms * 1000 ms / s\)).

The ST segment is another important morphological feature of the ECG wave since ST elevation and depression are both associated with heart dysfunction like acute myocardial ischemia or ST-elevation myocardial infarction (STEMI). Elevation or depression are calculated as the difference (in millimeters) between the J point (where the ST segment starts) and the PR segment. Finally, the T-wave reflects a repolarization of the contractile cells and is also associated with a range o heart conditions.

Pan-Thompkins algorithm

The Pan-Thompkins algorithm is widely used and can be used for real-time continuous QRS detection. The algorithm is based on analysis of slope, width, and amplitude of ECG using a series of filters. An ECG signal first goes through a bandpass filter, then a differentiator, a squaring operation, a moving window integrator, and finally adaptive thresholding and search-back to find the R-peak.

Raw ECG signals include muscle noise (from respiration), motion artifacts, the QRS complex, and P-T Waves. The band pass filter is designed to match the spectrum of QRS complex, attenuates muscle noise, 60Hz interference, baseline wander, and T-wave interference. Pass band of 5-15 Hz maximizes the QRS energy.

Filtering on waveforms can have negative effects. While low pass filters successfully reduce noise in ECG traces, they also reduce the QRS amplitude. High pass filters (e.g. low cutoff at 0.5 Hz) reduce baseline wander, but also introduce ST distortion. Using forward/backward filtering (highpass, reverse time, highpass, reverse time) removes most of the distortion introduced by high pass filters on the ST segment.

Here we use cascading filters combining a low pass filter and a high pass filter to mimic a bandpass filter. The filter attenuates the P and T waves (which peak at <5 Hz), which is a desired feature since the goal is to detect the QRS complex.

The filtered ECG is then differenced and squared to amplify the QRS complex. The derivative filter further suppresses low frequency components of P and T waves. Squaring makes the signal positive and enhances the derivatives by amplifying the high frequency QRS complex.

Next, the moving average filter over a 150ms window captures the duration of the QRS complex and gives us the integrated signal. This suppresses the smaller oscillations by smoothing out the residual high frequency components. Here we have to define an optimal window length for the moving average. Large windows merge the QRS and T waves together and small windows would produce several peaks at the QRS complex making it difficult to find the R-peak. In addition to detecting the QRS complex the moving average filter gives us the width of the QRS complex.

There are numerous heurestics for peak detection from the integrated signal (e.g. simple thresholding of the moving window integral). Pan-Thompkins proposes to use adpative thresholding and search-back to select a range of time values that correspond to QRS complexes by adapting to changes in ECG by computing running estimates of signal and noise peaks. Here instead I smooth the integrated signal with a gaussian filter to get the energy of the signal. The peak then corresponds to zero-crossings of the first difference where \(x[i+1] < x[i]\).

import numpy as np
import scipy.signal

ecg = np.loadtxt("ecg.txt")  # load the ECG signal
fs = 300  # Hz
tvec = np.arange(len(ecg)) / fs  # time vector

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)

The figure above visualizes each step of the beat detection algorith. The vertical red line are all the zero-crossings (zero_crossings), which lines up with the peak of the ECG signal.

Other beat detection algorithms

Additional resources

• Understanding ECG waves
• Philips DXL ECG Algorithm Physician Guide