Parallel use of a convolutional neural network and bagged tree ensemble for the classification of Holter ECG

We used a public dataset from the PhysioNet/CinC Challenge 2017—single-lead Holter ECG, acquired from a device (AliveCor Inc., CA, US) in 12-bit depth and with a sampling frequency of 300 Hz. The dataset contains files categorized by pathology into four classes. The expert labeled dataset was divided into a public (8528) and hidden (3658) part. We removed 345 files from the public part in cases in which we disagreed with the expert labeling. Therefore, a total of 8183 files was used for training and local testing. The hidden part of the dataset was used for remote tests on PhysioNet servers. Labeling corrected after the end of the 'Official Challenge Phase' was not used for training.

The ECG signal is transformed into several amplitude envelopes (envelograms) which are later used for QRS detection and feature extraction (figure 1(A)). Envelograms are computed using Fourier and Hilbert transformations in specific frequency ranges. Envelograms in lower frequencies (LF, 1–6 Hz, figure 2(C)) and middle frequencies (MF, 8–30 Hz, figure 2(B)) serve as sources for QRS detection (figure 1(B)—dots). The LF envelogram is more sensitive to T-waves, while the MF envelogram is more sensitive to QRS and this is used by specific conditions during QRS detection. These envelograms are also used for distinguishing ventricle contractions from regular QRS complexes using ratio of LF and MF envelograms [29]. The input ECG signal is also filtered into the 2–30 Hz band (CF signal) to remove baseline wandering and noise; the signal transformed in this way is prepared for QRS morphology clustering (figure 1(C)).

Zoom In Zoom Out Reset image size

For QRS detection (figure 1(B)), peaks were detected in the MF envelogram. Detection criteria were set to 'Descend' order and the minimal distance to other peaks was set to 150 ms (the maximal detectable heartrate is, therefore, 400 bpm). Each peak must then pass multiple checks. First, each peak is checked to see whether the MF signal maximum in a window  ±300 ms from the peak exceeds 1.5×  nominal peak value. This time window was set in respect to the maximal expected P-R distance and R-T (top of the T-wave) distance. If this check is not passed, it means that it is probably a P-wave or a case of a very sharp top of a T-wave or a noise; such peaks are removed from further processing. Next, the same time window around the peak is divided into three equal and consecutive parts (figure 2—R_A, R_B and R_C, each of a length of 200 ms) and the standard deviation inside these parts is compared. If the sum of the standard deviation of border regions R_A and R_C is higher than the standard deviation of R_B the peak is removed from further processing. Finally, each peak is checked for low-frequency activity on its borders using an LF envelogram, because we expect that low-frequency activity connected to any QRS complex should follow QRS due to a T-wave. Therefore, the sum of region  −250:0 ms from the tested peak (figure 2(C)-L) is tested against region 0:250 ms (figure 2(C)-R). If sum(L)  >  sum(R)  ×  1.5, the low-frequency activity in LF before the peak is much higher than after the peak and it is unlikely to be a QRS. Therefore, the peak is considered a false QRS and is removed from further processing. The time windows of 250 ms were set in respect of the capture of areas in LF envelopes related to the T-wave. However, in the case of a very fast pace it is possible that both T-waves prior and post the tested peak will be in the analyzed  ±250 ms area. For this reason, we had to add a 1.5×  multiplier to cover such a situation. Furthermore, during events such as ventricle tachycardia it may happen that the majority of the detected QRS are removed. Therefore, this step is omitted if more than 95% of all peaks are removed during this check. Finally, if the number of detected QRS is insufficient (<4), the process is aborted and the record is considered noisy.

Because the task is also to find other arrhythmias in addition to atrial fibrillation, we decided to cluster QRS complexes using their morphology (figure 1(C)). The reason for this is that atrial fibrillation, as well as the sinus rhythm, should contain a single QRS morphology, while premature ventricle contractions, for example, should result in additional morphological groups. Clustering is achieved using the Pearson correlation, computed in a  ±150 ms window centered at QRS complexes; correlation is computed using the ECG signal filtered in the 2–30 Hz band (CF). The size of the correlation window was set in respect to QRS duration, which is not expected to be longer than 300 ms.

The clustering process works as follows: the first QRS complex is assigned to the first morphological group. The next QRS complex is compared to all already assigned QRS complexes (called prototypes). If they correlate sufficiently (>0.9), the tested peak is assigned to the same morphological group as the prototype. It is typical of noisy signals that the number of detected morphological groups is very high, while signals containing a high-quality signal (even with a pathological content) should contain 2-3-4 morphological groups. Therefore, if an excessive number of morphological groups is found (>length of signal in seconds), the process is aborted and the record is considered noisy.

There is also a risk that the same QRS morphology will be registered multiple times, but with a time shift. To avoid this, an average shape for each morphology group is created and their central parts (±75 ms) are correlated with a shift (also  ±75 ms). If any two groups correlate well enough (>0.95), they are merged together; averaged shapes must be regenerated and the process continues until none of the groups can be merged with another.

Morphological groups are ordered by a count of QRS complexes. Therefore, the first morphological group (G1) should always describe a sinus rhythm (with the exception of bigeminy PVCs). G2 will usually contain the most common morphology of extra beats, G3 will contain a less common morphology, etc. Obviously, the presence of noise will create additional morphological groups, but these should not contain multiple QRS complexes because the noise is expected to be random.

Once morphological groups are defined, we can look for missing QRS complexes (which have not, for some reason, passed the initial checks). The search for missing QRS complexes works as follows: a cross-correlation function is generated from the averaged shape of each morphological group (containing more than one QRS complex) and the CF signal. Next, peaks are detected in this cross-correlation function (>0.94). Next, nominal values in the position of the detected peaks are compared to corresponding samples in averaged shapes of morphological groups. Finally, if a newly detected peak passes these checks and its position is further than 150 ms from a QRS already detected, the new QRS complex is defined and added to the specific morphological group.

Clustering based on the Pearson correlation may lead to incorrect results in cases in which shapes correlate sufficiently but have a different scale. In terms of ECG, beats represented by the same appearance of QRS, but with 0.1×  amplitude, probably do not belong to the same morphological group and are probably the result of noise. Therefore, the median sum of the region  ±150 ms around each QRS complex from group G1 is computed in an MF envelogram (Med_MF). This value serves as a base for comparison of each other QRS complex: this value of each QRS complex must be between Med_MF  ×  0.2 and Med_MF  ×  4 or it is excluded from the list of QRS complexes. Finally, morphological groups are reordered again using their number of QRS complexes.

A list of 43 extracted features for the bagged tree ensemble is shown in table 1. The first three features (table 1, #1–3) are not dependent on QRS detection and cover the hypothetical situation in which, for example, a block of ventricular fibrillation prohibits the QRS detector from detecting single beats. For this, we use a previously developed approach [29]. The longest block when the LF envelogram is dominant over the MF envelogram is evaluated in terms of length (s) (table 1 #1), MF to LF area ratio (table 1 #2) and the sum of the LF envelogram (table 1 #3) form additional features describing the same issue.

All other features used for the bagged tree ensemble are dependent on detected QRS complexes and their morphological group. Statistical description of the RR intervals forms seven features (table 1 #4–11); the number of peaks in consecutive RR values, as well as the number of RRs higher or lower than average RR (outside  ±0.1 tolerance), are additional features (table 1 #11–13).

Thanks to morphology clustering it is possible to quantify ratios of other morphological groups as well as the ratio between the most important groups G1 and G2 (table 1 #14–17). We also examined the whole sequence of RR intervals in a sliding 6-beat window, where seven features were extracted from the standard deviation and mean RR in this sliding window (table 1 #18–24).

Frequency components of all beats are described in five features (figure 3 and table 1 #25–29), where the ratio of LF and MF envelograms is sensitive to abnormal beats, specifically premature ventricle beats. This is also extracted for the average shape from G1, resulting in a specific feature (table 1 #36).

Zoom In Zoom Out Reset image size

Technically, most other arrhythmias as well as atrial fibrillation result in irregularities in RR intervals. Therefore, another feature (table 1 #30) describes the ratio of unusual RR intervals to all RR intervals (an RR change greater than 10% is considered an unusual interval). This is also covered by examining the first derivative of the RR series; five features (table 1 #31–35) statistically describe changes of RR intervals. The last feature describing changes in RR intervals (table 1 #40) is computed as the correlation of the whole RR series with its copy shifted by one RR interval.

The behavior of the P-wave, or its link to the QRS complex, is essential to stating pathology—namely atrial fibrillation, 2nd or 3rd degree AV-block, or atrial tachycardia. However, the detection of the P-wave is a difficult task because it may be buried under background noise. This is common in Holter recordings. Therefore, we decided not to detect the P-wave directly, but to analyze the region in which the P-wave should exist. We statistically describe the stability of this region preceding each QRS complex (figure 4(D)) in three features (table 1 #37–39).

Zoom In Zoom Out Reset image size

The last features (table 1 # 41–43) used for the bagged tree ensemble describe a number of extrema computed on the averaged shape from the most common morphology G1. Three regions are analyzed—preceding, surrounding and following the average QRS complex (figures 4(A)–(C)).

If any of these 43 features cannot be evaluated or if the mean RR interval (computed using the 1st morphological group only) is longer than 3 s, the process is aborted and the recording is considered noisy. Otherwise, these features are fed into the bagged tree ensemble. This condition correctly detects 21% (N  =  58) of noisy recordings and incorrectly detects 3% (N  =  9) in the training set.

A bagged tree ensemble was chosen because it has shown the best results of the tested machine learning methods (simple decision trees, shallow neural networks and supported vector machines with different kernels). 70% of the public dataset was used for training; the remaining 30% was used for testing. The tree ensemble consisted of 30 trees. The bagged tree ensemble was trained using the Machine Learning and Statistics toolbox in Matlab software (r2017a). The resulting confusion matrix is shown in table 2.

The presented method combines two—mostly independent—approaches. While the approach with a bagged tree ensemble relies on QRS detection, the convolutional neural network uses only filtered ECG data. The input for the CNN consists of 6 s-long (without zero padding) blocks of input ECG signal. These blocks are pre-processed with zero-phase bandpass digital filtering (Butterworth filter) and envelograms in ranges (1–5, 5–10, ..., 35–40 Hz) are estimated using the absolute value of the Hilbert transform. Furthermore, the RAW signal and computed envelograms are concatenated into a matrix that is used as an input image for CNN.

The CNN architecture consists of 13 layers. The input image of raw signal and envelograms (dimensions 9  ×  1800) is processed with the first convolutional layer with kernel dimension (9  ×  150) and 30 filters. Next, a non-linearity mapping layer with a ReLU activation function is used. Next, the extracted features are down-sampled with a max-pooling layer (1  ×  2). Subsequently, a second convolutional kernel (1  ×  15  ×  30) with 50 filters, ReLU and dropout (0.5) is used in order to extract more complex features. The resulting tensor is squashed into vector form and consequently a fully connected layer (15 neurons) with ReLU and dropout (0.5). Finally, the fully connected layer (three neurons) with softmax activations is attached and probabilities for each class are obtained.

The proposed CNN system is trained for classification into three groups, i.e. normal rhythm, atrial fibrillation and other arrhythmias. CNN was not trained for noisy recordings due to their low number. A gradient descent algorithm with momentum (100 training examples per minibatch) was used as a training method and cross-entropy error was used as a cost function. The training process was performed for 25 epochs. We have used a data augmentation technique (ECG lead inversion) in order to train CNN for the classification of signals with inverse polarity. On the other hand, we did not see any significant improvement after this augmentation and, therefore, this step should be reconsidered.

Because we used a short 6 s window for analysis by CNN, its output is a probability vector of being N, O or A (class X was not recognized by CNN due to the low number of X recordings). A typical situation for a record containing one premature ventricle beat (resulting in class O) is that most of the record points to class N, but in terms of the location of PVC it points to class O (figure 5). This means that quantifying simple maxima may be insufficient for encoding CNN results. Therefore, we decided to extract several statistical descriptors from CNN output vectors—${{\overrightarrow{{\rm CNN}}}_{N}}$, ${{\overrightarrow{{\rm CNN}}}_{A}}$ and ${{\overrightarrow{{\rm CNN}}}_{O}}$. These descriptors are used as input features for a shallow neural network. A list of these features is presented in table 3.

Zoom In Zoom Out Reset image size

Table 3. List of 17 features extracted from the convolutional neural network. These features are mostly statistical descriptors of probability vectors of each class N (normal), A (atrial fibrillation) and O (other arrhythmia). These features are fed into a shallow neural network.

#	Feature name	Description
1	cnns(1)	mean(${{\overrightarrow{{\rm CNN}}}_{N}}$)
2	cnns(2)	mean(${{\overrightarrow{{\rm CNN}}}_{A}}$)
3	cnns(3)	mean(${{\overrightarrow{{\rm CNN}}}_{O}}$)
4	cnns(4)	std(${{\overrightarrow{{\rm CNN}}}_{N}}$)
5	cnns(5)	std(${{\overrightarrow{{\rm CNN}}}_{A}}$)
6	cnns(6)	std(${{\overrightarrow{{\rm CNN}}}_{O}}$)
7	cnns(7)	max(${{\overrightarrow{{\rm CNN}}}_{N}}$)
8	cnns(8)	max(${{\overrightarrow{{\rm CNN}}}_{A}}$)
9	cnns(9)	max(${{\overrightarrow{{\rm CNN}}}_{O}}$)
10	cnns(10)	min(${{\overrightarrow{{\rm CNN}}}_{N}}$)
11	cnns(11)	min(${{\overrightarrow{{\rm CNN}}}_{A}}$)
12	cnns(12)	min(${{\overrightarrow{{\rm CNN}}}_{O}}$)
13	avgNsubMaxO	mean(${{\overrightarrow{{\rm CNN}}}_{N}}$)  −  max(${{\overrightarrow{{\rm CNN}}}_{O}}$)
14	avgNsubMaxA	mean(${{\overrightarrow{{\rm CNN}}}_{N}}$)  −  max(${{\overrightarrow{{\rm CNN}}}_{A}}$)
15	maxCodeNAO	Index of class with maximal value (1, 2, 3 for N, A, O)
16	maxCodeAvgNMaxAO	Index of class with maximal value; class N uses its average instead of maximum
17	multNprev	avgNsubMaxO (#13)  ×  avgNsubMaxA (#14)

The computed features are fed into a shallow neural network which points directly to one of the three classes N, A and O. This shallow neural network contains 17 neurons in an input layer, 10 neurons in a hidden layer and three neurons in an output layer. The output layer implements a softmax function. The estimated test confusion matrix (70/30 dataset division for training/testing) is shown in table 4; the output for specific classes is shown in figure 6.

Zoom In Zoom Out Reset image size

Finally, output values from CNN/NN are compared to thresholds. If any of them passes, it means that there is high credibility of the CNN/NN output and the given result is used. Otherwise, the result from the bagged tree ensemble is used. Thresholds (0.8 for classes A and N, 0.75 for class O) were obtained experimentally with regard to the resultant F1 score (figure 7). In the training set, these thresholds were passed in 6636 of 8183 cases meaning that the CNN/NN result was preferred in 81% of cases.

Zoom In Zoom Out Reset image size