An adaptive thresholding method for the wavelet based denoising of phonocardiogram signal

Highlights

• An adaptive threshold estimation method for wavelet based denoising of phonocardiography signal.

• The proposed method is adaptive to the level of noise present in the signal.

• New statistical parameter is proposed based on domain knowledge about the PCG signal.

• A non-linear mid threshold function is also proposed, and its parameters are optimized.

• The proposed method efficiently suppresses various types of noise and murmur sounds.

Abstract

Segmentation of the phonocardiography (PCG) signal into cardiac cycles is a primary task for the diagnosis of cardiovascular diseases. However, PCG is highly susceptible to noise, and extra sound called murmur may also be present in the PCG signal due to pathology. These components cause difficulties in the segmentation and therefore, segmentation is often preceded by the denoising of the PCG signal to emphasize the fundamental heart sounds S1 and S2, by removing these unwanted components. For the denoising of the PCG signal, discrete wavelet transform (DWT) based algorithms have shown good performance because such algorithms suppress in-band noise besides the out-of-band noise. Selection of threshold value and threshold function significantly affects the performance of these algorithms. In this paper, for threshold value estimation, an adaptive method based on statistical parameters of the given PCG signal is proposed. The statistical parameters are found to be highly effective for this purpose. We also propose a new threshold function, non-linear mid function, to address the issues of SNR and transients in the existing threshold functions, soft and hard. The proposed method is applied on a large number of PCG signals with additive white Gaussian noise, red noise, and pink noise. The Performance of the proposed method is also evaluated on the PCG signals recorded in real-life noisy scenarios and signals with murmur sound. The obtained results show that the proposed method is significantly superior to the competitive algorithms.

Introduction

Heart valvular diseases are manifested into heart sound signal much earlier than they cause serious damage to the heart [1]. It is because the genesis of the heart sound components is related to the movement of the heart valves [2]. Phonocardiography is the most widely used technique for the heart sound signal analysis because of its user-friendly features such as easy to use, timeless set-up process, and cost effectiveness [3], [4], [5]. It uses a sensor called electronic stethoscope, which acquires the heart sound signal and converts it into an electrical form using a microphone. Thus, analysis of the PCG signal is of paramount importance to diagnose heart valvular diseases at an early stage. However, availability of only a small number of experts in this field motivated researchers to develop algorithms for automatic analysis of the PCG signal.

Automatic analysis of the PCG signal is generally performed in two steps; segmentation and classification [6], [7]. In the segmentation step, fundamental heart sounds (FHS), S1 and S2, are identified and then the signal is segmented into systole and diastole periods [7]. Systole is the time duration from S1 to S2 and the diastole is the time duration from S2 to next S1 [6]. In the classification step, the PCG signal is classified as normal or abnormal having a particular disease. Thus, the primary task in the PCG signal analysis is its segmentation [7]. However, PCG signal is highly susceptible to various noises generated due to the motion of subject, subject's own speech and other's speech, movement of the stethoscope used to acquire the PCG signal, lung sounds, and ambient sources [4], [7], [8]. Furthermore, pathologies may cause the presence of extra sound in the PCG signal, called as murmur [6]. The presence of these components makes the segmentation task difficult. Therefore, to emphasize the S1 and S2, segmentation is often preceded by denoising to minimize the contamination level of noise and to remove the murmurs from the signal [9], [10]. However, the presence of murmurs in systole or diastole period and their time-frequency characteristics leads to a diagnosis of heart valvular diseases [11], [12]. For example, murmur due to aortic stenosis, mitral regurgitation and pulmonary stenosis occur in systole period and aortic regurgitation and mitral stenosis occur in diastolic period [6]. The features related to the murmurs can be extracted from the segmented cardiac cycles more efficiently as compared to the signal without segmented cardiac cycles [10]. Therefore, once the PCG signal is segmented into cardiac cycles, pathological features, if present, are extracted from these cycles [10].

Various denoising algorithms for the PCG signal in time domain and frequency domain have been proposed [3]. In the time domain, denoising algorithms have been proposed based on conventional filters such as Chebyshev IIR filter [13], adaptive noise canceller (ANC), and autocorrelation method [14]. The conventional filters are limited to suppress the noise which is out of the frequency band of the signal components. On the other hand, ANC based algorithms such as least mean square (LMS) [15], [16], suppress the noise in an adaptive manner and, hence, suppress in-band noise as well. The major drawback of the ANC algorithms is that they need reference (noise source) signal, which is not available in most cases of the real-life scenarios. Manikandan and Soman [14], proposed a computationally efficient denoising algorithm based on lag-1 autocorrelation method [14]. However, performance of these algorithms significantly degrades as the level of noise increases.

In the frequency domain based denoising algorithms, the time domain signal is first transformed into the frequency domain using a specific transform function such as Fourier transform and wavelet transform (WT), and then the transformed signal is processed. Analysis of the PCG signal in frequency domain provides the information about spectral characteristics of the components presented in the signal and, hence, more efficiency in noise removal can be obtained as compared to the time domain. Sanei et. al. proposed an approach to separate the murmurs from the PCG signal using singular spectrum analysis [11]. In another approach, murmurs were removed from the PCG signal using constrained tunable-Q wavelet transform [17]. However, both the algorithms require high computational time.

The most widely used method for the denoising of the PCG signal is based on the DWT [8], [18], [19] due to the fact that the DWT coefficients of the PCG signal components will be large and they will be confined to specific frequency band, while the coefficients for the noise components will have small amplitude and scattered in different frequency bands [20]. Thus, denoising can be achieved by suppressing the small coefficients. However, the performance of the DWT based denoising algorithm significantly depends on the choice of the parameters: (1) Mother wavelet, (2) Number of decomposition levels and the levels to be processed, (3) Threshold value, and (4) Threshold function [8].

For the denoising, mother wavelet should be orthogonal, which allows perfect reconstruction of the signal [5]. For the PCG signal denoising, various orthogonal wavelets have been suggested such as Coiflet [18], [21], Symlet [12] and Daubechies [2], [22]. In [4], Chourasia et al. developed a new wavelet for PCG signal of a foetus.

The second parameter is the number of decomposition levels. The number of decomposition levels should be selected precisely such that the useful signal components and unwanted components lie in different levels. The frequency range of each level depends on the sampling frequency of the signal. As the sampling frequency increases the frequency range of particular level increases [23]. Therefore, different choices for a number of decomposition levels are reported in the literature [5], [8]. After the decomposition, levels to be processed should be chosen appropriately. In literature, most of the algorithms processed all the decomposed levels [5], [8], [21], [24], which requires unnecessary high computation. In other approaches, the signal is reconstructed using coefficients at a few selected levels, while discarding others [2], [12], [22]. These approaches remove only out-of-band noise. To suppress the in-band noise, the levels associated with the signal components should also be processed [18]. Researchers have also proposed algorithms for the appropriate selection of the levels based on the energy and frequency range of the PCG signal [25], [26].

The third parameter, threshold value, plays a crucial role in DWT based denoising. A large value of threshold affects the useful signal components, while a low threshold value will be ineffective to suppress the unwanted signal components [9]. For the PCG signal, mostly used threshold estimation methods are ‘rigrsure’ [4], [18], [27], ‘heursure’ [5], [8], ‘sqtwolog’ [24], and ‘minimaxi’ [8], [21]. The ‘sqtwolog’ is a fixed form method and does not take into account the content of the signal, but only depends on the length of the signal [4]. It provides a threshold value larger than other methods and hence it may result into over thresholding. ‘minimaxi’ is also a fixed form threshold method, in which the threshold value is estimated such that the maximum risk of estimation error is minimized [5]. ‘rigrsure’ method determines a threshold value to minimize the Steins unbiased risk estimation (SURE). ‘rigrsure’ and ‘minimaxi’ methods estimates threshold value to minimize the risk estimation and results into low threshold value [28]. ‘heursure’ method selects one of the methods from the ‘sqtwolog’ and ‘rigrsure’ methods, based on the comparison between the SURE estimation and SNR [5]. Naseri and Homaeinezhad [22] devised a threshold estimation method based on the weighted variance of the noise while Kumar and Saha [29] calculated the threshold as the 20% of the weighted maximum energy of the coefficient vector. However, values of parameters used in these methods were obtained heuristically.

The fourth parameter is the threshold function, which defines the way to treat the wavelet coefficients using the estimated threshold value. Soft and hard are two existing threshold functions used extensively for the denoising of PCG signals [4], [18], [21]. In soft threshold function, the coefficients lower than the threshold are replaced by zeros while other coefficients get shrank by the threshold value [4]. In hard threshold function, the coefficients lower than the threshold are replaced by zeros, as in soft threshold, while larger coefficients remain unchanged. Hard threshold function may cause discontinuities in the reconstructed signal and make it oscillating [5]. In soft threshold function, shrinkage of the wavelet coefficients by threshold reduces the effect of singularities and transients that cannot be addressed by the hard threshold function [30]. However, hard threshold function produces larger SNR value than the soft threshold function. Zhao [27] proposed a generalized threshold function although it needs a selection of parameter, which controls the performance of the algorithm.

To address the issues discussed above related to threshold estimation and threshold function, we propose a new DWT based denoising algorithm for the PCG signals. We used ‘Coif-5’ wavelet as mother wavelet and performed five levels of decomposition of the PCG signal sampled at 2 kHz. When a signal with 2 kHz sampling frequency is decomposed, the 4th and 5th levels cover the frequency range 31–125 Hz, approximately [31], and hence cover most of the frequency range of the S1 and S2, which is 25–120 Hz [31]. Therefore, in the proposed algorithm, only these two levels are processed. Removal of the coefficients of lower detailed levels also removes the out-of-band noise. To further improve the performance of the denoising, we propose a novel adaptive threshold estimation method using statistical properties of the DWT coefficients. The proposed method uses the domain knowledge that the sum of the length of the S1 and S2 remains less than 25% of the length of a cardiac cycle [32], [33]. Therefore, a new parameter, med₇₅ is calculated instead of traditional median value. The med₇₅ represents the 75th percentile value in the sorted absolute values of a coefficient vector in ascending order. Further, to address the issue of threshold function, we also propose a new method called as ‘non-linear mid’ function for the PCG signal. Furthermore, its parameters are optimized using the genetic algorithm, to improve its performance for the PCG signals. Thus, the proposed algorithm adaptively shrinks the wavelet coefficients of the PCG signal.

Rest of the paper is organised as follow; The proposed denoising algorithm is described in Section 2. The obtained results of the proposed method and its comparison with state-of-the-art methods are discussed in Section 3. Section 4 concludes the result of the paper.