Research on Heart Sound Signal Denoising Algorithm Based on Variational Mode Decomposition and Wavelet Threshold

Qian Liu^1*, Yan Xu², Lianjun Zhang¹, Chunyan Liang¹

¹College of Computer Science and Technology, Shandong University of Technology, Zibo, China.
²Department of Cardiovascular Medicine, Zibo Hospital of Traditional Chinese Medicine, Zibo, China.

DOI: 10.4236/jcc.2021.910007                

Abstract

Heart sound signals are easy to introduce noise during the acquisition process, and traditional denoising algorithms always remove the characteristic information of the heart sound while removing the noise. The denoising effects in turn affect the subsequent diagnosis results. So an improved algorithm based on variational mode decomposition (VMD) and wavelet threshold method is proposed. First, the number of decomposition modes K of the VMD is selected by analyzing the average instantaneous frequency curve of the different decomposition values, and the noisy heart sound is decomposed into K modes by the VMD algorithm. Then, the modes that need to be retained are decided by the energy curve of each mode. Finally, wavelet threshold denoising method is performed on the retained modes. Experiment simulation results show that under different signal-to-noise ratio conditions, the proposed method can improve heart sounds’ ratio of signal to noise and reduce the root mean square error. Compared with traditional algorithms, it has good noise suppression capabilities under different noise levels.

Keywords

VMD, Wavelet Transform, Heart Sound

1. Introduction

The heart sound is a bioelectric signal collected under a strong noise background. The noise will influence the correct diagnosis and recognition of the heart sound. Therefore, the denoising effect will directly affect the subsequent diagnosis and treatment of cardiovascular diseases. So many researchers [1] have done plenty of researches on the heart sound denoising. Wavelet threshold denoising method [2] [3] [4] is often used to eliminate noise components of heart sounds, but the choice of threshold function [5] affects the denoising effect. The empirical mode decomposition (EMD) algorithm [6] can decompose heart sounds into a series of intrinsic mode functions (IMFs), but it is prone to mode mixing. Aiming at the shortcomings of EMD, the ensemble empirical mode decomposition (EEMD) algorithm [7] is proposed and it can solve the problem of mode mixing by adding noise. On the basis of EEMD, the complementary ensemble empirical mode decomposition (CEEMD) further solves the noise problem after adding noise. The article [8] uses the method combining CEEMD and wavelet packet to denoise the heart sounds. However, the above EMD algorithms are essentially a recursive decomposition of the signal in the time domain, and the denoising effect is not ideal. The VMD [9] algorithm was proposed in 2014. This method can realize the adaptive decomposition of each mode, which overcomes the mode mixing problem of EMD. However, this decomposition method needs to select appropriate parameters to achieve a better decomposition effect.

Based on the above problems, an improved algorithm combining VMD and wavelet threshold is proposed. The algorithm firstly chooses the decomposition value K of VMD according to the instantaneous frequency curve of the different decomposition numbers. According to the differences between noise energy and signal energy, the retained IMFs which are decomposed by VMD algorithm can be determined. Finally, the retained heart sound is denoised by the adaptive threshold wavelet transform.

2. Algorithms

2.1. VMD Algorithm

The VMD algorithm [9] is as follows:

1) Initially determine the number K of IMFs. For each mode, $k=1,2,\cdots$, calculate the relevant analytical signal through the Hilbert transform, so as to obtain the unilateral spectrum.

$\left[\delta \left(t\right)+\frac{j}{\pi t}\right]\ast u_k\left(t\right)$. (1)

2) Multiply each eigenmode by, mix it with the corresponding predicted center frequency, and thereby get the “baseband” signal.

$\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]{\text{e}}^{-j{\omega }_{k}t}$. (2)

3) Estimate the bandwidth by the Gaussian smoothness of the demodulated signal. The constrained variational problem is shown in Equation (3).

$\begin{array}{l}\underset{\left\{u_k\right\},\left\{{\omega }_k\right\}}{\min }\left\{{\displaystyle \underset{k=1}{\overset{K}{\sum }}{\Vert {\partial }_t\left[\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]{\text{e}}^{-j{\omega }_{k}t}\right]\Vert }_2^2}\right\}\\ \text{s}\text{.t}\text{.}{\displaystyle \underset{k=1}{\overset{K}{\sum }}{u}_k}=f\end{array}$ (3)

where $\left\{u_k\right\}$ is the set of all eigenmodes, $\left\{{\omega }_k\right\}$ is the set of center frequencies of the eigenmodes.

4) Good convergence is obtained through the quadratic penalty term α, and the constraints are strictly reconstructed through the Lagrangian multiplier $\lambda \left(t\right)$. The introduced augmented Lagrangian is shown in Equation (4).

$\begin{array}{l}L\left(\left\{u_k\right\}\left\{{\omega }_k\right\},\lambda \right)\\ \alpha {\displaystyle \underset{k=1}{\overset{K}{\sum }}{\Vert {\partial }_t\left[\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)\ast u_k\left(t\right)\right]{\text{e}}^{-j{\omega }_{k}t}\right]\Vert }_2^2}\\ +{\Vert f\left(t\right)-{\displaystyle \underset{k=1}{\overset{K}{\sum }}{u}_k\left(t\right)}\Vert }_2^2+\langle \lambda \left(t\right),f\left(t\right)-{\displaystyle \underset{k=1}{\overset{K}{\sum }}{u}_k\left(t\right)}\rangle \end{array}$ (4)

The saddle point of the augmented Lagrangian function is calculated by the alternate direction method of multiplers (ADMM), and the frequency domain update of each mode is obtained. The iterative expression of the mode components and the center frequency is:

${\stackrel{⌢}{u}}_k^{n+1}\left(\omega \right)=\frac{\stackrel{⌢}{f}\left(\omega \right)-{\displaystyle \underset{i\ne k}{\sum }{\stackrel{⌢}{u}}_i\left(\omega \right)+\frac{\stackrel{⌢}{\lambda }\left(\omega \right)}{2}}}{1+2\alpha {\left(\omega -{\omega }_k\right)}^2}$, (5)

${\omega }_k^{n+1}=\frac{{\displaystyle {\int }_0^{\infty }\omega {\left|{\stackrel{⌢}{u}}_k\left(\omega \right)\right|}^2\text{d}\omega }}{{\displaystyle {\int }_0^{\infty }{\left|{\stackrel{⌢}{u}}_k\left(\omega \right)\right|}^2\text{d}\omega }}$, (6)

${\stackrel{⌢}{\lambda }}^{n+1}\left(\omega \right)\leftarrow {\stackrel{⌢}{\lambda }}^n\left(\omega \right)+\tau \left(\stackrel{⌢}{f}\left(\omega \right)-{\displaystyle \underset{k}{\sum }{\stackrel{⌢}{u}}_k^{n+1}\left(\omega \right)}\right)$. (7)

where n is the number of iterations and $\tau$ represents the fidelity coefficient.

2.2. Wavelet Threshold

How to select the threshold and threshold function [10] is a difficult problem in wavelet transform. Threshold functions commonly used are hard threshold function and soft threshold function. The expression of the hard threshold function is:

${\stackrel{⌢}{w}}_{j,k}=\{\begin{array}{l}{w}_{j,k},\left|w_{j,k}\right|>\mu \\ 0,\left|w_{j,k}\right|\le \mu \end{array}$, (8)

In Equation (8), ${\stackrel{⌢}{w}}_{j,k}$ is the denoising form of $w_{j,k}$, $\mu$ is the corresponding threshold.

The expression of $\mu$ is:

$\mu ={\sigma }_j\sqrt{2\ln N}$, (9)

where N is the length of the signal, and ${\sigma }_j$ is the noise estimation variance of thejth layer. The variance of the noise estimate is as follows:

${\sigma }_j=\frac{median\left|w_{j,k}\right|}{0.6745}$. (10)

The expression of the soft threshold function is:

${\stackrel{⌢}{w}}_{j,k}=\{\begin{array}{l}sign\left(w_{j,k}\right)\left(\left|w_{j,k}\right|-\mu \right),\left|w_{j,k}\right|>\mu \\ 0,\left|w_{j,k}\right|\le \mu \end{array}$ (11)

The hard threshold causes discontinuity in denoising signal and although the soft threshold function will be smoother, it will cause the reconstructed signal to be close to the original signal not well. Based on this, a threshold function [11] is proposed and the expression is as follows:

${\stackrel{⌢}{w}}_{j,k}=\{\begin{array}{l}{w}_{j,k}-0.5sign\left(w_{j,k}\right)\frac{{\mu }^m}{{\left|w_{j,k}\right|}^{m-1}},\left|w_{j,k}\right|>\mu \\ 0.5sign\left(w_{j,k}\right)\frac{{\left|w_{j,k}\right|}^{m+1}}{{\mu }^m},\left|w_{j,k}\right|\le \mu \end{array}$ (12)

In Equation (12), m is a continuous real number which is greater than 1. While $m=1$, the function almost coincides with the soft threshold function. While $m>10$, the function is very close to the hard threshold function.

3. Improved Algorithm

Select the appropriate decomposition value of VMD can effectively solve the mode mixing problem. Make full use of the energy differences between noise and signal to retain the signal IMFs. Finally the adaptive wavelet threshold method is used to denoise the heart sounds.

3.1. Selection of Decomposition Value

The VMD algorithm is equivalent to an adaptive Wiener filter bank. Therefore, the decomposition effect of VMD is mainly affected by the value K. When K is too small, some important signal components in the original heart sound will be filtered, and it influences the effect of subsequent denoising. While K is too large, the center frequencies of adjacent mode components will be closer together, resulting in repeated mode or additional noise.

Based on the problems, the algorithm first analyzes the average instantaneous frequency curve [12] of each mode under the condition of different decomposition numbers. Figure 1 is the average instantaneous frequency curves of the AS (Aortic Stenosis) heart sound with different K, and its signal-to-noise ratio is 0 dB. It can be seen from the figure that whenK is 2 - 5, the curve is approximately a straight line, indicating that the mode component is approximately proportional to the average instantaneous frequency. However, the frequency of the heart sound is low, and if the average instantaneous frequency is far apart, it will cause too much noise components in the low-frequency mode. As the value of K increases, the mode decomposition is more adequate, and there might be cause the problem of mode mixing problem. Based on the above analysis, combined with Figure 1, the article selects K = 6 as the number of decompositions.

Figure 1. Average instantaneous frequency curve (AS).

After a large number of simulations, it is found that if the heart sound is a non-pathological heart sound, the change of instantaneous curve is not obvious with the increase of decomposition numbers, which is approximately a straight line. In this case, whenK is 8, compared with other decomposition numbers, the denoising effect is better and the signal-to-noise ratio is higher. Therefore, if the instantaneous frequency average curve is approximately a straight line at different decomposition numbers, the number of decompositions of VMD selects 8.

3.2. Energy Analysis

The frequency of normal heart sound is 20 Hz - 200 Hz, and the frequency of heart murmurs generally does not exceed 800 Hz. The VMD algorithm decomposes the heart sound into different modes which frequencies are from low to high, and the frequency spectrum is concentrated on the lower-order modes. So the signal energy is very large at the low-order components. That the energy based on signal and noise behaves differently in the IMFs indicates that the signal energy starts from strong to weak. Therefore, the point where the low-order energy drops rapidly is marked as the signal and noise boundary point. As the Figure 2, the AS is decomposed into 6 layers, and the energy of the first IMF is the largest, the energies of the second to sixth IMFs are almost indistinguishable. So the signal of the first mode is used as the signal for subsequent wavelet denoising.

3.3. The Algorithm Flow

In summary, the algorithm steps are as follows:

Step1: Decompose the noisy heart sound signal into different K, where K = 2 - 13. Analyze the average instantaneous frequency value curve, find the K value and record it as the decomposition value.

Figure 2. The energy curve.

Step2: Decompose the heart sound signal using VMD algorithm, where the K value is selected by step 1.

Step3: Retain the IMFs to be denoised based on the energy analysis method.

Step4: Select the appropriate wavelet to decompose the IMFs, and perform wavelet threshold according to the threshold function of Equation (12).

Step5: The IMFs component denoised by wavelet threshold is the heart sound signal after denoising.

The implement process is shown in Figure 3.

4. Simulation Results

The heart sound samples used in this article are from International Heart Sound Competition challenge2016. In the MATLAB 2016a simulation environment, this article uses the proposed algorithm and five traditional algorithms of CEEMDAN, EEMD, VMD and wavelet threshold to analyze and compare the denoising effects of heart sound under different signal-to-noise ratio conditions.

4.1. The Denoising Result under Different Noise Conditions

In order to intuitively reflect the denoising effect of the algorithms under different noise conditions, the heart sound (S4) with SNR from −5 dB to 15 dB is selected for denoising experiment. The denoising effect is shown in Table 1.

The mathematical expression of SNR is as follows:

$SNR=10\lg \frac{{\displaystyle \underset{n=1}{\overset{N}{\sum }}{x}^2\left(n\right)}}{{\displaystyle \underset{n=1}{\overset{N}{\sum }}{\left(\hat{x}\left(n\right)-x\left(n\right)\right)}^2}}$. (13)

Figure 3. The flow chart of the algorithm.

Table 1. The denoising result under different noise conditions.

The mathematical expression of RMSE is as follows:

$RMSE=\sqrt{\frac{1}{N}{\displaystyle \underset{n=1}{\overset{N}{\sum }}{\left(\hat{x}\left(n\right)-x\left(n\right)\right)}^2}}$. (14)

From Table 1, the denoising effect of the algorithm in this artical is still significantly proposed under the low signal-to-noise ratio. RMSE decreases with the increase of signal-to-noise ratio, indicating that the proposed algorithm retains the characteristics of the original heart sound signal while removing noise.

4.2. Analysis of Different Algorithms

The above 5 algorithms are used to denoise the different heart sound signals (NHS, ASD, AI, AS, S4) under the condition of 0dB. Algorithm 1 uses the denoising method based on this article, Algorithm 2 uses the combination of CEEMDAN and wavelet threshold method, Algorithm 3 uses the combination of EEMD and wavelet threshold method, Algorithm 4 uses the VMD method, Algorithm 5 uses wavelet threshold method. Based on the simulations of many experiments, the wavelet base in this article adopts “db10”, the wavelet decomposition level is 5, and the VMD penalty factor is 2000. The denoising results of the 5 methods are shown in Table 2.

From Table 2, whether NHS or pathological heart sounds ,the denoising effect of the proposed algorithm is better, and the signal-to-noise ratio is higher among the above 5 algorithms. For pure VMD, the signal-to-noise ratio has also been significantly improved, but the high-frequency noise is not filtered out. However, purely wavelet threshold denoising will filter out the heart sound information while filtering out the noise, resulting in a lower signal-to-noise ratio.

In order to further verify the effects of different denoising algorithms, this article uses the heart sound signals of NHS and AS as the experimental objects, and uses SNR and RMSE curves to analyze the denoising effects. The output results are shown in Figure 4. Under the same SNR environment, the wavelet

Table 2. The SNR value of different types of denoised heart sound.

(a)(b)(c)(d)

Figure 4. SNR and RMSE curve. (a) SNR curve (NHS); (b) RMSE curve (NHS); (c) SNR curve (AS); (d) RMSE curve (AS).

threshold has the worse denoising effect. The denoising effects of CEEMDAN and EEMD algorithms are relatively close. As the signal-to-noise ratio increases, the VMD becomes more effective. Obviously, compared to the other four algorithms, the proposed algorithm in this paper has a more significant increase in signal-to-noise ratio, a smaller root mean square error, and a better retention of the characteristics of heart sound signals.

Figure 5 is a graph showing the denoising effect of the proposed algorithm on different abnormal heart sounds. From the figure, it is known that the denoising effect is more obvious. All in all, the proposed denoising algorithms can filter out most of the noise and retain the basic characteristics of the original pathological heart sound, so the denoising effect is obvious.

(a)(b)(c)

Figure 5. Abnormal heart sound and the denoised singals. (a) AS and denoised signal; (b) ASD and denoised signal; (c) AI and denoised signal.

5. Conclusions

1) In the theoretical part, VMD algorithm is used to decompose the noisy heart sound signal. The article proposes to use the instantaneous frequency curve to select the decomposition number K, and use the differences between the signal energy and the noise energy to determine the final retained mode components. Finally, the retained signal is processed by wavelet adaptive threshold algorithm for denoising.

2) In the experimental part, compare the denoising effects of the algorithm in this article with the other 4 algorithms on normal heart sounds and abnormal heart sounds under different signal-to-noise ratio conditions. By analyzing, whether for normal or abnormal heart sounds, the denoising effect of the article is better than the other four traditional denoising algorithms, and it has obtained a higher signal-to-noise ratio, a lower root mean square error, and better to retention the characteristics of the heart sound signal.

The algorithm in this article achieves effective noise filtering, provides a good data basis for subsequent heart sound recognition and cardiovascular treatment, and also provides a way of thinking for future signal denoising methods.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Liu, F., Wang, Y.T. and Wang, Y.X. (2012) Research and Implementation of Heart Sound Denoising. Physics Procedia, 25, 777-785. https://doi.org/10.1016/j.phpro.2012.03.157
[2]	Jain, P.K. and Tiwari, A.K. (2017) An Adaptive Thresholding Method for the Wavelet Based Denoising of Phonocardiogram Signal. Biomedical Signal Processing and Control, 38, 388-399. https://doi.org/10.1016/j.bspc.2017.07.002
[3]	Misal, A. and Sinha, G.R. (2012) Denoising of PCG Signal by Using Wavelet Transforms. Advances in Computational Research, 4, 46-49. https://doi.org/10.26634/jcs.1.1.1740
[4]	Ali, M.N., El-Dahshan, E.S.A. and Yahia, A.H. (2017) Denoising of Heart Sound Signal Using Discrete Wavelet Transform. Circuits Systems & Signal Processing, 36, 4482-4497. https://doi.org/10.1007/s00034-017-0524-7
[5]	Zhou, K.L. and Liu, Y.Y. (2020) Heart Sound Denoising of New Threshold Wavelet Transform. Computer Engineering and Design, 41, 2476-2481. https://doi.org/10.16208/j.issn1000-7024.2020.09.012
[6]	Salman A.H., Ahmadi, N., Mengko, R., Langi, A.Z.R. and Mengko, T.L.R. (2016) Empirical Mode Decomposition (EMD) Based Denoising Method for Heart Sound Signal and Its Performance Analysis. International Journal of Electrical and Computer Engineering, 6, 2197-2204. https://doi.org/10.11591/ijece.v6i5.pp2197-2204
[7]	Gaci, S. (2016) A New Ensemble Empirical Mode Decomposition (EEMD) Denoising Method for Seismic Signals. Energy Procedia, 97, 84-91. https://doi.org/10.1016/j.egypro.2016.10.026
[8]	Dong, L.C., Guo, X.M. and Zheng, Y.N. (2019) Wavelet Packet De-Noising Algorithm for Heart Sound Signals Based on CEEMD. Journal of Vibration and Shock, 38, 192-198+222. https://doi.org/10.13465/j.cnki.jvs.2019.09.025
[9]	Dragomiretskiy, K. and Zosso, D. (2014) Variational Mode Decomposition. IEEE Transactions on Signal Processing, 62, 531-544. https://doi.org/10.1109/TSP.2013.2288675
[10]	Donoho, D.L. (1995) De-Noising by Soft-Thresholding. IEEE Transaction on Information Theory, 41, 613-627. https://doi.org/10.1109/18.382009
[11]	Wu, G.W., Wang, C.M., Bao, J.D., Chen, Y. and Hu, Y.P. (2014) A Wavelet Threshold De-Noising Algorithm Based on Adaptive Threshold Function. Journal of Electronics & Information Technology, 36, 1340-1347.
[12]	http://www.cnblogs.com/xingshansi/p/6511916.html