A New Underwater Acoustic Signal Denoising Technique Based on CEEMDAN, Mutual Information, Permutation Entropy, and Wavelet Threshold Denoising

Owing to the complexity of the ocean background noise, underwater acoustic signal denoising is one of the hotspot problems in the field of underwater acoustic signal processing. In this paper, we propose a new technique for underwater acoustic signal denoising based on complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN), mutual information (MI), permutation entropy (PE), and wavelet threshold denoising. CEEMDAN is an improved algorithm of empirical mode decomposition (EMD) and ensemble EMD (EEMD). First, CEEMDAN is employed to decompose noisy signals into many intrinsic mode functions (IMFs). IMFs can be divided into three parts: noise IMFs, noise-dominant IMFs, and real IMFs. Then, the noise IMFs can be identified on the basis of MIs of adjacent IMFs; the other two parts of IMFs can be distinguished based on the values of PE. Finally, noise IMFs were removed, and wavelet threshold denoising is applied to noise-dominant IMFs; we can obtain the final denoised signal by combining real IMFs and denoised noise-dominant IMFs. Simulation experiments were conducted by using simulated data, chaotic signals, and real underwater acoustic signals; the proposed denoising technique performs better than other existing denoising techniques, which is beneficial to the feature extraction of underwater acoustic signal.


Introduction
With the development of ocean scientific technology, the use and protection of the oceans have attracted more extensive attention. Because of the complexity of the marine environment and the time-varying nature of the underwater acoustic channel, it is more difficult to detect and reduce the noise of underwater acoustic signals [1,2]. Therefore, the research on underwater acoustic signal processing method and its application are very essential and important in the field of underwater acoustic. Underwater acoustic signals not only are nonlinear, non-stationary, and non-Gaussian, but also chaos and fractal, traditional signal processing methods based on the classic Fourier analysis are not suitable for underwater acoustic signals such as short-time Fourier transform, Fourier transform, Wigner-Ville, and wavelet transform [3,4]. Therefore, finding a suitable method is the key to analysis underwater acoustic signal. (1) Add white noise n i( t) to the target signal x(t): x i( t) = x(t) + n i( t), i = 1, 2, · · · , N (1) (2) Decompose x i( t) by EMD to obtain the first IMF c i1( t) and residual mode r i (t): · · · · · · c i1( t) r i (t) · · · · · · c N1( t) r N (t) (3) Obtain the first IMF of CEEMDAN by calculating the mean of c i1( t): (4) Obtain the residual mode of c 1( t): (5) Decompose white noise n i( t) by EMD: c n 1 1( t) c n 1 2( t) · · · c n 1 j( t) r n 1( t) c n 2 1( t) c n 2 2( t) · · · c n 2 j( t) r n 2( t) · · · · · · · · · · · · · · · c n i 1( t) c n i 2( t) · · · c n i j( t) r n i( t) · · · · · · · · · · · · · · · c n N 1( t) c n N 2( t) · · · c n N j( t) r n N( t) where c n i j( t) represents the j-th IMF of the i-th white noise, r n i( t) represents the residual mode of the i-th white noise. E j (s i (t)) is defined as a set which includes the j-th IMF of s i (t), E 1 (n i( t)) is expressed as: E 1 (n i( t)) = c n 1 1( t) c n 2 1( t) · · · c n i 1( t) · · · c n N( t) T (6) (6) Construct signal xnew 1 (t) and decompose it by EMD (only decompose the first IMF): xnew 1 (t) = r 1( t) +          c n 1 1( t) c n 2 1( t) · · · c n i 1( t) · · · c n N 1( t) · · · c r 1 n i 1( t) · · · c r 1 n N 1( t) Entropy 2018, 20, 563 4 of 23 (7) Obtain the second IMF c 2( t) and residual mode r 2( t) of CEEMDAN: (8) Obtain xnew j−1 (t) and repeat step (6) and (7), c j( t) and r j( t) are expressed as: (9) x(t) is expressed as: where L represents the number of IMF by CEEMDAN, r(t) represents the residual mode. The flow chart of CEEMDAN is designed in Figure 1.
In this paper, we chose the CEEMDAN algorithm for the following reasons: (1) CEEMDAN has better decomposition effect and lower computational cost than EEMD and CEEMD. (2) CEEMDAN is suitable for analyzing non-linear, non-stationary and non-Gaussian signals, in theory, it can decompose all signals. (3) CEEMDAN is self-adaptive and based on characteristic time scale of the data itself without basis function.

MI
For two discrete random variables X and Y, the MI can be defined as [27]: where p(x, y) is the joint probability distribution function of x and y, p(x) and p(y) are the marginal probability distribution functions of x and y, respectively. The MI of continuous random variables can be expressed as a double integral: )dxdy (16) In probability theory and information theory, the mutual information of two random variables represents a measure of the interdependence of variables. If X and Y are independent, I(X; Y) = 0. In addition, the symmetry of MI can be expressed as: Entropy 2018, 20, 563

of 23
Furthermore, the MI can also be expressed as: where H(X) and H(Y) are informationentropy, H(X|Y) and H(Y|X) are conditional entropy, H(X, Y) is joint entropy of X and Y.
Usually the MIs between noise IMFs are different from ones between non-noise IMFs. We take a ship signal as an example. The normalized ship signal is shown as shown in Figure 2, the sampling frequency and the number of sampling points are 44.1 kHz and 2000, respectively. The decomposition result of the ship signal by CEEMDAN is shown in Figure 3. The center frequency distribution of IMFs and the MIs of two neighboring IMFs are shown in Tables 1 and 2, where MI i represents the MI of IMF i and IMF i+1 . As shown in Tables 1 and 2, the center frequency decreases with the increase of IMF, the first three MI of IMFs are obviously less than the other ones of IMFs. According to the prior information of ship signal, its main frequency range is less than 5000Hz, the first three IMFs are noise IMFs, which is consistent with the judgment of MI. Therefore, we can use MI to identify noise IMFs in this paper, when the MI of IMF i and IMF i+1 increases, obviously more than the former MIs, the former i − 1 IMFs are considered as noise IMFs. Usually the MIs between noise IMFs are different from ones between non-noise IMFs. We take a ship signal as an example. The normalized ship signal is shown as shown in Figure 2, the sampling frequency and the number of sampling points are 44.1 kHz and 2000, respectively. The decomposition result of the ship signal by CEEMDAN is shown in Figure 3. The center frequency distribution of IMFs and the MIs of two neighboring IMFs are shown in Tables 1 and 2, where MIi represents the MI of IMFi and IMFi+1. As shown in Tables 1 and 2, the center frequency decreases with the increase of IMF, the first three MI of IMFs are obviously less than the other ones of IMFs. According to the prior information of ship signal, its main frequency range is less than 5000Hz, the first three IMFs are noise IMFs, which is consistent with the judgment of MI. Therefore, we can use MI to identify noise IMFs in this paper, when the MI of IMFi and IMFi+1 increases, obviously more than the former MIs, the former i − 1 IMFs are considered as noise IMFs.

PE
PE is proposed by Bandt in [24]. The brief process of PE is as follows [28]: (1) Reconstruct time series X = {x 1 , x 2 , · · · , x N }: where τ and m are the time lag and embedding dimension. (2) Rearrange each row vectorin ascending order: (3) Obtain a symbol-sequence for each row vector as: (4) Define PE as: where P j is the probability of one symbol-sequence.
(5) Define normalized PE as: More detail about PE was described previously [29]. In this study, we set m = 3 and τ = 1 according to the suggestiondescribed previously [30]. In a previous paper [21], PE is used to identify noisy IMF. Therefore, in this paper, we choose PE to identifynoise-dominant IMF.

Wavelet Threshold Denoising
Signal denoising is an important research direction of signal processing. The wavelet transform has multi-resolution characteristics. One-dimensional noisy time series can be expressed as follows [31]: where f (k) is original signal, e(k) isnoise signal, s(k) is noisy signal.
Assuming that e(k) is Gaussian white noise, f (k) is usually represented as a low-frequency signal in practical engineering applications. Therefore, we can use the following methods to reduce noise. The specific steps are as follows: (1) A proper wavelet basis function and decomposition level are selected to perform wavelet decomposition on the noisy signal. (2) Threshold is performed by selecting an appropriate threshold method for high frequency coefficients at different decomposition scales. (3) The low frequency coefficient of wavelet decomposition and the thresholdhigh frequency coefficient of different scales are used to reconstruct.
Wavelet thresholding with different thresholds usually has three methods: denoising by default threshold, denoising with specified threshold, and forcing threshold. Among them, denoising with a Entropy 2018, 20, 563 8 of 23 specified threshold is divided into soft threshold and hard threshold. In this paper, a soft threshold method is selected to estimate threshold.

Denoising Algorithm for Underwater Acoustic Signal
The proposed denoising algorithm for underwater acoustic signal is designed in Figure 4. The specific procedures are summarized as follows: (1) The underwater signal is decomposed by CEEMDAN, we can obtain a lot of IMFs, which contain noise IMFs, noise-dominant IMFs, and real IMFs.

Denoising Algorithm for Underwater Acoustic Signal
The proposed denoising algorithm for underwater acoustic signal is designed in Figure 4. The specific procedures are summarized as follows: (1) The underwater signal is decomposed by CEEMDAN, we can obtain a lot of IMFs, which contain noise IMFs, noise-dominant IMFs, and real IMFs.

CEEMDAN for Simulation Signal
Four kinds of simulation signals are selected for denoising, namely, Blocks, Bumps, Doppler, and Heavysine, as shown in Figure 5. The sampling frequency and data lengthare 1 kHz and 1024, respectively.
Taking the Blocks signal as an example, we can obtain the noisy Blocks signal with 0 dB signal-to-noise ratio (SNR) by adding Gaussian white noise. The time-domain waveform of the noisy Blocks signal with 0 dB is shown in Figure 6, and the decomposition result is shown in Figure  7. As shown in Figure 6, the Blocks signal has been completely submerged in noise. The noisy Blocks signal is decomposed using EMD, EEMD, and CEEMDAN.As shown in Figure 7, ten IMFs are obtained by three kinds of methods, however, there are some differences for different

CEEMDAN for Simulation Signal
Four kinds of simulation signals are selected for denoising, namely, Blocks, Bumps, Doppler, and Heavysine, as shown in Figure 5. The sampling frequency and data lengthare 1 kHz and 1024, respectively.
Taking the Blocks signal as an example, we can obtain the noisy Blocks signal with 0 dB signal-to-noise ratio (SNR) by adding Gaussian white noise. The time-domain waveform of the noisy Blocks signal with 0 dB is shown in Figure 6, and the decomposition result is shown in Figure 7. As shown in Figure 6, the Blocks signal has been completely submerged in noise. The noisy Blocks Entropy 2018, 20, 563 9 of 23 signal is decomposed using EMD, EEMD, and CEEMDAN.As shown in Figure 7, ten IMFs are obtained by three kinds of methods, however, there are some differences for different decomposition methods. IMF1 of each decomposition methods represent the shortest oscillation period, typically a noise component or the high frequency components.

Identifying Noise IMFs
In order to observe the effect of noise IMFs on denoising effect, we define the NWn signal as follows: where x(t) and N represent the noisy signal and the number of IMF by CEEMDAN, NWn is the restructured signal by removing the first N IMFs. For the noisy Blocks signal with 0 dB SNR, the six kinds of reconstructed signals are shown in Figure 8 using different decomposition methods.As shown in Figure 8, the noise IMFs is eliminated and the reconstructed signal becomes more smoothwith the increasing of n. When n is larger than a certain value, the non-noise IMF is eliminatedand the reconstructed signal is obviously different from the original signal. Therefore, how to identify noise IMFs is the key problem for denoising. MIs of two neighboring IMFs can expressed as: M n = MI(IMF n , IMF n+1 )(n = 1, 2, · · · , N − 1) (26) where M n represents the MI of IMF n and IMF n+1 . Usually, MI of two noise IMFs is obviously less than the MI of two non-noise IMFs. Therefore, when M n increases obviously, the former n − 1 IMFs can be judged as noise IMFs. For the noisy Blocks signal with 0 dB SNR, MIs of two neighboring IMFs by different decomposition methods are shown in Table 3. As shown in Table 3, M 4 is more than the former ones for EMD and EEMD, we can judge the first three IMFs as noise IMFs. Similarly, the first four IMFs are noise IMFs for CEEMDAN.

Identifying Noise-Dominant IMFs
Noise-dominant IMFs can be identified according to PEs of non-noise IMFs. For the noisy Blocks signal with 0 dB SNR, PEs of non-noise IMFs are shown in Table 4. As shown in Table 4, the PE of IMF5 is more than 0.5, IMF5 is the noise-dominant IMF for CEEMDAN; real IMFs are the last five IMFs.

Denoising for Noise-Dominant IMFs and Reconstruction
The wavelet soft-threshold denoising is applied to IMF5, the wavelet basis function and decomposition level are db4 and 4, respectively. The denoised Blocks signal is obtainedby reconstructing denoised IMF5 and real IMFs. The denoising results are shown in Figure 9. Denoising methods using MI combined with EMD, EEMD, and CEEMDAN are called EMD-MI, EEMD-MI, and CEEMDAN-MI, theproposed denoising method is calledCEEMDAN-MI-PE.  The parameters of different denoising methods are shown in Table 5. As shown in Table 5, theproposed denoising method has lower root mean square error (RMSE) and higher SNR, which outperforms other three denoising methods.

Wavelet Denoising
The wavelet soft-threshold denoising (WSTD) is applied to four kinds of noisy signals with different SNR, wavelet basis function is db4, decomposition level is from 1 to 6. WSTD results are shown in Table 6. As shown in Table 6, SNRs of the four kinds of signals increase with the increasing of decomposition levels. When the decomposition level increases to a certain value, the SNR reaches a maximum. For Doppler and Heavysine signals, when the decomposition level is 5, the denoising results are optimal. For Blocks and Bumps signals with different SNRs, the optimal denoising effects are distributed in different decomposition levels.  Table 7, where WSTD denoising results are optimal values in Table 6. All the results of SNRs and RMSEs are the mean of 500 simulations. As shown in Table 7, the CEEMDAN-MI is better than EMD-MI, EEMD-MI, and WSTD, the CEEMDAN-MI-PE has lower RMSE and higher SNR, which has a better performance than the other four denoising methods.

Denoising for Chaotic Signal
Underwater acoustic signals have the chaotic characteristic, a typical Lorenz chaotic system is used to test the effectiveness of the CEEMDAN-MI-PE denoising algorithm. The Lorenz system can be expressed as: where A is 10, B is 8/3, C is 28. The Runge-Kutta iteration method is used to calculate the x component with a step length of 0.01.The x component signal with a length of 2000 points is selected as Lorenz signal, and the Lorenz noisy signal with different SNR are obtained for CEEMDAN-MI-PE denoising.
Lorenz noisy and denoised signals with different SNRs and their chaotic attractor trajectories are shown in Figure 10. As shown in Figure 10, denoised Lorenz signals and their chaotic attractor trajectories by CEEMDAN-MI-PE are close to Lorenz signal and its attractor trajectory, the denoised chaotic attractor trajectories are more smooth and regular.
Denoising results of different SNR by CEEMDAN-MI-PE are shown in Table 8. As shown in Table 8, the SNR and RMSE are improved evidently, the proposed denoising method enhances the SNR more than 10 dB. Overall, the above results show that the CEEMDAN-MI-PE method is suitable for chaotic signals.

Denoising for Chaotic Signal
Underwater acoustic signals have the chaotic characteristic, a typical Lorenz chaotic system is used to test the effectiveness of the CEEMDAN-MI-PE denoising algorithm.
The Lorenz system can be expressed as: where A is 10, B is 8/3, C is 28.
The Runge-Kutta iteration method is used to calculate the x component with a step length of 0.01.The x component signal with a length of 2000 points is selected as Lorenz signal, and the Lorenz noisy signal with different SNR are obtained for CEEMDAN-MI-PE denoising.
Lorenz noisy and denoised signals with different SNRs and their chaotic attractor trajectories are shown in Figure 10. As shown in Figure 10, denoised Lorenz signals and their chaotic attractor trajectories by CEEMDAN-MI-PE are close to Lorenz signal and its attractor trajectory, the denoised chaotic attractor trajectories are more smooth and regular.
Denoising results of different SNR by CEEMDAN-MI-PE are shown in Table 8. As shown in Table 8, the SNR and RMSE are improved evidently, the proposed denoising method enhances the SNR more than 10 dB. Overall, the above results show that the CEEMDAN-MI-PE method is suitable for chaotic signals.

Denoising for Underwater Acoustic Signal
The CEEMDAN-MI-PE denoising is applied to three kinds of underwater acoustic signals, namely ship-1, ship-2, and ship-3. Three kinds of ship signals were recorded by calibratedomnidirectional hydrophones at a depth of 29 m in the South China Sea. During recording, there wereno observed disturbances from biological or man-made sources. The distance between the ship andhydrophone was about 1 km. The sampling frequency was set as 44.1 kHz. Ship signals and denoised ship signals and their attractor trajectories are shown in Figures 11-13. As shown in Figures 11-13, ocean background noiseis included in original ship signal, high frequency noise is removed effectively by CEEMDAN-MI-PE, denoised attractor trajectories of ship signals are more regular than original ones.    Denoising results of different ships by CEEMDAN-MI-PE are shown in Table 9. Two kinds of PE were used to evaluate the effect of denoising. PE can represent the complexity of time series. A new PE (NPE) was proposed in a previous paper [32], and is interpreted as the distance to noise, which shows a reverse trend to PE. As shown in Table 9, the PE after denoising is less than the one before denoising, which means that the complexity is reduced by denoising; the NPE after denoising is more than the one before denoising, which means that the distance to noise is increased by denoising. In summary, the above results show that the CEEMDAN-MI-PE method is effective and suitable for underwater acoustic signals.

Conclusions
To improve the denoising effect of underwater acoustic signal, a new denoising method is proposed based on CEEMDAN, MI, PE, and WSTD. CEEMDAN is used to decompose noisy signal into IMFs, noise IMFs, and noise-dominant IMFs which can be identified by MI and PE, WSTD is used for denoising noise-dominant IMFs. The innovations and conclusions of the proposed denoising method are as follows: (1) CEEMDAN, as an adaptive decomposition algorithm, is introduced for underwater acoustic signal denoising. (2) Compared with existing denoising methods, IMFs by CEEMDAN are divided into three parts (noise IMFs, noise-dominant IMFs, and real IMFs) for the first time. (3) Four kinds of signals (Blocks, Bumps, Doppler, and Heavysine) with different SNRs are denoised by EMD-MI, EEMD-MI, CEEMDAN-MI, CEEMDAN-MI-PE, and WSTD, the proposed denoising method has lower RMSE and higher SNR, which has a better performance. (4) For chaotic signals with different SNR and underwater acoustic signals, the CEEMDAN-MI-PE is also an effective denoising method, which is beneficial to the subsequent processing of underwater acoustic signals.