# Application of SN-EMD in Mode Feature Extraction of Ship Radiated Noise

Authors: Fang Niu; Juan Hui; Anbang Zhao; Yue Cheng; Yang Chen
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Due to the randomness of added noise, noise-assisted versions based on EMD (empirical mode decomposition) usually cause new “mode mixing” problem. In addition, these algorithms also have problems such as high time-consuming and large recovering error. For the reasons, a new method SN-EMD (Selective Noise-assisted EMD) is put forward in this paper. It determines whether to add noise as assistance by judging whether there is high frequency intermittent component contained in the signal or not. The new method was proved to have the optimal performance by comparing the performance parameters for evaluating the decomposition. In this paper, SN-EMD was used to decompose ship radiated noise. On account of the differences in the original information contained in each mode of radiated noise signals from different ship, we selected the first three modes for processing. Average instantaneous frequency, center frequency, energy density, and energy distribution ratio were extracted as mode feature of ship targets for classification and recognition. Spatial distribution of the feature quantities in three-dimensional space verified similarity of the same target and separability of different targets.


---

## Body

## 1. Introduction

EMD [1] is an adaptive signal processing method proposed by N. E. Huang, which is especially suitable for nonlinear and nonstationary signal. It can decompose the original signal in an adaptive manner in fast and slow oscillations based on the direct extraction of signals with different time scale. The original signal can be decomposed into a series of monocomponent signals and a component in a monotonic trend. These monocomponent signals are usually referred to as “intrinsic modal functions (IMFs)” (or modes). Although the modes obtained by EMD are local narrowband and can represent the actual physical process, therefore, analysis of instantaneous frequency by Hilbert transform can be performed for each mode.The real signal is complex, and intermittency of frequency algorithm causes mode mixing in EMD. Namely, signal components on different time scales are decomposed into a mode, or signal components at the same time scale are decomposed into different modes. In order to solve the problem of mode mixing, based on the EMD algorithm, a series of noise-assisted algorithms named Ensemble Empirical Mode Decomposition (EEMD) [2], Complementary Ensemble Empirical Mode Decomposition (CEEMD) [3], Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN) [4], and improved CEEMDAN [5] were proposed consecutively. Among them, the improved CEEMDAN has the best results. There is no problem that the number of modes is different due to the additive noise in EEMD and CEEMD, and the problem of false modes in CEEMDAN. However, the improved CEEMDAN still has aspects for improvement: (I) time-consuming issues: in order to reduce the residual error, it needs hundreds of decomposition over an ensemble of noisy copies of the original signal, so the algorithm is time-consuming; (II) the problem of mode mixing: for the reason of added noise, this algorithm can avoid decomposing the components on different time scales into the same mode, but there is still a second type of mode mixing phenomenon that the components at the same time scale are decomposed into different modes. To solve the two problems above, we introduce an improved algorithm SN-EMD.The decomposition algorithms have applications in many fields, such as electrocardiograph (ECG) signal denoising and remote monitoring of the respiration and heartbeat of a patient [5, 6], time-frequency analysis of seismic signals [7, 8], analysis of the regional features of temperature variation trends [9, 10], electrical filed such as power quality estimation [11, 12], mechanical fault diagnosis of roller bearings [13, 14], polymer electrolyte membrane fuel cell fault diagnosis in chemistry [15], visual enhancement of underwater images, and nonlinear feature extraction in underwater acoustics [16, 17]. Though the actually ship radiated noise is Gaussian-like in terms of probability distribution [18], it contains many definitive ingredients. The center frequency distribution of the modes of Gaussian noise usually obeys the simple bisection rule, while ship radiated noise does not [19, 20]. The difference of center frequency distribution of modes reflects the characteristics of different types of ship. In addition, the energy distribution of the mode of each order also has its own rules. As a result, the features such as instantaneous frequency and instantaneous amplitude extracted from each mode can be used as the basis for classification and identification of ship targets.This paper is organized as follows: Section2 introduces the existing empirical mode decomposition algorithms and the new algorithm SN-EMD. Section 3 compares the performance of these decomposition algorithms. In Section 4, instantaneous frequency and other factures are extracted from modes of real ship radiation noise data, and their separability is verified in a three-dimensional space. The conclusions are drawn in Section 5.


## 2. Methods

### 2.1. Existing Methods

#### 2.1.1. EMD

The EMD algorithm can be described as follows.Step 1.
Setk=0 and rk(t)=x(t).Step 2.
Find all the maxima (minima) ofrk(t). Interpolate between maxima (minima) of rk(t) to obtain the upper (lower) envelope emax(t)(emin(t)).Step 3.
Compute the mean envelopem(t)=(emax(t)+emin(t))/2.Step 4.
Compute the IMF candidateck+1(t)=rk(t)-m(t).Step 5.
Judge whetherck+1(t) is an IMF.(i)
Yes. Setimfk+1⁡(t)=ck+1(t) and save imfk+1⁡(t) as IMF of order k+1. Compute the residue rk+1(t)=x(t)-∑n=1k+1imfk+1⁡(t), do k=k+1, and treat rk(t) as input data in Step 2.(ii)
No. Treatck+1(t) as input data in Step 2.Step 6.
Continue until the final residuerk(t) satisfies some predefined stopping criterion.The original signalx(t) can be expressed as a sum of imfn⁡(t),(n=1,…,N).There is mode mixing problem in EMD decomposition, which makes it impossible for EMD to achieve separation of components of the signal at different frequency. The mode obtained by decomposition cannot reflect the instantaneous frequency characteristics of the signal and fundamentally loses the original intention of EMD algorithm. Meanwhile, the deficiency limits the engineering application of the algorithm.


#### 2.1.2. EEMD

In order to solve the problem of mode mixing, Wu, Z. [2] proposed the EEMD algorithm. EEMD defines the “true” modes as the average of the corresponding IMFs obtained from an ensemble of the original signal plus different realizations of finite variance white noise. On the basis of the characteristics of the EMD’s binary filter bank, it reduces mode mixing by filling the entire time-frequency space, with adding white noise to the signal.EEMD algorithm also has some drawbacks: (I) Different realizations of finite variance white noise may produce different numbers of modes; (II) the average mode may not meet the requirements of effective IMF; (III) when IMFs are used to reconstruct the signal, the residual noise cannot be ignored; (IV) ensemble average in EEMD usually needs operation for several hundred times, which is very time-consuming.


#### 2.1.3. CEEMD

Yeh, J.R. [3] proposed the CEEMD algorithm, which is different with EEMD in that the noise added to the original data is in pairs (a positive number and a negative number). Adding noise with a pair of opposite signs makes the residual noise basically cancel out each other and achieve the effect of “complementary.”


#### 2.1.4. CEEMDAN

CEEMDAN proposed in [4] is an important improvement to EEMD. The first-order mode using CEEMDAN algorithm is the same as that using EEMD, and other signal modes are the first-order mode decomposed from the residual signal plus white noise mode in the same order. Also, the coefficients of each noise mode are all adaptive, which are related to the residual signal and noise mode. The improvement of CEEMDAN reduces the reconstruction error and solves the problem that different realizations of signal plus noise might produce different number of modes.


#### 2.1.5. Improved CEEMDAN

In order to reduce the residual noise, Colominas, M.A. [5] made an improvement on CEEMDAN: estimating local mean instead mode. From the first-order mode, the noise added in the signal is mode component decomposed from white noise rather than the original white noise.


### 2.2. SN-EMD

When there are high frequency intermittent components contained in the signal, the EMD will lead to the first type of mode mixing problem that signal components on different time scales in one mode. In order to solve the problem, a series of noise-assisted versions based on EMD have been proposed. However, due to the randomness of the added noise, the noise-assisted algorithms may cause the second type of mode mixing phenomenon that the signal components on the same time scale are decomposed into different modes. In addition, these algorithms also have problems such as high time-consuming and large residual error. For the reasons, a new method SN-EMD (Selective Noise-assisted EMD) is proposed in this paper. Therefore, in order to avoid the phenomenon of mode mixing, it is necessary to judge whether high frequency intermittent component is contained before decomposition and then decide whether to add noise for decomposition or not.At the high frequency intermittent component of the signal, the interval time between the extrema of the signal envelope will decrease obviously. That is, at the beginning and the end of the intermittent component, the interval time between the extrema has a jump phenomenon. By detecting fluctuation of interval time between the extrema, whether there is high frequency intermittent component in the signal can be judged. In this paper, variance of interval time is used to judge whether there is intermittent component.LetM(·) be the operator which produces the local mean of the signal that is applied to. · indicates the action of averaging operation for I times instead of (1/I)∑i=1I/2(·). ω(i±)(t) is the realizations of additive white noise mode in pairs for the ith time.We propose the SN-EMD algorithm as follows.Step 1.
Setk=0 and rk(t)=x(t).Step 2.
Judge whether high frequency intermittent component is contained inrk(t).(i)
Yes. Calculate the local means ofrk(i±)(t)=rk(t)+βk[ωk+1(i±)(t)] by EMD to obtain the residue rk+1(t)=Mrk(i±)(t) and go to Step 3.(ii)
No. Get thek+1th mode imfk+1⁡(t) and the residue rk+1(t)=rk(t)-imfk+1(t) by EMD and go to Step 4.Step 3.
Get thek+1th mode imfk+1⁡(t)~=rk(t)-rk+1(t).Step 4.
Go toStep 2 and continue for nextk until residue satisfies predefined stopping criterion.Step 5.
Save all the modes fromSteps 2 and 3. The original signal x(t) can be expressed as a sum of imfn⁡(t),(n=1,…,N) and rN(t).Coefficientβk(βk=εkstd(rk(t))/std(ωk+1(i±)(t))) controls the amplitude of added noise. βk is chosen to obtain a desired SNR between the added noise and the residue to which the noise is added. εk is a constant not greater than 0.2. Figure 1 gives the flowchart of SN-EMD.Figure 1
SN-EMD algorithm flowchart.


## 2.1. Existing Methods

### 2.1.1. EMD

The EMD algorithm can be described as follows.Step 1.
Setk=0 and rk(t)=x(t).Step 2.
Find all the maxima (minima) ofrk(t). Interpolate between maxima (minima) of rk(t) to obtain the upper (lower) envelope emax(t)(emin(t)).Step 3.
Compute the mean envelopem(t)=(emax(t)+emin(t))/2.Step 4.
Compute the IMF candidateck+1(t)=rk(t)-m(t).Step 5.
Judge whetherck+1(t) is an IMF.(i)
Yes. Setimfk+1⁡(t)=ck+1(t) and save imfk+1⁡(t) as IMF of order k+1. Compute the residue rk+1(t)=x(t)-∑n=1k+1imfk+1⁡(t), do k=k+1, and treat rk(t) as input data in Step 2.(ii)
No. Treatck+1(t) as input data in Step 2.Step 6.
Continue until the final residuerk(t) satisfies some predefined stopping criterion.The original signalx(t) can be expressed as a sum of imfn⁡(t),(n=1,…,N).There is mode mixing problem in EMD decomposition, which makes it impossible for EMD to achieve separation of components of the signal at different frequency. The mode obtained by decomposition cannot reflect the instantaneous frequency characteristics of the signal and fundamentally loses the original intention of EMD algorithm. Meanwhile, the deficiency limits the engineering application of the algorithm.


### 2.1.2. EEMD

In order to solve the problem of mode mixing, Wu, Z. [2] proposed the EEMD algorithm. EEMD defines the “true” modes as the average of the corresponding IMFs obtained from an ensemble of the original signal plus different realizations of finite variance white noise. On the basis of the characteristics of the EMD’s binary filter bank, it reduces mode mixing by filling the entire time-frequency space, with adding white noise to the signal.EEMD algorithm also has some drawbacks: (I) Different realizations of finite variance white noise may produce different numbers of modes; (II) the average mode may not meet the requirements of effective IMF; (III) when IMFs are used to reconstruct the signal, the residual noise cannot be ignored; (IV) ensemble average in EEMD usually needs operation for several hundred times, which is very time-consuming.


### 2.1.3. CEEMD

Yeh, J.R. [3] proposed the CEEMD algorithm, which is different with EEMD in that the noise added to the original data is in pairs (a positive number and a negative number). Adding noise with a pair of opposite signs makes the residual noise basically cancel out each other and achieve the effect of “complementary.”


### 2.1.4. CEEMDAN

CEEMDAN proposed in [4] is an important improvement to EEMD. The first-order mode using CEEMDAN algorithm is the same as that using EEMD, and other signal modes are the first-order mode decomposed from the residual signal plus white noise mode in the same order. Also, the coefficients of each noise mode are all adaptive, which are related to the residual signal and noise mode. The improvement of CEEMDAN reduces the reconstruction error and solves the problem that different realizations of signal plus noise might produce different number of modes.


### 2.1.5. Improved CEEMDAN

In order to reduce the residual noise, Colominas, M.A. [5] made an improvement on CEEMDAN: estimating local mean instead mode. From the first-order mode, the noise added in the signal is mode component decomposed from white noise rather than the original white noise.


## 2.1.1. EMD

The EMD algorithm can be described as follows.Step 1.
Setk=0 and rk(t)=x(t).Step 2.
Find all the maxima (minima) ofrk(t). Interpolate between maxima (minima) of rk(t) to obtain the upper (lower) envelope emax(t)(emin(t)).Step 3.
Compute the mean envelopem(t)=(emax(t)+emin(t))/2.Step 4.
Compute the IMF candidateck+1(t)=rk(t)-m(t).Step 5.
Judge whetherck+1(t) is an IMF.(i)
Yes. Setimfk+1⁡(t)=ck+1(t) and save imfk+1⁡(t) as IMF of order k+1. Compute the residue rk+1(t)=x(t)-∑n=1k+1imfk+1⁡(t), do k=k+1, and treat rk(t) as input data in Step 2.(ii)
No. Treatck+1(t) as input data in Step 2.Step 6.
Continue until the final residuerk(t) satisfies some predefined stopping criterion.The original signalx(t) can be expressed as a sum of imfn⁡(t),(n=1,…,N).There is mode mixing problem in EMD decomposition, which makes it impossible for EMD to achieve separation of components of the signal at different frequency. The mode obtained by decomposition cannot reflect the instantaneous frequency characteristics of the signal and fundamentally loses the original intention of EMD algorithm. Meanwhile, the deficiency limits the engineering application of the algorithm.


## 2.1.2. EEMD

In order to solve the problem of mode mixing, Wu, Z. [2] proposed the EEMD algorithm. EEMD defines the “true” modes as the average of the corresponding IMFs obtained from an ensemble of the original signal plus different realizations of finite variance white noise. On the basis of the characteristics of the EMD’s binary filter bank, it reduces mode mixing by filling the entire time-frequency space, with adding white noise to the signal.EEMD algorithm also has some drawbacks: (I) Different realizations of finite variance white noise may produce different numbers of modes; (II) the average mode may not meet the requirements of effective IMF; (III) when IMFs are used to reconstruct the signal, the residual noise cannot be ignored; (IV) ensemble average in EEMD usually needs operation for several hundred times, which is very time-consuming.


## 2.1.3. CEEMD

Yeh, J.R. [3] proposed the CEEMD algorithm, which is different with EEMD in that the noise added to the original data is in pairs (a positive number and a negative number). Adding noise with a pair of opposite signs makes the residual noise basically cancel out each other and achieve the effect of “complementary.”


## 2.1.4. CEEMDAN

CEEMDAN proposed in [4] is an important improvement to EEMD. The first-order mode using CEEMDAN algorithm is the same as that using EEMD, and other signal modes are the first-order mode decomposed from the residual signal plus white noise mode in the same order. Also, the coefficients of each noise mode are all adaptive, which are related to the residual signal and noise mode. The improvement of CEEMDAN reduces the reconstruction error and solves the problem that different realizations of signal plus noise might produce different number of modes.


## 2.1.5. Improved CEEMDAN

In order to reduce the residual noise, Colominas, M.A. [5] made an improvement on CEEMDAN: estimating local mean instead mode. From the first-order mode, the noise added in the signal is mode component decomposed from white noise rather than the original white noise.


## 2.2. SN-EMD

When there are high frequency intermittent components contained in the signal, the EMD will lead to the first type of mode mixing problem that signal components on different time scales in one mode. In order to solve the problem, a series of noise-assisted versions based on EMD have been proposed. However, due to the randomness of the added noise, the noise-assisted algorithms may cause the second type of mode mixing phenomenon that the signal components on the same time scale are decomposed into different modes. In addition, these algorithms also have problems such as high time-consuming and large residual error. For the reasons, a new method SN-EMD (Selective Noise-assisted EMD) is proposed in this paper. Therefore, in order to avoid the phenomenon of mode mixing, it is necessary to judge whether high frequency intermittent component is contained before decomposition and then decide whether to add noise for decomposition or not.At the high frequency intermittent component of the signal, the interval time between the extrema of the signal envelope will decrease obviously. That is, at the beginning and the end of the intermittent component, the interval time between the extrema has a jump phenomenon. By detecting fluctuation of interval time between the extrema, whether there is high frequency intermittent component in the signal can be judged. In this paper, variance of interval time is used to judge whether there is intermittent component.LetM(·) be the operator which produces the local mean of the signal that is applied to. · indicates the action of averaging operation for I times instead of (1/I)∑i=1I/2(·). ω(i±)(t) is the realizations of additive white noise mode in pairs for the ith time.We propose the SN-EMD algorithm as follows.Step 1.
Setk=0 and rk(t)=x(t).Step 2.
Judge whether high frequency intermittent component is contained inrk(t).(i)
Yes. Calculate the local means ofrk(i±)(t)=rk(t)+βk[ωk+1(i±)(t)] by EMD to obtain the residue rk+1(t)=Mrk(i±)(t) and go to Step 3.(ii)
No. Get thek+1th mode imfk+1⁡(t) and the residue rk+1(t)=rk(t)-imfk+1(t) by EMD and go to Step 4.Step 3.
Get thek+1th mode imfk+1⁡(t)~=rk(t)-rk+1(t).Step 4.
Go toStep 2 and continue for nextk until residue satisfies predefined stopping criterion.Step 5.
Save all the modes fromSteps 2 and 3. The original signal x(t) can be expressed as a sum of imfn⁡(t),(n=1,…,N) and rN(t).Coefficientβk(βk=εkstd(rk(t))/std(ωk+1(i±)(t))) controls the amplitude of added noise. βk is chosen to obtain a desired SNR between the added noise and the residue to which the noise is added. εk is a constant not greater than 0.2. Figure 1 gives the flowchart of SN-EMD.Figure 1
SN-EMD algorithm flowchart.


## 3. Simulation and Performance Analysis

In this section, two artificial signals are decomposed by several decomposition algorithms in Section2. The performance of these algorithms is compared and analyzed based on the evaluating indicator given.


### 3.1. Quantitative Evaluation Indexes

In order to accurately describe the performance of the mode decomposition algorithm, five quantitative evaluation indexes for mode separation are introduced: Consuming Time (CT), Index of Orthogonality (IO), Similarity Index (SI), and Residual Error (RE).I O describes the orthogonality between different modes, which can examine whether there is mode mixing between every two modes. The ideal IMF components are completely orthogonal, while IO is equal to zero. Conversely, the larger the IO, the more severe the mode mixing of the two IMF components. In order to ensure IO positive, the definition given in [19] is modified to(1)IO=∑j=1N∑i=1N∑l=1Limfi⁡limfj⁡lx2li m f i ( l ) and imfj(l) are, respectively, the discrete sequences of theith and thejth mode (i≠j), and x(l) is in the discrete form of x(t). The following signals and modes are mostly represented in discrete form. The length of the signal sequence is L and the highest order of the IMF isN.S I describes the similarity between the each order IMF and the original signal. It is used to measure whether the mode decomposition algorithm truly separates the modes from the original signal. SI is defined as follows:(2)SI=1N∑n=1Ncxl,imfn⁡lThe larger theSI, the more similar the IMF and the original signal. x(l)~ indicates the mean value of the original input signal sequence. The correlation coefficient expression is as follows:(3)cxl,imfn⁡l=∑l=1Lxl-xl~imfn⁡l-imfn⁡l~∑l=1Lxl-xl~2∑l=1Limfn⁡l-imfn⁡l~2R E represents residual error of noise when recovering signal components. Several noise-assisted decomposition algorithms require that the residual noise be reduced by averaging for many times. The residual error when recovering x1(t) from artificial signal x(t) above can be expressed as(4)RE1=∑l=1Limf1⁡l-x1l2∑l=1Limf1⁡l2


### 3.2. Simulation

As an example, we propose here a classical mode mixing example. The input signal is an intermittent high frequency sine component plus a sustained intermediate frequency component and a low frequency component for 1 second. The artificial input signalx(t)=x1(t)+x2(t)+x3(t) is shown in Figure 2. Intermittent high frequency signal is expressed as x1(t)=0.2sin⁡(2π300t)(0.3≤t≤0.35, 0.6≤t≤0.65, 0.8≤t≤0.85); two sustained sine signals are x2(t)=sin⁡(2π65t)(0≤t≤1) and x3(t)=0.3sin⁡(2π15t)(0≤t≤1).Figure 2
Artificial signal.In order to determine whether there is high frequency intermittent component in the input signal, the maximum and minimum values of the signal are marked first, and the time interval between the extrema is calculated. The maximum and minimum intervals are normalized to calculate the respective variances. The average value of the two variances is used as the final judgment parameter. Considering the divergent effect of the boundary in decomposition, the first two extrema and the last two extrema are discarded, so the number of extreme points is usually no less than 6. In this paper, the threshold is calculated from decomposition of white noise. Mean variances of normalized interval time between the extrema of white noise and each residue are calculated. The mean value of the several mean variances multiplied by 1.5 is set to be the final threshold, which is equal to 0.032. When the judgment parameter is greater than 0.032, the signal is considered to contain intermittent component.Figure3 shows the time interval between the extrema of the signal when there is and there is no intermittent component x1(t), respectively. The judgment parameters are, respectively, 0.11 and 0.0019 when there is and there is no intermittent component. The threshold is proven to be effective.Figure 3
Interval time between (a) maxima and (b) minima.
(a) (b)Decompose the artificial signal separately using the six algorithms in Section2 and the IMFs are shown in Figures 4(a)–4(f), with ensemble sizes of I=100 and coefficient εn=0.05.Figure 4
Decomposition of artificial signal by (a) EMD, (b) EEMD, (c) CEEMD, (d) CEEMDAN, (e) improved CEEMDAN, and (f) SN-EMD.
(a) (b) (c) (d) (e) (f)From the comparison of Figures4(a)–4(f), it can be found that, in the first-order mode decomposed by EMD, the intermediate frequency component x2(t) is filled in the gap of high frequency component x1(t). In other words, using the EMD algorithm, the high frequency intermittent component cannot be extracted from the original signal. Frequency mixing appears when extracting low-frequency component x3(t) by EEMD and CEEMD; that is, a time scale signal component is decomposed into different modes. Not only mode mixing phenomenon but also false modes such as the second and fifth IMFs appear in decomposition by CEEMDAN as shown in Figure 4(d). Further, the improved CEEMDAN still has mode mixing problems. However, the SN-EMD algorithm proposed in this paper can accurately decompose the original signal into three components, solving the problem of mode mixing. The evaluation indexes CT, IO, and SI are counted and listed in Table 1. Of these, for the five noise-assisted algorithms, the ensemble average number of timesI is 100.Table 1
Performance statistics of mode composition algorithms.
EMD EEMD CEEMD CEEMDAN Improved CEEMDAN SN-EMD CT/s 0.23 13.49 12.97 36.04 29.42 1.91 IO 0.19 0.11 0.10 0.41 0.31 0.056 SI 0.26 0.32 0.32 0.28 0.34 0.37Comparing the evaluation indexes in Table1, it can be seen that the SN-EMD algorithm has the smallest IO and the largest SI, which shows the new algorithm has a better performance of extracting true modes from the signal. In addition, CTof SN-EMD is much smaller than other noise-assisted version.When using the noise-assisted method, due to added random noise, there is residual noise in each mode, which generally decreases with the increase of the average number of times. The ensemble average number of timesI is taken as 20, 100, 200, and 400, respectively, and we can get the average residual error of 100 calculations. For the reason of mode mixing, the true mode can be seen as the superimposition of mode components in the same time scale, which can be judged based on the orthogonality between two modes. The 3rd to 5th order modes, as decomposed by the EEMD algorithm in Figure 4(b), are in a time scale, which can be superimposed as the mode corresponding to x3(t). REi represents the residual error for recovering xi(t). Figures 5(a)–5(c) show the comparison of the residual error in recovering the three components of the original signal.Figure 5
Residual error when recovering (a)x1(t), (b) x2(t), and (c) x3(t).
(a) (b) (c)As can be seen from Figures5(a)–5(c), the residual noise error by SN-EMD when recovering the three components is less than that by the other four algorithms. Furthermore, only if the ensemble average number of times I is taken as 20, the residual noise error has been quite slow. So SN-EMD can reduce the average number of times and save the operation time. In a word, the new algorithm SN-EMD integrates the advantages of the existing decomposition methods.To research the decomposition with noise, another example is proposed. The input signal is composed of two sine components and noise. The artificial input signalx(t)=x1(t)+x2(t)+n(t)(0≤t≤1) is shown in Figure 6. The high frequency sine signal is expressed as x1(t)=sin⁡(2π65t), another sine signal is x2(t)=sin⁡(2π15t), and n(t) is filtered Gaussian noise by a 500 Hz low-pass filter.Figure 6
Artificial signal with noise.Decompose the artificial signal separately using the six algorithms and the modes are shown in Figures7(a)–7(f), with ensemble sizes of I=100 and coefficient εn=0.05.Figure 7
Decomposition of artificial signal with noise by (a) EMD, (b) EEMD, (c) CEEMD, (d) CEEMDAN, (e) improved CEEMDAN, and (f) SN-EMD.
(a) (b) (c) (d) (e) (f)Obviously, the noisen(t) leads to mode mixing problems marked in Figure 7(a) with red boxes. Furthermore, high frequency component and low frequency component alternately appear in. c2(t). From Figures 7(b)–7(e), it can be found that mode mixing, false mode extracting, and other problems cannot be avoided when using the four existing noise-assisted decomposition algorithms. However, the SN-EMD algorithm solves these problems as shown in Figure 7(f). The mean variances of normalized interval time between the maxima and minima of x(t) and ri(t)(i=1,…,4, for the number of maxima of r5(t) being less than 6) are 0.017, 0.076, 0.0071, 0.0017, and 0.024, respectively. Only 0.076 is greater than the threshold 0.032, so high frequency intermittent component firstly appears in r1(t). The interval time between the extrema of r1(t) is shown in Figure 8.Figure 8
Interval time between (a) maxima and (b) minima.
(a) (b)Comparing the evaluation indexes shown in Table2, it can be seen that the SN-EMD algorithm has a better performance of extracting true modes from the signal. In addition, the noise degrades the decomposition performance of each algorithm.Table 2
Performance statistics of mode composition algorithms.
EMD EEMD CEEMD CEEMDAN Improved CEEMDAN SN-EMD CT/s 0.29 15.47 15.14 38. 40 33.67 1.88 IO 0.38 0.31 0.27 0.64 0.44 0.011 SI 0.24 0.28 0.29 0.22 0.32 0.35The ensemble average number of timesI is taken as 20, 100, 200, and 400, respectively, and we can get the average residual error of 100 calculations. Figure 9 shows the comparison of the residual error when recovering the two signal components x1(t) and x2(t). In Figure 8, it can be seen that the residual noise error by SN-EMD when recovering the three components is less than that by the other four algorithms.Figure 9
Residual error when recovering (a)x1(t) and (b) x2(t).
(a) (b)


## 3.1. Quantitative Evaluation Indexes

In order to accurately describe the performance of the mode decomposition algorithm, five quantitative evaluation indexes for mode separation are introduced: Consuming Time (CT), Index of Orthogonality (IO), Similarity Index (SI), and Residual Error (RE).I O describes the orthogonality between different modes, which can examine whether there is mode mixing between every two modes. The ideal IMF components are completely orthogonal, while IO is equal to zero. Conversely, the larger the IO, the more severe the mode mixing of the two IMF components. In order to ensure IO positive, the definition given in [19] is modified to(1)IO=∑j=1N∑i=1N∑l=1Limfi⁡limfj⁡lx2li m f i ( l ) and imfj(l) are, respectively, the discrete sequences of theith and thejth mode (i≠j), and x(l) is in the discrete form of x(t). The following signals and modes are mostly represented in discrete form. The length of the signal sequence is L and the highest order of the IMF isN.S I describes the similarity between the each order IMF and the original signal. It is used to measure whether the mode decomposition algorithm truly separates the modes from the original signal. SI is defined as follows:(2)SI=1N∑n=1Ncxl,imfn⁡lThe larger theSI, the more similar the IMF and the original signal. x(l)~ indicates the mean value of the original input signal sequence. The correlation coefficient expression is as follows:(3)cxl,imfn⁡l=∑l=1Lxl-xl~imfn⁡l-imfn⁡l~∑l=1Lxl-xl~2∑l=1Limfn⁡l-imfn⁡l~2R E represents residual error of noise when recovering signal components. Several noise-assisted decomposition algorithms require that the residual noise be reduced by averaging for many times. The residual error when recovering x1(t) from artificial signal x(t) above can be expressed as(4)RE1=∑l=1Limf1⁡l-x1l2∑l=1Limf1⁡l2


## 3.2. Simulation

As an example, we propose here a classical mode mixing example. The input signal is an intermittent high frequency sine component plus a sustained intermediate frequency component and a low frequency component for 1 second. The artificial input signalx(t)=x1(t)+x2(t)+x3(t) is shown in Figure 2. Intermittent high frequency signal is expressed as x1(t)=0.2sin⁡(2π300t)(0.3≤t≤0.35, 0.6≤t≤0.65, 0.8≤t≤0.85); two sustained sine signals are x2(t)=sin⁡(2π65t)(0≤t≤1) and x3(t)=0.3sin⁡(2π15t)(0≤t≤1).Figure 2
Artificial signal.In order to determine whether there is high frequency intermittent component in the input signal, the maximum and minimum values of the signal are marked first, and the time interval between the extrema is calculated. The maximum and minimum intervals are normalized to calculate the respective variances. The average value of the two variances is used as the final judgment parameter. Considering the divergent effect of the boundary in decomposition, the first two extrema and the last two extrema are discarded, so the number of extreme points is usually no less than 6. In this paper, the threshold is calculated from decomposition of white noise. Mean variances of normalized interval time between the extrema of white noise and each residue are calculated. The mean value of the several mean variances multiplied by 1.5 is set to be the final threshold, which is equal to 0.032. When the judgment parameter is greater than 0.032, the signal is considered to contain intermittent component.Figure3 shows the time interval between the extrema of the signal when there is and there is no intermittent component x1(t), respectively. The judgment parameters are, respectively, 0.11 and 0.0019 when there is and there is no intermittent component. The threshold is proven to be effective.Figure 3
Interval time between (a) maxima and (b) minima.
(a) (b)Decompose the artificial signal separately using the six algorithms in Section2 and the IMFs are shown in Figures 4(a)–4(f), with ensemble sizes of I=100 and coefficient εn=0.05.Figure 4
Decomposition of artificial signal by (a) EMD, (b) EEMD, (c) CEEMD, (d) CEEMDAN, (e) improved CEEMDAN, and (f) SN-EMD.
(a) (b) (c) (d) (e) (f)From the comparison of Figures4(a)–4(f), it can be found that, in the first-order mode decomposed by EMD, the intermediate frequency component x2(t) is filled in the gap of high frequency component x1(t). In other words, using the EMD algorithm, the high frequency intermittent component cannot be extracted from the original signal. Frequency mixing appears when extracting low-frequency component x3(t) by EEMD and CEEMD; that is, a time scale signal component is decomposed into different modes. Not only mode mixing phenomenon but also false modes such as the second and fifth IMFs appear in decomposition by CEEMDAN as shown in Figure 4(d). Further, the improved CEEMDAN still has mode mixing problems. However, the SN-EMD algorithm proposed in this paper can accurately decompose the original signal into three components, solving the problem of mode mixing. The evaluation indexes CT, IO, and SI are counted and listed in Table 1. Of these, for the five noise-assisted algorithms, the ensemble average number of timesI is 100.Table 1
Performance statistics of mode composition algorithms.
EMD EEMD CEEMD CEEMDAN Improved CEEMDAN SN-EMD CT/s 0.23 13.49 12.97 36.04 29.42 1.91 IO 0.19 0.11 0.10 0.41 0.31 0.056 SI 0.26 0.32 0.32 0.28 0.34 0.37Comparing the evaluation indexes in Table1, it can be seen that the SN-EMD algorithm has the smallest IO and the largest SI, which shows the new algorithm has a better performance of extracting true modes from the signal. In addition, CTof SN-EMD is much smaller than other noise-assisted version.When using the noise-assisted method, due to added random noise, there is residual noise in each mode, which generally decreases with the increase of the average number of times. The ensemble average number of timesI is taken as 20, 100, 200, and 400, respectively, and we can get the average residual error of 100 calculations. For the reason of mode mixing, the true mode can be seen as the superimposition of mode components in the same time scale, which can be judged based on the orthogonality between two modes. The 3rd to 5th order modes, as decomposed by the EEMD algorithm in Figure 4(b), are in a time scale, which can be superimposed as the mode corresponding to x3(t). REi represents the residual error for recovering xi(t). Figures 5(a)–5(c) show the comparison of the residual error in recovering the three components of the original signal.Figure 5
Residual error when recovering (a)x1(t), (b) x2(t), and (c) x3(t).
(a) (b) (c)As can be seen from Figures5(a)–5(c), the residual noise error by SN-EMD when recovering the three components is less than that by the other four algorithms. Furthermore, only if the ensemble average number of times I is taken as 20, the residual noise error has been quite slow. So SN-EMD can reduce the average number of times and save the operation time. In a word, the new algorithm SN-EMD integrates the advantages of the existing decomposition methods.To research the decomposition with noise, another example is proposed. The input signal is composed of two sine components and noise. The artificial input signalx(t)=x1(t)+x2(t)+n(t)(0≤t≤1) is shown in Figure 6. The high frequency sine signal is expressed as x1(t)=sin⁡(2π65t), another sine signal is x2(t)=sin⁡(2π15t), and n(t) is filtered Gaussian noise by a 500 Hz low-pass filter.Figure 6
Artificial signal with noise.Decompose the artificial signal separately using the six algorithms and the modes are shown in Figures7(a)–7(f), with ensemble sizes of I=100 and coefficient εn=0.05.Figure 7
Decomposition of artificial signal with noise by (a) EMD, (b) EEMD, (c) CEEMD, (d) CEEMDAN, (e) improved CEEMDAN, and (f) SN-EMD.
(a) (b) (c) (d) (e) (f)Obviously, the noisen(t) leads to mode mixing problems marked in Figure 7(a) with red boxes. Furthermore, high frequency component and low frequency component alternately appear in. c2(t). From Figures 7(b)–7(e), it can be found that mode mixing, false mode extracting, and other problems cannot be avoided when using the four existing noise-assisted decomposition algorithms. However, the SN-EMD algorithm solves these problems as shown in Figure 7(f). The mean variances of normalized interval time between the maxima and minima of x(t) and ri(t)(i=1,…,4, for the number of maxima of r5(t) being less than 6) are 0.017, 0.076, 0.0071, 0.0017, and 0.024, respectively. Only 0.076 is greater than the threshold 0.032, so high frequency intermittent component firstly appears in r1(t). The interval time between the extrema of r1(t) is shown in Figure 8.Figure 8
Interval time between (a) maxima and (b) minima.
(a) (b)Comparing the evaluation indexes shown in Table2, it can be seen that the SN-EMD algorithm has a better performance of extracting true modes from the signal. In addition, the noise degrades the decomposition performance of each algorithm.Table 2
Performance statistics of mode composition algorithms.
EMD EEMD CEEMD CEEMDAN Improved CEEMDAN SN-EMD CT/s 0.29 15.47 15.14 38. 40 33.67 1.88 IO 0.38 0.31 0.27 0.64 0.44 0.011 SI 0.24 0.28 0.29 0.22 0.32 0.35The ensemble average number of timesI is taken as 20, 100, 200, and 400, respectively, and we can get the average residual error of 100 calculations. Figure 9 shows the comparison of the residual error when recovering the two signal components x1(t) and x2(t). In Figure 8, it can be seen that the residual noise error by SN-EMD when recovering the three components is less than that by the other four algorithms.Figure 9
Residual error when recovering (a)x1(t) and (b) x2(t).
(a) (b)


## 4. Experimental Results and Discussion

The mode decomposition of the ship radiated noise is to adaptively divide the original signal according to different time scales. The original information contained in the signal is also assigned to each mode when decomposing the signal. For different types of ships, the distribution of the original information contained in each mode is different. Based on these differences, time-frequency features can be extracted for each mode as a standard for classification and recognition. The flowchart of mode feature extraction from ship radiated noise based on SN-EMD is shown in Figure10.Figure 10
Flowchart of mode feature extraction.


### 4.1. Mode Feature

#### 4.1.1. Average Instantaneous Frequency

The average instantaneous frequency can exhibit the frequency distribution of a mode. The instantaneous frequencyfij at the moment of j can be achieved by taking the derivative of the instantaneous phase of theith mode obtained with Hilbert transform. Set K to be the number of sample points per unit time. So, the average instantaneous frequency per unit time is(5)fi¯=1K∑j=1Kfij


#### 4.1.2. Center Frequency

The center frequency is the average instantaneous frequency weighted by the instantaneous amplitude per unit time. The instantaneous amplitudeAij at the moment of j is calculated after Hilbert transforming the ith mode. Then the center frequency is defined as(6)fmi¯=∑j=1KAijfij∑j=1KAij


#### 4.1.3. Energy Density

The energy density reflects the average energy of the mode per unit time. The energy density is defined as(7)Ei¯=1K∑j=1KAij2


#### 4.1.4. Energy Distribution Ratio

The energy distribution ratioEpi describes the energy share of each mode. The ratio of the energy of the ith mode to the energy of the original signal is defined as the energy distribution ratio with the expression:(8)Epi=∑j=1KAij2∑i=1N∑j=1KAij2


### 4.2. Mode Feature Extraction

Two sets of ship radiated noise data are collected from a large ship named “Target A,” and other two sets are collected from a boat and a yacht, which are named “Target B” and “Target C,” respectively. Figure11 shows the first 8 modes decomposed by SN-EMD from the four sets of data of the three targets after 5 kHz low-pass filtering and normalization. The sampling rate is 48 kHz and the data length is 9000 (0.1875 seconds). The average instantaneous frequency, center frequency, energy density, and energy distribution ratio are extracted from the modes, whose separability is verified in three-dimensional vector space.Figure 11
Decomposition of ship radiated noise signal (a) Target A1; (b) Target A2; (c) Target B; (d) Target C.
(a) (b) (c) (d)The mean variances of normalized interval time between the maxima and minima of Target A1’s signal and first 7 residues are 0.019, 0.022, 0.066, 0.019, 0.018, 0.018, 0.024, and 0.027. So only the second-order residue contains intermittent component. The mean variances of other targets will not be listed here.It can be seen that the amplitude of each mode is greatly reduced after the fourth-order mode in Figure11, so the energy is mainly concentrated in the first four modes. Instantaneous phase curves over time of the first eight modes (c1~c8) of Target A’ radiated noise are shown in Figure 12. Calculate the derivative of the instantaneous phase and divide it by 2π, so as to get instantaneous frequency. The instantaneous phase of the first four modes almost linearly increases, indicating that the instantaneous frequency is relatively stable. Therefore, only the first four modes are processed below to extract feature quantities, which are shown by box-plot.Figure 12
Instantaneous phase of first 8 modes of Target A1’s radiated noise.
(a) (b)Instantaneous frequency and instantaneous amplitude are obtained for the first 4 modes. 30 average instantaneous frequencies and 30 central frequencies calculated by grouping and weighting average are shown in Figures13 and 14. The horizontal axis coordinates in the figures correspond to the mode orders.Figure 13
Statistics for average instantaneous frequency of first 4 modes.Figure 14
Statistics for center frequency of first 4 modes.From Figures13 and 14, it can be found that the average instantaneous frequency and center frequency of the same mode have no much difference. Target A1 has the mode frequency commensurate with Target A2. Comparing statistics of different targets, it can be seen that the first-order mode frequency of target B is higher than that of Target A and Target C about 400-500 Hz, and the frequency of the second- to the fourth-order modes of target C is higher than that of Target A and Target B. That is, mode frequency of different targets has different distribution.The energy density and energy distribution ratio of the first four modes of the four sets of data are counted as shown in Figures15 and 16. The statistical results indicate that the first three modes’ energies of Target A are not significantly different but are all higher than that of the fourth-order mode, while the mode energy of Target B is mainly concentrated in the first-order mode, with Target C’s mode energy gradually decreasing. Among the three targets, the energy distribution of Target A is the most uniform, and the distribution of energy of Target B is the most uneven.Figure 15
Statistics for energy density of first 4 modes.Figure 16
Statistics for energy distribution ratio of first 4 modes.


### 4.3. Separability Validation

Classifying and recognizing ship targets require that the mode feature of different targets should be distinguishable, while the mode feature of the same target is similar. In order to verify the separability of average instantaneous frequency, center frequency, energy density, and energy distribution ratio feature, feature quantities of the first 3 modes extracted from the four sets of data are displayed in the three-dimensional space. As shown in Figure17, the three-dimensional coordinates, respectively, stand for the feature quantities of first three modes.Figure 17
Feature quantities of first 3 modes distributed in three-dimensional space: (a) average instantaneous frequency; (b) center frequency; (c) energy density; (d) energy distribution ratio.
(a) (b) (c) (d)The feature quantities extracted from the two sets of data of Target A present mixed distribution in the same three-dimensional space range, while the feature quantities of the three targets mostly have different spatial distributions. This shows that the average instantaneous frequency, center frequency, energy density, and energy distribution ratio are easy to separate, and these feature quantities provide a basis for ship target classification and identification.


## 4.1. Mode Feature

### 4.1.1. Average Instantaneous Frequency

The average instantaneous frequency can exhibit the frequency distribution of a mode. The instantaneous frequencyfij at the moment of j can be achieved by taking the derivative of the instantaneous phase of theith mode obtained with Hilbert transform. Set K to be the number of sample points per unit time. So, the average instantaneous frequency per unit time is(5)fi¯=1K∑j=1Kfij


### 4.1.2. Center Frequency

The center frequency is the average instantaneous frequency weighted by the instantaneous amplitude per unit time. The instantaneous amplitudeAij at the moment of j is calculated after Hilbert transforming the ith mode. Then the center frequency is defined as(6)fmi¯=∑j=1KAijfij∑j=1KAij


### 4.1.3. Energy Density

The energy density reflects the average energy of the mode per unit time. The energy density is defined as(7)Ei¯=1K∑j=1KAij2


### 4.1.4. Energy Distribution Ratio

The energy distribution ratioEpi describes the energy share of each mode. The ratio of the energy of the ith mode to the energy of the original signal is defined as the energy distribution ratio with the expression:(8)Epi=∑j=1KAij2∑i=1N∑j=1KAij2


## 4.1.1. Average Instantaneous Frequency

The average instantaneous frequency can exhibit the frequency distribution of a mode. The instantaneous frequencyfij at the moment of j can be achieved by taking the derivative of the instantaneous phase of theith mode obtained with Hilbert transform. Set K to be the number of sample points per unit time. So, the average instantaneous frequency per unit time is(5)fi¯=1K∑j=1Kfij


## 4.1.2. Center Frequency

The center frequency is the average instantaneous frequency weighted by the instantaneous amplitude per unit time. The instantaneous amplitudeAij at the moment of j is calculated after Hilbert transforming the ith mode. Then the center frequency is defined as(6)fmi¯=∑j=1KAijfij∑j=1KAij


## 4.1.3. Energy Density

The energy density reflects the average energy of the mode per unit time. The energy density is defined as(7)Ei¯=1K∑j=1KAij2


## 4.1.4. Energy Distribution Ratio

The energy distribution ratioEpi describes the energy share of each mode. The ratio of the energy of the ith mode to the energy of the original signal is defined as the energy distribution ratio with the expression:(8)Epi=∑j=1KAij2∑i=1N∑j=1KAij2


## 4.2. Mode Feature Extraction

Two sets of ship radiated noise data are collected from a large ship named “Target A,” and other two sets are collected from a boat and a yacht, which are named “Target B” and “Target C,” respectively. Figure11 shows the first 8 modes decomposed by SN-EMD from the four sets of data of the three targets after 5 kHz low-pass filtering and normalization. The sampling rate is 48 kHz and the data length is 9000 (0.1875 seconds). The average instantaneous frequency, center frequency, energy density, and energy distribution ratio are extracted from the modes, whose separability is verified in three-dimensional vector space.Figure 11
Decomposition of ship radiated noise signal (a) Target A1; (b) Target A2; (c) Target B; (d) Target C.
(a) (b) (c) (d)The mean variances of normalized interval time between the maxima and minima of Target A1’s signal and first 7 residues are 0.019, 0.022, 0.066, 0.019, 0.018, 0.018, 0.024, and 0.027. So only the second-order residue contains intermittent component. The mean variances of other targets will not be listed here.It can be seen that the amplitude of each mode is greatly reduced after the fourth-order mode in Figure11, so the energy is mainly concentrated in the first four modes. Instantaneous phase curves over time of the first eight modes (c1~c8) of Target A’ radiated noise are shown in Figure 12. Calculate the derivative of the instantaneous phase and divide it by 2π, so as to get instantaneous frequency. The instantaneous phase of the first four modes almost linearly increases, indicating that the instantaneous frequency is relatively stable. Therefore, only the first four modes are processed below to extract feature quantities, which are shown by box-plot.Figure 12
Instantaneous phase of first 8 modes of Target A1’s radiated noise.
(a) (b)Instantaneous frequency and instantaneous amplitude are obtained for the first 4 modes. 30 average instantaneous frequencies and 30 central frequencies calculated by grouping and weighting average are shown in Figures13 and 14. The horizontal axis coordinates in the figures correspond to the mode orders.Figure 13
Statistics for average instantaneous frequency of first 4 modes.Figure 14
Statistics for center frequency of first 4 modes.From Figures13 and 14, it can be found that the average instantaneous frequency and center frequency of the same mode have no much difference. Target A1 has the mode frequency commensurate with Target A2. Comparing statistics of different targets, it can be seen that the first-order mode frequency of target B is higher than that of Target A and Target C about 400-500 Hz, and the frequency of the second- to the fourth-order modes of target C is higher than that of Target A and Target B. That is, mode frequency of different targets has different distribution.The energy density and energy distribution ratio of the first four modes of the four sets of data are counted as shown in Figures15 and 16. The statistical results indicate that the first three modes’ energies of Target A are not significantly different but are all higher than that of the fourth-order mode, while the mode energy of Target B is mainly concentrated in the first-order mode, with Target C’s mode energy gradually decreasing. Among the three targets, the energy distribution of Target A is the most uniform, and the distribution of energy of Target B is the most uneven.Figure 15
Statistics for energy density of first 4 modes.Figure 16
Statistics for energy distribution ratio of first 4 modes.


## 4.3. Separability Validation

Classifying and recognizing ship targets require that the mode feature of different targets should be distinguishable, while the mode feature of the same target is similar. In order to verify the separability of average instantaneous frequency, center frequency, energy density, and energy distribution ratio feature, feature quantities of the first 3 modes extracted from the four sets of data are displayed in the three-dimensional space. As shown in Figure17, the three-dimensional coordinates, respectively, stand for the feature quantities of first three modes.Figure 17
Feature quantities of first 3 modes distributed in three-dimensional space: (a) average instantaneous frequency; (b) center frequency; (c) energy density; (d) energy distribution ratio.
(a) (b) (c) (d)The feature quantities extracted from the two sets of data of Target A present mixed distribution in the same three-dimensional space range, while the feature quantities of the three targets mostly have different spatial distributions. This shows that the average instantaneous frequency, center frequency, energy density, and energy distribution ratio are easy to separate, and these feature quantities provide a basis for ship target classification and identification.


## 5. Conclusions

Based on the adaptive mode decomposition, the mode decomposition of ship radiated noise and mode feature extraction have been studied in this paper. Firstly, several existing mode decomposition algorithms were introduced, and the advantage and disadvantage were studied. Secondly, the SN-EMD algorithm was proposed to overcome the problems of the existing mode decomposition algorithms. Furthermore, the method of judging whether the signal contains intermittent component was given. By contrasting the performance parameters of these algorithms, the results showed that the new algorithm has better performance on resisting mode mixing and time-consuming than other noise-assisted algorithms, and residual error when recovering signal components and reconstruction error are also relatively minimal. Finally, based on the differences in the original information contained in each mode of different ships’ radiated noise, the first 3 modes were selected to extract the average instantaneous frequency, center frequency, energy density, and energy distribution ratio. In the three-dimensional space, distribution of feature quantities verifies that features of the same target are similar and the features of different targets are separable. As a result, the feature quantities can be seen as the basis for ship classification and identification.



---

*Source: 2184612-2018-12-03.xml*