Ischemia episode detection in ECG using kernel density estimation, support vector machine and feature selection

Materials

We used the European ST-T database from Physionet. European ST-T database has 90 records which are two-channel and each two hours in duration [18, 19]. Each record in this database has a different number of ST episodes. Overall there are 367 ischemic ST episodes in the database. Sampling frequency of each ECG data is 250 Hz.

We excluded 5 records because these had some problems. The records e0133, e0155, e0509, and e0611 had no ischemic ST episodes. The record e0163 had so limited ST episode whose length was just 31 seconds.

Removing baseline wandering

The ST segments in ECG can be strongly affected by baseline wandering [20]. Main causes of the baseline wandering are respiration and electrode impedance change due to perspiration [20, 21]. The frequency content of the baseline wandering is usually in a range below 0.5 Hz [20, 21].

We use discrete wavelet transform to remove baseline wandering in ECG. We transform signal vector into two sequences of coefficients, approximation and detail coefficients sequences [22]. We do this in each step in an iterative fashion, until we get an input signal whose length is smaller than the length of the filter which characterizes the wavelet. In our case, we used Daubechies8 wavelet with filter length of 8. The resulting approximation coefficient sequence becomes the input signal to the next discrete wavelet transform as shown in Figure 3(a) [22].

In each step, the coefficient sequence implies a band of frequencies. If the sampling frequency of a discrete ECG signal ecg(n) is x, we can determine a continuous and band-limited signal within frequency limits of $\left[0,\frac{x}{2}\right]$ by Nyquist sampling theorem [23]. Therefore if we have transformed the input signal ecg(n) into the approximation coefficient sequence h ₁(n) and detail coefficient sequence g ₁(n), then the frequency content of g ₁(n) is from $\frac{x}{4}$ to $\frac{x}{2}$, and the frequency content of h ₁(n) is below $\frac{x}{4}$. In this regard, if we have transformed the ecg(n) into the approximation coefficient sequence h _k (n), and the detail coefficient sequences g _k (n),g _k−1(n),· · ·,g ₁(n), the frequency contents of g _k (n),g _k−1(n),· · ·,g ₁(n) become $\left[\frac{x}{2^{k+1}},\frac{x}{2^k}\right],\left[\frac{x}{2^k},\frac{x}{2^{k-1}}\right],·\; ·\; ·,\left[\frac{x}{2^2},\frac{x}{2}\right]$ respectively [24, 25].

To remove baseline wandering, we should choose appropriate wavelet scale. We follow argument similar to that presented by Arvinti et al. except that they used stationary wavelet transform instead of its discrete counterpart [26]. We remove signal components whose frequency content is less than 1/2 Hz [20, 21]. If we have transformed the ECG signal ecg(n) into coefficient sequences h _k (n),g _k (n),g _k−1(n),· · ·,g ₁(n), the frequency contents of h _k (n) and g _k (n) become $\left[0,\frac{x}{2^{k+1}}\right]$ and $\left[\frac{x}{2^{k+1}},\frac{x}{2^k}\right]$ respectively, where x is the sampling frequency. If we choose k as $\frac{x}{2^{k+1}}\le \frac{1}{2}$ $k=⌈\underset{2}{log}x⌉$, the frequency content of the approximation coefficient sequence h _k (n) becomes less than 1/2 Hz. Thus, we assign zero sequence $\overrightarrow{0}\left(n\right)$ to all the detail coefficient sequences g _k (n),g _k−1(n),· · ·,g ₁(n), and calculate inverse transform of $h_k\left(n\right),\overrightarrow{0}\left(n\right),\overrightarrow{0}\left(n\right),·\; ·\; ·,\overrightarrow{0}\left(n\right)$ to form the baseline(n) in the bottom of Figure 3(b). If we subtract baseline(n) from ecg(n), we obtain the flattened signal like the one shown in Figure 3(c).

If we select a wrong wavelet scale k to find coefficient sequences of ecg(n), we obtain disappointing results. The flattened signal in Figure 3(c) is obtained when k is $⌈\underset{2}{log}250⌉=8$, where 250 is the sampling frequency expressed in Hz. When select k=4 to use h ₄(n),g ₄(n),g ₃(n),g ₂(n),g ₁(n), we obtain a plot in Figure 3(d). The middle waveform, baseline(n), resulted from the inverse discrete wavelet transform of $h_4\left(n\right),\overrightarrow{0}\left(n\right),\overrightarrow{0}\left(n\right),\overrightarrow{0}\left(n\right)\overrightarrow{0}\left(n\right)$. This middle waveform is too detailed, so the bottom waveform ecg(n)-baseline(n) was negatively affected. When we select k=12, see Figure 3(e), the bottom waveform was not different from the input waveform ecg(n).

We adopt a discrete wavelet transform to retain the details of the ECG waveform because filtering by some cut-off frequency can deteriorate the quality of the ECG waveforms [27].

Detecting QRS complexes

We have to select an appropriate wavelet scale to capture proper time positions of QRS complexes. We will deal with only the flattened ECG waveform ecg(n)-baseline(n) referred in the previous section. We will denote it as fecg(n).

We select the wavelet scale j which produces the largest drop of scor e _j −scor e _{j + 1}(j≥2). The bottom waveform in Figure 4(b) shows the time positions of QRS complexes when selecting this suitable wavelet scale.

After finding the locations of QRS complexes, we choose QRS onset and offset points in each QRS complex. We search QRS onset point in backward direction from a peak point in each QRS complex. We take the QRS onset point if the point is at the place of changing direction of rising and falling of fecg(n) twice. In the same way, we take the QRS offset point in forward direction from the peak point.

Algorithm 1 shows a process of removing baseline wandering and detecting QRS complexes.

Feature formation for classification problems

We deal with the flattened waveform, fecg(n), to obtain the values of the features. We take voltage level of QRS onset point as the reference from which we measure voltage deviation [2, 28]. We denote the mean value of electric potentials at QRS onset points as $\overline{\mathit{\text{fecg}}\left(\mathit{\text{QRS}}\text{onset}\right)}$. We consider this value as an effective zero voltage, so we measure voltage deviation from the $\overline{\mathit{\text{fecg}}\left(\mathit{\text{QRS}}\text{onset}\right)}$.

To form the first feature, we sum up all the voltage deviation from QRS offset point to T wave peak point as shown in Figure 5(a) and (b).

${\mathit{\text{feature}}}_1=\sum _{i=\mathit{\text{QRS}}\text{offset}}^{T\text{peak}}\left|\mathit{\text{fecg}}\left(i\right)-\overline{\mathit{\text{fecg}}\left(\mathrm{QRS}\text{onset}\right)}\right|$

(2)

We calculate these three feature values for each heart beat. Then we average these values in five successive beats, and arrange the three mean values as $\left({\mathit{\text{feature}}}_1,{\mathit{\text{feature}}}_2,{\mathit{\text{feature}}}_3\right)$.

Algorithm 2 shows the pseudo-code of computing feature values.

Classification by kernel density estimation

We classify a test point by examining posterior probabilities in which the test point belongs to two classes, normal or ischemic ST episode. We assume we have n _S points $\left\{{\mathbf{x}}_1^{\left(S\right)},{\mathbf{x}}_2^{\left(S\right)},·\; ·\; ·,{\mathbf{x}}_{n_S}^{\left(S\right)}\right\}$, and n _N points $\left\{{\mathbf{x}}_1^{\left(N\right)},{\mathbf{x}}_2^{\left(N\right)},·\; ·\; ·,{\mathbf{x}}_{n_N}^{\left(N\right)}\right\}$. The first and the second set designate training sets of ischemic ST episode and normal part, respectively. Each point is described by three components $\left({\mathit{\text{feature}}}_1,{\mathit{\text{feature}}}_2,{\mathit{\text{feature}}}_3\right)$.

The quantities $b_i^{\left(N\right)}$ and $b_i^{\left(S\right)}$ $\left(1\le i\le 3\right)$ are called kernel bandwidths. We calculate these bandwidths for each class (N or S) and component $\left(1\le i\le 3\right)$. These kernel bandwidths impact accuracy of kernel density estimation [32].

We choose half of the mean, $\frac{1}{2}{\mathit{\text{mean}}}_i^{\left(\mathrm{cl}\right)}$, as kernel bandwidth $b_i^{\left(\mathrm{cl}\right)}$ for each class cl (N or S), and component i $\left(1\le i\le 3\right)$.

Classification with the use of support vector machine

where the matrix H is expressed as $H_{\mathrm{ij}}\equiv t_i^{\left(\mathrm{cl}\right)}{t}_j^{\left(\mathrm{cl}\right)}K\left({\mathbf{x}}_i^{\left(\mathrm{cl}\right)},{\mathbf{x}}_j^{\left(\mathrm{cl}\right)}\right)=t_i^{\left(\mathrm{cl}\right)}{t}_j^{\left(\mathrm{cl}\right)}\varphi \left({\mathbf{x}}_i^{\left(\mathrm{cl}\right)}\right)·\varphi \left({\mathbf{x}}_j^{\left(\mathrm{cl}\right)}\right)=t_i^{\left(\mathrm{cl}\right)}{t}_j^{\left(\mathrm{cl}\right)}{e}^{-\frac{1}{3}{∥{\mathbf{x}}_i^{\left(\mathrm{cl}\right)}-{\mathbf{x}}_j^{\left(\mathrm{cl}\right)}∥}^2}$ [33]. When we classify a new pattern y, we examine decision function, $\mathrm{sgn}\left(\sum _{j=1}^{n_{\mathrm{cl}}}{t}_j^{\left(\mathrm{cl}\right)}{\alpha }_{j}K\left({\mathbf{x}}_j^{\left(\mathrm{cl}\right)},\mathbf{y}\right)+b\right)$. Whenever the input training set $\left\{{\mathbf{x}}_1^{\left(\mathrm{cl}\right)},{\mathbf{x}}_2^{\left(\mathrm{cl}\right)},·\; ·\; ·,{\mathbf{x}}_{n_{\mathrm{cl}}}^{\left(\mathrm{cl}\right)}\right\}$ was changed, we varied the parameter C to find its value which produced the highest classification rate.

Experiments setting

We used kernel density estimation and support vector machine methods to evaluate the proposed approach. We completed the experiment for each channel and record available in the European ST-T database. First, we trained the classifier based on a subset of ST episodes and normal ECG. Then we tested how well the feature values discriminated the two classes, ST episode and normal. When we formed the ST episode data, we used all the ischemic ST episodes except ST deviations data resulted from non-ischemic causes such as position related changes in the electrical axis of the heart. To preserve balance between ST episode and normal ECG data, we collected normal data from the beginning of each record as much as the amount of ST episode data. When dividing the data into training and test sets, we assigned one tenth of data to the training data, and the rest to the test data. In the cases of e0106 lead 0, e0110, e0136, e0170, e0304, e0601, and e0615 records, we constructed the training data of one third of all data and test data of two thirds because these records had much small ischemic ST episode data. To avoid ambiguous region between ischemic ST episode and normal ECG, we removed 10 seconds amount of ECG data from each side of the boundary.

When we classified a test set $\left\{{\mathbf{y}}_i\right\}$, four quantities were computed: true positive (TP), false negative (FN), false positive (FP), and true negative (TN). TP is a number of ischemic events correctly detected. FN is a number of erroneously rejected (missed) ischemic events. FP is a number of non-ischemic, that is, normal parts which the classifier erroneously detected as ischemic events. TN is a number of normal parts which our classifier correctly rejected as non-ischemic events [34]. These are numbers of corresponding y _i points which were obtained by averaging three feature values of successive five beats in Algorithm 2. The sensitivity and specificity are expressed in a usual fashion, Se=TP/(TP + FN) and Sp=TN/(TN + FP) respectively [6].

We tested the classifiers by counting how many ST episodes were correctly caught, out of 367 episodes in the 85 records of European ST-T database. For an interval of ischemic ST episode data, we formed n test points $\left\{{\mathbf{y}}_1,{\mathbf{y}}_2,·\; ·\; ·,{\mathbf{y}}_n\right\}$ from the data (Algorithm 2), and classified each test point and then counted numbers of two classes, “ischemic” and “normal”. If the number of class “ischemic” was larger than n/2, we declared the interval to be an ischemic ST episode. The experiments were completed for 367 ischemic ST episodes.

We compared the results of kernel density estimation (KDE) and support vector machine (SVM) methods with those formed by artificial neural network (ANN). The corresponding ANN classifier exhibits the following topology. The input layer has three nodes which accept featur e ₁ featur e ₂ and featur e ₃respectively. The output layer has two nodes which have target values (1, 0) and (0, 1) in the cases of “ischemia” and “normal” classes, respectively. We initialized bias weights as 0, and assigned random values between -1.0 and 1.0 to the weights of the network. The learning was carried out by running the backpropagation method [17] for 3000 iterations. We used a sigmoid activation function $1/\left(1+e^{-x}\right)$ and set learning rate 0.01. We adopted various topologies of hidden layers such as $3\to \left(5\right)\to \left(5\right)\to 2$ $3\to \left(6\right)\to 2$ $3\to \left(7\right)\to 2$ and $3\to \left(8\right)\to 2$ where the number in each parenthesis represents a number of nodes in the corresponding hidden layer. We used stochastic (incremental) gradient descent method to alleviate some drawbacks of the standard gradient descent method, see [17].

KDE with various kernels

Results for KDE, SVM and ANN with various wavelets

We examined the classifiers to find out how their performance depends on the mother wavelets which were used to produce training and test sets in Algorithm 1. We used 7 wavelets, Haar, Daubechies4, Daubechies8, Daubechies10, Coiflet6, Coiflet12 and Coiflet18 [22, 36]. The number forming a part of the name of each wavelet designates the length of filter which characterizes corresponding wavelet. Figure 6 shows selected shapes of wavelet functions except for the Haar wavelet which is given as

$\mathit{\text{Haar}}\left(t\right)=\left\{\begin{array}{c}1\left(0\le t\le 1/2\right)\\ -1\left(1/2\le t\le 1\right)\end{array}\right..$

(20)

Tables 3 5 and 6 show that the Daubechies8 and Daubechies10 wavelets give us superior results. The shapes of these two wavelets are similar to typical ECG waveforms [37, 38]. From now on, we use the Daubechies8 wavelet exclusively.

Effects of baseline wandering in ECG

Effects of simulated noise

We examined performance of the classifiers when we added simulated noise into the original ECG signal. We modeled the noise as the sum of wandering baseline and AC power line 60 Hz noise.

where a is an amplification factor and b is an angular frequency of the added baseline. Here Samp _Freq means sampling frequency which was 250 Hz in our case. We varied a from 0.1 to 1.0 in step of 0.1, and selected b to be equal to 2, 4 or 6.

Figure 7 shows the original ECG and its noise-impacted version. Tables 9, 10 and 11 show the experimental results for the noisy ECG signal. The first, second and third column in each b item represent the sensitivity, specificity and the “detect” respectively. The kernel bandwidth is set as $b_i^{\left(\mathrm{cl}\right)}={\mathit{\text{mean}}}_i^{\left(\mathrm{cl}\right)}/2$ for the KDE classifier. The layer composition of the ANN classifier was $3\to \left(7\right)\to 2$. The first column in each b item in Table 10 includes the trade-off parameter C which produced best results.

Comparison with others’ works

We used the Daubechies8 wavelet in Algorithm 1 to analyze the ECG waveform, and took the kernel bandwidths $b_i^{\left(\mathrm{cl}\right)}={\mathit{\text{mean}}}_i^{\left(\mathrm{cl}\right)}/2$ for the KDE classifier with Gaussian kernel. We used the SVM classifier with C=281.1.