Coronary heart disease (CHD)is the leadingcause of deathin the world. This disease can be asymptomaticfor a long time and over time can progress and result in death. Today electrocardiogram (ECG) can be done at home with the help of special equipment from CardioQvark. In this paper the possibility of CHD detection based on such ECwGass explored. Di erent approaches to the classification of such electrocardiograms werseurveyed. New algorithms and modifications to existing algorithms were proposed. A new method - the ensemble of di erent algorithms - has shown the best performance.
Coronary heart disease(CHD) [
At the initial stages of the disease, most people do not show any symptoms of this disease. It is very important to identify the CHD in time to slow the course of the disease and prevent the patient's death.
Traditionally, CHD can be detected with the help of specialists and a number of tests. It should be noted that these tests take a signi cant amount of patient's time, and also require high quali cation of the specialist who will conduct them. Since there are very few specialists and the number of potential patients is growing every year, the task of automatically determining CHD is extremely urgent at the present time. Currently, it is extremely relevant to create a device that will help determine the disease or its probability at home. Such devices will allow the person to be sent to a doctor in case of a high probability of having a CHD.
An electrocardiogram (ECG) { is a signal that re ects the electrical activity of the heart. ECG
is one of the most a ordable ways of diagnosing heart disease now due to its non-invasiveness
and low cost. Currently, there are many di erent studies that show that the ECG can be used
to determine CHD [
CardioQvark (project site: www.cardioqvark.ru) has created a device in a form of a smartphone case that allows you to make ECG measurements at home. The CardioQVARK device is a portable electrocardiograph in the form of a smartphone case (iPhone 5 / 5s / SE / 6 / 6s), allowing to register data of bio-electrical activity of the heart from the rst ECG lead and next leads: aVR, aVL, aVF, Vi (i = 1 ... 6) using the patient's cable. In this work, the possibility of CHD detection based on such ECGs from the rst lead was explored.
There are di erent approaches to ECG classi cation problem. Some of them were surveyed in this work. New algorithms and modi cations to existing algorithms were proposed. In order to improve the performance of the classi cation, an ensemble of 5 di erent methods was built. Each of these methods will be described below.
All research was conducted on the basis of the following medical centers: NGHCI "Semashko Central Clinical Hospital 2", Federal State Scienti c Institution "Petrovsky Russian Scienti c Center of Surgery", Federal State-Funded Health Care Institution "City Clinical Hospital 4 Health Care Moscow Department", Federal State-Funded Health Care Institution "Moscow Clinical Scienti c Center of Moscow Department Health Care", State Autonomous Health Care Institution of the Moscow Region "Clinical Center for Restorative Medicine and Rehabilitation".
A voluntary anonymous study included patients over 18 years of age. Annotated impersonal electrocardiograms (ECG) were recorded from the rst ECG-lead using a CardioQVARK cardio monitor. The duration of each recording was 5 minutes. The measurement was taken in the sitting position, with support for the back, hands on the knees or on the table. Data was collected in dynamics of 3-10 observations with an interval of at least 12 hours between measurements.
The sample that was used for this task consists of 1798 cardiograms. It contained 1055 cardiograms of healthy patients and 743 cardiograms of patients with CHD. The sampling frequency was 1000 Hz.
Signals were preprocessed before applying machine learning algorithms. For preprocessing, a
low-pass and high-pass Butterworth lters [
Below are descriptions of the algorithms that were used to build the ensemble.
The idea which was used for this algorithm was described in the article [
R-peaks are computed.
The intervals between two R-peaks (RR-intervals) are measured.
Each value of the RR-interval is converted to 60=RR.
The main idea of this method is a construction of various groups of features from HRV signal.
The rst group of features includes various entropic features, which indicate a measure of unpredictability in the signal. The following types of entropies are used: approximate entropy, sample entropy, and Shannon entropy. Each type of entropy is described in detail below.
Approximate entropy is calculated as follows. Here x = (x0; x1; :::; xN 1) is the HRV signal of length N .
The integer m and the real r are xed.
The values Cim(r) are calculated as follows: A set of vectors of the form xim = (xi; xi+1; :::; xi+m 1), where i 2 [0; N m] is composed.
Cim(r) = xkm : d xim; xkm
N
r; k 2 [0; N m + 1 m] ; where where the xim(a) component a of the vector xim.
The values m(r) are de ned as:
d xim; xkm = a2m[0;amx 1] xim(a)
xkm(a) ; m(r) = (N
M + 1) l
log(Cim(r)): N m+1
X i=1 ApEn(r) = m(r) m+1(r): In this paper, the approximate entropy was realized for m = 10 and r = 0:2std(x), std(x) standard deviation of the signal x.
The Sample entropy is calculated as follows.
Vectors xim of length m and vectors xim+1 of length m + 1 are formed similar to those that were formed in approximate entropy.
The values A and B are calculated as follows:
A(r) = (xkm+1; xlm+1) : d(xkm+1; xlm+1) r; 0 k l
N m 1 ; B(r) = (xkm; xlm) : d(xkm; xlm) r; 0 k l
N m ; where d is determined as in approximate entropy.
Sample entropy (SampEn) is de ned as
SampEn (r) = log
A(r) B(r)
: ShanEn = k X pf log pf ; f=1 In this paper, the sample entropy was determined for m = 10, r = 0:2std(x).
The Shannon entropy is calculated as where k is the number of di erent elements in the signal x, pf is the frequency of the element f in the signal x.
The following group of features is based on a recurrence plot. This plot shows the frequency and duration of repetitions in the signal. The element (i; j) of the given plot is de ned as 1, if jjxi xj jj < ", or as 0 otherwise, where x is the HRV signal.
Based on this plot, the following features are calculated, N the number of elements in the HRV signal, ld min the minimum length of the plot diagonal, lv min minimal length of the vertical line, ld max maximum length of the diagonal, lv max maximum length of the vertical line:
Density of points (REC)
REC =
The percentage of points that form the diagonal lines (DET)
The average length of the diagonals (Lmean) Entropy of diagonal lines (ENd) where pl is the frequency of diagonal lines of length l.
Entropy of vertical lines (ENv)
DET =
Plld=mldamxin lP (l)
PiN;j R(i; j)
; Lmean =
Plld=mldamxin lP (l) Pl=ld min P (l)
: ENd =
pl log pl; ld max
X l=ld min lv max
X l=lv min ENv = pvl log pvl; where pvl is the frequency of vertical lines of length l.
For the calculation of these features, the PyRQA library [
Another group of features that is used for this approach is a group of features based on the Poincare plot. The Poincare plot is constructed as follows: for the signal x = (x0; x1; :::; xN 1), the plot consists of the points
(x0; x1); (x1; x2); :::; (xi; xi+1) and so on. In this case, RR-intervals were used as a signal x. constructed:
The standard deviation of distances from the points of the plot to the line y = x. This feature describes the local variability of RR-intervals.
The standard deviation of distances from the points of the plot to the line y = x 2RRmean. RRmean is the average value of RR-intervals. This feature describes the long-term variability of RR-intervals.
The feature, which is based on the detrended uctuation analysis. This method allows determining the self-dependence of the signal. The following cumulative sum is de ned as: t xcumsum(t) = X(x(i) i=1 ); x is a signal consisting of RR-intervals, is the mean of the signal x.
The data is segmented with a window of size n. On each segment polynomial is found for the data, which most accurately represents it (usually linear). The union of all such polynomials forms a function x n(t), which is an approximation of the original function xcumsum(t). Then there is the following function: v
N F ( n) = tuu N1 Xt=1[xcumsum(t) x n(t)]2; where N is the length of the signal consisting of RR-intervals.
The feature is the slope of the line log F ( n) to log( n). More information about this
approach is written in [
The next feature that was used in this approach is the correlation of dimensions. This feature is a quantitative characteristic of the signal trajectory and is de ned as follows.
Vectors xim of length m are similar to those that were formed in approximate entropy are constructed.
The value g is de ned as: g(r) = (xkm; xlm) : d(xlm; xkm) r; 0 k l
N m 1 ; g is the number of pairs of vectors whose distance is less than or equal to r. The value C(r) is de ned as where N is the length of the signal x.
The correlation of dimensions (D2) is de ned as
C(r) = g(r) N 2
; D2 = lim r!0 log(r)
: log C(r)
In this paper, the implementation of this feature was computed using publicly available software Nonlinear measures for dynamical systems or nolds, version 0.3.2, which can be downloaded from (https://pypi.python.org/pypi/nolds).
The gradient boosting from the xgboost package [
This algorithm is a mixture of 3 di erent feature spaces, which were previously used in the classi cation of biomedical signals.
The rst group of features consists of the parameters of the Hjorth's parameters:
activity, mobility, complexity [
The second group of features consists of statistical signal features: mean, standard deviation, signal minimum, signal maximum, skew, kurtosis, selective quantiles of order: 0.1, 0.25, 0.5, 0.75, 0.9, sums and sums of signal values squares that are above / below certain values of quantiles: 0.1, 0.25, 0.5, 0.75, 0.9.
The next group of features was suggested by Uspenskiy for the disease detection by patient's
ECG [
A(n) T (n)
It is assumed that the values of A(n), T (n) are not important, but signs of their increments are. The method of signal encoding based on all possible signs of increments of these quantities is presented in Table 1.
After the code representation of the signal is received, the three-gram selection is performed. The feature space is the number of occurrences of each of the possible three-grams in a given code sequence derived from the signal.
The logistic regression from the scikit-learn package [
The idea used in this algorithm was described in the article [
The neighborhoods of the signal R-peaks are allocated: 200 points before R-peak and 500 after.
The averaged neighborhood is used as a feature space.
Neural network with the architecture described in Table 2 was used as the classi cation model for this algorithm.
In this work, the neural network was implemented using the libraries Theano [
The idea for this algorithm was described in paper [
Wavelet signal transformation is a convolution of the signal with functions (t), called wavelets. Such wavelet functions should possess speci c properties:
Z +1 (t)dt = 0
Z +1
j (t)j2dt < 1: 1 1
Wavelet transformation allows achieving signal compression, with good performance of
reproduced original signal [
All wavelet functions used in a wavelet transformation can be represented with prototype function (t) using scaling and shift. In case of discreet wavelet transformation all wavelet function can be written as:
In discrete wavelet transformation of m;n function can be separated into two parts. The
two types of functions correspond to coe cients of approximation and detail coe cients. In
the present work approximation coe cients are used for new representation of the signal, and
Daubechies wavelets [
After getting the coe cients of wavelet transformations, local segments are extracted from the signal. The segments are extracted by moving a window of certain length w with step s, with all elements within one window going into a separate segment. After such a procedure every signal is represented as a set of local segments.
All local segments in the training set are separated into k clusters using k-means algorithm. After this, every segment is replaced with the number of the cluster it belongs to. This way every signal is represented as a text of codewords, with each word representing the certain cluster.
In the implementation of described algorithm these parameter values were used: w = 100, s = 30, k = 200. K-means algorithm implementation from scikit-learn package was used.
In the train sample every local segment is replaced with the cluster it is the closest to. That means that for every local segment si in test sample cluster c it belongs to is determined using this formula: c = argmin d(bj ; si); d(bj ; si) = tuvu Xw(sik
bjk)2; k=1 where bj is a cluster center j, sik (bjk) is a kth element of local segment i (cluster j).
Using this approach and transforming the input signal into text it is possible to use natural language processing algorithms. The paper's authors that suggested such encoding used the bag of words as a feature space. Feature description is a number of occurrences of each code word in a speci c signal.
Features based on word2vec technology were used in order to extract dependencies between
local segments in the signal. This approach was suggested by Google and it allows to use
context-aware text processing, reducing the dimensions of the data [
Word2Vec model was trained using the length of embedding vector equals to 80. A mean of
all vectors in the signal was used as a feature. The model was trained using package gensim [
Bispectrum is a function of two variables f1 and f2 that specify the frequencies, expressed by
the following formula [
B(f1; f2) = X(f1)X(f2)X (f1 + f2);
where X(f ) is the Fourier transform of the signal, and X (f ) is the complex conjugate of
it. The signal bispectrum is usually calculated using a fast Fourier transform. A detailed
description of the algorithm for the bispectrum computation of the signal can be found in [
Based on the resulting matrix, the elements of which can be denoted as a(i; j), we can associate each signal with a certain image. This image is calculated as follows. A new matrix B = jjbi;j jj is calculated, the elements of which are equal to:
q b(i; j) =
Re2a(i; j) + Im2a(i; j); where Re denotes the real part of the complex number, and Im denotes the imaginary part of the complex number. The contour plot of matrix B is used as the image.
The authors of the article [
It was proposed to use neural networks for classi cation of such images. The architecture of the neural network is described in Table 3.
In this work, the neural network was implemented using the libraries Theano [
In order to increase classi cation performance, it was suggested to use ensembles of the algorithms. Several methods of ensembling are described in this section.
Given a set of algorithms A = (A1; A2; :::; An) and output a vector of predictions a = (a1; a2; :::; an) then the resulting answer a of the ensemble is equal to
a = mode(a1; a2; :::; an); where mode is a statistic, that is equal to the element which is most often encountered in the predictions. If there are several of them a random one of them is chosen.
The main idea of the EM-algorithm [
N size of available data, nikl whether the k-th algorithm has the answer l (l 2 f1; 2g) to the data i (i = 1:::N ),
jkl (j 2 f1; 2g) is the probability that the k th algorithm will return j when the true answer is l.
Tij = 1, if the true answer for the data i is j, otherwise, it equals 0. pj the probability of class j in the sample.
Step 1: Initialize the matrices majority.
Step 2: Recalculate the values of matrices
and pj : to the ideal case. T initialize the voting value for the jkl =
P Tij nikl i P P Tij nikl l i
P Tij pj = i
N Tij =
pj QkK=1 Ql2=1( jkl)nikl
Pq2=1 pq QkK=1 Ql2=1( qkl)nikl
Step 3. Recalculate Tij : Repeat steps 2 and 3 until the matrices stop changing.
At the end of this algorithm, we obtain the probability of the data belonging to each class in the matrix T . As an answer, the class is taken, the probability of belonging to which is the greatest.
Cross-validation was used to evaluate the performance of algorithm. To avoid over tting the ECGs of one patient did not fall simultaneously into the training and test set. The following performance criteria was introduced:
N i=1
N Pjn=i1 Itij=pij ;
1 X
ni
where tij is the true value of the target variable for the cardiogram j of the patient i, pij
the predicted value of the target variable for the cardiogram j of patient i, ni the number
of cardiograms of the patient i, N the number of patients, Itij=pij is the indicator function
which equals to 1 if tij equals to pij and equals to 0 otherwise. This criterion is called patient
performance. It allows us to evaluate how well the algorithm determines a person's disease by
any of his cardiograms. In addition, this performance criterion does not depend on the number
of cardiograms for each patient. ROC-AUC score [
The results of evaluations are shown in Table 4, where the rst column shows algorithms type or ensemble type. The algorithm based on wavelet transformation is included in two variants, with word2vec and without.
In the course of this paper, the following results were obtained. When CardioQvark equipment is used, it is possible to determine CHD with an accuracy of more than 0.81 for patient performance, with an accuracy greater than 0.77 for F-score and with an accuracy greater than 0.87 for ROC-AUC score. Word2vec can increase the performance of the classi cation method based on the wavelet transformation. Bispectrum can be used to classify CHD. The EM algorithm is applicable for ensemble and in this case, shows the best performance of classi cation for all selected performance criteria.
Bispectrum Wavelet transformation Wavelet transformation +word2vec R-peak's neighborhood The HRV signal 3 di erent feature spaces majority EM