An improved spline-local mean decomposition and its application to vibration analysis of rotating machinery with rub-impact fault

2.1. LMD algorithm

The fundamental objective of LMD algorithm is to decompose amplitude and frequency modulated signals into a small set of PFs. The LMD process essentially involves progressively separating a frequency modulated signal from an amplitude modulated signal [8]. This separation procedure can be briefly described as follows:

(1) Given a random signal $x\left(t\right)$, set the counter $k=1$.

(2) Identify all the local extrema of $x\left(t\right)$, and then respectively calculate the local mean $m_i$ and local magnitude $a_i$ between each two successive extrema $n_i$ and $n_{i+1}$ by:

(3) All the local means can be plotted as different straight line segments, and these line segments produce an original mean function. Then, this function is smoothed using moving averaging to form a smoothly varying continuous local mean function $m\left(t\right)$. An envelope estimate function $a\left(t\right)$ which also possesses the smoothly varying continuous property can be produced through implementing the same procedures on the local magnitudes.

(4) Subtract $m\left(t\right)$ from $x\left(t\right)$:

and subsequently, $h\left(t\right)$ is divided by $a\left(t\right)$ to achieve amplitude demodulation:

(5) If $s\left(t\right)$ is a purely frequency modulated signal, go to next step. Otherwise, treat $s\left(t\right)$ as $x\left(t\right)$ and repeat the steps (2)-(4) until $s\left(t\right)$ satisfies the prescribed criterion that $a\left(t\right)$ is equal to 1 at any time. In the meantime, the corresponding IA can be given by:

where $k$ denotes the number of the current PF, and $n$ denotes the number of iterations for computing the $k$th PF.

(6) Multiplying $s\left(t\right)$ by $a_k\left(t\right)$ gives the $k$th PF:

and the corresponding IF can be calculated by:

(7) Remove ${PF}_k\left(t\right)$ from $x\left(t\right)$:

(8) Treat $u_k\left(t\right)$ as $x\left(t\right)$ and repeat the steps (2)-(7). This iterative separation process will continue for $q$ times until $u_q\left(t\right)$ becomes a monotonic function or contains no more oscillations. Finally, several PFs and a residue are progressively produced, and accordingly the original signal $x\left(t\right)$ can be reconstructed by:

2.2. Deficiency of moving averaging process

LMD is a self-adaptive signal analysis method and can produce the entire time-frequency representation of a non-stationary signal. However, LMD has an inherent deficiency that the innermost loop in LMD is time-consuming and the smoothed results are inaccurate, especially for processing the complicated signals. There are two key factors leading to this troublesome problem. First, the smoothed results of moving averaging will be markedly different if the length of moving window is different; second, the moving averaging process needs to be repeated many times. In fact, after the local means and magnitudes are got by Eqs. (1) and (2), respectively, the calculation process of both local mean and envelope estimate functions can be described as follows [8]:

(1) Denote all the data points (i.e., the local means or magnitudes) as an initial data set $\left\{w\left(i\right);i=1,2,\dots ,N\right\}$, where $N$ is the number of data points.

(2) The $i$th updated data element is written as:

where $2h+1$ represents the window width (i.e., the length of moving window), and it should be an odd number, otherwise minus 1. Under five-point averaging mode, all the updated data elements can be given by:

If any two successive data elements are equal, treat the updated data series $\left\{w_s\left(i\right);i=1,2,\dots ,N\right\}$, as the initial data set $\left\{w\left(i\right);i=1,2,\dots ,N\right\}$, and repeat the above moving averaging process. This iterative computation process will be stopped until the updated data series becomes a smoothly varying continuous function.

Because the length of the moving window has a great influence on the convergence rate and computational accuracy of moving averaging algorithm, it is very significant to obtain an appropriate window width. Unfortunately, it is markedly difficult to identify an optimal window width for a complicated non-stationary signal, which may produce an unexpected result that LMD performs inefficiently and inaccurately in processing this kind of signals, especially in analyzing the mechanical vibration signals. On the other hand, in the original LMD method, computing the smoothly varying continuous local mean and envelope estimate functions must depend on the moving averaging algorithm.

3.1. Cubic spline interpolation

Given an interval $[a,\mathrm{}b]$ composed of $n$ ordered disjoint subintervals $[x_{i-1},x_i]$ with:

and the corresponding function values $y_i=f\left(x_i\right)$$(i=\mathrm{0,1},\dots ,n)$ if a function $S\left(x\right)$ satisfies the following conditions:

(1) On each subinterval $\left[x_{i-1},x_i\right]$$(i=\mathrm{0,1},\dots ,n)$$S\left(x\right)$ is a cubic polynomial;

(2) $\left(x_i\right)=y_i$$(i=\mathrm{0,1},\dots ,n)$;

(3) On the interval $[a,\mathrm{}b]$, if the second derivative of $S\left(x\right)$, $\ddot{S}\left(x\right)$, is a continuous function, $S\left(x\right)$ is called a cubic spline interpolation of the function $y=f\left(x\right)$ on the interval $[a,\mathrm{}b]$.

We hypothesize that the second derivative values at each $x_i$ are $\ddot{S}\left(x_i\right)=M_i$$\left(i=0,1,\dots ,n\right),$ where $M_i$ are $n+1$ unknown parameters. $S\left(x\right)$ is a piecewise cubic polynomial, so $\ddot{S}\left(\mathrm{x}\right)$ is a piecewise linear function and can be expressed as:

where $h_i=x_i-x_{i-1}$. Since the second derivative $\ddot{S}\left(\mathrm{x}\right)$ is a continuous function, we can operate the quadratic integration on Eq. (13), and obtain:

where $b_i$ and $c_i$ are integration constants. According to the interpolation conditions: $S\left(x_{i-1}\right)=y_{i-1}$ and $S\left(x_i\right)=y_i$, we can obtain:

that is:

Substituting Eq. (16) into Eq. (14), $S\left(x\right)$ on the subinterval $\left[x_{i-1},x_i\right]$$\left(i=\mathrm{0,1},\dots ,n\right)$, can be expressed as:

Subsequently, the first derivative of Eq. (17) can be derived:

Setting $x=x_i$, we can obtain the left derivative of $S\left(x\right)$ at $x_i$:

Similarly, setting $x=x_{i-1}$, we can obtain the right derivative of $S\left(x\right)$ at $x_{i-1}$:

and then:

Furthermore, according to the continuous property of $\dot{S}\left(x\right)$, we obtain:

Substituting Eqs. (19) and (21) into Eq. (22) gives:

Multiplying $6/\left(h_i+h_{i+1}\right)$ to Eq. (23) gives:

Setting:

Eq. (24) can be written as:

Since Eq. (25) contains $n-1$ equations, the $n+1$ unknown parameters $M_i$$(i=0,1,\dots ,n)$ still cannot be identified. In order to solve Eq. (25), we need to add other two equations.

Assuming that the first derivative values at endpoints $a$ and $b$ are $\dot{f}\left(a\right)$ and $\dot{f}\left(b\right)$, which are known constants, and setting ${\dot{y}}_0=\dot{f}\left(a\right)$ and ${\dot{y}}_n=\dot{f}\left(b\right)$, we obtain:

Substituting Eqs. (19) and (21) into Eq. (26) gives:

i.e.:

where:

Combining Eqs. (25) and (28) results in:

and Eq. (29) can be further expressed as:

Solving Eq. (30) by the pursuit method based on LU decomposition, and then substituting the obtained $M_i$$(i=0,1,\dots ,n)$ into Eq. (17), a cubic spline interpolation can be produced.

Compared with the moving averaging process, the cubic spline interpolation needs to be performed only once in the process of computing local mean and envelope estimate functions, while the moving averaging process needs to be repeated many times to achieve the smoothing objective. Moreover, cubic spline interpolation possesses low interpolation error and excellent convergence property. Therefore, cubic spline interpolation consumes less time and obtains more accurate local mean and envelope estimate functions than moving averaging process.

3.2. Improved scheme for computing local mean and envelope estimate functions

In order to overcome the weakness of the LMD method, a new SLMD method is developed in this part. On the basis of the analyses given in subsections 2.2 and 3.1, the moving averaging process in the original LMD algorithm is replaced with the following procedures to produce local mean and envelope estimate functions.

(1) After all extrema of the signal are identified, respectively connect all the local maxima and minima with two different cubic spline lines to form the upper envelope $E_u\left(t\right)$ and lower envelope $E_l\left(t\right)$.

(2) The local mean function $m\left(t\right)$ and local envelope estimate function $a\left(t\right)$ can be given by:

(3) Continue the rest of the procedures in the original LMD algorithm.

In fact, the above improvements only substitute the cubic spline interpolation for the moving averaging process, but it avoids a lowly efficient iterative computation process. Although some boundary distortions, overshoots and undershoots exist in the envelope estimate functions [5], overshoot and undershoot phenomena will be greatly weakened due to the arithmetic means computed by Eq. (31). As for the boundary distortions, a signal extending technique based on self-adaptive waveform matching is proposed to eliminate this problem.

3.3. Signal extending based on self-adaptive waveform matching

Signal extending is a very effective scheme for overcoming the boundary distortion in the process of signal processing and analysis [18, 19]. There are several signal extending methods, but a most suitable extending method needs to be specified according to the characteristics of the vibration signal. Usually, the vibration of a running rotor is a quasi-periodic oscillation behavior [20], so the similar waveforms are bound to arise in the entire vibration signal. Based on this characteristic, a self-adaptive signal extending method is developed to process the raw vibration signal. Here, only left extending is mainly described because of the same principle for right extending, and the left extending process can be depicted as follows [21, 22]:

(1) Suppose a given signal $\left(t\right)$$(t=1,2,\dots ,L)$, where $L$ is the number of sampled points.

(2) Identify all the local extrema of $x\left(t\right)$, which are denoted by $n_i$$(i=1,2,\dots ,T)$. The corresponding time series of these extrema in $x\left(t\right)$ is denoted by $t_i$$(i=1,2,\dots ,T)$ and the data points between $x\left(1\right)$ and $n_2$ are named characteristic waveform (denoted by $p\left(t\right)$), in which the number of data points is denoted by $k$.

(3) Find out the time interval $\left[\left(t_i-t_1\right)+1,t_i+(k-t_1)\right]$$(i=2,3,\dots ,T-1)$ of each matching waveform (denoted by $p_i\left(t\right)$) in $x\left(t\right)$ via the rest of extrema with the same local characteristic as $n_1$ one by one, and compute the corresponding waveform matching error:

where $i$ represents the sequence number of the extrema.

(4) The minimal waveform matching error is denoted by:

and the extremum point with the same local characteristic as $n_1$ in this matching waveform is denoted by $n_e$. The data points between $n_1$ and $n_e$ are copied and put on the left of $n_1$. Thus, left extending of the sampled vibration signal is accomplished. In this way, right extending also can be implemented after the left extended signal is reversed. Finally, reversing the right extended signal results in a boundary extended signal with the normal sequence.

When subsections 3.2 and 3.3 are orderly combined, an improved SLMD method is formed and its flowchart is shown in Fig. 1. It is well known that the LMD algorithm has a three-layer cyclic structure, while it can be clearly seen that the improved SLMD algorithm is only with a two-layer cyclic structure, so improved SLMD is relatively simple. For convenience, improved SLMD is stilled named SLMD in the following texts.

Fig. 1Flowchart of the improved YSLMD algorithm

4.1. Iteration number and time consumption

Since the cyclic moving averaging process is replaced with the relatively simple cubic spline interpolation, the iteration number and time consumption of decomposing the given signal $x\left(t\right)$ via SLMD should be different from those via LMD. And therefore, their respective iteration number and time consumption are compared in this part.

Fig. 3Decomposed results of the simulative signal: a) via LMD; b) via SLMD

(a)

(b)

Although SLMD is different from LMD in a few computation steps, both they decompose $x\left(t\right)$ into two PFs and a residue, and we can hardly find, only by visual observation, any obvious differences between the two decomposed results shown in Fig. 3. However, when $x\left(t\right)$ was decomposed using LMD and SLMD, respectively, we could perceive the difference that the running time of LMD is obviously longer than that of SLMD. The actual iteration number and time consumption of decomposing $x\left(t\right)$ via the two methods are shown in Table 1.

From Table 1, it can be seen that the iteration number of computing the PFs and the time consumption of operating the whole decomposition process of SLMD are less and shorter than those of LMD, respectively. These results indicate that SLMD actually decrease the iteration number and running time of the decomposition process, which demonstrates that the improvements of SLMD really plays an important role in improving the efficiency of computing the local mean and envelope estimate functions.

4.2. Errors of each PF and residue

Two sets of similar decomposed results respectively shown in Figs. 3(a), (b) are compared and investigated in this part. For each ${PF}_i\left(t\right)$, we subtract it from its corresponding real component $x_i\left(t\right)$ to form a component error:

Consequently, four different error curves can be produced via Eq. (35). They are classified into two groups and shown in Fig. 4. The results obviously show that there actually are some errors between the PFs and their corresponding real components. Moreover, compared with LMD, SLMD only produces relatively small errors. In other words, each PF derived by SLMD is more accurate than that derived by LMD. As for the two residues, nothing has been done to them and they are shown in Fig. 5. It can be seen that the residue produced by SLMD is closer to a zero vector than that produced by LMD.

Fig. 4Errors of the two sets of PFs: a) e1(t); b) e2(t)

(a)

(b)

Fig. 5Two residue u2(t)

In fact, since SLMD can produce more accurate PFs than LMD and the residue is equal to the original signal subtracting the corresponding PFs, the residue resulting from SLMD should be smaller than that resulting from LMD. Table 2 shows the integrals of the absolute errors and residues produced by LMD and SLMD, respectively. Each integral value is computed by:

where $e\left(t\right)$ represents the errors $e_1\left(t\right)$, $e_2\left(t\right)$ and the residue $u_2\left(t\right)$. From Table 2, we can see that the integral values of three error functions (including the residue) resulting from SLMD are obviously smaller than those resulting from LMD. These results indicate that SLMD can more accurately decompose the non-stationary signal than LMD.

4.3. Accuracy of IF and IA

IF and IA are the important features of a non-stationary signal, and they can be obtained along with LMD and SLMD processes. Fig. 6 shows the IF of ${PF}_1\left(t\right)$ and Fig. 7 shows the IA of ${PF}_1\left(t\right)$; here, ${PF}_1\left(t\right)$ represents two different PF components derived by LMD and SLMD, respectively.

Fig. 6IF of PF1(t)

(a)

(b)

Fig. 7IA of PF1(t)

(a)

(b)

The two IF curves are both distorted in the whole range, but the distortions of IF in Fig. 6(a) is much severer than that in Fig. 6(b). This phenomenon illustrates that IF actually becomes more accurate because of the improvements in SLMD. Moreover, the two IA curves plotted in Fig. 7 show that the IA computed by LMD also suffers from some obvious distortions, while the IA computed by SLMD is very accurate and there is almost no distortion in the whole range. This result indicates that the SLMD method actually improves the accuracy of IA. Through comparing the two IFs and IAs, respectively, we find that SLMD actually can produce more accurate instantaneous characteristics of a non-stationary signal than LMD.

4.4. Energy error

In accordance with the basic principle of signal decomposition, any signal $x\left(t\right)$ can be expressed as:

where ${\alpha }_k\left(t\right)$ is the coefficient, and $g_k\left(t\right)$ is the basic function. So, the energy of $x\left(t\right)$ can be determined by:

If $g_k\left(t\right)$$(k=\mathrm{1,2},\dots ,q)$ is an orthogonal basis set, the energy $E_x$ can be written as:

On the other hand, respective PFs derived by LMD and SLMD are almost orthogonal to each other, i.e., any two different PFs satisfy the following condition [23]:

where $〈\bullet ,\bullet 〉$ represents the inner product operator. Therefore, each component ${\alpha }_k\left(t\right)g_k\left(t\right)$ is essentially a PF, i.e.:

Substituting Eq. (41) into Eq. (39) gives:

It should be noted that the residue $u_q\left(t\right)$ is neglected in Eq. (42). Generally, the better orthogonality a group of PFs possess, the smaller the error existing in Eq. (42) becomes. In other words, the error between two sides of the approximately equal sign “$\approx$” in Eq. (42) can reflect the performance of LMD and SLMD to some extent. In fact, when a non-stationary signal is decomposed by any method, the signal energy computed via the decomposed components will change more or less. Therefore, the energies of the three groups of signal sets are calculated through Eq. (42), respectively. For clarity and convenience, the energy of the original signal is denoted by $E_0$, and the energies of the PFs derived by LMD and SLMD are denoted by $E_1$ and $E_2$, respectively. The three computed energy values are shown in Table 3. From Table 3, it can be easily found that the value of $E_0$ is the largest and the value of $E_1$ is the smallest. Although $E_1$ and $E_2$ are not equal to $E_0$, both they are very close to $E_0$. However, the energy error between $E_1$ and $E_0$ is 0.0872, while the energy error between $E_2$ and $E_0$ is only 0.0100. The result indicates that the PFs derived by SLMD make the error in Eq. (42) smaller than those derived by LMD, which verifies that SLMD actually can produce more accurate PFs than LMD.

It is clearly demonstrated through the above comparative study in four aspects that SLMD actually outperforms LMD in computational accuracy and efficiency. The improvements of SLMD play an important role in computing the local envelope and mean functions, and consequently, the derived PFs, IA and IF all become more accurate. SLMD not only reduces the running time of decomposing a signal, but also makes the instantaneous characteristics of a no-stationary signal more accurate. Therefore, SLMD indeed has a more powerful ability in processing the complicated non-stationary signal than LMD.