Improved Empirical Wavelet Transform for Compound Weak Bearing Fault Diagnosis with Acoustic Signals

Most of the current research on the diagnosis of rolling bearing faults is based on vibration signals. However, the location and number of sensors are often limited in some special cases. Thus, a small number of non-contact microphone sensors are a suboptimal choice, but it will result in some problems, e.g., underdetermined compound fault detection from a low signal-to-noise ratio (SNR) acoustic signal. Empirical wavelet transform (EWT) is a signal processing algorithm that has a dimension-increasing characteristic, and is beneficial for solving the underdetermined problem with few microphone sensors. However, there remain some critical problems to be solved for EWT, especially the determination of signal mode numbers, high-frequency modulation and boundary detection. To solve these problems, this paper proposes an improved empirical wavelet transform strategy for compound weak bearing fault diagnosis with acoustic signals. First, a novel envelope demodulation-based EWT (DEWT) is developed to overcome the high frequency modulation, based on which a source number estimation method with singular value decomposition (SVD) is then presented for the extraction of the correct boundary from a low SNR acoustic signal. Finally, the new fault diagnosis scheme that utilizes DEWT and SVD is compared with traditional methods, and the advantages of the proposed method in weak bearing compound fault diagnosis with a single-channel, low SNR, variable speed acoustic signal, are verified.


Introduction
Rotating machinery is a type of equipment that is widely used in various industries, the safety and reliable operation of which have become increasingly important [1]. As a significant component of rotating machinery, rolling bearings often operate under very harsh environments, such as high temperatures, high speeds and heavy loads, and are thus prone to failure [2]. Therefore, it is of great significance to make timely and accurate diagnoses of rolling bearing failures [3]. In most cases, rolling bearing faults are often compound because of the mutual contact between various components, and multiple fault features are superimposed and interfere with each other, thereby complicating the fault features and creating difficulties in bearing fault diagnosis [4,5]. Bearing fault diagnosis is often based upon vibration signals collected by acceleration sensors, but the sensor positions and numbers are limited in some special cases, e.g., nuclear power equipment monitoring [6].
When acoustic sensors can be positioned in a non-intrusive and non-contacting way, the use of a small number of microphone sensors to collect device status signals is an appropriate option for such special cases. However, acoustic signals are easily interfered with environmental noise, which makes it difficult to diagnose compound weak fault features [7].
In recent decades, researchers have proposed a number of methods for weak bearing fault diagnosis, a popular one among which is the adaptive signal decomposition method [8]. Empirical wavelet transform (EWT) is a commonly used method that mainly solves already-existing problems in the widely-used empirical mode decomposition (EMD) method, e.g., mode aliasing, endpoint effects, and unclear physical meaning [9,10]. EWT first segments the Fourier spectrum of a signal, and then constructs a series of empirical wavelet-based filter banks to extract the intrinsic mode of the original signal. EWT has been extensively used in many fields, including bearing fault diagnosis [11][12][13]. Li [14] applied EWT to mechanical fault diagnosis. Li [15] proposed a bearing feature extraction and classification method based on EWT and compared it with ensemble empirical mode decomposition (EEMD), whose findings indicated that EWT has more adaptive feature extraction and fault classification characteristics than EEMD. Chen [16] made a fault diagnosis of a generator bearing for a wind turbine by using EWT with vibration signals. Kedadouche [17] applied EWT to the bearing fault diagnosis by combining EWT and operational modal analysis (OMA), which improved the detection accuracy of fault diagnosis. Yan [18] carried out research on bearing fault diagnosis based on order analysis and EWT, revealing that the combination of order analysis and EWT can be used to accurately identify bearing fault features under non-steady working conditions. However, EWT still has some problems when dealing with compound weak bearing faults that must be solved [19]. If the input signal is composed of two chirps that overlap in both time and frequency domains, EWT will not be able to separate them. Additionally, the adaptive detection of the number of signal modes is also a problem that must be further studied [9]. EWT boundary detection is based upon the Fourier spectrum, and because the characteristics of weak faults are not obvious, and the noise of the acoustic signal is relatively high, it is difficult to effectively identify the fault impact. In addition, due to the mutual influence of bearing faults, multiple fault features are superimposed and interfere with each other, thereby complicating the fault features and bringing difficulties to bearing fault diagnosis [20]. The bearing signal also exhibits a modulation phenomenon, and thus the bearing fault frequency cannot be accurately detected in the Fourier spectrum, especially in the case of compound faults [21]. Hence, boundary detection is not well performed in the Fourier spectrum. Therefore, EWT is unable to correctly determine the number of separation modes, and thus cannot effectively separate the intrinsic modes from the measured bearing fault signal.
To solve these problems, an improved EWT strategy is proposed in this paper. First, a novel envelope demodulation-based EWT (DEWT) is proposed. The signal is decomposed by DEWT, each component is treated as a virtual observation signal, and all components constitute a new multi-channel virtual observation signal. The multi-channel virtual observation signal is then input into singular value decomposition (SVD). The mode number is determined by a method of adjacent singular value difference, according to which the number of DEWT components is selected, and DEWT is performed once more. Finally, the bearing health status is determined via spectrum analysis. The main contributions of this paper are as follows.
(1) An envelope demodulation-based strategy for EWT boundary detection is developed. The high-frequency modulation problem of the bearing vibration signal is solved, and the improved EWT can be better used in the field of rolling bearing signal processing. (2) A source number estimation method based on the DEWT and SVD is proposed to determine the number of decomposed modes, which solves the problem of compound bearing fault diagnosis with single-channel acoustic signals.
The remainder of this paper is organized as follows. The basic theory of the proposed method is briefly introduced in Section 2. Section 3 describes the proposed method. Experimental validation and method comparison are presented in Section 4. Finally, the conclusions are given in Section 5.

Empirical Wavelet Transform
The empirical wavelet transform (EWT) divides the Fourier spectrum of a signal and constructs a wavelet filter bank to extract the intrinsic modes, and can be divided into three main components.
1. Fourier spectrum segmentation of a signal Assuming that the Fourier support [0, π] is segmented into contiguous N segments, then a total number of N + 1 boundaries is needed. According to the local maximum values, the Fourier spectrum of a signal is divided into N segments. Denote ω n to be the limits between each segments (where ω 0 = 0 and ω n = π). Each segment is represented as Λ n = [ω n−1 , ω n ]. A transition phase T n of width 2τ n is then defined, as shown in Figure 1.

Empirical Wavelet Transform
The empirical wavelet transform (EWT) divides the Fourier spectrum of a signal and constructs a wavelet filter bank to extract the intrinsic modes, and can be divided into three main components.

Fourier spectrum segmentation of a signal
Assuming that the Fourier support [ ] 0,π is segmented into contiguous N segments, then a total number of 1 N + boundaries is needed. According to the local maximum values, the Fourier spectrum of a signal is divided into N segments. Denote n ω to be the limits between each segments (where 0 =0 ω and n ω π = ). Each segment is represented as The ratio γ in Equations (1) and (2) is restricted to a small value as 0 ≤ γ ≤ min n [(ω n+1 − ω n )/(ω n+1 + ω n )] to ensure the empirical scaling function and the empirical wavelets are a tight frame of L 2 (R). ω n : τ n = γω n , 0 < γ < 1. The function β(x) is an arbitrary C k ([0, 1]) function, the most used β(x) is defined as: 3. Calculation of approximation coefficients and the detail coefficients The approximation coefficients w ε f (0, t) and the detail coefficients w ε f (n, t) can be calculated by Equations (3) and (4), respectively.

Adjacent Singular Value Difference Method
The basic principle of the singular value decomposition (SVD) method is as follows: the covariance matrix R s of n signal sources s(t) is a square matrix of n*n, and the number of linearly independent columns (rows) is the number of uncorrelated sources (NIS) of s(t). In general, the covariance matrix R s cannot be directly obtained (because the source signal s(t) is unknown), but the observation signal covariance matrix R x can be directly obtained. Therefore, the SVD method is conducted through R x to estimate the number of unrelated sources.
Then, a set of adjacent singular values Λ i,j = λ i − λ j , i = 1, 2, · · · , M − 1, j = i + 1 is obtained. When one or both of the λ i and λ j values are singular values of fault signals, Λ i,j is larger. When both λ i and λ j values are singular values of noise signals, Λ i,j is relatively smaller. The difference of the adjacent singular values is detected, the minimum difference corresponds to the interface between the fault signal and noise subspace, and the number of sources can then accordingly be estimated.

Envelope Demodulation-Based Empirical Wavelet Transform
The analytical signal g(t) of the original signal x(t) is obtained by Equations (7) and (8): The amplitude of g(t): A(t) is the envelope of the original real signal x(t).
A Fourier transform is then performed on Equation (7): Appl. Sci. 2020, 10, 682 5 of 16 where F represents a Fourier transform of the function, and sgn() is a symbolic function.
Equation (10) shows that ∧ X( f ) is obtained by phase shifting in the frequency domain by X( f ), delaying π/2 in the positive frequency domain and leading π/2 in the negative frequency domain. The obtained demodulated signal is taken as the input. Based on the empirical scaling function and the empirical wavelet expressions presented in Section 2.1, the approximation coefficients w ε f (0, t) and the detail coefficients w ε f (n, t) are obtained, and are respectively represented by Equations (11) and (12).
where ∧ φ 1 (ω) and ∧ ψ(ω) are defined by Equations (1) and (2), respectively. The schematic diagram of the proposed method is illustrated in Figure 2. We are assuming that the Fourier support [0, π] is segmented into contiguous N segments. The red dashed lines in the figure represent the boundary detected by EWT, and the solid blue lines represent the spectrum contour. Due to the existence of high-frequency modulation, most of the boundaries detected by traditional EWT are located in the bearing resonance frequency band, which leads to modal aliasing in EWT. The proposed method demodulates the compound fault feature information through resonance demodulation. This boundary detection mechanism solves the abovementioned problems well.

Adaptive Estimation of the Modes Number
Because the object signal is measured from a single channel acoustic sensor in this study, the dimension can be increased by DEWT. The specific steps are as follows: Suppose you have the observation signal ( ) x t of the sensor; ( ) x t is decomposed adaptively by DEWT to obtain the intrinsic mode, represented by matrix where c represents the intrinsic mode;

Adaptive Estimation of the Modes Number
Because the object signal is measured from a single channel acoustic sensor in this study, the dimension can be increased by DEWT. The specific steps are as follows: Suppose you have the observation signal x(t) of the sensor; x(t) is decomposed adaptively by DEWT to obtain the intrinsic mode, represented by matrix x i = [c 11 , · · · , c 1n ], where c represents the intrinsic mode; Combine x(t) and x i into a virtual multi-channel sensor observation signal x oi ; Calculate the covariance matrix R i of x oi ; Singular value decomposition is carried out on R i , by which n singular values are obtained. They are then arranged in descending order; The difference of adjacent singular values is obtained, and the number of sources is estimated according to the adjacent singular value difference method.
A schematic diagram of the method proposed in this section is illustrated in Figure 3, in which Λ 1,2 to Λ 6,7 represent the differences between two adjacent singular values. When the difference between two adjacent singular values is as close to 0 as possible, the number of modes is the point abscissa value minus 1. For example, in this diagram, the value of Λ 4,5 is close to 0, and the subsequent values are also close to 0, which means that the number of modes is 3.

Improved EWT Strategy
The method proposed in Section 3.1 improves the boundary detection mechanism of EWT, and a new method for determining the number of modes of EWT was proposed in Section 3.2. Through the improvement of the EWT boundary detection mechanism and mode number estimation method, they can be better applied to the compound weak fault diagnosis of the rolling bearings. Experiments using a variable speed acoustic signal were conducted to verify the effectiveness of the method. The reason that using the variable speed signal is that the early weak fault of the bearing is more easily exposed during the variable speed phase. The acoustic signal was used because an acoustic sensor can simultaneously monitor multiple parts of the device, and can thus overcome the limited number of sensor installations. Because the bearing speed slowly increases linearly, angular resampling was used in this study to conduct preprocessing. The flow chart of the method proposed in this paper is presented in Figure 4.

Improved EWT Strategy
The method proposed in Section 3.1 improves the boundary detection mechanism of EWT, and a new method for determining the number of modes of EWT was proposed in Section 3.2. Through the improvement of the EWT boundary detection mechanism and mode number estimation method, they can be better applied to the compound weak fault diagnosis of the rolling bearings. Experiments using a variable speed acoustic signal were conducted to verify the effectiveness of the method. The reason that using the variable speed signal is that the early weak fault of the bearing is more easily exposed during the variable speed phase. The acoustic signal was used because an acoustic sensor can simultaneously monitor multiple parts of the device, and can thus overcome the limited number of sensor installations. Because the bearing speed slowly increases linearly, angular resampling was used in this study to conduct preprocessing. The flow chart of the method proposed in this paper is presented in Figure 4.

Test Rig
The acoustic signal was collected using a microphone sensor from a rolling bearing fault simulation test rig with an acceleration process. The test rig was driven by an electric motor, and the power was transmitted to the rotating shaft through flexible coupling. The fault bearing was installed at the farthest end of the motor. The locations of the fault bearing and the microphone sensor are shown in Figure 5.

Test Rig
The acoustic signal was collected using a microphone sensor from a rolling bearing fault simulation test rig with an acceleration process. The test rig was driven by an electric motor, and the power was transmitted to the rotating shaft through flexible coupling. The fault bearing was installed at the farthest end of the motor. The locations of the fault bearing and the microphone sensor are shown in Figure 5.
The faulty bearing is depicted in Figure 6. Use wire cutting to preset faults on bearings. The shape of the fault was a through-slot with a width and depth of 0.1 mm. The signal sampling frequency was 20 kHz, and the rotational speed was an acceleration process from 0 to 1328.2 r/min.

The Proposed Method
The time domain diagram of the collected original signal is presented in Figure 7. From the signal waveform, it can be concluded that the bearing state experienced three stages: static, linear acceleration and uniform speed. In the linear acceleration phase, the time-domain waveform was stable without obvious periodic impact. It can also be found that the amplitude of the signal was low, and contained some noise interferences. The linear acceleration signal required in this study was obtained from the collected signal. Through the signal waveform, the 40,000 th to 100,000 th points of the acceleration process were selected, and the instantaneous speed at each time point was calculated by the speed pulse signal collected by the tachometer. The equation of the accelerating speed was fitted as follows: 5.0003 6.1888 V t = × + , where V is the instantaneous speed (Hz) and t is the sampling time point (s). The faulty bearing is depicted in Figure 6. Use wire cutting to preset faults on bearings. The shape of the fault was a through-slot with a width and depth of 0.1 mm. The signal sampling frequency was 20 kHz, and the rotational speed was an acceleration process from 0 to 1328.2 r/min.

Test Rig
The acoustic signal was collected using a microphone sensor from a rolling bearing fault simulation test rig with an acceleration process. The test rig was driven by an electric motor, and the power was transmitted to the rotating shaft through flexible coupling. The fault bearing was installed at the farthest end of the motor. The locations of the fault bearing and the microphone sensor are shown in Figure 5.
The faulty bearing is depicted in Figure 6. Use wire cutting to preset faults on bearings. The shape of the fault was a through-slot with a width and depth of 0.1 mm. The signal sampling frequency was 20 kHz, and the rotational speed was an acceleration process from 0 to 1328.2 r/min.

The Proposed Method
The time domain diagram of the collected original signal is presented in Figure 7. From the signal waveform, it can be concluded that the bearing state experienced three stages: static, linear acceleration and uniform speed. In the linear acceleration phase, the time-domain waveform was stable without obvious periodic impact. It can also be found that the amplitude of the signal was low, and contained some noise interferences. The linear acceleration signal required in this study was obtained from the collected signal. Through the signal waveform, the 40,000 th to 100,000 th points of the acceleration process were selected, and the instantaneous speed at each time point was calculated by the speed pulse signal collected by the tachometer. The equation of the accelerating speed was fitted as follows: 5.0003 6.1888 V t = × + , where V is the instantaneous speed (Hz) and t is the sampling time point (s).

The Proposed Method
The time domain diagram of the collected original signal is presented in Figure 7. From the signal waveform, it can be concluded that the bearing state experienced three stages: static, linear acceleration and uniform speed. In the linear acceleration phase, the time-domain waveform was stable without obvious periodic impact. It can also be found that the amplitude of the signal was low, and contained some noise interferences. The linear acceleration signal required in this study was obtained from the collected signal. Through the signal waveform, the 40,000th to 100,000th points of the acceleration process were selected, and the instantaneous speed at each time point was calculated by the speed pulse signal collected by the tachometer. The equation of the accelerating speed was fitted as follows: V = 5.0003 × t + 6.1888, where V is the instantaneous speed (Hz) and t is the sampling time point (s).    Because the original signal-to-noise ratio (SNR) of the variable speed signal was very low, as shown in Figure 9, the extraction of the fault characteristic order after angular domain resampling was difficult. It can be determined from the envelope order spectrum that the amplitude of the signal was very small. There was only a low-frequency component in the low-frequency part, which does not accurately reflect the state of the bearing. The proposed method was then used to extract the  Table 1.  Table 1.   Because the original signal-to-noise ratio (SNR) of the variable speed signal was very low, as shown in Figure 9, the extraction of the fault characteristic order after angular domain resampling was difficult. It can be determined from the envelope order spectrum that the amplitude of the signal was very small. There was only a low-frequency component in the low-frequency part, which does not accurately reflect the state of the bearing. The proposed method was then used to extract the  Because the original signal-to-noise ratio (SNR) of the variable speed signal was very low, as shown in Figure 9, the extraction of the fault characteristic order after angular domain resampling was difficult. It can be determined from the envelope order spectrum that the amplitude of the signal was very small. There was only a low-frequency component in the low-frequency part, which does not accurately reflect the state of the bearing. The proposed method was then used to extract the bearing fault characteristic order under such complicated conditions, and the superiority of the proposed method was proven.
In the first step, the DEWT-SVD was used to estimate the source number. After the signal was decomposed by DEWT, 10 components were obtained. These 10 components were combined with the original signal to form a virtual multi-channel observation signal. The virtual observation signal was then processed by SVD, and 11 singular values were obtained, as presented in Table 2. The adjacent singular value difference is presented in Table 3. Figure 10 illustrates the broken line diagram of the adjacent singular value difference. It is evident from Tables 2 and 3 and Figure 10 that the singular value change after λ4 was very small, as was the difference between λ4 and λ5, which was about equal to 0. Hence, it can be concluded that the number of source signals of the compound signal was 3, which is consistent with the three fault sources of the inner ring, outer ring, and rolling element fault that had been artificially processed prior. The results of this step indicate that the proposed DEWT-SVD source number estimation method exhibits superior performance.    bearing fault characteristic order under such complicated conditions, and the superiority of the proposed method was proven.
In the first step, the DEWT-SVD was used to estimate the source number. After the signal was decomposed by DEWT, 10 components were obtained. These 10 components were combined with the original signal to form a virtual multi-channel observation signal. The virtual observation signal was then processed by SVD, and 11 singular values were obtained, as presented in Table 2. The adjacent singular value difference is presented in Table 3. Figure 10 illustrates the broken line diagram of the adjacent singular value difference. It is evident from Tables 2 and 3 and Figure 10 that the singular value change after λ4 was very small, as was the difference between λ4 and λ5, which was about equal to 0. Hence, it can be concluded that the number of source signals of the compound signal was 3, which is consistent with the three fault sources of the inner ring, outer ring, and rolling element fault that had been artificially processed prior. The results of this step indicate that the proposed DEWT-SVD source number estimation method exhibits superior performance.  It can be clearly seen from Figure 10 that the value of Number 4 is close to 0, which is also an obvious inflection point. At the same time, due to the low signal-to-noise of the bearing acoustic signal, there is a large noise interference in the low-frequency region, so set N = 4. The angle domain waveform and envelope order spectrogram of each component state after decomposition are presented in Figures 11 and 12, respectively. It can be clearly seen from Figure 10 that the value of Number 4 is close to 0, which is also an obvious inflection point. At the same time, due to the low signal-to-noise of the bearing acoustic signal, there is a large noise interference in the low-frequency region, so set N = 4. The angle domain waveform and envelope order spectrogram of each component state after decomposition are presented in Figures 11 and 12, respectively. By analyzing the time domain of each component of DEWT, the obvious periodic impact can be observed. The corresponding order of fault can also be clearly determined by the envelope order spectrum of each component. The second order of the rotation frequency can be clearly found in component 1, and the obvious fault characteristic order corresponding to the outer ring, rolling body and inner ring can be respectively found in components 2, 3 and 4. The experimental results demonstrate that the proposed method can effectively extract the fault characteristic order in a complex environment. To prove the superiority of the proposed method as compared to traditional methods, the traditional EWT and EMD methods were used for comparison.

Comparison with Conventional EWT
First, by using the same treatment method, the adjacent singular value difference graph of EWT-SVD is presented in Figure 13. Because there is no clear zero value point, the number of original signals cannot be clearly estimated from the figure.
methods, the traditional EWT and EMD methods were used for comparison.

Comparison with Conventional EWT
First, by using the same treatment method, the adjacent singular value difference graph of EWT-SVD is presented in Figure 13. Because there is no clear zero value point, the number of original signals cannot be clearly estimated from the figure.
For better comparison with the method proposed in this paper, the number of separated modes was also set to four. EWT was used to decompose the original signal in the same way. As a result, no obvious cyclical impact is apparent in the angle domain spectrum (Figure 14), and no fault features can be found in the corresponding envelope order spectra (Figure 15) of the modes extracted by EWT. Figure 13. A line graph of adjacent singular value differences obtained by the EWT-SVD method.
These results demonstrate the problems that EWT presents in dealing with the signals mentioned herein. Due to the existence of high-frequency modulation phenomena, the EWT method cannot detect the correct boundary, and thus cannot extract the fault features. Therefore, it is evident from Figure 15 that the components still contain redundant information after EWT decomposition.  For better comparison with the method proposed in this paper, the number of separated modes was also set to four. EWT was used to decompose the original signal in the same way. As a result, no obvious cyclical impact is apparent in the angle domain spectrum (Figure 14), and no fault features can be found in the corresponding envelope order spectra (Figure 15) of the modes extracted by EWT.     These results demonstrate the problems that EWT presents in dealing with the signals mentioned herein. Due to the existence of high-frequency modulation phenomena, the EWT method cannot detect the correct boundary, and thus cannot extract the fault features. Therefore, it is evident from Figure 15 that the components still contain redundant information after EWT decomposition.

Comparison with EMD
The adjacent singular value difference graph of EMD-SVD is presented in Figure 16, from which it is clear that the EMD-SVD method incorrectly estimated the number of sources as two. 3. It seems that the font size of the caption of Figure 16 is different from those of the other figures, please make a double-check. EMD was used to decompose the original signal in the same way. As a result, no fault features can be found in the corresponding envelope order spectra of the modes extracted by EMD, as shown in Figure 17.  The results in this section clearly reflect the difficulties of EMD, such as the mode aliasing problem (which is evident from Figure 17; each component contains almost the same information), which is also a popular research topic by scholars.
The results presented in this section demonstrate that neither EWT nor EMD can effectively extract the fault characteristic order in such a complicated environment. The problems revealed in the experimental results are in line with the authors' expectations, and prove the correctness of the experiment. The superiority of the method proposed in this paper is therefore well verified.

Conclusions
An improved EWT strategy for the compound weak fault diagnosis of rolling bearings is proposed, and the EWT is improved from two aspects. First, DEWT is proposed to solve the high-frequency modulation phenomenon, which effectively improves its boundary detection capability. Then a source number estimation method based on the DEWT and SVD is proposed to adaptively determine the number of the modes. Through the DEWT and mode number estimation method, the method proposed in this paper can be better applied to the composite weak fault diagnosis of rolling bearings.
A fault simulation experiment of rolling bearing with variable speed is performed, and the acoustic signal of a weak compound fault is collected. Through preliminary analysis of experimental data, the fault characteristics are very weak and contain huge noise interference. Experimental results prove that the proposed method has excellent effects, even on variable speed acoustic signals with extremely low signal-to-noise ratio.