Noise Reduction Method of Underwater Acoustic Signals Based on Uniform Phase Empirical Mode Decomposition, Amplitude-Aware Permutation Entropy, and Pearson Correlation Coefficient

Noise reduction of underwater acoustic signals is of great significance in the fields of military and ocean exploration. Based on the adaptive decomposition characteristic of uniform phase empirical mode decomposition (UPEMD), a noise reduction method for underwater acoustic signals is proposed, which combines amplitude-aware permutation entropy (AAPE) and Pearson correlation coefficient (PCC). UPEMD is a recently proposed improved empirical mode decomposition (EMD) algorithm that alleviates the mode splitting and residual noise effects of EMD. AAPE is a tool to quantify the information content of nonlinear time series. Unlike permutation entropy (PE), AAPE can reflect the amplitude information on time series. Firstly, the original signal is decomposed into a series of intrinsic mode functions (IMFs) by UPEMD. The AAPE of each IMF is calculated. The modes are separated into high-frequency IMFs and low-frequency IMFs, and all low-frequency IMFs are determined as useful IMFs (UIMFs). Then, the PCC between the high-frequency IMF with the smallest AAPE and the original signal is calculated. If PCC is greater than the threshold, the IMF is also determined as a UIMF. Finally, all UIMFs are reconstructed and the denoised signal is obtained. Chaotic signals with different signal-to-noise ratios (SNRs) are used for denoising experiments. Compared with EMD and extreme-point symmetric mode decomposition (ESMD), the proposed method has higher SNR and smaller root mean square error (RMSE). The proposed method is applied to noise reduction of real underwater acoustic signals. The results show that the method can further eliminate noise and the chaotic attractors are smoother and clearer.


Introduction
The processing and analysis of underwater acoustic signals play a significant role in the area of military and marine exploration. Due to the influence of the accuracy of measuring instrument and the interference of marine environment, underwater acoustic signals are often contaminated by various noises. The true physical characteristics of underwater acoustic signals are masked, consequently, it is necessary to denoise [1,2]. Underwater acoustic signals are nonlinear, non-Gaussian, and non-stationary chaotic signals [3]. In the time-frequency domain, it has strong broadband properties and pseudo-randomness. This causes the frequency band of the signal to partially or completely overlap with the frequency band of the noise. Therefore, the conventional filtering method is not suitable for noise reduction of underwater acoustic signals [4,5]. In order to effectively suppress additive noise, wavelet analysis, local projection, and singular spectral analysis have been applied to underwater acoustic signals processing, but these methods have certain limitations to different degrees [6,7]. As an example, the wavelet transform has the limits of the selection of the appropriate wavelet basis, decomposition level, and wavelet threshold; the local projection method has the limits of the selection of neighborhood radius; the singular spectral analysis method needs to solve the problem of principal component selection.
Empirical mode decomposition (EMD) [8] provides a new idea of underwater acoustic signals processing. It can decompose the original signal into a series of intrinsic mode functions (IMFs). Each oscillation mode of the original signal represented an IMF. However, EMD has some problems, such as mode mixing, residual noise, and endpoint effects. To overcome these shortcomings, a multitude of methods have been put forward including the masking signal EMD (MS-EMD) [9], ensemble empirical mode decomposition (EEMD) [10], extreme-point symmetric mode decomposition (ESMD) [11], complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) [12], nonlinear mode decomposition (NMD) [13], and improved empirical mode decomposition (IEMD) [14]. Both MS-EMD and EEMD are based on a similar idea: adding certain disturbances to the original signal in order to obtain a perturbed signal with a more uniform distribution of extrema. For example, MS-EMD adds a sinusoidal signal, and EEMD adds Gaussian white noise. Based on EMD, ESMD optimizes the external envelope interpolation method to internal pole-symmetric interpolation, thus, it improves the mode mixing problem. Wang et al. proposed the uniform phase empirical mode decomposition (UPEMD) algorithm [15], which selected sinusoidal functions with uniform phase distribution as perturbed signals, and significantly alleviated the mode splitting and residual noise effects of MS-EMD and EEMD simultaneously.
Shannon entropy is used to quantify the information content of time series. In recent years, researchers have also come up with the approximate entropy (ApEn) [16], fuzzy entropy (FE) [17], sample entropy (SE) [18], and permutation entropy (PE) [19]. PE is based on permutation patterns or order relations among values of a signal, and it has the advantages of simple calculation and relative robustness to observational and dynamical noise. Although the calculation efficiency is high, PE does not consider the amplitude information of the signal. To overcome this problem, many researchers have proposed modified algorithms, such as weighted-permutation entropy (WPE) [20], permutation min-entropy (PME) [21], and amplitude-aware permutation entropy (AAPE) [22]. AAPE is put forward in 2016, which takes into account the amplitude information of the signal and has a wider application prospect.
Recently, an ocean of noise reduction methods based on signal decomposition algorithms and Shannon entropy have been developed and used in different fields, such as acoustic signal [23,24], hydropower unit vibration signal [25], bearing vibration signal [26], medical signal [27], wind speed prediction [28], and so on. Xiao et al. [29] proposed a fault denoising and feature extraction method of rolling bearing based on NMD and continuous wavelet transform (CWT). Zhan et al. [30] used CEEMDAN and FE to denoise the shock signal. The relevant modes (noisy signal modes and useful signal modes) can be identified by the difference between the FE of the new waveform and the next adjacent new waveform. Li et al. [31] proposed a noise reduction method of underwater acoustic signals based on CEEMDAN, mutual information, PE and CWT. Ma et al. [32] proposed EEG signal denoising method based on EEMD and cosine similarity (CS). The similarity between IMF and the original signal is calculated by CS. Then, the IMFs are divided into noise-dominated and signal-dominated part by finding the latter location of the first minimal value of CS curve. Bi et al. [33] used EEMD, robust independent component analysis (RobustICA) and CWT to study the blind source separation and noise source identification of gasoline engines. As far as we know, UPEMD and AAPE have not been used to solve the problem of underwater acoustic signals.
In this paper, a noise reduction method of underwater acoustic signals based on UPEMD, AAPE and Pearson correlation coefficient (PCC) is proposed. We used UPEMD to decompose the original signal into a series of IMFs. According to the threshold of AAPE, high-frequency IMFs and low-frequency IMFs can be distinguished effectively. It is considered that all low-frequency IMFs are useful IMFs (UIMFs), and then PCC is used to determine whether the high-frequency IMF with the smallest AAPE is a UIMF. Finally, the noise reduction can be realized by reconstructing all UIMFs. The structure of this paper is as follows: Section 2 is the basic theory of UPEMD, AAPE, and PCC; in Section 3, the proposed noise reduction method of underwater acoustic signals and evaluation criteria is introduced; in Sections 4 and 5, the proposed method is applied to chaotic signal and underwater acoustic signals respectively. The conclusion is given in Section 6.

Uniform Phase Empirical Mode Decomposition (UPEMD)
The main purpose of the EMD algorithm is to decompose the original signal f (t) into a series of IMFs that characterize the time scale, which can be expressed below: The specific steps of EMD are briefly summarized as follows: Step 1. Connect the local maxima/minima of x(t) to obtain the upper/lower envelope using the cubic spline.
Step 2. Derive the local mean of envelope, m(t), by averaging the upper and lower envelopes.
Step 4. If h(t) satisfies some predefined stoppage criteria [8], h(t) is assigned as an IMF noted as c m (t) where m is the IMF index. Otherwise set x(t) = h(t) and repeat Step 1 to Step 3.
Step 6. Set x(t) = r m (t) and repeat Step 1 to Step 5 to extract the next IMF.
Note that the fixed sifting number n s is one of the stoppage criteria of EMD [8].  [9]. The masking EMD uses a sinusoid signal w(t) as the assisted disturbance with its frequency f w no less than that of the highest frequency component of data. Then EMD is applied to decompose the signal x(t) perturbed by w(t) into two IMFs. The same process is repeated by applying EMD to the signal perturbed by −w(t). Finally, the average of the two resultant sets of IMFs is identified as the final result [9]. Let E m (·) be an operator, which produces the mth IMF decomposed by EMD.
The specific steps of MS-EMD are summarized as follows: Step 1. A masking signal is constructed according to the frequency information of the original signal x(t): Step 2. Compute c + The disturbance amplitude ε and phase θ of the assisted sinusoid signal can be preset according to the relevant rules of the MS-EMD. The specific method is described in [9].

The Two-Level EMD (2L-UPEMD)
MS-EMD uses two sinusoidal signals to eliminate residual noise, but the effect is underdeveloped. Wang et al. [15] pointed out that residual noise can be minimized by searching for all possible phases. Based on this idea, 2L-UPEMD is proposed. Denote the number of phases as n p with n p ∈ N. Let these n p phases be uniformly distributed in the interval [0, 2π]. Then the phase θ k in the k th realization in Equation (2) is calculated as: The specific steps of 2L-UPEMD are summarized as follows: Step 1. Assign T w (= 1/ f w ), ε and n p .
Step 4. Repeat Step 2 to Step 3 for k = 1 to n p .
Step 5. Obtain the resultant IMF1 and IMF2 as c m (t) = 1 n p ∑ k c k,m (t).
The observation formula shows that when n p = 2, 2L-UPEMD is MS-EMD.

The Multi-Level UPEMD
As mentioned in Section 2.1.2, 2L-UPEMD can decompose the signal into two IMFs. The IMF1 is identified as the resultant IMF1, and IMF2 is identified as the new signal. The same method is applied recursively to extract the resultant IMFs at lower frequencies [15].
Assuming a time series of n samples in length. The masking frequency f w is predetermined by taking the dyadic property of EMD that acts as an adaptive dyadic filter bank for the decomposition of the white noise with n s = 10 [10]. The number of resultant IMFs is approximately equal to n im f = log 2 n. The period with index m is determined as T w = 2 m , m = 1, 2, . . . , log 2 n. Let U m (·) be an operator, which produces the mth IMF decomposed by 2L-UPEMD.
The specific steps of multi-level UPEMD are summarized as follows: Step 1. Assign n p , set n im f = log 2 n and initial residue: r 0 (t) = x(t).
Step 3. Perform the 2L-UPEMD to obtain the IMF c m (t), that is c m (t) = U 1 r m−1 ; n p , ε m , (T w ) m , n s .
Step 5. Repeat Step 2 to Step 5 for m = 1 to n im f to extract all IMFs.
Note that the UPEMD used in this paper is the multi-level UPEMD. The disturbance amplitude ε 0 is usually empirically chosen over the range of ε 0 ≈ 0.1∼1.0, n p is the maximum number of phases allowed in each IMF. According to the default values recommended in [15], ε 0 and n p are set to 0.2 and 8, respectively. The Matlab codes of UPEMD are available at http://in.ncu.edu.tw/mzlo/drLo.html.
PE only considers ordinal structure of patterns, the amplitude of each sample is discarded [34]. 26, 8} produce the same ordinal pattern π k = {2, 0, 1}, thus, their PE values are equal. Azami et al. proposed the amplitude-aware permutation entropy (AAPE) for the above problems. Briefly, the repetition probability of each pattern π k is estimated by considering its relative frequency, as well as the average absolute (AA) and relative amplitudes (RA) associated with the vectors Y d . For a specific vector Y d (i), the amplitudes can be expressed as follows: Then, the occurrence probability of π k can be expressed as: where α is an adjusting coefficient, for α ∈ [0, 1], it makes the AAPE algorithm more flexible. In this paper, we set α = 0.5 according to the suggestion in [22]. Finally, the normalized AAPE is computed as follows: The value of d is crucial to the calculation of PE and AAPE. The larger d, the more reliable results are usually obtained. But if d is too large, it takes longer to calculate the entropy. In this paper, based on the length of time series analyzed, we set d = 5.
In order to more intuitively verify and compare the performances of the PE and the AAPE, some quintessential simulation signals are given. The simulation signals f 1 , f 2 , and f 3 can be determined by: where s 1 and s 2 are two components of f 1 , f 2 , and f 3 , respectively, and the sampling frequency is 1000 Hz. The PEs and AAPEs of the three simulated signals are listed in Table 1. As shown in Table 1, for three different simulated signals, the PEs are completely equal and the AAPEs are different. The results of this study indicate that PE cannot distinguish three signals in this case accordingly, and AAPE can better measure time series complexity.

Pearson Correlation Coefficient (PCC)
Pearson correlation coefficient (PCC) is an indicator used to quantify the linear relation between two stochastic variables, X and Y. In order to determine whether the high-frequency IMF with the smallest AAPE is a UIMF, this paper uses PCC to calculate the similarity between the IMF and the original signal. Then the PCC can be expressed as follows: where µ X is the mean of the variable X, µ Y is the mean of the variable Y, σ X is the standard deviation of X and σ Y is the standard deviation of Y. The correlation coefficient can range between −1 and +1 and reveals two aspects about the linear relation between variables: its strength and direction. The direction is determined by the sign, which indicates whether the variables are positive correlated or anti-correlated [35]. The relationship between PCC and correlation is given in Table 2.

The Proposed Noise Reduction Method
The flow chart of the proposed noise reduction method (namely UPEMD-AAPE-PCC) is shown in Figure 1.
The specific steps of the proposed method are as follows: Step 1. Decompose the original signal using UPEMD.
Step 2. Calculate the AAPE of each IMF.
Step 3. Determine the threshold of AAPE. These IMFs, where AAPE is less than a given threshold, are determined as low-frequency IMF (all low-frequency IMFs to be UIMFs). The remaining IMFs are determined as high-frequency IMFs. When the embedded dimension of AAPE is 5, it is appropriate to set the threshold to 0.35, which will be proved in Section 4.1.
Step 4. Calculate PCC between the high-frequency IMF with the smallest AAPE and the original signal. If PCC is greater than 0.4, the IMF is also determined as a UIMF.
Step 5. Reconstruct all UIMFs. After the reconstruction, the process of noise reduction is completed.

The Proposed Noise Reduction Method
The flow chart of the proposed noise reduction method (namely UPEMD-AAPE-PCC) is shown in Figure 1. The specific steps of the proposed method are as follows: Step 1. Decompose the original signal using UPEMD.
Step 2. Calculate the AAPE of each IMF.
Step 3. Determine the threshold of AAPE. These IMFs, where AAPE is less than a given threshold, are determined as low-frequency IMF (all low-frequency IMFs to be UIMFs). The remaining IMFs are determined as high-frequency IMFs. When the embedded dimension of AAPE is 5, it is appropriate to set the threshold to 0.35, which will be proved in Section 4.1.
Step 4. Calculate PCC between the high-frequency IMF with the smallest AAPE and the original signal. If PCC is greater than 0.4, the IMF is also determined as a UIMF.
Step 5. Reconstruct all UIMFs. After the reconstruction, the process of noise reduction is completed.

Evaluation Criteria for Chaotic Signal Noise Reduction
In order to evaluate the performance of the proposed noise reduction method, the proposed method is applied to the denoising experiment of chaotic signal. For the sake of evaluating the performance of these methods, the signal-to-noise ratio (SNR) and root mean square error (RMSE) are adopted and defined as follows [36]:

Evaluation Criteria for Chaotic Signal Noise Reduction
In order to evaluate the performance of the proposed noise reduction method, the proposed method is applied to the denoising experiment of chaotic signal. For the sake of evaluating the performance of these methods, the signal-to-noise ratio (SNR) and root mean square error (RMSE) are adopted and defined as follows [36]: where x andx are the noiseless signal and the denoised signal; N represents the length of the signal.

Evaluation Criteria for Underwater Acoustic Signals Noise Reduction
For the underwater acoustic signals without prior knowledge, the dynamical behavior of denoised signal is analyzed with noise intensity [1], correlation dimension [37] and spatial-dependence recurrence sample entropy (SdrSampEn) [38] to show the validity of the proposed method. Since the original signal contains a multitude of noises, these feature parameters will be larger than the denoised signal, and we use these feature parameters to evaluate the noise reduction effect. Assuming a time series of N samples in length, that is x(n) = {x(1), x(2), . . . , x(N)}. Its noise intensity can be approximated by the standard deviation, i.e.,: x(n) represents the mean of the time series.

Correlation Dimension
The given time series of a dynamic system x k = x(k∆t), (k = 1, 2, . . . , N) is embedded into a m-dimensional space to fully expose the information contained in the time series [37]. For this phase space, it is converted into a set of vectors as follows: where τ = k∆t is the time delay, ∆t is the sampling interval, X n is the N vector in m-dimensional space, n = 1, 2, . . . , M, and M = N − (m − 1)τ is the total number of vectors in reconstructed m-dimensional phase space. Thus, the reconstructed phase space is defined as: First we calculated the distance between all the vectors in the phase space. We define the phase point X i , which is chosen arbitrarily from the M data points, as the reference phase point and calculated the spatial distance d ij between X i and the other M − 1 phase points: where X i and X j are the two points in the phase space reconstruction. After calculating all the points, the correlation integral can be obtained as: where r > 0, and 1 ≤ i ≤ j ≤ M. C(r) represent the percentage of the number of point pairs in the phase space whose radius is smaller than r, and θ(x) is the Heaviside function, which is defined as: Therefore, the correlation dimension D 2 is calculated by the correlation integral C(r) and the radius r with: D 2 is regarded as the slope of the ln C(r)-ln r line segments when r is in a relatively small value and the correlation of ln C(r) with ln r is indicated approximately straight. With a growing m, D 2 increases and becomes saturated gradually. It is worth noting that the calculation of the correlation dimension is more reliable when the length of the data is large enough [39].

Spatial-Dependence Recurrence Sample Entropy (SdrSampEn)
SdrSampEn is based on the mathematical principle of sample entropy and enables the capture of sequential information of a time series in the context of spatial dependence provided by the binarylevel co-occurrence matrix of a recurrence plot. Pham et al. [38] used SdrSampEn to measure the uncertainty of Lorenz signal. At the same time, it is pointed out that SdrSampEn have a better discriminative ability in measuring the irregularity of chaotic time series. Detailed calculation of SdrSampEn can be seen in [38].

The Chaotic Signal Denoising Experiment
Underwater acoustic signals have obvious chaotic characteristics. Three quintessential chaotic signals are selected for simulation experiments, including Lorenz signal, Rossler signal, and Duffing signal. In order to show the superiority of the method, EMD-AAPE-PCC, EMSD-AAPE-PCC and the proposed method are used for comparison.

Choice the Threshold of AAPE
After decomposing the original signal into a series of IMFs, and determining whether they are UIMFs in accordance with AAPE, we need to set the corresponding threshold. Therefore, it is indispensable to study the threshold of AAPE. For the sake of brevity, this section selects the Lorenz signal to analyze.
The Lorenz system can be expressed as: As can be seen in Figure 2, EMD and ESMD decompose the original signal into 9 IMFs, and UPEMD decomposes the original signal into 10 IMFs. The AAPE of each IMF is listed in Table 3. In order to choose the appropriate AAPE threshold, we set the threshold from 0.1 to 1 and the increment is 0.01. When the AAPE is less than this value, the corresponding IMFs are added to obtain a recombined signal. Finally, the SNR of the recombined signal is shown in Figure 3. As indicated in Figure 3, as the AAPE increases, the SNR first increases to a certain extent and then decreases. When the AAPE is around 0.3, the SNR of the three decomposition algorithms reaches the maximum for the first time. As AAPE continues to increase, the SNR of UPEMD begins to decline around 0.4, and the SNR of EMD and ESMD begins to decrease around 0.5. After testing, we found that the threshold selection presented similar results for the Duffing signal and the Rossler signal, which are not listed here. For convenience, this paper sets the AAPE threshold to 0.35. As can be seen in Figure 2, EMD and ESMD decompose the original signal into 9 IMFs, and UPEMD decomposes the original signal into 10 IMFs. The AAPE of each IMF is listed in Table 3. In order to choose the appropriate AAPE threshold, we set the threshold from 0.1 to 1 and the increment is 0.01. When the AAPE is less than this value, the corresponding IMFs are added to obtain a recombined signal. Finally, the SNR of the recombined signal is shown in Figure 3.   As can be seen in Figure 2, EMD and ESMD decompose the original signal into 9 IMFs, and UPEMD decomposes the original signal into 10 IMFs. The AAPE of each IMF is listed in Table 3. In order to choose the appropriate AAPE threshold, we set the threshold from 0.1 to 1 and the increment is 0.01. When the AAPE is less than this value, the corresponding IMFs are added to obtain a recombined signal. Finally, the SNR of the recombined signal is shown in Figure 3.

Denoising for Noisy Chaotic Signal
This section also selects the Lorenz signal to analyze, and the test results of the Rossler signal and the Duffing signal will be given later. The x component signal with a length of 2048 points is selected as a chaotic signal, and the signal is added the Gaussian white noise with different SNRs. The time-domain waveform before and after noise reduction for noisy Lorenz signals when SNR = 10 dB are shown in Figure 4, and their phase diagrams are shown in Figure 5. The noise reduction results for chaotic signal (including Lorenz signal, Rossler signal and Duffing signal with SNR = 5 dB, 10 dB, 15 dB, and 20 dB are listed in Table 4.  As indicated in Figure 3, as the AAPE increases, the SNR first increases to a certain extent and then decreases. When the AAPE is around 0.3, the SNR of the three decomposition algorithms reaches the maximum for the first time. As AAPE continues to increase, the SNR of UPEMD begins to decline around 0.4, and the SNR of EMD and ESMD begins to decrease around 0.5. After testing, we found that the threshold selection presented similar results for the Duffing signal and the Rossler signal, which are not listed here. For convenience, this paper sets the AAPE threshold to 0.35.

Denoising for Noisy Chaotic Signal
This section also selects the Lorenz signal to analyze, and the test results of the Rossler signal and the Duffing signal will be given later. The x component signal with a length of 2048 points is selected as a chaotic signal, and the signal is added the Gaussian white noise with different SNRs. The time-domain waveform before and after noise reduction for noisy Lorenz signals when SNR = 10 dB are shown in Figure 4, and their phase diagrams are shown in Figure 5. The noise reduction results for chaotic signal (including Lorenz signal, Rossler signal and Duffing signal with SNR = 5 dB, 10 dB, 15 dB, and 20 dB are listed in Table 4.   Figure 4 revealed that there are still some noises based on the EMD and ESMD methods, the noise reduction result of UPEMD has the highest similarity with the noiseless signal. In order to better prove the noise reduction effect of the three methods, we compare the phase diagrams before and after noise reduction. As presented in Figure 5a, the phase diagram of the noisy Lorenz signal is annihilated by noise and appears as a disorderly pseudo-random feature. As can be seen from Figure 5d, the phase diagram has a smoother track after noise reduction using UPEMD. Compared with EMD and ESMD, we can see that the denoised phase diagram of UPEMD is closest to the noiseless signal. Table 4 presents the comparison of the SNR and the RMSE of the above three methods. Compared with EMD-AAPE-PCC and ESMD-AAPE-PCC, the proposed method further enhances the SNR and reduces the RMSE. It is obvious that the method in this paper has improving the performance.   Figure 4 revealed that there are still some noises based on the EMD and ESMD methods, the noise reduction result of UPEMD has the highest similarity with the noiseless signal. In order to better prove the noise reduction effect of the three methods, we compare the phase diagrams before and after noise reduction. As presented in Figure 5a, the phase diagram of the noisy Lorenz signal

Data Collection
In this paper, three different types of ship radiated noise are selected as sample data, namely the Ship-I, the Ship-II, and the Ship-III. The ship radiated noise are measured in the South China Sea. Each type of underwater acoustic signals has 100 sample data. Each sample length is 2048 points and sampling interval is 0.05 ms. Before the noise reduction experiment, the sample data have been filtered, normalized, and sampled.

Denoising for Underwater Acoustic Signals
For the sake of brevity, this section selects the Ship-I to analyze. The time-domain waveform of the Ship-I is shown in Figure 6a. The decomposition result of UPEMD for the Ship-I is presented in Figure 6b. The AAPEs and PCCs of these components are indicated in Figure 7.
have been filtered, normalized, and sampled.

Denoising for Underwater Acoustic Signals
For the sake of brevity, this section selects the Ship-I to analyze. The time-domain waveform of the Ship-I is shown in Figure 6a. The decomposition result of UPEMD for the Ship-I is presented in Figure 6b. The AAPEs and PCCs of these components are indicated in Figure 7.  From the AAPE curve in Figure 7, the proposed method determines the high-frequency IMFs and low-frequency IMFs by AAPE, and the AAPE decreases as the IMF order increases. The AAPEs of the last six IMFs (IMF5~IMF10) are all less than 0.35, and they are identified as low-frequency IMFs. From the PCC curve, we find the high-frequency IMF with the smallest AAPE, which is IMF4. It can be found that the PCC between IMF4 and the original signal is greater than 0.4, consequently, Amplitude (normalized) AAPE and PCC

Denoising for Underwater Acoustic Signals
For the sake of brevity, this section selects the Ship-I to analyze. The time-domain waveform of the Ship-I is shown in Figure 6a. The decomposition result of UPEMD for the Ship-I is presented in Figure 6b. The AAPEs and PCCs of these components are indicated in Figure 7.  From the AAPE curve in Figure 7, the proposed method determines the high-frequency IMFs and low-frequency IMFs by AAPE, and the AAPE decreases as the IMF order increases. The AAPEs of the last six IMFs (IMF5~IMF10) are all less than 0.35, and they are identified as low-frequency IMFs. From the PCC curve, we find the high-frequency IMF with the smallest AAPE, which is IMF4. It can be found that the PCC between IMF4 and the original signal is greater than 0.4, consequently, Amplitude (normalized) AAPE and PCC From the AAPE curve in Figure 7, the proposed method determines the high-frequency IMFs and low-frequency IMFs by AAPE, and the AAPE decreases as the IMF order increases. The AAPEs of the last six IMFs (IMF5~IMF10) are all less than 0.35, and they are identified as low-frequency IMFs. From the PCC curve, we find the high-frequency IMF with the smallest AAPE, which is IMF4. It can be found that the PCC between IMF4 and the original signal is greater than 0.4, consequently, IMF4 is also determined as a UIMF. The time-domain waveforms before and after noise reduction of the Ship-I are shown in Figure 8. The phase diagrams before and after noise reduction are shown in Figure 9. IMF4 is also determined as a UIMF. The time-domain waveforms before and after noise reduction of the Ship-I are shown in Figure 8. The phase diagrams before and after noise reduction are shown in Figure 9.   Using the same method, we perform noise reduction on the Ship-II and the Ship-III. The timedomain waveforms before and after noise reduction are presented in Figures 10 and 11, respectively. The phase diagrams before and after noise reduction are presented in Figure 12 and Figure 13, respectively. In order to measure the effect of the noise reduction method for real underwater acoustic signals, we calculate some feature parameters and make a comparison before and after noise reduction, as listed in Table 5. Using the same method, we perform noise reduction on the Ship-II and the Ship-III. The time-domain waveforms before and after noise reduction are presented in Figures 10 and 11, respectively. The phase diagrams before and after noise reduction are presented in Figures 12 and 13, respectively. In order to measure the effect of the noise reduction method for real underwater acoustic signals, we calculate some feature parameters and make a comparison before and after noise reduction, as listed in Table 5. Using the same method, we perform noise reduction on the Ship-II and the Ship-III. The timedomain waveforms before and after noise reduction are presented in Figures 10 and 11, respectively. The phase diagrams before and after noise reduction are presented in Figure 12 and Figure 13, respectively. In order to measure the effect of the noise reduction method for real underwater acoustic signals, we calculate some feature parameters and make a comparison before and after noise reduction, as listed in Table 5.  (a) (b) Figure 11. The time-domain waveform before and after noise reduction for the Ship-III. (a) The time-domain waveform before noise reduction; and (b) the time-domain waveform before noise reduction (select 300 points).
(a) (b)  It can be seen from Figures 8, 10 and 11 that underwater acoustic signals before noise reduction is full of a great deal of noise, and it is completely infeasible to distinguish the trend of the waveform. After noise reduction, the noises are well suppressed. Figures 9, 12 and 13 indicated that the phase diagrams of underwater acoustic signals are obliterated by noise before noise reduction, and fractal characteristics are almost not observed. After noise reduction, the proposed method can roughly restore the phase diagram of the pure attractor. The denoised signal has a smoother orbit, and the attractor has better geometric regularity. As demonstrated in Table 5, the denoised signal has smaller noise intensity, correlation dimension, and SdrSampEn, indicating that the method removes most of the noises in the underwater acoustic signals, further reflecting the effectiveness of the proposed method.

Conclusions
In this paper, a new noise reduction method for underwater acoustic signals based on UPEMD, AAPE, and PCC is proposed. The proposed method uses UPEMD to decompose the original signal, distinguishes the high-frequency and low-frequency IMFs by AAPE, judges whether the high-frequency IMF with the smallest AAPE is a UIMF according to PCC, and finally reconstructs all UIMFs. The main findings in this paper are highlighted as follows: (1) UPEMD, as a new adaptive decomposition algorithm, is first used in the noise reduction of underwater acoustic signals. (2) A simulation experiment shows that AAPE can reflect the amplitude information of the signal compared with PE. Consequently, AAPE is used to measure the complexity of IMF in this paper. (3) Quantitative comparisons based on the noisy chaotic signals demonstrate that the proposed method performs better than the EMD-AAPE-PCC and the ESMD-AAPE-PCC method by providing lower RMSE and higher SNR value. (4) Through the noise reduction experiments of three types of underwater acoustic signals, it is proved that the proposed method can further eliminate the noise and recover the true dynamic characteristics of the chaotic signal more clearly, which lays a foundation for the study of the detection, feature extraction and classification of underwater acoustic signals.