Improving the QRS detection for one-channel ECG sensor

Abstract

We analyzed several QRS detection algorithms in order to build a quality industrial beat detector, intended for a small, wearable, one channel electrocardiogram sensor with a sampling rate of 125 Hz, and analog-to-digital conversion of 10 bits. The research was a lengthy process that included building several hundred rules to cope with the QRS detection problems and finding an optimal threshold value for several parameters. We obtained 99.90% QRS sensitivity and 99.90% QRS positive predictive rate measured on the first channel of rescaled and resampled MIT-BIH Arrhythmia ECG database. Even more so, our solution works better than the algorithms for the original signals with a sampling rate of 360 Hz and analog-to-digital conversion of 11 bits.

Keywords

QRS detection beat detector electrocardiogram wearable sensors rescale resample optimization

1. Introduction

Advances in the Internet of Things (IoT) field have encouraged researchers to intensify their focus on Electrocardiogram (ECG) processing, especially for wearable devices. This field has become a popular research topic in biomedical engineering [1]. Our primary focus is, on building a quality industrial QRS detector for a wearable ECG sensor that is significantly better than existing robust algorithms for QRS detection, which are inconvenient for mobile devices with limited resources [2].

The performance of the QRS detector is evaluated by calculating how many real QRS peaks are found (QRS sensitivity, denoted as $Q_{\textit{SE}}$ ), and how many of those detected QRS peaks are real beats (QRS positive predictive rate, denoted as $Q_{+P}$ ). The default testing database was the Massachusetts Institute of Technology – Beth Israel Hospital Arrhythmia database (MITDB) [3] with 48 records of 30 minute ECG measurements. The original signals are sampled with 360 Hz and 11 bit AD conversion. Our research target is a QRS detector for signals using a sampling frequency of 125 Hz and 10 bit AD conversion.

One of the most cited papers for QRS detection built for small devices with limited resources is the Pan and Tomkins algorithm [4]. Its robustness lies in the fact that it runs fast enough to be used in real time and can cope with noisy signals. However, the performance of this algorithm depends on bit resolution in the AD conversion. In our case, when using one channel ECG sensor and smaller sampling frequencies, its performance was not satisfactory, especially for signals with smaller amplitudes. It was not a good solution for us because of these factors.

Another alternative is the Hamilton algorithm [5]. Compared to the Pan and Tompkins algorithm, it is quite similar but uses different filters and decision rules. It is a stable solution, but, still is unable to cope with small amplitudes, or variations in consecutive amplitude levels, especially, when using a smaller bit resolution in AD conversion.

Physionet.org [6] is a very comprehensive resource where one can find several QRS detection algorithms, including Wavedet, gqrs, wqrs, and sqrs. They represent simple and fast algorithms demanding a small number of resources and obtain high sensitivity and positive predictive rate values. However, they lack the beat classification, and the obtained sensitivity and the positive predictive rate are also considered to be lower than the demands of a quality industrial QRS detector for a wearable one channel ECG sensor.

After these initial efforts, the attention of researchers gradually focused on developing more sophisticated QRS detection algorithms, including Machine Learning and other methods, as described in Section 2 (related work). Although some of the new approaches achieved better performance, they generally require computationally intensive algorithms, not suitable for smartphones that collect continuous ECG data from wearable sensors.

In our research, we improve Hamilton’s algorithm [5] in order to make it efficient for industrial application. The improvement was a rather long process due to the exponential nature of the effort to improve the algorithm. The closer you are to the margin of 100%, the more effort is needed for a very small improvement of the performance. We introduced several hundred rules to cope with the identified problems in QRS detection, and several thousand tests to tune parameters, and threshold values for identified solutions. Some threshold values obtained good performance on some test file while performing badly on others. When we fine-tuned some parameters, it so happened that some of the rules did not work on other test datasets, which was even more challenging.

Finally, the results proved that we have realized our goals. We reached higher values than all other published results, though we were working with almost three times less sampling frequency and half of the AD conversion bit resolution.

The structure of this paper is organized as follows. Section 2 discusses the related work in the area of QRS detection algorithms and Section 3 presents the background. Analysis of problems are presented in Section 4, and our approach to improving the algorithms in Section 5. The results from the experiments are evaluated in Section 6 and compared with other solutions. Section 7 is devoted to conclusions and future work.

2. Related work

QRS detection algorithms generally follow the same routine [7], starting with Digital Signal Processing (DSP) filtering that eliminates noise and baseline wander. Then the output is matched to a set of thresholds. The algorithms mainly differ in the way they calculate the threshold values and implement decision rules, which is another layer implemented after thresholding [8].

Long ECG recordings are generally occupied with noise in the form of a subtle deviation of the heart rhythm. In addition, signal quality changes to alternating changes in the wave amplitude. Unless noise is eliminated, detection of such QRS complexes becomes difficult, and decreases the accuracy [9]. It is important to state that any algorithm designed for quality industrial QRS detection must be adaptable to any type of noise.

There are several survey papers [9, 10] which give a comprehensive overview of popular QRS detection methods. In summary, published QRS detection algorithms are based on the following techniques:

•
Differentiation (derivation), where the difference between the current and previous samples is calculated, as a way of identifying the slope, and then it is compared to a given threshold value including the Pan & Tompkins algorithm [4], Hamilton’s algorithm [5], or other relevant approaches [11, 12, 13];
•
Pure DSP algorithms, where fundamental DSP filters with different characteristics are combined to produce a bandpass filter, with the aim to eliminate noise, and filter the signal so that a threshold will determine the beats [14, 15, 16, 17, 18];
•
Pattern Recognition algorithms, where the signal data is matched with predefined patterns and a waveform is detected in case of similarity within given constraints of the amplitude and slopes [19, 20, 21, 22, 23];
•
Neural Network, Multilayer Perceptron (MLP), Radial Basis Function (RBF), and Learning Vector Quantization (LVO) networks are used to adaptively predict the location of the next peak [24, 25, 26, 27, 21];
•
Digital Wavelet Transformation (DWT), where the signal is decomposed to a certain level of scales, and then recomposed, which effectively reduces noise. Then a threshold is applied to select proper peaks [28, 29, 30, 31, 32];
•
Genetic Algorithms are used to optimize the preprocessing polynomial filter. The ECG signal is compared to an adaptive threshold and the parameters are optimized with a genetic optimization approach [33];
•
Hidden Markov Models (HMM), used to train the probability function varying according to the hidden Markov chain, and then the model predicts the current state, which could be QRS complex. P and T waves can also be computed [34, 35, 36];
•
Hilbert Transform, where Hilbert transform of ECG signal is calculated by Fast Fourier Transform (FFT), and that is used for calculating the signal envelope [37, 38, 39]; and
•
Phasor Transformation, where each ECG sample is converted into a phasor to correctly manage P, and T waves, by definition having lower amplitudes than an R-peak with low computational cost, and then compares it against a threshold [40].

High performance QRS detector directly affects the amount and the quality of valuable information on ECG. QRS detection is the initial step for further ECG analysis.
3. Background

In this section, we will explain the evaluation metrics and give an overview of the original Hamilton’s QRS detector algorithm.

3.1 Performance measures

The benchmarks used in our testing methodology are the same used in the IEC 60601-2-47 standard for particular requirements for the safety, including an essential performance of ambulatory electrocardiographic systems, and ANSI/AAMI EC57:2012 for Testing and Reporting Performance Results of Cardiac Rhythm and ST Segment Measurement Algorithms. These standards use the MIT-BIH ECG arrhythmia database [3], and the American Heart Association’s (AHA database) [41].

MITDB contains half-hour ECG recordings for 48 anonymized persons, and only 44 records exclude those that contain paced beats. These recordings are publicly available on the physionet.org web site [6]. The recording frequency is 360 samples per second, per channel, with a 11-bit resolution. Even though each recording contains two-channels, we used the first channel, identified as ML II, in most of the records.

In addition, we follow the requirements according to the standard IEC 60601-2-47:2012 for medical electrical equipment, particularly, requirements for essential performance of ambulatory electrocardiographic systems. According to these requirements, any calculated peak is considered as detected if it is at most 150 ms away from the real beat.

A detected QRS is denoted to be True Positive (TP), if the QRS detector has found a QRS closer than 150 ms from the one which is annotated. A False Negative (FN) is a missed QRS, or if the QRS detector has found a QRS outside the 150 ms perimeter, while a False Positive (FP) is an erroneously detected QRS (extra found).

The commonly used performance measures are sensitivity and positive predictive rate, calculated by Eq. (1).

$\displaystyle\textit{SE}=\frac{\textit{TP}}{\textit{TP}+\textit{FN}}+P=\frac{% \textit{TP}}{\textit{TP}+\textit{FP}}$ (1)

In addition, to find an optimal value of a parameter, we provided a lot of test experiments, and calculate the number of Total Errors by Eq. (2), as a sum of errors denoted by FP and FN. The smaller the errors are, the better performance is achieved.

$\displaystyle\textit{Errors}=\textit{FP}+\textit{FN}$ (2)

3.2 Analysis of Hamilton’s algorithm

EP Limited’s open source software for arrhythmia detection serves as a basis for this research [42, 5]. It has a complete C-code implementation of Hamilton’s algorithm, with three different detectors, and a simple beat classifier. Two of the detectors are for general-purpose, whereas the third one is for environments with a small amount of memory.

Algorithmic details are theoretically provided in their original work [5]. Figure 1 represents the conceptual level for the two phases, and the high-level steps conducted for each of them.

Figure 1.

Hamilton’s approach for QRS detection.

The algorithm starts with a low pass filter (LPF) with a cutoff frequency of 16 Hz. Then, the signal is passed through a high pass filter (HPF) of 8 Hz. Both filters have the effect of a band pass filter (BPF). Then, the slope of the signal is calculated with a differentiation method ( $\frac{d[i]}{dt}$ ), followed by the calculation of an absolute value (ABS). The final step (AVG) consists of calculating a time average over an 80 ms moving window.

After eliminating the noise in the DSP filtering phase, the algorithm continues with the peak detection phase. It already has two thresholds for the AVG signal, classified as a

•

static threshold with a fixed value, and

•

dynamic adaptive threshold (DAT) which is affected by the amplitudes of the latest peaks.

The original algorithm sets the static threshold (STHR) at value $\textit{MIN\_PEAK\_AMP}=$ 7. A general rule of thumb is that a lower value of the static threshold will find more peaks, but also detect lots of artifacts. On the other hand, a higher static threshold value yields fewer peaks, but a smaller number of artifacts.

When a new local peak is found, the dynamic adaptive threshold is calculated by taking the mean values for real peaks and noise peaks into account. The mean value is computed by Eq. (3).

$\displaystyle\textit{mean }=\frac{\sum_{n=1}^{8}X_{n}}{8}$ (3)

Let the mean value for real beats and noise peaks be denoted as qmean and nmean respectively, and also TH be the constant multiplier (with a default value of 0.3125), then the DAT is calculated by Eq. (4).

$\displaystyle\textit{DAT}=\textit{nmean}+(\textit{qmean-nmean})*\textit{TH}$ (4)

When a new local maximum is detected by calculating an AVG value, both of the thresholds are compared to this value. If the value is higher than the static threshold, then it is considered to be a potential peak, otherwise, it is an artifact. It is classified as a noise peak if the calculated value is lower than the dynamic peak, and as a real beat if it is higher than the dynamic adaptive threshold.

Figure 2 presents both dynamic and static thresholds. Note, that detected peaks A1, A2, $\ldots$ , A7 are categorized as artifacts (smaller than the static threshold), and R1, R2, R3 and R4 as real beats. N1, N2, and N3 are considered noise peaks since the local maxima of each label are smaller than the dynamically calculated threshold.

Figure 2.

Detecting artifacts, noise and real peaks based on values of static and dynamic thresholds in the original Hamilton’s algorithm presented on signal extract over MITDB record 124 (1046 sec).

4. Identification of performance issues

A discrepancy in peak amplitudes may introduce bad detection. We have identified two cases when this happens:

•
a sequence of low-amplitude peaks after an isolated high-amplitude peak;
•
an isolated low-amplitude peak after a sequence of high-amplitude peaks.

Furthermore, apart from fine-tuning the threshold value, and expecting lower performance on cropped signals, we analyzed the segments where the algorithm showed lower sensitivity, and specificity even though the signal was not contaminated with noise. The conclusion was that lower performance was obtained for specific segments, and the problems can be classified as:

•
a mixture of low and high-amplitude peaks;
•
artifact elimination; and
•
wrong R-peak location.

4.1 Bad detection of low-amplitude peaks

Figure 3 presents the case when an isolated high-amplitude peak is preceded and followed by a sequence of low-amplitude peaks. An 8-second segment extract of MITDB record 114 is displayed including the original signal and outputs after executing each of the processing steps BPF, ABS, and AVG.

Figure 3.

Signal extracts of executing the Hamilton algorithm over the MITDB record 114 (240 sec) a) Original ECG signal; b) Output after bandpass filtering; c) Output after differentiation and absolute calculation; d) Output after average over an 80 ms window.

Figure 3d identifies the static and dynamic thresholds and shows the case where the beats between the two high amplitude are considered as artifacts, although they should be real QRS beats.

The reason for bad detection of low-amplitude peaks after high-amplitude peaks is primarily due to the high level of static threshold. Even if one makes a correction by decreasing the static threshold value to include these peaks, there will still be a problem regardless of the fact that the peak will be treated as a candidate, and the dynamic threshold check will be applied. This is because the high-amplitude peak will increase the dynamic threshold value caused by the calculation of the mean, and so the peaks will be classified as noise peaks.

4.2 Isolated peaks after sequences of high-amplitude peaks

High-amplitude peaks directly affect the calculation of dynamic adaptive threshold. It is recalculated each time a new local maximum is found with an amplitude higher than the static threshold. In this case, it is considered as a potential peak, and the values lower than the dynamic threshold are considered as a noise peak, while the others as real QRS peaks. The original algorithm buffers the latest eight peak amplitudes and calculates a DAT value by Eq. (4).

Figure 4.

Signal extracts of executing the Hamilton algorithm over the MITDB record 201 (424 sec) a) Original ECG signal; b) Output after bandpass filtering; c) Output after differentiation and absolute calculation; d) Output after average over an 80 ms window.

An extensive analysis of the MITDB record 201 shows too many misses, especially in cases of aberrated atrial premature beats (classified as a beat) as illustrated in Fig. 4. Two of the beats highlighted by default cannot be captured due to the dynamic adaptive threshold and mean calculation, since most of the latest beats have a high amplitude. In this case, neither the static nor the dynamic adaptive threshold will work. The example is highlighted in Fig. 5 by peaks $C$ and $D$ , which should be classified as QRS peaks, but they are detected as artifacts because their value is lower than the static threshold.

Figure 5.

Static and dynamic threshold values on the output after average over an 80 ms window on MITDB record 201 (424 sec).

4.3 Classification of artifacts

Dynamic adaptive thresholding identifies noise and real peaks. Although in most cases, the dynamic adaptive threshold reacts properly, there are still cases where a noise peak is incorrectly calculated as real. Correct classification of artifacts is of primary importance for a quality industrial QRS detector.

Increasing the static threshold directly decreases the number of artifacts, however, this drastically increases the missed beats. This also applies to the dynamic adaptive threshold. To find an optimum of the static and dynamic peaks means to search for a comprise that would reach high values of both sensitivity, and positive predictive rate.

An example is illustrated in Fig. 5. Peaks labeled as $A$ , $C$ , $D$ , $E$ and $G$ are low amplitude peaks. However, peaks labelled as $B$ and $F$ are artifacts. With the default threshold, $A$ and $E$ peaks are considered as candidates for QRS, whereas the rest are considered as artifacts. For this particular case, decreasing the static threshold to 3 would catch all of the real peaks, but, will consider artifact peak $F$ as a peak. On the other hand, keeping the static threshold at 4 will only detect $G$ as a peak, and the other peaks will remain artifacts again.

4.4 Calculation of R-peak location

Figure 6.

Signal extracts and R-peak detection over the MITDB record 201 (426.4 sec): a) Original signal and local peaks $A$ , $B$ , $C$ , $D$ and $E$ ; b) Output after bandpass filter; c) Output after average over an 80 ms window; d) Real and calculated R-peak locations with constant delay introduced by the filter.

One of the issues in executing the original Hamilton algorithm is the proper detection of a QRS peak. Fig. 6 illustrates such a case, where local maxima are labeled with $A$ , $B$ , $C$ , $D$ , and $E$ . Another peak appears in the output after bandpass filter denoted by $F$ , as seen in Fig. 6b). The proximity of marked peaks $B$ , $C$ , and $F$ cause two local peaks on the output of 80 ms average window, marked as $B$ and $C$ on Fig. 6c. When static thresholding is applied to the time average over an 80 ms window, $C$ , $F$ , and $D$ peaks are detected to be artifacts, whereas $A$ , $B$ , and $E$ are identified as potential peaks. Since the dynamic threshold is below the static, these beats are classified as real.

Note that the filter makes a constant delay, which is deducted from the location of peak values (as displayed on the signal average output). It correctly determines properly the previous QRS peak $A$ , and the next QRS peak $E$ , but makes a mistake in determining the peak $B$ .

The original Hamilton’s algorithm detects the location of R-peaks to be at points $A^{C}$ , $B^{C}$ , and $E^{C}$ , even though they have different real peak locations $A^{R}$ , $B^{R}$ , and $E^{R}$ .

Hamilton’s algorithm detects a peak based on a calculated amount of delay [42]. Although filters produce a fixed delay, the author has pointed out that the detection delay can easily vary from 395 ms to 1 sec, depending on the heart rate and detection rule. It is important to note that the search back detection method [42]can produce a fixed delay if the search back algorithm does not report any local peaks. Thus, we can conclude that the original Hamilton’s algorithm provides the best likely position of the R-peak, however, it is not always exact. This problem can particularly increase the total number of FP’s. This particular case is observed in almost all MITDB records. Even though the difference is not so big, there are cases where the difference between the real and detected location is higher.

5. Algorithm improvement

Increasing the algorithm performance is directly related to fine-tuning thresholds and algorithm improvements.

5.1 Improving the detection of low-amplitude peaks

A careful analysis shows that the step (ABS) executed prior to time average, will not be able to cope with the bad performance in detection of low-amplitude peaks. This is especially crucial in cases when the signal is a mixture of one high amplitude beat, and then followed by several beats with low amplitudes. We used the idea introduced in the Pan Tompkins algorithm [4] to square the signal, instead of calculating the absolute value.

Figure 7 presents a case where a combination of a square mode and the optimized static threshold will improve the detection of low energy peaks. The peaks labelled as $A$ , $C$ , $E$ , $F$ , $H$ , and $I$ are real R peaks, whereas $B$ and $D$ are artifacts. The original algorithm, which uses the calculation of an absolute value and average along with the static threshold, is not able to classify $F$ , and $I$ as real peaks. However, the square mode and signal average in combination with a new (smaller) static threshold, is able to detect that there is sufficient energy for a potential peak. The square average signal also marks $B$ , and $D$ as potential peaks. Such peaks can be reduced with the dynamic threshold or by introducing rules for artifact detection.

Nevertheless, there are side effects, especially in the calculation of the dynamic adaptive threshold. This threshold increases due to increased amplitudes, and it becomes slightly difficult to adapt to sudden changes in the amplitudes.

Figure 7.

Signal extracts of executing the Hamilton’s algorithm over the MITDB record 201 (424 sec): a) Original signal; b) Output after average over a 80 ms window using the original Hamilton’s algorithm; c) Output after square average over a 80 ms window with the optimized static threshold.

This operation behaves as an important amplifier especially if it followed by a calculation of an average over an 80 ms moving window. It shows that this copes better with the identified problems. However, this is not enough, since this algorithm cannot perform with static threshold values, and needs dynamic calculation by other rules.

We have conducted several tests to experiment with threshold values STHR from 2 to 50, to find the optimal threshold value for which the number of errors is minimal. The left part of Fig. 13 presents the false detections for the conducted experiment. The performance of the algorithm gradually decreases as the threshold increases, mainly due to the high values of FP. The best performance is obtained for $\textit{STHR}=$ 2. We noticed a decreased number of FN, which gives an idea for how to get better performance if we decrease the number of FPs through other methods.

5.2 Improving the calculation of the R-peak location

Since the Hamilton algorithm only gives an approximation where the peak is, one way to reduce the number of FP’s that occur due to the determination of a proper R peak location, is to search the real peak in near proximity. The idea is to search the most convenient local maximum by analyzing the noise-eliminated output after bandpass filtering.

Figure 8.

Proper calculation of the R-peak location on the MITDB record 201 (426.4 sec): a) Calculated R-peak location $B^{C}$ and the search intervals over the output of bandpass filter; b) The re-calculated R-peak location $B^{RC}$ on the actual signal.

Figure 9.

False detections improving the calculation of the R-peak location.

The original Hamilton’s algorithm will determine the best approximate for the location of the R-peak. This is the starting point for the search of the local maximum, within a range, found SearchL ms to the left, and SearchR ms to the right of the approximated R-peak location. Once the local maximum is found on the output of the bandpass filter, we continue to find the local maximum on the actual signal, though the range for searching will be limited to 48 ms.

Figure 8 illustrates the basic steps of this improvement algorithm in the example presented in Fig. 6. The corresponding search segments are marked on the original signal. The original Hamilton’s algorithm detects the local peak $B^{C}$ , and our improvement algorithm finds $B^{RC}$ to be the real location.

We realized another experiment to locate the optimal values for SearchL, and SearchR. False detections for thresholding values are plotted on a surface graph presented in Fig. 9, for different values for SearchL, and SearchR, using the static threshold value $\textit{STHR}=$ 4. The best results are obtained when $\textit{SearchL}=$ 160 ms, and $\textit{SearchR}=$ 120 ms.

5.3 Improving the detection of low-amplitude peaks after sequences of high-amplitude peaks

Two functions, denoted as mean Eq. (3) and thresh Eq. (4), play a crucial role in the original QRS detection algorithm.

We proposed a change in the mean function in order to alleviate the effect of high amplitude complex proceeded by a significantly lower one. Instead of calculating the mean of the last 8 peaks, our algorithm considers only half of this value when the peak amplitude is higher than the dynamic threshold (DTHR) and the exact value in all other times, as defined by Eq. (5). This prevents a linear increase in the threshold especially for high amplitude signals and solves the identified problem.

$\displaystyle\textit{mean}=\frac{\sum_{n=1}^{8}\begin{cases}X_{n},&\textit{if % $X_{n}<\textit{DTHR}$}.\\ \frac{X_{n}}{2},&\textit{otherwise}.\end{cases}}{8}$ (5)

In addition, we changed the thresh method. Previously, the calculated threshold was multiplied by the constant $\textit{TH}=$ 0.3125. Since the square mode filter is used, we updated the multiplication constant by its square, i.e. $\textit{TH}*\textit{TH}$ . Thus the new DAT calculation is defined by Eq. (6). Both these interventions enabled the detection of such beats.

$\displaystyle\textit{DAT}=\textit{nmean}+(\textit{qmean-nmean})*\textit{TH}^{2}$ (6)

A good performance is achieved on both cases with:

•

a sequence of low-amplitude peaks after an isolated high-amplitude peak;

•

an isolated low-amplitude peak after a sequence of high-amplitude peaks.

Figure 10 illustrates the improvement idea for the example presented in Fig. 3, where a sequence of low-amplitude peaks is followed by a high-amplitude peak, which will increase the dynamic threshold to a value where all consequent low-amplitude peaks are marked as noise peaks.

The figure demonstrates the output of an average of an 80 ms window realized on squared (not absolute) values along with the static, and dynamic thresholds. Note that the detected potential peaks and their absolute values presented in Fig. 3 are smaller than the original dynamic threshold value.

Figure 10.

Effect of a dynamic threshold to detect peaks after average of squared values over an 80 ms window over the MITDB record 114 (240 sec).

Table 1

Parameter list for artifact detection rules

Parameter	Description
C	Current detected peak
P	Previous detected beat
CH	Current peak time average height
PH	Previous beat time average height
RR	Current beat to peak interval in ms
TA $x$	Time in ms optimizing A $x$ , $x\in\{0,1,2\}$
THRA $x$	Parameter optimizing A $x$ , $x\in\{1,2\}$
MS $x x x$	$x x x$ ms interval

Figure 11.

Effect of a dynamic threshold to detect peaks after average of squared values over an 80 ms window over the MITDB record 201 (424 sec).

The effect of applying the new way to calculate the dynamic threshold can be also observed in Fig. 11, presenting isolated, low-amplitude peaks after a sequence of high-amplitude peaks. The new minimum threshold detects 17 candidate peaks, whereas the original algorithm with default dynamic threshold is not able to capture the peaks labeled as $C$ , $D$ and $G$ . The newly optimized threshold is able to capture all 15 peaks correctly, and also classify $B$ , and $F$ peaks as noise.

We have conducted a lot of experiments to determine the optimized value for the $D T H R$ threshold value. The test cases included testing the threshold values of 100 up to 400. The middle part of Fig. 13 presents false detection as a function of the threshold values. Threshold values near 200 are the best candidates for the most effective performance.

5.4 Improvement of artifact elimination

In the final step, we introduce a new Classification phase. The aim of this is to classify whether the calculated beat is real or an artifact. Three important decision rules decide whether the peak is an artifact. If none of these rules are satisfied, the beat is calculated as a real peak.

From the preliminary analysis, we observed that artifacts generally follow a real beat and are closer than 320 ms away. The second important issue is that an artifact obviously has lower energy when compared to the previously detected beat. The original Hamilton’s algorithm eliminates artifacts closer than 195 ms. Our findings show that this value can also be optimized. Table 1 describes some parameters used in our optimization approaches. We introduce the following optimization rules for artifact elimination:

$\displaystyle C=\textit{Artifact}\ \text{if}$

A0:
$\textit{RR}\leqslant\textit{TA0}=\textit{MS250}$
A1:
$\textit{RR}\leqslant\textit{TA1}=\textit{MS260}\ \&\ \textit{PH/CH}>\textit{% THRA1}$
A2:
$\textit{RR}\leqslant\textit{TA2}=\textit{MS320}\ \&\ \textit{PH/CH}>=\textit{% THRA2}$

Figure 12.
False detections of optimization approaches for low-amplitude peaks (left); sequences of high-amplitude peaks (middle), beat rate impact (right).

Figure 13.
False detections of optimization approaches A0 (left), A1 (middle) and A2 (right).

Figure 14.
Signal extracts and outputs after squared average over 80 ms window of executing our algorithm over MITDB: a) signal and d) output for A0 type artifact in record 103 (1304.4 sec); b) signal and e) output for A1 type artifact in record 124 (413.3 sec); c) signal and f) output for A2 type artifact in record 101 (132.2 sec).

Examples of different types of detected artifacts are presented in Fig. 14. The identified segments demonstrate that the detected peaks are artifacts, since their energy is higher than static, and dynamic thresholds, but satisfies one of the A0, A1, and A2 rules. Otherwise, if none of these rules are satisfied, then the detected beat $C$ is not an artifact.

The results of the experiment using the A0 optimization approach, are presented in the left part of Fig. 13. A threshold value of $\textit{TA0}=\textit{MS250}=$ 250 ms yields promising results, due to the lowest level of false detections.

The impact of threshold values on the A1 approach is demonstrated in the middle part of Fig. 13. The $x$ -axis denotes values which are multiplied by 100 and the optimized value 80 for THRA1 corresponds to 80/100 $=$ 0.8.

The right part of Fig. 13 shows how the threshold parameter impacts the performance of the A2 optimization method. The $x$ -axis denotes values that are multiplied by 100, and the optimized value 250 corresponds to 250/100 $=$ 2.5 for THRA2.
5.5 Beat rate impact on artifacts

Our experiments have shown that the beat rate affects the intervals in the A0-A2 optimization approaches. Let $\textit{RR}_{\textit{avg}}$ be the average value of the last R to R intervals, and $f_{s}$ the sampling frequency. Then beat rate BPM is measured by beats per minute, and calculated by Eq. (7).

$\displaystyle\textit{BPM }=\frac{60}{(\textit{RR}_{\textit{avg}}/f_{s})}=\frac% {60f_{s}}{\textit{RR}_{\textit{avg}}}$ (7)

To classify whether a peak is an artifact, we have set three time constants $\textit{TA0}=$ 250 ms, $\textit{TA1}=$ 260 ms and $\textit{TA2}=$ 320 ms. Usually, the peaks closer than 250 ms are considered as peaks (corresponds to a beat rate higher of 240 BPM). The second and third thresholds correspondingly check if the peak is closer than 260 ms (corresponds to 230 BPM) or 320 ms (corresponds to 188 BPM).

However, our analysis has shown that premature beats might appear closer than these values. This is why we introduced a scaling factor to the previous improvement, and use the time thresholds calculated by Eq. (8), which are multiplied by the scaling factor BPM_BASE and heart rate BPM.

$\displaystyle\textit{RR}\leqslant\textit{TAx}\frac{\textit{BPM\_BASE}}{\textit% { BPM }}x\in\{0,1,2\}$ (8)

Table 2

Comparison of algorithm performance over MIT-BIH Arrhythmia database

MIT-BIH Arrhythmia database								All 48 records			No paced records (44)
Algorithm	$f_{S}$ (Hz)	Bits	Scale (mV)	Total peaks	TP	FP	FN	Tot err	$Q_{\textit{SE}}$	$Q_{+P}$	Tot err	$Q_{\textit{SE}}$	$Q_{+P}$
Our Work	125	10	6	109494	109382	110	112	222	99.90	99.90	194	99.91	99.90
Ghaffari [43]	360	11	10	109428	109327	129	101	230	99.91	99.88	$N/A$	$N/A$	$N/A$
Bahouraa[29]	250	11	10	116137	109625	133	174	307	99.83	99.88	303	99.82	99.88
Elgendi [44]	360	11	10	109985	109775	82	247	329	99.78	99.92	322	99.76	99.92
Martinez [40]	360	11	10	109428	109111	35	317	352	99.71	99.97	$N/A$	$N/A$	$N/A$
Martinez [31]	360	11	5	109428	109208	153	220	373	99.80	99.86	$N/A$	$N/A$	$N/A$
Cvikl [45]	250	$N/A$	$N/A$	109494	109294	200	200	400	99.82	99.82	373	99.81	99.82
Chiarugi [46]	360	11	10	109494	109228	210	266	476	99.76	99.81	443	99.75	99.81
Lee [47]	$N/A$	$N/A$	$N/A$	109481	109146	137	335	472	99.69	99.87	459	99.68	99.87
Zidemal [48]	360	$N/A$	$N/A$	109494	109101	193	393	586	99.64	99.82	540	99.64	99.83
Hamilton [5]	360	11	10	109267	108927	248	340	588	99.69	99.77	569	99.68	99.76
Choi [49]	360	11	$N/A$	109494	109118	218	376	594	99.66	99.80	561	99.65	99.79
GQRS [6]	360	11	10	109494	109196	302	298	600	99.73	99.72	562	99.72	99.72
Christov [50]	360	$N/A$	5	109855	109615	386	288	674	99.74	99.65	670	99.72	99.62
Arzeno [10]	360	11	10	$N/A$	109099	405	354	759	99.68	99.63	$N/A$	$N/A$	$N/A$
Pan Tompkins[4]	200	11	10	116137	109532	507	277	784	99.75	99.54	771	99.73	99.50
Paoletti [51]	360	11	10	109809	109430	565	379	944	99.65	99.49	924	99.64	99.45
Poli [33]	120	11	10	109963	109522	545	441	986	99.60	99.51	$N/A$	$N/A$	$N/A$
Elgendi [52]	360	11	10	109493	109397	97	1715	1812	98.31	99.92	1798	98.33	99.91
Chouakri [53]	360	11	10	109488	108043	3068	1446	4514	98.68	97.24	4155	98.63	97.29

The right part of Fig. 13 shows how the scaling factor BPM_BASE impacts the performance. The $x$ -axis shows the values of the scaling factor BPM_BASE in a range from 40 to 130, with increments of 2, and $y$ -axis the number of errors.

We observe that a value of $\textit{BPM\_BASE}=$ 90 minimizes the errors.

6. Evaluation and discussion

The following parameter combination achieves the best overall performance: $\textit{STHR}=$ 2, $\textit{SearchL}=$ 152 ms, $\textit{SearchR}=$ 56 ms, $\textit{DTHR}=$ 200, $\textit{TA0}=$ 250 ms, $\textit{TA1}=$ 260 ms, $\textit{TA2}=$ 320 ms, $\textit{THRA1}=$ 0.8, $\textit{THRA2}=$ 2.5, and $\textit{BPM\_BASE}=$ 90.

Table 2 gives an overview of the obtained results and shows that our algorithm has reached better-combined values of sensitivity and positive predictive rate.

One can face several problems comparing any QRS detection method to other published papers as summarized below:

•
no source code provided to check other approaches;
•
no info about positive predictive rate; or
•
no info about the number of errors.

When analyzing the performance, only a small number of papers give information on the achieved positive predictive rate and they usually target achieving higher sensitivity values. However, it is very easy to achieve a higher sensitivity value and capture most of the results you would like to include in your algorithm by relaxing the constraints on the optimization parameters, but, at the same time, this will produce many extra generated peaks that do not represent a QRS peak. This is why it is very important to address both the sensitivity and positive predictive to evaluate performance.

This means that there is no direct comparison method. To cope with this problem, we have analyzed the number of errors as a performance measure (as defined by Errors in Eq. (2)). Although this performance measure can be achieved by calculation the sum of FP and FN, it can also be calculated through the harmonic mean (HM) of the QRS sensitivity and positive predictive rate by Eq. (9).

$\displaystyle\textit{Errors}=\textit{TotalQRS }\Bigg{(}\frac{1}{Q_{\textit{SE% }}}+\frac{1}{Q_{+P}}-2\Bigg{)}$ (9)

We used Eq. (10) to evaluate this relation. In addition, we assumed that the number of extra detected (false negative) peaks is a lot smaller than the number of correctly detected QRS peaks, i.e. $\textit{TP}>>\textit{FN}$ , leading to $\textit{TotalQRS}\approx\textit{TP}$ .

$\displaystyle\frac{1}{Q_{\textit{SE}}}+\frac{1}{Q_{+P}}=\frac{\textit{TP}+% \textit{FN}}{\textit{TP}}+\frac{\textit{TP}+\textit{FP}}{\textit{TP}}=2+\frac{% \textit{FP}+\textit{FN}}{\textit{TP}}$ (10)

A number of discrepancies were found while analyzing related work. Table 1 of [4] shows that the number of errors is 782 although it is 784. The calculation of total beats is the most unambiguous. For example, TP $+$ FN is greater than TB in [52]. Lee publishes two papers [54, 47] providing correspondingly 109486 and 109481 total beats. The former one has 6 additional beats, which are from files 118, 201, 205, 220, 221* and 233, whereas the latter one has 1 additional beat at record 114.

We found that various authors have used a different number of detected peaks. They should use the total number of detected beats, since the total number of peaks includes also non-beat annotations, for example, locations where there is a rhythm change. That is why we use the total number of annotated beats in the MIT-BIH Arrhythmia database 109494.

Table 2 also shows that researchers focusing on QRS detection algorithms mostly tend to use the original sampling frequencies, and resolution of reference ECG databases. Running the algorithm on 125 Hz means that nearly 3 times less data is feed to the QRS detector, thus the execution times are deduced accordingly. Moreover, the adjustments that we proposed increased the performance metrics, which in total yields a better QRS detector intended for a small wearable one channel ECG sensor.
7. Conclusion

We have introduced several methods to improve the QRS detection in differentiation-based algorithms. Even though our approach is demonstrated on Hamilton’s QRS detection algorithm, it can also be implemented in other algorithms. The results show superior performance over other published results. The improvements were efficiently built in on an industrial QRS detector, for a wearable ECG monitor, where the sampling frequency is 125 Hz, with 10-bit resolution of the AD converter.

Tuning the threshold values might increase the performance, but one needs to develop new relations to generate the thresholds, such as finding

•

how many previous beats might be analyzed to estimate the mean value,

•

which peaks will be classified as QRS beats since they are very similar in shape, and

•

the impact of the beat rate on classifying the artifacts.

For this purpose, we have introduced several new rules for 1) threshold calculation to better classify QRS peaks; 2) elimination of QRS-like artifacts, and 3) artifact elimination based on beat rate.

Due to rescaling, loose ECG contacts, and noise caused by muscles, 62 beats cannot be detected by analyzing the first ECG channel, without the analysis of a second channel. This gives a higher performance of the QRS sensitivity to 99.96%, and our algorithm has reached a 99.90% QRS sensitivity, with 99.90% positive predictive rate when all records are analyzed in the MITDB, and 99.91% QRS sensitivity for 44 records without paced beats.

Generally, the published papers on QRS detection algorithms do not offer their source codes, and only some of them are validated with referenced ECG databases, such as the MIT-BIH Arrhythmia database. Most of the algorithms give a brief description of design issues without implementation details, addressing only theoretical related results. This is why one cannot directly compare results. Different approaches generally do not achieve results as the ones achieved in a real implementation.

This study can facilitate possible QRS detection algorithms to consider rescaling, and resampling on reference ECG databases in return for better performances. Our findings show that the adjusted version of Hamilton’s QRS detection algorithm yields better results with 125 Hz data on nearly three times shorter execution times.

Our future work will be in the direction of classifying QRS detection errors, and eventually eliminating them, as well as for the purpose of modeling the dependence on sampling rate, and bit resolution. Additionally, we are considering building a quality industrial heart beats classifier, with higher performance, on top of the detection algorithm.

Footnotes

Conflict of interest

None to report.

References

Böhm

Brüser

Leonhardt

. A novel BCG sensor-array for unobtrusive cardiac monitoring. Acta Polytechnica. 2013; 53(6).

Álvarez

Penín

AJM

Sobrino

XAV

. A comparison of three QRS detection algorithms over a public database. Procedia Technology. 2013; 9: 1159-1165.

Moody

Mark

. The impact of the MIT-BIH arrhythmia database. IEEE Engineering in Medicine and Biology Magazine. 2001; 20(3): 45-50.

Pan

Tompkins

. A real-time QRS detection algorithm. IEEE transactions on biomedical engineering. 1985; (3): 230-236.

Hamilton

Tompkins

. Quantitative investigation of QRS detection rules using the MIT-BIH arrhythmia database. IEEE transactions on biomedical engineering. 1986; 12): 1157-1165.

Goldberger

Amaral

Glass

Hausdorff

Ivanov

Mark

, et al. Physiobank, physiotoolkit, and physionet. Circulation. 2000; 101(23): e215-e220.

Pahlm

Sörnmo

. Software QRS detection in ambulatory monitoring: A review. Medical and Biological Engineering and Computing. 1984; 22(4): 289-297.

Friesen

Jannett

Jadallah

Yates

Quint

Nagle

. A comparison of the noise sensitivity of nine QRS detection algorithms. IEEE Transactions on biomedical engineering. 1990; 37(1): 85-98.

Kohler

Hennig

Orglmeister

. The principles of software QRS detection. IEEE Engineering in Medicine and Biology Magazine. 2002; 21(1): 42-57.

10.

Arzeno

Deng

Poon

. Analysis of first-derivative based QRS detection algorithms. IEEE Transactions on Biomedical Engineering. 2008; 55(2): 478-484.

11.

Ahlstrom

Tompkins

. Automated high-speed analysis of Holter tapes with microcomputers. IEEE Transactions on Biomedical Engineering. 1983; (10): 651-657.

12.

Fraden

Neuman

. QRS wave detection. Medical and Biological Engineering and computing. 1980; 18(2): 125-132.

13.

Morizet-Mahoudeaux

Moreau

Quarante

. Simple microprocessor-based system for on-line ECG arrhythmia analysis. Medical and Biological Engineering and Computing. 1981; 19(4): 497-500.

14.

Afonso

Tompkins

Nguyen

Luo

. ECG beat detection using filter banks. IEEE transactions on biomedical engineering. 1999; 46(2): 192-202.

15.

Borjesson

Pahlm

Sornmo

Nygards

. Adaptive QRS detection based on maximum a posteriori estimation. IEEE Transactions on Biomedical Engineering. 1982; (5): 341-351.

16.

Engelse

Zeelenberg

. A single scan algorithm for QRS-detection and feature extraction. Computers in Cardiology. 1979; 6(1979): 37-42.

17.

Fancott

Wong

. A minicomputer system for direct high speed analysis of cardiac arrhythmia in 24 h ambulatory ECG tape recordings. IEEE Transactions on Biomedical Engineering. 1980; (12): 685-693.

18.

Okada

. A digital filter for the QRS complex detection. IEEE Transactions on Biomedical Engineering. 1979; (12): 700-703.

19.

Gusev

Ristovski

Guseva

. Pattern recognition of a digital ECG. International Conference on ICT Innovations. Springer; 2016; p. 93-102.

20.

Ciaccio

Dunn

Akay

. Biosignal pattern recognition and interpretation systems. IEEE Engineering in Medicine and Biology Magazine. 1993; 12(3): 89-95.

21.

Maglaveras

Stamkopoulos

Diamantaras

Pappas

Strintzis

. ECG pattern recognition and classification using non-linear transformations and neural networks: a review. International journal of medical informatics. 1998; 52(1-3): 191-208.

22.

Steinberg

Abraham

Caceres

. Pattern recognition in the clinical electrocardiogram. IRE Transactions on Bio-Medical Electronics. 1962; 9(1): 23-30.

23.

Trahanias

Skordalakis

. Syntactic pattern recognition of the ECG. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1990; 12(7): 648-657.

24.

Xue

Tompkins

. Neural-network-based adaptive matched filtering for QRS detection. IEEE Transactions on Biomedical Engineering. 1992; 39(4): 317-329.

25.

Barro

Fernandez-Delgado

Vila-Sobrino

Regueiro

Sanchez

. Classifying multichannel ECG patterns with an adaptive neural network. IEEE Engineering in Medicine and Biology Magazine. 1998; 17(1): 45-55.

26.

Dokur

Ölmez

Yazgan

Ersoy

. Detection of ECG waveforms by neural networks. Medical Engineering and physics. 1997; 19(8): 738-741.

27.

Tompkins

Urrusti

Afonso

. Applications of artificial neural networks for ECG signal detection and classification. Journal of electrocardiology. 1993; 26: 66-73.

28.

Zheng

Tai

. Detection of ECG characteristic points using wavelet transforms. IEEE Transactions on biomedical Engineering. 1995; 42(1): 21-28.

29.

Bahoura

Hassani

Hubin

. DSP implementation of wavelet transform for real time ECG wave forms detection and heart rate analysis. Computer methods and programs in biomedicine. 1997; 52(1): 35-44.

30.

Shambi

Tandon

Bhatt

. Using wavelet transforms for ECG characterization. IEEE Engineering in Medicine and Biology. 1997; p. 77-83.

31.

Martínez

Almeida

Olmos

Rocha

Laguna

. A wavelet-based ECG delineator: evaluation on standard databases. IEEE transactions on biomedical engineering. 2004; 51(4): 570-581.

32.

Milchevski

Gusev

. Improved pipelined wavelet implementation for filtering ECG signals. Pattern Recognition Letters. 2017; 95: 85-90.

33.

Poli

Cagnoni

Valli

. Genetic design of optimum linear and nonlinear QRS detectors. IEEE Transactions on Biomedical Engineering. 1995; 42(11): 1137-1141.

34.

Coast

Stern

Cano

Briller

. An approach to cardiac arrhythmia analysis using hidden Markov models. IEEE Transactions on biomedical Engineering. 1990; 37(9): 826-836.

35.

Andreão

Dorizzi

Boudy

. ECG signal analysis through hidden Markov models. IEEE Transactions on Biomedical engineering. 2006; 53(8): 1541-1549.

36.

Cost

Cano

. QRS detection based on hidden Markov modeling. In: Engineering in Medicine and Biology Society, 1989. Images of the Twenty-First Century, Proceedings of the Annual International Conference of the IEEE Engineering in. IEEE; 1989. p. 34-35.

37.

Song-Kai

Jian-Tao

Jun-Rong

. The real-time detection of QRS-complex using the envelope of ECG. In: Engineering in Medicine and Biology Society, 1988. Proceedings of the Annual International Conference of the IEEE. IEEE; 1988. p. 38s.

38.

Nygårds

Sörnmo

. Delineation of the QRS complex using the envelope of the ECG. Medical and Biological Engineering and Computing. 1983; 21(5): 538-547.

39.

Benitez

Gaydecki

Zaidi

Fitzpatrick

. The use of the Hilbert transform in ECG signal analysis. Computers in biology and medicine. 2001; 31(5): 399-406.

40.

Martínez

Alcaraz

Rieta

. Application of the phasor transform for automatic delineation of single-lead ECG fiducial points. Physiological measurement. 2010; 31(11): 1467.

41.

http://www.heartorg.AmericanHealthAssociation,AHAdatabase.

42.

Hamilton

. Open Source ECG Analysis Software Documentation; 2002.

43.

Ghaffari

Homaeinezhad

Akraminia

Atarod

Daevaeiha

. A robust wavelet-based multi-lead electrocardiogram delineation algorithm. Medical Engineering and Physics. 2009; 31(10): 1219-1227.

44.

Elgendi

Mohamed

Ward

. Efficient ECG Compression and QRS Detection for E-Health Applications. Scientific Reports. 2017; 7(1): 459.

45.

Cvikl

Jager

Zemva

. Hardware implementation of a modified delay-coordinate mapping-based QRS complex detection algorithm. EURASIP Journal on Applied Signal Processing. 2007; 2007(1): 104-104.

46.

Chiarugi

Sakkalis

Emmanouilidou

Krontiris

Varanini

Tollis

. Adaptive threshold QRS detector with best channel selection based on a noise rating system. In: Computers in Cardiology, IEEE; 2007, p. 157-160.

47.

Lee

Jeong

Yoon

Lee

. A simple real-time QRS detection algorithm. In: Engineering in Medicine and Biology Society, 1996. Bridging Disciplines for Biomedicine. Proceedings of the 18th Annual International Conference of the IEEE. vol. 4. IEEE; 1996. p. 1396-1398.

48.

Zidelmal

Amirou

Adnane

Belouchrani

. QRS detection based on wavelet coefficients. Computer Methods and Programs in Biomedicine. 2012; 107(3): 490-496.

49.

Choi

Adnane

Lee

Jang

Jiang

Park

. Development of ECG beat segmentation method by combining lowpass filter and irregular R-R interval checkup strategy. Expert Systems with Applications. 2010; 37(7): 5208-5218.

50.

Christov

. Real time electrocardiogram QRS detection using combined adaptive threshold. Biomedical engineering online. 2004; 3(1): 28.

51.

Paoletti

Marchesi

. Discovering dangerous patterns in long-term ambulatory ECG recordings using a fast QRS detection algorithm and explorative data analysis. Computer Methods and programs in biomedicine. 2006; 82(1): 20-30.

52.

Elgendi

Jonkman

DeBoer

. Frequency bands effects on QRS detection. Pan. 2010; 5: 15 Hz.

53.

Chouakri

Bereksi-Reguig

Taleb-Ahmed

. QRS complex detection based on multi wavelet packet decomposition. Applied Mathematics and Computation. 2011; 217(23): 9508-9525.

54.

Lee

Kim

Lee

. A real time QRS detection using delay-coordinate mapping for the microcontroller implementation. Annals of biomedical Engineering. 2002; 30(9): 1140-1151.