QRS peak detection#
Detection of cardiac cycles from ECG signals
# setup
import sys
import numpy as np
import scipy.signal as sp
from matplotlib import pyplot as plt
import wfdb
Download records#
Identify and download records in the MIMIC III Waveform Database
# Select the first record
selected_record = '3000063_0013'
database_name = 'mimic3wdb/1.0/30/3000063/'
print("Selected record: {}".format(selected_record))
# load data from this record
start_seconds = 20
no_seconds_to_load = 10
fs = 125
record_data = wfdb.rdrecord(record_name = selected_record, sampfrom = fs*start_seconds, sampto = fs*(start_seconds + no_seconds_to_load), pn_dir = database_name)
print("{} seconds of data loaded from: {}".format(no_seconds_to_load, selected_record))
# Plot the data loaded from this record
title_text = "First record from " + database_name + " (record: " + selected_record + ")"
wfdb.plot_wfdb(record=record_data, title=title_text, time_units='seconds')
Selected record: 3000063_0013
10 seconds of data loaded from: 3000063_0013
Separate records#
Separate ECG, PPG and ABP signals
# ECG, PPG and ABP extraction
h = 0
for i in record_data.sig_name:
if i == 'II':
print('ECG at position ' + str(h))
ecg = record_data.p_signal[:,h]
elif i == 'PLETH':
print('PPG at position ' + str(h))
ppg = record_data.p_signal[:,h]
elif i == 'ABP':
print('ABP at position ' + str(h))
abp = record_data.p_signal[:,h]
else:
print('Other signal (' + i + ') at position ' + str(h))
h = h + 1
t = np.arange(start_seconds,(start_seconds + no_seconds_to_load),1.0/fs)
fig, (ax1,ax2,ax3) = plt.subplots(3, 1, sharex = True, sharey = False)
ax1.plot(t, ecg)
ax1.set(xlabel = '', ylabel = 'ECG')
ax2.plot(t, ppg)
ax2.set(xlabel = '', ylabel = 'PPG')
ax3.plot(t, abp)
ax3.set(xlabel = 'Time [s]', ylabel = 'ABP')
ECG at position 0
ABP at position 1
PPG at position 2
[Text(0.5, 0, 'Time [s]'), Text(0, 0.5, 'ABP')]
# Remove NaNs
x, = np.where(np.isnan(ecg))
if len(x) != 0:
print('NaNs in ECG: ' + str(x))
for i in x:
ecg[i] = ecg[i - 1]
x, = np.where(np.isnan(ppg))
if len(x) != 0:
print('NaNs in PPG: ' + str(x))
for i in x:
ppg[i] = ppg[i - 1]
x, = np.where(np.isnan(abp))
if len(x) != 0:
print('NaNs in ABP: ' + str(x))
for i in x:
abp[i] = abp[i - 1]
Filter data#
# Filter ECG
sos_ecg = sp.butter(10, [0.5, 40], btype = 'bp', analog = False, output = 'sos', fs = fs)
ecg_ff = sp.sosfiltfilt(sos_ecg, ecg)
# Filter PPG
sos_ppg = sp.butter(10, [0.2, 10], btype = 'bp', analog = False, output = 'sos', fs = fs)
ppg_ff = sp.sosfiltfilt(sos_ppg, ppg)
# Filter ABP
sos_abp = sp.butter(10, 10, btype = 'low', analog = False, output = 'sos', fs = fs)
abp_ff = sp.sosfiltfilt(sos_abp, abp)
fig, (ax1,ax2,ax3) = plt.subplots(3, 1, sharex = False, sharey = False)
fig.suptitle('Filtered signals')
ax1.plot(t, ecg, color = 'red')
ax1.plot(t, ecg_ff, color = 'orange')
ax1.set_xlabel('Time [s]')
ax1.set_ylabel('ECG [V]')
ax2.plot(t, ppg, color = 'red')
ax2.plot(t, ppg_ff, color = 'orange')
ax2.set_xlabel('Time [s]')
ax2.set_ylabel('PPG [V]')
ax3.plot(t, abp, color = 'red')
ax3.plot(t, abp_ff, color = 'orange')
ax3.set_xlabel('Time [s]')
ax3.set_ylabel('ABP [V]')
plt.subplots_adjust(left = 0.1,
bottom = 0.1,
right = 0.9,
top = 0.9,
wspace = 0.4,
hspace = 0.4)
#plt.show()
QRS detection#
Define QRS detector functions#
# QRS detection function
def qrs_detect(ecg,fs,w):
"""
Description: QRS peak detection and correction
Inputs: ecg, array with ECG signal [user defined units]
fs, sampling rate of signal [Hz]
w, window length for analysis [s]
Outputs: qrs, positions of R peaks [number of samples]
n_int, number of intervals of length w in the signal
Libraries: NumPy (as np), SciPy (Signal, as sp), Matplotlib (PyPlot, as plt)
Version: 1.0 - June 2022
Developed by: Elisa Mejía-Mejía
City, University of London
"""
# Normalization of ECG
ecg_n = (ecg - np.min(ecg))/(np.max(ecg) - np.min(ecg))
ecg_n = ecg_n - np.mean(ecg_n)
# Peak detection in windows of length w
n_int = np.floor(len(ecg)/(w*fs))
for i in range(int(n_int)):
start = i*fs*w
stop = (i + 1)*fs*w - 1
#print('Start: ' + str(start) + ', stop: ' + str(stop))
aux = ecg_n[range(start,stop)]
locs, amps, nlocs = rpeakdetect(aux,fs,0.2)
locs = locs + start
if i == 0:
qrs = locs
else:
qrs = np.append(qrs,locs)
if n_int*fs*w != len(ecg_n):
start = stop + 1
stop = len(ecg_n)
aux = ecg_n[range(start,stop)]
if len(aux) > 20:
locs, amps, nlocs = rpeakdetect(aux,fs,0.2)
locs = locs + start
qrs = np.append(qrs,locs)
# Peak correction
# (1) Not peaks
med_rr = np.median(np.diff(qrs))
med_amp = np.median(ecg[qrs])
pos_for = 0
pos_rev = 0
for i in range(len(qrs)):
cur_amp = ecg[qrs[i]]
if i == 1: # Peak is at the first position of the signal
next_amp = ecg[qrs[i] + 1]
cond = cur_amp >= next_amp
#print(cur_amp, next_amp, cond)
elif i == len(ecg): # Peak is at the last position of the signal
prev_amp = ecg[qrs[i] - 1]
cond = cur_amp >= prev_amp
#print(cur_amp, prev_amp, cond)
else:
next_amp = ecg[qrs[i] + 1]
prev_amp = ecg[qrs[i] - 1]
cond = cur_amp >= prev_amp and cur_amp >= next_amp
#print(cur_amp, prev_amp, next_amp, cond)
if not(cond):
# Forward search
j = qrs[i] + 1
if i == len(qrs) - 1:
stop = len(ecg)
else:
stop = qrs[i + 1]
while j < stop:
cur_amp = ecg[j]
next_amp = ecg[j + 1]
prev_amp = ecg[j - 1]
if not(cur_amp > next_amp and cur_amp > prev_amp):
j = j + 1
else:
if ecg[j] >= 0.5*med_amp:
pos_for = j
j = stop
else:
j = j + 1
# Reverse search
j = qrs[i] - 1
if i == 0:
stop = 1
else:
stop = qrs[i - 1]
while j > stop:
cur_amp = ecg[j]
next_amp = ecg[j + 1]
prev_amp = ecg[j - 1]
if not(cur_amp > next_amp and cur_amp > prev_amp):
j = j - 1
else:
if ecg[j] >= 0.5*med_amp:
pos_rev = j
j = stop
else:
j = j - 1
# Selection of peak
if pos_for != 0 and pos_rev != 0:
pos = [pos_for, pos_rev]
if i - 1 == 0:
dif = - pos
else:
dif = qrs[i - 1] - pos
dif_med = np.abs(dif - med_rr)
#print(dif_med)
min_dif = np.min(dif_med)
min_dif, = np.where(dif_med == min_dif)
#print(min_dif,pos(min_ind))
qrs[i] = pos(min_ind)
# (2) Low peaks
med_amp = np.median(ecg[qrs])
#print(len(qrs))
result, = np.where(ecg[qrs] < 0.5*med_amp)
if len(result) > 0:
#print(result)
qrs = np.delete(np.array(qrs),result)
#print(len(qrs))
# (3) Too close or too far
rr = np.diff(qrs)
avg_rr = np.mean(rr)
med_rr = np.median(rr)
q3_rr, q1_rr = np.percentile(rr,[75, 25])
iqr_rr = q3_rr - q1_rr
#print(rr, avg_rr, med_rr, q3_rr, q1_rr, iqr_rr)
if (med_rr < 0.5*fs or avg_rr < 0.5*fs) and (iqr_rr > 0.2*fs):
result, = np.where(rr <= med_rr)
qrs = np.delete(qrs, result)
med_rr = np.median(rr)
# (a) Too close
ind, = np.where(rr < 0.5*med_rr)
if len(ind) != 0:
i = 1
while i <= len(ind):
ind = no.sort(np.append(ind, ind[i] + 1))
#print(ind)
i = i + 2
qrs = np.delete(qrs, ind)
# (b) Too far
rr = np.diff(qrs)
ind, = np.where(rr > 1.5*med_rr)
th_amp = 0.5*np.median(ecg[qrs])
#print(ind, th_amp)
if len(ind) != 0:
i = 0
while len(ind) != 0:
#ind = np.round(ind)
#qrs = np.round(qrs)
#print(i, ind[i], len(qrs))
if qrs[ind[i]] == qrs[-1]:
aux = ecg[qrs[ind[i]]:len(ecg)]
else:
aux = ecg[qrs[ind[i]]:qrs[ind[i] + 1]]
locs, _ = sp.find_peaks(aux)
ind_del, = np.where(aux[locs] < th_amp)
locs = np.delete(locs, ind_del)
if len(locs) == 0:
new = med_rr + qrs[ind[i]] - 1
elif len(locs) == 1:
new = locs + qrs[ind[i]] - 1
else:
aux_dif = np.abs(locs - med_rr)
min_dif = np.min(aux_dif)
ind_new, = np.where(aux_dif == min_dif)
new = locs[ind_new] + qrs[ind[i]] - 1
qrs = np.sort(np.append(qrs, new))
rr = np.diff(qrs)
ind, = np.where(rr > 1.5*med_rr)
rr = np.diff(qrs)
ind_del, = np.where(rr < 0.3*fs)
rr = np.delete(rr, ind_del)
#fig = plt.figure()
#plt.plot(ecg, color = 'black')
#plt.scatter(qrs,ecg[qrs],marker = 'o',color = 'red')
return qrs, n_int
# R peak detect
def rpeakdetect(ecg, fs, th):
"""
Description: QRS peak detection based on rpeakdetect2.m by G. Clifford
A batch QRS detector based upon that of Pan, Hamilton and Tompkins:
J. Pan & W. Tompkins. A real-time QRS detection algorithm. IEEE Trans Biomed Eng, vol. BME-32 NO. 3. 1985.
P. Hamilton & W. Tompkins. Quantitative Investigation of QRS Detection Rules Using the MIT/BIH Arrythmia
Database. IEEE Trans Biomed Eng, vol. BME-33, NO. 12.1986.
Inputs: ecg, array with ECG signal [user defined units]
fs, sampling rate of signal [Hz]
th, threshold for peaks in integrated waveform - Default: 0.2
Outputs: locs, positions of R peaks [number of samples]
amps, amplitudes of R peaks [user defined units]
nlocs, number of R peaks found
Libraries: NumPy (as np), SciPy (Signal, as sp), Matplotlib (PyPlot, as plt)
Version: 1.0 - June 2022
Developed by: Elisa Mejía-Mejía
City, University of London
"""
# Create time array
t = np.divide(range(0,len(ecg) - 1), fs)
# Preprocessing:
# (1) Filter data
sos = sp.butter(6, 40, btype = 'low', analog = False, output = 'sos', fs = fs)
ecg_f = sp.sosfiltfilt(sos, ecg)
# (2) Differentiate ECG
ecg_d = np.diff(ecg_f)
# (3) Square ECG
ecg_s = ecg_d*ecg_d
# (4) Integrate data
if fs >= 256:
d = np.ones(int(np.round(7*fs/256)))
else:
d = np.ones(21)
ecg_fir = sp.lfilter(d, 1, ecg_s)
ecg_med = sp.medfilt(ecg_fir, kernel_size = 9)
# Remove filter delay:
delay = np.ceil(len(d)/2)
ecg_med = ecg_med[int(delay):len(ecg_med)]
# Find peaks:
# (1) Find highest bumps in data
start_int = round((len(ecg) - 1)/4)
stop_int = round(3*(len(ecg) - 1)/4)
max_h = np.max(ecg_med[start_int:stop_int])
# (2) Determine left and right boundaries
ecg_med = np.insert(ecg_med,0,0)
ecg_med = np.append(ecg_med,0)
n_left = 0
n_right = 0
for i in range(len(ecg_med)):
if i > 0 and i < len(ecg_med):
if ecg_med[i] > (th*max_h) and ecg_med[i - 1] < (th*max_h): # left boundary
if n_left == 0:
left_bound = i - 1
n_left = 1
else:
left_bound = np.append(left_bound,i - 1)
if ecg_med[i] > (th*max_h) and ecg_med[i + 1] < (th*max_h): # right boundary
if n_right == 0:
right_bound = i - 1
n_right = 1
else:
right_bound = np.append(right_bound,i - 1)
#print(left_bound)
#print(right_bound)
# (4) Look through all possibilities
if left_bound[0] > right_bound[0]:
right_bound = np.delete(right_bound,0)
if left_bound[-1] > right_bound[-1]:
left_bound = np.delete(left_bound,0)
nlocs = 0
for i in range(len(left_bound)):
#print(i)
lb = left_bound[i]
rb, = np.where(np.array(right_bound) > lb)
if len(rb) != 0:
rb = right_bound[rb[0]]
#print(lb, rb)
amp = np.max(ecg[lb:rb])
pos, = np.where(np.array(ecg[lb:rb]) == np.amax(ecg[lb:rb]))
pos = pos[0]
pos = int(pos + lb)
#print(pos, amp)
if nlocs == 0:
locs = pos
amps = ecg[pos]
else:
locs = np.append(locs, pos)
amps = np.append(amps, ecg[pos])
nlocs = nlocs + 1
#print(locs)
#fig = plt.figure()
#plt.plot(ecg)
#plt.plot(ecg_f)
#plt.plot(ecg_d)
#plt.plot(ecg_s)
#plt.plot(ecg_fir)
#plt.plot(ecg_med)
#plt.scatter(locs, amps, marker = 'o')
return locs, amps, nlocs,
# Detect cardiac cycles
qrs, n_int = qrs_detect(ecg_ff,fs,5)
fig = plt.figure()
plt.title('QRS detection')
plt.plot(t, ecg_ff, color = 'black')
plt.scatter(start_seconds + qrs/fs, ecg_ff[qrs], color = 'red', marker = 'o')
<matplotlib.collections.PathCollection at 0x7fb7ad332950>
