Generalized adaptive mode decomposition for nonstationary signal analysis of rotating machinery: Principle and applications

Mechanical Systems and Signal Processing 136 (2020) 106530

Contents lists available at ScienceDirect

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

Generalized adaptive mode decomposition for nonstationary signal analysis of rotating machinery: Principle and applications Zhipeng Feng a, Xinnan Yu a, Dong Zhang a,⇑, Ming Liang b a b

School of Mechanical Engineering, University of Science and Technology Beijing, Beijing 100083, China Department of Mechanical Engineering, University of Ottawa, Ottawa, ON K1N 6N5, Canada

a r t i c l e

i n f o

Article history: Received 9 June 2019 Received in revised form 13 November 2019 Accepted 16 November 2019

Keywords: Adaptive mode decomposition Mono-component Generalized demodulation Time-frequency representation Fault diagnosis

a b s t r a c t Rotating machinery signals are usually intricate and often nonstationary. The analysis of such signals is further complicated because of the spectral overlaps between constituent frequency components under time-varying running conditions. Most conventional signal processing methods are unable to effectively extract the embedded meaningful information. Adaptive mode decomposition (AMD) methods are flexible to describe complex signals, free from the limitations inherent in conventional basis expansion and thus more suited to unveiling the underlying physical nature. However, as those methods work in a way similar to filter banks, i.e. via rectangular partition of time-frequency plane, they are effective only for signals with constituent frequency components that are almost constant and parallel to each other. As such, it would be inappropriate to apply such methods when signal constituent frequency components vary considerably over time and particularly overlap in frequency domain. To address this issue, this paper proposes a general framework by exploiting the unique capability of generalized demodulation to transform an arbitrary time-varying instantaneous frequency into a constant frequency. In doing so, many state-of-the-art AMD methods (including empirical mode composition, local mean decomposition, intrinsic time-scale decomposition, local characteristic scale decomposition, Hilbert vibration decomposition, empirical wavelet transform, variational mode decomposition, and adaptive local iterative filtering) can be extended to analyze highly nonstationary complex signals with spectral overlaps. The principle and advantage of generalized adaptive mode decomposition (GAMD) are illustrated through numerical simulation. The GAMD’s performance in monitoring different rotating machinery components (including a planetary gearbox, an aircraft engine rolling bearing and a hydraulic turbine rotor) is demonstrated and validated under time-varying conditions. Ó 2019 Elsevier Ltd. All rights reserved.

Abbreviations: TFA, time-frequency analysis; TFR, time-frequency representation; STFT, short-time Fourier transform; CWT, continuous wavelet transform; WVD, Wigner-Ville distribution; IGD, iterative generalized demodulation; AMD, adaptive mode decomposition; GAMD, generalized adaptive mode decomposition; EMD, empirical mode decomposition; GEMD, generalized empirical mode decomposition; LMD, local mean decomposition; GLMD, generalized local mean decomposition; ITD, intrinsic time scale decomposition; GITD, generalized intrinsic time scale decomposition; LCD, local characteristic scale decomposition; GLCD, generalized local characteristic scale decomposition; HVD, Hilbert vibration decomposition; GHVD, generalized Hilbert vibration decomposition; EWT, empirical wavelet transform; GEWT, generalized empirical wavelet transform; VMD, variational mode decomposition; GVMD, generalized variational mode decomposition; ALIF, adaptive local iterative filtering; GALIF, generalized adaptive local iterative filtering. ⇑ Corresponding author. E-mail address: [email protected] (D. Zhang). https://doi.org/10.1016/j.ymssp.2019.106530 0888-3270/Ó 2019 Elsevier Ltd. All rights reserved.

2

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

1. Introduction Rotating machinery condition monitoring and fault diagnosis significantly rely on detecting the presence of characteristic frequencies in dynamic signals and/or monitoring changes in their magnitudes. In practice, the dynamic signal of rotating machinery is often dominated by multiple frequency components generated from various mechanical parts, such as rotors, shafts, gears and bearings. These frequencies are essentially the integer or fractional harmonics of the rotating frequency. Under nonstationary running conditions, particularly during time-varying speed processes, they change over time, following the rotating speed profile. In this situation, dynamic signals feature strong nonstationarity. To accurately reveal the frequency contents and time-varying pattern of nonstationary signals, it is highly desirable to develop effective signal analysis methods [1]. Time-frequency analysis (TFA) methods depict signal frequency composition and time-varying pattern in a joint timefrequency domain, and are widely used to analyze nonstationary signals. However, conventional TFA methods suffer from some drawbacks. For example, linear time-frequency representations (TFRs) including short-time Fourier transform (STFT) and continuous wavelet transform (CWT) are essentially an integral transform with a predefined basis or window function, subject to Heisenberg uncertainty principle and thus limited time-frequency resolution. Bilinear TFRs all originate from the Wigner-Ville distribution (WVD) and involve a double integral and quadratic term of a signal, therefore suffering from both outer interferences due to cross term between multi-components and inner interferences existent with nonlinearly timevarying frequency mono-components [2–5]. To improve time-frequency readability, reassignment methods [6] reallocate diluted energy to the gravity center of time-frequency distributions, and synchrosqueezing transform [7] compresses diluted energy in frequency/scale domain aided by instantaneous frequency estimation. However, these post-processing methods may fail when signal components are closely adjacent to each other on a time-frequency plane. Adaptive mode decomposition (AMD) methods are truly data-driven and posterior. They require neither any a priori basis function to match signal characteristics, nor a double integral and quadratic term of a signal. With an AMD method, a signal is first decomposed into mono-components according to its local characteristics. Then, the instantaneous frequency and amplitude envelope of each mono-component are calculated based on the Hilbert transform. Next, TFR of each monocomponent is obtained by allocating its energy along the corresponding instantaneous frequency. Finally, a high quality TFR can be constructed as a linear superposition of all constituent components’ TFR. Such a TFR not only has high timefrequency resolution but also is free from both inner and outer interferences [8–17]. This eliminates the two main hurdles present in the conventional TFAs and makes AMD methods flexibly fit to transient features and capable of highlighting local characteristics of signals, hence enhancing time-frequency readability. As such, AMD has attracted increasingly more attention from researchers in related fields. In 1998, Huang et al. [9] proposed empirical mode decomposition (EMD) to decompose a nonstationary signal into constituent mono-components for instantaneous frequency calculation, and then construct TFR for nonstationary signal analysis. The EMD separates intrinsic mode functions (IMFs) in a frequency order from high to low, and exhibits a dyadic filter bank property when applied to white noise [10]. However, it has some shortcomings such as the lack of mathematical formulation, susceptibility to mode mixing under singularities, instability under noise interferences, and over/under fitting due to cubic spline interpolation [8–10]. To address these issues, researchers have proposed some new methods. For instance, the local mean decomposition (LMD) by Smith [11] in 2006, the intrinsic time scale decomposition (ITD) by Frei and Osorio [12] in 2007, the local characteristic scale decomposition (LCD) by Cheng et al. [13] in 2013, following the same mono-component sifting framework of EMD to solve the mode mixing and over/under fitting problems; the Hilbert vibration decomposition (HVD) by Feldman [14] in 2006, the empirical wavelet transform (EWT) by Gilles [15] in 2013, the variational mode decomposition (VMD) by Dragomiretskiy and Zosso [16] in 2014, and the adaptive local iterative filtering (ALIF) by Cicone et al [17] in 2016, exploit the amplitude modulation and frequency modulation (AM-FM) property of intrinsic mode function (IMFs) or the filter bank nature of EMD to obtain a rigorous mathematical formulation and better robustness to noise. These algorithms improve AMD. However, all these AMD methods share the same filter bank property, and therefore are subject to a narrow band assumption where the frequencies of signal components are required to be near constant. This property makes them only suitable to analyze signals without spectral overlap, for example those composed of almost constant frequencies. However, when signal frequencies vary excessively over time and overlap with each other in frequency domain, these methods no longer perform satisfactorily. This is because of their filter bank partition nature in the time-frequency plane. Conventional filtration essentially divides a time-frequency plane into rectangles along time axis within a specific pass band, and the pass region of a conventional filter cannot adapt to and track the time trajectory of signal frequency components. For signals with spectral overlaps, filtration will split the instantaneous frequency trajectories on the time-frequency plane, and lead to interferences from neighboring components or energy leakage and information loss, thus unable to obtain true constituent monocomponents. Take Fig. 1(a) for example (see Eq. (12) in Section 2.3 for detailed explanation on this numerical simulation signal), no filtration can effectively separate any frequency component intact, neither the wide band pass filter within the dashed lines, nor the narrow one within the dot dashed lines. Suppose the middle frequency component is to be separated. Within a wide pass band, for instance the one defined by the dashed lines, although the middle component can be separated, part of either the upper or lower component will also be separated together with the middle one, resulting in interferences.

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

3

Fig. 1. Filtration rectangularly divides time-frequency plane.

Within a narrow pass band, for example the one defined by dot dashed lines, the filtration will split the target middle component, causing information loss. The middle component can be separated without any interferences from either upper or lower component but only the beginning part is separated, while the ending part is lost. This problem will cause deficiency of AMD methods. Generalized demodulation can transform any time-varying frequency component into a constant frequency component by multiplying a demodulation function (a conjugate complex exponential function of each component’s instantaneous phase). This makes the target component separable from other neighboring ones through filtration and finally recoverable by multiplying an inverse demodulation function (a complex exponential function of each component’s instantaneous phase) [18]. Feng and his colleagues [19–22] proposed an iterative generalized demodulation (IGD) method, and successfully extracted vibration features of a hydraulic turbine and gear fault features of planetary gearboxes in joint time-frequency domain. Shi et al. [23] further extended it to rolling bearing fault diagnosis. It should be noted, however, that although the IGD can describe the time-frequency feature of major frequency components, it relies on conventional filtering to separate the components of interest. As a result, the filtered components cannot be guaranteed to be mono-components because the passband inevitably contains background noise. Consequently, the estimated instantaneous frequency based on the Hilbert transform is error-prone. The above discussions lead us to address an important open issue: the deficiency of AMD caused by spectral overlaps. To this end, this paper will exploit the unique capability of generalized demodulation in transforming arbitrary time-varying instantaneous frequency into a constant one. The instantaneous frequency estimation errors inherent in the IGD will be eliminated by applying AMD to separate mono-component subject to background noise. The above goal can be achieved as follows. Firstly, for a target frequency component, design a suitable generalized demodulation phase vector, referring to the instantaneous frequency ridge of target component. Then, perform generalized demodulation on the signal by multiplying with the designed generalized demodulation phase vector. This will transform the target frequency component into an almost constant frequency component, without spectral overlap with other components, thus being separable from other frequency components through time-frequency partition like filters. Next, apply AMD methods to the generalized demodulated signal, and separate the target mono-component. Finally, recover the true frequency component through inverse generalized demodulation, i.e. multiply the target mono-component with the inverse generalized demodulation phase vector. Once the true constituent components are separated, accurate instantaneous frequency calculation and high quality TFR are ready to be generated. Compared with the angular resampling in order tracking to convert time-varying frequencies (in time domain) into constant ones (in angular domain), the generalized demodulation is free of the error due to interpolation and/or decimation in angular resampling, and retains accurately the time information essential to time-frequency analysis. This generalization overcomes the time-frequency smearing problem inherent with AMD due to spectral overlaps, and avoids noise influences on filtered component inherent with IGD. This proposed method enables AMD to effectively analyze arbitrary complex signals even with spectral overlaps, and more applicable in rotating machinery fault diagnosis applications. Hereafter, this paper is organized as follows. Section 2 introduces the principle of generalized adaptive mode decomposition (GAMD) methods. Section 3 illustrates their performance via numerical simulation. Section 4 validates their capability through vibration signal analysis of a planetary gearbox lab experiment, an aircraft rolling bearing ground test and a hydraulic turbine on-site measurement. Section 5 draws conclusions.

2. Generalized adaptive mode decomposition In this section, the principle of generalized demodulation is introduced first, then those of AMD methods are revisited briefly, and thereafter GAMD is presented.

4

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

2.1. Generalized demodulation Olhede and Walden [18] proposed generalized demodulation for time-frequency analysis of nonstationary signals. Through generalized demodulation, a time-varying frequency component can be mapped into any desired shape on a time-frequency plane, thus becoming separable from other neighboring components through filtration. Generalized demodulation is essentially an extension of the modulation or frequency shift property of the Fourier transform from time-invariant case to time-variant case. Multiplying a signal xðtÞ by a demodulation function expf  j½tðtÞ  x0 tg, and applying the Fourier transform, yields

Z

XðxÞ ¼

1

1

Z

xðtÞexpf  j½tðtÞ  x0 tgexp½  jxtdt ¼

1

1

xðtÞexpf  j½ðx  x0 Þt þ tðtÞg dt

ð1Þ

where t and x stand for the time and angular frequency respectively, tðtÞ is a demodulation phase function, x0 the base 2

frequency, and j ¼ 1. This is defined as generalized Fourier transform. Correspondingly, the inverse generalized Fourier transform is defined as

Z

xðtÞ ¼

1

1

Z

XðxÞexpf j½tðtÞ  x0 tgexpðjxtÞdx ¼ expf j½tðtÞ  x0 tg

1

1

XðxÞexpðjxtÞdx

ð2Þ

where expf j½tðtÞ  x0 tg is an inverse demodulation function. According to Eq. (1), if xðtÞ ¼ exp½jtðtÞ, then XðxÞ ¼ dðx  x0 Þ, where dðÞ is the Dirac delta function. This means that a signal with a time-varying instantaneous frequency t_ ðtÞ is mapped to a constant frequency point x ¼ x0 after applying the generalized Fourier transform. On the time-frequency plane, a constant frequency component concentrates on a line parallel to the time axis but vertical to the frequency axis, being separable from others via filtering. If the instantaneous frequency trajectory of a signal component on the time-frequency plane is specified by a function xðtÞ ¼ t_ ðtÞ, and the signal component is to be mapped into a passband around the base frequency x0 on the time-frequency plane, then a demodulation phase R function is needed to approximate the time-varying phase tðtÞ  x0 t ¼ xðtÞdt. After generalized demodulation, the target component will not overlap with other components in frequency domain, thus becoming separable through AMD because of their mono-component decomposition capability. For a separated constant frequency component expðjx0 tÞ, multiplying it by an inverse demodulation function expf j½tðtÞ  x0 tg recovers the original time-varying frequency component [18]

xðtÞ ¼

R1

1

dðx  x0 Þexpf j½tðtÞ  x0 tgexpðjxtÞdx R1 dðx  x0 Þexp½jðx  x0 Þtdx 1

¼ exp½jtðtÞ

ð3Þ

¼ exp½jtðtÞ This is inverse generalized demodulation. Both the demodulation and inverse demodulation functions are unit complex exponential, therefore the generalized demodulation and inverse generalized demodulation do not affect the amplitude envelopes of signals. As such, the result generated from the above generalized demodulation, separation, and inverse generalized demodulation retains the true instantaneous frequency and amplitude envelope of the original frequency component. 2.2. Adaptive mode decomposition Instantaneous frequency is fundamental and is the key to describing the time-varying pattern of each frequency component and revealing the underlying physical nature of nonstationary complex signals. It is meaningful for single frequency component only, and therefore is defined based on mono-component. However, most real world signals are of multicomponents. In order to decompose complex multi-component signals into their constituent mono-components, various AMD methods including EMD, LMD, ITD, LCD, HVD, EWT, VMD and ALIF have been proposed. These AMD methods share almost the same filter bank nature, and some follow a similar algorithmic framework. For example, in EMD, LMD, ITD, LCD and ALIF, local minima and maxima in a signal are detected, and an instantaneous mean (a ‘‘low-pass” centerline) is estimated based on these local extrema in local characteristic time scale. Then, the instantaneous mean is removed, and the high frequency component is separated as a mono-component. Such procedure is recursively applied on the remaining ‘‘low-pass” centerline, until all the mono-components are separated. For a real signal xðtÞ, their algorithm is as follows [8-13,17]. (1) Initialize parameters: Set iteration index i ¼ 1, residual signal r 0 ðtÞ ¼ xðtÞ. (2) Extract the ith mono-component ci ðtÞ: (2.1) Let k ¼ 0, and the proto-mode function hik ðtÞ ¼ ri1 ðtÞ. (2.2) Find the local minima and the local maxima of hik ðtÞ.

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

5

(2.3) Calculate the instantaneous mean mik ðtÞ of hik ðtÞ according to the local extrema. (2.4) Update the proto-mode function hik ðtÞ ¼ hik ðtÞ  mik ðtÞ. (2.5) If hik ðtÞ satisfies the stop criterion for mono-component sifting, then set the ith mono-component ci ðtÞ ¼ hik ðtÞ. Otherwise, let k ¼ k þ 1, return to step (2.2). (3) Update the residual signal r i ðtÞ ¼ ri1 ðtÞ  ci ðtÞ. (4) If ri ðtÞ satisfies the stop criterion for adaptive mode decomposition, terminate the decomposition. Otherwise, let i ¼ i þ 1, return to step (2). For the inner loop, i.e., steps (2.1)–(2.6), a Cauchy-type criterion is usually used. When the standard deviation of two consecutive proto-mode functions reaches a predefined threshold (usually set between 0.2 and 0.3) stop the mono-component sifting iteration, i.e., steps (2.1)–(2.6). For the outer loop, namely steps (2)–(4), when the residual signal r i ðtÞ becomes a trend, i.e., it has no more than one local extremum, stop the AMD iteration, i.e., steps (2)–(4). However, the instantaneous mean estimation approaches of EMD, LMD, ITD, LCD and ALIF are different. In EMD, the upper and lower envelopes of residual signal are derived by interpolating to local minima and maxima using cubic splines respectively, and the instantaneous mean is set to the local average of upper and lower envelopes [9]. In LMD, the instantaneous mean is estimated as the smoothed mean of consecutive local minima and maxima [11]. The instantaneous mean in ITD is obtained through piecewise linear fitting to consecutive local minima and maxima [12]. For LCD, the mirror of local minima (maxima) between consecutive local maxima (minima) is calculated via linear fitting first, then the local average of these local extrema and their mirrors is deduced, and finally the instantaneous mean is obtained by cubic spline interpolating to these local averages [13]. In ALIF, the instantaneous mean is derived by low pass filtering the residual signal [17]. EWT and VMD exploit the AM-FM feature and spectral characteristics of mono-components [15,16]. EWT is inspired by the compact Fourier support of AM-FM mono-components. It segments the signal Fourier spectrum first and then separates mono-components within each Fourier support through wavelet filtering, where wavelets are defined empirically by the Fourier supports rather than prescribed by a fixed dilation factor. On each Fourier support ½xn1 ; xn , the empirical scaling function and empirical wavelet are defined as band pass filters

8 jxj  xn  sn 1; > > > < ^ ðxÞ ¼ cosfp b½ 1 ðjxj  xn þ sn Þg; xn  sn  jxj  xn þ sn / n 2sn 2 > > > : 0; other

^ ðxÞ ¼ w n

8 1; xn þ sn  jxj  xnþ1  snþ1 > > > > > > p 1 > < cosf 2 b½2snþ1 ðjxj  xnþ1 þ snþ1 Þg; xnþ1  snþ1  jxj  xnþ1 þ snþ1 > > sinfp2 b½21sn ðjxj  xn þ sn Þg; > > > > > : 0;

xn  sn  jxj  xn þ sn

ð4aÞ

ð4bÞ

other

respectively, where sn ¼ cxn , when c < min½ðxnþ1  xn Þ=ðxnþ1 þ xn Þ, and bðxÞ ¼ x4 ð35  84x þ 70x2  20x3 Þ is an arbitrary n

C k ð½0; 1Þ function. Given the empirical scaling and wavelet functions, the empirical mode on each Fourier support can be constructed through empirical wavelet transform and inverse transform [15]. VMD shares the same nature as EWT, and is also motivated by the narrow-band properties of AM-FM IMF definition. It assumes each AM-FM mono-component concentrate around its center frequency xi , which is set as the gravity center of ith AM-FM mono-component’s spectrum. Accordingly, it separates AM-FM mono-components by minimizing their bandwidth

min

fci ðtÞg;fxi g

2   X X @  1   ci ðtÞ ¼ xðtÞ @t ½dðtÞ þ j pt   ci ðtÞ expðjxi tÞ s:t:; 2

i

ð5Þ

i

where dðÞ is the Dirac delta function, and the symbol  denotes the convolution operator [16]. HVD separates mono-components following a descending amplitude order. It obtains the reference instantaneous frequency x1 ðtÞ of the largest amplitude component through low pass filtering the signal instantaneous frequency, and then estimates the amplitude envelope via synchronous detection and low pass filtering [14]. Through synchronous detection, the in-phase and quadrature parts are obtained as

R R Ai ðtÞcos½ xi ðtÞdt þ ui cos½ x1 ðtÞdt i R ¼ 12 Ai ðtÞðcosui þ cosf ½xi ðtÞ þ x1 ðtÞdt þ ui gÞ

xi¼r ðtÞ ¼

P

ð6aÞ

6

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

^xi¼r ðtÞ ¼

 Z  1 ½xi ðtÞ þ x1 ðtÞdt þ ui Ai ðtÞ sinui þ sin 2

ð6bÞ

respectively, where Ai ðtÞ, xi ðtÞ and ui are the instantaneous amplitude, frequency and phase of the ith component [14]. Low pass filtering the in-phase and the quadrature parts, removes the fast oscillating components and yields the instantaneous amplitude and instantaneous phase [14]

(

hxi¼r ðtÞi ¼

1 A ðtÞcos 2 i

0;

ui ; xi ¼ x1 ^ ; hxi¼r ðtÞi ¼ other

(

1 A ðtÞsin 2 i

ui ; xi ¼ x1

0;

ð7aÞ

other

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h^xi¼r ðtÞi 2 Ai ðtÞ ¼ 2 hxi¼r ðtÞi2 þ h^xi¼r ðtÞi ; ui ¼ arctan hxi¼r ðtÞi

ð7bÞ

R Then, the largest amplitude component can be reconstructed as x1 ðtÞ ¼ a1 ðtÞcos½ x1 ðtÞdt. Repeating iteratively the synchronous detection and low pass filtering will separate mono-components one by one. In summary, all these AMD methods share filter bank properties (for example, EMD, LMD, ITD, LCD and ALIF), involve signal filtration (for instance ALIF and HVD), or rely on Fourier spectral analysis (for example, EWT and VMD). Therefore, they require signals to be of (near) constant frequencies or free of spectral overlaps. 2.3. Generalized adaptive mode decomposition Generalized demodulation with a proper demodulation function can effectively avert overlap of a target signal constituent component with other components in the frequency domain, thus making it easily separable from others. Most rotating machinery signals are mainly composed of rotating frequency harmonics. The higher the harmonic order, the faster the corresponding instantaneous frequency will change. Most AMD methods (except EWT and VMD) separate monocomponents in descending order of frequency, and hence the first few separated ones are more accurate because data fitting errors tend to accumulate in those separated afterwards. As such, we separate constituent mono-components following such an order. For a given signal xðtÞ, the proposed GAMD procedure is summarized as follows, and shown as Fig. 2. (1) Initialization: Set the iteration number i ¼ 1, and the residual signal r i1 ðtÞ ¼ xðtÞ. (2) Generalized demodulation: (2.1) Determine signal components first by visual inspection of conventional TFRs (for example STFT spectrogram or CWT scalogram, and reassigned TFR for the sake of better time-frequency readability) or by adaptive surrogate tests [22], detect time-frequency ridges respectively through maxima tracking f ðtÞ ¼ arg max jTFRðt; f Þj, and then estimate their f

instantaneous frequencies t_ ðtÞ ¼ f ðtÞ by polynomial fitting to their respective time-frequency ridge. Usually, such analyses on raw signal can roughly estimate the number of components to be separated and the highest harmonic order. (2.2) Design the generalized demodulation phase function, referencing to the highest harmonic component identified in the residual signal. For the target highest frequency component t_ i ðtÞ, construct its corresponding demodulation phase function ti ðtÞ  x0 t by integrating its instantaneous frequency with respect to time. (Note: The base frequency x0 is usually set to a higher frequency, e.g., 0.4 times of the sampling frequency in this paper. After generalized demodulation, all the frequency components will be frequency shifted by x0 . The target time-varying frequency will be transformed into a constant frequency, and it will not overlap with other components. x0 acts as a signal carrier frequency, and a higher carrier frequency guarantees accurate estimation of amplitude envelope and instantaneous frequency.) (2.3) Perform generalized demodulation on the residual signal with the designed demodulation phase function ri1 ðtÞexpfj½ti ðtÞ  x0 tg. For the signal in Fig. 1(a), after a generalized demodulation with the highest frequency component as reference, its true TFR is shown in Fig. 1(b). The highest frequency component is transformed into a constant or near constant base frequency x0 without spectral overlap with the middle and lower ones, thus facilitating component separation through AMD. (3) Adaptive mode decomposition: Apply AMD to the generalized demodulated residual signal, and separate the target highest frequency component ci ðtÞ. Among the obtained mono-components, select the first one with an instantaneous frequency fluctuating around the base frequency x0 as the true one ci ðtÞ for further analysis. (Note: As explained earlier, since most AMD methods (except EWT and VMD) separate mono-components in a descending order of frequency, the first few separated ones are more accurate, while the later ones suffer from data fitting errors accumulated). (4) Inverse generalized demodulation: After separating a target frequency component, recover the original frequency component by multiplying the separated target component with an inverse generalized demodulation phase vector xi ðtÞ ¼ ci ðtÞexpfj½ti ðtÞ  x0 tg.

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

7

Fig. 2. GAMD analysis procedure.

(5) Residual signal update: Subtracting the recovered frequency component from the previous residual, yields the new residual signal ri ðtÞ ¼ r i1 ðtÞ  xi ðtÞ. (6) Stop condition check: This can be done by checking the number of components separated or the residual signal energy. If the number of separated components i does not reach a predefined value, then let the iteration number i ¼ i þ 1, and repeat steps (2) through (5). Otherwise, terminate iteration. This stop condition is suited for the case where all the frequency components are estimated through conventional TFA of the raw signal. Alternatively, calculate the ratio R R of residual signal energy to raw signal energy c ¼ r2i ðtÞdt= x2 ðtÞdt. If the ratio is larger than a predefined threshold, for example 0.0001, then let the iteration number i ¼ i þ 1, and repeat steps (2) through (5). Otherwise, terminate iteration. This stop condition is suited for the case where only the time-frequency ridge of one single or a few dominant components are detected through TFA of the raw signal. In this situation, after the detected dominant components are separated and subtracted, others will become dominant in the updated residual signal and therefore detectable. Applying the separation procedure on the residual signal iteratively, the multiple components will be separated one by one. Given all the separated constituent components, a high quality TFR of raw signal can be generated as follows. For each separated mono-component xi ðtÞ, calculate its amplitude envelope

Ai ðtÞ ¼ jxi ðtÞ þ jHT½xi ðtÞj ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2i ðtÞ þ HT2 ½xi ðtÞ

ð8Þ

and its instantaneous frequency

Table 1 AM-FM parameters. Component i

ai

ci

1 2 3

0.008 0.01 0.012

20 30 40

8

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

xi ðtÞ ¼

d HT½xi ðtÞ arctan dt xi ðtÞ

ð9Þ

via the Hilbert transform based analytic signal approach (Note: In coding implementation, the derivative can be calculated by a finite difference method, such as symmetric and asymmetric difference). Then construct its Hilbert spectrum in joint time-frequency domain

TFRi ðt; xÞ ¼ A2i ðtÞd½x  xi ðtÞ

ð10Þ

Superposing the Hilbert spectrum of all mono-components yields the TFR of raw signal

TFRðt; xÞ ¼

X

TFRi ðt; xÞ

ð11Þ

i

This proposed framework generalizes AMD methods to highly time-varying signals with spectral overlaps between constituent frequency components. Moreover, it addresses the issue inherent in IGD due to noise interferences with the components separated through conventional filters. As the TFR via Hilbert transform based analytic signal approach emphasizes local characteristics of signals, it naturally provides fine time-frequency resolution. Additionally, it does not

Fig. 3. Simulated signal and its TFRs based on conventional methods and IGD.

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

9

Fig. 3 (continued)

require double integral of signals and is hence free of both outer and inner interferences. These advantages make the proposed GAMD method well suited to the analysis of nonstationary complex multi-component signals. In the aspects of mono-component sifting and instantaneous parameter estimation, GAMD shares almost the same properties of AMD. In applications, users can refer to the pros and cons of various AMD methods summarized in [7], to select a proper GAMD method. For simplicity, hereafter we name the (iterative) GAMD versions of EMD, LMD, ITD, LCD, HVD, EWT, VMD, and ALIF as GEMD, GLMD, GITD, GLCD, GHVD, GEWT, GVMD, and GALIF respectively. 3. Numerical simulation To illustrate GAMD methods and their advantages over conventional TFA and AMD methods, an AM-FM signal is numerically simulated and analyzed. The synthetic signal consists of three AM-FM components (i.e. a quadratic chirp multiplied by a quadratic function), and a white Gaussian noise nðtÞ at signal-to-noise ratio of 0 dB

xðtÞ ¼

3 X

2

ai f ðtÞcosf2p½ci

Z f ðtÞdtg þ nðtÞ

ð12Þ

i¼1

where f ðtÞ ¼  12t2  6t þ 20, t ¼ 0; 0:0005;    ; 0:8 s, i.e. the sampling frequency is 2000 Hz, and other parameter values are listed in Table 1. This synthetic signal simulates the AM-FM characteristics of rotating machinery vibration signals during time-varying speed processes and the three constituent frequency components overlap in frequency domain. Fig. 3 shows the signal and its TFA results via conventional TFRs and IGD. As shown in Fig. 3(c) through (i), the STFT, CWT and smoothed pseudo WVD suffer from low time-frequency resolution. Their reassigned versions and synchrosqueezing wavelet transform can improve time-frequency resolution, but in the meantime generate some unwanted ripples around true time-frequency ridges and speckles due to noise interferences (see Fig. 3(f) through (i)). These shortcomings make the above conventional TFA methods unsuitable to analyze highly complex nonstationary signals. Because the three constituent frequency components overlap in frequency domain, no conventional filtration can separate any component and keep it intact. As shown in Fig. 1(a), before generalized demodulation, in a wide pass band within the dashed lines, integrity of the target middle component is kept, but interferences from upper and lower components are

10

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

involved; in a narrow pass band within the dot dashed lines, only the ending part is separated, while the beginning part is lost, leading to information loss. Since all the AMD methods share a filter bank property, the spectral overlap in frequency domain makes every original AMD method generate a blurred TFR full of meaningless speckles, and fail to reveal the true time-frequency structure, as shown in the left column in Fig. 4 (i.e., Fig. 4(a1) through (h1)).

Fig. 4. TFRs by AMD (*1, left column) and GAMD (*2, right column).

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

Fig. 4 (continued)

11

12

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530 Table 2 Comparison of computational time spent by AMD and GAMD (s). AMD

Time

GAMD

Time

EMD LMD ITD LCD HVD EWT VMD ALIF

0.286223 0.282275 0.293441 0.286973 1.360714 0.269055 0.467156 0.526033

GEMD GLMD GITD GLCD GHVD GEWT GVMD GALIF

0.302868 0.292986 0.300998 0.299544 3.559485 0.314175 1.386789 1.567320

Fig. 5. Planetary gearbox test rig: (1) motor, (2) tachometer, (3) fixed-shaft gearbox, (4) planetary gearbox stage 1, (5) planetary gearbox stage 2, (6) accelerometers, and (7) magnetic powder brake.

Table 3 Gearbox configuration. Number of gear teeth

Input Intermediate Output

Fixed-shaft gearbox

Number of gear teeth

Stage 1

Stage 2

32 80 –

– 40 72

Ring Planet Sun

Note: The number in parentheses indicates the number of planet gears.

Fig. 6. Stage 1 sun gear wear.

Planetary gearbox Stage 1

Stage 2

100 40(4) 20

100 36(4) 28

13

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530 Table 4 Planetary gearbox characteristic frequencies. Frequencies

Meshing

Sun gear rotating

Planet gear rotating

Planet carrier rotating

Sun gear fault

Planet gear fault

Ring gear fault

Notation

f m ðtÞ

f s ðtÞ

f p ðtÞ

f r ðtÞ

100 27 f d ðtÞ 175 216 f d ðtÞ

ðrÞ f p ðtÞ 1 18 f d ðtÞ 7 486 f d ðtÞ

f c ðtÞ

Stage 1

ðrÞ f s ðtÞ 2 9 f d ðtÞ 1 27 f d ðtÞ

1 27 f d ðtÞ 7 864 f d ðtÞ

20 27 f d ðtÞ 175 1512 f d ðtÞ

5 54 f d ðtÞ 175 7776 f d ðtÞ

4 27 f d ðtÞ 7 216 f d ðtÞ

Stage 2

After proper generalized demodulation, each target frequency component is transformed into a constant frequency, thus is separable through filtration, as shown in Fig. 1(b) for the highest frequency component case. The instantaneous frequency is estimated by quadratic polynomial fitting to each time-frequency ridge in the CWT scalogram as f i ðtÞ ¼ ci ð12t2  6t þ 20Þ. By subtracting a base frequency f 0 ¼ 800 Hz from each estimated instantaneous frequency, then integrating over time and multiplying by 2p, yields the forward generalized demodulation phase function ti ðtÞ ¼ 2p½ci ð4t3  3t2 þ 20tÞ  800t. Correspondingly, the inverse generalized demodulation phase function 

ti ðtÞ ¼ ti ðtÞ ¼ 2p½ci ð4t3  3t2 þ 20tÞ  800t. Using these parameters, the IGD based TFR is obtained, as shown in Fig. 3 (j). The time-frequency readability is improved. However, IGD relies on filtration to separate target components, background

Fig. 7. Baseline signals and GAMD TFRs.

14

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

Fig. 7 (continued)

noise also exists in the filter passband, and causes errors in instantaneous frequency estimation via the Hilbert transform. This can be seen in Fig. 3(j), in the ending part from 0.6 s to 0.8 s, the noise interferences cause all the three timefrequency ridges to ripple around and deviate from the true instantaneous frequencies. The proposed GAMD framework can address the issues due to spectral overlaps and noise in passband. Using the same forward and inverse generalized demodulation phases as in IGD, the GAMD based TFRs are obtained, as displayed in the right column in Fig. 4, (a2) through (h2). In contrast to the TFRs via AMD, after enhancement by GAMD, the TFRs are free of outer and inner interferences inherent in bilinear time-frequency distributions, and small ripples in instantaneous frequencies in IGD due to noise interferences. More importantly, they also lead to fine time-frequency resolution. Therefore, the three true time-varying frequencies are clearly revealed, and their time-varying amplitudes are also accurately presented. The timefrequency ridges coincide with true instantaneous frequencies. In this paper, the parameters for each AMD method are set to the default values recommended by its respective inventors [8–17], and are utilized for the corresponding GAMD version. In EMD, LMD, ITD, LCD, ALIF, and their GAMD versions, the Cauchy-type stop criterion is utilized in the mono-component sifting process, where the threshold is set to 0.25. For all AMD and GAMD methods, the number of components to be separated is set to the number of time-frequency ridges manually detected in the CWT scalogram of raw signal, where Morlet wavelet with a half length equal to the square root of signal

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

15

length is used. For example, for the numerical simulated signal analysis in this section, the number of components to be separated is set to 3, according to the time-frequency ridges detected in the CWT scalogram, Fig. 3(d). The computational time of GAMD depends on that of AMD methods and the number of components to be separated. Suppose the signal is composed of K components, and the computational time for AMD of raw signal is O. Then the AMD will be repeated for K times. Each time one component is separated. The total computational time will be KO to the extreme (Note: As the components are separated one by one, the new residual signal becomes simpler, and therefore the computational time for AMD of the new residual signal will become shorter). This is confirmed by the numerical simulated signal analysis. See Table 2 for the comparison of computational time spent by AMD and GAMD, where the Matlab program runs on a laptop with an Intel Core i7 2.7 GHz Dual CPU and 8 GB RAM.

4. Applications to rotating machinery signal analysis In this section, we demonstrate that the proposed method is not only effective in monitoring one type of rotating machinery component but also versatile in detecting faults for different rotating machinery components under time-varying conditions. This is carried out by analyzing data from three different rotating machinery components from different sites: a) gear condition signals collected from a planetary gearbox, b) rolling bearing signals from an aircraft engine ground test, and c) hydraulic turbine rotor signals on-site measured during a transient process. The details are presented in the following subsections.

Fig. 8. Sun gear wear signals and GAMD TFRs.

16

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

Fig. 8 (continued)

4.1. Planetary gearbox fault diagnosis Planetary gearbox vibration signals feature multiple modulations due to: (1) the effects of both amplitude modulation (AM) and frequency modulation (FM) of gear faults on meshing vibrations, and (2) an additional AM effect on meshing vibrations caused by the time-varying vibration propagation paths from rotating sun-planet and planet-ring gear pair meshing locations to vibration sensors fixed on gearbox casing. These modulation effects lead to intricate sidebands around gear meshing frequency harmonics. Under time-varying speeds, these sidebands change over time as they are proportional to rotating speed. Gear fault diagnosis relies on detecting the presence of these sidebands and monitoring their changes in magnitude [24–26]. Therefore, reliable time-frequency analysis is the key to effective identification of these intricate timevarying sidebands.

17

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

Fig. 9. Civil aircraft engine gearbox schematic.

Table 5 Roller bearing fault frequencies with reference to shaft L5 rotation.

Fault frequency

Outer race

Inner race

Roller

Cage

7.759

10.241

3.556

0.431

4.1.1. Experiment settings Fig. 5 shows the planetary gearbox test rig at the University of Ottawa. It consists of an induction motor, a two-stage fixed-shaft gearbox, a two-stage planetary gearbox and a magnetic powder brake. All the gearboxes are used as reducers. For the planetary gearbox, the ring gear is fixed, while the sun gear and planet carrier are rotatable. The input shaft is connected to the sun gear, and the planet carrier to the output shaft. Table 3 lists the gearbox configuration parameters. Two tests are conducted: (1) baseline case when all gears are healthy, and (2) sun gear wear case when the stage 1 sun gear has wear distributed on every gear tooth while all the other gears are healthy. Fig. 6 shows the gear fault picture. During the tests, the induction motor speed f d ðtÞ changes sinusoidally, and vibration signals are collected from accelerometers fixed on the planetary gearbox casing. According to the gearbox configuration, the characteristic frequencies are calculated, as listed in Table 4.

4.1.2. Signal analysis Figs. 7 and 8 display the baseline and sun gear wear signal analysis results obtained using GAMD respectively. In the baseline case, Fig. 7(d) through (k), the time-varying sidebands appear around the gear meshing frequency, and they mainly correspond to the sun gear rotating frequency or the drive motor rotating frequency. For example, the meshing frequency minus ðrÞ

ðrÞ

three times the sun gear rotating frequency f m ðtÞ  3f s ðtÞ, plus the sun gear rotating frequency f m ðtÞ þ f s ðtÞ, and plus six and seven time the drive motor rotating frequency f m ðtÞ þ ð6; 7Þf d ðtÞ. Additionally, the drive motor rotating frequency harmonics, such as ð2; 6; 7Þf d ðtÞ are also observed. In contrast, in the sun gear wear case, Fig. 8(d) through (k), in addition to the drive motor rotating frequency harmonics ðrÞ

ð2; 6; 7Þf d ðtÞ and the meshing frequency plus the sun gear rotating frequency f m ðtÞ þ f s ðtÞ, a few new time-varying sidebands appear around the gear meshing frequency. These new sidebands are related to the stage 1 sun gear fault frequency, for instance, the meshing frequency minus the sun gear fault frequency f m ðtÞ  f s ðtÞ, plus twice the sun gear fault frequency f m ðtÞ þ 2f s ðtÞ. These detected features imply a stage 1 sun gear fault, which is consistent with the actual stage 1 sun gear condition.

4.2. Rolling bearing fault diagnosis Rolling bearing faults lead to amplitude modulation on resonance vibration at a modulating frequency equal to the bearing fault frequency. Rolling bearing fault frequency is proportional to shaft rotating frequency. Under time-varying speed conditions, the AM modulating frequency changes over time, following the time-varying profile of shaft rotating frequency.

18

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

Since the AM part (amplitude envelope) of vibration signals contains the key modulating frequency information, we identify the time-varying bearing fault frequency through time-frequency analysis of amplitude envelope. 4.2.1. Test settings Fig. 9 shows the schematic of a civil aircraft engine gearbox. The engine has two main shafts and an accessory gearbox with pumps, filters, alternators and starter. The accessory gearbox is linked to the high pressure shaft (labelled HP in the figure) by a radial drive shaft (labelled RDS) and a horizontal drive shaft (labelled HDS). In the gearbox, the roller bearing supporting shaft L5 has spalling on outer race, in a sector of 32° along all the functional line of the race, and with a depth of approximately 0.1 mm. Table 5 lists the bearing fault frequencies with reference to the shaft L5 rotational frequency. A tachometer is mounted on shaft L4 and has a resolution of 44 pulses per revolution. An accelerometer is located on the gearbox flange close to shaft L5. During a slow acceleration from idle to full power in a ground test, the vibration and tachometer signals are collected at a sampling frequency of 50 kHz. The signal amplitude is normalized to a unitary maximum amplitude [27]. 4.2.2. Signal analysis To demonstrate the capability of proposed method, we analyze the amplitude envelope of accelerometer signal from 150 s to 162 s approximately when the engine speed varies substantially and the constituent frequency components overlap in frequency domain. Fig. 10 shows the TFRs derived by GAMD methods. They clearly reveal the frequency composition and their time-varying pattern within the frequency band of [7.3, 7.7] kHz, even though the frequencies are highly time-varying. All the major time-frequency ridges in this frequency band link to four times the outer race fault frequency 4f o , and its sum or difference combinations with the integer or fractional harmonics of cage rotating frequency kf c , such as 4f o  f c , 4f o  2f c , 4f o þ 3f c and 4f o  ð1=2Þf c . They form sidebands around four times the outer race fault frequency 4f o with a spacing equal to

Fig. 10. Rolling bearing fault signals and GAMD TFRs.

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

19

Fig. 10 (continued)

the cage rotating frequency f c or its half ð1=2Þf c . Such features indicate the outer race fault and the cage unbalance, which are consistent with the actual test settings. 4.3. Hydraulic turbine vibration signal analysis Hydraulic turbine rotor vibration signals are mainly composed of rotor rotating frequency and its harmonics. These frequency components are key indicators for hydraulic turbine health monitoring. Usually, hydraulic turbines run under low speeds but hydraulic force is volatile. Particularly, hydraulic turbines often switch between different states. This makes on-site measured vibration signals nonstationary [28,29]. Effective identification of frequency contents and their timevarying pattern is an important topic for hydraulic turbine condition monitoring. 4.3.1. On-site measurement settings Fig. 11 shows the schematic of a vertical axial-flow hydraulic turbine in a hydroelectric power plant and the on-site vibration measurement point configuration. Table 6 lists the main parameters of the hydraulic turbine.

20

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

Fig. 11. Hydraulic turbine schematic: 1 Upper guide bearing, 2 Generator rotor, 3 Thrust bearing, 4 Main shaft, 5 Water turbine guide bearing, 6 Spiral case, 7 Guide vane, 8 Turbine runner, 9 Draft tube. Measurement point configuration: V1 At upper guide bearing, V2 At thrust bearing, V3 At water turbine guide bearing.

The rotor radial displacement signals are measured at the upper guide, thrust, and water turbine bearings, and are sampled at 16 Hz during a shut-down transient process without power output. The signals last 50 s. 4.3.2. Signal analysis Fig. 12 shows the time-frequency representations based on GAMD methods. All the prominent components and their time-varying pattern are clearly identified in fine resolution. The rotor rotating frequency f r ðtÞ dominates in all the TFRs and its harmonics up to the fifth order ð2; 3; 4; 5Þf r ðtÞ are also pronounced. Additionally, an almost constant frequency f n around 4.5 Hz appears at relatively high magnitude during the beginning interval from 0 to 15 s. The nonstationary signal time-frequency structure is clearly revealed, though the amplitudes of the higher harmonic components are much lower than that of the fundamental rotating frequency, and even crossings exist between the constant frequency f n and the fifth rotating frequency harmonic 5f r ðtÞ. These findings are helpful to unveil the dynamic properties and evaluate the hydraulic turbine health condition. 5. Conclusions Rotating machinery signals are usually dominated by rotating frequency and its integer or fractional harmonics. These frequencies and their magnitudes are essential to condition monitoring and fault diagnosis. Under time-varying conditions, rotating machinery signals are complex and nonstationary, and their constituent frequencies can even overlap in frequency domain. A reliable time-frequency analysis method is the key to revealing the frequency contents and their time-varying pattern under such circumstances. To this end, effective mono-component decomposition is necessary for accurately calculating instantaneous frequency and constructing TFR of fine resolution. Various AMD methods (including empirical mode composition, local mean decomposition, intrinsic time-scale decomposition, local characteristic scale decomposition, Hilbert vibration decomposition, empirical wavelet transform, variational mode decomposition, and adaptive local iterative filtering) are effective in separating mono-components from signals composed of constant or nearly constant frequencies. However, they often fail when signals are highly nonstationary so that their constituent components overlap in frequency domain due to the filter bank nature (i.e. rectangular partition of time-frequency plane) of these AMD methods. Generalized demodulation is capable of transforming a time-varying instantaneous frequency into a nearly constant frequency, and making it separable from other frequency components even through filters. By exploiting this unique capability of generalized demodulation, a GAMD framework has been proposed in this paper to avoid the key difficulty present in the original AMD methods due to spectral overlaps. Accordingly, AMD methods have been generalized under such a framework to effectively analyze nonstationary signals even with spectral overlaps. The GAMD based TFRs are of fine resolution and are free from both outer interferences of cross-terms and inner interferences of auto-terms (with nonlinearly time-varying frequency). These advantages make the GAMD based time-frequency analysis very effective in the analysis of complex nonstationary signals, especially those signals of rotating machinery under time-varying speed conditions and even with spectral overlaps. The proposed framework is illustrated via numerical simulation, and validated through signal analysis of typical rotating machinery, including a planetary gearbox, an aircraft engine and a hydraulic turbine.

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530 Table 6 Hydraulic turbine parameters. Model

Rated power

Rated speed

ZD(F23)-LH-700

75 MW

88 rpm (1.467 Hz)

Fig. 12. Hydraulic turbine rotor signals and GAMD TFRs.

21

22

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530

Fig. 12 (continued)

CRediT authorship contribution statement Zhipeng Feng: Conceptualization, Methodology, Supervision, Writing - original draft. Xinnan Yu: Formal analysis, Software. Dong Zhang: Formal analysis, Software. Ming Liang: Writing - review & editing. Acknowledgements This work is supported by the National Key R&D Program of China (2018YFC0810500), the National Natural Science Foundation of China (51875034), and the Fundamental Research Funds for the Central Universities (FRF-TP-19-008B1). Comments and suggestions from the reviewers are appreciated very much. References [1] Z. Feng, M. Liang, F. Chu, Recent advances in time-frequency analysis methods for machinery fault diagnosis: a review with application examples, Mech. Syst. Sig. Process. 38 (2013) 165–205. [2] B. Boashash, Time-Frequency Signal Analysis and Processing: A Comprehensive Reference, second ed., Academic Press, London, 2016. [3] Z. Peng, F. Chu, Application of the wavelet transform in machine condition monitoring and fault diagnostics: a review with bibliography, Mech. Syst. Sig. Process. 18 (2) (2004) 199–221. [4] R. Yan, R.X. Gao, X. Chen, Wavelets for fault diagnosis of rotary machines: a review with applications, Signal Process. 96 (2014) 1–15. [5] J. Chen, Z. Li, J. Pan, G. Chen, Y. Zi, J. Yuan, B. Chen, Z. He, Wavelet transform based on inner product in fault diagnosis of rotating machinery: a review, Mech. Syst. Sig. Process. 70–71 (2016) 1–35. [6] F. Auger, P. Flandrin, Improving the readability of time–frequency and time–scale representations by the reassignment method, IEEE Trans. Signal Process. 43 (5) (1995) 1068–1089. [7] I. Daubechies, J.F. Lu, H.T. Wu, Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool, Appl. Comput. Harmon. Anal. 30 (2011) 243–261. [8] Z. Feng, D. Zhang, M.J. Zuo, Adaptive mode decomposition methods and their applications in signal analysis for machinery fault diagnosis: a review with examples, IEEE Access 5 (1) (2017) 24301–24331. [9] 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, Proceedings of the Royal Society of London, Series A, 1998, Vol. 454: 903–995. [10] Z. Wu, N.E. Huang, A study of the characteristics of white noise using the empirical mode decomposition method, Proceedings of the Royal Society of London, Series A, 2004, Vol. 460: pp. 1597–1611. [11] S. Smith, The local mean decomposition and its application to EEG perception data, J. R. Soc. Interface 2 (2005) 443–454. [12] M.G. Frei, I. Osorio, Intrinsic time-scale decomposition: time-frequency-energy analysis and real-time filtering of non-stationary signals, Proceedings of the Royal Society of London, Series A, 2007, 463: pp. 321–342.

Z. Feng et al. / Mechanical Systems and Signal Processing 136 (2020) 106530 [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]

23

J. Zheng, J. Cheng, Y. Yang, A rolling bearing fault diagnosis approach based on LCD and fuzzy entropy, Mech. Mach. Theory 70 (2013) 441–453. M. Feldman, Time-varying vibration decomposition and analysis based on the Hilbert transform, J. Sound Vib. 295 (2006) 518–530. J. Gilles, Empirical wavelet transform, IEEE Trans. Signal Process. 61 (16) (2013) 3999–4010. K. Dragomiretskiy, D. Zosso, Variational mode decomposition, IEEE Trans. Signal Process. 62 (2014) 531–544. A. Cicone, J. Liu, H. Zhou, Adaptive local iterative filtering for signal decomposition and instantaneous frequency analysis, Appl. Comput. Harmon. Anal. 41 (2016) 384–411. S. Olhede, A.T. Walden, A generalized demodulation approach to time-frequency projections for multi component signals, Proc. R. Soc. A 461 (2005) 2159–2179. Z. Feng, F. Chu, M.J. Zuo, Time-frequency analysis of time-varying modulated signals based on improved energy separation by iterative generalized demodulation, J. Sound Vib. 330 (2011) 1225–1243. Z. Feng, X. Chen, M. Liang, Iterative generalized synchrosqueezing transform for fault diagnosis of wind turbine planetary gearbox under nonstationary conditions, Mech. Syst. Sig. Process. 52–53 (2015) 360–375. X. Chen, Z. Feng, Iterative generalized time-frequency reassignment for planetary gearbox fault diagnosis under nonstationary conditions, Mech. Syst. Sig. Process. 80 (2016) 429–444. Z. Feng, X. Chen, Adaptive iterative generalized demodulation for nonstationary complex signal analysis: principle and application in rotating machinery fault diagnosis, Mech. Syst. Sig. Process. 110 (2018) 1–27. J. Shi, M. Liang, D. Necsulescu, Y. Guan, Generalized stepwise demodulation transform and synchrosqueezing for time-frequency frequency analysis and bearing fault diagnosis, J. Sound Vib. 368 (2016) 202–222. Z. Feng, M.J. Zuo, Vibration signal models for fault diagnosis of planetary gearboxes, J. Sound Vib. 331 (2012) 4919–4939. Z. Feng, M. Liang, Y. Zhang, S. Hou, Fault diagnosis for wind turbine planetary gearboxes via demodulation analysis based on ensemble empirical mode decomposition and energy separation, Renew. Energy 47 (47) (2012) 112–126. Z. Feng, M. Liang, Fault diagnosis of wind turbine planetary gearbox under nonstationary conditions via adaptive optimal kernel time-frequency analysis, Renew. Energy 66 (2014) 468–477. J. Antoni, J. Griffaton, H. Andre, et al, Feedback on the surveillance 8 challenge vibration-based diagnosis of a Safran aircraft engine, Mech. Syst. Sig. Process. 97 (2017) 112–144. Z. Feng, F. Chu, Nonstationary vibration signal analysis of hydroturbine based on adaptive chirplet decomposition, Struct. Health Monitor. 6 (4) (2007) 265–279. H. Ohashi, Vibration and Oscillation of Hydraulic Machinery, Avebury Technical, England, 1991.