# Instantaneous Amplitude-Frequency Feature Extraction for Rotor Fault Based on BEMD and Hilbert Transform

Authors: Chuanjin Huang; Haijun Song; Wenping Lei; Zhanya Niu; Yajun Meng
Journal: Shock and Vibration
 (2019)

---

## Abstract

The vibration signals propagating in different directions from rotating machines can contain a variety of characteristic information. A novel feature extraction method based on bivariate empirical mode decomposition (BEMD) for rotor is proposed to comprehensively extract the fault features. In this work, the number of signal projection directions is determined through simulation, and the energy end condition based on the energy threshold is increased using BEMD to enhance the decomposition quality. Mixed vibration signals are generated along two orthogonal directions. Then, the acquired vibration signal can be decomposed into several intrinsic mode functions (IMFs) at the rotational speed using the BEMD method. Furthermore, the instantaneous frequency and instantaneous amplitude of the real signals and the imaginary part of the IMF signals are obtained using the Hilbert transform. The fault features along two and three dimensions can be investigated, providing more comprehensive information to aid in the fault diagnosis of rotor. Experimental results on oil film oscillation, the oil whirl, the bistability of the rotor, and looseness and rotor rubbing composite fault indicate the effectiveness of the proposed method.


---

## Body

## 1. Introduction

The features of vibration signals form the basis of fault diagnosis of rotating machinery. A variety of common fault feature extraction methods, such as the wavelet transform (WT) [1], empirical mode decomposition (EMD) [2], and local mean decomposition (LMD) [3, 4], are widely used for such diagnosis. In the WT, the wavelet function must be preset, which leads to limitations in the self-adaptability of the WT. The EMD and LMD methods adaptively decompose the signals into a series of intrinsic mode functions (IMFs) and product functions (PFs), respectively, according to the characteristics of the signal envelope. However, these traditional methods are suitable only for processing real-value data on a single channel [5].The fault vibration characteristics from two orthogonal sensors, which are placed in various locations on the machine, have been shown to be important for diagnosing faults. To improve the reliability of the diagnostic results, the feature extraction method based on multichannel signal features is becoming increasingly more popular.The orbit of each harmonic from the vibration signals is an ellipse when the rotor runs in a steady state. The fault feature extraction methods based on homology information, such as the holospectrum [6], full spectrum [7], and full vector spectrum [8], have been proposed in previous reports. The vibration signals in the orthogonal direction can form an ellipse. Because the homologous information incorporates the orthogonal direction of the vibration signal features, the diagnostic results are more comprehensive and accurate [9]. Fault features can be described by the homologous information through round or elliptical information. Unfortunately, in most cases, the circle or ellipse information can be extracted using only a Fourier transform. Therefore, the homologous information method is applicable primarily to the analysis of stationary signals.The EMD or LMD methods can locally smooth the nonlinear signal. Homologous information technology combined with EMD and LMD [10, 11] has been proposed to extract the features of the nonlinear signal. However, the EMD and LMD methods decompose the signals from different directional sensors individually and separately. The circle or ellipse information is calculated according to the one-to-one mode principle [11]. When the EMD or LMD is used to decompose the orthogonal vibration signals, the number of decomposition result is nonuniform, which causes the information fusion to become increasingly challenging [12]. In addition, the EMD or LMD method is relatively sensitive to noise, which causes the IMF group or PF group frequency to be inconsistent with other IMF or PF groups.To improve the capacity of EMD and LMD, these methods have been extended recently to process multivariate signals. Some examples include bivariate empirical mode decomposition (BEMD) [13], complex local mean decomposition (CLMD) [14], trivariate EMD (TEMD) [15], and multivariate EMD (MEMD) [16]. However, these extended methods were not specifically developed for machine fault feature extraction. In BEMD, the local mean of the bivariate signal is calculated by projecting the signal to a number of directions. Following the same idea, the TEMD and MEMD were proposed, respectively, for three-dimensional and n-dimensional signals. The CLMD estimates the local mean by projecting only a bivariate signal onto the x- and y-axes. The bivariate EMD has been employed to detect wind turbine mechanical and electrical faults by decomposing the fault signals [5]. Here, the electrical signal was analyzed without considering the noise and the number of BEMD projections. The associated experiments indicated that the number of projection directions affects the decomposition results, particularly when noise is included. In addition, additional false components were identified in [5]. The MEMD has also been employed to fault diagnosis of rolling bearing [17]. Methods based on BEMD or MEMD are used to analyze and extract signal characteristics in multichannels. However, the joint information between multiple sensors was not considered in [5, 17]. This method, combined with the CLMD or MEMD and full spectrum, has been proposed to obtain joint information among multichannel signals [12, 18]. The methods based on CLMD or BEMD simultaneously decompose the signals into multichannels, thus ensuring that the number of decomposition results is the same and is readily incorporated. However, the instantaneous feature cannot be extracted due to the elliptical information that was generated using the conventional Fourier transform.The signals are analyzed as the superposition of slow and fast oscillations in the EMD, and the bivariate extension and bivariate signals, e.g., the orthogonal vibrations, are analyzed as the superposition of IMFs at the rotational speed in the BEMD. To extract the instantaneous vibration characteristics and the joint information among the multichannel signals, the fault feature extraction method based on the BEMD is proposed to decompose the multicomponent rotation signal into a series of single-component rotation signals. The instantaneous amplitude and instantaneous frequency of the IMFs are further obtained using the Hilbert transform (HT).


## 2. Brief Description and Improvement of BEMD

### 2.1. Brief Description of BEMD

The fundamental concept of BEMD is that bivariate signals are made up of slow rotation signals and fast rotation signals superimposed on the slow rotation signals [13]. For a mixed signal zt, the decomposition process based on the BEMD is as follows [5]:Step 1: determine the number of projectionsN and calculate the projection directions:(1)φn=2nπN,n∈1,N.Step 2: project the complex-valued signalzt on the directions φn:(2)pφnt=Ree−jφnzt,j=−1.Step 3: extract all local maxima ofpφnt:tin,pφntin. Herein, i indicates the number order of individual local maxima.Step 4: interpolate the settin,pφntin by spline interpolation to obtain the tangent along the direction φn:eφnt.Step 5: repeat 2–4 until the tangents in allN projections are obtained.Step 6: compute the mean of all tangents:(3)mt=1N∑n=1Neφnt.Step 7: subtractmt from zt to obtain ht, i.e.,(4)ht=zt−mt.Step 8: perform sifting process whether the stopping criterion similar to the one proposed in [19] is met; however, ht and mt are bivariate signals. If not, regard ht as original signal and repeat Steps 2–7 until the stopping criterion is met.Step 9: record the obtained IMF and remove it fromzt, i.e.,(5)c1t=ht,r1t=zt−c1t.Step 10: taker1t as the original signal and repeat the above calculation until the second IMF c2t is obtained. The remainder is then calculated as follows:(6)r2t=r1t−c2t.Step 11: iterate the previous calculations until acquiring all IMFs contained inzt.After the above process,zt can be expressed as follows:(7)zt=∑k=1Kckt+rKt,where K represents the total number of IMFs.During the decomposition process, the bivariate rotating signal is required to rotate around the zero point. The rotating machinery often revolves around a central movement, which coincides with the requirements of the BEMD. The present study takes the typical Jeffcott rigid rotor-bearing system dynamics equation to structure the complexz (z = x + jy) with x and y representing the vibration signals collected by two orthogonal sensors. The obtained characteristics curve is shown in Figure 1. The BEMD analysis of the zt with the parameter N of 4 is shown in Figure 2, which shows that the target signal is decomposed into bivariate rotation components at the rotation speed, allowing investigation of the instantaneous amplitude-frequency (IAF) characteristics of the main components.Figure 1
The simulated signal ofz: (a) the three-dimensional time-domain waveform; (b) the two-dimensional plots.
(a)
(b)Figure 2
(a)zt and its four directions of the tangent and tangential mean; (b) the decomposition results of zt based on the BEMD; (c) the planes of c1 and c2.
(a)
(b)
(c)


### 2.2. Improved BEMD Method

In the original BEMD algorithm, the loop cutoff condition is such that all IMF components are obtained from the original signal, which leads to additional IMF components and increases the computing time of BEMD. In most cases, the fault characteristics are primarily contained in the IMFs with a higher energy, and the vibration faults corresponding to the rotational modes are limited. In this study, the end condition of the BEMD based on the energy threshold is proposed based on the reasons mentioned above. A ratioλ is set, and the ratio between the signal to be decomposed and the energy of the original signal is less than a specific value that serves as a criterion to stop the BEMD algorithm. λ is calculated using the following formula:(8)∑L=1nabsrL2∑L=1nabszL2≤λ,L=1,2,…,n.


## 2.1. Brief Description of BEMD

The fundamental concept of BEMD is that bivariate signals are made up of slow rotation signals and fast rotation signals superimposed on the slow rotation signals [13]. For a mixed signal zt, the decomposition process based on the BEMD is as follows [5]:Step 1: determine the number of projectionsN and calculate the projection directions:(1)φn=2nπN,n∈1,N.Step 2: project the complex-valued signalzt on the directions φn:(2)pφnt=Ree−jφnzt,j=−1.Step 3: extract all local maxima ofpφnt:tin,pφntin. Herein, i indicates the number order of individual local maxima.Step 4: interpolate the settin,pφntin by spline interpolation to obtain the tangent along the direction φn:eφnt.Step 5: repeat 2–4 until the tangents in allN projections are obtained.Step 6: compute the mean of all tangents:(3)mt=1N∑n=1Neφnt.Step 7: subtractmt from zt to obtain ht, i.e.,(4)ht=zt−mt.Step 8: perform sifting process whether the stopping criterion similar to the one proposed in [19] is met; however, ht and mt are bivariate signals. If not, regard ht as original signal and repeat Steps 2–7 until the stopping criterion is met.Step 9: record the obtained IMF and remove it fromzt, i.e.,(5)c1t=ht,r1t=zt−c1t.Step 10: taker1t as the original signal and repeat the above calculation until the second IMF c2t is obtained. The remainder is then calculated as follows:(6)r2t=r1t−c2t.Step 11: iterate the previous calculations until acquiring all IMFs contained inzt.After the above process,zt can be expressed as follows:(7)zt=∑k=1Kckt+rKt,where K represents the total number of IMFs.During the decomposition process, the bivariate rotating signal is required to rotate around the zero point. The rotating machinery often revolves around a central movement, which coincides with the requirements of the BEMD. The present study takes the typical Jeffcott rigid rotor-bearing system dynamics equation to structure the complexz (z = x + jy) with x and y representing the vibration signals collected by two orthogonal sensors. The obtained characteristics curve is shown in Figure 1. The BEMD analysis of the zt with the parameter N of 4 is shown in Figure 2, which shows that the target signal is decomposed into bivariate rotation components at the rotation speed, allowing investigation of the instantaneous amplitude-frequency (IAF) characteristics of the main components.Figure 1
The simulated signal ofz: (a) the three-dimensional time-domain waveform; (b) the two-dimensional plots.
(a)
(b)Figure 2
(a)zt and its four directions of the tangent and tangential mean; (b) the decomposition results of zt based on the BEMD; (c) the planes of c1 and c2.
(a)
(b)
(c)


## 2.2. Improved BEMD Method

In the original BEMD algorithm, the loop cutoff condition is such that all IMF components are obtained from the original signal, which leads to additional IMF components and increases the computing time of BEMD. In most cases, the fault characteristics are primarily contained in the IMFs with a higher energy, and the vibration faults corresponding to the rotational modes are limited. In this study, the end condition of the BEMD based on the energy threshold is proposed based on the reasons mentioned above. A ratioλ is set, and the ratio between the signal to be decomposed and the energy of the original signal is less than a specific value that serves as a criterion to stop the BEMD algorithm. λ is calculated using the following formula:(8)∑L=1nabsrL2∑L=1nabszL2≤λ,L=1,2,…,n.


## 3. IAF Feature Extraction of the Bivariate Rotation Signal

As noted in Section2, BEMD decomposes the complex-valued signal of the multiple components into the complex-valued signals of a single component. ci is a complex-valued signal with ci=cxi+jcyi. cxi represents the horizontal component of the vibration signal, and cyi represents the vertical component of the vibration signal. The HT is a classical method to obtain the IAF of the signal. The real and the imaginary components of ci are transformed with the HT to obtain the IAF. The following formulae are established:(9)cxit¯=1π∫−∞+∞cxiτt−τdτ,cyit¯=1π∫−∞+∞cyiτt−τdτ.The corresponding resolution signal expressions are(10)xit=cxit+j×cxit¯=axitejΦxit,yit=cyit+j×cyit¯=ayitejΦyit,where axi and ayi are the instantaneous amplitudes of the analytic signals of xi and yi, respectively, Φxi and Φyi are the phase functions of xi and yi, respectively, and(11)Φxit=arctancxit¯cxit,Φxit=arctancxit¯cxit.The instantaneous frequency function is derived from the phase function:(12)fxit=12π×dΦxitdt,fyit=12π×dΦyitdt.We defineaxit and fxit as the instantaneous amplitude and the instantaneous frequency of cxi, respectively, and define ayit and fyit as the instantaneous amplitude and the instantaneous frequency of cyi, respectively.


## 4. Experiment Verification

### 4.1. Analysis of the Oil Film Oscillation Signal

The rotor oil film oscillation test devices are shown in Figure3. The signals are collected using eddy current sensors along the orthogonal direction to the axis. The first-order critical speed is found to be in the range of 3,200 to 3,400 rpm using the resonance phenomena by increasing the speed of the rotor test rig whose speed is increased from 0 to 8,331 rpm and then decreased to 0. Certain oil film vortex and oil film oscillation phenomena occur during the experiment. The data recording device uses the INV306 collector, which has a 16-bit 4-way parallel AD converter. The sampling frequency of the collector was set to 2,048 Hz in the following experiments.Figure 3
The schematic diagram and the experimental apparatus.The oil film oscillation signals alongx and y (where the rotor speed is 6,480 rpm) in the horizontal and vertical directions are collected by eddy current sensors. x and y and their Fourier spectra are shown in Figure 4, where X = 108 Hz is the frequency of the fundamental frequency signal. Let z = x + jy. The three-dimensional time-domain waveform and the two-dimensional plots of z are shown in Figure 5. Figures 4(b) and 4(d) indicate that the amplitude of the 0.48X component is more prominent, particularly in the vertical direction and even exceeds the amplitude of the fundamental frequency signal. Figure 5 indicates that the amplitude of the oscillation signals is variable but that the Fourier spectra cannot distinguish this change.Figure 4
The oil film oscillation signal: (a) the horizontal and vertical direction signalsx; (b) the Fourier spectrum of signal x; (c) the vertical direction signals y; (d) the Fourier spectrum of signal y.
(a)
(b)
(c)
(d)Figure 5
The oil film oscillation signal ofz: (a) the three-dimensional time-domain waveform; (b) the two-dimensional plots.
(a)
(b)The EMD is applied to decompose the signalsx and y, as shown in Figure 6. The S1 and S2 columns represent the decomposition results of x and y, respectively. The number of IMFs generated from signal x is 10, but there are only 8 from signal y. The mismatch in the number of IMFs leads to challenges in information fusion. In addition, modal aliasing occurs in IMF1 from signal x. The mode aliasing obscures the physical meaning of the IMF components and affects the subsequent analysis.Figure 6
The decomposition results of signalsx and y using EMD.The decomposition results of the oil film oscillation signal using BEMD (N=4) are shown in Figure 7. The signal z is decomposed into eight complex rotational components, c1–c8. After decomposition, the IMF numbers from x are the same as those from y. The noise component c1 is more easily separated from signal z relative to the case of IMF1 from signal x. c2 and c3 are considered fundamental and oscillatory components, respectively; the frequency component of c2 is 108 Hz, and the frequency of c3 is 52 Hz. The two trajectories composed of c2 and c3 generated by the BEMD method are shown in Figures 8(a) and 8(b), respectively, while the corresponding ones of IMF2 and IMF3 generated by EMD method are shown in Figures 8(c) and 8(d), respectively. Figure 8 indicates that the orbit arrangements obtained by the BEMD method are superior to those of the EMD method.Figure 7
Decomposition results of the oil film oscillation signal using BEMD (N=4).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 8
The orbits made up of (a) the real and imaginary parts ofc2, (b) real and imaginary parts of c3, (c) IMFx2 and IMFy2, and (d) IMFx3 and IMFy3.
(a)
(b)
(c)
(d)Relative to the EMD method, the decomposition effect is improved when using BEMD to decompose the oil film oscillation signal. However, there are still two problems when using BEMD to decompose the oil film oscillation signal. One is the occurrence of modal aliasing in the IMFs, and the other is increased false components. From Figure6, the modal aliasing occurs at the end of c4, which causes the value of c3 to become small at the end. After decomposition, eight IMF components and one residual r are obtained. However, most of the IMF components contain useless information. In most cases, the fault characteristics are often included in higher-energy components. The projection directions and the energy end condition based on the energy threshold are increased in the improved BEMD method to enhance the decomposition quality.The decomposition results of the oil film oscillation signal are shown in Figure9 using the improved BEMD method (N=16, λ = 0.05). The oil film oscillation signal is decomposed into four complex rotational components c1–c4 and the remaining component r. c1 is considered as the noise component. According to the Fourier analysis, c2 is the fundamental signal with a frequency of 108 Hz, and c3 is the oscillation signal with a frequency of 52 Hz. It is clear that the number of IMFs is reduced, and the fundamental frequency component c2 and the oscillation component c3 are more accurately extracted from the original signal. In particular, the value of the c4 end becomes smaller, and the value of the c3 end becomes larger, relative to the values of c3 and c4 shown in Figure 7. Two trajectories made up of c2 and c3 are shown in Figure 10. From c2 and c3 three-dimensional time domain and the trajectories, the fundamental frequency signal c2 has a small elliptic amplitude change. c3 is a large series of amplitude conversion elliptical compositions that cause a significant oscillation. Relative to the orbit of c3 in Figure 6, the orbit of c3 in Figure 10 is improved, particularly in the center region.Figure 9
Decomposition results of the oil film oscillation signal using the improved BEMD method (N=16, λ = 0.05).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 10
The orbits were made up ofc2 (a) and c3 (b) obtained with the improved BEMD method.
(a)
(b)Increasing the number of signal projection directions results in an increase in the number of projection signals. Then, the tangent mean, which is obtained by interpolating the local maximum of the projected signals with a spline interpolation, is more accurate. However, it is meaningless to continue to increase the number of projection directions when the tangent mean is accurately fitted. It is found that the signal decomposition results ofN=16 are almost the same as N=1024 when considering the complex rotation components separated by the original signal. However, the calculation time of the BEMD algorithm with N=1024 is greatly increased.The HT is applied to the real and imaginary parts ofc2 and c3, respectively, and the instantaneous frequency and instantaneous amplitude are obtained, as shown in Figure 11, where ax2 and ax3 represent the instantaneous amplitude of the real parts of c2 and c3, respectively, ay2 and ay3 represent the instantaneous amplitude of the imaginary parts of c2 and c3, respectively, fx2 and fx3 represent the instantaneous frequency of the real parts of c2 and c3, respectively, and fy2 and fy3 represent the instantaneous frequency of the imaginary parts of c2 and c3, respectively.Figure 11
(a) The instantaneous amplitude of the real part and imaginary part ofc2; (b) the instantaneous frequency of the real part and imaginary part of c2; (c) the instantaneous amplitude of the real part and imaginary part of c3; (d) the instantaneous frequency of the real part and imaginary part of c3.
(a)
(b)
(c)
(d)In Figure11, ax3 is much larger than ay3, and their phases are separated by nearly 180°, which shows that the amplitude and the phase of the oil film oscillation signal in different directions can vary. fx2 and fy2, and fx3 and fy3, are approximately the same. ax2 and ay2 show little change, and their phases are the same. Since BEMD is a bivariate extension of EMD, like EMD, BEMD also has an endpoint effect. There are some fluctuations in the instantaneous amplitude and instantaneous frequency due to the end effect and the edge effect of the Hilbert transform. The three-dimensional time domain of ax2 and ay2 and of ax3 and ay3 is shown in Figure 12. The amplitude range of c3 is much larger than that of c2. It is can be inferred that the main component of oscillation is c3 with the frequency of 52 Hz in the oil film oscillation signal.Figure 12
The 3D time domain of instantaneous amplitude was made up of (a)ax2 and ay2 and (b) ax3 and ay3.
(a)
(b)The decomposition results of the oil film oscillation signal are shown in Figure13 using the CLMD method proposed in reference [14]. The oil film oscillation signal is decomposed into four complex product functions cpf1–cpf4. The noise component is not separated from the oil film oscillation signal using the CLMD method. Moreover, the single component fundamental frequency signal and the oil film oscillation signal were not successfully separated. One of the possible reasons is that in the CLMD algorithm, the complex signal is only projected onto the x-axis and the y-axis, unlike BEMD, which projected on multiple directions. In addition, CLMD is a bivariate extension of LMD. LMD used a moving average algorithm when fitting the signal envelope, which can filter noise to a certain extent. The signals other than the noise component c1 in Figure 9 are added to obtain a filtered oil film oscillation signal, which is then decomposed by the CLMD method, and the first two decomposed results are shown in Figure 13. It is seen that the single component fundamental frequency signal and the single component oil film oscillation signal are separated.Figure 13
The first two decomposed results of the filtered signal based on the CLMD method.
(a)
(b)
(c)
(d)cpf1 and cpf2 consisted of real part signals and imaginary part signals, both of which were composed of the product of the envelope signal and the pure frequency modulation function. The envelope signal is the instantaneous amplitude of the signal. The corresponding instantaneous frequency was obtained by deriving the inverse function of the cosine pure frequency modulation function. The instantaneous amplitude and instantaneous frequency curves are shown in Figure 14. The instantaneous amplitude and frequency obtained by the CLMD method are smoother than in Figure 12. The reason is mainly that the CLMD method uses the moving average filtering algorithm to obtain the signal envelope curve. However, this is the result of using the CLMD algorithm after filtering out noise with BEMD. If the BEMD algorithm is not used for filtering noise, the instantaneous amplitude and instantaneous frequency curves of the single component were not obtained by the CLMD method.Figure 14
(a) The instantaneous amplitude of the real part and imaginary part of cpf1; (b) the instantaneous frequency of the real part and imaginary part of cpf1; (c) the instantaneous amplitude of the real part and imaginary part of cpf2; (d) the instantaneous frequency of the real part and imaginary part of cpf2.
(a)
(b)
(c)
(d)


### 4.2. Analysis of Oil Whirl Signal Based on the Improved BEMD Method

In the method similar to that presented in Section4.1, the oil whirl signals of the rotor test rig with a speed parameter of 4,320 rpm are collected by two orthogonal sensors, as shown in Figure 15. The figure shows the typical whirl phenomenon of large circles with embedded smaller ones. The decomposition results based on the improved BEMD method (N=16, λ = 0.05) are shown in Figure 16.Figure 15
The oil whirl signalz: (a) the 3D time-domain wave of z; (b) the 2D plane of z; (c) the Fourier spectrum of Real[z]; (d) the Fourier spectrum of Imag[z].
(a)
(b)
(c)
(d)Figure 16
The decomposition results of the oil whirl signal based on the improved BEMD method.
(a)
(b)
(c)
(d)
(e)
(f)Only three IMFs appear in Figure16, and the single components c2 and c3 and the noise component c1 are successfully separated from the original signal. The HT is applied to the real and imaginary parts of c2 and c3, respectively, and the instantaneous frequency and instantaneous amplitude are obtained, as shown in Figure 17. The three-dimensional time domain of ax2 and ay2 and of ax3 and ay3 is shown in Figure 18. Figure 17 indicates that the frequency of c2 is approximately twice the frequency of c3 and that ax2 is larger than ay2. In addition, ay3 is slightly larger than ax3, but the range of change for ax3 is greater than the range for ay3. It is inferred that c2 is the fundamental frequency signal and that c3 is the half-frequency signal in the oil whirl signal.Figure 17
The instantaneous amplitude and frequency ofc2 and c3 from the oil whirl signal obtained by the HT: (a) the instantaneous amplitude of the real part and imaginary part of c2; (b) the instantaneous frequency of the real part and imaginary part of c2; (c) the instantaneous amplitude of the real part and imaginary part of c3; (d) the instantaneous frequency of the real part and imaginary part of c3.
(a)
(b)
(c)
(d)Figure 18
The 3D time domain of the instantaneous amplitude of ax2 and ay2 (a) and ax3 and ay3 (b) from the oil whirl signal.
(a)
(b)


### 4.3. Analysis of Looseness and Rotor Rubbing Composite Fault Signal Based on the Improved BEMD Method

Loose and rotor rubbing composite faults are set on the test equipment shown in Figure3 in Section 4.1. Loose fault is set on the nondrive end of the motor, and the distance between the plastic rod and the shaft is fixed near the sensor on the left side of the disk. As the rotor speed increases, the vibration increases and the rubbing fault occurs, which is stable at around 1,700 r/min. The composite fault signals are collected by two orthogonal sensors, as shown in Figure 19. The decomposition results based on the improved BEMD method are shown in Figure 20.Figure 19
The composite fault signalz: (a) the 3D time domain wave of z; (b) the 2D plane of z; (c) the Fourier spectrum of Real[z]; (d) the Fourier spectrum of Imag[z].
(a)
(b)
(c)
(d)Figure 20
The decomposition results of the composite fault signal based on the improved BEMD method.
(a)
(b)
(c)
(d)
(e)Figures19(c) and 19(d) indicate that the signal component mainly contain 1X, 2X (X = 28 Hz), and a frequency modulated signal generated due to time-varying stiffness. This phenomenon is similar to that described reference [20]. Four IMFs appear in Figure 20, and the single components 1X, 2X signals and the FM signal c2 are successfully separated from the original signal. The HT is applied to the real and imaginary parts of c2, c3, and c4, respectively, and the instantaneous frequency and instantaneous amplitude are obtained, as shown in Figure 21. There are some fluctuations in the frequencies of c2 and c3, but these fluctuations are different from the random fluctuations in the above cases. They have obvious regularity and are characteristic of FM signals. The frequency modulation characteristics of c2 are more obvious than those of c3. In addition, by observing the instantaneous amplitude of c2, it is seen that c2 is still an amplitude modulation signal.Figure 21
The instantaneous amplitude and frequency ofc2, c3, and c4 obtained by the HT: (a) the instantaneous amplitude of the real part and imaginary part of c2; (b) the instantaneous frequency of the real part and imaginary part of c2; (c) the instantaneous amplitude of the real part and imaginary part of c3; (d) the instantaneous frequency of the real part and imaginary part of c3; (e) the instantaneous amplitude of the real part and imaginary part of c4; (f) the instantaneous frequency of the real part and imaginary part of c4.
(a)
(b)
(c)
(d)
(e)
(f)In order to further verify the correctness of the instantaneous amplitude-frequency characteristics of the proposed method, the real and imaginary parts of the composite fault signalz are analyzed separately using synchrosqueezed wavelet transforms (SWT) proposed in reference [21]. The results are shown in Figure 22. It is seen that the time-frequency representations of the composite fault signal z also include the AM-FM signal and the 1X signal, which proves the correctness of the proposed method. Compared with the SWT method, the instantaneous amplitude-frequency characteristics acquired by the HT method are relatively straightforward.Figure 22
The results of the composite fault signalz based on SWT: the time-frequency representation of (a) the real part of z and (b) the imaginary part of z.
(a)
(b)


### 4.4. The Bistable Behavior Analysis of the Fan Rotor Based on BEMD

The bistability of the rotor is a nonlinear behavior of the rotor-bearing system, which is the state in which the rotor jumps from one stable state to another, forming a step. The bivariate signal of the bistable behavior is composed of two signals collected by two displacement sensors from orthogonal locations on the experimental devices in literature [22], as shown in Figure 23. Literature [22] shows that the cause of the bistable behavior remains to be further explored. This paper uses this case to illustrate the feasibility of BEMD to analyze nonstationary signals.Figure 23
The schematic diagram of the experimental apparatus.Thex and y signals in the horizontal and vertical directions of the left and right bearings, respectively, from the fan rotors are collected with four displacement sensors. Letting z = x + jy, the time and frequency domain plots of z are shown in Figure 24, where the fan rotor speed is 5500 rpm, the sampling frequency is 2000 Hz, and the number of sampling points is 1024. The left and right columns, respectively, show the time and frequency domain plots of the vibration signals from the left and right bearings of the fan rotor. Bistable behavior arises in the fan rotor, and the amplitudes of the vibration signals vary significantly in different positions and directions. Further studies are required to explain the causes of this bistability. The present study focuses on extracting the bistable behavioral signal characteristics to verify the feasibility of the proposed method.Figure 24
The time and frequency domain plots of the bistable behavior signals.The decomposition results of the bistable behavioral signals based on the improved BEMD method are shown in Figure25. c1, c2, c3, and r are separated in order from z using the improved BEMD method. c1 shows a random arrangement and is considered the high-frequency noise signal. c2 is considered to represent the extracted bistable behavior signals. The HT is applied to the real and imaginary parts of c2 to obtain the instantaneous amplitude and frequency of c2 from the left and right columns from Figure 25, as shown in Figure 26. Figure 27 shows the three-dimensional time domain of ax2 and ay2 from the left and right columns, respectively. Figure 26 shows that the vibration signal amplitude on the left side of the fan decreases from large to small, opposite of the behavior of the right. The horizontal vibration signal amplitude on the left side of the fan is larger than that of the vertical direction signal, opposite of the right. This result validates that the vibration signals from different directions or positions are different when the fan produces bistable behavior. In addition, the time of the bistable behavior can be determined according to the jump point of the amplitude or frequency.Figure 25
The decomposition results of the bistable behavior signals based on the improved BEMD method.Figure 26
The instantaneous amplitude and frequency ofc2 from the (a) left and (b) right columns.
(a)
(b)Figure 27
The three-dimensional time domain ofax2 and ay2.


## 4.1. Analysis of the Oil Film Oscillation Signal

The rotor oil film oscillation test devices are shown in Figure3. The signals are collected using eddy current sensors along the orthogonal direction to the axis. The first-order critical speed is found to be in the range of 3,200 to 3,400 rpm using the resonance phenomena by increasing the speed of the rotor test rig whose speed is increased from 0 to 8,331 rpm and then decreased to 0. Certain oil film vortex and oil film oscillation phenomena occur during the experiment. The data recording device uses the INV306 collector, which has a 16-bit 4-way parallel AD converter. The sampling frequency of the collector was set to 2,048 Hz in the following experiments.Figure 3
The schematic diagram and the experimental apparatus.The oil film oscillation signals alongx and y (where the rotor speed is 6,480 rpm) in the horizontal and vertical directions are collected by eddy current sensors. x and y and their Fourier spectra are shown in Figure 4, where X = 108 Hz is the frequency of the fundamental frequency signal. Let z = x + jy. The three-dimensional time-domain waveform and the two-dimensional plots of z are shown in Figure 5. Figures 4(b) and 4(d) indicate that the amplitude of the 0.48X component is more prominent, particularly in the vertical direction and even exceeds the amplitude of the fundamental frequency signal. Figure 5 indicates that the amplitude of the oscillation signals is variable but that the Fourier spectra cannot distinguish this change.Figure 4
The oil film oscillation signal: (a) the horizontal and vertical direction signalsx; (b) the Fourier spectrum of signal x; (c) the vertical direction signals y; (d) the Fourier spectrum of signal y.
(a)
(b)
(c)
(d)Figure 5
The oil film oscillation signal ofz: (a) the three-dimensional time-domain waveform; (b) the two-dimensional plots.
(a)
(b)The EMD is applied to decompose the signalsx and y, as shown in Figure 6. The S1 and S2 columns represent the decomposition results of x and y, respectively. The number of IMFs generated from signal x is 10, but there are only 8 from signal y. The mismatch in the number of IMFs leads to challenges in information fusion. In addition, modal aliasing occurs in IMF1 from signal x. The mode aliasing obscures the physical meaning of the IMF components and affects the subsequent analysis.Figure 6
The decomposition results of signalsx and y using EMD.The decomposition results of the oil film oscillation signal using BEMD (N=4) are shown in Figure 7. The signal z is decomposed into eight complex rotational components, c1–c8. After decomposition, the IMF numbers from x are the same as those from y. The noise component c1 is more easily separated from signal z relative to the case of IMF1 from signal x. c2 and c3 are considered fundamental and oscillatory components, respectively; the frequency component of c2 is 108 Hz, and the frequency of c3 is 52 Hz. The two trajectories composed of c2 and c3 generated by the BEMD method are shown in Figures 8(a) and 8(b), respectively, while the corresponding ones of IMF2 and IMF3 generated by EMD method are shown in Figures 8(c) and 8(d), respectively. Figure 8 indicates that the orbit arrangements obtained by the BEMD method are superior to those of the EMD method.Figure 7
Decomposition results of the oil film oscillation signal using BEMD (N=4).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 8
The orbits made up of (a) the real and imaginary parts ofc2, (b) real and imaginary parts of c3, (c) IMFx2 and IMFy2, and (d) IMFx3 and IMFy3.
(a)
(b)
(c)
(d)Relative to the EMD method, the decomposition effect is improved when using BEMD to decompose the oil film oscillation signal. However, there are still two problems when using BEMD to decompose the oil film oscillation signal. One is the occurrence of modal aliasing in the IMFs, and the other is increased false components. From Figure6, the modal aliasing occurs at the end of c4, which causes the value of c3 to become small at the end. After decomposition, eight IMF components and one residual r are obtained. However, most of the IMF components contain useless information. In most cases, the fault characteristics are often included in higher-energy components. The projection directions and the energy end condition based on the energy threshold are increased in the improved BEMD method to enhance the decomposition quality.The decomposition results of the oil film oscillation signal are shown in Figure9 using the improved BEMD method (N=16, λ = 0.05). The oil film oscillation signal is decomposed into four complex rotational components c1–c4 and the remaining component r. c1 is considered as the noise component. According to the Fourier analysis, c2 is the fundamental signal with a frequency of 108 Hz, and c3 is the oscillation signal with a frequency of 52 Hz. It is clear that the number of IMFs is reduced, and the fundamental frequency component c2 and the oscillation component c3 are more accurately extracted from the original signal. In particular, the value of the c4 end becomes smaller, and the value of the c3 end becomes larger, relative to the values of c3 and c4 shown in Figure 7. Two trajectories made up of c2 and c3 are shown in Figure 10. From c2 and c3 three-dimensional time domain and the trajectories, the fundamental frequency signal c2 has a small elliptic amplitude change. c3 is a large series of amplitude conversion elliptical compositions that cause a significant oscillation. Relative to the orbit of c3 in Figure 6, the orbit of c3 in Figure 10 is improved, particularly in the center region.Figure 9
Decomposition results of the oil film oscillation signal using the improved BEMD method (N=16, λ = 0.05).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)Figure 10
The orbits were made up ofc2 (a) and c3 (b) obtained with the improved BEMD method.
(a)
(b)Increasing the number of signal projection directions results in an increase in the number of projection signals. Then, the tangent mean, which is obtained by interpolating the local maximum of the projected signals with a spline interpolation, is more accurate. However, it is meaningless to continue to increase the number of projection directions when the tangent mean is accurately fitted. It is found that the signal decomposition results ofN=16 are almost the same as N=1024 when considering the complex rotation components separated by the original signal. However, the calculation time of the BEMD algorithm with N=1024 is greatly increased.The HT is applied to the real and imaginary parts ofc2 and c3, respectively, and the instantaneous frequency and instantaneous amplitude are obtained, as shown in Figure 11, where ax2 and ax3 represent the instantaneous amplitude of the real parts of c2 and c3, respectively, ay2 and ay3 represent the instantaneous amplitude of the imaginary parts of c2 and c3, respectively, fx2 and fx3 represent the instantaneous frequency of the real parts of c2 and c3, respectively, and fy2 and fy3 represent the instantaneous frequency of the imaginary parts of c2 and c3, respectively.Figure 11
(a) The instantaneous amplitude of the real part and imaginary part ofc2; (b) the instantaneous frequency of the real part and imaginary part of c2; (c) the instantaneous amplitude of the real part and imaginary part of c3; (d) the instantaneous frequency of the real part and imaginary part of c3.
(a)
(b)
(c)
(d)In Figure11, ax3 is much larger than ay3, and their phases are separated by nearly 180°, which shows that the amplitude and the phase of the oil film oscillation signal in different directions can vary. fx2 and fy2, and fx3 and fy3, are approximately the same. ax2 and ay2 show little change, and their phases are the same. Since BEMD is a bivariate extension of EMD, like EMD, BEMD also has an endpoint effect. There are some fluctuations in the instantaneous amplitude and instantaneous frequency due to the end effect and the edge effect of the Hilbert transform. The three-dimensional time domain of ax2 and ay2 and of ax3 and ay3 is shown in Figure 12. The amplitude range of c3 is much larger than that of c2. It is can be inferred that the main component of oscillation is c3 with the frequency of 52 Hz in the oil film oscillation signal.Figure 12
The 3D time domain of instantaneous amplitude was made up of (a)ax2 and ay2 and (b) ax3 and ay3.
(a)
(b)The decomposition results of the oil film oscillation signal are shown in Figure13 using the CLMD method proposed in reference [14]. The oil film oscillation signal is decomposed into four complex product functions cpf1–cpf4. The noise component is not separated from the oil film oscillation signal using the CLMD method. Moreover, the single component fundamental frequency signal and the oil film oscillation signal were not successfully separated. One of the possible reasons is that in the CLMD algorithm, the complex signal is only projected onto the x-axis and the y-axis, unlike BEMD, which projected on multiple directions. In addition, CLMD is a bivariate extension of LMD. LMD used a moving average algorithm when fitting the signal envelope, which can filter noise to a certain extent. The signals other than the noise component c1 in Figure 9 are added to obtain a filtered oil film oscillation signal, which is then decomposed by the CLMD method, and the first two decomposed results are shown in Figure 13. It is seen that the single component fundamental frequency signal and the single component oil film oscillation signal are separated.Figure 13
The first two decomposed results of the filtered signal based on the CLMD method.
(a)
(b)
(c)
(d)cpf1 and cpf2 consisted of real part signals and imaginary part signals, both of which were composed of the product of the envelope signal and the pure frequency modulation function. The envelope signal is the instantaneous amplitude of the signal. The corresponding instantaneous frequency was obtained by deriving the inverse function of the cosine pure frequency modulation function. The instantaneous amplitude and instantaneous frequency curves are shown in Figure 14. The instantaneous amplitude and frequency obtained by the CLMD method are smoother than in Figure 12. The reason is mainly that the CLMD method uses the moving average filtering algorithm to obtain the signal envelope curve. However, this is the result of using the CLMD algorithm after filtering out noise with BEMD. If the BEMD algorithm is not used for filtering noise, the instantaneous amplitude and instantaneous frequency curves of the single component were not obtained by the CLMD method.Figure 14
(a) The instantaneous amplitude of the real part and imaginary part of cpf1; (b) the instantaneous frequency of the real part and imaginary part of cpf1; (c) the instantaneous amplitude of the real part and imaginary part of cpf2; (d) the instantaneous frequency of the real part and imaginary part of cpf2.
(a)
(b)
(c)
(d)


## 4.2. Analysis of Oil Whirl Signal Based on the Improved BEMD Method

In the method similar to that presented in Section4.1, the oil whirl signals of the rotor test rig with a speed parameter of 4,320 rpm are collected by two orthogonal sensors, as shown in Figure 15. The figure shows the typical whirl phenomenon of large circles with embedded smaller ones. The decomposition results based on the improved BEMD method (N=16, λ = 0.05) are shown in Figure 16.Figure 15
The oil whirl signalz: (a) the 3D time-domain wave of z; (b) the 2D plane of z; (c) the Fourier spectrum of Real[z]; (d) the Fourier spectrum of Imag[z].
(a)
(b)
(c)
(d)Figure 16
The decomposition results of the oil whirl signal based on the improved BEMD method.
(a)
(b)
(c)
(d)
(e)
(f)Only three IMFs appear in Figure16, and the single components c2 and c3 and the noise component c1 are successfully separated from the original signal. The HT is applied to the real and imaginary parts of c2 and c3, respectively, and the instantaneous frequency and instantaneous amplitude are obtained, as shown in Figure 17. The three-dimensional time domain of ax2 and ay2 and of ax3 and ay3 is shown in Figure 18. Figure 17 indicates that the frequency of c2 is approximately twice the frequency of c3 and that ax2 is larger than ay2. In addition, ay3 is slightly larger than ax3, but the range of change for ax3 is greater than the range for ay3. It is inferred that c2 is the fundamental frequency signal and that c3 is the half-frequency signal in the oil whirl signal.Figure 17
The instantaneous amplitude and frequency ofc2 and c3 from the oil whirl signal obtained by the HT: (a) the instantaneous amplitude of the real part and imaginary part of c2; (b) the instantaneous frequency of the real part and imaginary part of c2; (c) the instantaneous amplitude of the real part and imaginary part of c3; (d) the instantaneous frequency of the real part and imaginary part of c3.
(a)
(b)
(c)
(d)Figure 18
The 3D time domain of the instantaneous amplitude of ax2 and ay2 (a) and ax3 and ay3 (b) from the oil whirl signal.
(a)
(b)


## 4.3. Analysis of Looseness and Rotor Rubbing Composite Fault Signal Based on the Improved BEMD Method

Loose and rotor rubbing composite faults are set on the test equipment shown in Figure3 in Section 4.1. Loose fault is set on the nondrive end of the motor, and the distance between the plastic rod and the shaft is fixed near the sensor on the left side of the disk. As the rotor speed increases, the vibration increases and the rubbing fault occurs, which is stable at around 1,700 r/min. The composite fault signals are collected by two orthogonal sensors, as shown in Figure 19. The decomposition results based on the improved BEMD method are shown in Figure 20.Figure 19
The composite fault signalz: (a) the 3D time domain wave of z; (b) the 2D plane of z; (c) the Fourier spectrum of Real[z]; (d) the Fourier spectrum of Imag[z].
(a)
(b)
(c)
(d)Figure 20
The decomposition results of the composite fault signal based on the improved BEMD method.
(a)
(b)
(c)
(d)
(e)Figures19(c) and 19(d) indicate that the signal component mainly contain 1X, 2X (X = 28 Hz), and a frequency modulated signal generated due to time-varying stiffness. This phenomenon is similar to that described reference [20]. Four IMFs appear in Figure 20, and the single components 1X, 2X signals and the FM signal c2 are successfully separated from the original signal. The HT is applied to the real and imaginary parts of c2, c3, and c4, respectively, and the instantaneous frequency and instantaneous amplitude are obtained, as shown in Figure 21. There are some fluctuations in the frequencies of c2 and c3, but these fluctuations are different from the random fluctuations in the above cases. They have obvious regularity and are characteristic of FM signals. The frequency modulation characteristics of c2 are more obvious than those of c3. In addition, by observing the instantaneous amplitude of c2, it is seen that c2 is still an amplitude modulation signal.Figure 21
The instantaneous amplitude and frequency ofc2, c3, and c4 obtained by the HT: (a) the instantaneous amplitude of the real part and imaginary part of c2; (b) the instantaneous frequency of the real part and imaginary part of c2; (c) the instantaneous amplitude of the real part and imaginary part of c3; (d) the instantaneous frequency of the real part and imaginary part of c3; (e) the instantaneous amplitude of the real part and imaginary part of c4; (f) the instantaneous frequency of the real part and imaginary part of c4.
(a)
(b)
(c)
(d)
(e)
(f)In order to further verify the correctness of the instantaneous amplitude-frequency characteristics of the proposed method, the real and imaginary parts of the composite fault signalz are analyzed separately using synchrosqueezed wavelet transforms (SWT) proposed in reference [21]. The results are shown in Figure 22. It is seen that the time-frequency representations of the composite fault signal z also include the AM-FM signal and the 1X signal, which proves the correctness of the proposed method. Compared with the SWT method, the instantaneous amplitude-frequency characteristics acquired by the HT method are relatively straightforward.Figure 22
The results of the composite fault signalz based on SWT: the time-frequency representation of (a) the real part of z and (b) the imaginary part of z.
(a)
(b)


## 4.4. The Bistable Behavior Analysis of the Fan Rotor Based on BEMD

The bistability of the rotor is a nonlinear behavior of the rotor-bearing system, which is the state in which the rotor jumps from one stable state to another, forming a step. The bivariate signal of the bistable behavior is composed of two signals collected by two displacement sensors from orthogonal locations on the experimental devices in literature [22], as shown in Figure 23. Literature [22] shows that the cause of the bistable behavior remains to be further explored. This paper uses this case to illustrate the feasibility of BEMD to analyze nonstationary signals.Figure 23
The schematic diagram of the experimental apparatus.Thex and y signals in the horizontal and vertical directions of the left and right bearings, respectively, from the fan rotors are collected with four displacement sensors. Letting z = x + jy, the time and frequency domain plots of z are shown in Figure 24, where the fan rotor speed is 5500 rpm, the sampling frequency is 2000 Hz, and the number of sampling points is 1024. The left and right columns, respectively, show the time and frequency domain plots of the vibration signals from the left and right bearings of the fan rotor. Bistable behavior arises in the fan rotor, and the amplitudes of the vibration signals vary significantly in different positions and directions. Further studies are required to explain the causes of this bistability. The present study focuses on extracting the bistable behavioral signal characteristics to verify the feasibility of the proposed method.Figure 24
The time and frequency domain plots of the bistable behavior signals.The decomposition results of the bistable behavioral signals based on the improved BEMD method are shown in Figure25. c1, c2, c3, and r are separated in order from z using the improved BEMD method. c1 shows a random arrangement and is considered the high-frequency noise signal. c2 is considered to represent the extracted bistable behavior signals. The HT is applied to the real and imaginary parts of c2 to obtain the instantaneous amplitude and frequency of c2 from the left and right columns from Figure 25, as shown in Figure 26. Figure 27 shows the three-dimensional time domain of ax2 and ay2 from the left and right columns, respectively. Figure 26 shows that the vibration signal amplitude on the left side of the fan decreases from large to small, opposite of the behavior of the right. The horizontal vibration signal amplitude on the left side of the fan is larger than that of the vertical direction signal, opposite of the right. This result validates that the vibration signals from different directions or positions are different when the fan produces bistable behavior. In addition, the time of the bistable behavior can be determined according to the jump point of the amplitude or frequency.Figure 25
The decomposition results of the bistable behavior signals based on the improved BEMD method.Figure 26
The instantaneous amplitude and frequency ofc2 from the (a) left and (b) right columns.
(a)
(b)Figure 27
The three-dimensional time domain ofax2 and ay2.


## 5. Discussion

The BEMD algorithm decomposes two orthogonal directions of vibration signals as a complex signal, which is a two-dimensional digital signal processing method, thus ensuring that the real and imaginary parts have the same decomposition scale. Similar to EMD, the envelope mean is critical for the decomposition effect of BEMD, but the envelope mean in BEMD is three-dimensional. If the number of projection directions of the complex signal in three-dimensional space is larger, the corresponding envelope signal is also more. Thus, the envelope mean value is more accurate and the BEMD decomposition effect is better. Increasing the number of projection directions can improve modal aliasing. Like EMD, BEMD also produces false components when decomposing signals. Generally speaking, the energy of the false components is low, and these low-energy false components do not contain fault characteristic information, and the introduction of the energy threshold criterion in the termination condition can increase the decomposition speed of the BEMD.The experimental results show that there is a certain difference in the existence of vibration signal characteristics in different directions when rotating machinery fails. In addition, when the number of projection directions is increased, the decomposition speed of BEMD will decrease.


## 6. Conclusions

We use BEMD and HT to extract the instantaneous amplitude-frequency features of rotor faults. A bivariate instantaneous feature extraction method based on the improved BEMD method and the HT is investigated, which extends the fault feature extraction technology to two dimensions. The BEMD method is suitable to analyze the complex, multicomponent bivariate signals. The main single-component bivariate signals are separated from the multicomponent bivariate signals of the fan rotor bistability for the oil film oscillation and the oil film vortex using the improved BEMD method. For the single-component bivariate signal, the HT is used to obtain the corresponding instantaneous amplitude and frequency characteristics. The proposed method can examine the detailed information of a single rotation component.



---

*Source: 1639139-2019-03-12.xml*