The Rational Spline Interpolation Based-LOD Method and Its Application to Rotating Machinery Fault Diagnosis

Local oscillatory-characteristic decomposition (LOD) is a relatively new self-adaptive time-frequency analysis methodology. The method, based on local oscillatory characteristics of the signal itself uses three mathematical operations such as differential, coordinate domain transform, and piecewise linear transform to decompose the multi-component signal into a series of mono-oscillation components (MOCs), which is very suitable for processing multi-component signals. However, in the LOD method, the computational efficiency and real-time processing performance of the algorithm can be significantly improved by the use of piecewise linear transformation, but the MOC component lacks smoothness, resulting in distortion. In order to overcome the disadvantages mentioned above, the rational spline function that spline shape can be adjusted and controlled is introduced into the LOD method instead of the piecewise linear transformation, and the rational spline-local oscillatory-characteristic decomposition (RS-LOD) method is proposed in this paper. Based on the detailed illustration of the principle of RS-LOD method, the RS-LOD, LOD, and empirical mode decomposition (EMD) are compared and analyzed by simulation signals. The results show that the RS-LOD method can significantly improve the problem of poor smoothness of the MOC component in the original LOD method. Moreover, the RS-LOD method is applied to the fault feature extraction of rotating machinery for the multi-component modulation characteristics of rotating machinery fault vibration signals. The analysis results of the rolling bearing and fan gearbox fault vibration signals show that the RS-LOD method can effectively extract the fault feature of the rotating mechanical vibration signals.


Introduction
Rotating machinery, such as gears and rolling bearings, is an important component of mechanical equipment whose fault will have a serious impact on the safe operation of the mechanical device. Therefore, it is of great significance to realize real-time online monitoring of the running condition of rotating machinery [1]. When rotating machinery breaks down, its vibration signals are usually nonlinear and non-stationary. Consequently, traditional time-domain or frequency-domain analysis methods are not suitable for the direct analysis of these signals [2]. The time-frequency analysis method is often used to extract the fault features of the vibration signals of rotating machinery because it can provide both the time-domain and frequency-domain information of signals, and has been adopted which increases the amplitude ratio of the high-and low-frequency components of the signal and amplifies the local oscillation characteristics of the signal, thereby effectively eliminating the occurrence of modal mixing; (2) coordinate domain transformation, which transforms the signal from the original coordinate domain to the sawtooth domain. By increasing the number of sampling points close to the local extreme points of the signal without changing the total number of sampling points, it improves the decomposition accuracy of LOD based on the local extreme points; and (3) piecewise linear transformation, which can process the data between two adjacent extreme points of the signal in parallel without simultaneously knowing the complete local extreme point information of the signal, thereby improving the calculation efficiency and reducing the endpoint effect. A study [21] described the decomposition process of LOD in detail, where LOD was used to decompose a multi-component AM-FM simulation signal x(t): x 1 (t) = (1 + 0.6cos(8πt)) × cos(180πt + 3cos(8πt)), (1) x 2 (t) = sin(25πt), x 3 (t) = sin(2πt), where t ∈ [0, 5], setting sampling frequency 1000 Hz and the simulation signal consists of a modulated signal and two sine signals. The time domain waveforms of x(t) and its three components x 1 (t), x 2 (t), and x 3 (t) are shown in Figure 1. This simulation signal x(t) is decomposed by LOD method and the results are presented in Figure 2. As shown in Figure 2, the original signal x(t) is decomposed into three MOC components and a monotonic residue, and the three MOC components are corresponding to the x 1 (t), x 2 (t), and x 3 (t), respectively. Therefore, LOD, as a self-adaptive time-frequency analysis method based on the characteristics of the signal itself, can reflect the essential characteristics of the signal. However, as illustrated in Figure 3, the locally enlarged view of MOC 3 (t) shows that the waveform of MOC 3 (t) is evidently jagged and exhibits poor smoothness. The main reason is that the use of piecewise linear transformation leads to an unsmooth local mean curve, which causes the gradual accumulation of errors with an increasing number of iterations, thereby affecting the decomposition accuracy. The above decomposition error can be eliminated or reduced if a suitable interpolation method is identified to replace piecewise linear transformation and to obtain a local mean curve with enhanced smoothness. Previous studies [25,26] proposed the adoption of the rational spline interpolation instead of cubic spline interpolation and moving average to calculate local mean curves in EMD and LMD, respectively, which achieved good results. Accordingly, this study intended to introduce the rational spline interpolation to replace the piecewise linear transformation, so as to improve the poor smoothness of the LOD decomposition.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 19 the occurrence of modal mixing; (2) coordinate domain transformation, which transforms the signal from the original coordinate domain to the sawtooth domain. By increasing the number of sampling points close to the local extreme points of the signal without changing the total number of sampling points, it improves the decomposition accuracy of LOD based on the local extreme points; and (3) piecewise linear transformation, which can process the data between two adjacent extreme points of the signal in parallel without simultaneously knowing the complete local extreme point information of the signal, thereby improving the calculation efficiency and reducing the endpoint effect. A study [21] described the decomposition process of LOD in detail, where LOD was used to decompose a multi-component AM-FM simulation signal x(t): x1(t) = (1 + 0.6cos(8πt)) × cos(180πt + 3cos(8πt)), where t ∈ [0, 5], setting sampling frequency 1000 Hz and the simulation signal consists of a modulated signal and two sine signals. The time domain waveforms of x(t) and its three components x1(t), x2(t), and x3(t) are shown in Figure 1. This simulation signal x(t) is decomposed by LOD method and the results are presented in Figure 2. As shown in Figure 2, the original signal x(t) is decomposed into three MOC components and a monotonic residue, and the three MOC components are corresponding to the x1(t), x2(t), and x3(t), respectively. Therefore, LOD, as a self-adaptive time-frequency analysis method based on the characteristics of the signal itself, can reflect the essential characteristics of the signal. However, as illustrated in Figure 3, the locally enlarged view of MOC3(t) shows that the waveform of MOC3(t) is evidently jagged and exhibits poor smoothness. The main reason is that the use of piecewise linear transformation leads to an unsmooth local mean curve, which causes the gradual accumulation of errors with an increasing number of iterations, thereby affecting the decomposition accuracy. The above decomposition error can be eliminated or reduced if a suitable interpolation method is identified to replace piecewise linear transformation and to obtain a local mean curve with enhanced smoothness. Previous studies [25,26] proposed the adoption of the rational spline interpolation instead of cubic spline interpolation and moving average to calculate local mean curves in EMD and LMD, respectively, which achieved good results. Accordingly, this study intended to introduce the rational spline interpolation to replace the piecewise linear transformation, so as to improve the poor smoothness of the LOD decomposition.

Rational Spline Interpolation
Since the piecewise linear transformation lacks enough smoothness, a new alternative interpolation method named rational spline interpolation is introduced in this paper. When using the rational spline interpolation, the shape of the fitting curve will change as the pole parameter p varies. Hence, it is easy for us to determine which curve is more suitable to fit the discrete points. Here the definition of rational spline interpolation function is outlined in some detail [24]. Then in the fragment interval [xi, xi+1] (i = 1, 2, …, n − 1), the rational spline can be described as follows: In formula (5), p is the pole parameter used to control the degree of tightness of the spline curve. Ai, Bi, Ci and Di are the coefficients that make s(x) have continuous first and second derivatives, which can be described as the second derivative value at the node:

Rational Spline Interpolation
Since the piecewise linear transformation lacks enough smoothness, a new alternative interpolation method named rational spline interpolation is introduced in this paper. When using the rational spline interpolation, the shape of the fitting curve will change as the pole parameter p varies. Hence, it is easy for us to determine which curve is more suitable to fit the discrete points. Here the definition of rational spline interpolation function is outlined in some detail [24]. Then in the fragment interval [xi, xi+1] (i = 1, 2, …, n − 1), the rational spline can be described as follows: In formula (5), p is the pole parameter used to control the degree of tightness of the spline curve. Ai, Bi, Ci and Di are the coefficients that make s(x) have continuous first and second derivatives, which can be described as the second derivative value at the node:

Rational Spline Interpolation
Since the piecewise linear transformation lacks enough smoothness, a new alternative interpolation method named rational spline interpolation is introduced in this paper. When using the rational spline interpolation, the shape of the fitting curve will change as the pole parameter p varies. Hence, it is easy for us to determine which curve is more suitable to fit the discrete points. Here the definition of rational spline interpolation function is outlined in some detail [24]. Then in the fragment interval [x i , x i+1 ] (i = 1, 2, . . . , n − 1), the rational spline can be described as follows: In Formula (5), p is the pole parameter used to control the degree of tightness of the spline curve. A i , B i , C i and D i are the coefficients that make s(x) have continuous first and second derivatives, which can be described as the second derivative value at the node: Definition 2. The second derivative at the node can be obtained by the following linear equations: . Equation set (10) contains n − 2 equations, plus boundary conditions s"(x 1 ) = m 1 and s"(x n ) = m n , can form a system of n linear equations with n unknowns. When p ≥ −1, the system of equations is a strictly diagonally dominant matric, and y i " (i = 1, 2, . . . , n) can be solved by chasing method. Then the results of A i , B i , C i , D i , and s i (x) in the interval [x i , x i+1 ] can be obtained by taking y i " into Formulas (6)- (9).
From Formula (5), it can be seen that the pole parameter p can control the shape of the spline curve s(x). When p = −1, s(x) becomes a quadratic spline function showed relatively larger volatility and generally did not apply to strive for fitted curve of discrete points; when p = 0, q = 6, it can be seen from Equation (10) that s(x) is a cubic spline function, so the cubic spline function is a special case of the rational spline function. In addition, as p increases, the last two terms of Formula (5) will gradually decrease, and s(x) will continue to tighten. Furthermore, when p tends to infinity, the last two terms of Formula (5) will become zero, and s(x) will become a piecewise linear interpolation function. Therefore, the piecewise linear transform is also a special case of rational spline interpolation. From the above, we can know that, by adjusting the pole parameter p, the desired fitting curve shape can be obtained theoretically, and the decomposition error can be reduced.
The shape of the local mean curve can be changed by adjusting the pole parameter p, as shown in Figure 4. It can be seen that with the change of p, the degree of tightness and smoothness of the local mean curve has also changed accordingly. Thus, the result and precision of subsequent decomposition will be affected.
Definition 2. The second derivative at the node can be obtained by the following linear equations: where, di = (yi+1 − yi)/hi , q = 2(p 2 + 3p + 3). Equation set (10) contains n − 2 equations, plus boundary conditions s"(x1) = m1 and s"(xn) = mn, can form a system of n linear equations with n unknowns. When p ≥ −1, the system of equations is a strictly diagonally dominant matric, and yi" (i = 1, 2, …, n) can be solved by chasing method. Then the results of Ai, Bi, Ci, Di, and si(x) in the interval [xi, xi+1] can be obtained by taking yi" into formulas (6)-(9).
From formula (5), it can be seen that the pole parameter p can control the shape of the spline curve s(x). When p = −1, s(x) becomes a quadratic spline function showed relatively larger volatility and generally did not apply to strive for fitted curve of discrete points; when p = 0, q = 6, it can be seen from equation (10) that s(x) is a cubic spline function, so the cubic spline function is a special case of the rational spline function. In addition, as p increases, the last two terms of formula (5) will gradually decrease, and s(x) will continue to tighten. Furthermore, when p tends to infinity, the last two terms of formula (5) will become zero, and s(x) will become a piecewise linear interpolation function. Therefore, the piecewise linear transform is also a special case of rational spline interpolation. From the above, we can know that, by adjusting the pole parameter p, the desired fitting curve shape can be obtained theoretically, and the decomposition error can be reduced.
The shape of the local mean curve can be changed by adjusting the pole parameter p, as shown in Figure 4. It can be seen that with the change of p, the degree of tightness and smoothness of the local mean curve has also changed accordingly. Thus, the result and precision of subsequent decomposition will be affected.

Rational Spline Based-LOD Method
Based on the above discussion about the interpolation method, it is not difficult for us to conclude that rational spline interpolation is available for precisely fitting discrete points. Furthermore, the mean function algorithm plays a crucial role in the decomposing results with respect to LOD. Hence, the rational spline-based LOD method (RS-LOD) is proposed in this paper. For a given a signal x(t), the algorithmic procedure of the RS-LOD can be described as follows:

Rational Spline Based-LOD Method
Based on the above discussion about the interpolation method, it is not difficult for us to conclude that rational spline interpolation is available for precisely fitting discrete points. Furthermore, the mean function algorithm plays a crucial role in the decomposing results with respect to LOD. Hence, the rational spline-based LOD method (RS-LOD) is proposed in this paper. For a given a signal x(t), the algorithmic procedure of the RS-LOD can be described as follows: (i) Finding out all local extrema X k , k = 1, 2, . . . , M and corresponding time values τ k , k = 1, 2, . . . , M of original signal x(t). In order to reduce the decomposed error, the function with respect to coordinate domain transform is defined by performing a transformation from original data From the Formulas (11) and (12), we can see that the coordinate transformation does not change the amplitude of the signal (the direction of ordinate), but only compresses or expands the time axis (the direction of abscissa). Hence, the data point density of different signal positions can be changed without changing the amount of data, and thus to highlight the local characteristic information.
(ii) Making a differential operation for the original signal x(t) to get x (t) and finding out the time values τ k , k = 1, 2, . . . , N corresponding to the extrema of x (t). At the same time, finding out function values X k , k = 1, 2, . . . , N corresponding to the time values τ k in original signal x(t).
Then the local mean function m 1 (u) can be derived by using rational spline function to fit the (τ k , X k ) in sawtooth domain. (iii) Subtract the local mean function m 1 (u) from the sawtooth domain function s(u) of the original signal, and then the high-frequency oscillatory function c 1 (u) can be obtained: (iv) Restore the coordinates of c 1 (u) from the sawtooth domain to the original data domain c 1 (t) using the inverse transformation. The inverse transformation formula can be written as: Ideally, c 1 (t) is the first component MOC 1 (t) of original signal x(t) when c 1 (t) is a mono-component signal of instantaneous frequency with physical meaning.
(v) If c 1 (t) is not a mono-component signal of the instantaneous frequency with a physical meaning, we should repeat the above procedure m times until c m (t) is a mono-component signal of instantaneous frequency with physical meaning. In other words, in this situation, c m (t) will be the first component MOC 1 (t) of the original signal. (vi) This derived MOC 1 (t) is then subtracted from the original signal x(t), obtaining a new function r 1 (t). Regard r 1 (t) as a new original signal and repeat the procedure (i)~(v) n times until r n (t) is a monotonic function. Thus, x(t) may be decomposed into the sum of n MOCs and residue r n (t). i.e., The flowchart of the RS-LOD method is illustrated in Figure 5.

Application to Simulation Signal
In order to verify the effectiveness of the RS-LOD method and compare it with other methods, the same simulation signal is analyzed using RS-LOD, EMD, and LOD methods, respectively. In a comparison of the analysis processing, the end effects are not solved. It is to say that end effects may appear in all methods. Therefore, the severity of end effects can be used as an assessing index to compare the three methods. Additionally, the standard deviation (SD) criterion [7] is adopted as the iterative stopping condition, that is, by calculating the standard deviation of the results of two consecutive iterations, and comparing with the predetermined threshold to determine the iterative termination. The predetermined threshold is 0.5 in this paper.
Here we consider a multi-component AM-FM simulation signal x(t): where t ∈ [0, 6], the setting sampling frequency 1000 Hz, and the simulation signal consists of a modulated signal and two sine signals. The time-domain waveforms of x(t) and its three components x 1 (t), x 2 (t), and x 3 (t) are shown in Figure 6.
where t ∈ [0, 6], the setting sampling frequency 1000 Hz, and the simulation signal consists of a modulated signal and two sine signals. The time-domain waveforms of x(t) and its three components x1(t), x2(t), and x3(t) are shown in Figure 6. The simulation signal x(t) is decomposed by using the LOD, EMD, and RS-LOD method and the decomposition results are presented in Figures 7-9, respectively. It can be seen from Figures 7-9 that all three methods can decompose the original signal from high frequency to low frequency in turn, and the decomposed three components correspond to the x1(t), x2(t), and x3(t) of the original signal, respectively. However, it is also obvious found that there exist some differences in the results of their decomposition. In Figure 7, the decomposition result of the LOD is serrated, and the smoothness is poor (marked by a red circle). The RS-LOD decomposition result shown in Figure 9 (set p = 5) does not show this phenomenon, and because of the good smoothness, the cumulative error is small; thus the amplitude of the residue is obviously smaller than that of LOD. In Figure 8, due to the cubic spline function used by EMD, the end effects are obvious [20], and the amplitude of residue is also larger than that of RS-LOD because of the error of the end effects. In summary, we can see that the decomposition result of the RS-LOD is better than that of LOD and EMD as a whole.  The simulation signal x(t) is decomposed by using the LOD, EMD, and RS-LOD method and the decomposition results are presented in Figures 7-9, respectively. It can be seen from Figures 7-9 that all three methods can decompose the original signal from high frequency to low frequency in turn, and the decomposed three components correspond to the x 1 (t), x 2 (t), and x 3 (t) of the original signal, respectively. However, it is also obvious found that there exist some differences in the results of their decomposition. In Figure 7, the decomposition result of the LOD is serrated, and the smoothness is poor (marked by a red circle). The RS-LOD decomposition result shown in Figure 9 (set p = 5) does not show this phenomenon, and because of the good smoothness, the cumulative error is small; thus the amplitude of the residue is obviously smaller than that of LOD. In Figure 8, due to the cubic spline function used by EMD, the end effects are obvious [20], and the amplitude of residue is also larger than that of RS-LOD because of the error of the end effects. In summary, we can see that the decomposition result of the RS-LOD is better than that of LOD and EMD as a whole. Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 19    In addition, in order to quantitatively compare the decomposition results of the three methods, the assessing indicators, including the index of orthogonality (IO) [7], the number of sifting times, and the consuming time are used to evaluate the results. The three algorithms have been run 30 times using MATLAB R2017b to above simulation signal x(t) on the same computer, respectively, and all using SD < 0.5 as iterative stopping condition. The average values of three assessing indicators are calculated and counted, and the results are shown in Table 1.
It can be seen from Table 1, the IO value, which reflects the decomposition effects, the RS-LOD is the best of the three. Comparing the LOD method which uses linear transformation completely, the consuming time of the RS-LOD is a little longer, but they are both obviously shorter than that of the EMD method. Therefore, through a quantitative comparison, the RS-LOD method also has some advantages in terms of decomposition effects and computational efficiency.

Application to Rotating Machinery Fault Diagnosis
When gear or rolling bearing is operated with local faulty, the measured vibration signal generally is a multi-component AM-FM signal, and the high-frequency vibration components, such as natural frequency vibration and mesh frequency vibration, are modulated by shock pulse caused by local faults [27,28]. The RS-LOD method can decompose the multi-component AM-FM signal into a series of single-component signals, while performing the Hilbert transform on each singlecomponent signal can obtain an envelope signal with a physical meaning. Theoretically, the modulation information caused by fault can be extracted from the spectrum of these envelope signals. Therefore, RS-LOD is very suitable for extracting characteristic information from fault vibration signals of the rolling bearing and gear.
The main steps of extracting fault characteristics based on the RS-LOD method can be summarized as follows: Firstly, the vibration signal is decomposed into a serious of the MOC components by RS-LOD method. Secondly, the envelope signal of MOC component with fault characteristics is obtained by the Hilbert transform. Finally, the frequency spectrum of the envelope signal is analyzed for effectively extracting the fault characteristic frequency. In addition, in order to quantitatively compare the decomposition results of the three methods, the assessing indicators, including the index of orthogonality (IO) [7], the number of sifting times, and the consuming time are used to evaluate the results. The three algorithms have been run 30 times using MATLAB R2017b to above simulation signal x(t) on the same computer, respectively, and all using SD < 0.5 as iterative stopping condition. The average values of three assessing indicators are calculated and counted, and the results are shown in Table 1. It can be seen from Table 1, the IO value, which reflects the decomposition effects, the RS-LOD is the best of the three. Comparing the LOD method which uses linear transformation completely, the consuming time of the RS-LOD is a little longer, but they are both obviously shorter than that of the EMD method. Therefore, through a quantitative comparison, the RS-LOD method also has some advantages in terms of decomposition effects and computational efficiency.

Application to Rotating Machinery Fault Diagnosis
When gear or rolling bearing is operated with local faulty, the measured vibration signal generally is a multi-component AM-FM signal, and the high-frequency vibration components, such as natural frequency vibration and mesh frequency vibration, are modulated by shock pulse caused by local faults [27,28]. The RS-LOD method can decompose the multi-component AM-FM signal into a series of single-component signals, while performing the Hilbert transform on each single-component signal can obtain an envelope signal with a physical meaning. Theoretically, the modulation information caused by fault can be extracted from the spectrum of these envelope signals. Therefore, RS-LOD is very suitable for extracting characteristic information from fault vibration signals of the rolling bearing and gear.
The main steps of extracting fault characteristics based on the RS-LOD method can be summarized as follows: Firstly, the vibration signal is decomposed into a serious of the MOC components by RS-LOD method. Secondly, the envelope signal of MOC component with fault characteristics is obtained by the Hilbert transform. Finally, the frequency spectrum of the envelope signal is analyzed for effectively extracting the fault characteristic frequency.

Application to Rolling Bearing Fault Diagnosis
In order to verify the proposed method, the rolling bearing failure experiment was carried out on a rotating machinery fault test bench. The top view of the experimental device is shown in Figure 10. The adjustable speed motor in this system is a DC servo motor with a power of 600 W. The rolling bearings are installed on four bearing chocks, and both driving gear and driven gear are spur gears. The moment of the inertia of the load is 0.03 kg·m 2 . In this experiment, the tested bearing is the 6311-type rolling bearing, which was formed the inner-race and outer-race faults by laser cutting technology, and the slot width and depth are 0.15 mm and 0.13 mm, respectively. The vibration signals are collected by two piezoelectric acceleration transducers, which are horizontally and vertically fixed to the bearing housing. Figure 11 shows the time-domain waveform of vibration acceleration signal of the rolling bearing with outer-race fault. The sampling frequency is 4096 Hz, and the rotating frequency is f r = 25 Hz, and by calculating the outer-race fault characteristic frequency is f o = 76 Hz.

Application to Rolling Bearing Fault Diagnosis
In order to verify the proposed method, the rolling bearing failure experiment was carried out on a rotating machinery fault test bench. The top view of the experimental device is shown in Figure  10. The adjustable speed motor in this system is a DC servo motor with a power of 600 W. The rolling bearings are installed on four bearing chocks, and both driving gear and driven gear are spur gears. The moment of the inertia of the load is 0.03 kg·m 2 . In this experiment, the tested bearing is the 6311type rolling bearing, which was formed the inner-race and outer-race faults by laser cutting technology, and the slot width and depth are 0.15 mm and 0.13 mm, respectively. The vibration signals are collected by two piezoelectric acceleration transducers, which are horizontally and vertically fixed to the bearing housing. Figure 11 shows the time-domain waveform of vibration acceleration signal of the rolling bearing with outer-race fault. The sampling frequency is 4096 Hz, and the rotating frequency is fr = 25 Hz, and by calculating the outer-race fault characteristic frequency is fo = 76 Hz  The vibration signal shown in Figure 11 is decomposed by the RS-LOD method. The value of the pole parameter p will affect decomposition results and decomposition time; therefore, we need to determine the value of p firstly. Table 2 shows the number of sifts required to identify each MOC. It can be seen that the number of sifts required for each MOC component and the total number of sifts for per decomposition for a range of p values. When p takes relatively smaller value at a certain range, the total number of sifts required is less, and the number of MOCs obtained is also less. From these two standards, when p takes the range of 0 to 10, it is possible to get higher decomposition efficiency.

Application to Rolling Bearing Fault Diagnosis
In order to verify the proposed method, the rolling bearing failure experiment was carried out on a rotating machinery fault test bench. The top view of the experimental device is shown in Figure  10. The adjustable speed motor in this system is a DC servo motor with a power of 600 W. The rolling bearings are installed on four bearing chocks, and both driving gear and driven gear are spur gears. The moment of the inertia of the load is 0.03 kg·m 2 . In this experiment, the tested bearing is the 6311type rolling bearing, which was formed the inner-race and outer-race faults by laser cutting technology, and the slot width and depth are 0.15 mm and 0.13 mm, respectively. The vibration signals are collected by two piezoelectric acceleration transducers, which are horizontally and vertically fixed to the bearing housing. Figure 11 shows the time-domain waveform of vibration acceleration signal of the rolling bearing with outer-race fault. The sampling frequency is 4096 Hz, and the rotating frequency is fr = 25 Hz, and by calculating the outer-race fault characteristic frequency is fo = 76 Hz  The vibration signal shown in Figure 11 is decomposed by the RS-LOD method. The value of the pole parameter p will affect decomposition results and decomposition time; therefore, we need to determine the value of p firstly. Table 2 shows the number of sifts required to identify each MOC. It can be seen that the number of sifts required for each MOC component and the total number of sifts for per decomposition for a range of p values. When p takes relatively smaller value at a certain range, the total number of sifts required is less, and the number of MOCs obtained is also less. From these two standards, when p takes the range of 0 to 10, it is possible to get higher decomposition efficiency. The vibration signal shown in Figure 11 is decomposed by the RS-LOD method. The value of the pole parameter p will affect decomposition results and decomposition time; therefore, we need to determine the value of p firstly. Table 2 shows the number of sifts required to identify each MOC. It can be seen that the number of sifts required for each MOC component and the total number of sifts for per decomposition for a range of p values. When p takes relatively smaller value at a certain range, the total number of sifts required is less, and the number of MOCs obtained is also less. From these two standards, when p takes the range of 0 to 10, it is possible to get higher decomposition efficiency.  0  5  7  6  7  ----25  0.5  5  6  6  8  8  ---33  1  5  7  6  7  8  8  --41  2  5  7  6  6  8  7  --39  5  5  7  6  6  8  7  8  -47  10  5  7  7  6  7  7  7  -46  20  5  7  8  8  8  8  7  8  59  50  5  7  8  7  7  8  8  8  58   Table 3 lists the IO value and decomposition time for different p values. It can be seen that with the increase of p, the IO value decreases gradually, and the decomposition time increases during the fluctuation. However, when p is greater than 10, the decrease of IO value is not obvious and undergoes gradual stabilization, but the decomposition time increases obviously. Therefore, consider the decomposition effects and decomposition time to determine the p value is 10 in this decomposition.  When p = 10, the decomposition results are shown in Figure 12, and the high-frequency components MOC 1 (t) and MOC 2 (t) are shown in Figure 13. It can be seen that seven MOC components and one residue are obtained from high-frequency to low-frequency, and the decomposition effects are ideal. According to the rolling bearing fault vibration mechanism, the fault vibration signal of the rolling bearing is the modulation signal with the natural frequency of a rolling bearing system as the carrier frequency and the fault characteristic frequency as the modulation frequency. Therefore, the working condition of the rolling bearing can be judged by the envelope spectrum of the fault vibration signal. Additionally, the natural frequency of the rolling bearing system usually exists in the high-frequency components. Hence, by applying the Hilbert transform to the MOC 1 (t), and the obtained Hilbert envelope spectrum, as shown in Figure 14, it can be easily seen that, except for the spectral line of rotating frequency (25 Hz), there are also obvious spectral lines of the outer-race fault characteristic frequency (76 Hz) and its double (152 Hz) can be found in the Hilbert envelope spectrum. Therefore, the fault of the outer race of the rolling bearing can be judged, which is consistent with the actual situation, and shows the effectiveness of the proposed method.     Figure 15 shows the structure of the No.3 fan gearbox of a petrochemical company, in which the teeth number of gears 1, 2, 3, and 4 are 11, 31, 23, and 53, respectively. The rotational speed of shaft I is nI = 980 r/min, namely, the rotating frequency of shaft I is about fI = 16.3 Hz, the rotating frequency of shaft II is about fII = 5.8 Hz, and the rotating frequency of shaft III is about fIII = 2.5 Hz. One day, the online monitoring system found that the vibration of the no.3 fan gearbox was abnormal. The root means square value of the vibration acceleration signal (on monitoring point A shown in Figure 15) was relatively bigger than the threshold value, but the fan running was not affected, and the date was close to the normal shutdown maintenance date, so keep observation. During the subsequent shutdown maintenance, the fan gearbox was examined. It was found that four gears of gearbox had different degrees of wear and tooth surface glue, especially the wear of gear 1 and gear 2 were very serious, which can be shown in Figure 16. In order to verify that Hilbert envelope spectrum based on RS-LOD can accurately identify the fault condition of the gearbox, a vibration acceleration signal (the sampling frequency is 1024 Hz) measured on the same day is selected for analysis, and which time-domain waveform and amplitude spectrum are as shown in Figures 17 and 18. By preliminary analysis, the main frequency components include the mesh frequency of gear 1 and gear 2, the mesh frequency of gear 3 and gear 4, and many other unknown frequencies. Moreover, the sidebands of the gear mesh frequency are very difficult to distinguish. Therefore, it is hard to estimate whether and where the fault has occurred in the fr fo 2fo Figure 14. The Hilbert envelope spectrum of the MOC 1 (t) obtained by the RS-LOD. Figure 15 shows the structure of the No.3 fan gearbox of a petrochemical company, in which the teeth number of gears 1, 2, 3, and 4 are 11, 31, 23, and 53, respectively. The rotational speed of shaft I is n I = 980 r/min, namely, the rotating frequency of shaft I is about f I = 16.3 Hz, the rotating frequency of shaft II is about f II = 5.8 Hz, and the rotating frequency of shaft III is about f III = 2.5 Hz. One day, the online monitoring system found that the vibration of the no.3 fan gearbox was abnormal. The root means square value of the vibration acceleration signal (on monitoring point A shown in Figure 15) was relatively bigger than the threshold value, but the fan running was not affected, and the date was close to the normal shutdown maintenance date, so keep observation. During the subsequent shutdown maintenance, the fan gearbox was examined. It was found that four gears of gearbox had different degrees of wear and tooth surface glue, especially the wear of gear 1 and gear 2 were very serious, which can be shown in Figure 16.  Figure 15 shows the structure of the No.3 fan gearbox of a petrochemical company, in which the teeth number of gears 1, 2, 3, and 4 are 11, 31, 23, and 53, respectively. The rotational speed of shaft I is nI = 980 r/min, namely, the rotating frequency of shaft I is about fI = 16.3 Hz, the rotating frequency of shaft II is about fII = 5.8 Hz, and the rotating frequency of shaft III is about fIII = 2.5 Hz. One day, the online monitoring system found that the vibration of the no.3 fan gearbox was abnormal. The root means square value of the vibration acceleration signal (on monitoring point A shown in Figure 15) was relatively bigger than the threshold value, but the fan running was not affected, and the date was close to the normal shutdown maintenance date, so keep observation. During the subsequent shutdown maintenance, the fan gearbox was examined. It was found that four gears of gearbox had different degrees of wear and tooth surface glue, especially the wear of gear 1 and gear 2 were very serious, which can be shown in Figure 16. In order to verify that Hilbert envelope spectrum based on RS-LOD can accurately identify the fault condition of the gearbox, a vibration acceleration signal (the sampling frequency is 1024 Hz) measured on the same day is selected for analysis, and which time-domain waveform and amplitude spectrum are as shown in Figures 17 and 18. By preliminary analysis, the main frequency components include the mesh frequency of gear 1 and gear 2, the mesh frequency of gear 3 and gear 4, and many other unknown frequencies. Moreover, the sidebands of the gear mesh frequency are very difficult to distinguish. Therefore, it is hard to estimate whether and where the fault has occurred in the fr fo 2fo Figure 15. The no.3 fan gearbox structure schematic diagram.

Application to Fan Gearbox Fault Diagnosis
In order to verify that Hilbert envelope spectrum based on RS-LOD can accurately identify the fault condition of the gearbox, a vibration acceleration signal (the sampling frequency is 1024 Hz) measured on the same day is selected for analysis, and which time-domain waveform and amplitude spectrum are as shown in Figures 17 and 18. By preliminary analysis, the main frequency components include the mesh frequency of gear 1 and gear 2, the mesh frequency of gear 3 and gear 4, and many other unknown frequencies. Moreover, the sidebands of the gear mesh frequency are very difficult to distinguish. Therefore, it is hard to estimate whether and where the fault has occurred in the gearbox from the amplitude spectrum in Figure 18. Further, the envelope demodulation technology based on the RS-LOD is used to identify the condition of the gearbox. According to the manner in Section 5.1, we can determine the value of p in this decomposition is 9. A total of ten MOC components are obtained by RS-LOD decomposition, and the first five MOC components are shown in Figure 19. It can be seen that the decomposition effects are good, and there are no obvious distortions. The first three MOC components are further selected for Hilbert envelope analysis, and the results are shown in Figure 20, respectively.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 19 gearbox from the amplitude spectrum in Figure 18. Further, the envelope demodulation technology based on the RS-LOD is used to identify the condition of the gearbox. According to the manner in Section 5.1, we can determine the value of p in this decomposition is 9. A total of ten MOC components are obtained by RS-LOD decomposition, and the first five MOC components are shown in Figure 19. It can be seen that the decomposition effects are good, and there are no obvious distortions. The first three MOC components are further selected for Hilbert envelope analysis, and the results are shown in Figure 20, respectively.     gearbox from the amplitude spectrum in Figure 18. Further, the envelope demodulation technology based on the RS-LOD is used to identify the condition of the gearbox. According to the manner in Section 5.1, we can determine the value of p in this decomposition is 9. A total of ten MOC components are obtained by RS-LOD decomposition, and the first five MOC components are shown in Figure 19. It can be seen that the decomposition effects are good, and there are no obvious distortions. The first three MOC components are further selected for Hilbert envelope analysis, and the results are shown in Figure 20, respectively.     gearbox from the amplitude spectrum in Figure 18. Further, the envelope demodulation technology based on the RS-LOD is used to identify the condition of the gearbox. According to the manner in Section 5.1, we can determine the value of p in this decomposition is 9. A total of ten MOC components are obtained by RS-LOD decomposition, and the first five MOC components are shown in Figure 19. It can be seen that the decomposition effects are good, and there are no obvious distortions. The first three MOC components are further selected for Hilbert envelope analysis, and the results are shown in Figure 20, respectively.     In Figure 20a, it can be seen that there is a spectral line with high amplitude at the double frequency of rotating frequency fI. In Figure 20b, there are several obvious spectral lines, which at the 1, 2, 3, 4 times of rotating frequency fI, and at the 1, 2 times of rotating frequency fII. In Figure 20c. It can be seen at the rotating frequency fI and 1, 2 times of rotating frequency fIII there are clear spectral lines. According to the distribution and amplitude of above spectral lines, it can be determined that there are certain degree faults in all four gears on three shafts, and the faults on gear 1 and gear 2 are more serious, while the faults on gear 3 and gear 4 are relatively slight. It is consistent with the actual gearbox check, which confirmed the effectiveness of the proposed diagnosis method.  In Figure 20a, it can be seen that there is a spectral line with high amplitude at the double frequency of rotating frequency fI. In Figure 20b, there are several obvious spectral lines, which at the 1, 2, 3, 4 times of rotating frequency fI, and at the 1, 2 times of rotating frequency fII. In Figure 20c. It can be seen at the rotating frequency fI and 1, 2 times of rotating frequency fIII there are clear spectral lines. According to the distribution and amplitude of above spectral lines, it can be determined that there are certain degree faults in all four gears on three shafts, and the faults on gear 1 and gear 2 are more serious, while the faults on gear 3 and gear 4 are relatively slight. It is consistent with the actual gearbox check, which confirmed the effectiveness of the proposed diagnosis method.

Discussion
Through the analysis of simulation signals and experimental signals, we know that the RS-LOD method inherits the strengths of the LOD method, such as fast calculation speed, no obvious end effects, stable decomposition results, and much slighter mode mixing phenomena etc. At the same time, the rational spline function with a pole parameter p is introduced to overcome the problems of MOC component poor smoothness and distortion. For different signals, the optimal local mean curve can be obtained flexibly by adjusting the pole parameter p, reducing the rounding error, and obtaining the optimal component. However, the optimal selection range of pole parameter p of rational spline function is estimated through multiple numerical experiments, i.e., the selection of f III 2f III f I Figure 20. The Hilbert envelope spectrum of the MOC components: (a) Hilbert envelope spectrum of MOC 1 (t); (b) Hilbert envelope spectrum of MOC 2 (t); (c) Hilbert envelope spectrum of MOC 3 (t).
In Figure 20a, it can be seen that there is a spectral line with high amplitude at the double frequency of rotating frequency f I . In Figure 20b, there are several obvious spectral lines, which at the 1, 2, 3, 4 times of rotating frequency f I , and at the 1, 2 times of rotating frequency f II . In Figure 20c. It can be seen at the rotating frequency f I and 1, 2 times of rotating frequency f III there are clear spectral lines. According to the distribution and amplitude of above spectral lines, it can be determined that there are certain degree faults in all four gears on three shafts, and the faults on gear 1 and gear 2 are more serious, while the faults on gear 3 and gear 4 are relatively slight. It is consistent with the actual gearbox check, which confirmed the effectiveness of the proposed diagnosis method.

Discussion
Through the analysis of simulation signals and experimental signals, we know that the RS-LOD method inherits the strengths of the LOD method, such as fast calculation speed, no obvious end effects, stable decomposition results, and much slighter mode mixing phenomena etc. At the same time, the rational spline function with a pole parameter p is introduced to overcome the problems of MOC component poor smoothness and distortion. For different signals, the optimal local mean curve can be obtained flexibly by adjusting the pole parameter p, reducing the rounding error, and obtaining the optimal component. However, the optimal selection range of pole parameter p of rational spline function is estimated through multiple numerical experiments, i.e., the selection of pole parameter p is subjective and lack of adaptivity. Therefore, how to determine the optimal pole parameter p automatically is a problem that needs further study. In addition, other interpolation methods [29,30] can also be considered to improve the decomposition effects in the future.

Conclusions
In order to solve the problems of MOC component poor smoothness and distortion caused by the piecewise linear transformation in the LOD method, the rational spline interpolation is used to replace the piecewise linear transformation and the RS-LOD method is proposed in this paper. Then, the RS-LOD method is used for rotating machinery fault feature extraction. The following conclusions are obtained: 1.
The analysis of the simulation signal shows that although the consuming time of the RS-LOD method is a little longer than the original LOD method, it solves the problems of MOC component poor smoothness and distortion in the original LOD method, and overall analysis effects are also better than that of EMD and LOD method.

2.
The analyses of fault vibration signals of rolling bearing and fan gearbox show that the envelope spectrum based on RS-LOD can effectively extract the characteristic fault information of the corresponding vibration signal, which provides a new effective way to extract the vibration signal feature.