The work presented here applies deep learning to the task of automated cardiac auscultation, i.e. recognizing abnormalities in heart sounds. We describe an automated heart sound classification algorithm that combines the use of time-frequency heat map representations with a deep convolutional neural network (CNN). Given the cost-sensitive nature of misclassification, our CNN architecture is trained using a modified loss function that directly optimizes the trade-off between sensitivity and specificity. We evaluated our algorithm at the 2016 PhysioNet Computing in Cardiology challenge where the objective was to accurately classify normal and abnormal heart sounds from single, short, potentially noisy recordings. Our entry to the challenge achieved a final specificity of 0.95, sensitivity of 0.73 and overall score of 0.84. We achieved the greatest specificity score out of all challenge entries and, using just a single CNN, our algorithm differed in overall score by only 0.02 compared to the top place finisher, which used an ensemble approach.
Advances in deep learning [LeCun et al., 2015] are being made at a rapid pace, in part due to challenges such as ILSVRC – the ImageNet Large-Scale Visual Recognition Challenge [Russakovsky et al., 2015]. Successive improvements in deep neural network architectures have resulted in computer vision systems that are better able to recognize and classify objects in images [Lin et al., 2013; Szegedy et al., 2015] and winning ILSVRC entries [Szegedy et al., 2014; He et al., 2015]. While a large focus of deep learning has been on automated analysis of image and text data, advances are also increasingly being seen in areas that require processing other input modalities. One such area is the medical domain where inputs into a deep learning system could be physiologic time series data. An increasing number of large scale challenges in the medical domain, such as [Kaggle, 2014] and [Kaggle, 2015] have also resulted in improvements to deep learning architectures [Liang and Hu, 2015].
PhysioNet [Goldberger et al., 2000] has held a Computing in Cardiology Challenge since 2000 that requires participants to automatically analyze physiologic time series data. The 2016 challenge [Clifford et al., 2016] asked participants to perform automated analysis of phonocardiogram (PCG) waveforms, i.e. heart sound data collected using digital stethoscopes. The objective of the challenge was to accurately classify normal and abnormal heart sounds. Recordings were collected from both healthy individuals, as well as those with heart disease, including heart valve disease and coronary artery disease. A PCG plot showing the recording of the (normal) sounds made by the heart is given in Figure 1.
Heart disease remains the leading cause of death globally, resulting in more people dying every year due to cardiovascular disease compared to any other cause of death [World Health Organization, 2017]. Successful automated PCG analysis can serve as a useful diagnostic tool to help determine whether an individual should be referred on for expert diagnosis, particularly in areas where access to clinicians and medical care is limited.
In this work, we present an algorithm that accepts PCG waveforms as input and uses a deep convolutional neural network architecture to classify inputs as either normal or abnormal using the following steps:
semi-Markov model is used to segment incoming heart sound instances into shorter segments beginning at the start of each heartbeat, i.e. the S1 heart sound.
Mel-frequency cepstral coefficients, one dimensional time series input segments are converted into twodimensional spectrograms (heat maps) that capture the time-frequency distribution of signal energy.
A convolutional neural network is trained to perform automatic feature extraction and distinguish between normal and abnormal heat maps.
The contributions of this work are as follows: 1. We introduce a deep convolutional neural network architecture designed to automatically analyze physiologic time series data for the purposes of identifying abnormalities in heart sounds. 2. Given the cost-sensitive nature of misclassification, we describe a novel loss function used to train the above network that directly optimizes the sensitivity and specificity trade-off. 3. We present results from the 2016 PhysioNet Computing in Cardiology Challenge where we evaluated our algorithm and achieved a Top 10 place finish out of 48 teams who submitted a total of 348 entries.
The remainder of this paper is organized as follows. In Section 2, we discuss related works, including historical approaches to automated heart sound analysis and deep learning approaches that process physiologic time series input data. Section 3 introduces our approach and details each step of the algorithm. Section 4 further describes the modified costsensitive loss function used to trade-off the sensitivity and specificity of the network’s predictions, followed by Section 5, which details the network training decisions and parameters. Section 6 presents results from the 2016 PhysioNet Computing in Cardiology Challenge and in Section 7 we provide a final discussion and end with conclusions in Section 8. 2
Before the 2016 PhysioNet Computing in Cardiology Challenge there were no existing approaches (to the authors’ knowledge) that applied the tools and techniques of “deep learning” to the automated analysis of heart sounds [Liu et al., 2016]. Previous approaches relied upon a combination of feature extraction routines input into classic supervised machine learning classifiers. Features extracted from heart cycles in the time and frequency domains, as well as wavelet features, time-frequency and complexity-based features were input into artificial neural networks [De Vos and Blanckenberg, 2007; Ug˘uz, 2012a; 2012b; Sepehri et al., 2008; Bhatikar et al., 2005] and support vector machines [Maglogiannis et al., 2009; Ari et al., 2010; Zheng et al., 2015] for classification. Previous works have also employed Hidden Markov Models for both segmenting PCG signals into the fundamental heart sounds [Springer et al., 2014; 2016], as well as classifying normal and abnormal instances [Wang et al., 2007; Sarac¸oglu, 2012].
While there have been many previous efforts applied to automated heart sound analysis, gauging the success of historical approaches has been somewhat difficult, due to differences in dataset quality, number of recordings available for training and testing algorithms, recorded signal lengths and the environment in which data was collected (e.g. clinical vs. non-clinical settings). Moreover, some existing works have not performed appropriate train-test data splits and have reported results on training or validation data, which is highly likely to produce optimistic results due to overfitting [Liu et al., 2016]. In this work, we report results from the 2016 PhysioNet Computing in Cardiology Challenge, which evaluated entries on a large hidden test-set that was not made publicly available. To reduce overfitting, no recordings from the same subject were included in both the training and the test set and a variety of both clean and noisy PCG recordings, which exhibited very poor signal quality, were included to encourage the development of accurate and robust algorithms.
The work presented in this paper, is one of the first attempts at applying deep learning to the task of heart sound data analysis. However, there have been recent efforts to apply deep learning approaches to other types of physiological time series analysis tasks. An early work that applied deep learning to the domain of psychophysiology is described in [Mart´ınez et al., 2013]. They advocate the use of preference deep learning for recognizing affect from physiological inputs such as skin conductance and blood volume pulse within a game-based user study. The authors argue against the use of manual ad-hoc feature extraction and selection in affective modeling, as this limits the creativity of attribute design to the researcher. One difference between the work of [Mart´ınez et al., 2013] and ours is that they perform an initial unsupervised pre-training step using stacked convolutional denoising auto-encoders, whereas our network does not require this step and is instead trained in a supervised fashion end-to-end.
Similar deep learning efforts that process physiologic time series have also been applied to the problems of epileptic seizure prediction [Mirowski et al., 2008] and human activity recognition [Hammerla et al., 2016]. 3
Recall from Section 1 that our approach consists of three general steps: segmentation, transformation and classification. Each is described in detail below. 3.1
The main goal of segmentation is to ensure that incoming time series inputs are appropriately aligned before attempting to perform classification. We first segment the incoming heart sound instances into shorter segments and locate the beginning of each heartbeat, i.e. the S1 heart sound. A logistic regression hidden semi-Markov model [Springer et al., 2016] is used to predict the most likely sequence of heart sound states (S1 ! Systole ! S2 ! Diastole) by incorporating information about expected state durations.
Once the S1 heart sound has been identified, a time segment of length, T , is extracted. Segment extraction can either be overlapping or non-overlapping. Our final model used a segment length of, T = 3 seconds, and we chose to use overlapping segments as this led to performance improvements during initial training and validation. 3.2
Each segment is transformed from a one-dimensional time series signal into a two-dimensional heat map that captures the time-frequency distribution of signal energy. We chose to use Mel Frequency Cepstral Coefficents [Davis and Mermelstein, 1980] to perform this transformation, as MFCCs capture features from audio data that more closely resembles how human beings perceive loudness and pitch. MFCCs are commonly used as a feature type in automatic speech recognition [Godino-Llorente and Gomez-Vilda, 2004].
We apply the following steps to achieve the transformation: 1. Given an input segment of length, T , and sampling rate, ⌫ , select a window length, `, and step size, , and extract overlapping sliding windows, si(n), from the input time series segment, where i 2 [1, ⌅ T ⇧] is the window index and n 2 [1, `⌫ ] is the sample index. We chose a window length of 0.025 seconds and a step size of 0.01 seconds. 2. Compute the Discrete Fourier transform for each window. (1) (2) log transformation as sound volume is not perceived on a linear scale.
K c⇤i,j = log(X dj,kPi(k)) k=1 (3) We used a filterbank consisting of J = 26 filters, where frequency ranges were derived using the Mel scale that maps actual measured frequencies, f , to values that better match how humans perceive pitch, M (f ) = 1125 ln(1 + 7f00 ). 4. Finally, apply a Discrete Cosine Transform to decorrelate the log filterbank energies, which are correlated due to overlapping windows in the Mel filterbank.
ci,j =
J j=1 X c⇤i,j cos h k(2i
1)⇡ i 2J , k = 1 . . . J (4) The result is a collection of cepstral coefficients, ci,j for window, i. For i = 1 . . . ⌅ T ⇧, ci,j can be stacked together to give a time-frequency heat map that captures changes in signal energy over heart sound segments. Figure 2 illustrates two example heat maps (one derived from a normal heart sound input and the other from an abnormal input), where ci,j is the MFCC value (represented by color) at location, i, on the horizontal axis and, j, on the (inverted) vertical axis. 3.3
The result of transforming the original one-dimensional timeseries into a two-dimensional time-frequency representation is that now each heart sound segment can be processed as an image, where energy values over time can be visualized as a heat map (see Figure 2). Convolutional neural networks are a natural choice for training image classifiers, given their ability to automatically learn appropriate convolutional filters. Therefore, we chose to train a convolutional neural network architecture using heat maps as inputs.
Decisions about the number of filters to apply and their sizes, as well as how many layers and their types to include in the network were made by a combination of initial manual Si(k) = `⌫ X si(n)h(n)ei2⇡n `k⌫ n=1 where k 2 [1, K], K is the length of the DFT and h(n) is a hamming window of length N . The power spectral estimate for window, i, is then given by (2).
1
Pi(k) = N |Si(k)|2 3. Apply a filterbank of, j 2 [1, J ], triangular band-pass filters, dj,1...K , to the power spectral estimates, Pi(k), and sum the energies in each filter together. Include a exploration by the authors, followed by employing a random search over a limited range of network architecture parameters. Figure 3 depicts the network architecture of a convolutional neural network that accepts as input a single channel 6x300 MFCC heat map and outputs a binary classification, predicting whether the input segment represents a normal or abnormal heart sound.
The first convolutional layer learns 64 2x20 kernels, using same-padding. This is followed by applying a 1x20 maxpooling filter, using a horizontal stride of 5, which has the effect of reducing each of the 64 feature maps to a dimension of 6x60. A second convolutional layer applies 64 2x10 kernels over the previous layer, once again using same padding. This is again followed by a max-pooling operation using a filter size of 1x4 and a stride of 2, further reducing each feature map to a dimension of 6x30. At this stage in the architecture a flattening operation is applied that unrolls each of the 64 6x30 feature maps into a single dimensional vector of size 11,520. This feature vector is fed into a first fully connected layer consisting of 1024 hidden units, followed by a second layer of 512 hidden units and finally a binary classification output. 4
The loss function of the network was altered from a standard softmax cross entropy loss function to instead directly tradeoff between sensitivity and specificity.
Given unnormalized log-probabilities, y = W x + b, from a classifier consisting of weight matrix, W , and bias b. The softmax function: s(y(1)) s(y(2)) . .
. s(y(n)) 3 7 6 7 , Y ⇤ = 6 7 6 5 4 2 y⇤ (1) y⇤ (2) . .
. y⇤ (n) 3 where s(yi(j)), refers to the ith entry of row j and Y ⇤ is the corresponding one hot encoded matrix of actual class labels.
For the binary class labels of normal (y0⇤ ) and abnormal (y1⇤ ), we define the mask matrices, YNn and YAa, where entries within each matrix are softmax prediction values extracted 8 s(y(j))2 Y , as follows: (5)
6 Y = 6 6 4 2 YNn = YAa = 8 s(y0(j)), ><>> >>>: 0, 8 s(y1(j)), >>>< >:>> 0, where y0⇤(j) = 1 and arg max {s(y(j))} = arg max{y⇤ (j)} otherwise otherwise where y1⇤(j) = 1 and arg max {s(y(j))} = arg max{y⇤ (j)}
We then define softmax sensitivity, Se, and specificity, Sp, as follows:
Se = X YA(ja) , j Y A⇤a(j)
Sp = X YN(jn) j Y N⇤(nj) The final loss function we wish to minimize is given in (7).
LSeSp =
(Se + Sp) + R (W ) where R (W ) is a regularization parameter and routine, respectively. (6) (7) s(yi) =
eyi P eyj j gives probability predictions P (yi|x; W, b) for the class at index, i, for input x.
Consider, L2 256 Adam Optimization L2 regularization was computed for each of the fully connected layers’ weight and bias matrices and applied to the loss function. Dropout was applied within both fully connected layers. Table 1 shows the values of hyper-parameters chosen by performing a random search through parameter space, as well as a list of other network training choices, including weight updates and use of regularization. Adam optimization [Kingma and Ba, 2014] was used to perform weight updates. Models were trained on a single NVIDIA GPU with between 4 – 6 GB of memory. A mini-batch size of 256 was selected to satisfy the memory constraints of the GPU. The overall dataset used within the PhysioNet Computing in Cardiology Challenge was provided by the challenge organizers and consisted of eight heart sound databases collected from seven countries over a period of more than a decade [Clifford et al., 2016]. In total 4,430 recordings were taken from 1,072 subjects, resulting in 30 hours of heart sound recordings. From this total dataset, 1,277 heart sound recordings from 308 subjects were removed to be used as heldout test data for evaluating challenge submissions. The test dataset was not made publicly available and challengers were only allowed to make 15 submissions, in total, to the challenge server to evaluate their models on a small 20% subset of the hidden dataset, before final results were computed. The number of allowed submissions was limited to avoid the issue of participants implicitly overfitting their models on the hidden test dataset.
From the 3153 publicly available PCG waveforms supplied by the challenge organizers, the authors set aside a further 301 instances to be used as a local held-out test-set to gauge model performance before making a submission to the challenge server. The remaining instances were used to train initial models. Models were trained on the overlapping 3-second MFCC segments extracted from the remaining 2852 PCG waveforms. This resulted in approximately 90,000 MFCC heat maps, which were split into a training (⇠ 75, 000 instances) and validation set (⇠ 15, 000 instances). This training and validation set was unbalanced, consisting of approximately 80% normal segments and 20% abnormal segments. Training was performed on the unbalanced dataset and no attempt was made to compensate for this class imbalance.
Given that each model was trained on 3-second MFCC heat map segments, it was necessary to stitch together a collection of predictions to classify a single full instance. The simple strategy of averaging each class’s prediction probability was employed and the class with the greatest probability was selected as the final prediction. 6
Equations (8) and (9) show the modified sensitivity and specificity scoring metrics that were used to assess the submitted entries to the 2016 PhysioNet Computing in Cardiology Challenge [Clifford et al., 2016]. Uppercase symbols reflect the true class label, which could either be (A)bnormal, or (N )ormal. Lowercase symbols refer to a classifier’s predicted output where, once again, a is abnormal, n is normal and q is a prediction of unsure. A subscript of 1 (e.g. Aa1, N a1) refers to heart sound instances that were considered good signal quality by the challenge organizers and a subscript of 2 (e.g. An2, N n2) refers to heart sound instances that were considered poor signal quality by challenge organizers. Finally, the weights used to calculate sensitivity, wa1 and wa2, capture the percentages of good signal quality and poor signal quality recordings in all abnormal recordings. Correspondingly for specificity, the weights wn1 and wn2 are the proportion of good signal quality and poor signal quality recordings in all normal recordings. Overall, scores are given by Se+Sp . 2 Se = Sp =
wa1 · Aa1 Aa1 + Aq1 + An1
wn1 · N n1 N a1 + N q1 + N n1 + + wa2 · (Aa2 + Aq2) , Aa2 + Aq2 + An2 wn2 · (N n2 + N q2) N a2 + N q2 + N n2 (8) (9)
Table 2 shows a selected subset of the results for the 2016 PhysioNet Computing in Cardiology Challenge. For each selected entry, sensitivity, specificity and overall scores are shown, as well as the entry’s final ranking and a brief description of its approach. In total, 348 entries were submitted by 48 teams. Our entry, as described by the algorithm presented in this paper, was ranked 8th with a sensitivity of 0.7278 and specificity of 0.9521, giving an overall score of 0.8399. The top entry to the competition achieved sensitivity of 0.9424, specificity of 0.7781 for an overall score of 0.8602. Also included in Table 2 is the result of a benchmark entry that was supplied by the challenge organizers, which ranked 43rd overall, with a sensitivity of 0.6545 and specificity of 0.7569, for an overall score of 0.7057. 7
Table 2 shows that the overall scores for the top entries to the PhysioNet Computing in Cardiology challenge were very close. In particular, our entry, which achieved an 8th place ranking, had a difference in score of only 0.02, compared to the top place finisher. For our entry, the overall score of 0.8399 was achieved using a single convolutional neural network, whereas other top place finishers achieved strong classification accuracies using an ensemble of classifiers. Improvements in performance have often been witnessed using an ensemble of networks or separate classifiers and we leave this for future work/improvement. For practical purposes, a diagnostic tool that relies on only a single network, as opposed to a large ensemble, has the advantage of limiting the amount of computational resources required for classification. Deployment of such a diagnostic tool on platforms that impose restricted computational budgets, e.g mobile-based, could perhaps benefit from such a trade-off between accuracy and computational cost.
Another point of interest is that our entry to the PhysioNet Computing in Cardiology challenge achieved the greatest specificity score (0.9521) out of all challenge entries. However, the network architecture produced a lower sensitivity score (0.7278). Once again, considering the practical result of deploying a diagnostic tool that relied upon our algorithm, this would likely result in a system with few false positives, but at the expense of misclassifying some abnormal instances. Final decisions about the trade-off between sensitivity and specificity would require further consideration of the exact conditions and context of the deployment environment.
A final point of discussion and area of future improvement is that the approach presented was limited to binary decision outputs, i.e. either normal or abnormal heart sounds. An architecture that also considered signal quality as an output would likely result in performance improvement. 8
Conclusion The work presented here is one of the first to apply deep convolutional neural networks to the task of automated heart sound classification for recognizing normal and abnormal heart sounds. We have presented a novel algorithm that combines a CNN architecture with MFCC heat maps that capture the time-frequency distribution of signal energy. The network was trained to automatically distinguish between normal and abnormal heat map inputs and it was designed to optimize a loss function that directly considers the trade-off between sensitivity and specificity. We evaluated the approach by submitting our algorithm as an entry to the 2016 PhysioNet Computing in Cardiology Challenge. The challenge required the creation of accurate and robust algorithms that could deal with heart sounds that exhibit very poor signal quality. Overall, our entry to the challenge achieved a Top-10 place finish out of 48 teams who submitted 348 entries. Moreover, using just a single CNN, our algorithm differed by a score of at most 0.02 compared to other top place finishers, all of which used an ensemble approach of some kind.