The electrocardiogram (ECG) is a quasi-periodic biomedical signal that provides information about cardiac muscle electrical activities. One cardiac cycle in a typical ECG signal is delineated by arrangements of P, the QRS complex, T waves, and PQ and ST segments. Correct R-peak detection is the first and most critical step in almost all ECG analysis methods. The R-peak is the highest and only positive peak within the QRS complex, reflecting the ventricular depolarization of the heart’s electrical activity. Precise detection of the R-peak location plays a critical role in obtaining the morphology of the QRS complex and revealing the location of other ECG fiducial points. Furthermore, R-peak localization serves as the basis for automated determination of the heart rate, which is a significant criterion for heart arrhythmia diagnoses such as premature atrial contraction, tachycardia, and bradycardia. Many other diseases can also be diagnosed in a non-invasive way using R-peak detection due to the relationship between heart rate variability and several physiological systems (e.g., vasomotor, respiratory, central nervous, and thermoregulatory).
Various approaches have been proposed in the literature for detecting R-peaks in an ECG signal [beraza2017comparative]. Typically, these methods consist of two main steps: pre-processing and detection. In the pre-processing step, the algorithm attempts to eliminate the noise and artifacts and to highlight the relevant sections of the ECG [castells2015simple, sharma2017qrs]
. In the second step, various methods are used to locate R-peaks based on the result of the pre-processing step, and then other waves are detected by defining a set of heuristic rules[hou2018real]. However, these approaches suffer from some critical drawbacks that limit their performance in practical applications. First, in real-time data processing and ambulatory care settings, where the collected data are highly noisy, preprocessing-based algorithms are less effective. Second, these algorithms can fail to detect R-peaks in some determinant morphological patterns resulting from certain life-threatening heart arrhythmias due to the time-varying morphology of the QRS complex. Incorrect detection of R-peaks can affect the correct identification of subsequent waves.
The R-peak detection step can be generally accomplished either by implementing a threshold-based technique or by employing an independent threshold technique. The amplitude of the peak and time duration between two consecutive R peaks (i.e., the RR interval) are typically used to determine a suitable threshold [song2015new]. A constant threshold is only efficient for detecting R-peaks within records with normal morphological patterns. Therefore, recent studies have employed adaptive thresholds, for which there is no need to determine the threshold experimentally. In [nayak2019optimally] and [sahoo2016noising], the Hilbert transform with an adaptive thresholding technique was utilized to detect R-peaks. Some threshold-based techniques with other criteria have also been used to specify the threshold. In [hou2018real], an adaptive threshold concerning the geometric angle between two consecutive samples of the ECG signal was defined. The performance of the threshold-based technique is highly dependent on the selection of initial parameters; hence, it can lead to a significantly higher number of false beats. Therefore, independent threshold techniques are more desirable than the threshold-based technique.
Most of the state-of-the-art methods for R-peak detection are based on wavelet transform [park2017r, rakshit2017efficient, thirrunavukkarasu2019detection], simple mathematical operations [nayak2019optimally, PanTompkins1985, gutierrez2015novel]
, hidden Markov models, and machine learning. Wavelet transform is a suitable approach for considering the non-stationary behavior of the ECG signal. However, considering the various shapes of the QRS complex, it is difficult to select the optimal mother wavelet or find the required threshold in the detection step of the wavelet transform. Additionally, discrete wavelet transform fails to provide reliable results in a short-recording duration. Mathematical operation-based algorithms have a low computational cost, which is more appropriate for real-time applications and large dataset analysis. However, achieving high performance when the signal-to-noise ratio is high remains challenging for these algorithms. Hidden Markov models are also widely used in ECG segmentation because they are powerful tools for considering the temporal dependency among the waveforms[akhbari2016ecg, monroy2020hidden, huque2019hmm]. The majority of the studies on machine learning-based methods have utilized sparse signal processing to represent an approximation of the nonlinear ECG signal using sparsity constraints [ning2013ecg, zifan2006automated, parviziomran2019optimization, tavakoli2013decoding, monti2018adaptive, tibshirani1996regression]
. Some studies have also applied deep learning techniques to detect the ECG waveforms considering its high performance in various classification tasks[xiang2018automatic, tavakoli2017learning]. However, the caveat with deep learning-based approaches is that they need large-scale datasets for the training phase and often suffer from the imbalanced class problem [8683140, mousavi2020single].
In this paper, we propose a new class of graphical models based on optimal changepoint detection models, named the graph-constrained changepoint detection (GCCD) model, to locate R-peaks in the ECG signal. A changepoint detection model identifies abrupt changes in data when a property of the time series changes. In the non-stationary ECG signal, ECG waves can also be considered as abrupt up or down changes over time during the heart cycle. We exploit the model introduced by Hocking et al. [hocking2017log, hocking2018generalized], in which a graph-based optimal changepoint detection model was used for detecting abrupt changes in the genomics data. In their work, they propose a new class of functional pruning algorithms with log-linear time complexity in the amount of data, which is capable of handling the large datasets that are common to ECG analysis.
Only a few studies in the literature have applied changepoint detection models for cardiac analysis. Gold et al. [gold2018doubly]
adopted a changepoint detection method based on Bayesian inference to extract the onset of the QRS complex over a small time window containing just one QRS complex. In[qi2014novel], a changepoint detection approach based on the Haar wavelet and Kolmogorov-Smirnov statistic was applied to find normal and abnormal ECG segments within the assembled ECG samples from different ECG datasets. Sinn et al. [sinn2012detecting] analyzed heart rate changes in ECG recordings by detecting abrupt changes in the ordinal pattern distributions, which are used to represent the order structure of a time series. Some studies have also applied changepoint detection models to investigate sleep problems by analyzing heart rate variability in the ECG signal during sleep [staudacher2005new, ma2020multiple].
To the best of our knowledge, this is the first study in which changepoint detection models have been proposed to detect ECG fiducial points in long records of ECG signals. In this novel framework, prior biological knowledge about the expected sequence of changes can be specified in a constraint graph. Then, functional pruning dynamic programming algorithms can compute the globally optimal model (mean, changes, and hidden states) in fast log-linear time. We furthermore propose a new algorithm for learning the graph structure using labeled ECG data. Therefore, the main contributions of this study are:
A new class of graphical models based on optimal changepoint detection models to detect R-peak positions in the ECG signal. The proposed method does not require any noise removal preprocessing step as it uses the sparsity of changepoints to detect abrupt changes.
A new algorithm to learn the graph structure and parameters using labeled ECG data. Thus, the model’s performance is no longer dependent on an expert to encode prior knowledge into the constraint graph.
Comparison of the learned graphs with the manually constructed graphs in terms of graph structure and detection accuracy. Results demonstrate that there can exist different optimum graph structures for one subject, and the proposed graph learning algorithm can find global optima depending on the initial graph structure.
The rest of the paper is organized as follows. In the next section, we describe the proposed model for R-peak detection in the ECG signal. We explain the GCCD model in Section II-A and the constraint graph in Section II-B. Section II-B also defines the manual graph and the proposed graph learning algorithm. Section III provides a description of the dataset used in this study and a discussion of the results as well as a comparison between the performance of the manually defined graphs and learned graphs. Finally, Section IV summarizes this research work and its contributions.
The proposed method treats ECG wave detection as a changepoint detection problem for a non-stationary ECG signal. It extracts the R-peaks in the raw ECG signal by representing the periodic non-stationary ECG signal as a piecewise locally stationary time series with constant mean values (i.e., each piece is the mean of one segment of datapoints). The model takes a raw ECG signal and a constraint graph as inputs and computes the onset/offset and the mean of desired segments (i.e., hidden states). Then, the center of each state is associated with the location of a peak. The constraint graph allows the incorporation of prior knowledge into the model and regularizes the model. Figure 1 illustrates an overview of the proposed algorithm in the detection of R-peak positions in the ECG signal. It is worth re-emphasizing that the model takes the raw ECG signal as the input, without applying any preprocessing step, as it leverages the sparsity of changepoints to denoise the signal and to detect abrupt changes.
The constraint graph, which encodes the expected sequence of changes in the ECG signal, can be defined manually by an expert or automatically from the data. In the following sections, we describe the details of various parts of the proposed model.
Ii-a Graph-Constrained Changepoint Detection Model
ECG fiducial points detection can be defined as the problem of finding abrupt changes over one cardiac cycle caused by changes in statistical characteristics. From this point of view, a proper changepoint detection algorithm can be employed to detect ECG waves in a fast and effective way. We applied the optimal changepoint detection model introduced in [hocking2017log] to localize R-peak positions in the ECG signal. In this model, prior biological knowledge about the expected sequence of changes is incorporated into the model as a graph constraint. Then, a dynamic programming algorithm using functional pruning computes the globally optimal model (mean, changes, and hidden states) in fast log-linear time.
We assumed a directed graph as the constraint graph, where the vertex set represents the hidden states/segments (not necessarily a waveform), and the edge set represents the expected changes between the states/segments. Each edge incorporates the following associated prior knowledge about the expected sequences of changes:
The source and target are vertices/states for a changepoint from to .
A non-negative penalty constant is the cost of changepoint .
A constraint function defines the possible mean values before and after each changepoint . If is the mean before the changepoint and is the mean after the changepoint, then the constraint is . These functions can be used to constrain the direction (up or down) and/or the magnitude of the change (greater/less than a certain amount).
Mathematically, given the input signal and the directed graph , the problem of finding changepoints , segment means , and hidden states can be solved using the following optimization problem:
The changepoints can be assigned to any of the pre-defined edges (). Consequently, indicates no change with zero cost, . Function (1) consists of a data-fitting term and a model complexity term [amini2019iterative, amini2019iterative2]. represents the negative log-likelihood of each datapoint, and is a non-negative penalty on each changepoint. In other words, regularizes the number of predicted changepoints/segments by the model so that a largeralso encodes the expected up/down change and the least amplitude gap between the mean of two states. When there is no change , Constraint (2) forces the model to stay in the current state with no change in mean . However, when there is a change , Constraint (3) imposes a change in the mean implied by the constraint function as well as a change in the state
. An open-source implementation of the Generalized Functional Prunining Optimal Partitioning (GFPOP) algorithm is available in C++ code inside an R package named GFPOP on GitHub[GFPOP].
Ii-B Constraint Graph
The constraint graph in the optimization problem of Equation (1) encodes prior biological knowledge about the expected sequences of changes within one cardiac cycle. It can be designed manually by an expert or be learned from the data by the model. The two following subsections detail both the manual and learning-based designs.
Ii-B1 Manual Graph Definition
To manually define the constraint graph , we took into account the possible morphological categories for the ECG waves (i.e., P, QRS, and T waves) and the overall morphological properties of the signal in each record [fotoohinasab2020graph]. An expected hidden state/segment in the signal is characterized as a node in the constraint graph, and the required conditions for transition between states are encoded in the edges. The required conditions are determined based on the expected minimum amplitude difference of two successive states and the polarity of each transition (i.e., up/down).
The caveat with the manual definition of the constraint graph is that it can be inefficient for ECG signal analysis considering the various morphological patterns for each waveform. Furthermore, the model’s performance depends on the expert knowledge encoded into the constraint graph. In the next subsection, we explain the proposed graph learning algorithm for learning the constraint graph using the R-peak labels provided by the gold standard.
Ii-B2 Constraint Graph Learning
To automate the R-peak detection task, we modified the previous model by learning the constraint graph from the data (see the dashed part in Figure 1). In this new framework, the proposed model takes the raw signal and an initial graph structure as inputs and yields the desired outputs, including the onset/offset and the mean of segments specified in the nodes of the learned constraint graph [fotoohinasab2021graph]. Here, the model architecture is comprised of two stages: training and detection. The training step tries to heuristically find an optimum graph structure by which the label errors in the training set are minimized (the block named “Graph Learning Algorithm” in Figure 1). The detection step then extracts the R-peaks in the raw ECG record constrained to the graph learned in the previous step (the block named “Changepoint Detection Model” in Figure 1).
The novelty of this new structure lies in the training step, which is comparable to the previous model in Section II-B1. The main idea of the training step is to automatically discover the desired topology of the constraint graph and the information about the edges from the data. As described in Section II-A, each edge contains the following information: (1) the expected up/down change in the segment means, (2) the least amplitude gap between the means of two states, and (3) a non-negative penalty imposed by the edge transition. Suppose that the initial graph for each record is denoted as , where and are the corresponding graph node and edge sets, respectively. Each node in the set represents initial hidden states in the model. Each edge in the set represents a transition between two consecutive hidden states (i.e., a changepoint from the source to the target in section II-A) and also contains initial values for parameters of , , and , which are the initial type, the initial gap between two states, and the initial penalty, respectively. Figure 2a shows the simple initial graph used for the optimization process. It should be noted that the initial edge information was chosen based on the overall results obtained from the manual definition of the constraint graph.
A sketch of the proposed graph learning algorithm is summarised as Algorithm 1. The greedy graph search algorithm starts with the initial graph and iteratively optimizes the graph structure and edge parameters to find a graph that maximizes the accuracy regarding the provided labels. At the -th iteration, the function finds the graph candidate set using the editing candidates for each edge of the output graph from the previous iteration . In this study, the algorithm considers 11 editing candidates per edge to optimize the graph topology and the three edge parameters. For example, in the iteration , if the parent graph (i.e., ) has two edges, the graph candidate set will have no more than 22 members . These editing candidates include three types of adding a node, two types of deleting a node, one type of adding two nodes, changing the type of the abrupt change, and increasing or decreasing the penalty and gap corresponding to an edge. We believe all morphological patterns of the ECG waves can be constructed using these editing candidates. Figure 2 illustrates the graph editing candidates related to the edge with an up change.
Ii-C Computational Complexity
As can be seen in Algorithm 1, the time complexity of the GCCD algorithm is theoretically proportional to the number of graph candidates at each iteration (Line 10) and the number of required iterations to achieve an optimum graph with minimum label errors (Line 5). There are three main aspects that characterize the time complexity of the algorithm:
Given a record with data samples and a graph candidate with vertices and edges, the time complexity to detect R-peaks (Lines 11–15) is in the worst case (pathological simulated data) and in the average case (typical in real data). Also note that since we consider only graphs with a single circular path, , and the time complexity is further reduced to (for average case/non-pathological data).
Considering graph edit candidates in the iteration , the time complexity to compute all the models is (where is the time complexity of solving for optimal model parameters given a single graph). It should be noted that the number of graph candidates in the iteration depends on , which is the graph from the previous iteration (Line 10). The time complexity to compute the label error given labels is , which can be effectively ignored from the overall time complexity as this task is fast.
Finally, iterating over iterations to obtain the graph with the minimum label error (Line 5) causes the overall time complexity of the algorithm to be , where is the time to solve for a single graph, and is the number of edit candidates considered in each iteration.
Iii Experimental Studies
We applied the well-known MIT-BIH Arrhythmia (MIT-BIH-AR) database to evaluate the GCCD model. This database contains 48 ECG recordings taken from 47 subjects. Each record’s duration is 30 min, and each recording is sampled at 360 Hz with a resolution of 200 samples over a 10 mV range [moody2001impact, goldberger2000physiobank]. Each recording consists of two ambulatory ECG channels from the modified lead II (MLII) and one of the leads V1, V2, V4, or V5. In this study, all 48 records with one MLII or V5 lead were used to evaluate the algorithm. The database has been annotated with both RR intervals and heartbeat class information by two or more expert cardiologists independently.
Iii-B Results and Discussion
This section presents a comprehensive discussion of the results obtained by the proposed model and a detailed comparison between the manually defined graphs and the learned graphs. We also provide some suggestions for the future development of the GCCD model.
Figure 3 illustrates an example of the model’s performance with a manually defined constraint graph in the R-peak detection task for a window of Record 230 of the MIT-BIH-AR dataset. However, as mentioned in Section II-B1, the performance of the model using manually defined graphs depends on an expert with prior knowledge. Furthermore, manual annotation by an expert is time consuming and expensive. To address this issue, we proposed a new graph learning algorithm that searches for a locally optimal constraint graph using a greedy scheme on the labeled ECG data. Regarding the various morphological patterns for the ECG signal, the proposed graph learning algorithm can relax the model from the manual definition of the constraint graph for each record.
We adopted the intra-patient paradigm to train the constraint graph to address the intra-patient variation in ECG morphologies. Thus, the training and testing sets were generated by randomly splitting the intra-samples for each record with an approximate ratio of . We used a -fold cross-validation approach to evaluate the model performance with a size of . More specifically, we divided the intra-sample data into folds so that each trial used four folds to train the model and one fold for validation.
Figures 4 and 5 show representative examples of the R-peak detection task performed by the model integrated with the graph learning algorithm for two records from the MIT-BIH-AR database. These figures illustrate how the proposed graph learning algorithm iteratively edits the graph structure to yield a model with maximum accuracy in detecting R-peaks. We initialized the constraint graphs using the graph structure in Figure 2a with the initial values of and for Record 107 and and for Record 219. It should be noted that the initial edge information was assigned based on the overall results derived from the manually defined graphs in all experiments. However, graph candidates 2f and 2g can adjust the parameters and for the optimum values. For these two examples, we chose the initial edge information so that all the training steps could be completely displayed. Label errors are omitted from Figures 5a–5c to reduce clutter in the figures. The red part of the graph structure in each iteration presents the chosen editing candidate in the current iteration over the graph in the previous iteration. More interestingly, Figure 5 demonstrates the model’s capability to detect R-peaks in the presence of a baseline wandering artifact, which is a typical artifact in the ECG signal. Baseline wandering can change the shape of the QRS complex and thereby causes incorrect detection of the R-peak. The performance of the Pan and Tompkins [PanTompkins1985] algorithm, algorithms derived from the ECG signal slope [sabzevari2018ultra], and methods based on wavelet transform are highly dependent on the removal of this artifact. Figure 6 shows the test result for this record over two different time windows of data. Figure 7 illustrates the training progress for these two records, where the Y-axis shows the sum of false negative and false positive error rates. Indeed, the training progress curve reflects the number of label errors produced by the model in each iteration given the provided labels for the training and validation sets. It is worth mentioning that the proposed graph learning algorithm avoids possible overfitting issues as it tries to extract the morphology of the ECG signal that contains multiple various morphological patterns.
The proposed graph learning algorithm employs a greedy search scheme to select the best performing graph in terms of detection accuracy (see Section II-B2). Therefore, the performance of the algorithm depends heavily on the initial graph structure and will likely lead to local optima. Figure 8 compares the training progress for Record 106 of the MIT-BIH-AR database initialized with the two simple (see Fig. 2a) and complex graph structures (i.e., a graph with eight nodes representing the morphology of a normal ECG signal). Figure 9 also presents a comparison of the final selected graphs and their performances for a window of this record. As these figures show, the model initialized with the complex graph structure can achieve higher accuracy (i.e., a lower number of label errors) in a lower number of iterations than the model initialized with the simple graph structure.
The investigation of the experimental results shows that the greedy graph search algorithm can achieve optimal performance for the model trained with the manually defined graphs, although its performance is affected by the initial graph. We noticed that for most of the records from the MIT-BIH-AR database, the learned graphs could reach the performance of the manually defined graphs but with different graph structures. This means that the GCCD model can obtain global optima using various initialization structures, which will likely lead to different final graph structures. Figure 10 compares the constraint graph structures defined manually vs. those learned automatically using the initial graph structure in Figure 2a for Record 100 of the MIT-BIH-AR dataset. As shown in this figure, the manually defined graph and the learned graph both achieved the optimal performance but with different graph structures. We also noticed that for some records from the MIT-BIH-AR database, the graph learning algorithm chose the same structure as the manually defined graph structure. Figure 11 shows the model performance using the graph learning algorithm for Record 232 from the MIT-BIH-AR dataset, for which the manually defined constraint graph and the learned graph had the same structures.
Different metrics were adopted to evaluate the performance of the proposed model with both the manual and learning-based graph designs. These metrics included sensitivity (), positive predictivity rate (), and detection error rate (), which are calculated by:
where is true positives, is false positives, is false negatives, and is true negatives. Table I presents the performance of the proposed model regarding both the manually defined and learning graphs against other state-of-the-art methods for R-peak detection (QRS complex). As shown in the table, the proposed algorithm achieved Sen = %99.76, PPR = %99.68, and DER = 0.55 for the manual definition of the constraint graph and Sen = %99.64, PPR = %99.71, and DER = 0.19 for the learning constraint graph using the MIT-BIH-AR database. Note that the model constrained to the manually defined graphs outperformed the model combined with the graph learning algorithm because in the latter, the model’s performance was largely dependent on the initial graph structure.
|Method||Sen (%)||PPR (%)||DER (%)|
|Park et al. [park2017r]||99.93||99.91||0.163|
|Sharma and Sunkaria [sharma2017qrs]||99.50||99.56||0.93|
|Castells-Rufas and Carrabina [castells2015simple]||99.43||99.67||0.88|
|GCCD Model with Manual Definition of the Constraint Graph||99.76||99.68||0.55|
|GCCD Model with Learning of the Constraint Graph||99.64||99.71||0.19|
ECG recordings in the MIT-BIH-AR database were chosen to challenge the R-peak detection task because they represent a wide variety of QRS morphologies with real-world variability. Our proposed model yielded outstanding results when detecting R-peaks in these tricky records. Records 103, 104, 105, 108, 111, 112, 116, 200, 201, 203, 205, 208, 210, 217, 219, 222, and 228 are comprised of abrupt changes in ECG morphology, and they are severely affected by noise and artifacts. Figure 5 shows the capability of the model to detect R-peaks in the presence of baseline wandering noise. We re-emphasize that these comparable results were obtained without applying any preprocessing operations, as opposed to other methods in the literature. Records 108, 113, 117, 201, 202, 203, 213, 219, 222, 223, 231, and 232 contain many peaks with unusual amplitudes. Small-amplitude R-peaks or high-amplitude P- and T-peaks embedded in high-amplitude QRS complexes can lead to high FN and FP errors in the R-peak detection task. As a representative example, Figure 4 illustrates the efficiency of the GCCD model in R-peak detection for Record 117, which contains many beats with high-amplitude T-peaks.
The experimental results obtained using the proposed model justify changepoint detection models as a potential approach to extract ECG fiducial points. In this study, we demonstrated the capability of the GCCD model in locating R-peaks within various morphological patterns of ECG. The proposed greedy graph search algorithm can potentially detect ECG waves other than the R wave (i.e., P, Q, S, and T waves) by considering corresponding prior knowledge of the graph editing candidates. We noticed that in Records 114, 200, 203, 207, and 210, the Sen and PPR values were less than 99%. These records contain multiple different morphological patterns, including negative QRS complexes, and Records 200 and 203 have several QRS complexes with ventricular arrhythmias. The constraint graph for these records involves learning a graph with more than one optimum graph path. Learning a multi-path constraint graph is also required to detect all ECG waves due to the various morphological patterns of each wave incorporated into the graph. The other point that should be considered here is that the GCCD model estimates the ECG signal using a Gaussian function. A modified model with a multi-Gaussian fitting method can drastically improve the ECG-related changepoint detection task.
Future work should focus on developing the proposed model with a multi-Gaussian fitting and a multi-path graph learning algorithm. Incorporating these modifications into the proposed model could provide a promising platform for evolving new graph-based tools to detect and classify heart arrhythmias. A multi-path graph learning algorithm could reveal the morphology of the ECG signal (time duration, amplitude, and direction of each wave) in each cardiac cycle. Subsequently, new graph-based features could be extracted from the constraint graph path for an ECG cycle to classify heartbeats.
The accurate delineation of R-peaks in the ECG signal plays a crucial role in most automated ECG analysis tools. This paper proposed a novel graphical model based on changepoint detection techniques for detecting R-peaks within a non-stationary ECG signal. The proposed model was highly successful at detecting R-peaks in noisy ECG data without applying any preprocessing steps. To our knowledge, this is the first time that a changepoint detection model has been applied for ECG fiducial points detection. In this new framework, prior biological knowledge about the expected sequences of changes was incorporated into the model using a graph. We defined the constraint graph manually and automatically using a proposed greedy graph search algorithm. Using the proposed graph learning algorithm, the initial graph structure can develop into a structure containing edge parameters with maximum detection accuracy for a record. The experimental results provided in this paper demonstrate that the GCCD model can be a promising approach for detecting ECG waves and developing new graph-based tools for further ECG analysis. The proposed graphical model approach can be advanced by learning a multi-path constraint graph and fitting a multi-Gaussian curve model to the ECG signal, which should be considered in future studies.
This material is based on work supported by the National Science Foundation under Grant Number 1657260. Research reported in this publication was supported by the National Institute on Minority Health and Health Disparities of the National Institutes of Health under Award Number U54MD012388.