Enhancement of Fault Feature Extraction from Displacement Signals by Suppressing Severe End Distortions via Sinusoidal Wave Reduction

Abstract

Displacement signals, acquired by eddy current sensors, are extensively used in condition monitoring and health prognosis of electromechanical equipment. Owing to its sensitivity to low frequency components, the displacement signal often contains sinusoidal waves of high amplitudes. If the digitization of the sinusoidal wave does not satisfy the condition of full period sampling, an effect of severe end distortion (SED), in the form of impulsive features, is likely to occur because of boundary extensions in discrete wavelet decompositions. The SED effect will complicate the extraction of weak fault features if it is left untreated. In this paper, we investigate the mechanism of the SED effect using theories based on Fourier analysis and wavelet analysis. To enhance feature extraction performance from displacement signals in the presence of strong sinusoidal waves, a novel method, based on the Fourier basis and a compound wavelet dictionary, is proposed. In the procedure, ratio-based spectrum correction methods, using the rectangle window as well as the Hanning window, are employed to obtain an optimized reduction of strong sinusoidal waves. The residual signal is further decomposed by the compound wavelet dictionary which consists of dyadic wavelet packets and implicit wavelet packets. It was verified through numerical simulations that the reconstructed signal in each wavelet subspace can avoid severe end distortions. The proposed method was applied to case studies of an experimental test with rub impact fault and an engineering test with blade crack fault. The analysis results demonstrate the proposed method can effectively suppress the SED effect in displacement signal analysis, and therefore enhance the performance of wavelet analysis in extracting weak fault features.

1. Introduction

The rotating machinery covers a wide range of engineering applications in the modern industry [1,2]. Owing to the long term service of rotating machinery under hash working conditions, faults are likely to happen on mechanical components, which may lead to economic loss as well as catastrophic accidents [3]. Vibration measurement plays an important role in detecting mechanical faults at an early stage, because it carries abundant information regarding the failure of rotating machines [1,3,4,5,6].

On the other hand, due to complex structures of rotating machinery, the acquired vibration signals often appear as mixtures of multiple vibration modes [7,8]. As it is difficult to identify incipient non-stationary features using methods based on pure time domain or pure frequency domain, signal decompositions are usually required decompose the original measurement [9,10]. This challenging problem attracts attentions from researchers in both academic and engineering fields [11]. In the literature, the multiscale analysis is very popular in studies related to mechanical fault diagnosis and structural health monitoring.

Techniques related to the multiscale analysis can be roughly categorized into two types. The first type is based on predetermined bases, for example the wavelet transform. Wavelet transform, derived from the classical short time Fourier transform, enables the separation of different modes into a number of localized time-frequency subspaces. Currently there are many variations of wavelet transform which have been successfully used to conduct condition monitoring [11,12]. Hemmati extracted bearing characteristic frequencies from the raw acoustic emission signals masked by background noise using wavelet packet transform [13]. Cao employed complex-valued wavelets and convolutional neural network for identifying cutting tool wear states in milling process [14]. Chen extracted gearbox fault signatures based on rational dilation wavelet transform [15]. He conducted the study of automatic feature extraction of induction motors based on parameterized overcomplete wavelet dictionaries [16]. Multiple wavelet dictionaries can also be jointly used to realize sparse representation of incipient fault features [17,18,19,20]. The second type of multiscale analysis is the adaptive decomposition driven by the input data [20,21]. Sun employed empirical mode decomposition to remove stochastic noises in the dynamic measurement of cutting temperatures [22]. Lei investigated early rub-impact fault diagnosis of a heavy oil catalytic cracking machine set using the ensemble empirical mode decomposition (EEMD) [23]. Wang analyzed weak fault features from displacement signals via local mean decomposition-based demodulation (LMD) approaches [24]. LMD was also combined with fast kurtogram to diagnose gearbox faults and bearing faults [25]. In Recent years, the ideal of sparsity is also introduced to adaptive signal decomposition to develop the method of variational mode decomposition (VMD) [26,27]. Recently, many improvements have been added to enhance VMD in the field of condition monitoring [28,29,30,31].

Distortions of end samples in reconstructed signals of multiscale analysis is referred to as the end effect. The existence of end effect has been confirmed in discrete wavelet decompositions. Because we want to maintain the length of the reconstructed signal in each wavelet subspace, boundary extensions are necessary for convolutional operations. The boundary extensions of zero padding, periodic extension, and symmetric extension may create artificial singularities at ends of the analyzed signal. In this paper, we report a special phenomenon of severe end effect in wavelet analysis of the displacement signal. The displacement is explained as the double integral of the acceleration, and the operator of integration assigns high weights for slow-varying contents. Therefore, low frequency contents, such as sinusoidal waves related to the working frequency of the mechanical system, always have extremely high amplitudes. In the practice, we have found that sever distortions on boundary samples are very common in wavelet decompositions of the displacement signal. Such distortions often complicate the process of fault feature extractions. In the past, to suppress the effect of severe end distortion (SED), a few samples on both ends of the reconstructed signals were discarded.

In this paper, we present an insightful study to reveal the cause of the SED effect. It is revealed that if the sampling of the sinusoidal wave does not satisfy the full period sampling (FPS) condition, energy leakages of sinusoidal waves would incur severe distortions on boundary samples, in the form of impulsive features. We present an alternative way to alleviate this negative effect without having to truncate signal lengths. Ratio-based spectrum corrections, serving as postprocessing of the fast Fourier transform (FFT), are used to construct compensations signals of sinusoidal waves of high amplitudes. Spectrum correction methods based on the rectangle windows and the Hanning windows [32] were compared to obtain an optimized estimation of harmonic information of sinusoidal waves. The residual signal is further decomposed by a compound wavelet dictionary developed from the dual tree complex wavelet analysis [33]. This wavelet dictionary, consisting of a dyadic part and a non-dyadic part, provides a more comprehensive time-scaling representation for non-stationary fault features. The statistical indicator of kurtosis is used to select wavelet subspaces of high impulsiveness. The proposed method was applied to case studies of an experimental setup with rub-impact fault and an engineering machine with blade crack fault. Analyzing results showed the extraction of fault features within the actual displacement measurement are significantly enhanced.

The rest of the paper is outlined as below. In Section 2, the FPS condition and the two types of spectrum correction methods are introduced. In Section 3, the construction of the compound wavelet dictionary is explained. In Section 4, the mechanism behind the SED effect during wavelet analysis is investigated. In Section 5, the proposed method is described. In Section 6 and Section 7, two case studies are shown in detail. Finally, some concluding remarks are given in Section 8.

2. Signal Sampling and Discrete Fourier Spectrum of a Sinusoidal Wave

2.1. Fast Fourier Transform of Discrter Digitized Signals

In this paper, a sinusoidal wave is referred to as a simple harmonic wave expressed as

$$sh{w}_c(t)=A_c\cos (2\pi f_{c}t+{\varphi }_c),$$

(1)

where the variables $A_c$, $f_c$ and ${\varphi }_c$ stand for the amplitude, the frequency, and the initial phase at the referential time instant. Let ${\{x(n)\}}_{n=0,\Delta t\dots ,(L-1)\Delta t}$ be a record of signal, of length $L$, sampled from an analog sinusoidal $x(t)$. As shown in Figure 1, it provides a sampling of the associated dynamic process at the time instants of $0,1/f_s,\dots ,(L-1)/f_s$, where $f_s$ denotes the sampling frequency.

The Fourier transform of $x(t)$ is a Dirac impulse function $0.5A_c\delta (f-f_c)+0.5A_c\delta (f+f_c)$ in the frequency domain. Due to the segmentation of finite length, the discrete time Fourier transform for ${\{x(t)\}}_{0\le t<L/f_s}$ becomes a wide band signal. If the shape of the window function is rectangular, the spectrum can be represented as

$$X(f)=DTFT\{x(n)\}=\frac{A_c}{2}·\frac{\sin (\pi T(f-f_c))}{\pi (f-f_c)}{e}^{-\mathrm{j}(\pi Tf-{\varphi }_c)}+\frac{A_c}{2}\frac{\sin (\pi T(f+f_c))}{\pi (f+f_c)}{e}^{-\mathrm{j}(\pi Tf+{\varphi }_c)},$$

(2)

where $T=L/f_s$ and $\mathrm{j}=\sqrt{-1}$ denote the imaginary number. While for ${\{x(n)\}}_{n=0,\dots ,L-1}$, the fast Fourier transform(FFT) can be regarded as a discrete sampling of $X(f)$ at the frequencies of $k/T$, where $k=0,\dots ,L-1$.

In the frequency analysis, the sampling parameters, such as $L$, $f_c$ and $f_s$, have a significant impact on the shape of FFT spectrum of the analyzed signal. We first introduce the definition of full period sampling (FPS) as below.

Condition of Full Period Sampling: If there is a sinusoidal component in a digital signal of length $L$ sampled at the frequency $f_s$, let the harmonic information of the sinusoidal wave be $A_c$, $f_c$, and ${\varphi }_c$. If $N{f}_c/f_s\in {\mathbb{Z}}^{+}$, the sinusoidal wave is said to be full period sampled.

For convenience of discussion, we only discuss the spectrum within the interval of $[0,f_s/2)$. For a sinusoidal wave satisfying the FPS condition, its FFT is also of narrow band. Otherwise it will be recognized as wide band signal in the Frequency domain. As shown in Figure 2, it provides that FFT spectra of sinusoidal waves were affected by the condition of FPS.

2.2. Ratio-Based Spectrum Correction for High-Precision Estimation of Harmonic Information

For a sinusoidal wave not satisfying the FPS condition, it is not possible to retrieve accurate harmonic information from any single spectral bin from the FFT spectrum. To cure this problem, ratio-based spectrum correction (RBSC) methods can be used to obtain high precision estimations of these information using a number of neighboring spectral bins. RBSC can be regarded as post-processing of FFT spectra of signals. Generally, the rectangle window and the Hanning window can be utilized. The former is more suitable for relatively low signal-to-noise ratio (SNR) situations while the latter is more suitable for high SNR situations [32].

2.2.1. Spectrum Correction Based on the Rectangle Window

In the frequency domain, the rectangle window is narrow in main lobe width but slow in side lobe decaying rate. A rectangular window of length $L$ is expressed as

$$w_{\mathrm{rec}}(n)=1forn=0,1,2\dots ,L-1.$$

(3)

The Fourier transform of $w_{\mathrm{rec}}(n)$ is

$$W_{\mathrm{rec}}(w)=\frac{\sin (Lw/2)}{\sin (w/2)}{e}^{j(N-1)w/2},$$

(4)

where the modulus function is $W_{\mathrm{rec}}(w)=\sin (Lw/2)/\sin (w/2)$. As stated above, the frequency resolution of FFT is $\Delta w=2\pi /L$. The angular frequency associated with the $k^{th}$ spectral bin in the FFT spectrum is $w_k=2\pi k/N$. Therefore, the modulus function is approximated as

$$T(k)=\frac{\sin (\pi k)}{\sin (\pi k/L)}\approx \frac{L\sin (\pi k)}{\pi k}.$$

(5)

This approximation is feasible when the value of $L$ is adequately large. The two neighboring spectral bins located within the main lobe of sinusoidal component are used to estimate the harmonic information, as shown in Figure 3. Let the integers $\tilde{k}$ and $\tilde{k}+1$ be the indexes of main lobe spectral bins, $k^{\prime }$ be the index of the rectangular window vertex, and $\Delta k$ be the normalized error between $\tilde{k}$ and $k^{\prime }$. It is easy to know that $-1<\Delta k<0$. It can be derived from the modulus function that

$$\Delta k·T(\Delta k)+(\Delta k+1)T(\Delta k+1)=0.$$

(6)

Therefore, the ratio function can be defined as

$$v=\frac{T(\tilde{k})}{T(\tilde{k}+1)}.$$

(7)

The normalized error $\Delta k$ is as

$$\Delta k=-\frac{1}{1+v}=-\frac{T(\tilde{k}+1)}{T(\tilde{k})+T(\tilde{k}+1)}.$$

(8)

In conclusion, the amplitude, the frequency and the phase can be rectified using the following formula

$$\{\begin{array}{l}{A}_c=\frac{\pi \Delta k·T(k)}{\sin (\pi \Delta k)}\\ f_c=(k+\Delta k)\frac{f_s}{L}\\ {\varphi }_c=\arctan \frac{\mathrm{Im}(k^{\prime })}{\mathrm{Re}(k^{\prime })}-\pi \Delta k\end{array},$$

(9)

where $\mathrm{Re}(·)$ and $\mathrm{Im}(·)$ stand the real part and the imaginary part of the ${k^{\prime }}^{th}$ spectral bin.

2.2.2. Spectrum Correction Based on the Hanning Window

A Hanning window of length $L$ is defined using the mathematical prototype

$$w_{\mathrm{han}}(n)=\frac{1}{2}-\frac{1}{2}\cos \left(\frac{2\pi n}{L}\right)forn=0,1,2\dots ,L-1.$$

(10)

Its Fourier transform is written as

$$W_{\mathrm{han}}(w)=\left\{\frac{1}{2}{W}_{\mathrm{rec}}(w)+\frac{1}{4}\left[W_{\mathrm{rec}}\left(w-\frac{2\pi }{L}\right)+W_{\mathrm{rec}}\left(w+\frac{2\pi }{L}\right)\right]\right\}{e}^{-jNw/2}.$$

(11)

The modulus function of the above equation is

$$T(k)=\frac{\sin (\pi k)}{2\pi k}·\frac{1}{1-k^2}.$$

(12)

Different from the case of rectangle window, the Hanning window is wider in the main lobe but extremely fast in decaying of side lobes. In order to estimate the harmonic information, two neighboring main lobe spectral bins with largest amplitude are also employed, as shown in Figure 4. Let the integers $\tilde{k}$ and $\tilde{k}+1$ be the indexes of two main lobe spectral bins with largest magnitudes, and $\tilde{k}$ be the Hanning window vertex. The normalized error $\Delta k=\tilde{k}-k^{\prime }$. It is easy to know that $-1<\Delta k<0$. It can be derived from the modulus function that

$$(\Delta k-1)T(\tilde{k})+(\Delta k+2)T(\tilde{k}+1)=0.$$

(13)

The ratio function can be calculated as

$$v=\frac{T(\tilde{k})}{T(\tilde{k}+1)}=\frac{\Delta k+2}{1-\Delta k},$$

(14)

and the value of $\Delta k$ is computed as

$$\Delta k=\frac{v-2}{v+1}.$$

(15)

In conclusion, the estimated harmonic information based on the Hanning window is shown as

$$\{\begin{array}{l}{A}_c=\frac{\pi \Delta k·T(\tilde{k})}{\sin (\pi \Delta k)}·2[1-{(\Delta k)}^2]\\ f_c=(k+\Delta k)\frac{f_s}{L}\\ {\varphi }_c=\arctan \frac{\mathrm{Im}(k^{\prime })}{\mathrm{Re}(k^{\prime })}-\pi \Delta k\end{array}.$$

(16)

3. Compound Wavelet Dictionaries Based on Complex-Valued Wavelet Bases

Wavelet transform is an effective tool to characterize non-stationary contents in vibration signals. In discrete wavelet transform, the input signal is decomposed into the sum of a few fast decaying time-frequency atoms. A dual-tree complex wavelet basis (DTCWB) possesses two scaling functions and two wavelet functions, which form two approximate Hilbert transform pairs. DTCWB is equipped with many attractive properties such as approximate analytic, approximate shift-invariance and approximate linear phase. On the other hand, the redundancy factor for DTCWB is 2, which allows fast computational efficiency.

In this paper, wavelet packet decomposition (WPD) based on the DTCWB is utilized in combination with its post-processed subspaces to form a compound wavelet dictionary. As shown in Figure 5, a dyadic part and a non-dyadic part are both employed for dynamic signal decomposition in order to ensure a more comprehensive investigation of time-scale representation of input signals.

3.1. Dyadic Wavelet Packet Decompositon Based on DTCWB

The structure for dual tree complex wavelet packet decomposition, shown in Figure 6a, can be interpreted as two separated conventional wavelet packet transform (WPT) without data flow in between. The filterbank structure of a traditional WPT is also plotted in Figure 6b. A complex-valued scaling function ${\phi }^{ℂ}(t)={\phi }^{\Re e}(t)+\mathrm{j}·{\phi }^{\Im m}(t)$ and a complex-valued wavelet function ${\psi }^{ℂ}(t)={\psi }^{\Re e}(t)+\mathrm{j}·{\psi }^{\Im m}(t)$ are employed in the decompositions. The scaling functions and the wavelet functions in the two parts form approximate Hilbert transform pairs, shown as

$$\{\begin{array}{l}{\psi }^{\Im m}(t)\approx \mathscr{H}\left\{{\psi }^{\Re e}(t)\right\}\\ {\phi }^{\Im m}(t)\approx \mathscr{H}\left\{{\phi }^{\Re e}(t)\right\}\end{array},$$

(17)

where $\mathscr{H}\{·\}$ denotes the Hilbert transform, ‘$\Re e$’ means the functions in the real tree and ‘$\Im \mathrm{m}$’ means the functions in the imaginary tree.

The complex-valued wavelet coefficients $d_{j,k}^{ℂ}(n)=d_{j,k}^{\Re e}(n)+\mathrm{j}{d}_{j,k}^{\Im m}(n)$ are used to reconstruct each wavelet packet in the time domain $r_{j,k}(n)$. A gray code mapping should be used to rearrange the elements in the set of $\{r_{j,k}(n)\}$ such that wavelet packets in $\{w{p}_{j,k}(n)\}$ are listed according to the ascending order of central analysis frequency value. This process is shown in Figure 7.

3.2. Construction of Implicit Wavelet Packets

Although dyadic WPD can realize multiscale analysis of input signals, each wavelet packet has a unique central frequency. In order to enhance the feature extraction ability, implicit wavelet packets (IWPs) can be derived based on the reconstructed dyadic wavelet packets. The signal of $iw{p}_{j,k}(n)$ is obtained by summing two adjacent wavelet packets, shown as

$$iw{p}_{j,k}(n)=w{p}_{i,2k}(n)+w{p}_{i,2k}(n)fork=1,\dots ,2^{j-1}-1.$$

(18)

The dyadic wavelet packets and the IWPs are interleaved to offer a more comprehensive time-scale representation of input signals, which is shown in Figure 8.

3.3. Boundary Extensions in Discrete Wavelet Decomposition

Because discrete wavelet analysis is conducted based on convolutions between the input signal and specific time-frequency atoms, boundary extensions are required such that the filtered signal has the same length with the original signal. Generally, the strategies of zero padding, periodic extension and symmetric extension can be used. Owing the operations of boundary extensions, artificial singularity is often introduced, which leads to the phenomenon of end distortions. In this paper, we focus on investigating end effects produced by the periodic extension during wavelet decomposition.

4. Propagation of Strong Sinusoidal Waves in Multiscale Decompositions

In this section, we will show that the vibration measurement of displacement signal is more sensitive to low frequency components. This often leads to sinusoidal waves of extremely high amplitudes in the sampled signal.

4.1. Sensitivity of Displacement Signal to Low-Frequency Components

There are three types of dynamic signals that can be used for vibration-based condition monitoring and fault diagnosis. Let $\ddot{x}(t)$, $\dot{x}(t)$, $x(t)$ be the acceleration signal, the velocity signal and the displacement signal, respectively, collected from the same sensor point. According the theory of calculus, the three types of signals are related as below.

$$\mathrm{Acceleration}\ddot{x}(t)\underset{\frac{d}{dt}}{\overset{{\displaystyle \int ()dt}}{\rightleftharpoons }}Velocity\dot{x}(t)\underset{\frac{d}{dt}}{\overset{{\displaystyle \int ()dt}}{\rightleftharpoons }}Displacementx(t)$$

(19)

Using Fourier transform, a record of acceleration signal $\ddot{x}(t)$, of finite length, can be represented as the sum of a few sinusoidal waves.

$$\ddot{x}(t)={\displaystyle \sum _{k=0}^{L-1}{A}_k\cos ({\omega }_{k}t+{\varphi }_k)},$$

(20)

where ${\omega }_k=2\pi f_k$. The associated velocity signal and displacement signal can be expressed as

$$\{\begin{array}{l}velocitysignal\dot{x}(t)={\displaystyle \int \ddot{x}\mathrm{d}t}={\displaystyle \sum _{k=0}^{L-1}\frac{A_k}{{\omega }_k}\sin ({\omega }_{k}t+{\varphi }_k)}\\ displacementsignalx(t)={\displaystyle \int \dot{x}\mathrm{d}t}=-{\displaystyle \sum _{k=0}^{L-1}\frac{A_k}{{\omega }_k^2}\cos ({\omega }_{k}t+{\varphi }_k)}\end{array}$$

(21)

For ${\omega }_i<{\omega }_j$, the weights in the integral operations satisfy that ${\omega }_i^{-1}>{\omega }_j^{-1}$ and ${\omega }_i^{-2}>{\omega }_j^{-2}$. It means that the integration operator assigns greater weights to low frequency components and smaller weights to high frequency components. In multiple integrations, the components in the low frequency range will be substantially enhanced. These statements can be illustrated using the following example. We simulate an acceleration signal $\ddot{x}(t)$ as a combination of three sinusoidal waves.

$$\ddot{x}(t)=4{\pi }^{2}100\cos (2\pi 100t+{\varphi }_1)+4{\pi }^{2}500\cos (2\pi 500t+{\varphi }_2)+4{\pi }^{2}900\cos (2\pi 900t+{\varphi }_3)$$

(22)

The dominant component is the sinusoidal wave of 900 Hz. The information regarding its associated velocity signal and displacement signal is listed in

Table 1. The amplitude of each sinusoidal component in the three types of signals.

Acceleration Signal	Frequency	100 Hz	500 Hz	900 Hz
	Amplitude	4π2100	4π2500	4π2900
Velocity Signal	Frequency	100 Hz	500 Hz	900 Hz
	Amplitude	2π	2π	2π
Displacement Signal	Frequency	100 Hz	500 Hz	900 Hz
	Amplitude	1/100	1/500	1/900

. In $\dot{x}(t)$, the three sinusoidal waves have the same amplitude; while in $x(t)$, the sinusoidal wave of 100 Hz becomes the dominant component with the largest amplitude. The information in Table 1 is also illustrated by Figure 9.

4.2. Explanations on Mechanism of the SED Effect

In this subsection, we attempt to give an intuitive explanation on the cause of the SED effect. Assuming that we have two sinusoidal waves of the same amplitude, the Fourier spectra of the two signals are displayed in Figure 10a. As stated in Section 2.1, the sinusoidal wave not satisfying the FPS condition will become a signal of wide band in the Fourier spectrum. On the other hand, for the sinusoidal wave perfectly satisfying the FPS condition, the leakage components are extremely weak in energy.

Because a wavelet function can be interpreted as a digital filter of band-pass (Figure 10b,c), the phenomenon of energy leakage due to the FPS condition will significantly affect the decomposition results. In Figure 10b,c, we illustrate a wavelet function whose theoretical passing band does not include the frequency of the two sinusoidal waves. For the sinusoidal wave not satisfying the FPS condition, the analysis result is a signal with a concentrated energy area in the frequency domain (Figure 10d). The corresponding temporal waveform will become very impulsive on both ends of the signal (Figure 10e). For the sinusoidal wave satisfying the FPS condition, the energy leakage, caused by numerical computation, is also negligible in the frequency domain of the decomposition result (Figure 10f). Therefore, the temporal waveform does not exhibit impulsive features (Figure 10g).

The above explanations have demonstrated the reason for the SED effect in wavelet decomposition. Usually it is not possible for all the sinusoidal waves in the dynamic signal to satisfy the FPS condition simultaneously. Since the SED effect is caused by sinusoidal waves not satisfying the FPS condition, a feasible way to address this problem can be a preprocessing to remove such sinusoidal waves prior to the wavelet decomposition.

4.3. Numerical Simulations

To evaluate the side effect caused by sinusoidal waves not satisfying the FPS condition, we simulate three sinusoidal waves with identical harmonic information of $A_c=1$, $f_c=30\mathrm{Hz}$ and ${\varphi }_c=\pi /6$. The three signals are of different lengths. The sinusoidal wave in $x_1(t)$ satisfies the FPS condition, while $x_2(t)$ and $x_3(t)$ do not. Let the sampling frequency of these signals be 1000 Hz. Their FFT spectra are displayed in Figure 11a.

$$\begin{array}{l}{x}_1(t)=\cos \left(2\pi 30t+\frac{\pi }{6}\right)t\in \left[0,1\right]\\ x_2(t)=\cos \left(2\pi 30t+\frac{\pi }{6}\right)t\in \left[0,\frac{1002}{1000}\right]\\ x_3(t)=\cos \left(2\pi 30t+\frac{\pi }{6}\right)t\in \left[0,\frac{1016}{1000}\right]\end{array}$$

(23)

Two-stage wavelet packet decompositions are performed on these signals. The sinusoidal wave of 30 Hz is expected to be included in $w{p}_{2,1}$ whose theoretical passing band is [0, 125] Hz. We only displayed the reconstructed wavelet packets of $w{p}_{2,2}$, $w{p}_{2,3}$ and $w{p}_{2,4}$, which are shown in Figure 11b–d. Due to satisfying of the FPS condition, the reconstructed signals of high frequency subspaces are small in energies. The maximal absolute values of these wavelet packets are $3\times {10}^{-3}$, $2\times {10}^{-12}$ and $2\times {10}^{-12}$. This is caused by non-ideal filtering properties of wavelet atoms and numerical computational errors. While for $x_2(t)$ and $x_3(t)$, which do not satisfy the FPS condition, strong artificial impulses are found at both ends of $w{p}_{2,2}$, $w{p}_{2,3}$ and $w{p}_{2,4}$.

To investigate the impact of the initial phase in wavelet decomposition of sinusoidal waves, we change the value of ${\varphi }_c$ in Equation (23) continuously in the range $[0,2\pi ]$. The maximal absolute values of $w{p}_{2,2}(t)$, $w{p}_{2,3}(t)$, $w{p}_{2,4}(t)$ are also calculated and are plotted in Figure 12. It is seen, for sinusoidal waves satisfying FPS condition, no matter what the initial phase is, the maximal absolute value of the related wavelet packets ($w{p}_{2,2}(t)$, $w{p}_{2,3}(t)$, $w{p}_{2,4}(t)$) is smaller than $3\times {10}^{-3}$. Recalling that the amplitude of the simulated signal is 1, such results indicate that no severe boundary distortion effect is produced. Such results are obtained because the periodic extension do not create artificial singularities on both ends of the signal $x_1(t)$. While, for the maximal absolute values of wavelet packets decomposed from $x_2(t)$ and $x_3(t)$, the values can be as high as 0.25 and 0.45, which indicates the existence of SED effect.

5. The Proposed Method for Processing Displacement Signals

It is shown in Section 4 that strong sinusoidal waves can be very hazardous for weak fault feature extraction if they are left untreated. Because sinusoidal waves of high amplitudes are so common in analysis of displacement signals, we propose a novel approach based on the combination of ratio-based spectrum correction and DTCWB-based discrete wavelet analysis. The procedure of the approach is mainly divided into three steps.

Step 1. Reduction of Strong Sinusoidal Waves

Investigate FFT spectra of input displacement signals and search spectral bins of high amplitudes. Because the SNR of the signal is not known, ratio spectrum corrections based on the two types of windows should be attempted. Compensation signals are constructed using the corrected harmonic information. A better sinusoidal wave reduction can be obtained by minimizing the energies within the neighborhoods around the spectral bins.

Step 2. Multiscale Analysis of the Residual Signal

Wavelet packet decomposition based on DTCWB is performed on the residual signal. IWPs are derived from post-processing of the reconstructed dyadic wavelet packets.

Step 3. Investigation of information contained in the compound dictionary.

The statistical indicator of kurtosis is employed to evaluate the impulsiveness of the subspace signals. FFT and Hilbert transform are utilized to extract periodic impulsive features of the signals.

A schematic flow chart of the proposed method is illustrated in Figure 13.

6. Case Study of Experiment Test

In order to show the effectiveness of proposed method in actual analysis of displacement signals, an experiment is used in this section. In the test, a fault of rub impact was simulated. The photograph of the test rig is shown in Figure 14. A direct current motor was used to drive the shaft, on which two rotor mass were fixed. Two pairs of eddy current sensors, installed perpendicularly, were used to acquire displacement signals. A rub screw was carefully installed to act as the source of rub impact. In the test, the rotor system was operated at the speed of 2200 revolutions per minute (36.67 Hz). The rub screw contacted with the rotating shaft once in each revolution of the shaft.

The displacement signals were sampled at the frequency of 2000 Hz. The sampling length of each record of signal is 1024. The time domain waveform of one collected signal is shown in Figure 15, where a sinusoidal wave can be observed in the time domain. The sinusoidal component is related to the shaft working frequency (SWF).

The FFT spectrum of the signal is shown in Figure 16. In its zoom-in plot, it is found that the energy of the spectral bin related to the SWF is so strong that higher order harmonic tones of SWF in the higher frequency range are submerged by its slow decaying side lobes. On the other hand, the side lodes related to SWF also appear in the frequency range of [500, 1000] Hz. To confirm the existence of the rub impact fault, we compute the instantaneous frequencies (IFs) of the signal based on Hilbert transform. From the curve displayed in Figure 17, repetitive IF changes are recognized at the interval of 0.0273 s (36.63 Hz).

In order to extract the fault features related to the rub impact fault, the proposed method is used to analyze the signal. The two spectrum correction methods, mentioned in Section 3.2, are employed to estimate the harmonic information related to the SWF. The zoom-in FFT spectra of the residual signals are shown in Figure 18. The corrected harmonic information of two spectrum correction methods are listed in

Table 2. Harmonic information derived from spectrum correction methods.

Spectrum Correction Method	Amplitude	Frequency [Hz]	Phase [rad]
Rectangle window-based spectrum correction	101.5120	36.7102	−0.0198
Hanning window-based spectrum correction	101.5309	36.7110	−0.0158

. The compensation result using the Hanning window is smaller in energy in the neighborhood of the SWF. It actually indicates that the original signal is of high SNR.

The time domain waveform of the optimized residual signal and its FFT spectrum are shown in Figure 19. The absolute maximal value of the signal is about 5, which is much smaller compared with that of the original signal. On the other hand, features of decaying oscillation and impulsiveness are observed. In Figure 19b, the higher order tones related to SWF are significantly enhanced.

We calculate the kurtosis of the wavelet packets and the implicit wavelet packets generated from the processed displacement signal in Figure 20. The wavelet packet of $w{p}_{3,2}(n)$, whose theoretical passing band is [125, 250] Hz, is found to have the largest value of kurtosis. The kurtosis related to this wavelet packet is 4.1667.

To investigate the fault signatures, the time domain waveform of the reconstructed signal is plotted in Figure 21a. A group of impulsive transients, spaced at the interval of 0.0273 s, are detected. In its FFT spectrum, Figure 21b, a few sidebands spaced at 36.71 Hz are observed in the theoretical passing band of the selected wavelet subspace. Moreover, the energy leakage due to non-FHS condition of SWF is successfully suppressed, compared with the zoom-in spectrum in Figure 16. The Hilbert envelope spectrum of the wavelet packet is shown in Figure 21c. The SWF with a few harmonic tones can be found.

The instantaneous curve of the selected wavelet subspace is shown in Figure 22. It is seen that the curve has a periodicity associated with SWF. However, there are two sharp IF changes in each period of shaft rotation. It indicates that the extracted signal may still consist multiple modes. By checking the nearby wavelet subspaces, we selected two WPs, $w{p}_{4,3}(n)$ and $w{p}_{4,4}(n)$, from the fourth WPD stage. Time domain waveforms of the selected wavelet subspaces, shown in Figure 23a,c, reveal that either of the two WPs has a group of impulsive transients occurred at the interval of SWF. In Figure 23b,d, the IF curves of WPs show only one sharp IF change in each period of shaft rotation. The ringing frequencies of the impulses between the IF changes stay comparatively stable. Therefore, we can conclude that vibration modes of more complete physical meaning are extracted in the frequency-scale plane of the compound dictionary.

Figure 24a shows the analysis results using wavelet analysis without reduction of strong harmonic waves. The wavelet subspace with the highest kurtosis is $w{p}_{2,1}(n)$. The value of this WP is more than 350. While $w{p}_{3,2}(n)$ also has a high of kurtosis, which is about 65. In the time domain waveform (Figure 24b), two impulses appear at both ends of the reconstruction signals. In Figure 24c, the even spaced spectral bins are submerges by leakage components of SWF. The same phenomenon can also be found in the Hilbert envelope spectrum (Figure 24d). Comparing Figure 21b with Figure 24c, it is verified that the additional interfering components in the frequency domain is caused by the SWF. Because the sampling of the SWF does not satisfy the FPS condition, the leakage components are also propagated in the signal decomposition process and cause the impulse on both ends of the reconstructed signal.

As a comparison, we also decompose the displacement signal using EEMD. Eight intrinsic mode functions are shown in Figure 25. Although modes of high oscillation can be separated from the sinusoidal wave, the information is not so evident. Among these functions, IMF1 and IMF2 are two modes containing high frequency contents. However, in these two modes, no explicit periodicity of impulsive transients can be directly identified. The zoom-in plot of the IMF2 is shown in Figure 26a and distortions at end samples are also found. In the FFT spectrum (Figure 26b), a few side bands are observed with more interfering noises and it causes the corruption of impulsive transients in the Hilbert envelope spectrum (Figure 26c). The IFs of IMF2, shown in Figure 26d, is entirely irregular because the continuity of the curve is poor.

7. Case Study of an Engineering Application

In this section, an engineering application is used to demonstrate the effectiveness of the proposed method. Figure 27 provides the schematic diagram of the investigated rotating machine.

The booster fan was installed in a blast furnace gas power plant. An alternated current motor was used to drive the centrifugal compressor. After eight years of service, an accident occurred on the compressor. In the accident, a blade fell off from the compressor during operation. The blade of failure knocked against the volute and broke through it.

In the power plant, vibrations of the shaft in Figure 27 were regularly collected and stored in the data server. Although the signals were recorded, the on-site staff did not have the sophisticated skills to interpret them. To investigate the potential fault signatures caused by the blade crack before the accident, we studied the displacement measurements using the proposed method. The sampling frequency and sampling length of the signal were 2560 Hz and 5120 respectively. The rotation speed of the booster fan was controlled at 993 revolutions per minute. The SWF is calculated at 16.55 Hz.

The time domain waveform of the signal is shown in Figure 28. The sinusoidal wave related to the SWF is again the most dominant component. The FFT spectrum of the signal is shown in Figure 29. We mark the higher harmonic tones of SWF by SWF2, SWF3 and SWF4 in the Figure. The spectrums of SWF2, SWF3, SWF4 are significantly affected by side lobes of SWF in the frequency domain. The IF curve of the original displacement signal is displayed in Figure 30. Although there is some vague information indicating a periodic change of the curve, it is still not easy to recognize it explicitly.

To reduce the strong sinusoidal waves, ratio based spectrum correction methods were applied on the signal. Figure 31 show the FFT spectra of the residual signals. Because this engineering signal has more noises in the low frequency range, these two spectrum correction methods attain similar results. For reductions of SWF, SWF2, and SWF3, the Hanning window based method is better. While, for reduction of SWF4, the rectangular window based method is better due to relatively low SNR condition in the signal.

In this case, either spectrum correction method can be selected. The corresponding residual signal and its FFT spectrum are shown in Figure 32. There are still trend terms in the low frequency range but their energies are relatively small. Besides, the high frequency components are significantly enhanced.

The wavelet analysis was applied to decompose the residual signal. Kurtosis values of the wavelet packets and IWPs generated from the residual signal are calculated and shown in Figure 33.

The wavelet packet of $iw{p}_{4,4}(n)$, whose theoretical passing band is [560, 720] Hz, is found to have the largest value of kurtosis. The kurtosis related to this wavelet packet is 16.157. The time domain waveform of $iw{p}_{4,4}(n)$ is shown in Figure 34a, in which a group of periodic transients spaced at the interval are found. The occurrence frequency of the features is calculated at 0.0602 s (16.613 Hz). This is much closed to the value of the SWF of the booster fan. In the FFT spectrum (Figure 34b), a few sidebands spaced at SWF can be observed. In the Hilbert envelope spectrum (Figure 34c), the associated spectral bins can also be clearly detected. The IF curve of the wavelet space is plotted in Figure 35, in which the ringing frequency of the impulsive transients stay relatively stable. The above results show the proposed method can successfully extract vibration mode related to the blade crack fault.

As comparison, the original displacement signal is decomposed directly by the wavelet analysis. The distribution of kurtosis is given in Figure 36a. The wavelet packet $w{p}_{2,1}$ has the largest kurtosis of almost 900. The time domain waveform of the implicit wavelet $iw{p}_{4,4}(n)$ is investigated and two impulses are found on both ends of the signal Figure 36b. Comparing the FFT spectrum in Figure 36c with that in Figure 34b, it is observed the frequency counterparts of the impulses are submerged by a flat side lobe caused by the leakage of the sinusoidal waves related to SWF. In Figure 36c, it is noticed that the frequency range of the side lobe perfectly coincides with the theoretical passing band of the implicit wavelet packet ([560, 720] Hz). This has verified that leakage components due to SWF not satisfying the FPS condition has incurred the impulsive features in the wavelet packets. In the Hilbert envelope spectrum (Figure 36d), no spectral bins related to SWF can be seen. In the results by the proposed method, the sinusoidal waves related the SWF of the system are suppressed, therefore the associated leakage components are very small in energy and do not cause severe distortions.

The decomposition results of the original displacement signal by EEMD are given in Figure 37. The nonstationary features are mainly found in IMF1, IMF2 and IMF3. However, it is unable to identify periodic impulsive transients in these signals.

8. Conclusions

Aiming at suppressing the SED effect in multiscale analysis of displacement signals, an approach based on the combination of spectrum correction and the compound wavelet dictionary is proposed. The major findings are concluded as below.

(1) The displacement measurement is very sensitive to components in the low frequency range. In actual dynamic measurements, sinusoidal waves are usually large in amplitude and therefore become the dominant components in the displacement signal. If the digitization of sinusoidal wave does not satisfy the condition of FPS, a severe problem of energy leakage emerges in the Frequency spectrum of the signal. The leakage components will be propagated during multiscale decompositions and cause the impulsive features on both ends of wavelet packets. Theoretical explanations and numerical simulations were used to show the presence of strong sinusoidal waves not satisfying the FPS condition will significantly intensify the SED effect.

(2) Because strong sinusoidal waves can be a major source of the SED effect, spectrum correction methods are employed to construct optimized compensation signals. The rectangular window-based method and the Hanning windows-based method are combined to derive an optimized estimation of the harmonic information of sinusoidal waves. After reductions of the strong sinusoidal waves in the original measurements, the SED effect has been effectively suppressed. A compound wavelet dictionary, consisting dyadic wavelet packets and implicit wavelet packets, is utilized to decompose the residual signal. The wavelet dictionary provides a more comprehensive time-scale representation for nonstationary features. The preprocessing of sinusoidal wave reduction and the compound wavelet dictionary are combined to enhance the extraction of periodic impulsive features in displacement signals.

(3) The performance of the proposed method is investigated using an experimental test and an engineering case study. In the two investigated case studies, the original signals possess strong sinusoidal waves related to the shaft working frequencies of the systems. In the analysis, wavelet decompositions without preprocessing of sinusoidal wave reduction inevitably caused severe distortions on both end samples in resultant wavelet packets. Many interfering contents irrelevant to the actual machinery faults are also of high kurtosis values. By using the proposed method, the artificial impulsive features in the wavelet packets are successfully suppressed. The analysis results are also compared with those of the EEMD. It is validated the proposed method can effectively suppress the SED effect without having to truncate some samples on decomposed signals.

In this paper, we have focused on discussing the effectiveness of the proposed method when periodic extension is adopted in wavelet decompositions. However, it can be verified the proposed method works well no matter what kind of boundary extension is selected.