Detection and Classification of Abnormities of First Heart Sound Using Empirical Wavelet Transform

It is expected that an automatic detection and classification algorithm for the abnormities of first heart sound (S1) can realize computer artificial intelligence diagnosis of some relative cardiovascular disease. Few studies have focused on the detection and classification of the abnormities of S1 and given out in detail the essential differences between abnormal and normal S1. This work applied Empirical Wavelet Transform (EWT) to decompose S1 and extracted the instantaneous frequency (IF) of mitral component (M1) and tricuspid component (T1) by using Hilbert Transform. Firstly, the heart sound signal is preprocessed following these processes: filtering, resampling, normalization and segmentation. Secondly, S1 is decomposed into several modes based on EWT. First two maximal points with a distance greater than 20Hz in Fourier Spectrum of S1 are selected and the nearest minimal points on both sides of the maximal points are found out as the boundaries for segmentation of the spectrum. S1 is decomposed into 5 modes and every mode’s IF are calculated through Hilbert transformation. At last, a k-mean cluster algorithm is applied to cluster the IF of different modes. TD and  $A_{peak\_{}ratio}$  are calculated for decision tree classifier and S1s are divided into three categories: normal S1, S1 with abnormal split and S1with abnormal amplitude change. When the proposed method is applied to detect normal S1, Se=94.6%, Pp=98.6% and Oa=93.3%; When it is applied to detect S1 with abnormal split, Se=92.6%, Pp=96.9% and Oa=90%; When it is applied to detect S1 with abnormal amplitude change, Se=94.4%, Pp=95.7% and Oa=90.6%; Comparison experiments are carried out between the proposed method and HVD method. The results show Oa of the proposed method is higher than HVD method when detecting the three different S1s.


I. INTRODUCTION
Heart sounds reflect the most primitive mechanical movement of heart. With the advent of electronic stethoscope, the real-time storage and playback of phonocardiogram (PCG) are realized [1]. The analysis of heart sound signal can not only overcome the limitations of human hearing sensitivity and experience, but also discover more The associate editor coordinating the review of this manuscript and approving it for publication was Liangxiu Han . pathological information [2]. It is expected to be a noninvasive, low-cost and convenient diagnostic method through analyzing of heart sound signal [2]- [3].
So far, many time-frequency analytical methods were applied to analyze various types of heart sounds, but there are few methods to analyze the abnormities of S1. Abnormities of S1 usually include abnormal split and abnormal amplitude changes. The classical and generally accepted theory for the production of S1 is that it is associated with the closure of the mitral (M1) and tricuspid (T1) valves [4]. M1 is the first VOLUME 7, 2019 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ audible component of S1. This component normally occurs before T1 in time domain [5]. It has great significance to diagnose the abnormities of S1 because inactivity of calcified atrioventricular valves and conduction delay or systolic delay are associated with abnormities of S1. However, it is usually difficult to distinguish abnormal S1 from normal S1 and other easily confused additional heart sounds and murmurs, such as the atrial consonant before S1 and the highfrequency murmurs in the pre-contraction of ventricle [6]. Clinical judgment of abnormal S1 is completely subjective and comparative without clear objective criteria. It is expected that an automatic detection and classification algorithm of the abnormities of S1 can realize computer intelligence diagnosis of some relative cardiovascular disease automatically. Heart sound signal is non-stationary and non-linear signal [7]. At present, the time-frequency analysis method is generally used to analyze the heart sound signal, such as Short-time Fourier Transform (STFT), wavelet transform (WT), Wigner distribution (WD), empirical mode decomposition (EMD), etc. [8]. STFT cannot perfectly detect the two internals components of S1, WD also cannot accurately detect the components M1 and T1 and seems them to be one component [9]. Wavelet transform lacks adaptability, and the decomposition effect of signals mainly depends on the choice of wavelet basis [10]. The WD has a disadvantage because it shows cross-terms in its response [11]. For discrimination M1 and T1, decomposition of S1 is a major concern. In recent years, Huang et al. have put forward an Empirical Mode decomposition method (EMD) [12]- [13]. EMD is fully adaptive, but it lacks mathematical theoretical basis, and has problems of endpoint effect, over-envelope, under-envelope, false modes. Shovan Barma and Bo-Wei Chen decomposed S2 using Hilbert vibration decomposition (HVD) [9]. However, the mean relative error shows a very low value in contrast to the others in discriminating of S2 split. In 2013, French scholar Gilles puts forward the empirical wavelet transform (EWT) combining the advantages of EMD [14]. In this method, an orthogonal wavelet filter is constructed in each interval to extract the inherent modes with tight support of the Fourier spectrum though adaptive segmenting the signal's Fourier spectrum. The extracted modes are transformed into IF through Hilbert transformation. EWT algorithm is simple, visual and intuitive, which can avoid unnecessary and incorrect modes. In 2017, V. Nivitha Varghees has successfully applied the EWT algorithm to the precise segmentation and classification of various heart sounds and murmurs [15]. However, the premise of obtaining the correct IF is that each segmentation contains only one spectral peak point. If it contains two or more, the IF obtained by EWT algorithm is wrong. Therefore, how to apply EWT flexibly to divide spectrum and decompose signal determines whether the time frequency analysis results are correct or not. This paper proposed a new method to decompose S1 based on EWT although S1's Fourier spectrum is not simple and have too much number of peaks. In fact, S1 contains multiple frequency components, and M1 and T1 are the two main frequency components. We focus on extraction of the IF of M1 and T1 and pay less attention to other frequency components. Therefore, maximum points in the spectrum of S1 were sorted and the first two discontinuous maximum points with a distance greater than 20Hz were extracted. Then, the minimum points closest to these two maximum points were selected as the boundaries for segmentation of the spectrum of S1.
As a matter of fact, there is little research on detection of abnormities of S1. Most research were focused on S2 split. Meanwhile, there is no research on classification of abnormities of S1. This work not only detects S1 split but also detect abnormal amplitude change of S1. The proposed algorithm can automatically clearly show how the maximal IF of M1 and T1 varies with time and amplitude in a threedimensional scatter diagram. Since the moment corresponds to the maximum of the IF of M1 is the moment that the mitral valve closes, and the moment corresponds to the maximum of the IF of T1 is the moment that the tricuspid valve closes, the difference between the two moments (TD) can be used as an objective criterion to evaluate S1 split. As is known to all, the criterion to evaluate S1 split is TD >0.04ms. The ratio between the maximal amplitude of IF of T1 to the maximal amplitude of the IF of M1 (A peak_ratio ) will be used as an objective criterion to evaluate S1 with abnormal amplitude change because S1 with abnormal amplitude changes are always accompanied by abnormal amplitude changes in the IF of M1 or T1.However, from now on, there is few criterions for detection of S1 with abnormal amplitude change. This work intends to use the statistical method to give out the exact objective criterion for evaluation of S1 with abnormal amplitude change through a large number of analysis of the IF's amplitude of M1 and T1 of different S1. Finally, the performance of the proposed algorithm was evaluated using another database including 150 normal S1s, 70 S1s with abnormal split and 75 S1s with abnormal amplitude change.
This paper is organized as follows: Section II introduces the overview of the proposed method. Section III then describe details of the proposed method. In section IV, the proposed method was applied to analyze a small database and the key characteristic parameters (TD and A peak_ratio ) are extracted for detection and classification of abnormities of S1. Finally, a bigger database was used to evaluate the performance of the proposed method. At the same time, comparison experiments were carried out between HVD method and the proposed method. Section V is the conclusion of this study.

II. SYSTTEM OVERVIEW
Some of the heart sound data used in this paper are collected from website and others are collected by the MEMS electronic stethoscope made in our laboratory and 3M littmann 3200 electronic stethoscope [1]. The proposed detection method of abnormities of S1 includes several steps: preprocessing, decomposition of S1 based on EWT, discrimination M1 and T1 from IF and classification of S1 into three classes. After accurately extracting the S1 of each heart sound signal using segmental locating algorithm, this paper mainly focuses on the time-frequency analysis of S1. At first, in order to find M1 and T1, S1 was decomposed into several modes based on EWT. This includes Fourier transform of signal, determining boundaries points, segmentation of Fourier spectrum, decomposition of S1 using Meyer's wavelets filters, Hilbert transform of every mode and drawing a three-dimensional map of IF. Then we designed a clustering algorithm to classify IF in three-dimensional diagram and select the classes which contain the maximal points in Fourier spectrum. In most cases, the modes of the maximum points are in the second and fourth modes. So, the IF of T1 and M1 are corresponding to the second and fourth modes' IF. Finally, some characteristic parameters are calculated based on the three-dimensional figure, such as the time difference of M1and T1 (TD) and the ratio of maximum amplitude of T1 and M1 A peak_ratio . Then S1 was classified into three classes using decision tree classifier: normal S1, S1 with abnormal split, S1 with abnormal amplitude change. Figure 1 illustrates a simplified block diagram of the proposed method for detection of abnormities of S1 based on EWT.
Every heart sound signal will be preprocessed, including resampling, denoising, normalization and segmentation. The purpose of resampling is to reduce the amount of data and improve the speed of data processing. All the data are resampled at a frequency of 1,200 Hz. If heart sound including murmurs or noise, segmentation will become very difficult. So, EWT algorithm is firstly used to separate pure heart sounds (5∼100Hz) from murmurs and noise (>100Hz) [15]. The pure heart sounds are the denoised heart sounds which are prepared for further segmentation to extract S1. Due to the different amplification coefficients adopted by different stethoscope, the amplitude difference of heart sound data collected from a same patient is very large. In order to improve the robustness of the algorithm, pure heart sound data consistency needs to be normalized by using (1). The envelopes of several continuous cycles of pure heart sounds are obtained through Shannon Envelope Extraction. Then the two-threshold segmental localization algorithm is used to localize S1, systole, S2, and diastole accurately in every cycle. Then every S1 of every pure heart sound cycle can be easily extracted for further analyzed. EWT algorithm can be used not only to denoise heart sounds but also to analyze S1 in this paper. It will become very powerful and convenient if you can use it appropriately.
In the next section, how to use EWT to analyze S1 will be emphatically described.

III. PROPOSED METHOD
A. DECOMPOSITION OF S1 BASED ON EWT EWT algorithm divides the Fourier spectrum of a signal to construct some wavelet filter banks and extracts the AM-FM components of the signal with tight support spectrum through orthogonal empirical wavelet transform (EWT). Assuming the signal consists of N AM-FM single components with the spectrum range normalized to [0, ], in order to extract all the individual components, we need to divide the spectrum [0, ] into N continuous segments. In addition to the boundary points 0 and , N-1 boundary points must be identified [11]. In the EWT method, the N-1 boundaries are determined according to the distribution of the maximum points of the spectrum. The continuous interval in Fourier spectrum of a signal can be expressed as [w n−1 , w n ], (n=1, 2, . . . , N), the amplitude is dimensionless and the segmented spectrum is shown in Figure 2. The normalized range of the signal's spectrum is [0, π ], w n is the center of the interval and its width is 2τ n . The shadow part in Figure 2 is the transition region of each segmentation interval, which satisfies (2): Gills only proved that EWT algorithm can decompose the signal into right and meaningful modes when the signal's Fourier spectrum is simple and include several maximum points. However, when the signal's spectrum is complicated and contain too much maximum points, he didn't give out clear methods how to divide the spectrum can ensure getting meaningful decomposition of signal and right instantaneous frequency (IF). V. Nivitha Varghees segmented the spectrum using boundary frequencies (10Hz, 150Hz, 600Hz) [15]. Thus, the spectrum of heart sound will be segmented into 4 IMF components. Then heart sound signal and murmurs (a) Normal S1; (b) Abnormal S1 split; (c) S1 with mitral stenosis; (d) S1 with incomplete tricuspid closure. are divided. However, her method is not suitable for decomposing S1 because S1's spectrum dose not contain these boundary frequencies completely. The frequency of S1 varies from 0-100Hz and its spectrum is not simple. How to divide S1's spectrum using EWT seems to become a major concern in this work.
As is known to all, S1 contains several frequency components, M1 and T1 are the two major frequency components [16]. The time difference between M1 and T1 (TD) is greater than 0.02s directly cause S1 split. The amplitude changes of M1 and T1's IF causes S1's amplitude change [17]. Chinese scholar S. C. Du used Hilbert-Huang transform to analyze the IF of M1 and T1. His statistical results show that the IF of M1 are always greater than T1 over 20Hz [18]. Hence, the abnormities of S1 can be detected and classified through analyzing the changes of M1 and T1.

1) SEGMENTATION OF FOURIER SPECTRUM OF S1
Gills has shown that the segmentations including only one maximum points of the spectrum can ensure getting right IF. Although, S1's frequency varies from 0-100Hz, its Fourier spectrum contains a lot of maximum points. In order to avoid a lot of tedious spectrum segmentations and find the right IF of M1 and T1 accurately and quickly, the peak points in the spectrum were sorted firstly. Then the first two maximum points with a distance greater than 20Hz were selected. Thirdly, the nearest minimum points on both sides of the maximum points were selected as the boundaries points for dividing the spectrum. The process of the segmentation of Fourier spectrum in different S1 were illustrated in Fig.3.

2) DECOMPOSTION OF S1 USING EWT
After the boundaries for the spectrum's segmentation are determined, empirical scale functionφ n (w), and empirical wavelets functionψ n (w) can be obtained according to Meyer's wavelets [11]. The approximate coefficients and detail coefficients obtained by the empirical wavelet transform (EWT) can be expressed as (3) and (4): . Decomposition of the normal S1. (a) Normal S1's decomposition; (b) Abnormal S1 split; (c) S1 with mitral stenosis; (d) S1 with incomplete tricuspid closure.
The original signal can be reconstructed as (5): where, the symbol denotes convolution;Ŵ f 1 (0,w) and W f 2 (n, w) are the Fourier transform of W f 1 (0,t) and W f 2 (n, t) respectively. Finally, the signal x (t) is decomposed into the sum of N single component signals as shown in (6).
Then we can get the different S1's decomposition results using EWT as shown in Fig. 4. They were decomposed into 5 modes because the number of boundaries points for dividing their Fourier spectrum are always 4. In most cases, the spectrum is divided into 5 intervals, T1 and M1 are in the two segments that contain the maximum points. Fig.3(a)-(d) all belong to this situation. So, mode 2 and mode 4 are the waveform of T1 and M1 in time domain because frequency of T1 is usually lower than M1. The decomposition results of different S1 in Fig. 3 were shown in Fig. 4(a)-(d).

3) HILBERT TRANSFORM AND INSTANTANEOUS FREQUENCY
Time-frequency representation is useful to have the information of all components or modes summarized in a single domain. We can apply the same process by evaluating the Hilbert transform of each filter output [11]. Hilbert transform is defined as (7): where, h (t) = 1 πt , * means convolution operation. x (t) is each filter's output andx (t) is its results of Hilbert transform. An analytical signal is constructed withx (t) as the imaginary part, denoted as (8)-(10): This analytical signal's instantaneous frequency can be expressed as (11): After all the modes are evaluated using Hilbert transform, the analytical signal can be expressed as (12): The equation (12) reveals how the signal's IF varies with time and its amplitude. Therefore, we can use (12) to plot a three-dimensional diagram which can clearly shows how the amplitude of IF of each mode change with time. However, when a three-dimensional scatter image is presented, it is impossible to distinguish every dot representing the IF which mode it belongs to. Fig. 5 shows the distribution of a normal S1's IF. Next, a k-mean cluster algorithm was applied to the cluster analysis of IF of different modes.

B. DISCRIMINATION M1 AND T1
Each scatter shown in Fig. 5 include an IF's three parameters: frequency, time and amplitude. We can see the two continuous scatter distributions which look like two sinusoidal waves in Fig. 5. Then a k-mean cluster algorithm was applied in automatically distinguishing these distributions of instantaneous frequency. The number of modes is the number of clusters. The mean value of the discrete instantaneous frequency of every mode will be used as the initial vector of the cluster algorithm. This will greatly reduce the number of iterations of the cluster algorithm. Fig. 6 shows four different S1's IF's distribution after analyzing and processing by the cluster algorithm. According to the segmentation of Fourier spectrum in Fig. 3, the results of clustering are right. The number of clusters is 5. The second cluster and the fourth cluster are consecutive while others are chaotic. Fig.6 further proves that when using EWT algorithm for decomposition of signal, only by ensuring that each segmentation contains one maximum point is IF of each mode correct and meaningful.
According to the clustering results of in Fig. 6(a), the time (td1) at the maximal amplitude of black cluster is earlier than the time (td2) at the maximal amplitude of green cluster (td1). Because M1 is the first audible major component in normal S1 and the frequency of M1 is higher than T1 in normal S1, the black cluster is the distribution of M1's IF and the green cluster is the distribution of T1's IF as shown in Fig. 6(a). Although, Fig.6(b)-(d) are the distribution of abnormal S1's IF and the amplitude and the time of the maximal point of the black clusters and the green clusters will change, the position relationship between M1 and T1 on the frequency domain axis will not change greatly. Therefore, the black cluster is the IF of M1 and T1's is corresponding to green cluster in Fig.6(b)-(d).

C. CLASSIFICATION OF S1
Since the IF of M1 and T1 are discriminated, some characteristic parameters for detection and classification of abnormities of S1 can be calculated. Because the human ear can distinguish splitting sounds with a time difference of greater than 0.02s, it takes a time difference of greater than 0.03s between M1 and T1 as the criterion for S1 split or not. In fact, there are two types of S1 split: physiological split and pathological split. About 70%-85% of healthy people showed interval S1 (0.02-0.03s), and this kind of S1 split is mostly normal physiological split. If the interval between the two components (M1 and T1) is greater than 0.04s, abnormal S1 split occurs.
The criterion for detection and classification of S1 split are mainly according to (13). (13) td1 is the time at the maximal amplitude of black cluster (M1) and td2 is the time at the maximal amplitude of green cluster (T1). S1 split is caused by the delay of a component (M1 or T1), which can be caused by electrical or mechanical delay. Electrical delay is common in bundle branch block. Mechanical delay is mainly mitral valve or tricuspid valve closing delay. In all, the final symptom of S1 split is that the mitral or tricuspid valve closes later than usual. In normal S1, TD is greater than 0 and is smaller than 0.03. If TD<0, it must be mitral valve closure later than tricuspid valve closure, and if the absolute value of TD is greater than 0.04, abnormal S1 split occurs. If TD>0 and the absolute value of TD is greater than 0.04, this kind of abnormal S1 split is caused by the late closure of T1.
Although the loudness of S1 is related to the amplitude of S1 in time domain, the clinical auscultation judgment of the loudness of S1is completely subjective and comparative, and there is no clear objective standard. In addition to shape and thickness of chest wall and lung volume, S1's loudness is affected by the following important factors: the integrity of the atrioventricular valve, the activity of the atrioventricular valve, the position of the atrioventricular valve, and the rate of ventricular contraction. The maximum amplitudes of the FIGURE 6. Distribution of IF of different S1 after using cluster algorithm. (a) Normal S1; (b) Abnormal S1 split; (c) S1 with mitral stenosis; (d) S1 with incomplete tricuspid closure.
M1's and T1's IF, as well as the mitral and tricuspid closing times, are clearly visible in the three-dimensional diagram. Taking the S1 with mitral stenosis as an example, it can be seen that the amplitude of M1 is significantly weakened from Fig. 6(c). Tricuspid valve closed incomplete also appeared the phenomenon of the amplitude of T1's IF increase and M1's weakened. Generally speaking, the factors causing the abnormal change of S1's amplitude are due to the change of M1 or T1's amplitude. The ratio of the maximum value of T1's amplitude to the maximum value of M1's amplitude A peak_ratio can almost determine whether there is abnormal amplitude change in S1. Finally, the abnormal change of S1's amplitude is detected by (14) According to the statistical results of Table 1, it is found that TD of normal S1 is from 0.02s to 0.03s and the A peak_ratio is about from 0.5 to 2, which means the maximal amplitude of M1 and T1 has not too much difference in normal S1. The TD of abnormal S1 split is greater than 0.04s. No matter S1 with mitral stenosis or S1 with incomplete tricuspid closure or S1 with complete atrioventricular block, they all infect the amplitude and the occurrence time of M1 and T1's IF. Finally, the TD and A peak_ratio will change greatly. In Table 1, some cases of S1 with mitral stenosis have the phenomenon of abnormal S1 split since the TD is greater than 0.04s. Some cases' A peak_ratio is higher than 6 while some cases' A peak_ratio is greater than 20. That means in these cases, M1's amplitude decreases or T1's amplitude increases. And the results of M1 or T1's amplitude change are usually an abnormal change in the amplitude of S1.

IV. RESULTS
To evaluate the performance of proposed method for detection and classification of S1, a wide variety of normal and abnormal heart sound data are taken from some databases, including Michgan Heart sounds Database, 2012 Latest Database and Heart Sounds Database built by 3M littmann. The test heart sounds include normal S1, abnormal Split S1, S1 with mitral stenosis, S1 with mitral regurgitation, S1 with tricuspid regurgitation, S1 with incomplete tricuspid closure and S1 with incomplete mitral closure.
Finally, the performance of the proposed method is evaluated using sensitivity (Se), positive predictivity (Pp) and overall accuracy (Oa), which are defined as (15), (16), (17), where TP denote true positive, FP denote false positive and FN denote false negative.
In this study, the performance of the proposed method is compared with the HVD method (the decomposition of signal using HVD method and smoothed pseudo Wigner-Ville distribution for discrimination A2 and P2). The classification results of Se, Pp and Oa using HVD method and proposed method are summarized in Table 2. From the classification results, it is found that the Oa of proposed method is greatly higher than HVD method. At the same time, we demonstrate the accuracy and robustness of the proposed method in detection and classification of abnormities of S1.

V. CONCLUSION
This work applied EWT algorithm in the detection and classification of abnormities of S1. The fundamental concept of this method is to discriminate maximal IF of M1-T1 components by decomposing S1 based on the segmentation of the Fourier spectrum of S1. At first, because M1 and T1 are the two major frequency components of S1, which are generated by the mitral and tricuspid valve closing, first two maximal points with a distance greater than 20Hz in Fourier spectrum of S1 were selected and the nearest minimal points on both sides of the maximal points were found out as the boundaries for segmentation of the spectrum. Then S1 was decomposed into 5 modes using Meyer's wavelets filters based on the 4 boundaries and the IF of every mode was calculated by using Hilbert Transformation. Secondly, the threedimensional scatter diagram of IF of each mode changing with time and amplitude clearly shows two continuous curves which correspond to the IF of M1 and T1. In order to extract the two continuous curves from numerous scatters of IF, a kmean cluster algorithm was applied to cluster the scatters. According to our accumulated knowledge before, in normal S1, M1 appears earlier than T1, and the maximal IF of T1 is usually higher than M1. We can deduce that the green continuous curve in the three-dimensional clustering scatter diagram is the maximal IF of T1. And the black continuous curve is the maximal IF of M1. As can be seen from the three-dimensional diagram, the corresponding frequency of these two continuous curves is constant, and the amplitude of two IF increases at first and then decreases with time.
The peak of the curves corresponds to the moment when the mitral valve and tricuspid valve close. The value of maximal amplitude of two curves also directly reflects the amplitude generated by mitral valve and tricuspid valve closing, which in turn directly affects the amplitude of S1. Therefore, TD and A peak _ratio are used as characteristic parameters to detect and classify S1. Finally, S1s are classified into three classes: normal S1, S1 with abnormal split and S1 with amplitude change. When the proposed method is applied to detect normal S1, Se=94.6%, Pp=98.6% and Oa=93.3%; When the proposed method is applied to detect S1 with abnormal split, Se=92.6%, Pp=96.9% and Oa=90%; When the proposed method is applied to detect S1 with abnormal amplitude change, Se=94.4%, Pp=95.7% and Oa=90.6%; Comparison experiment were carried out between the proposed method and HVD method. The results show Oa of the proposed method is much higher than HVD method when detecting the three different S1.
HAIXIA LI is currently pursuing the Ph.D. degree in instruments and electronic engineering with the North University of China, Taiyuan, China. She is also a Lecturer with the Department of Electronics, Xinzhou Teachers University. Her research interests include design of micro-electro-mechanical systems (MEMS) sensitive units and heart sound signal processing.