QRS Detection Based on Medical Knowledge and Cascades of Moving Average Filters

: Heartbeat detection is the ﬁrst step in automatic analysis of the electrocardiogram (ECG). For mobile and wearable devices, the detection process should be both accurate and computationally efﬁcient. In this paper, we present a QRS detection algorithm based on moving average ﬁlters, which affords a simple yet robust signal processing technique. The decision logic considers the rhythmic and morphological features of the QRS complex. QRS enhancing is performed with channel-speciﬁc moving average cascades selected from a pool of derivative systems we designed. We measured the effectiveness of our algorithm on well-known benchmark databases, reporting F1 scores, sensitivity on abnormal beats and processing time. We also evaluated the performances of other available detectors for a direct comparison with the same criteria. The algorithm we propose achieved satisfying performances on par with or higher than the other QRS detectors. Despite the performances we report are not the highest that have been published so far, our approach to QRS detection enhances computational efﬁciency while maintaining high accuracy. In the algorithm we present reduction of noise and undesired waves is performed by applying to the signal two different moving average cascades (MACs) of different length. The synchronized outputs of the two MACs are then subtracted. The ﬁrst MAC consists of the cascade of a 20 ms moving average with a 8 ms one. Thus, we refer to it as MAC(0.02, 0.008). At a sampling frequency of 250 Hz, the impulsive response of MAC(0.02, 0.008) ms long. This length reduces narrow peaks of noise, slightly impacting QRS peaks. The second MAC is MAC(0.12, 0.024). For a sampling frequency of 250 Hz, this ﬁlter has a 140 ms long impulsive response interval. The output of MAC(0.12, 0.024) tracks baseline wandering partially, T Thus, their amplitude at output of the synchronized difference of the two MACs is reduced. resulting impulsive response of noise reduction


Introduction
Electrocardiography is the recording of the electrical activity originated by the heart in its duty cycle. The analysis of both rhythm and shape of the characteristic PQRST pattern can reveal altered conditions of the cardiac system. Software-based analysis of the electrocardiogram is a valuable tool, even more so with the increased presence of portable and wearable systems. The first and foremost step of ECG software analysis is heartbeat detection, specifically the ventricular systole which is represented by the combination of the Q, R, and S characteristic waves, the QRS complex. In the ECG analysis pipeline the accuracy of heartbeat detection determines the result quality of the following stages. Thus, this first step should be as accurate as possible. The main limiting factors in QRS detection accuracy are undesired noise, prominent T waves, and abnormal (non-sinus) beats.
Scientific literature on heartbeat detection is rich, and several approaches have been proposed over time. The most notable and referenced is Pan and Tompkins' algorithm [1], to which newer approaches are still compared to. Their work is based on differentiation, rectification and integration of the filtered ECG signal and a decision threshold calculated from the estimates of signal and noise levels. Several techniques have been applied to QRS detection [2], like filter banks, neural networks, wavelets, and other means of signal decomposition. Recently, Phukpattaranont [3] introduced quadratic filtering to QRS detection that attempts a thorough frequency characterization of the ECG. Wavelet-based processing has been a mainstay of ECG processing and heartbeat detection [4][5][6][7]. Pal et al. [8] and Hossain et al. [9] proposed methods based on empirical mode decomposition while Bashar et al. [10] proposed an algorithm consisting in a set of rules applied on variable frequency complex demodulation, a time-frequency representation. Another time-frequency representation was applied in QRS detection by Zidelmal et al. [11], the S-Transform. to our knowledge. We did not edit the code of the detection routines of other authors in any way.
Authors of QRS detection literature often report performances on a single database. This is the case for many of the works mentioned above. While the printed results may be excellent, the generalizability of the proposed methods to real world data remains indeterminate. For this reason, we included several databases in the evaluation of our work. When optimizing an algorithm on several data sources, overperforming on a single data source may considerably reduce the accuracy achieved on the other databases.
In this paper, we propose a fast, accurate approach to QRS detection based on combining simple unweighted moving averages. In particular, the basic elements of our ECG signal processing are cascades of moving averages (MACs). Our algorithm enhances the QRS complex by selecting a MAC from a pool of wide, derivative MACs that are inherently resistant to noise. We compare the performance of the QRS detector with those of three public algorithms, both in terms of accuracy and execution time, on the same calculator and the same initial conditions. The algorithm we propose aims to be an efficient, accurate, and generalizable method for QRS detection that improves on the concept of derivative transforms.
Our software was developed with MATLAB (The Mathworks Inc., Natick, MA, USA). A preliminary version of this work has been reported in [27]. The algorithm has since then been further optimized, and one of the algorithms we considered for comparison has been changed to Kim's detection approach, as it performed better. In this paper, we also elaborate with greater detail on the decision logic. Additionally, we report results on additional databases: the INCART, T Wave Alternans, and MIT-BIH Noise Stress Test databases.

Data
The algorithms we describe were tested on five public datasets available on the Physionet repository [28]: the MIT-BIH Arrhythmia (here abbreviated as MIT-BIH) database, the European ST-T (Eu ST-T) database, INCART database, the MIT-BIH Noise Stress Test (NST) database, the T Wave Alternans (TWA) database and the QT database. Table 1 is a summary of their features. These databases are often the common ground of performance validation for studies on detection of ECG characteristic waves. Considering multiple databases entails the advantage of testing the algorithms on a wide range of possible scenarios influenced by noise and physiological variability. The MIT-BIH Arrhythmia database [24] is a collection of 48 records sampled at 360 Hz. Each record is 30 min long and comprises two channels. The first channel is a modified lead II. The second channel is often akin to lead V 1 and occasionally to V 2 , V 4 , or V 5 . We report the performances for QRS detection without excluding any record.
The European ST-T Database [29] is a collection of ambulatory ECGs developed for the study of ST tract and T wave changes in myocardial ischemia and other abnormalities. The collection consists of 90, 2-h long excerpts, sampled at 250 Hz. The first channel is usually a modified limb lead (I or III), where the QRS complex is evident. The second channel is usually a precordial lead, where the ST episodes are prominent (V 1-5 ).
The St. Petersburg INCART 12-lead Arrhythmia Database [28] includes 75 excerpts sampled at 257 Hz, 30 min long each. The excerpts were extracted from 32 Holter recordings of different patients testing for coronary artery disease. The included annotations were first produced by an algorithm and then corrected manually. Each recording in the database is a full 12-lead recording.
The MIT-BIH Noise Stress test database [30] combines two clean recordings of the MIT-BIH Arrhythmia database with excerpts of noise. Noise sources include baseline wandering, miographic noise and electrode motion artifacts. The resulting collection consists of 12 half-hour ECG records, sampled at 360 Hz, with two channels each.
The T Wave Alternans database [31] is composed by 100 records, each with a number of channels that can vary from 2 to 12, sampled at 500 Hz. The recordings have been obtained from both healthy patients and patients with risk factors for sudden cardiac death. Additionally, some records are synthesized to display calibrated amounts of T wave alternans. Each record is about 2 min long.
The QT Database [32] is an anthological database composed by 105 ECG excerpts, extracted from seven other public databases. However, only 82 records are accompanied by the '.atr' files containing the beat annotations. Each record is 15 min long, is resampled at 250 Hz and includes two channels. A total of 15 records of the QT Database are derived from the MIT-BIH Arrhythmia Database, while 33 of them are extracted from the Eu ST-T database. In our study, the annotations for the sel232 record of the QT Database were replaced by the corresponding ones from record 232 of the MIT-BIH Arrhythmia Database.

Methods
ECG signal analysis has often been approached with methods and mindsets pertaining to the field of telecommunications-i.e., frequency characterized signals and linear systems. As a consequence, signal analysis experts who approached ECG signal analysis were trained in that area and therefore applied methods and techniques of the frequency domain to the ECG signal.
It should be noted that the information of physiological and pathophysiological interest in the ECG signal is concentrated in localized and separable events in time. The P wave is associated with atrial depolarization, the QRS complex with ventricular depolarization and the T wave with ventricular repolarization. On the contrary, in the frequency domain, there is significant overlap of this information and it is not separable. For this reason substantial frequency-based ECG signal analysis is unjustified. There are, however, a few instances of specific frequency components in the ECG like powerline disturbance at 50 or 60 Hz and atrial or ventricular fibrillation. Therefore, we designed the signal preprocessing and QRS detection by taking into account the time domain properties of the QRS complex with respect to other events. In addition, to achieve a small processing time, we used combinations of moving average cascades (MACs). A MAC is composed by applying two or more unweighted moving averages in series. When two moving averages of the same length are combined in cascade, the resulting system will exhibit a triangular impulsive response. If two moving averages of different length are combined, the resulting impulsive response will take a trapezoidal shape. As the number of moving averages in cascade increases, the shape of the resulting impulsive response will get smoother, wider, and increasingly similar to the Gaussian bell function. This effect is shown in Figure 1. These cascades of moving averages are equivalent to linear FIR filters. However, we remark that they are designed in time, by taking advantage of the knowledge of the pattern of the target ECG event. Moreover, their application does not require any convolution, which is a time-consuming operation.
We built two types of linear systems by combining MACs. The first linear system consists in the difference of two centered MACs of different length (Section 2.2.1), which was used to enhance events characterized by mono-directional changes whose time duration is comprised between the lengths of the two MACs. The second linear system (Section 2.2.2) consists in the difference of two time-shifted MACs and was designed to enhance events We built two types of linear systems by combining MACs. The first linear syst consists in the difference of two centered MACs of different length (Section 2.2.1), wh was used to enhance events characterized by mono-directional changes whose time du tion is comprised between the lengths of the two MACs. The second linear system (Sect 2.2.2) consists in the difference of two time-shifted MACs and was designed to enhan events whose ascending or descending fronts' duration is greater than the shift of the t MACs. This second system enhances QRS complexes and has a derivative behavior.
The proposed QRS detection algorithm is articulated in three steps: noise reduction i.e., non QRS signal changes; enhancing QRS events with a derivative transform; detect heartbeats. The absolute value of the derivative transform will also be referred to as feature signal. Figure 2 shows a diagram of the algorithm. This algorithm does n resample the signal. The filters automatically scale at the sampling frequency of the EC record that is being processed.  The proposed QRS detection algorithm is articulated in three steps: noise reductioni.e., non QRS signal changes; enhancing QRS events with a derivative transform; detecting heartbeats. The absolute value of the derivative transform will also be referred to as the feature signal. Figure 2 shows a diagram of the algorithm. This algorithm does not resample the signal. The filters automatically scale at the sampling frequency of the ECG record that is being processed. We built two types of linear systems by combining MACs. The first linear system consists in the difference of two centered MACs of different length (Section 2.2.1), which was used to enhance events characterized by mono-directional changes whose time duration is comprised between the lengths of the two MACs. The second linear system (Section 2.2.2) consists in the difference of two time-shifted MACs and was designed to enhance events whose ascending or descending fronts' duration is greater than the shift of the two MACs. This second system enhances QRS complexes and has a derivative behavior.
The proposed QRS detection algorithm is articulated in three steps: noise reductioni.e., non QRS signal changes; enhancing QRS events with a derivative transform; detecting heartbeats. The absolute value of the derivative transform will also be referred to as the feature signal. Figure 2 shows a diagram of the algorithm. This algorithm does not resample the signal. The filters automatically scale at the sampling frequency of the ECG record that is being processed.

Noise Reduction
The target of this step is to reduce the effects of two non-QRS components on the signal. The first one is undesired noise, which is a ubiquitous aspect of biosignals. In ECG, noise can be observed as baseline wander, power line interference, spike noise, and myographic noise. The other undesired aspect in heartbeat detection are prominent T waves, as they could be mistaken for a QRS complex. Thus, the first step of the algorithm reduces the amplitude of these two sources of error.
In the algorithm we present reduction of noise and undesired waves is performed by applying to the signal two different moving average cascades (MACs) of different length. The synchronized outputs of the two MACs are then subtracted. The first MAC consists of the cascade of a 20 ms moving average with a 8 ms one. Thus, we refer to it as MAC(0.02, 0.008). At a sampling frequency of 250 Hz, the impulsive response of MAC(0.02, 0.008) is 24 ms long. This length reduces narrow peaks of noise, while only slightly impacting QRS peaks. The second MAC is MAC(0.12, 0.024). For a sampling frequency of 250 Hz, this filter has a 140 ms long impulsive response interval. The output of MAC(0.12, 0.024) tracks baseline wandering and, partially, T waves. Thus, their amplitude at the output of the synchronized difference of the two MACs is reduced. Figure 3 shows the resulting impulsive response of the noise reduction filter.
noise can be observed as baseline wander, power line interference, spike noise, and myo graphic noise. The other undesired aspect in heartbeat detection are prominent T wave as they could be mistaken for a QRS complex. Thus, the first step of the algorithm reduce the amplitude of these two sources of error.
In the algorithm we present reduction of noise and undesired waves is performed b applying to the signal two different moving average cascades (MACs) of different length The synchronized outputs of the two MACs are then subtracted. The first MAC consis of the cascade of a 20 ms moving average with a 8 ms one. Thus, we refer to it as MAC(0.0 0.008). At a sampling frequency of 250 Hz, the impulsive response of MAC(0.02, 0.008) 24 ms long. This length reduces narrow peaks of noise, while only slightly impacting QR peaks. The second MAC is MAC(0.12, 0.024). For a sampling frequency of 250 Hz, th filter has a 140 ms long impulsive response interval. The output of MAC(0.12, 0.024) track baseline wandering and, partially, T waves. Thus, their amplitude at the output of th synchronized difference of the two MACs is reduced. Figure 3 shows the resulting impu sive response of the noise reduction filter. The design of this QRS detector is based on examination of the shape, time of occu rence and duration of the waves that represent ventricular depolarization. These wave are distinguished from artifacts and the P, T characteristic waves that concern othe events. The QRS complex is typically 70-110 ms long, with two to four rising and fallin fronts. Considering the largest amplitude wave (either Q R or S), its duration can var from 30-70 ms. Therefore, a MAC(0.02, 0.008) system that performs a moving averag over approximately 24 ms reduces the fast variations caused by noise without consistentl reducing the QRS peaks. The MAC(0.12, 0.024) system is designed to follow the slow base line variations but not the fast QRS variations. In fact, this system averages over 120 ms o signal and does not track the variations of the QRS waves.

Derivative Filter Selection
The next step of the algorithm further enhances the QRS complex with respect to th other waves. The QRS complex features a high derivative which involves an interval o about 100 ms. Thus, a derivative transform enhances the prominent characteristic of th The design of this QRS detector is based on examination of the shape, time of occurrence and duration of the waves that represent ventricular depolarization. These waves are distinguished from artifacts and the P, T characteristic waves that concern other events. The QRS complex is typically 70-110 ms long, with two to four rising and falling fronts. Considering the largest amplitude wave (either Q R or S), its duration can vary from 30-70 ms. Therefore, a MAC(0.02, 0.008) system that performs a moving average over approximately 24 ms reduces the fast variations caused by noise without consistently reducing the QRS peaks. The MAC(0.12, 0.024) system is designed to follow the slow baseline variations but not the fast QRS variations. In fact, this system averages over 120 ms of signal and does not track the variations of the QRS waves.

Derivative Filter Selection
The next step of the algorithm further enhances the QRS complex with respect to the other waves. The QRS complex features a high derivative which involves an interval of about 100 ms. Thus, a derivative transform enhances the prominent characteristic of the QRS complex. However, the degree of enhancement will depend on how the derivative matches the QRS shape. Hence, a specific derivative transform may better enhance the QRS complex in some ECG recordings than others. For this reason, the authors designed a dictionary of 18 different derivative filters. Each filter results from the subtraction between two MACs, involving about 8-30 ms of the signal each and shifted by 8-30 ms. These derivative filters will be referred to as dMACs. A few examples of such derivative filters (dMACs) are shown in Figure 4. filters (dMACs) are shown in Figure 4.
These derivative filters are designed to match the rising or falling fronts of the Q complex. Thus, the output of these filter is highest when the length of the longest (ris or falling) front of the QRS complex matches the main front of the derivative filter's i pulsive response. This characteristic makes these filters inherently specific for the tar signal, and less susceptible to noise. Since the QRS complexes do not have a perfectly fix duration and shape in each ECG recording, different shapes of derivative filters may p form better than others in different recordings. The dMAC that better enhances the QRS complexes of a specific record is autom cally chosen for each channel in the initialization phase of the algorithm. This choice performed by maximization of the signal quality index (SQI) defined as where mD is the trimmed mean of the derivative signal maxima computed in success windows of 1.8 s, aimed at representing the QRS-specific high derivative values. In 1 at least one heartbeat is expected. The trimmed mean discards undesired outliers. T term mD is the trimmed mean of the derivative signal maxima computed in success windows of 0.1 s that accounts for the background noise. Most of the 0.1 s long windo of signal will not include a heartbeat. The trimmed mean discards the few high values t relate to the windows containing a heartbeat. When the input of the SQI is an ECG o derivative of an ECG, the term mD estimating the QRS derivative is high, resulting i high signal quality index. When the input is noise, mD and mD have a similar va and the SQI is approximately unitary. The SQI is calculated in a selected 15 s interval each of the 18 derivative filters we designed. The derivative filter that featured the high SQI value is chosen as the derivative filter for the remainder of the record. The SQI previously described was defined from basic fuzzy considerations on EC Thus, the values of all the parameters involved in SQI definition are not critical. T These derivative filters are designed to match the rising or falling fronts of the QRS complex. Thus, the output of these filter is highest when the length of the longest (rising or falling) front of the QRS complex matches the main front of the derivative filter's impulsive response. This characteristic makes these filters inherently specific for the target signal, and less susceptible to noise. Since the QRS complexes do not have a perfectly fixed duration and shape in each ECG recording, different shapes of derivative filters may perform better than others in different recordings.
The dMAC that better enhances the QRS complexes of a specific record is automatically chosen for each channel in the initialization phase of the algorithm. This choice is performed by maximization of the signal quality index (SQI) defined as where mD s is the trimmed mean of the derivative signal maxima computed in successive windows of 1.8 s, aimed at representing the QRS-specific high derivative values. In 1.8 s at least one heartbeat is expected. The trimmed mean discards undesired outliers. The term mD n is the trimmed mean of the derivative signal maxima computed in successive windows of 0.1 s that accounts for the background noise. Most of the 0.1 s long windows of signal will not include a heartbeat. The trimmed mean discards the few high values that relate to the windows containing a heartbeat. When the input of the SQI is an ECG or a derivative of an ECG, the term mD s estimating the QRS derivative is high, resulting in a high signal quality index. When the input is noise, mD s and mD n have a similar value and the SQI is approximately unitary. The SQI is calculated in a selected 15 s interval for each of the 18 derivative filters we designed. The derivative filter that featured the highest SQI value is chosen as the derivative filter for the remainder of the record. The SQI previously described was defined from basic fuzzy considerations on ECG. Thus, the values of all the parameters involved in SQI definition are not critical. The window length for mD s can range between 1 and 2 s. The window length for mD n can range between 0.05 and 0.3 s. The constant k s can be omitted, while k n has a very small value and should be included to avoid occasional divisions by zero. Figure 5 shows an example of the derivative filters selected with SQI maximization for different records of the MIT-BIH Arrhythmia database. Each row refers to the first channel of a different record: top row refers to record 100, middle row refers to record 200 and the bottom row refers to record 203. Narrow derivative filters are selected to enhance narrow QRS complexes (top row), while wide derivative filters are selected for records that display wider QRS complexes (middle and bottom row). window length for mD can range between 1 and 2 s. The window length for mD can range between 0.05 and 0.3 s. The constant k can be omitted, while k has a very small value and should be included to avoid occasional divisions by zero. Figure 5 shows an example of the derivative filters selected with SQI maximization for different records of the MIT-BIH Arrhythmia database. Each row refers to the first channel of a different record: top row refers to record 100, middle row refers to record 200 and the bottom row refers to record 203. Narrow derivative filters are selected to enhance narrow QRS complexes (top row), while wide derivative filters are selected for records that display wider QRS complexes (middle and bottom row).  (1)).
The 15 s interval where the derivative filter is chosen is, ideally, the beginning of the record. However, it is common to find flatlines or transitory interference at the very beginning of an ECG record. For this reason, the selected 15 s interval is identified in the initialization phase of the algorithm using a signal quality index similar to the previous one that includes an additional artifact-specific term. The index is computed for successive intervals of 15 s using a non-specific raw dMAC. In this process, the algorithm looks for an interval of acceptable ECG quality. The duration of the initialization phase, 15 s, is an arbitrary number. It should be high enough to capture at least five mD windows, but The 15 s interval where the derivative filter is chosen is, ideally, the beginning of the record. However, it is common to find flatlines or transitory interference at the very beginning of an ECG record. For this reason, the selected 15 s interval is identified in the initialization phase of the algorithm using a signal quality index similar to the previous one that includes an additional artifact-specific term. The index is computed for successive intervals of 15 s using a non-specific raw dMAC. In this process, the algorithm looks for an interval of acceptable ECG quality. The duration of the initialization phase, 15 s, is an arbitrary number. It should be high enough to capture at least five mD s windows, but low enough to avoid the computational burden of an unnecessarily long SQI maximization phase.

Decision Logic
The decision logic of the QRS detection algorithm compares a dynamic threshold T with the absolute value of the derivative signal (the feature signal). Figure 6 shows the result of the processing stages described above for a sample ECG excerpt, from input signal to absolute derivative.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 9 low enough to avoid the computational burden of an unnecessarily long SQI maxim tion phase.

Decision Logic
The decision logic of the QRS detection algorithm compares a dynamic thresho with the absolute value of the derivative signal (the feature signal). Figure 6 shows result of the processing stages described above for a sample ECG excerpt, from input nal to absolute derivative. When the absolute derivative crosses the detection T, a detection is registered. first defined in the initialization phase, and then updated in a bounded fashion depend on its initial value, after each successful QRS detection. Let D ̅ 0 be the trimmed mea the derivative maxima during the initialization phase. Let M be the maximum valu the derivative signal corresponding to the last successful detection. Let us consider D ̅ a variable that reflects the average of the local derivative maxima at detection n. Sud rising spikes in M could interfere with the detection process, and therefore changes t are limited The threshold T is updated in the following fashion: The threshold T is updated in the following fashion: where c is a control parameter operating on the detection threshold. The term D n estimates the local average of the QRS complex derivative maxima at detection n. The value 0.97 that multiplies D n determines the memory factor of recursive updating. Additionally, the 0.2 D 0 limit prevents D n from diminishing to zero while the 2.5 D 0 limit prevents excessive overshoots. After each successful detection, the value of T is set to 2.4 times its new value to prevent detecting more than once the same heartbeat. This value decreases as the time elapsed from the last QRS detection increases. Occasionally, the threshold T n+1 could get excessively high and prevent detection of the following QRS complexes. For example, when a premature ventricular contraction is followed by normal sinus beats of lower amplitude. To avoid this occurrence, the threshold value T continues to decrease as the time with no detections increases. These temporal thresholds are dependent on the mean interval between two consecutive R peaks, the mean RR. As the mean RR decreases, the threshold evolves faster. Figure 7 shows the decreasing dynamic of the detection threshold for a sample mean RR interval of 1000 ms. The gray area denotes a refractory period, where no detection is registered.
where c is a control parameter operating on the detection threshold. The term D ̅ estimates the local average of the QRS complex derivative maxima at detection n. The value 0.97 that multiplies D ̅ determines the memory factor of recursive updating. Additionally, the 0.2 D ̅ 0 limit prevents D ̅ from diminishing to zero while the 2.5 D ̅ 0 limit prevents excessive overshoots. After each successful detection, the value of T is set to 2.4 times its new value to prevent detecting more than once the same heartbeat. This value decreases as the time elapsed from the last QRS detection increases. Occasionally, the threshold T +1 could get excessively high and prevent detection of the following QRS complexes. For example, when a premature ventricular contraction is followed by normal sinus beats of lower amplitude. To avoid this occurrence, the threshold value T continues to decrease as the time with no detections increases. These temporal thresholds are dependent on the mean interval between two consecutive R peaks, the mean RR. As the mean RR decreases, the threshold evolves faster. Figure 7 shows the decreasing dynamic of the detection threshold for a sample mean RR interval of 1000 ms. The gray area denotes a refractory period, where no detection is registered. Figure 7. Decreasing dynamic of the detection threshold. At time = 0, a QRS detection has just occurred. Initially, the threshold T ′ decreases from 2.4 T +1 to T +1 . During the refractory period, portrayed by the gray area, no detection is registered. As the time elapsed without a detection increases, T ′ decreases from T +1 to 0.6 T +1 , then plateaus again.
The time signature or fiducial point of the detected QRS complexes is placed on the maximum signed derivative of the QRS complex. During the initialization phase, the algorithm chooses the positive or negative derivative depending on their maximum amplitude. This choice is kept throughout the record.

Results
The detection performances are reported as F1 scores, which is a global metric of performance. F1 is the harmonic mean of Sensitivity (Se) and Positive Predictive Value (+P) Figure 7. Decreasing dynamic of the detection threshold. At time = 0, a QRS detection has just occurred. Initially, the threshold T decreases from 2.4 T n+1 to T n+1 . During the refractory period, portrayed by the gray area, no detection is registered. As the time elapsed without a detection increases, T decreases from T n+1 to 0.6 T n+1 , then plateaus again.
The time signature or fiducial point of the detected QRS complexes is placed on the maximum signed derivative of the QRS complex. During the initialization phase, the algorithm chooses the positive or negative derivative depending on their maximum amplitude. This choice is kept throughout the record.

Results
The detection performances are reported as F1 scores, which is a global metric of performance. F1 is the harmonic mean of Sensitivity (Se) and Positive Predictive Value (+P) Se = TP TP + FN (5) where TP are the true positives, FP are the false positives and FN are the false negatives of QRS detection. Figure 8 shows a strip of ECG signal, its respective absolute derivative, and detection markers for a strip of signal where true positives, false positives, and false negatives are registered. This is an example of a critical instance of signal. The first false positive of Figure 8 is caused by an artifact that closely resembles a QRS complex. This misdetection causes the algorithm to enter its refractory period, thus the next true heartbeat is missed resulting in a false negative. The second false negative is a heartbeat with very low amplitude which is challenging to spot even by eye. As this beat is missed, the adaptive threshold decreases until the second false positive is detected, despite its relatively low amplitude. While these instances of error are not common, they represent the challenging aspects of QRS detection.

R PEER REVIEW
11 of 16 Figure 8 shows a strip of ECG signal, its respective absolute derivative, and detection markers for a strip of signal where true positives, false positives, and false negatives are registered. This is an example of a critical instance of signal. The first false positive of Figure 8 is caused by an artifact that closely resembles a QRS complex. This misdetection causes the algorithm to enter its refractory period, thus the next true heartbeat is missed resulting in a false negative. The second false negative is a heartbeat with very low amplitude which is challenging to spot even by eye. As this beat is missed, the adaptive threshold decreases until the second false positive is detected, despite its relatively low amplitude. While these instances of error are not common, they represent the challenging aspects of QRS detection. Heartbeats can be divided in two groups: (1) normal beats, which are originated from the dominant pacemaker, the sinus node and (2) ectopic beats, which are generated by secondary (ectopic) pacemakers. Ectopic beats are often premature and can be atrial, supraventricular or ventricular in origin. Sometimes, sinus beats and ventricular beats collide in a fusion beat. Ectopic beats are expression of cardiac automaticity disorders and can be the hallmark of life-threatening conditions like ventricular tachycardia. The shape of these heartbeats on the ECG can differ substantially from normal beats, and are thus more challenging to detect. This is especially true for beats of ventricular origin. Moreover, conduction disorders like bundle branch blocks can also significantly alter the shape of the normal, healthy ECG pattern. Since bundle branch blocks and ventricular beats are clinically relevant for diagnosis and therapy, we also report sensitivity to these abnormal beats (Se-A) where TP A and FN A are true positive and false negatives of confusion matrix computed only by considering the heartbeats annotated as: bundle branch block, premature ventricular contractions, R-on-T premature ventricular contractions, fusion beats and ventricular escape beats.
Comparing algorithm results to annotations requires the definition of an acceptance window. This is a two-sided interval of appropriate size, centered around annotated beats. When a detected beat falls inside the acceptance window of an annotated beat it is considered a true positive. The following QRS detection results are obtained with a window of 300 ms (+/− 150 ms). This is the same window size used in Physionet's bxb annotation comparator.
To evaluate the performances of the QRS detector we introduced we implemented and operated three publicly available detectors. They are Behar's jqrs [33], Kim's [25] and Heartbeats can be divided in two groups: (1) normal beats, which are originated from the dominant pacemaker, the sinus node and (2) ectopic beats, which are generated by secondary (ectopic) pacemakers. Ectopic beats are often premature and can be atrial, supraventricular or ventricular in origin. Sometimes, sinus beats and ventricular beats collide in a fusion beat. Ectopic beats are expression of cardiac automaticity disorders and can be the hallmark of life-threatening conditions like ventricular tachycardia. The shape of these heartbeats on the ECG can differ substantially from normal beats, and are thus more challenging to detect. This is especially true for beats of ventricular origin. Moreover, conduction disorders like bundle branch blocks can also significantly alter the shape of the normal, healthy ECG pattern. Since bundle branch blocks and ventricular beats are clinically relevant for diagnosis and therapy, we also report sensitivity to these abnormal beats (Se-A) where TP A and FN A are true positive and false negatives of confusion matrix computed only by considering the heartbeats annotated as: bundle branch block, premature ventricular contractions, R-on-T premature ventricular contractions, fusion beats and ventricular escape beats.
Comparing algorithm results to annotations requires the definition of an acceptance window. This is a two-sided interval of appropriate size, centered around annotated beats. When a detected beat falls inside the acceptance window of an annotated beat it is considered a true positive. The following QRS detection results are obtained with a window of 300 ms (+/− 150 ms). This is the same window size used in Physionet's bxb annotation comparator.
To evaluate the performances of the QRS detector we introduced we implemented and operated three publicly available detectors. They are Behar's jqrs [33], Kim's [25] and Sedghamiz's [34] detectors. Behar's QRS detector features a Mexican hat pre-processing filter, maximization of the signal's energy in local windows, a fixed threshold and a backward search phase for missed heartbeats. In Kim's approach the QRS detection is performed on the signal's energy as well, by checking both the quantity and the dynamic of the QRS complex' energy with adaptive thresholds. A valid heartbeat detection requires both criteria to be satisfied. Sedghamiz's algorithm is structured on two adaptive thresholds that distinguish between signal and noise peaks, a decision logic that corrects multiple detections due to noise or T waves and a backward search procedure.
These algorithms were made available by their authors. We did not edit their code, but we did optimize their control parameter(s). Usually, each QRS detector is controlled by one or more control parameters that operate on the decision thresholds. The results we report were calculated using the control parameters that yielded the best performances for each algorithm.
Each of these programs was run in MATLAB, on the same calculator (2013 notebook with I7-dual core-cpu@2.00 GHz, 6 Gb RAM), with the same background conditions, at the same time of day, each one in succession. Table 2 reports the gross and average F1 score achieved by the QRS detectors on all the considered. No record was excluded. Table 3 shows the gross and average Se-A results. The T Wave Alternans database is not feature in this table, as it lacks annotations for non-sinus beats. The highest F1 and Se-A scores for each channel of each database and for the whole dataset are highlighted in bold. Regarding the computational cost of the algorithms, in Table 4 we show that our approach was faster than the other QRS detectors we considered. The comparison was made with the same background conditions. Each detector was run on all the databases we introduced, in succession. We report the average execution time. Our moving averagebased approach took on average three times less to process an hour of ECG recording than Sedghamiz's algorithm: 0.1137 s against 0.3388 s, effectively resulting in a reduction of 0.1137/0.3388 = 33.56% of required CPU time. The processing time of Kim's algorithm was unexpectedly high.

Discussion
In this paper, we presented a QRS complex detection algorithm based on the knowledge of the electrocardiographic phenomena. The algorithm was evaluated on well-known databases and the performances we report are comparable to the best ones published in the literature. We considered several databases, in order to avoid overfitting on a single data source. It should be noted that some published approaches reach higher performances on specific databases, but caution should be exercised in direct comparison between published results. Different papers may adopt different subsets of records from the same databases and different acceptance windows. For this reason, the performances of the algorithm we proposed was compared to three available QRS detectors. As a performance metric, we computed the F1 scores as well as execution speed. We also reported sensitivity on ventricular beats and beats where bundle branch blocks were present, which we defined as Se-A, as they are more challenging to detect. As stated above, we chose to compare available QRS detectors instead of implementing other methods from literature, which was beyond the scope of this work. Additionally, rewriting a method is a developing effort itself. The resulting algorithm is not likely to reach the same accuracy and/or computational efficiency reported by the authors, as algorithms employ a set of specific parameters and rules that the authors have fine-tuned but do not end up in the final paper. We considered available QRS detectors for a direct comparison, so that results are reported for the same set of records, with the same acceptance rule for true positives. We chose the QRS detectors of Behar [33], Sedghamiz [34], and Kim [25] as they were performing best among those we found available to our knowledge.
We evaluated the QRS detection algorithm on the MIT-BIH Arrhythmia, EU ST-T, INCART, MIT-BIH Noise Stress Test, T Wave Alternans, and QT databases. They are popular databases in the topic of automated analysis, and especially for heartbeat detection. Our algorithm outperformed the others, reaching 99.18% gross F1 and 94.81% gross Se-A. Despite this, Kim's detector obtained a higher F1 score on the first channel of the MIT-BIH database. Detection accuracy for the second leads of the MITH.BIH Arrhythmia, Eu ST-T, and NST databases are lower for each algorithm. This happens as channel 2 is often chosen to highlight specific ECG features, like ST tract elevation, while the QRS complex is prominent in channel 1.
The Se-A score was lower than the overall sensitivity, as anticipated. This is due to the different shape of this subset of beats, as the slopes of the QRS waves are low and so are the corresponding derivatives. For our algorithm the worst-case scenario of Se-V is 86.16%, for the first channel of the INCART database. This particular result is, however, substantially different from the other Se-V scores we report on the other databases (Table 4), or even the second channel of the same database. This may be caused by the annotation criteria on abnormal beats for channel 1 of the INCART database. For channel 2 of the INCART database and the other databases, sensitivity on abnormal beats lies between 94% and 99%. This means that up to 1 in every 20 abnormal beats is not detected correctly. Future work on this algorithm could address these instances of detection failure. Since identification of abnormal beats is one of the main goals in Holter exams and other lengthier ECG recordings, we feel that sensitivity on abnormal beats is a key characteristic of any QRS detector.
To highlight the light computational cost of the proposed QRS detection algorithm, we measured its execution time. On average, 0.1137 s were required to process an hour of ECG recording on our test machine, which is 33.56% of the time required by the quickest of the other programs. This result confirms the low processing power required by the algorithm we propose. While the detection accuracy of any method remains the fundamental performance metric, low execution time is nonetheless relevant for low-power devices such as wearable and portable devices.
In our algorithm, the derivative filter that maximizes the signal quality index is chosen on a selected 15 s signal interval at the beginning of the record. On longer records, e.g., 24-hlong Holter ECG, the derivative that maximizes the quality index could change throughout the record. Posture changes of the patient can alter the relative position of the heart and the electrodes, thus changing the amplitude and shape of the signal in each lead. To tackle this effect, the derivative filter selection could be repeated at regular intervals. As an alternative, the signal quality index of each derivative filter could be updated iteratively. This would remove the need to maximize the signal quality index multiple times. A derivative filter updating rule is a possible future development of this algorithm.
In the existing literature, the performance reported in some papers achieved higher scores than our algorithm. For example, Cai and Hu [16] reported higher F1 scores on the first channel of the MIT-BIH Arrhythmia (99.95% vs. 99.84%), the MIT-BIH Noise Stress Test (99.53% vs. 91.71%) and QT (99.98% vs. 99.92%) databases using convolutional neural networks. Machine learning, and specifically deep learning, have been applied with great effect to many biomedical topics, and rightly so. Yet, these approaches come at a cost in both in terms of processing time and memory requirements. Moreover, machine learning algorithms crucially depend on the training data and must satisfy specific regulatory demands to be implemented in commercial devices. The scope of our algorithm is instead to propose an accurate algorithm with an emphasis on execution efficiency and a clear, explainable decision logic.
The literature on QRS detection is rich. Many approaches have been proposed over time, and results are often reported using heterogeneous criteria. Moreover, database annotations are another source of confusion in results comparison, as human observers occasionally disagree on annotation protocols. For this reason, the question of finding the 'best' QRS detection approach has no simple, unambiguous answer and comparison of published results may not be equitable.
Elgendi [8] developed a QRS detection approach that aimed to be efficient and accurate, with moving averages, considering multiple databases, and is thus closer to the scope of our algorithm. Elgendi reported the performances for only the first channel of several databases, including the databases we considered. The comparison with our work is favorable for the MIT-BIH Arrhythmia (99.82% vs. 99.84%) and QT (99.81% vs. 99.92%) databases. Conversely, the performances of the algorithm we propose are lower on the INCART (98.05% vs. 96.51%), MIT-BIT Noise Stress Test (92.75% vs. 91.71%), and T Wave Alternans (99.00% vs. 98.20%) databases. A future comparison of execution time between these two algorithms would be of interest.
Kim [25] reported performances as high as 99.91% F1 on the MIT-BIH database, which have been confirmed in our implementation of their code (99.88% F1) and are higher than our score (99.84% F1). Despite this, on all the other databases we included their algorithm reached lower performances. This may possibly be an effect of overfitting an algorithm on a small pool of data. For this reason, comparing implemented algorithms on multiple databases can offer greater insight on the characteristics of a QRS detection algorithm. This task, however, is easier when code is shared by its authors. When an algorithm must be re-written from scratch following the pertinent literature, its performances may vary significantly from the published ones.
The algorithm we present inherits the legacy left by the work of Carlo Marchesi and Alessandro Taddei [35]. At the time, fast QRS detection was a necessity as calculators had limited processing power, which had to be optimized. Today, fast QRS detection remains a relevant topic for all applications involving battery-driven or portable devices.

Conclusions
In this paper, we presented a knowledge-based QRS detection algorithm that emphasizes both computational efficiency and accuracy. This algorithm features a preprocessing stage in the time domain, a decision logic based on medical knowledge, and an adaptive derivative behavior. While more accurate approaches exist in literature, this algorithm could be a useful tool in performing on-line ECG heartbeat detection in limited CPU environments.  Data Availability Statement: Each database we included in this work is freely accessible on the Physionet online repository, https://www.physionet.org/about/database/ (accessed on 28 July 2021). The MATLAB code of the QRS detector algorithm proposed in this paper can be received by contacting any of the authors.

Conflicts of Interest:
The authors declare no conflict of interest.