A compound interpolation envelope local mean decomposition and its application for fault diagnosis of reciprocating compressors

Mechanical Systems and Signal Processing 110 (2018) 273–295

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

A compound interpolation envelope local mean decomposition and its application for fault diagnosis of reciprocating compressors Zhao Haiyang a,b, Wang Jindong a,⇑, Jay Lee b, Li Ying a a b

Mechanical Science and Engineering Institute, Northeast Petroleum University, Daqing 163318, China NSF I/UCR Center for Intelligent Maintenance Systems, Department of Mechanical Engineering, University of Cincinnati, Cincinnati, OH 45221, USA

a r t i c l e

i n f o

Article history: Received 29 November 2017 Received in revised form 12 February 2018 Accepted 15 March 2018

Keywords: LMD MPCHI CSI Reciprocating compressor Fault diagnosis Bearing clearance

a b s t r a c t Local mean decomposition (LMD) is a self-adaptive analysis method which can decompose a signal into a set of product function (PF) components, and the construction of envelope function plays an important role in the accuracy of its PF components. According to the local nonstationary characteristics of vibration signals, a compound interpolation envelope (CIE) LMD was proposed through a novel envelope construction method. By defining an nonstationary coefficient to evaluate the local nonstationary characteristics of vibration signals, an compound envelope construction method which use Monotonic Piecewise Cubic Hermite Interpolation (MPCHI) for nonstationary part and Cubic spline interpolation (CSI) for smooth part was proposed, and an CIE LMD algorithm was give based on the novel envelope construction method. A numerical example simulation was conducted to verify the performance of CIE LMD, results indicated that CIE LMD outperforms three other methods. The CIE LMD was employed to diagnose the oversized bearing clearance fault in reciprocating compressor, and the envelope frequency spectrum of PF component gives a more significant peak of fault frequency than that of original signal, which further indicates that this proposed method is competent for the diagnosis of reciprocating compressor oversized bearing clearance fault. Ó 2018 Elsevier Ltd. All rights reserved.

1. Introduction In recent years, advancements in science and technology have led to the rapidly development of mechanical systems, and machineries have become more complicated in structure and functions. Therefore, machinery fault diagnosis is becoming increasingly important in the field of process monitoring due to greater demands for mechanical systems that provide higher performance, safety, and reliability [1]. Many sorts of modern machinery are often composed of mechanical, electrical and hydraulic systems. Their dynamic responses are a complex mixture of various dynamic phenomena, including mechanical vibrations, electrical oscillations, hydraulic fluctuations, their coupling effects, and the dynamic responses to external environmental excitations [2]. As a result, machinery data becomes more complicated. Feature extraction is the sequential task to extract intrinsic information of a targeted machinery and transform it into an understandable structure for further usage.

⇑ Corresponding author at: Mechanical Science and Engineering Institute, Northeast Petroleum University, Daqing 163318, China. E-mail addresses: [email protected] (Z. Haiyang), [email protected] (W. Jindong). https://doi.org/10.1016/j.ymssp.2018.03.035 0888-3270/Ó 2018 Elsevier Ltd. All rights reserved.

274

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

The key of this process is to find an effective signal analysis method to extract useful information implicitly existing in the machinery field data [3]. Reciprocating compressors are widely used in the gas storage, gas transport and petrochemical industries, and their operational safety has become the significant subject of the intensive investigation in recent years [4–6]. As a common used data source, vibration signals are easy to acquire and contain abundant operational state information [7,8], so it is highly effective in health monitoring and fault diagnosis of machinery and equipment such as gearbox and bearing under certain conditions [9,10]. However, the vibration signals of reciprocating compressor have characteristics of nonlinearity, nonstationarity and multi-component coupling due to the factors such as clearance, nonlinear stiffness of bearings and the unbalanced and timevarying forces of components, and those traditional techniques successful applied in rotating machinery fault diagnosis field, may fail to extract effective features from vibration signals of reciprocating compressors. Therefore, an effective and accurate fault diagnosis method is required to help increase the reliability and reduce the maintenance costs of reciprocating compressors [11]. Time–frequency analysis can identify the signal frequency components, reveals their time variant features, and is an effective tool to extract machinery health information contained in nonstationary signals. In recent years, a series selfadaptive time–frequency analysis method were proposed, such as empirical mode decomposition (EMD) [12], local mean decomposition (LMD) [13] and intrinsic time scale decomposition (ITD) [14]. They are completely signal-driven and adaptive non-parametric approaches which extract the intrinsic oscillating mode by data fitting or data smoothing, and were intensive applied in the fault diagnosis process of machinery with signal of nonstationary characteristics. Among them, LMD is a self-adaptive time–frequency analysis method proposed by Smith, which is proven to be better than EMD at inhibiting mode mixing, end effect, computation cost and retaining original signal information of frequency and envelope [15]. However, the end effect and mode mixing problems are two main limitations of LMD. The two limitations can be alleviated by proper parameter selection on boundary condition, envelope function estimation, and sifting stopping criterion [3]. Most published papers aiming at improving performance of LMD can be grouped into three categories corresponding to the above three parameters. Among them, envelope estimation are the main focuses during the past ten years. The construction of envelope estimation function in LMD algorithm plays an important role in the accuracy of the PF components it extracts. In the original LMD algorithm, envelope function are calculated by using the moving averaging method according to the local extreme points. While the moving averaging method may cause phase error of the functions after several iterations, and thus has an adverse effect on the accuracy of PF components [16]. Cubic spline interpolation gives an interpolating polynomial that is smoother and has smaller error than some other interpolating polynomials, and it is a widely used method in many fields, Qiao [17,18] has developed a regularized cubic B-spline collocation method to identify the impact force time history, which overcomes the deficiency of the ill-posed problem. Xue [19] has introduced a modified Hermitian cubic spline wavelets on interval to avoid the boundary problem, results indicate that this method has a higher precision and costs less time. Inspired by the empirical mode decomposition algorithm, Hu and Ren [20] has presented a novel approach to calculate the envelope estimation function using the envelopes of local maximum points and local minimum points by cubic spline interpolation. This approach can not only improve the computational efficiency, but also avoid the selection of the step size for moving averaging method. While, similar to EMD, the envelopes obtained by CSI may occur overshoot and undershoot problems due to its continuous second derivative for vibration signals with strong nonstationary characteristic. Several interpolation methods have been further proposed to improve the performance of envelope estimation function. For example, the B-spline and some higher order polynomial interpolation methods are used to replace the CSI technique [21,22]. Monotonic Piecewise Cubic Hermite Interpolation (MPCHI) is an alternative interpolation method, and its derivative is only first order continuous and can be set monotonically [23,24]. Therefore, to calculate the envelope estimation function by using MPCHI can avoid the overshoot and undershoot problem of CSI envelope for strong nonstationary signals. Li [25] and author [26] have proposed a MPCHI based LMD method, and this improves the accuracy of the PF components it extracts from signal of strong nonstationary characteristics. While, after a detailed observation we found the amplitude of strong nonstationary signal doesn’t vary significantly at any time, and it is strong impulsive sometimes, but stationary at most other time. Therefore, to acquire a better decomposition accuracy, the envelope should have distinctive characteristics for different parts of signal, and the envelope should be smooth for stationary part to ensure the stability, but be flexible for strong nonstationary part to avoid the overshoot and undershoot distortion. However, in the traditional envelope interpolation method the envelope can only be constructed using one interpolation method. CSI is a smooth interpolation method due to its continuous second derivative, while MPCHI is only first order continuous and can be set monotonically, and it is more flexible compared to CSI. Hence, to build the envelope with CSI for the stationary part and MPCHI for the strong nonstationary part of signal, this paper propose a novel compound interpolation envelope method using distinctive interpolation method for different parts of signal. In this research, we are interested in a compound interpolation envelope (CIE) method for LMD algorithm based on MPCHI and CSI and its application for fault diagnosis of reciprocating compressor oversized bearing clearance. In Section 2 we review the LMD, CSI and MPCHI method. In Section 3 we explain the compound interpolation envelope method based on MPCHI and CSI. In Section 4 we describe the numerical simulation of CIE LMD. Feature diagnosis process for the reciprocating compressor oversized bearing clearance states is presented in Section 5. Finally, conclusion is outlined in Section 6.

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

275

2. The review of LMD, CSI and MPCHI method 2.1. The cubic spline enveloped LMD LMD method decomposes a complicated signal into a set of product functions, each of which is the product of an envelope signal and a purely frequency modulated signal. In the traditional LMD algorithm, the local mean function and envelope function are calculated by using the moving averaging method according to the local extrema points, while this algorithm calculates the two functions using the upper and lower envelopes of local maximum points and local minimum points derived by CSI [22]. For any signal x(t), it can be decomposed as follows: (1) Determine all the local maximum points and the local minimum points of the original signal x(t), and then calculate the upper envelope Eu and the lower envelope El of extrema by using CSI. (2) The continuous local mean function m11(t) and the envelope function a11(t) are given by the upper envelope Eu(t) and the lower envelope El(t) as

m11 ðtÞ ¼ ðEu ðtÞ þ El ðtÞÞ=2

ð1Þ

a11 ðtÞ ¼ jEu ðtÞ  El ðtÞj=2

ð2Þ

(3) The local mean function m11(t) is subtracted from the original data x(t).

h11 ¼ xðtÞ  m11

ð3Þ

h11(t) is then divided by a11(t).

s11 ðtÞ ¼ h11 ðtÞ=a11 ðtÞ

ð4Þ

The envelope a12(t) of s11(t) can then be calculated. If a12(t) – 1, the procedure needs to be repeated for s11(t). A local mean m12(t) is calculated for s11(t), subtracted from s11(t), and the resulting function amplitude demodulated using a12(t). This iteration process continues n times until a purely frequency modulated signal s1n(t) is obtained. In practice, a variation d can be determined in advance. If 1  d  a1(n+1)(t)  1 + d and 1  s1n(t)  1, then iterative process would be stopped [15]. (4) The corresponding envelope a1(t) is given by

a1 ðtÞ ¼ a11 ðtÞa12 ðtÞ    a1n ðtÞ ¼

n Y a1q ðtÞ

ð5Þ

q¼1

Multiplying envelope signal a1(t) by the purely frequency modulated signal s1n(t), the first product function PF1 (t) of the original signal can be obtained.

PF 1 ðtÞ ¼ a1 ðtÞs1n ðtÞ

ð6Þ

(5) Subtract the first PF component PF1(t) from the original signal x(t) and we have a new signal u1(t), which becomes the new original signal and the whole of the above procedure is repeated k times until uk becomes a monotonic function.

8 u1 ðtÞ ¼ xðtÞ  PF 1 ðtÞ > > > > < u2 ðtÞ ¼ u1 ðtÞ  PF 2 ðtÞ .. > > . > > : uk ðtÞ ¼ uk1 ðtÞ  PF k ðtÞ

ð7Þ

Thus, the original signal x(t) was decomposed into k PF components and a monotonic function uk

xðtÞ ¼

k X PF k ðtÞ  uk ðtÞ

ð8Þ

p¼1

where uk is the residual term, k is the number of the PF components. 2.2. Cubic spline interpolation (CSI) Cubic spline interpolation is a special case for Spline interpolation, and it gives an interpolating polynomial that is smoother and has smaller error than some other interpolating polynomials such as Lagrange polynomial and Newton polynomial. Given an interval [a,b] composed of n ordered disjoint subintervals [xi1, xi] with: a = x0 < x1 <    < xn 1 < xn = b, and the corresponding function values yi = () (i = 0,1, . . ., n), if a function S(x) satisfies the following conditions [27]:

276

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

(1) On each subinterval [xi1, xi] (i = 0,1, . . ., n), S(x) is a cubic polynomial; (2) (xi) = yi (i = 0,1, . . ., n); (3) On the interval [a, b], if the second derivative of (x), S(x), is a continuous function, S(x) is called a cubic spline interpolation of the function y = f(x) on the interval [a, b]. The essential idea is to fit a piecewise function of the form

SðxÞ ¼

8 s1 ðxÞ > > > > < s2 ðxÞ > > > > :

x1 6 x 6 x 2 x2 6 x 6 x 3 .. .

ð9Þ

sn1 ðxÞ xn1 6 x 6 xn

where si is a third degree polynomial defined by

si ðxÞ ¼ ai ðx  xi Þ3 þ bi ðx  xi Þ2 þ ci ðx  xi Þ þ di

ð10Þ

for i = 1, 2, . . ., n  1. The first and second derivatives of these n1 equations are fundamental to this process, and they are

s0i ðxÞ ¼ 3ai ðx  xi Þ2 þ 2bi ðx  xi Þ þ ci s00i ðxÞ ¼ 6ai ðx  xi Þ þ 2bi

ð11Þ

for i = 1, 2, . . ., n1. To determine this cubic spline S(x), we need to determine ai, bi, ci, di for each i by:

Sðxi Þyi ; Sðxi  0Þ ¼ Sðxi þ 0Þ;

ði ¼ 2; . . . ; n  1Þ

S0 ðxi  0Þ ¼ S0 ðxi þ 0Þ;

ði ¼ 2; . . . ; n  1Þ

S00 ðxi  0Þ ¼ S00 ðxi þ 0Þ;

ð12Þ

ði ¼ 2; . . . ; n  1Þ

We can see that there are n + n + (n  1) + (n  1) = 4n  2 conditions, but we need to determine 4n efficients, so usually we add two boundary conditions to solve this problem. There are three types of common boundary conditions: I. First derivatives at the endpoints are given as:

S0 ðx1 Þ ¼ y01 ;

S0 ðxn Þ ¼ y0n

ð13Þ

II. Second derivatives at the endpoints are given as:

S00 ðx1 Þ ¼ y001 ;

S00 ðxn Þ ¼ y00n

ð14Þ

The special case y001 = yn‘‘ = 0 is called natural or simple boundary conditions. III. When the exact function f(x) is a periodic function with period, S(x) is a periodic function with period b-a too. Thus

S0 ðx1 þ 0Þ ¼ S0 ðxn  0Þ;

S00 ðx þ 0Þ ¼ S00 ðxn  0Þ

ð15Þ

The spline functions S(x) satisfying this type of boundary condition are called periodic splines. The Fritsch-Carlson RPN 15A data which is a monotone increasing example have been used to form the CSI envelope [28], and this envelope was given in Fig. 1. We can see that the envelope are not monotone in some intervals, and it will cause the overshoot and undershoot problem. 2.3. Monotonic Piecewise Cubic Hermite Interpolation (MPCHI) The performance of MPCHI mainly depends on its first derivatives of interpolation points. Through a suitable definition of the first derivative-value, it can preserve the monotonicity of envelopes, which can avoid the overshoot and undershoot problem of CSI envelops. The MPCHI can be defined as follow. Let (xi, Yi, di), i = 0,. . ., n be given real data, where a = x0 < x1 < . . . < xn = b is a partition of the interval [a, b], while yi and di, i = 0,. . .,n are assigned function values and derivative-values respectively. Further let hi, Di, Dyi, i = 0,. . .,n, be defined by

hi ¼ xiþ1  xi ;

Dyi ¼ yiþ1  yi ; Di ¼ Dyi =hi :

ð16Þ

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

277

0 Fig. 1. The CSI envelope for simulation example.

We suppose that the data are monotonic, i.e. Dyi P 08i or Dyi 6 08i and



di ¼ diþ1 ¼ 0 for Di ¼ 0 ; sgnðdi Þ ¼ sgnðdiþ1 Þ ¼ sgnðDi Þ for Di –0

ð17Þ

where we assume the convention sgn(0) = sgn(Di ). The piecewise cubic function SðxÞ 2 C 1 ½a; b interpolating the given values, i.e. Sðxi Þ ¼ yi and S0 ðxi Þ ¼ di , i = 0,. . .,n, is defined as

SðxÞ  Si ðxÞ ¼

ðdi þ diþ1  2Di Þ 2 hi

ðx  xi Þ3 þ

ð2di  diþ1 þ 3Di Þ ðx  xi Þ2 þ di ðx þ xi Þ þ yi hi

ð18Þ

for x e [xi,xi+1]. It is known that Eq. (17) is a necessary but not sufficient condition for the monotonicity of Si ðxÞ and a complete characterization of the monotonicity of Si ðxÞ can be given in terms of the ratios Di and Diþ1 [29,30].

dðiÞ ¼

3Dimin Dimax

Dimax þ 2Dimin

ð19Þ

where Dimax ¼ maxðDi ; Diþ1 Þ, Dimin ¼ minðDi ; Diþ1 Þ At the boundaries we use the second-order uncentered parabolic method, and the parabolic algorithms is

dð0Þ ¼

ð2h0 þ h1 ÞDi  h0 D1 h0 þ h1

ð20Þ

for the left end point,

dðnÞ ¼

ð2hn1 þ hn2 ÞDn1  hn1 Dn2 hn1 þ hn2

ð21Þ

for the right end point. The Fritsch-Carlson RPN 15A data in Fig. 1 was also used to form the MPCHI envelope, and this envelope was given in Fig. 2. We can see that the envelope are monotone in all intervals, and it can avoid the overshoot and undershoot problem in CSI envelope. 3. Compound interpolation envelope (CIE) LMD algorithm 3.1. Comparison between CSI and MPCHI envelops As mentioned in the introduction section, the envelope should be smooth for the stationary part of signal to ensure the stability, but be flexible for strong nonstationary part to avoid the overshoot and undershoot distortion. CSI is a smooth interpolation method due to its continuous second derivative, while MPCHI is more flexible compared with CSI for its only first order continuous derivative and monotonicity. To illustrate this phenomenon, we employ a strong nonstationary simulation

278

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

Fig. 2. The MPCHI envelope for simulation example.

signal to compare the envelopes of CSI and MPCHI, and this signal which contains a periodic impact signal and a sinusoidal signal was defined as follow:

8 > < > :

yðtÞ ¼ y1 ðtÞ þ y2 ðtÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  n2 t

ð22Þ

y1 ðtÞ ¼ yn enxn t sin xn

y2 ðtÞ ¼ ym sin xm t

where yn , n, xn are the amplitude, damping coefficient and natural frequency of periodic impact signal respectively, and yn was defined as 5, n as 0.02, xn as 10,000. ym , xm are the amplitude and natural frequency of sinusoidal signal respectively, and ym was defined as 5, xm as 100. The time t was defined as 0:0.000001:0.025 s, and the impact signal repeats every 0.005 s. The simulation signal was shown in Fig. 3, and we can see that there are five impact waves in the signal. At the beginning of each impact wave, the amplitude changes significantly, and it is the strong nonstationary part. While the amplitude of following part decreases smoothly, and it is the smooth part. To compare the behaviors of two envelope construction method, the upper and lower envelops of this signal were created by CSI and MPCHI respectively, and they were both shown in Fig. 4 with the original signal. Some local parts were zoomed in on the top of this figure to show the details, and we find that the envelope of CSI occurs overshoot and undershoot problem in the nonstationary part due to its smooth characteristics, while the envelop of MPCHI can avoid this problem for its flexibility. For a further analysis on the envelope accuracy for the nonstationary and smooth part, the two local mean functions were gained from the upper and lower envelopes created by the two methods according to Eq. (1) respectively, then they were compared with the theoretic local mean function which is the y2 ðtÞ in Eq. (22), and the absolute value errors for both methods were shown in Fig. 5. From Fig. 5 we can observe that the two errors both reach the maximum at the beginning point of

8 6

smooth part

nonstationary part

Amplitude

4 2 0 -2 -4 -6 -8

0

0.005

0.01

0.015

Time(s) Fig. 3. The strong nonstationary simulation signal.

0.02

0.025

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

279

Fig. 4. The overshoot and undershoot problem in the strong nonstationary simulation signal.

impact wave, while on the area near the beginning point in the local zoom, that is the nonstationary part, the error of CSI are much bigger than MPCHI. On the contrary, in the smooth part, that is the area away from the two beginning point of impact wave from each side, the errors of MPCHI are much higher than CSI. Based on above analysis, we can conclude that the MPCHI performances well than CSI in nonstationary part, and CSI is more suitable than MPCHI in the smooth part. Therefore, it is necessary to propose an envelope construction approach to combine these two methods. 3.2. Two key issues of compound interpolation envelope method In this section we will propose a compound interpolation envelope method by using CSI and MPCHI, while there are two key issues to be further discussed. One of the issue is to judge the nonstationary characteristic of the local signal, that is, to determine which part should use CSI or MPCHI to build the envelopes. Another issue is how to connect the end point between the envelopes built by two different methods for the smoothness of this envelope. When a local signal is strong nonstationary, envelop built by CSI will occurs overshoot and undershoot problem, and we use MPCHI envelop to replace it. To judge which part should use CSI or MPCHI for building the envelopes, we must give a standard to evaluate the nonstationary characteristic of the local signal. After a detailed observation, we find that when the slope Di of adjacent interpolation points varies greatly, for example Di is more than two times Diþ1 or less than an half of Diþ1 , the CSI built envelop will occurs overshoot or undershoot problem, and some examples are shown in Fig. 6. Therefore, we employ the ratio between adjacent slopes Di and Diþ1 as the standard to evaluate the nonstationary characteristic of the local signal []:

 Zi ¼

jDiþ1 =Di j for jDiþ1 j P jDi j jDi =Diþ1 j for jDiþ1 j < jDi j

ð23Þ

where Z i is the nonstationary coefficient. If Z i is higher than the given threshold value, the local envelope will be built by MPCHI, otherwise it will be created by CSI. After determined the method to build the envelope of local signal, the next issue is to connect the end point between adjacent envelopes and keep this new envelope smooth. CSI is a second derivative continuous function, while MPCHI is only first derivative continuous. As discussed in Section 2.2, CSI has three methods to determine the boundary conditions of envelope end point, one of them is to set the first derivative manually in condition I, while the first derivative of MPCHI in envelope end point has to calculate according to Eqs. (20) and (21) in Section 2.3. In order to keep the envelope smooth, the first derivative must be the same at end point. Therefore, we employ the value of MPCHI as the common first derivative at this connecting end point, and assign it to CSI according to Eq. (13) as the boundary conditions during the calculate process. 3.3. The steps and flowchart of compound interpolation envelope algorithm The main steps for the compound interpolation envelope construction method is: first finding out the nonstationary points and building the local MPCHI envelopes, then obtaining the envelope for the other interval by CSI method, finally combining the local envelopes together. For any signal x(t), it can be created as follows:

280

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

Fig. 5. The absolute value errors of CSI and MPCHI envelope for the simulation signal.

Fig. 6. The nonstationary coefficient of local signal.

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

281

x(t) Determine the local extrema series ni

Calculate slopes and nonstationary coefficient

Find out the position Lj of nonstationary points

Define [Lj-1, Lj+1] as the local MPCHI interval

Y Is local MPCHI interval empty?

N Create the local MPCHI envelopes

Find out the local CSI interval [Nk-1, Nk+1]

Y Is local CSI interval empty?

N Create the local CSI envelopes Connect the local MPCHI envelops and local CSI envelops

Fig. 7. The flowchart of compound interpolation envelope.

(1) Determine the local extrema series ni (i = 1,. . ., w) of the original signal x(t) (the extrema series ni stands for either the local maximum points nmaxi or the local minimum points nmini); (2) For one of the extrema series ni (the local maximum points nmaxi or the local minimum points nmini), first calculate its slopes Di (i = 1,. . ., w  1) according to Eq. (19), then get its nonstationary coefficient Z i (i = 2,. . ., w  2) based on Eq. (23);

282

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

(3) Find out those nonstationary points of which the nonstationary coefficient Z i is higher than the given threshold value, then record their position as Lj (j = 1,. . ., l); (4) Define [Lj  1, Lj + 1] as the local MPCHI interval (Because the MPCHI needs more than one interpolation point, here we define the point n(Lj) plus the previous point n(Lj  1) and following point n(Lj + 1) as the local MPCHI interpolation points); if the position of two adjacent nonstationary points position Li and Lj+1 are less than five extrema points(counts by i), then connect the two local MPCHI interval together as [Lj  1, Lj+1 + 1]; if the distance from the position of nonstationary points Li to the left or right end of extrema series (i = 1 or w) is less than three extrema points(counts by i), then the local MPCHI interval begins or ends from the left or right end of extrema series, that is [1, Lj + 1] or [Lj  1, w]; (5) If the local MPCHI interval [Lj  1, Lj + 1] is empty, create the envelope by using CSI for the whole extrema series ni (i = 1,. . ., w), else create envelopes for all the intervals [Lj  1, Lj + 1] (j = 1,. . ., l) by using MPCHI, then record the first derivative of left and right end points for each local envelope; (6) Find out the local CSI interval [Nupk, Nlowk] (k = 1,. . ., p) by remove the local MPCHI interval [Lj  1, Lj + 1] from the whole position series of ni (i = 1,. . ., w), if the local CSI interval [Nupk, Nlowk] is empty, create the envelope by using MCHPI for the whole extrema series ni (i = 1,. . ., w); (7) Create envelopes for the local CSI interval [Nupk, Nlowk] (k = 1,. . ., p) by using the first derivative of end points from MPCHI envelope as the boundary conditions; (8) Connect the local MPCHI envelops and local CSI envelops in turn to form the whole envelop for the extrema series ni (i = 1,. . ., w) of the original signal x(t); (9) Repeat the steps (2)–(8) for envelop of the other extrema series. Based on the above steps, a flowchart of compound envelope interpolation method is shown in Fig. 7. We use this method instead of CSI to build the upper and lower envelopes in step(1) of LMD algorithm given in Section 2.1, and follow the rest steps of it, the CIE LMD algorithm is proposed[20].

4. Simulated signal analysis In this section, we compare the proposed CIE LMD method with the original LMD, CSI LMD and MPCHI LMD to validate its effectiveness by using the above mentioned strong nonstationary numerical example in Eq. (22). The simulated signal x(t) was decomposed by the four LMD methods in turn, and the results were shown from Figs. 8–11, respectively. The threshold value of nonstationary coefficient in CIE LMD is 2. Obviously, four methods all succeed in extracting two PF components from signal x(t), and they are very similar in general. These two components correspond to the impulsive signal and the sinusoidal signal x1(t) in results of all methods. Furthermore, we use the Mean Squared Error (MSE), iterations of PF component and total CPU time to assess their performance, and the details are given in Table 1. The simulated signal contains a periodic impact and a sinusoidal signal, and the obtained PF components correspond to these two parts. Therefore, the Mean Squared Error (MSE) between PF component and its corresponding signal is a best indicator to reflect the performance of proposed method. We can know that the proposed method has the minimal MSE for both PF components than other three methods in Table 1, while the CSI LMD method has the maximal MSE because the overshoot and undershoot problem occurred at the strong nonstationary point.

5

PF1

0 -5 0

0.005

0.01

0.015

0.02

0.025

0.005

0.01

0.015

0.02

0.025

0.005

0.01

0.015

0.02

0.025

5

PF2

0 -5 0 5

u2

0 -5 0

Time [s] Fig. 8. The PF components of original LMD for simulation signal.

283

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

5

PF1

0 -5 0

0.005

0.01

0.015

0.02

0.025

0.005

0.01

0.015

0.02

0.025

0.005

0.01

0.015

0.02

0.025

5

PF2

0 -5 0 5

u2

0 -5 0

Time [s] Fig. 9. The PF components of CSI LMD for simulation signal.

Table 1 Performance of decomposition results for simulated signal with different LMD methods. Method

PF1

Original LMD CSI LMD MPCHI LMD CIE LMD

PF2

Time/s

MSE

Iterations

MSE

Iterations

0.0061 0.0093 0.0053 0.0047

5 4 3 2

0.0224 0.0496 0.0198 0.0121

5 8 4 3

0.732 0.665 0.465 0.608

5

PF1

0 -5 0

0.005

0.01

0.015

0.02

0.02

0.005

0.01

0.015

0.02

0.02

0.005

0.01

0.015

0.02

0.02

5

PF2

0 -5 0 5

u2

0 -5 0

Time [s] Fig. 10. The PF components of MPCHI LMD for simulation signal.

In general, a relatively small number of sifting iterations may result in some benefits in a single sifting process. This is because of not only the thing of efficiency, but also the thing that it may reduce error accumulation effects in a sifting process and lead to a high quality of modulation information mining [3]. The proposed method also outperforms than other methods due to its local self-adaptive envelopes for a more accurate local mean function in the sifting process. Furthermore, total CPU time reflects an overall efficiency, and a smaller value indicates a more efficient process. The whole decomposition time of the CIE LMD is 0.628 s, which is slower than 0.465 s of MPCHI LMD, and similar with 0.665 s of CSI LMD. The MPCHI is a high efficiency method than CSI in CSI LMD and moving average (MA) in original LMD, and MPCHI LMD has a less sifting iterations than CSI LMD and LMD, so it is the efficient LMD method. Though the CIE LMD

284

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

5

PF1

0 -5 0

0.005

0.01

0.015

0.02

0.02

0.005

0.01

0.015

0.02

0.02

0.005

0.01

0.015

0.02

0.02

5

PF2

0 -5 0 5

u2

0 -5 0

Time [s] Fig. 11. The PF components of CIE LMD for simulation signal.

has the least sifting iterations, it employs a nonstationary points searching process and builds the envelope using MPCHI and CSI separately, so it is slower than MPCHI LMD. The CSI in CSI LMD and MA in original LMD use more sifting iterations, hence they are slower than the other two methods. Based on above analysis, though CIE LMD is not the most efficient method, it out performs in the MSE and sifting iterations of PF components for a better accuracy. 5. Application for fault diagnosis of reciprocating compressor 5.1. Oversized bearing clearance state experiment Reciprocating compressor is one of the machineries widely used to compress gas in chemical industry, and their operational safety is a significant subject of the intensive investigation. A two-stage 2D12 type reciprocating compressor with a 500 kW shaft power, 240 mm piston stroke and 496 rpm motor rotation speed is studied in this paper, and its structural diagram was shown in Fig. 12(a). As shown in Fig. 12(b), the double acting pistons in this two-stage compressor were driven by a crank-connecting rod mechanism, and the journal bearings which used to assemble the connecting rod to crankshaft or to crosshead pin often occur oversized clearance fault due to tolerances and defects in manufacturing process or wearing after a certain working period, and this may cause the vibratory running condition which can stop the compressor. Thus, an effective measure to diagnose the clearance state of bearings is necessary in reducing the maintenance costs of the equipment. Vibration based measurement and analysis technique is an effective measure to monitor and diagnose fault state of machinery due to the abundant operational state information embodied in vibration signal. An experiment on clearance states of bearing between the first stage crankshaft pin and connecting rod was conducted in a real working condition. Three bearing clearance states, a 0.1 mm normal clearance state, two 0.25 mm and 0.35 mm oversized clearance states simulated with worn bearing bushing, were tested respectively by an ICP acceleration sensor placed on the top of crosshead guide surface to collect vibration signal. The position of sensor was marked with a red square in Fig. 12(a) and (c), and the sample frequency is 5 KHz. 5.2. Analysis of vibration signals in fault states The tested vibration acceleration signals in three clearance states were shown from Figs. 13–15 respectively. The signal of normal state in Fig. 13 is chaotic, while shocks appear in Figs. 14 and 15, and with the growth of clearance size the amplitude of shocks increases significantly. A proper clearance size can keep the shaft or journal contact with the metal sleeve or shell continuously. So, as we can see from Fig. 16, there are no obvious shocks in the envelope frequency spectrum for normal bearing clearance state. After a long work period, the clearance increases gradually due to tolerances and defects in manufacturing process or wearing. The oversized clearance destroyed the continuous oil film in bearing, then journal and sleeve impact directly and excite strong contact forces. They usually occurs two impacts for they change the movement direction twice in one crank revolution, and this will cause a double rotation frequency peak in the envelope frequency spectrum of vibration acceleration shown in Figs. 17 and 18. However, due to the complex structure and the nonlinear factors such as clearance, stiffness and the unbalanced and time-varying forces of components, the vibration signals of reciprocating compressor may present many noisy frequencies except the periodical clearance induced frequency, which disturbs the accurate

285

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

Motor

Cylinder of first stage

Crosshead

second stage 3

Sensor 2

Cylinder of

Crosshead

Crankcase

Right

Front

(a) The structural draw of reciprocating compressor Sliding bearings

(b) The two-stage reciprocating compressor transmission mechanism. Sensor

Sensor Sensor

( ) The test bench of reciprocating compressor Fig. 12. A two-stage double-acting reciprocating compressor of type 2D12.

recognition of fault information. So a signal pretreat method is necessary for the accurate diagnosis of bearing clearance state from the collected vibration acceleration signal.

5.3. CIE LMD decomposition of the vibration signal According to the theory of structural dynamics, the vibration acceleration of a structure can be considered as the superposition of different natural frequencies modulated by the excitation forces. While LMD method can decompose a multicomponent signal into a series of PF components, where each component is the product of an envelope signal and a purely

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

2

Acceleration [m/s ]

286

Times [s]

2

Acceleration [m/s ]

Fig. 13. The vibration acceleration in normal bearing clearance state.

Times [s]

2

Acceleration [m/s ]

Fig. 14. The vibration acceleration in 0.25 bearing clearance state.

Times [s] Fig. 15. The vibration acceleration in 0.35 mm bearing clearance state.

2

Acceleration [m/s ]

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

Frequency [Hz]

Acceleration [m/s 2]

Fig. 16. The envelope frequency spectrum of vibration acceleration in normal bearing clearance state.

Frequency [Hz]

2

Acceleration [m/s ]

Fig. 17. The envelope frequency spectrum of vibration acceleration in 0.25 mm bearing clearance state.

Frequency [Hz] Fig. 18. The envelope frequency spectrum of vibration acceleration in 0.35 mm bearing clearance state.

287

288

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

frequency-modulated signal. Therefore, LMD is an appropriate method to extract fault information from the vibration signal of reciprocating compressor. However, due to the strong nonstationary characteristics of reciprocating compressor vibration signal, the common envelope construction methods may lose its accuracy in sifting process of PF components. This novel LMD method proposed a compound envelope approach according to the local characteristics of signal to improve the accuracy of envelope function and local mean function. The CIE LMD was employed to decompose the vibration signals of two oversized clearance states in Figs. 14 and 15 using 2 as threshold value of nonstationary coefficient. To verify the performance of the proposed LMD, we also employ the MPCHI LMD, CSI LMD and original LMD to decompose these two vibration signals.

5.4. Results and discussion Usually the first several PF components contain the main fault information, so we give the first three PF components of four LMD methods for each oversized clearance state. The decomposition results of 0.25 mm oversized clearance state were illustrated from Figs. 19–22. We can know that the first PF component of four LMD are more dominated with the strong periodical shocks, and the PF components of CSI LMD have a mutation due to the overshoot and undershoot problem of CSI envelope. We also give the envelope frequency spectrum of the first PF component for four LMD methods from Figs. 23–26. We

PF1

PF2

PF3

Time [s] Fig. 19. Decomposition results of CIE LMD for vibration signal in 0.25 mm bearing clearance state.

PF1

PF2

PF3

Time [s] Fig. 20. Decomposition results of MPCHI LMD for vibration signal in 0.25 mm bearing clearance state.

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

PF1

PF2

PF3

Time [s] Fig. 21. Decomposition results of CSI LMD for vibration signal in 0.25 mm bearing clearance state.

PF1

PF2

PF3

Time [s]

Acceleration [m/s 2]

Fig. 22. Decomposition results of Original LMD for vibration signal in 0.25 mm bearing clearance state.

Frequency [Hz] Fig. 23. The envelope frequency spectrum of the first PF component with CIE LMD in 0.35 mm bearing clearance state.

289

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

Acceleration [m/s 2]

290

Frequency [Hz] Fig. 24. The envelope frequency spectrum of the first PF component with MPCHI LMD in 0.35 mm bearing clearance state.

0.8 0.7

Acceleration [m/s2]

0.6 0.5 0.4 0.3 0.2 0.1 0 0

20

40

60

80

100

120

140

160

180

200

Frequency [Hz] Fig. 25. The envelope frequency spectrum of the first PF component with CSI LMD in 0.25 mm bearing clearance state.

can see a clearly peak in double rotation frequency followed by multi rotation frequency, and the noise frequency are further suppressed compared with the envelope frequency spectrum in Fig. 17, and it is because some undesired noise are removed in the process of LMD decomposition. Furthermore, in these four envelope frequency spectrums, CIE LMD not only has the highest peak value in double rotation frequency, but also suppress the noise frequency effectively. This is because the novel envelope construction method can self-adaptive with the local characteristic of signal which indeed improves the performance of PF components it extracted. The decomposition results of 0.35 mm oversized clearance state were illustrated from Figs. 27–30, and the first PF component of four LMD are dominated with the stronger periodical shocks compared with those of 0.25 mm oversized clearance state. The PF components of CSI LMD have a more serious mutation due to the stronger impulsive local characteristics of signal. The envelope frequency spectrum of the first PF component for four LMD are shown from Figs. 31–34. Even though the envelope frequency spectrums in Fig. 18 has a peak in double rotation frequency, but it is followed by many strong noise frequency. The envelope frequency spectrum of PF component of four LMD method have a more significant peak, and CIE LMD also has the highest peak value in double rotation frequency and suppress the noise frequency effectively, which can improve recognition accuracy of fault state. (see Figs. 31–34) For further analysis, a detailed comparison between the four methods on iterations of PF component, Average Index of Orthogonality (IOav) and Index of Energy Conservation (IEC) were given in Tables 2 and 3 for two oversized clearance fault states.

Acceleration [m/s2]

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

Frequency [Hz] Fig. 26. The envelope frequency spectrum of the first PF component with Original LMD in 0.25 mm bearing clearance state.

PF1

PF2

PF3

Time [s] Fig. 27. Decomposition results of CIE LMD for vibration signal in 0.35 mm bearing clearance state.

PF1

PF2

PF3

Time [s] Fig. 28. Decomposition results of MPCHI LMD for vibration signal in 0.35 mm bearing clearance state.

291

292

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

PF1

PF2

PF3

Time [s] Fig. 29. Decomposition results of CSI LMD for vibration signal in 0.35 mm bearing clearance state.

PF1

PF2

PF3

Time [s]

2

Acceleration [m/s ]

Fig. 30. Decomposition results of Original LMD for vibration signal in 0.35 mm bearing clearance state.

Frequency [Hz] Fig. 31. The envelope frequency spectrum of the first PF component with CIE LMD in 0.35 mm bearing clearance state.

2

Acceleration [m/s ]

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

Frequency [Hz]

2

Acceleration [m/s ]

Fig. 32. The envelope frequency spectrum of the first PF component with MPCHI LMD in 0.35 mm bearing clearance state.

Frequency [Hz]

2

Acceleration [m/s ]

Fig. 33. The envelope frequency spectrum of the first PF component with CSI LMD in 0.35 mm bearing clearance state.

Frequency [Hz] Fig. 34. The envelope frequency spectrum of the first PF component with Original LMD in 0.35 mm bearing clearance state.

293

294

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

Table 2 Results performance for 0.25 mm oversized clearance state with different LMD methods. Method

Iterations

Original LMD CSI LMD MPCHI LMD CIE LMD

PF1

PF2

PF3

10 14 8 6

9 15 7 5

11 17 6 4

IOav e

m

0.267 0.313 0.178 0.121

0.932 0.915 0.978 0.992

IOav e

IEC

0.244 0.313 0.189 0.154

0.924 0.908 0.950 0.994

Table 3 Results performance for 0.35 mm oversized clearance state with different LMD methods. Method

Iterations

Original LMD CSI LMD MPCHI LMD CIE LMD

PF1

PF2

PF3

9 14 7 5

9 12 6 5

10 15 6 4

As mentioned in Section 3, a relatively small number of sifting iterations may reduce error accumulation effects and lead to a high quality of modulation information mining. All the three PF components of the CIE LMD have the least sifting iterations due to its local self-adaptive envelopes for two oversized clearance fault states. We can see that the second and third PF component of CSI LMD in Figs. 21 and 29 both have a mutation, and they have the largest number of sifting iterations. This is due to the overshoot and undershoot problems of such strongly nonstationary vibrational signals enveloped by cubic spline curves, thus reducing the accuracy of the local mean function and the envelope function. In theory, every two PF components should be orthogonal, so the Index of Orthogonality (IO) between the two PF components should be equal to zero [31]. While errors actually exist in the local mean function, the orthogonality between two PF components is relative. So a small value of IO indicates a high orthogonality. To compare all PF components of a LMD method, the average of IO (IOav) is defined as follows:

IOav ¼ Av erage

jhpfj ; pfk ij kpfj k2 kpfk k2

j–k;

ð24Þ

where x(t) is the original signal, pfj(t) is PF component. From another perspective, if the PF components of the LMD decomposition are orthogonal, the energy of this signal before and after decomposition should be conservative, and IEC equals one. Usually the more accurate the PF components are, the closer to one the IEC is [32]. The IEC are defined as

PT Pn t¼0

IEC ¼ PT

j¼1 jpfj ðtÞj

t¼0 jxðtÞ

2

 un ðtÞj2

ð25Þ

where x(t) is the original signal, pfj(t) is PF component and un(t) is the residual component. Based on the results as summarized in Tables 2 and 3, we can see that the CIE LMD indeed give a better results compare with other three methods for two oversized clearance fault states. The proposed LMD performs especially well on the energy conservation test, which indicates that the total energy deviation amongst the components is the smallest. 6. Conclusions Based on the local characteristics of strong nonstationary vibration signals, a CIE LMD was proposed through a novel envelope construction method, and this LMD method improved the fault recognition accuracy significantly for the oversized bearing clearance fault of reciprocating compressor. Through a numerical simulation example analysis, we found that the MPCHI envelope performances well than CSI envelope in nonstationary part, and CSI envelope is more suitable than MPCHI envelope in the smooth part. A nonstationary coefficient was defined to evaluate the local nonstationary characteristics of vibration signals, then a compound interpolation envelope construction method which use MPCHI for nonstationary part and CSI for smooth part was proposed. A numerical simulation example was conducted to verify the performance of CIE LMD, results indicated that CIE LMD is superior on MSE and sifting iterations of PF components, though it is not the most efficient method. The CIE LMD was employed to analyze the vibration signals for two oversized bearing clearance fault of reciprocating compressor, and it outperforms other three LMD methods on the sifting iterations, IOav and IEC. Furthermore, the envelope frequency spectrum of PF component presents a more significant peak of fault frequency.

Z. Haiyang et al. / Mechanical Systems and Signal Processing 110 (2018) 273–295

295

Acknowledgements This work was supported by China Scholarship Council (201608230308), General Financial Grant from The China Postdoctoral Science Foundation (2015M581423), Natural Science Foundation of Heilongjiang Province in China (E2015037, E2016009), and Natural Science Foundation in China (51505079). References [1] X. Chen, S. Wang, B. Qiao, et al, Basic research on machinery fault diagnostics: past, present, and future trends, Front. Mech. Eng. (2017) 1–28. [2] Z. Feng, M. Liang, F. Chu, Recent advances in time–frequency analysis methods for machinery fault diagnosis: a review with application examples, Mech. Syst. Signal Process. 38 (2013) 165–205. [3] Z.L. Liu, M.J. Zuo, Y.Q. Jin, et al, Improved local mean decomposition for modulation information mining and its application to machinery fault diagnosis, J. Sound Vib. 397 (9) (2017) 266–281. [4] M. Elhaj, F. Gub, A.D. Ballb, Numerical simulation and experimental study of a two-stage reciprocating compressor for condition monitoring, Mech. Syst. Signal Prcess. 22 (2008) 374–389. [5] Y. Tang, Q. Liu, J. Jing, et al, A framework for identification of maintenance significant items in reliability centered maintenance, Energy 118 (2017) 1295–1303. [6] H.Y. Zhao, M.Q. Xu, J.D. Wang, Y.B. Li, A Parameters optimization method for planar joint clearance model and its application for dynamics simulation of reciprocating compressor, J. Sound Vib. 334 (5) (2015) 416–433. [7] J.B. Ali, N. Fnaiech, L. Saidi, et al, Application of empirical mode decomposition and artificial neural network for automatic bearing fault diagnosis based on vibration signals, Appl. Acoust. 89 (2015) 16–27. [8] B.S. Yang, W.W. Hwang, D.J. Kim, A.C. Tan, Condition classification of small reciprocating compressor for refrigerators using artificial neural networks and support vector machines, Mech. Syst. Signal Process. 19 (2) (2005) 371–390. [9] Y. Yang, H.Y. Pan, L. Ma, J.S. Cheng, A fault diagnosis approach for roller bearing based on improved intrinsic timescale decomposition de-noising and kriging-variable predictive model-based class discriminate, J. Vib. Control (2014) 1–16. [10] Z.P. Feng, M. Liang, Complex signal analysis for wind turbine planetary gearbox fault diagnosis via iterative atomic decomposition thresholding, J. Sound Vib. 333 (2014) 5196–5211. [11] H.Y. Zhao, M.Q. Xu, J.D. Wang, Local mean decomposition based on rational hermite interpolation and its application for fault diagnosis of reciprocating compressor, J. Mech. Eng-En. 51 (1) (2015) 83–89. [12] N.E. Huang, Z. Shen, S.R. Long, et al, The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis, Proc. Roy. Soc. London, Ser. A. 454 (1998) 903–995. [13] J.S. Smith, The local mean decomposition and its application to EEG perception data, J. Roy. Soc. Interface 2 (2005) 443–454. [14] M. Frei, I. Osorio, Intrinsic time-scale decomposition: time–frequency–energy analysis and real-time filtering of non-stationary signals, Proc. Roy Soc. 463 (2007) 321–342. [15] J.S. Cheng, Y. Yang, A rotating machinery fault diagnosis method based on local mean decomposition, Digit. Signal Proces. 22 (2012) 356–366. [16] K. Zhang, J.S. Cheng, Y. Yang, The local mean decomposition method based on rational spline and its application, J. Vib. Eng. 24 (2009) 97–103. [17] B. Qiao, X. Chen, X. Xue, et al, The application of cubic B-spline collocation method in impact force identification, Mech. Syst. Signal Process. 64 (2015) 413–427. [18] B. Qiao, X. Zhang, X. Luo, et al, A force identification method using cubic B-spline scaling functions, J. Sound Vib. 337 (2015) 28–44. [19] X. Xue, X. Zhang, B. Li, et al, Modified Hermitian cubic spline wavelet on interval finite element for wave propagation and load identification, Finite Elem. Anal. Des. 91 (2014) 48–58. [20] J.S. Hu, S.X. Yang, D.Q. Ren, Spline-based local mean decomposition method for vibration signal, J. Data Acqu. Proc. 24 (2009) 78–83. [21] M.D. Wang, L.B. Zhang, L. Wei, Local mean decomposition method based on B-spline interpolation, J. Vib. Shock 29 (1) (2010) 73–78. [22] L.F. Deng, R.Z. Zhao, An improved spline-local mean decomposition and its application to vibration analysis of rotating machinery with rub-impact fault, J. Vibroeng. 16 (1) (2014) 414–433. [23] J.L. Merrien, P. Sablonnière, Rational splines for Hermite interpolation with shape constraints, Comput. Aided Geometric Des. 30 (2013) 296–309. [24] R.J. Cripps, M.Z. Hussain, C1 monotone cubic Hermite interpolant, Appl. Math. Lett. 25 (2012) 1161–1165. [25] L.H. Li, L.H. Guo, T. Chen, Iris recognition based on PCHIP-LMD, Opt. Precis. Eng. 21 (01) (2013) 197–206. [26] H.Y. Zhao, J.D. Wang, H. Han, Y.Q. Gao, A feature extraction method based on HLMD and MFE for bearing clearance fault of reciprocating compressor, Measurement 89 (2016) 34–43. [27] S.A. Dyer, X. He, Cubic-spline interpolation: part 2, Instrum. Meas. Mag. IEEE 4 (2001) 34–36. [28] James M. Hyman, Accurate monotonicity preserving cubic interpolation, SIAM J. Sci. Stat. Comput. 4 (4) (1983) 645–654. [29] M.G. Gasparo, R. Morandi, Pieeewise cubic monotone interpolation with assigned slopes, Computing 46 (1991) 355–365. [30] F.N. Fritsch, R.E. Carlson, Monotone piecewise cubic interpolation, Soc. Industr. Appl. Math. 17 (2) (1980) 238–247. [31] D.R. Huang, A sufficient condition of monotone cubic splines, Comput. Math. 2 (1982) 214–215. [32] Q.H. Chen, N.S. Huang, S. Riemenschneider, A B-spline approach for empirical mode decompositions, Adv. Comput. Math. 24 (2006) 171–195.