Event Detection - R Peaks (ECG)
Difficulty Level:
Tags detect☁ecg☁r-peaks

One of the distinctive characteristics of the electrocardiographic (ECG) signals is their periodicity, something that is not common in physiological terms.

Due to this characteristic, the study of cardiac cycle variability (heart rate variability) defines an important segment of ECG analysis.

However, heart rate variability analysis is dependent of the detection of ECG R peaks, which is the main topic of the present Jupyter Notebook . This task can be achieved by applying the Pan-Tompkins algorithm, translated to Python paradigm by Raja Selvaraj


1 - Importation of the needed packages and definition of auxiliary functions

In [1]:
# biosignalsnotebooks python package
import biosignalsnotebooks as bsnb

# Scientific packages
from numpy import linspace, diff, zeros_like, arange, array

2 - Load of acquired ECG data

In [2]:
# Load of data
data, header = bsnb.load_signal("ecg_4000_Hz", get_header=True)

3 - Identification of the device mac-address and the channel used during acquisition

In [3]:
channel = list(data.keys())[0]
In [4]:
print ("Channel: " + str(channel))
Channel: CH1

4 - Storage of sampling rate and acquired data inside variables

In [5]:
# Sampling rate.
sr = header["sampling rate"]

# Signal Samples.
signal = data[channel]

4.1 - A time-axis should be generated in order to a more intuitive interpretation of the final results

In [6]:
time = bsnb.generate_time(signal)

5 - Simplification of ECG signal (isolation of abrupt transitions)
5.1 - Step 1 of Pan-Tompkins Algorithm - ECG Filtering (Bandpass between 5 and 15 Hz)

In [7]:
# Step 1 of Pan-Tompkins Algorithm
filtered_signal = bsnb.detect._ecg_band_pass_filter(signal, sr)
In [8]:
bsnb.plot_ecg_pan_tompkins_steps(time, signal, filtered_signal, sr, titles=["Original Signal", "Post Filtering Signal"])

5.2 - Step 2 of Pan-Tompkins Algorithm - ECG Differentiation

In [9]:
# Step 2 of Pan-Tompkins Algorithm
differentiated_signal = diff(filtered_signal)
In [10]:
bsnb.plot_ecg_pan_tompkins_steps(time, filtered_signal, differentiated_signal, sr, titles=["Post Filtering Signal", "Post Differentiation Signal"])

5.3 - Step 3 of Pan-Tompkins Algorithm - ECG Rectification

In [11]:
# Step 3 of Pan-Tompkins Algorithm
squared_signal = differentiated_signal * differentiated_signal
In [12]:
bsnb.plot_ecg_pan_tompkins_steps(time, filtered_signal, squared_signal, sr, titles=["Post Differentiation Signal", "Post Rectification Signal"])

5.4 - Step 4 of Pan-Tompkins Algorithm - ECG Integration ( Moving window integration )
5.4.1 - Definition of the samples number inside the moving average window

In [13]:
nbr_sampls_int_wind = int(0.080 * sr)

5.4.2 - Initialisation of the variable that will contain the integrated signal samples

In [14]:
integrated_signal = zeros_like(squared_signal)

5.4.3 - Determination of a cumulative version of "squared_signal"

In the cumulative version of the signal under analysis, his sample value $i$ will be sum of all values included between entry 0 and entry $i$ of the studied signal (in our case "squared_signals").

In [15]:
cumulative_sum = squared_signal.cumsum()

5.4.4 - Estimation of the area/integral below the curve that defines the "squared_signal"

Implicitly, with the current procedure, "squared_signal" is divided into multiple rectangles with fixed width (equal 1 sample) and height determined by the sample value under analysis .

In [16]:
integrated_signal[nbr_sampls_int_wind:] = (cumulative_sum[nbr_sampls_int_wind:] - 
                                           cumulative_sum[:-nbr_sampls_int_wind]) / nbr_sampls_int_wind
integrated_signal[:nbr_sampls_int_wind] = cumulative_sum[:nbr_sampls_int_wind] / arange(1, nbr_sampls_int_wind + 1)
In [17]:
bsnb.plot_ecg_pan_tompkins_steps(time, squared_signal, integrated_signal, sr, titles=["Post Rectification Signal", "Post Integration Signal"])

6 - Application of the ECG simplified version to a sequence of peak identification steps

This task is achieved by using a threshold system, with this threshold being dynamically adjusted along the acquisition, like originally proposed by Pan and Tompkins

6.1 - Initialisation of the R peak detection algorithm

In [18]:
rr_buffer, signal_peak_1, noise_peak_1, threshold = bsnb.detect._buffer_ini(integrated_signal, sr)

6.2 - Detection of possible and probable R peaks

Each sample $i$ of our signal is a possible peak if his value is greater than the ones at sample $i-1$ and $i+1$. On the other hand, a probable peak is a possible peak that meets some criteria defined by Pan and Tompkins at their original article and synthesised inside _detect_peaks function of biosignalsnotebooks Python package.

In [19]:
probable_peaks, possible_peaks= bsnb.detect._detects_peaks(integrated_signal, sr)

Entry "0" of the returned result contains the list of probable peaks, while entry "1" refers to the possible_peaks.

6.3 - Identification of definitive R peaks

Taken into consideration the list of previously detected probable peaks, a set of additional criteria were defined by Pan and Tompkins, in order to exclude peaks from the list of probable peaks.

In [20]:
definitive_peaks = bsnb.detect._checkup(probable_peaks, integrated_signal, sr, rr_buffer, signal_peak_1, noise_peak_1, threshold)

# Conversion to integer type.
definitive_peaks = array(list(map(int, definitive_peaks)))

6.5 - Correcting step

Due to the multiple pre-processing stages there is a small lag in the determined peak positions, which needs to be corrected !

In [21]:
map_integers = definitive_peaks - 40 * (sr / 1000)
definitive_peaks_reph = array(list(map(int, map_integers)))

7 - Evolution of the detected peaks (from possible to definitive R peaks)

In [22]:
bsnb.plot_ecg_pan_tompkins_peaks(time, signal, integrated_signal, sr, possible_peaks, probable_peaks, definitive_peaks)

This procedure can be automatically done by detect_r_peaks function in detect module of opensignalstools package

In [23]:
detected_peaks = bsnb.detect_r_peaks(signal, sr, time_units=True, plot_result=True)

As demonstrated, R peak detection algorithm contains a great number of stages, obviously there are alternative and simpler methodologies. However, to obtain more precise results, a bigger complexity becomes mandatory !

Now you can understand better the complexity behind the heart beat monitoring systems but you can also program it.

We hope that you have enjoyed this guide. biosignalsnotebooks is an environment in continuous expansion, so don"t stop your journey and learn more with the remaining Notebooks !

In [24]:
from biosignalsnotebooks.__notebook_support__ import css_style_apply
css_style_apply()
.................... CSS Style Applied to Jupyter Notebook .........................
Out[24]: