Fault feature extraction method for rolling bearing based on MVMD and complex Fourier transform

. The vibration signals caused by rolling bearing defects in different directions may be different, and the fault diagnosis based on single channel vibration signals may be made incorrectly, and the observation results may be understood wrong. To avoid it, a new rolling bearing fault feature extraction method based on multivariate variational mode decomposition (MVMD) and complex Fourier transform (CFT) were proposed. First, the orthogonally sampled vibration signals were combined into a multivariate signal, and the multivariate signal was decomposed into several intrinsic mode functions (IMFs) using the MVMD. As per this method, a unified mathematical model was used to model vibration signals in two directions, ensuring that fault features were decomposed to the same level. Finally, the CFT was applied to fuse the envelope signals in two directions in order to obtain a clearer and comprehensive amplitude-frequency feature. Simulation and test results verify the feasibility and superiority of the proposed method.


Introduction
Rolling bearing is an important component of rotating machinery. However, because of complex working environments, wear, fatigue, corrosion, overloading and many other factors, local defects can occur in rolling bearing. If left unchecked, these defects may damage the whole system. Therefore, it is very important to monitor and diagnose the condition of the rolling bearing [1][2][3].
Since the vibration signals collected by sensors installed in rolling bearing housings contain a wealth of information about the machine functionality, the vibration analysis is the preferred method for diagnosing local defects in rolling bearings. When a local defect occurs in a rolling bearing, the periodic impact caused by the defect will excite the resonance frequency of the bearing and adjacent components [4]. However, there are various types of rolling bearing defects, such as rolling element defects, outer race defects, inner race defects and their composite defects. Therefore, the frequency distribution of the vibration signals caused by rolling defects is more complicated. Some are in high-frequency bands, some are in mid-frequency bands, and some are even in low-frequency bands, or in the combinations of two or even three bands. However, there is a problem of correct determination of the frequency band to locate the fault signal. In addition, the energy of the vibration signals caused by a weak defect is low, and the signals collected by sensors inevitably contain noise components, which makes the analysis of the vibration signals more difficult [5].
The usage of the wavelet transform (WT) to analyze vibration signals is one of the main ways to monitor rolling bearings in an early stage [6]. Because the detection effect based on WT is not only restricted by the Heisenberg uncertainty principle but also depends on the selection of the basis function and the decomposition scale, this technique cannot guarantee an optimal detection result. With the development of digital signal analysis and processing technology, some data-driven signal analysis methods with good adaptability have been proposed, such as empirical Different from displacement signals, the track formed by acceleration signals does not have rotation characteristics. The analysis of acceleration signals with FVS lost its physical significance. The integration method of the signal characteristics in the orthogonal direction of rolling bearings needs further discussion.
Based on the definition of multivariate modulation oscillation signals, Rehan proposed a multivariate VMD (MVMD) method to extract the set of bandwidth modes that contain inherent multivariate modulation oscillations in multivariate input signals, and the sum of the bandwidths of all signal modes residing in input data channels is the smallest [30]. In the MVMD, the same mathematical model is used to model multi-channel signals to ensure that fault features are decomposed to the same layer. The MVMD was originally used to analyse the EEG signals. Unlike the Fourier transform, the CFT no longer satisfies the conjugate symmetry. From the characteristic frequency and conjugate characteristic frequency, the one with the larger amplitude is selected to characterize fault features, so that the obtained fault features are clearer. A new method for the fault feature extraction of rolling bearings based on MVMD and CFT is proposed.
The main contributions of this paper can be summarized as follows: 1) There is no mode aliasing and mode splitting phenomena when applying MVMD to decompose the muti-channel vibration signals.
2) The fault features in the multi-channel vibration signals are decomposed by MVMD into the multiple modulated oscillation signals of the same layer to facilitate the subsequent feature information fusion.
3) The fault features of the vibration signals of the two channels can be effectively fused from the feature layer by CFT. 4) Selection of the fault characteristic frequency or the conjugate fault characteristic frequency from the complex Fourier spectrum according to the phase difference can highlight the fault feature.
The remainder of this paper is organized as follows: Section 2 introduces the theoretical background of the method, including MVMD and CFT. Section 3 describes the wholediagnostic process of the proposed method for roll bearing. Section 4 consists of two parts: experimental setup and experimental data analysis. The method feasibility is verified by two real cases and compared with other methods. Section 5 contains discussion, which mainly involves the influence on the decomposition results when the parameters of MVMD take different values. Section 6 contains the conclusions and prospects for future works.

Multivariate modulated oscillations
To express conveniently a multi-oscillation signal, the amplitude modulation and frequency modulation signal is defined in the form of vector { ( )} in the real number field: where ( ) and ∅ ( ) are the amplitude function and the phase function, respectively, and is the number of channels. Hilbert operator is applied to Eq. (1), and the corresponding analytic vector of ( ) is obtained. The analytic vector ( ) is as follows [30]: The amplitude function ( ) and the phase function ∅ ( ) can be obtained by the following modulo and phase operations: The original real-value signal vector ( ) can be obtained from the real part of ( ): For vibration signals of rolling bearing defects collected by 2D or 3D sensors, so one or more components with the same frequency may be present. Therefore, a simple multivariate analytic signal model for ( ) that assumes a single common component among all data channels is proposed:

Multivariate VMD
MVMD is an extension of the VMD method. For more details regarding the VMD, see [11]. As a generalization of the VMD algorithm in multidimensional space, the main goal of MVMD is to extract the predefined multivariate modulation oscillation signals ( ) containing C channel data from the input data ( ). MVMD is used to decompose a multivariate signal ( ) into the sum of multiple modulated oscillation signals ( ): where is the number of multivariable modulated oscillation signals.
The purpose of MVMD is to extract a set of multivariable modulation oscillations { ( )} from the input data. Additionally, the sum of bandwidths of the extracted sets is the smallest, and the original signals can be reconstructed accurately. To this end, Eq. (2) is used to represent an analytical vector, and the bandwidth of ( ) can be estimated by the L2 norm of the gradient function of ( ). Therefore, the multivariable function that needs to be optimized in MVMD is expressed by the following formula [30]: The symbol represents the partial derivative operation on time.
Another characteristic of Eq. (8) is that a single frequency component is used for harmonic mixing of the entire vector ( ). Therefore, the design provides the possibility to find the multivariate oscillation with a single common frequency component in all channels in ( ), and this causes a multivariate modulation oscillation model defined in Eq. (6). The bandwidth of the modulation multivariable oscillation is estimated by shifting the single-side spectrum of all channels of ( ) to and taking the Frobenius norm of the matrix. It is convenient to express as follows: where , ( ) represents an analytical modulation signal corresponding to channel number and mode number . In contrast to the vector signal in Eq. (8), , ( ) is a single component complex value signal. The constrained optimization problem of MVMD is as follows: where ∑ , ( ) = ( ), = 1,2, ⋯ , , { , ( )} and { } represent the set of all modes and their central frequencies, respectively.
The corresponding augmented Lagrangian function becomes: where is the penalty factor. In the MVMD, the alternating direction multiplier (ADMM) is used to solve the above unconstrained optimization problem. The ADMM method transforms the complex optimization problem in Eq. (11) into several simpler sub-optimization problems through Eq. (12)(13)(14): where is the time step. Eq. (12) can be equivalent to the following formula: One of these sub-problems is depicted in the updated diagram of the original VMD. In the MVMD, the relationship is updated using the following pattern [12] in the Fourier domain: For the renewal of the central frequency of Eq. (13), considering that the last two terms of the Lagrangian function Eq. (11) do not depend on , the related problems are simplified as follows: By using Plancherel's theorem of the inner product of the function in the time domain and frequency domain, the optimization process can be carried out more easily in the frequency domain. In the Fourier domain, the equivalence problem of Eq. (17) becomes: By setting the first derivative of the quadratic function to zero and then performing some simple algebraic operations, the above sum of the quadratic function is minimized, and the following relations are obtained: It can be seen from Eq. (19) that when updating of each mode during MVMD, the power spectrum contribution of all channels is considered.

Complex Fourier transform
Fourier transform is performed on the basis of a complex sequence signal ( ) = ( ) + j ( ) ( is an integer and ∈ [1, ], is the sequence length, j = −1) to obtain ( ): where FFT stands for fast Fourier transform.
Using the Fourier transform in the finite-length sequence of real and even functions on the /2 symmetry and the Fourier transform of the pure imaginary odd-function finite-length sequence on the /2 anti-symmetric and the linear properties of the Fourier transform, it is possible to get: where ( ) and ( ) are the Fourier transforms of the cosine signals ( ) and ( ) respectively, = 1, 2,…, /2. Because the amplitude of the Fourier transform of the complex sequence ( ) is related to | ( )|, if it is represented by a trigonometric function, it can be found that | ( )| is related to the initial phases of and .
The sampling frequency is equal to 800 Hz. The characteristic frequency amplitude obtained by analyzing ( ) with the phase difference of its imaginary and real parts based on the CFT and FVS methods, respectively, is shown in Fig. 1. In Fig. 1, when ∈ (− , 0), the amplitude of the characteristic frequency based on CFT is larger. When ∈ (0, ), the amplitude of the conjugate characteristic frequency ( − ) based on CFT is larger. Therefore, the characteristic frequency of the complex signal can be selected according to the phase difference to characterize the fault.  21) is applied to to obtain a complex Fourier spectrum, and to select the fault feature frequency according to the phase difference between the imaginary part and the real part of .
The key point of the whole procedure is selection. Generally, the vibration signals caused by rolling bearing defects have four components: fundamental frequencies, defect frequencies, natural vibration frequencies and ultrasonic frequency. Therefore, when the MVMD method is used to analyze vibration signals caused by rolling bearing defects, is set to 4.

Experimental Setup
In this paper, the test bench shown in [31] is used. As described in [31], the experiment rig consists of a rotor-bearings system driven by a 1 HP triple-phase asynchronous motor, as shown in Fig. 2(a). Multiple rotating speeds, which fluctuation range is almost 60 rpm, can be monitored by a motor speed controller. The vibration data were obtained from 3D-accelerometer mounted on a rotor-bearing system. Different kinds of bearing faults are simulated by replacing bearing 1 with a rolling element defective bearing or compound defective bearing having roller defect and outer race defect. The defect on the outer race and ball of bearings is shown in Fig. 2(c) and (d). All channels of vibration data from the 3D-accelerometer were acquired simultaneously using a LMS SCADAS Mobile data acquisition system in Fig. 2 Fig. 3(a) and (b), the original vibration signals and have no obvious amplitude modulation characteristics, and the fault feature information is submerged by noise. In Fig. 3(g) and Fig. 3(h), the and contain very obvious modulation characteristics, but the amplitude of is larger than that of . In Fig. 4, the central frequency of IMFs fluctuated greatly at the early stage of evolution, but it entered a stable state after 24 evolutions.

Fig. 4. Evolution of IMFk center-frequencies
To illustrate the frequency band structure of a vibration signal excited due to the bearing defect, the reason for taking = 4, and the advantage of using the MVMD to analyze multi-channel signals, the Fourier spectrums of the original signals and the Fourier spectrums of their decomposition results based on MVMD are shown in Fig. 5. As generally accepted, the Fourier spectrums of the rolling bearing fault vibration signals contain four subbands. Therefore, it is appropriate to set to 4 when using MVMD to decompose the rolling bearing fault vibration signals. In addition, and have similar spectral structures and are arranged in order of frequency from low to high, and there are no modal aliasing and modal splitting between the subbands of and the subbands of . The fault features obtained by using the MVMD method are decomposed into the same layer, which is beneficial to the subsequent fault feature fusion. Hilbert operations are respectively performed on and , and the envelope functions and are obtained. The Fourier spectrums of and are shown in Fig. 6.  Fig. 6, the rolling element fault feature frequency = 101 Hz is in the mid-frequency subband 3. In addition, the amplitude of is much larger than that of . Let = + j , then the phase difference spectrum will be between and . So it is possible to depict the complex Fourier spectrum of and full vector spectrum of as shown in Fig. 7. Because the phase difference between and is between -and 0, the amplitude of the feature frequency , which is 3.563, in the complex Fourier spectrum of is larger than that of the conjugate frequency ( -), which is 3.456. According to Fig. 6 and Fig. 7, it can be seen that the amplitude of is the largest in the complex Fourier spectrums, i.e. 0.3563, followed by 0.3506 in the full vector spectrum, and it is the smallest in the Fourier spectrum of , i.e. 0.3223. The Fourier spectrums of the first three-order IMFs and obtained by the BEMD method are shown in Fig. 8. The Fourier spectrums of first three-order product functions (PFs) and obtained by the CLMD method are shown in Fig. 9. The envelope spectrums of the IMFs and PFs are shown in Fig. 10 and Fig. 11, respectively.
Obviously, the IMFs obtained by the BEMD method or the PFs obtained by the CLMD method induce modal aliasing and mode splitting, which will cause the extracted fault features to be indistinct. For example, Fig. 8(a) shows that the mid-frequency subband is almost submerged by the dominant high-frequency subband, and mode splitting occurs in the mid-frequency subband. As a result, the amplitude of the rolling element fault feature frequency 2 in Fig. 10(a) is very small and almost submerged. Fig. 8(b) demonstrates mode aliasing in , and the mid-frequency subbands containing fault features are decomposed into and , obviously the mode splitting occurs. As a result, the amplitude of the fault feature frequency in Fig. 10(b) and Fig. 10(d) is reduced. when applying the CLMD method to analyze rolling element defects, the similar situation as described above occurs, as shown in Fig. 11.   Fig. 8. Fourier spectrum of first three-order IMFs and obtained by BEMD method Fig. 9. Fourier spectrum of first three-order PFs and obtained by CLMD method In Fig. 10 and Fig. 11, the fault features on the same layer are inconsistent. This means the information fusion becomes increasingly challenging. Comparing the envelope spectrums of IMFs and PFs, the amplitude of the fault characteristic frequency of the vibration signal in the vertical direction is obviously larger than that in the horizontal direction. It is very likely to make a wrong diagnosis based on only the signal feature in the horizontal direction. The full vector envelope spectrums of and are obtained by applying the FVS to fuse the envelope signals of and ( = 1, 2) as shown in Fig. 12. The complex envelope spectrums of are obtained by applying the CFT to fuse the envelope signals of and as shown in Fig. 13. The maximum amplitudes of the rolling element fault feature frequency obtained by the above methods are shown in Table 1. The amplitude of obtained by the MVMD method is much larger than that obtained by the BEMD or CLMD methods. The amplitude obtained by the proposed MVMD+CFT method is larger than that obtained by the MVMD+FVS method.

Composite defects
The time-domain waveform and Fourier spectrum of the composite fault vibration signals with rolling element defects and outer race defects are shown in Fig. 14 This property of MVMD provides the possibility to determine where the fault signature is located. In Fig. 15, the frequency centers of the subbands in the spectrums of and are approximately the same that is very beneficial to the subsequent fault feature information fusion.  The HT is applied to and respectively, and the envelope spectrums of and are shown in Fig. 16. The rolling element defect feature frequency = 101 Hz is still in the mid-frequency subband 3. The outer race defect feature frequency = 160 Hz is in lowfrequency subband 1, subband 2 and in high-frequency subband 4. Compared with Fig. 6, it can be seen that when only the rolling element defect occurs, the fault feature amplitude in the vertical direction is larger than that in the horizontal direction. But when a composite defect occurs, the situation changes cardinally to the reverse. Let = + , where and are the envelope signals of and respectively. Then the phase difference spectrums of and shown in Fig. 17, and the complex Fourier spectrums and the full vector spectrums of shown in Fig. 18 and in Fig. 19 are actual. In Fig. 18, the larger amplitude is selected from the characteristic frequency and conjugate characteristic frequency to reveal the characteristics of composite defects according to the phase difference between and . Taking and as two examples, In Fig. 17(c) and (d), the phase difference between and at the ball defect characteristic frequency is 6.481° and the phase difference between and at the outer race defect characteristic frequency is 82.31°, which belongs to (0, ). According to 2.3, when the phase difference is between (0, ), the magnitude of the conjugate characteristic frequency is larger than that of the characteristic frequency. In Fig. 18, the magnitude of the ball defect characteristic frequency is 0.2053, while the magnitude of the conjugate characteristic frequency is 0.2228. The magnitude of the outer race fault characteristic frequency is 0.238, while the magnitude of the conjugate characteristic frequency ( -) is 0.5825. In Fig. 19, the amplitude of in the full vector spectrum of is only 0.4102. In Fig. 16, the amplitude of in the envelope spectrum of is larger than that of , i.e. 0.4094, but much smaller than that in the full vector spectrums or in the complex Fourier spectrums. The BEMD is used to decompose the composite defect signals into a series of IMFs and , where = 1, 2, 3, 4. The Fourier specturms of the first four and are shown in Fig. 20. Their envelope specturms are shown in Fig. 21. Obviously, there are modal aliasing and modal splitting in the spectrum of IMFs shown in Fig. 20. The envelope spectrums of and have significant differences as shown in Fig. 21. The outer race fault feature frequency is in the envelope spectrums of and , and the rolling element fault feature frequency is in . In the envelope spectrum of , each of them has an outer race fault feature frequency , but without rolling element fault feature frequency . Obviously, the diagnosis based only on the envelope spectrum of is easy to induce errors and inaccuracies.   The full vector envelope spectrums of the IMFs obtained by BEMD are shown in Fig. 22. The complex envelope spectrums of obtained by CLMD are shown in Fig. 23. In Fig. 22 and Fig. 23, the rolling element fault features are small and difficult to be identified.
The maximum amplitudes of the compound defect feature frequencies and obtained by the above methods are shown in Table 2. It can be seen from Table 2 that even if compound defects are analyzed, the amplitude of defect features extracted by the proposed method is clearer. Using fast spectral kurtosis based on 4th order statistics, it is possible demodulate the rolling element defect feature frequency and its frequency multiplication from the signal , and the outer race defect feature frequency and its frequency multiplication from the signal . The demodulation result is shown in Fig. 24. Obviously, only the rolling element defect is diagnosed according to the envelope spectrum of , and the only outer race defect is diagnosed according to the envelope spectrum of .

Discussion
Similar to BEMD and CLMD, the MVMD separates multi-channel signals into a series of IMFs in the order of frequency. However, the IMFs obtained by the BEMD or CLMD methods are arranged in order of frequency from the highest to the lowest, while the IMFs obtained by the MVMD method are the opposite. In the envelope spectrum of high-order IMFs, the amplitude of the fault feature frequencies is created as per the selected signal decomposition method with the distribution of the original fault signals in the frequency band. When the fault signals exist only in one frequency band, as shown in section 4.2, the rolling element defect signals are located only in the mid-frequency band. The amplitude of the fault feature frequency obtained by the MVMD method is much higher than that of the BEMD method or the CLMD method. The reason is that the mode splitting occurs when the fault signals are decomposed by the BEMD method or the CLMD method, which causes the fault signals to be decomposed into different layers, resulting in a reduction in amplitude. When the fault signals are distributed in different frequency bands, such as composite defect signals mentioned in section 4.3, the MVMD can decompose this type of signals into the corresponding frequency bands, while the BEMD has modal aliasing, and the separated high-order IMF contains fault signals in almost all frequency bands, so the amplitude of fault feature frequency obtained by the BEMD is relatively higher.   Bw= 3200Hz, f c =11200Hz The MVMD has a solid mathematical foundation. Because a unified mathematical model is used for multi-channel signals. In the MVMD, it ensures that the same fault features residing in multi-channel signals are decomposed to the same layer to lay a good foundation for subsequent multi-channel information fusion. The test proves that the vibration signals excited by rolling bearing defects are more complicated. Vibration signals from different directions may inform about different defects. The diagnosis based on the information of vibration signals which come from a single direction can induce wrong conclusions and misunderstanding. The distribution of fault vibration signals is also quite different. For example, the rolling element fault feature frequencies are distributed in the mid-frequency band, and the outer race fault feature frequencies are widely distributed, and almost all frequency bands exist.
There are three benefits of analyzing multi-channel rolling bearing fault signals with the MVMD and CFT. Firstly, the signals in the same central frequencies band are decomposed into the same layer, which is conducive to information fusion. Secondly, compared with the BEMD or CLMD method, there is no mode splitting or mode aliasing in the decomposition results based on the MVMD. Thirdly, the amplitude of fault feature frequency obtained by the CFT is larger than that obtained by the FVS. The fault characteristics obtained by using the proposed method are more obvious, especially when analyzing composite faults.
As a newly emerging algorithm, the MVMD also has some shortcomings. For example, the decomposition levels and end conditions of MVMD need to be determined in advance, thus reducing the adaptability of the MVMD. Additionally, as demonstrated by this paper, is determined based on the Fourier spectrum of original signals. The affects the decomposition effect of the MVMD. is lager, so the tightness of the decomposition result is better, but the fault information is easy to lose. With a smaller , the tightness of the decomposition result is worse, and the modal aliasing is prone to occur, although the fault information will not be lost. Experiments have revealed that when is between 200 and 2000, the decomposition results are roughly the same.

Conclusions
The MVMD is introduced into the field of rolling bearing fault diagnosis. The multi-channel data analysis method based on the MVMD and the multi-channel data fusion method based on the CFT are investigated. The main conclusions are as follows: 1) The same fault features are decomposed into the same layer, which facilitates the information fusion of multi-channel signals.
2) There are no mode aliasing and mode splitting which is easy to occur when using the BEMD or CLMD to analyze vibration signals.
3) By using the CFT to fuse the envelope characteristics of multi-channel vibration signals, the obtained spectrum characteristics are clearer and more comprehensive. 4) When ∈ (− , 0), the amplitude of the characteristic frequency based on CFT is larger. When ∈ (0, ), the amplitude of the conjugate characteristic frequency ( -) based on CFT is larger.
In a future work, the following situations should be well addressed: (1) Research on the distribution law of rolling bearing defect features; (2) Explaining the reasons for the signal difference from different rolling channels in terms of rolling bearing defects; (3) Adaptive determination of the number of MVMD layers and end conditions.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.