Noise Reduction Method of Underwater Acoustic Signals Based on CEEMDAN, Effort-To-Compress Complexity, Refined Composite Multiscale Dispersion Entropy and Wavelet Threshold Denoising

Owing to the problems that imperfect decomposition process of empirical mode decomposition (EMD) denoising algorithm and poor self-adaptability, it will be extremely difficult to reduce the noise of signal. In this paper, a noise reduction method of underwater acoustic signal denoising based on complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN), effort-to-compress complexity (ETC), refined composite multiscale dispersion entropy (RCMDE) and wavelet threshold denoising is proposed. Firstly, the original signal is decomposed into several IMFs by CEEMDAN and noise IMFs can be identified according to the ETC of IMFs. Then, calculating the RCMDE of remaining IMFs, these IMFs are divided into three kinds of IMFs by RCMDE, namely noise-dominant IMFs, real signal-dominant IMFs, real IMFs. Finally, noise IMFs are removed, wavelet soft threshold denoising is applied to noise-dominant IMFs and real signal-dominant IMFs. The denoised signal can be obtained by combining the real IMFs with the denoised IMFs after wavelet soft threshold denoising. Chaotic signals with different signal-to-noise ratio (SNR) are used for denoising experiments by comparing with EMD_MSE_WSTD and EEMD_DE_WSTD, it shows that the proposed algorithm has higher SNR and smaller root mean square error (RMSE). In order to further verify the effectiveness of the proposed method, which is applied to noise reduction of real underwater acoustic signals. The results show that the denoised underwater acoustic signals not only eliminate noise interference also restore the topological structure of the chaotic attractors more clearly, which lays a foundation for the further processing of underwater acoustic signals.


Introduction
The underwater acoustic signals processing is one of the most active subjects in modern information fields [1,2]. Underwater acoustic signal is a non-linear, non-Gaussian, non-stationary chaotic signal, it is easily effected by the other targets, marine environment and various equipment in the process of acquisition and transmission. These problems will inevitably cause some noise for received signal, which is extremely harmful to target detection, location, classification and recognition [3,4]. And self-characteristic of signal is easily ignored in the traditional underwater acoustic signals denoising methods. Therefore, the more effective technology has a wide application prospect in underwater acoustic signal processing [5][6][7].
Because of white noise ε i (t) involved in the decomposition in each experiment is different, so the residual signals are different. The residual signals are defined as:

CEEMDAN
CEEMDAN is an improved algorithm of EEMD. The key technique of this algorithm is that adds adaptive white noise and calculates the specific residuals to obtain corresponding IMFs. E n ( * ) is defined as the n-th modal component by EMD. In this paper, I MF n (t) represents the n-th modal component obtained by CEEMDAN. the specific steps of CEEMDAN are summarized as follows [28]: Step 1: Perform N experiments on the signal x i (t), the first modal I MF 1 (t) is defined as: Step 2: Calculate the first residual, r 1 (t) = x(t) − I MF 1 (t).
Step 3: The i(i = 1, . . . , N) experiments are conducted continuously r 1 (t) + E 1 ε i (t) is decomposed until the first modal component of EMD is obtained. Calculate the second modal component I MF 2 (t): Step 4: Calculate the n-th residual signal and the (n + 1)-th modal component for the remaining phases according to the calculation process of step 3. r n (t) = r n−1 (t) − I MF n (t) (5) Step 5: Continue to execute step 4 until the residual signal is no longer decomposed this criterion is that the number of extreme points of residual signal is less than or equal to 2. the final residual signal is as follows: Therefore, original signal will be eventually decomposed as: Entropy 2019, 21, 11 4 of 21 This decomposition process is complete and adaptive from the algorithm implementation steps of CEEMDAN. This method not only accurately reorganize original signal by adding different white noise also restore some characteristics of EMD while solving the mode mixing problem. The flow chart of CEEMDAN [29] is designed in Figure 1. This decomposition process is complete and adaptive from the algorithm implementation steps of CEEMDAN. This method not only accurately reorganize original signal by adding different white noise also restore some characteristics of EMD while solving the mode mixing problem. The flow chart of CEEMDAN [29] is designed in Figure 1.

Effort-To-Compress Complexity (ETC)
Shannon entropy has been widely used to characterize the complexity of time series caused by stochastic processes such as the chaotic dynamic systems. However, it does not perform well on noise-intensive non-stationary time series. By the time, the measure of compression complexity becomes a great alternative solution [30]. At present, the two methods have been used to measure the complexity of signals, namely Lempel-Ziv test (LZ) and effort-to-compress test (ETC). Many researches have confirmed that LZ and ETC are superior to Shannon entropy in accurately characterizing the dynamic complexity of non-linear dynamic systems. However, the ETC has more different complexity value than LZ and Shannon entropy, it can achieve more sophisticated resolution [31]. Therefore, the application of ETC in non-linear, non-stationary underwater acoustic signals should be promising.
ETC is used to compress a known series by an algorithm, which named non-sequential recursive pair substitution (NSRPS) [32], the specific steps of ETC are summarized as follows: Step 1: Input time series is converted to symbol series so as to compress them by compression calculation, where s 1 , . . . , s n are defined as input series with the number of n.
Step 2: The first iteration is performed when the number of a symbol occurrences is greatest, the series will be replaced by a new symbol. For example, the symbol of "11010010" can be transformed into "12202" caused by "10" has the largest number of occurrences compared with "00", "01" and "11" in the iteration.
Step 3: Complete the second iteration continuously, in which "12202" is converted to "3202". In fact, the frequency of occurrence of all symbols is no difference, we choose to replace "12".
Step 4: The remaining series s 2 , . . . , s n are respectively iterated according to above algorithm until the length of string is reduced to 1 or the series becomes a constant series. Thus, the change of series "11010010" is expressed as: 11010010 → 12202 → 3202 → 402 → 52 → 6 .
Step 5: The value of complexity by ETC is obtained according to execution times of this algorithm. where p is defined as the value of ETC, N represents the number of algorithm that required to convert an input series into a constant series by NSRPS.
where L represents the length of symbol series, N is a non-negative integer from 0 to L − 1, attention should also be paid to:

Dispersion Entropy (DE)
DE is a non-linear dynamic analysis method that characterizes the complexity and irregularity of time series, the algorithm is based on the mapping of normal distribution function. Thus, the expectation and standard deviation of data should be considered, the calculation steps of DE are summarized as follows [33]: Step 1: Define time series is x = x j , j = 1, 2, . . . , N , x is mapped to y = y j , j = 1, 2, . . . , N according to normal distribution function, where y j ∈ (0, 1). and the normal distribution function y j is defined as: where µ and σ respectively represent expectation and standard deviation of time series. Step 2: The y is mapped to the range of [1, 2, . . . , c] by linear transformation.
where R and c respectively represent integer function and the number of categories.
Step 3: Calculate the embedded vector z m,c i : where i = 1, 2, . . . , N − (m − 1)d, m and d respectively represent embedding dimensions and time delays.
Step 4: The dispersion pattern is defined as: Step 5: For each dispersion patterns, relative frequency is defined as follows: where Step 6: According to the definition of Shannon, the DE of original time series is defined as: The parameters selection has been described in Reference [34]. the value of embedding dimensions and categories should be appropriate, m is usually taken as 2 or 3, c is taken as an integer from 3 to 8, d is 1, the length of input time series should be greater than 2000.
From the results of DE, the larger the value of DE, the higher the irregularity of time series, conversely, the irregularity is lower. From the process of algorithm, when all possible dispersion patterns obtain equal probability value, the complexity of signal is the highest, DE obtains the maximum value ln(c m ). If there is a value of p π v 0 ,v 1 ,...,v m−1 and it is not zero, it shows that the lower the complexity of time series, the smaller the value of DE [35].

Multiscale Dispersion Entropy (MDE)
MDE is an interval method based on DE for measuring the complexity and regularity of time series. and the specific steps of MDE are summarized as follows [36,37]: Step 1: For initial time series x, when embedding dimension m and similar tolerance r are respectively determined, the k-th coarse-grained time series with the scale factor τ can be constructed. w τ k is defined as: x a , k = 1, 2, . . . , N/τ where τ is a positive integer. For each scale factor, sample data is divided into several sequences with length N/τ.
Step 2: The MDE with the change of scale factor can be obtained according to the coarse-grained time series, the DE of sample data for each scale factor is defined as: The MDE algorithm fulfills multiscale transformation by equidistant segmentation and averaging of original data. Although its calculation process is simple and fast, there are some relationships between the data after segmentation, which can easily lead to the loss of information.

Refined Composite Multiscale Dispersion Entropy (RCMDE)
In order to solve the above problems, original data is pre-processed on the basis of MDE and then improve the algorithmic process. the specific steps of RCMDE are summarized as follows: Step 1: For original time series x, length( * ) represents the length of signal. the k-th of coarse-grained series is w τ k = w τ k,1 , w τ k,2 , . . . , w τ k,j , . . . , w τ k,N , where k = 1, 2, . . . , τ, w τ k,j is defined as: Step 2: For each scale factor, RCMDE is defined as follows: where p π v 0 ,v 1 ,...,v m−1 represents the average probability of dispersion pattern w τ k , p π v 0 ,v 1 ,...,v m−1 is defined as follows: In this paper, the relevant characteristics of RCMDE is verified by analyzing the simulated signals.  In order to solve the above problems, original data is pre-processed on the basis of MDE and then improve the algorithmic process. the specific steps of RCMDE are summarized as follows: Step 1: For original time series x , ( ) length * represents the length of signal. the k -th of coarse-grained series is Step 2: For each scale factor, RCMDE is defined as follows: where ( ) represents the average probability of dispersion pattern k w τ , ( ) is defined as follows: In this paper, the relevant characteristics of RCMDE is verified by analyzing the simulated signals. And the mean value and standard deviation figures of Gaussian white noise (Noise 1) and   As can be seen from Figure 2, the overall trend of this three methods is basically the same. It is found that the entropy of Noise 1 is larger than that of Noise 2 on the low scale, the entropy of Noise 1 monotonously decreases with the increase of scale. Which indicates that the degree of irregularity of Gaussian noise is higher, the main information is on the low scale. However, the curve change of Noise 2 is not obvious, which indicates that internal structure of Noise 2 is more complex, its main information is not on the low scale. As shown in Figure 2c, the curve of RCMDE is smoother and more stable. Therefore, it has been verified that the RCMDE is more suitable for analyzing these two signals, the stability and accuracy are higher.

Wavelet Threshold Denoising
An effective signal noise reduction method can play a vital role in the field of signal processing. Wavelet analysis developed from Fourier analysis is a new time-frequency analysis tool, which has favorable time-frequency localized and multi-resolution properties. Wavelet analysis has been As can be seen from Figure 2, the overall trend of this three methods is basically the same. It is found that the entropy of Noise 1 is larger than that of Noise 2 on the low scale, the entropy of Noise 1 monotonously decreases with the increase of scale. Which indicates that the degree of irregularity of Gaussian noise is higher, the main information is on the low scale. However, the curve change of Noise 2 is not obvious, which indicates that internal structure of Noise 2 is more complex, its main information is not on the low scale. As shown in Figure 2c, the curve of RCMDE is smoother and more stable. Therefore, it has been verified that the RCMDE is more suitable for analyzing these two signals, the stability and accuracy are higher.

Wavelet Threshold Denoising
An effective signal noise reduction method can play a vital role in the field of signal processing. Wavelet analysis developed from Fourier analysis is a new time-frequency analysis tool, which has favorable time-frequency localized and multi-resolution properties. Wavelet analysis has been widely applied in signal processing field [38][39][40]. The specific steps of wavelet transform are as follows: We suppose that the mathematical expression of one dimensional signal with noise is f (t) = s(t) + δe(t), where t = 0, 1, . . . , (n − 1). s(t), e(t), f (t) and δ are respectively defined as real signal, Gaussian noise, noise signal and the correlation coefficient of noise.
Step 1: A proper wavelet basis function and decomposition level are selected to perform wavelet decomposition on the noisy signal f (t).
Step 2: For the high frequency coefficients obtained by wavelet decomposition, the thresholds are estimated according to the appropriate threshold selection criteria.
Step 3: The high frequency coefficients in different decomposition scales are quantified by thresholds.
Step 4: The low frequency coefficients of wavelet decomposition and high frequency coefficients after processing are reconstructed to obtain denoised signals.
There are many wavelet basis functions and threshold selection criteria in the wavelet analysis, the db4 wavelet basis function is used in this paper. Owing to soft threshold denoising can flexibly overcome the discontinuity of hard threshold among many threshold estimation methods, which has been widely applied to signal processing fields. Therefore, the wavelet soft threshold denoising (WSTD) is applied to this paper.

The Proposed Noise Reduction Algorithm
In this paper, a noise reduction method based on CEEMDAN, effort-to-compress complexity, refined composite multiscale dispersion entropy and wavelet soft threshold denoising is proposed. The flow chart of the proposed algorithm is designed in Figure 3.
The specific steps of the proposed algorithm are as follows: Step 1: An input signal is decomposed into several IMFs by CEEMDAN and arranged from high frequency to low frequency in turn.
Step 2: Calculating the ETC of all IMFs. If the ETC is greater than or equal to threshold p, this IMF will be determined as noise IMF. And set p as 0. 85 after multiple experiments.
Step 3: Calculating the RCMDE of remaining IMFs. If the RCMDE is greater than or equal to threshold q, this IMF is defined as noise-dominant IMFs. If the RCMDE is greater than or equal to threshold r, this IMF is the real signal-dominant IMFs. and the remaining IMFs are judged as the real IMFs, where q is 1.85, r is 1.10.
Step 4: Wavelet soft threshold denoising is applied to noise-dominant IMFs and real signal-dominant IMFs. Owing to the noise signals contain different degrees of noise, the optimal denoising effect is distributed in different decomposition levels. In this paper, wavelet basis function is db4, decomposition level is from one to six, the optimal decomposition level is determined according to signal noise ratio (SNR) and root mean square error(RMSE) of denoised IMFs.
Step 5: In order to ensure more effective denoising effect, the noise IMFs are abandoned. We can obtain the final denoised signal by combining real IMFs and denoised IMFs obtained by wavelet soft threshold. In steps 2 and 3, the algorithm of determining p, q, r is as follows: Step 1: Many researches have confirmed that the noise of signal mainly exists in the high frequency component. In this paper, the IMF1 is a high frequency component by modal decomposition. Thus, the threshold q is set to 0.85 according to the ETC of IMFs.
Step 2: We suppose that the RCMDE of the remaining IMFs has a value range of [ ] Step 3: Firstly, the 1 w is determined to be a fixed value. And then, the value of 2 w is continuously adjusted until the signal-to-noise ratio (SNR) of the reconstructed sequences Z is the largest and the root mean square error (RMSE) of Z is the smallest. Finally, the RCMDE of 2 w is defined as the final value of q .
Step 4: The value of 1 w is continuously adjusted until the SNR of the reconstructed sequences X is the largest and the RMSE of X is the smallest according to the algorithm of Step 3. Therefore, In steps 2 and 3, the algorithm of determining p, q, r is as follows: Step 1: Many researches have confirmed that the noise of signal mainly exists in the high frequency component. In this paper, the IMF1 is a high frequency component by modal decomposition. Thus, the threshold q is set to 0.85 according to the ETC of IMFs.
Step 2: We suppose that the RCMDE of the remaining IMFs has a value range of [m, n], the average value is w. Take two thresholds in [m, n], namely w 1 and w 2 (w 1 < w 2 ), these IMFs are divided into three parts, namely X(m ≤ RCMDE < w 1 ), Y(w 1 ≤ RCMDE ≤ w 2 ), Z(w 2 < RCMDE ≤ n).
Step 3: Firstly, the w 1 is determined to be a fixed value. And then, the value of w 2 is continuously adjusted until the signal-to-noise ratio (SNR) of the reconstructed sequences Z is the largest and the root mean square error (RMSE) of Z is the smallest. Finally, the RCMDE of w 2 is defined as the final value of q.
Step 4: The value of w 1 is continuously adjusted until the SNR of the reconstructed sequences X is the largest and the RMSE of X is the smallest according to the algorithm of Step 3. Therefore, the RCMDE of w 1 is can be defined as the final value of r. After repeated experiments, the values of q and r are roughly determined as 1.85 and 1.10.

Evaluation Method of Chaotic Time Series
In order to evaluate chaotic time series more conveniently, many scholars have proposed some evaluation methods, such as signal-to-noise ratio (SNR), root mean square error (RMSE), correlation dimension, Lyapunov exponent and noise intensity [41,42] and so forth. Thus, in this paper, these evaluation methods are used to evaluate the effect of the proposed noise reduction method.

Signal-To-Noise Ratio (SNR)
Signal-to-noise ratio shows an energy relationship between signal and noise. The higher SNR, the more useful information and the less noise of signal. Therefore, the SNR is a very intuitive method to evaluate the effect of denoised signal by analyzing whether the SNR is improved. The definition of SNR is defined as follows: where x,x and * respectively indicate the noise signal, denoised signal and norm.

Root Mean Square Error (RMSE)
Root mean square error shows the difference between denoised signal and original signal in numerical. and the smaller RMSE, the better noise reduction effect. The RMSE is defined as follows: where length( * ) represents the length of signal.

Correlation Dimension
Fractal dimension is an important parameter to quantitatively analyze the chaotic attractor, which is applied to describe the nonlinear behavior of system. Correlation dimension is a branch of fractal dimension, it has been widely used in signal processing because of simple calculation. In 1983, Grassberger and Procacca proposed the GP algorithm for calculating the correlation dimension of time series [43]. For the time series {x 1 , x 2 , . . . , x n }, let the embedding dimension of reconstructed phase space is m, the delayed sampling is applied to a series with a delay τ. The reconstructed phase space is as follows:  where N = n − (m − 1)τ. For the reconstructed dynamical system, strange attractors are composed of y i = x i , x i+τ , x i+2τ , . . . , x i+(m−1)τ . For any two vectors of y i and y j in phase space, the distance between them are as follow: Suppose there are N vectors in the reconstructed phase space, the correlation integral is defined as: where H(x) is Heaviside unit function.
The correlation integral C(r) has the following relationship with r when r → 0 : where D represents correlation dimension, it can be obtained by calculating D = ln(C(r)) ln(r) .

Noise Intensity
For the time series {x(n)}, whose noise intensity is approximated by standard deviation σ.
where n = 0, 1, . . . , N, x(n) represents the mean of time series, it is found that the smaller noise intensity, the better the noise reduction effect.

Lyapunov Exponent
The Lyapunov exponent judges the chaotic characteristics of the system based on the presence or absence of diffusion motion characteristics of the phase trajectory. Lyapunov exponent is defined as: where n represents the number of iterations, f (x) is the differential equation of the dynamic system and x is the distance of the neighboring points. The positive and negative magnitudes of Lyapunov λ i respectively represent the degree of divergence or convergence of adjacent trajectories in the phase space. Normally, we only need to calculate the maximum Lyapunov exponent. If the maximum Lyapunov exponent is a positive number, we can determine that there is chaotic component in the system. In this paper, the maximum Lyapunov exponent is used to quantitatively analyze the phase space attractors of the signals.

The Chaotic Signal Denoising Experiment
In this section, the Chens model is selected for simulation experiments and added Gaussian white noise with different SNR as input signals. In order to verify the noise reduction effect of the proposed algorithm, two combined noise reduction methods are chosen to compare with CEEMDAN_ETC_RCMDE_WSTD. They are EMD_MSE_WSTD and EEMD_DE_WSTD, the first two methods divide IMFs into two reconstructed series, namely noise-dominant IMFs, real signal-dominant IMFs. and we can obtain the final denoised signal by combining real signal-dominant IMFs. The Chens system is expressed as:       It can be seen from Table 1 that the results of the denoised signals are improved by the three methods. However, the SNR of the proposed algorithm is higher, the RMSE is lower. As shown in Figure 4, the clarity and similarity of time-domain waveform by the proposed algorithm are the highest. Which not only achieves noise reduction but also restores the most of useful information. In Figure 5, although the three methods reduce the noise interference on a certain degree, the geometry of attractor obtained by the proposed algorithm has stronger regularity and higher clarity. In order to quantitatively analyze the phase space attractors of the Chens signal before and after noise reduction, the maximum Lyapunov exponent, correlation dimension and noise intensity before and after noise reduction can be calculated. The Characteristic parameters before and after noise reduction for Chens signal with SNR is 10 dB are shown in Table 2.  It can be seen from Table 1 that the results of the denoised signals are improved by the three methods. However, the SNR of the proposed algorithm is higher, the RMSE is lower. As shown in Figure 4, the clarity and similarity of time-domain waveform by the proposed algorithm are the highest. Which not only achieves noise reduction but also restores the most of useful information. In Figure 5, although the three methods reduce the noise interference on a certain degree, the geometry of attractor obtained by the proposed algorithm has stronger regularity and higher clarity. In order to quantitatively analyze the phase space attractors of the Chens signal before and after noise reduction, the maximum Lyapunov exponent, correlation dimension and noise intensity before and after noise reduction can be calculated. The Characteristic parameters before and after noise reduction for Chens signal with SNR is 10 dB are shown in Table 2. It can be seen from Table 2 that after the above three methods are used to denoise the Chens signal, the above characteristic parameters of the chens signal are significantly improved compared with before the noise reduction, the improvement of the CEEMDAN_ETC_RCMDE_WSTD is most obvious. Thus, it shows that the noise reduction effect of the proposed algorithm is better than the other two methods.

The Underwater Acoustic Signals Denoising Experiment
In order to further verify the effectiveness of this proposed algorithm for chaotic signals. The data used in this paper are three different types of real underwater acoustic signals measured by calibrated omnidirectional hydrophone in the south China sea, namely the Ship-1, Ship-2 and Ship-3. Each type of underwater acoustic signals has 100 sample data. Each sample length is 2048 points and sampling interval is 0.05 ms. The sample data have been filtered, normalized and sampled before the noise reduction experiment, the decomposition results by CEEMDAN are shown in Figure 6. It can be seen from Table 2 that after the above three methods are used to denoise the Chens signal, the above characteristic parameters of the chens signal are significantly improved compared with before the noise reduction, the improvement of the CEEMDAN_ETC_RCMDE_WSTD is most obvious. Thus, it shows that the noise reduction effect of the proposed algorithm is better than the other two methods.

The Underwater Acoustic Signals Denoising Experiment
In order to further verify the effectiveness of this proposed algorithm for chaotic signals. The data used in this paper are three different types of real underwater acoustic signals measured by calibrated omnidirectional hydrophone in the south China sea, namely the Ship-1, Ship-2 and Ship-3. Each type of underwater acoustic signals has 100 sample data. Each sample length is 2048 points and sampling interval is 0.05 ms. The sample data have been filtered, normalized and sampled before the noise reduction experiment, the decomposition results by CEEMDAN are shown in Figure 6. It can be seen from Figure 6 that the three types of underwater acoustic signals are decomposed into several IMFs. The different time scale components are not included in single IMF, the same scale component does not appear in different IMFs. It shows that CEEMDAN does not exhibit mode mixing and boundary effects when applied to underwater acoustic signals, which will make more sense for subsequent processing of underwater acoustic signals. The reconstructed series of underwater acoustic signals are shown in Table 3. Table 3. The reconstruction series of underwater acoustic signals.

Underwater acoustic signals
Noise IMFs

Noise-dominant IMFs
Real signal-dominant  IMFs  Real IMFs   Ship-1  IMF1  IMF2, IMF3  IMF4, IMF5, IMF6, IMF7 IMF8, IMF9, IMF10  Ship-2  IMF1  IMF2, IMF3, IMF4  IMF5, IMF6, IMF7, IMF8  IMF9, IMF10   Ship-3  IMF1  IMF2, IMF3  IMF4, IMF5, IMF6  IMF7, IMF8, IMF9,  IMF10 It can be seen from Table 3 that Ship-1, Ship-2 and Ship-3 are respectively divided into four parts. There is no problem that an IMF is repeatedly defined or undefined. It shows that the proposed algorithm meets the requirements of the IMF. Which will greatly contribute to the noise reduction of the underwater acoustic signals. The time-domain waveform of underwater acoustic signals and phase space attractors after noise reduction are respectively shown in Figures 7-9. It can be seen from Figure 6 that the three types of underwater acoustic signals are decomposed into several IMFs. The different time scale components are not included in single IMF, the same scale component does not appear in different IMFs. It shows that CEEMDAN does not exhibit mode mixing and boundary effects when applied to underwater acoustic signals, which will make more sense for subsequent processing of underwater acoustic signals. The reconstructed series of underwater acoustic signals are shown in Table 3. Table 3. The reconstruction series of underwater acoustic signals.

Signal-Dominant IMFs
Real IMFs   Ship-1  IMF1  IMF2, IMF3  IMF4, IMF5, IMF6, IMF7  IMF8, IMF9, IMF10  Ship-2  IMF1  IMF2, IMF3, IMF4  IMF5, IMF6, IMF7, IMF8  IMF9, IMF10   Ship-3  IMF1  IMF2, IMF3  IMF4, IMF5, IMF6  IMF7, IMF8, IMF9,  IMF10 It can be seen from Table 3 that Ship-1, Ship-2 and Ship-3 are respectively divided into four parts. There is no problem that an IMF is repeatedly defined or undefined. It shows that the proposed algorithm meets the requirements of the IMF. Which will greatly contribute to the noise reduction of the underwater acoustic signals. The time-domain waveform of underwater acoustic signals and phase space attractors after noise reduction are respectively shown in  As shown in Figures 7a, 8a and 9a, the time domain waveform before noise reduction is full of noise, some useful information and the change of time-domain waveform cannot be distinguished. It can be seen from Figures 7b, 8b and 9b that the noise of the underwater acoustic signals are well As shown in Figures 7a, 8a and 9a, the time domain waveform before noise reduction is full of noise, some useful information and the change of time-domain waveform cannot be distinguished. It can be seen from Figures 7b, 8b and 9b that the noise of the underwater acoustic signals are well As shown in Figures 7a, 8a and 9a, the time domain waveform before noise reduction is full of noise, some useful information and the change of time-domain waveform cannot be distinguished.
It can be seen from Figures 7b, 8b and 9b that the noise of the underwater acoustic signals are well suppressed and the waveform change of the denoised signals are clearer by comparing 300 points before and after noise reduction. In addition, we can also determine whether the noise is effectively removed by comparing the chaotic attractors of the underwater acoustic signals before and after noise reduction. Because the degree of damage of the attractor self-similar structure is determined by the noise intensity. The greater the noise of signal, the weaker the regularity of attractor trajectory and the lower the self-similarity. It can be seen from Figure 7c,d, Figure 8c,d and Figure 9c,d that the regularity of denoised signals are stronger, the self-similarity are higher. It shows that the proposed algorithm can reduce the noise interference to a large extent.
In order to further quantitatively describe the effectiveness of the proposed algorithm by calculating the correlation dimension, noise intensity, PE and RCMDE for noise signals and denoised signals. The results before and after noise reduction are shown in Table 4. As shown in Table 4, the change of correlation dimension, noise intensity, PE and RCMDE are smaller than original signals. However, the change of CEEMDAN_ETC_RCMDE_WSTD is the most obvious, which indicates that the noise is greatly suppressed, the complexity is greatly reduced, Therefore, it is shown that the proposed algorithm can not only effectively remove most of the noise also the chaotic characteristics of underwater acoustic signals is greatly improved. Which will have great advantages in processing actual underwater acoustic signals.

Conclusions
In order to solve the problem that inaccurate discrimination of IMFs because of imperfect decomposition process of EMD denoising algorithm and poor self-adaptability, a noise reduction method of underwater acoustic signal denoising based on CEEMDAN, combining ETC, RCMDE and wavelet threshold denoising is proposed. The innovations and conclusions of the proposed denoising method are as follows: (1) CEEMDAN, as an adaptive decomposition algorithm based on EEMD, is introduced for underwater acoustic signal denoising. which has great development potential in the field of non-linear signal processing.
(2) Compared with existing denoising methods, the IMFs by CEEMDAN are divided into four parts (noise IMFs, noise-dominant IMFs, real signal-dominant IMFs and real IMFs) for the first time.
(3) The RCMDE is better than MSE and DE in analyzing the complexity of chaotic signals, is introduced for underwater acoustic signal denoising. Thus, the RCMDE will have greater potential in chaotic signals processing.
(4) The proposed method is applied to Chens model and three different types of real underwater acoustic signals. The proposed method is compared with EMD_MSE_WSTD and EEMD_DE_WSTD, making qualitative and quantitative analysis for denoised signals. The results show that the proposed algorithm can not only reduce the noise interference to a large extent also obtain more regular and clear chaotic attractors. Which will play an important role in researching the physical characteristics of underwater acoustic signals based on chaos theory.