Full-oscillatory components decomposition from noisy machining vibration signals by minimizing the Q-factor variation

Abstract

Generally, the machining vibration frequency spectrum is dominated by the tooth cutting frequency and its harmonics, the part structure and its natural frequency, and the spindle-tool subsystem natural frequency, exhibiting full-oscillatory behaviour. In order to identify the machining status, especially for those thin-walled workpiece machining, the on-machine detected monitoring signals with noise should be decomposed precisely. Actually, the signals’ inherent characteristics, such as the Q-factor, could be employed. In this article, decomposition of the full-oscillatory components from noisy machining vibration signals by minimizing the Q-factor variation is presented. The Q-factor will be calculated using quadratic interpolation of linear prediction coefficients. On this basis, the measured signals can be decomposed into high-, low- and residual-oscillatory signal components using the sparsity-enabled signal analysis. Furthermore, the signal decomposition process is repeated iteratively until the minimization of the Q-factor variation. Finally, the simulation and the thin-walled machining experiments were designed. From comparison of the signal decomposition results with the wavelet packet transform (WPT), it was shown that the signal decomposition accuracy and reliability using the proposed strategy has been improved significantly.

Keywords

Full-oscillatory machining vibration monitoring Q-factor signal decomposition

Introduction

High-performance parts with thin-walled structures to be machined have been widely employed in aerospace fields. Unfortunately, machining vibration is always inevitable due to its significant low rigidity (Quintana and Ciurana, 2011), which imposes a negative effect on the productivity including poor surface quality and low dimensional accuracy. Furthermore, it generates excessive noise and aggravates cutting-tool wear. To maintain a high level of machining efficiency and reliability, it is necessary to identify and predict the machining status using the on-machine detected monitoring signals. Generally, the vibration frequency spectrum is dominated by the tooth cutting frequency and its harmonics, the part structure and its natural frequency, and the spindle-tool subsystem natural frequency (Davies and Balachandran, 2000; Kolluru and Axinte, 2013). In fact, the measured machining vibration signals mainly contain three main components: chatter-related signals characterized by low oscillation, periodic cutting excitation signals characterized by high oscillation and noise signals (i.e. other extraneous signals with stochastic oscillation). Thus, it essentially exhibits full-oscillatory behaviours. Generally, the original frequency and/or time domain signals without decomposition are directly used to estimate and predict the current machining status. However, only the rough machining status can be obtained. The reliability for confirming the machining vibration diagnosis is very poor and sometimes leads to error results. Decomposing the detected signals and finding the status-sensitive components can supply an efficient solution. However, it is difficult to decompose the full-oscillatory signals into different components and extract the periodic transient components related to machining chatter.

As the key point of machining stability monitoring, the chatter feature extraction from vibration signals has been studied in the last decades. In a traditional way, frequency-based analysis and filtering have been widely employed as a fundamental tool in signal processing. Furthermore, some advanced signal process techniques have been developed for the fault feature extraction in different fault diagnosis areas, such as short-time Fourier transformation (STFT) (Klein et al., 2001), wavelet transformation (WT) (Peng and Chu, 2004), ensemble empirical mode decomposition (EEMD) (Cao et al., 2015; Lei et al., 2013), blind source separation (BSS) (Guo et al., 2014) and tunable Q-factor wavelet transform (TQWT) (Selesnick, 2011a). Moreover, the aforementioned methods have also been combined together systematically (Bin et al., 2012; Han et al., 2015; Peng et al., 2005; Shao et al., 2011) or improved (Hoell and Omenzetter, 2016) to make full use of advantages and avoid the disadvantages. In these methods, the TQWT performs excellent signal analysis especially for signals with different oscillatory behaviours by adjusting the Q-factor, and redundancy to match the oscillatory behaviour of the wavelet basis with that of the signal of interest. For the signal that is a mixture of multiple signal components with different oscillatory properties, Selesnick (2011b) proposed a new sparsity-enabled signal analysis method with two distinct Q-factor transforms simultaneously. Cai et al. (2013) proposed a new fault feature extraction technique for gearboxes using sparsity-enabled signal analysis method and investigated the effect of two parameters pertinent to morphological component analysis (MCA) and split augmented Lagrangian shrinkage algorithm (SALSA), the Lagrange multiplier and the penalty parameter. A novel method for the original signal of a rolling bearing’s early weak fault was decomposed by EEMD first and the intrinsic mode function (IMF) with the biggest kurtosis index value selected and handled by the TQWT subsequently is proposed (Wang et al., 2014). Lastly, the envelope demodulation method was applied on the low Q-factor transient impact component and a satisfactory extraction result was obtained. He et al. (2013) proposed a denoising method based on the tunable Q-factor using neighbouring coefficients. Zhang et al. (2014) proposed an application of wavelet transform with a tunable Q-factor to analysis of non-stationary harmonics.

However, the determination of the Q-factor lacks quantitative criteria, which has increased the instability of the decomposition result and limited the TQWT application. In the research mentioned above, the determination of the Q-factor is based on two strategies. One, counting the oscillation number of the periodic transient signals becomes a way of determining the Q-factor. In a sense, the Q-factor determines the oscillation times of the wavelet waveform (the wavelet with Q-factor of Q sustained about 2Q times oscillations; Selesnick, 2011b). However, in engineering, the waveform of the chatter components is unknown and the Q-factor is rarely determined in advance. As an alternative, for calculating convenience, the Q-factor has been determined as a compromise between calculating efficiency and frequency resolution in a high-frequency band. The value of the Q-factor reflects the frequency resolution and the high-oscillatory components need a higher frequency resolution than the low-oscillatory components. For the low-frequency band, the frequency resolution can be preserved regardless of the value the Q-factor and redundancy, but for the high-frequency band, higher frequency resolution can be achieved by increasing Q-factor and redundancy. As a result, more decomposition levels are needed to cover a wide frequency range, which will in turn reduce the calculating efficiency inevitably. The trade-off between frequency resolution and computational efficiency limits the performance of the TQWT. Luo et al. (2013) proposed a kurtosis-guided adaptive demodulation technique for bearing fault detection based on TQWT. From two steps of merging of the band-pass filters of the TQWT, the optimal filter for the Hilbert demodulation analysis is obtained. However, only the maximum allowable decompositions for a given (Q, R) combination have been used in this method and the Q-factor was given by experience.

Aiming at this problem, a comprehensive strategy by combining the sparsity-enabled signal analysis method and precise calculation of Q-factor is developed to decompose full-oscillatory components from noisy machining vibration signals. Based on the resonance peak calculating method, centre frequency and band width, which determines the value of the Q-factor, can be precisely determined utilizing quadratic interpolation of linear prediction coefficients (LPC). Then, the sparsity-enabled signal analysis method is used iteratively to decompose full-oscillatory components from noisy machining chatter signals into periodic cutting components and transient components, until the minimization of the Q-factor variation.

The paper is organized as follows. The next section describes the principle of TQWT, and the two main aspects of the sparsity-enabled signal analysis, MCA and SALSA is introduced. An iterative algorithm is designed, and simulation study and machining experiments are conducted.

Principle of the TQWT

The TQWT can be realized through a series of band-pass filters. The analysis and synthesis filter banks with J-stage wavelet transform are illustrated in Figure 1 (Cai et al., 2013). TQWT can be implemented by applying inverse over-sampling filter banks with real-value sampling parameters. As shown in Figure 1, the transform consists of a sequence of two-channel filter banks, with the low-pass output of each filter bank being the input of the successive filter bank. By applying two channel filter banks iteratively on the low-pass channel iteratively followed by the low-pass scaling parameter $α$ and the high-pass scaling parameter $β$ , the TQWT can decompose the signal into J+1 sub-bands: the high-pass filter output signal of each filter bank and the low-pass filter output signal of the final filter bank. The relationship between the scaling parameters ( $α$ and $β$ ) and parameters of TQWT (the Q-factor Q and the redundancy r) can be represented as (Selesnick, 2011a)

β = \frac{2}{Q + 1} α = 1 - \frac{β}{r}

(1)

The frequency responses $H_{0} (ω)$ of low-pass filter and $G_{0} (ω)$ of high-pass filter are defined as (Selesnick, 2011a)

H_{0} (ω) = {\begin{matrix} 1 \\ θ (\frac{ω + (β - 1) π)}{α + β - 1}) \\ 0 \end{matrix} \begin{matrix} | ω | \leq (1 - β) π \\ (1 - β) π < | ω | < α π \\ α π \leq | ω | < π \end{matrix}

(2)

G_{0} (ω) = {\begin{matrix} 0 \\ θ (\frac{α π - ω}{α + β - 1}) \\ 1 \end{matrix} \begin{matrix} | ω | \leq (1 - β) π \\ (1 - β) π < | ω | < α π \\ α π \leq | ω | < π \end{matrix}

(3)

where $θ (ω) = \frac{1}{2} (1 + \cos ω) \sqrt{2 - \cos ω}, | ω | \leq π$ is the Daubechies frequency response with two vanishing moments. Obviously, the low-pass filter $H_{0} (ω)$ and high-pass filter $G_{0} (ω)$ meet the requirement of ${| H_{0} (ω) |}^{2} + {| G_{0} (ω) |}^{2} = 1$ for perfect reconstruction. The parameters $α$ and $β$ determine the property of frequency response function of low-pass and high-pass filter banks, which meet the requirements of the following (Selesnick, 2011a).

{\begin{matrix} 0 < α \leq 1 \\ 0 < β \leq 1 \\ α + β > 1 \end{matrix}

(4)

The first two formulas ensure that the wavelet transform will not be overly redundant and the last one ensures the strictly oversampled for perfect signal reconstruction and filter banks accurate response.

The TQWT can be implemented with the Radix-2 Fast Fourier Transforms. With the defined linear matrix $Φ$ , the coefficients can be calculated with the equation $w = Φ X$ and represented as $w = {w^{{1}}, w^{{2}}, \dots, w^{{J + 1}}}$ , where the wavelet coefficients $w^{{j}}$ is the set $[d_{j} (1), d_{j} (1), \dots, d_{j} (n)]$ for $j = 1, 2, \dots, J$ and the scaling coefficients $w^{{J + 1}}$ is the set $[a_{J + 1} (1), a_{J + 1} (2), \dots, a_{J + 1} (n)]$ .

Sparsity-enabled signal analysis method

The sparsity-enabled signal analysis method is based on signal resonance, rather than on frequency or scale. In this method, MCA is used to construct the objective function and decompose the signal non-linearly, and a SALSA is used to minimize the objective function of MCA. Finally, the signal components with different oscillatory properties can be obtained.

Morphological component analysis (MCA)

For the full-oscillatory signal s, which is the sum of multiple oscillatory components and noise component, $s = s_{H} + s_{L} + s_{n}$ , the high-oscillatory component $s_{H}$ and low-oscillatory component $s_{L}$ can be determined by performing MCA to obtain the sparse coefficient vectors $w_{H}$ and $w_{L}$ with the transforms $Φ_{H}^{*}$ and $Φ_{L}^{*}$ such that

s = Φ_{H}^{*} w_{H} + Φ_{L}^{*} w_{L} + n

(5)

where n is the error component that cannot be expressed with inverse TQWT ( $Φ_{H}^{*}$ and $Φ_{L}^{*}$ ). Following this, the basis pursuit denoising can be further expressed as

{w_{H}^{o p t}, w_{L}^{o p t}} = \underset{w_{H}, w_{L}}{\arg \min} \frac{1}{2} {‖ x - Φ_{H}^{*} w_{H} - Φ_{L}^{*} w_{L} ‖}_{2}^{2} + λ_{H} {‖ w_{H} ‖}_{1} + λ_{L} {‖ w_{L} ‖}_{1}

(6)

where $λ = [λ_{H}, λ_{L}]^{T} = [θ_{H} ‖ Ψ_{H}^{{k}} ‖_{2}, θ_{L} ‖ Ψ_{L}^{{k}} ‖_{2}]$ is the Lagrange multiplier that controls the trade-off between sparsity of ${w_{H}, w_{L}}$ and reconstruction fidelity. Considering that the energy of the wavelet in different decomposition level is different, Equation (6) can be overwritten with a great flexible form:

{w_{H}^{o p t}, w_{L}^{o p t}} = \underset{w_{H}, w_{L}}{\arg \min} \frac{1}{2} {‖ x - Φ_{H}^{*} w_{H} - Φ_{L}^{*} w_{L} ‖}_{2}^{2} + \sum_{i = 1}^{J_{1} + 1} λ_{H, i} {‖ w_{H, i} ‖}_{1} + \sum_{k = 1}^{J_{2} + 1} λ_{L, k} {‖ w_{L, k} ‖}_{1}

(7)

where the Lagrange multipliers ${λ_{H, i}}^{T} = [η_{H, i} ‖ Ψ_{H}^{{i}} ‖_{2}]^{T}$ and ${λ_{L, k}}^{T} = [η_{L, k} ‖ Ψ_{L, k}^{{k}} ‖_{2}]^{T}$ are proportional to the $l_{2}$ -norm of the wavelet coefficients $w_{H}^{i}$ and $w_{L}^{k}$ , namely $‖ Ψ_{H}^{{i}} ‖_{2}$ and $‖ Ψ_{L}^{{k}} ‖_{2}$ . The scalar parameters $η_{H}$ and $η_{L}$ are associated with the signal and the estimation error. For the energy of certain components, which is expected to be negligible, the corresponding scalar parameter of $θ_{H}$ or $θ_{L}$ is introduced to set a higher value in order to minimize the energy of certain components for the optimization problem expressed in Equation (7). Then the objective function can be expressed as (Cai et al., 2013)

O (w_{H}, w_{L}) = \frac{1}{2} {‖ x - Φ_{H}^{*} w_{H} - Φ_{L}^{*} w_{L} ‖}_{2}^{2} + θ_{H} \sum_{i = 1}^{J_{1} + 1} λ_{H, i} {‖ w_{H, i} ‖}_{1} + θ_{L} \sum_{k = 1}^{J_{2} + 1} λ_{L, k} {‖ w_{L, k} ‖}_{1}

(8)

Split augmented Lagrangian shrinkage algorithm

For the mathematical optimization problem of the form in Equation (8), the SALSA proved to be an effective way to solve the optimization problem. The SALSA contains two aspects: variable splitting and an augmented Lagrangian algorithm. Variable splitting mainly solves the optimization problem of the sum of two or more functions. By variable splitting, one of the multiple variables can be substituted and the unconstrained optimization problem, which is hard to solve, can be transformed into a constrained optimization problem, which is easy to solve. The optimization of the objective function in Equation (8) can be written in the following form:

min_{w, u} (f_{1} (u) + f_{2} (w)) such that w = u

(9)

where $f_{1} (u) = \sum_{i = 1}^{J_{1} + 1} λ_{H, i} {‖ w_{H, i} ‖}_{1} + \sum_{k = 1}^{J_{2} + 1} λ_{L, k} {‖ w_{L, k} ‖}_{1} and f_{2} (w) = \frac{1}{2} ‖ x - Φ_{H}^{*} w_{H} - Φ_{L}^{*} w_{L} ‖_{2}^{2}$ .

With the alternating Lagrangian algorithm (ALM), more particularly, the alternating direction method of multipliers (ADMM), the constrained optimization problem in Equation (9) can be iteratively solved, as shown in Table 1.

Table 1.

Alternating direction method of multipliers.

① Initialize:

μ > 0, d

②

u \in \arg min λ_{H} ‖ u_{H} ‖_{1} + λ_{L} ‖ u_{L} ‖_{1} + \frac{μ}{2} ‖ u - w - d ‖_{2}^{2}

③

w \in \arg min \frac{1}{2} ‖ x - Φ_{H}^{*} w_{H} - Φ_{L}^{*} w_{L} ‖_{2}^{2} + \frac{μ}{2} ‖ u - w - d ‖_{2}^{2}

④

d = d - (u - w)

⑤ Repeated

Absolutely, the solution to the minimization problem for ② in Table 1 can be regarded as an $l_{1}$ -norm regularized denoising problem and the optical u of $λ ‖ u ‖_{1} + \frac{u}{2} ‖ u - y ‖_{2}^{2}$ can be obtained from the soft-threshold rule with threshold $soft (y, λ / μ)$ , which is defined as

s o f t (y, λ / μ) = y \times \max (0, 1 - λ / (μ | y |)), y \in ℂ, λ \in R_{+}

(10)

The minimization problem for ③ in Table 1 can be considered a constrained least squares problem to which the solution can be expressed explicitly in matrix form. According to the analysis above, the algorithm can be expressed as shown in Table 2 by simplifying the computational redundancy and slightly rearranging the equations.

Table 2.

Split augmented Lagrangian shrinkage algorithm.

① Initialize:

k = 0, μ > 0, d_{0}

and

w_{0}

②

u^{(k + 1)} \leftarrow soft (w^{(k)} + d^{(k)}, λ / μ) - d^{(k)}

③

d^{(k + 1)} \leftarrow \frac{1}{μ + 2} Φ (x - Φ^{*} u^{(k + 1)})

④

w^{(k + 1)} \leftarrow d^{(k + 1)} + u^{(k + 1)}

⑤ Repeated

Iterative algorithm

Q-factor precise calculation

The basic principle of the linear prediction analysis is that the current sampling value of the signal can be approximated with a linear combination of some original sampling points. By minimizing the quadratic sum of the difference between the actual sampling value and the predicted value, the only set of prediction coefficients can be determined. Except for the linear prediction function, the LPC can also offer an excellent signal model. The properties of the signal model, including centre frequencies and frequency bandwidth, can be described by the LPC.

According to the principle of the linear predictive analysis, the signal can be represented as the output of a signal model with the input $u (n)$ and the properties of the signal can be described with the parameter of the signal model. When the output signal $x (n)$ is deterministic, the input signal $u (n)$ is a unit pulse series and when $x (n)$ is a stochastic signal, the input signal $u (n)$ is a white noise sequence.

The transfer function of the model can be represented as the following:

H (z) = G \frac{1 + \sum_{l = 1}^{q} b_{l} z^{- l}}{1 - \sum_{k = 1}^{p} a_{k} z^{- k}}

(11)

where p and q are the orders of the model, and $b_{l}$ and $a_{k}$ are the model parameters. G is the gain coefficient of the model. Considering the calculating convenience in which the parameter estimation of all-pole model is a process of solving linear equations, the model can be represented as the auto-regressive model (AR model) by supposing the parameter of $b_{l}$ to be zero. In the actual calculation, the ARMA model, namely the pole-zero model, can also be regarded as the all-pole model when the order of p is high enough.

For the AR model, the signal can be represented with the difference equation as follows:

x (n) = \sum_{k = 1}^{p} a_{k} x (n - k) + Gu (n)

(12)

Equation (12) indicates that the current output signal derives from the current input signal and linear combination of some original output signals. Then the transfer function of all-pole model can be formulated as follows:

H (z) = \frac{G}{1 - \sum_{k = 1}^{p} a_{k} z^{- k}}

(13)

The basic problem of the linear prediction is the determination of the predictor coefficients ${a_{k}}$ and minimization of the prediction error in certain criterion. The model parameter of ${a_{k}; k = 1, 2, \dots p}$ can be obtained through the Levinson–Durbin algorithm. By substituting $z^{- 1} = \exp (- j ω T)$ or $z^{- 1} = \exp (- j 2 π f / f_{s})$ ( $f_{s}$ is the sampling frequency) into the transfer function and modulus operation, the module value of the power spectrum denoted by $P (f)$ can be represented in Equation (14). Then the power spectrum response can be obtained through the FFT of the prediction coefficients.

P (f) = {| H (f) |}^{2} = \frac{G^{2}}{{| 1 - \sum_{k = 1}^{p} a_{k} \exp (- j 2 π kf / fs) |}^{2}}

(14)

For every peak in the power spectrum response curve shown in Figure 2, the centre frequency and corresponding bandwidth can be calculated with the parabolic interpolation method. Supposing the frequency of a certain local peak to be $m \times Δ f$ ( $Δ f$ is the frequency space in the spectrum) and the adjacent frequencies to be $(m - 1) \times Δ f$ and $(m + 1) \times Δ f$ , the responding power spectrum amplitudes are $P_{(m - 1)}$ , $P_{(m)}$ and $P_{(m + 1)}$ . Then the precise centre frequency $f_{i}$ and band width $B w_{i}$ can be calculated with the quadratic equations of $a_{i} λ^{2} + b_{i} λ + c_{i}$ .

For calculating convenience, the original point of the frequency axis is moved to m and the local peak frequency $m \times Δ f$ is set to zero. The power spectrum amplitude is divided with equal frequency space $Δ f$ . The power spectrum value $P_{(m - 1)}$ , $P_{(m)}$ and $P_{(m + 1)}$ responding to $- Δ f$ , 0 and $+ Δ f$ can be represented in the following formulation:

{\begin{matrix} P_{(m - 1)} = a_{i} Δ f^{2} - b_{i} Δ f + c_{i} \\ P_{(m)} = c_{i} \\ P_{(m + 1)} = a_{i} Δ f^{2} + b_{i} Δ f + c_{i} \end{matrix}

(15)

The coefficients of $a_{i} λ^{2} + b_{i} λ + c_{i}$ can be obtained as $a_{i} = (P_{(m + 1)} + P_{(m - 1)} - 2 \cdot P_{(m)}) / Δ f^{2}$ , $b_{i} = (P_{(m + 1)} + P_{(m - 1)}) / Δ f^{2}$ and $c_{i} = P_{(m)}$ .The local maximum of the power spectrum can be calculated through the derivation of the quadratic equation $\frac{d}{d λ} (a_{i} λ^{2} + b_{i} λ + c_{i}) = 0$ . Considering the actual frequency space $Δ f$ , the centre frequency $f_{i}$ can be represented as

f_{i} = λ_{pi} Δ f + m Δ f = (- b_{i} / 2 a_{i} + m) Δ f

(16)

where $λ_{pi} = - b_{i} / 2 a_{i}$ .

The responding power spectrum value $P_{p}$ is

P_{si} = a_{i} λ_{pi}^{2} + b_{i} λ_{pi} + c_{i} = \frac{{b_{i}}^{2}}{4 a_{i}} - \frac{{b_{i}}^{2}}{2 a_{i}} + c_{i} = - \frac{{b_{i}}^{2}}{4 a_{i}} + c_{i}

(17)

The bandwidth $B w_{i}$ refers to the distance between two frequency points on which the maximum of the power spectrum amplitude reduced by 3 dB, as shown in Figure 3, which can be represented by the following formula.

B w_{i} = - \frac{\sqrt{{b_{i}}^{2} - 4 a_{i} (c_{i} - 0.5 P_{si})}}{a_{i}} \times Δ f

(18)

The power spectrum curve after quadratic interpolation of LPC is shown in Figure 3. As shown in Figure 3, there are four centre frequencies, $f_{1}$ , $f_{2}$ , $f_{3}$ and $f_{4}$ , and four band widths $B w_{1}$ , $B w_{2}$ , $B w_{3}$ and $B w_{4}$ , respectively. Each one corresponds to a certain spectral line.

Algorithm design

In this paper, a novel signal decomposition method based on joint application of the sparsity-enabled signal analysis method and precise calculation of Q-factor is employed to extract chatter components from the noisy measurement signals. The value of the Q-factor can be updated with quadratic interpolation of LPC, which is obtained from the signal model founded by the measured signal.

The flow chart of the proposed feature extraction technique is illustrated in Figure 4.

Step 1: With the obtained vibration signal $S_{m}$ , the signal model is constructed and LPC are obtained, and the centre frequencies and band width can be obtained with quadratic interpolation of LPC. Then the Q-factors can be obtained.

Step 2: Substituting the obtained Q-factors into TQWT and with MCA and SALSA, the original vibration signal are separated into high-oscillatory, low-oscillatory and the residual components.

Step 3: Summing the low-oscillatory and the residual components as the input signal of the new iterative operation, the new signal model is constructed and the new Q-factors are updated.

Step 4: The iterative operation will continue until the difference between the newly updated Q-factors and the previous one is limited within a certain degree. Then the final high- and low-oscillatory components can be obtained.

Validations

Simulation study

In this section, a simulation signal is constructed to illustrate the effectiveness of the proposed method for signal decomposition, which consists of amplitude-modulated signal as the high-oscillatory components and periodic single-sided transient signal as the low-oscillatory components. Furthermore, the noise is added into the signal to simulate the noisy environment. To illustrate further the effectiveness of the proposed method, the correlation coefficients between the simulation signal and decomposed signal have been calculated.

Considering the characteristics of the signals measured in the thin-walled machining process, the simulation signal $s (t)$ can be formulated as

{\begin{matrix} s (t) = s_{1} (t) + s_{2} (t) + s_{3} (t) = (1 + 0.07 \cos (2 π f_{m} t_{2})) \cos (2 π f_{c} t_{2}) + e^{- 40 t_{1}} \sin (2 π f_{0} t_{2}) + A \times N \\ t_{1} = mod ((0 : 1023) / 1024, 1 / 10), t_{2} = (0 : 1023) / 1024 \end{matrix}

(19)

where the first part of the simulation signal $s_{1} (t)$ is the amplitude-modulated signal with the frequencies $f_{m} = 15$ Hz and $f_{c} = 100$ Hz, and the second part $s_{2} (t)$ represents the periodic single-sided transient signal with the frequency $f_{0} = 230$ Hz. N is the white noise. The third part of the simulation signal is the white noise weighted by $A = 0.3$ . The sampling frequency is 1024 Hz and the length of the simulation signal is 1024. The waveform of the simulation signal is shown in Figure 5, in which the amplitude-modulated part $s_{1} (t)$ , periodic single-sided transient part $s_{2} (t)$ and noisy part are illustrated from (a) to (d).

The proposed method is applied to decompose the simulation signals s(t) into the amplitude-modulated component and periodic transient component. In order to guarantee the model of the signal to be all-pole model, the order of the linear prediction coefficient p is set to be 12, and the penalty parameter $μ = 2$ and number of iteration $Nit = 200$ . The threshold values of the stop criterion for the iterative operation are set to be $ε_{H} = 0.2$ and $ε_{L} = 0.2$ . The set of parameters ${a_{i}, i = 1, 2, \dots, 13}$ in seven iterations are obtained through the Levinson–Durbin algorithm. With Fourier transformation of the parameter of ${a_{i}, i = 1, 2, \dots, 13}$ in each iteration, the power spectrum response curves after quadratic interpolation in seven iterations are shown in Figure 6.

As shown in Figure 6, there are four main centre frequencies and corresponding band widths. The Q-factors corresponding to the centre frequencies of 100 and 230 Hz reflect the main oscillatory behaviours of the simulation signal and are chosen as the Q-factors of the TQWT. Furthermore, in order to increase the overlap between adjacent subbands and preserve the localization of the certain characteristic of the signal, the redundancies of the TQWT in every iteration are all set to be $r_{1} = 3$ and $r_{2} = 5$ , respectively. The decomposition levels are set to be large enough that the frequency response of the wavelets can cover the main frequency components of 100 and 230 Hz. Moreover, it should not be more than the maximum decomposition level. As to the parameter of MCA and SALSA, $α$ and $β$ are set to be different values with consideration of the tuning energy between different components. The corresponding parameters including Q-factors, redundancy, decomposition level, $α$ and $β$ are listed in Table 3.

Table 3.

Parameters of six iterations.

i	$Q_{H, i}$	$Q_{L, i}$	$r_{H, i}$	$r_{L, i}$	$J_{H, i}$	$J_{L, i}$	$α_{i}$	$β_{i}$
1	39.92	7.63	3	5	106	17	0.62	0.38
2	9.31	0.82	3	5	5	21	0.3	0.54
3	4.86	0.80	3	5	5	11	0.3	0.55
4	3.82	0.81	3	5	5	9	0.3	0.56
5	3.69	0.81	3	5	5	9	0.33	0.55
6	3.24	0.81	3	5	5	8	0.33	0.54
7	3.12	0.80	3	5	5	8	0.28	0.57

As shown in Table 3, the value of $Q_{H, i}$ decreases continuously and there is a plunge between the first iteration and the second iteration, and likewise $Q_{L, i}$ . This is because the amplitude-modulated components corresponding to high-oscillatory components have been separated in the first iteration and the remaining amplitude-modulated components play a small role in the remaining signal. From the second iteration to the seventh iteration, the periodic transient components become the high-oscillatory components and are separated iteratively, which induce the continuous decrease of ${Q_{H, i}, i = 2, 3, \dots 7}$ , but the pulse signal corresponding to ${Q_{L, i}, i = 2, 3, \dots 7}$ remains and the value stays around 0.8. Based on the definition of Q-factor, ${Q_{L, i}, i = 2, 3, \dots 7}$ is set to be 1.

The correlative coefficients of periodic transient signal between the simulated and decomposed shown in Figure 7 validate the effectiveness of the proposed signal separation method. With the increase of iteration, the correlative coefficients increase continuously and change a little at the sixth and seventh iterations. Figure 7 shows the final decomposed waveform and frequency response of the simulated signal, which show that the periodic transient components have been separated effectively.

However, significant difference has been found between the decomposition results of the simulated signal and the given signal with the proposed method, e.g. the amplitude of decomposed periodic transient component is half of the simulated signal (Figure 8). The main reasons for the difference are that the Daubechies wavelet basis function is symmetrical and double-sided regardless of the value of the chosen Q-factor, but the given signal is single-sided. The linear combination of the double-sided wavelet basis function chosen based on the Q-factor has limited ability to express the single-sided signal. Through the proposed method, the given single-sided signal is redistributed and becomes double-sided. Therefore, the amplitude of the decomposed signal becomes half of the simulated signal.

Thin-walled parts machining experiments

The proposed method was adopted for machining chatter feature extraction of thin-walled parts machining. Furthermore, the decomposition results and transient feature extraction ability of the proposed method is compared with that of the wavelet packer transform (WPT) to demonstrate further the effectiveness of the proposed method.

The experiment concerns with the machining chatter of thin-walled parts. The machining process was carried on a thin wall open geometry of a cantilever of a 6-mm thick aluminium workpiece overhanging for a length of 41 mm using a Φ10-mm, four-flute milling cutter, as shown in Figure 9, on a Mikron HSM 500 vertical machining centre with a 4.2-kW, 42000 rpm spindle. Experimental trials have been carried out with axial (a_p=0.5 mm) and radial (a_e=5 mm) depth of cut while maintaining a constant cutting speed (n=4000 rpm) and feed rate of 400 mm/min. The bottom middle side of the workpiece was mounted on a jaw vice and the jaw vice is mounted on a dynamometer (Kistler 9253B23) so that the cutting force is aligned along the X-direction of the dynamometer. Furthermore, mounting an eddy current displacement sensor (Keyence EX-V10) on the exactly opposite side of the thin-walled parts with a magnetic base enables the displacement signal and cutting force signal to be compared and studied simultaneously. The synchronous acquisition of the displacement and force signal proceeds through signal conditioner and the data acquisition system, a National Instruments PCI-MIO-16E at 3000 Hz sampling rate with in-house software DEWESoft 6.

The measured displacement signal and cutting force signal with length of 500 000 are shown in Figure 10 and divided into three parts based on the cutting location. As shown in Figure 10, chatter happens at the beginning and end of the machining process. Undergoing a process from a stable state to a chattering state, the displacement signal of III becomes large and exhibits more oscillatory behaviours compared with the signal in II. The cutting force signal is reduced at the end of the machining process due to the deformation of the workpiece and becomes more oscillatory because of the chatter. Figure 11 shows the part displacement and force signal corresponding to II and III in Figure 10, which indicates that the chatter signal in III includes more pulses and exhibit more oscillatory behaviours than the signal in II. The signal spectrums of the part displacement and force signal are shown in Figure 12. In spite of the obvious difference of the waveform between the stable signal and chatter signal, it is difficult to demonstrate specific characteristic features from the time domain waveform of the displacement and force signals. Moreover, the frequency spectrums of displacement and force signal take in the same dominant frequencies between the stable signals and chatter signals. Both in the time waveform and frequency spectrum, there is no direct and completely useful information about the machining chatter.

The change in the spectrum for signal in the stable state and chattering state gives an insight of the chatter characteristics separation. The spectrum in Figures 12(a) and (b) show that the stable signal are distributed over a wide frequency range with a relative low amplitude, but the spectrum of chatter signal in Figures 12(c) and (d) demonstrate narrow frequency band characteristics and the signal energy ranges mainly around the cutting tool frequency (266 Hz for displacement signal) or the frequency doubled of the cutting tool frequency (532 Hz for force signal). Furthermore, the appearance of component corresponding to frequency of 254 Hz becomes the significant difference between the stable state and chattering state. This change in the spectrum results from the change in the oscillatory properties of the signal and can be quantified with Q-factor. Considering the more oscillatory behaviours of the chatter signals, the proposed method for signal separation can be used to analyse the measured signal and extract chatter characteristics.

The measured vibration signals in Figure 11(c) are used to evaluate the effectiveness of the proposed signal separation method. Similarly, for the simulation signal, the order of the signal model p is set to be 200 to ensure the established model of signal be an all-pole model. The parameter of ${a_{i}, i = 1, 2, \dots, 201}$ is obtained through the Levinson–Durbin algorithm. The penalty parameter $μ$ of the SALSA is set to be 2 and number of iterations for SALSA is 200. The corresponding parameters including Q-factors, redundancy, decomposition level, $α$ and $β$ are listed in Table 4. According to method shown in Figure 4, the power spectrum curves of four iterations are shown in Figure 13.

Table 4.

Parameters of four iterations.

i	$Q_{H, i}$	$Q_{L, i}$	$r_{H, i}$	$r_{L, i}$	$J_{H, i}$	$J_{L, i}$	$α_{i}$	$β_{i}$
1	88.78	186.54	3	3	264	486	0.1	0.1
2	103.56	121.06	3	3	308	316	0.1	0.1
3	102.01	65.82	3	3	303	173	0.1	0.1
4	90.86	–	3	3	270	–	0.1	0.1

As shown in Figure 13, the Q-factors corresponding to the centre frequencies of 266 and 254 Hz can reflect the main oscillatory behaviours and are chosen as the Q-factors of the TQWT. A group of ${Q_{L, i}, i = 1, 2, 3, 4}$ decreasing gradually with the iterative calculation can be obtained. The frequency 266 Hz disappeared in the fourth iteration. As a result, the cutting vibration information corresponding to 266 Hz was extracted in the third iteration.

The final decomposition results, including machining chatter component, cutter rotation component and residual component, are shown from Figures 14(a)–(c) and the corresponding frequency spectrums are shown in Figures 15(a)–(c). The decomposition results show that the cutting component and the chatter component have been separated effectively.

For the purpose of comparison, the same signal is processed by WPT to demonstrate the effectiveness of the proposed method. The decomposition level is set to be 4 and the sym8 wavelet of Sym wavelet series is chosen as the basis function. The main separated components and corresponding frequency spectrums are shown in Figures 16 and 17. As shown in Figure 17(a), the low-band frequency component has been separated effectively, but the components corresponding to frequency of 254 and 266 Hz shown in Figure 17(b) have not been separated effectively and they are still mixed up. WPT fails to extract the machining chatter component from the raw signal.

Conclusions

In this article, decomposition of the full-oscillatory components from noisy machining vibration signals by minimizing the Q-factor variation has been developed. A precise calculation method of Q-factors based on quadratic interpolation of LPC was proposed in this paper. An iterative algorithm by combining the sparsity-enabled signal analysis method and precise calculation of Q-factor was designed to decompose full-oscillatory components from noisy machining vibration signals. Numerical simulation and practical machining vibration procedure was performed to validate the effectiveness of the proposed method. Comparative research between the proposed approach and WPT was performed. The results indicated that the proposed approach exhibited a more effective extraction ability.

The component directly related to the chatter has been separated from the mixed signal. With further analysis of the separated component, more precise information such as the centre frequency, band width, spectral kurtosis, crest factor and so on will be obtained, which will deepen the understanding of chatter in thin-walled processing. Secondly, the mechanical part, which induces the chatter, will be identified and the further action will be proceeded to avoid the chatter. All this work will be performed in the future study.

Footnotes

Funding

This work is supported by National Basic Research Program Funding Agency of China (Grant No. 2014CB046604) and the National Natural Science Funding Agency of China (Grant No. 51305062).

Conflict of Interests

The authors declare that there is no conflict of interest.

References

Bin

Gao

et al . (2012) Early fault diagnosis of rotating machinery based on wavelet packets – empirical mode decomposition feature extraction and neural network. Mechanical Systems and Signal Processing 27: 696–711.

Cai

Chen

(2013) Sparsity-enabled signal decomposition using tunable Q-factor wavelet transform for fault feature extraction of gearbox. Mechanical Systems and Signal Processing 41: 34–53.

Cao

Zhou

Chen

(2015) Chatter identification in end milling process based on EEMD and nonlinear dimensionless indicators. International Journal of Machine Tools & Manufacture 92: 52–59.

Davies

Balachandran

(2000) Impact dynamics in milling of thin wall structures. Nonlinear Dynamics 22: 375–392.

Guo

et al . (2014) Envelope extraction based dimension reduction for independent component analysis in fault diagnosis of rolling element bearing. Journal of Sound and Vibration 333: 2983–2994.

Han

Guo

et al . (2015) Feature extraction method of bearing AE signal based on improved FAST-ICA and wavelet packet energy. Mechanical Systems and Signal Processing 62–63: 91–99.

Chen

et al . (2013) Tunable Q-factor wavelet transform denoising with neighboring coefficients and its application to rotating machinery fault diagnosis. Science China Technological Sciences 56: 1956–1965.

Hoell

Omenzetter

(2016) Optimal selection of autoregressive model coefficients for early damage detectability with an application to wind turbine blades. Mechanical Systems and Signal Processing 70–71: 557–577.

Klein

Ingman

Braun

(2001) Non-stationary signals: phase-energy approach – theory and simulations. Mechanical Systems and Signal Processing 15(6): 1061–1089.

10.

Kolluru

Axinte

(2013) Coupled interaction of dynamic responses of tool and workpiece in thin wall milling. Journal of Materials Processing Technology 213: 1565–1574.

11.

Lei

Lin

et al . (2013) A review on empirical mode decomposition in fault diagnosis of rotating machinery. Mechanical Systems and Signal Processing 35(1–2): 108–126.

12.

Luo

Liang

(2013) A kurtosis-guided adaptive demodulation technique for bearing fault detection based on tunable-Q wavelet transform. Measurement Science and Technology 24: 055009.

13.

Peng

Chu

(2004) Application of the wavelet transform in machine condition monitoring and fault diagnostics: a review with bibliography. Mechanical Systems and Signal Processing 18: 199–221.

14.

Peng

Tseb

Chu

(2005) An improved Hilbert–Huang transform and its application in vibration signal analysis. Journal of Sound and Vibration 286: 187–205.

15.

Quintana

Ciurana

(2011) Chatter in machining processes: a review. International Journal of Machine Tools & Manufacture 51: 363–376.

16.

Selesnick

(2011a) Wavelet transform with tunable Q-factor. IEEE Transactions on Signal Processing 59(8): 3560–3575.

17.

Selesnick

(2011b) Resonance-based signal decomposition: a new sparsity-enabled signal analysis method. Signal Process 91(12): 2793–2809.

18.

Shao

Shi

(2011) Power signal separation in milling process based on wavelet transform and independent component analysis. International Journal of Machine Tools & Manufacture 51: 701–710.

19.

Wang

Chen

Dong

(2014) Feature extraction of rolling bearing’s early weak fault based on EEMD and tunable Q-factor wavelet transform. Mechanical Systems and Signal Processing 48: 103–119.

20.

Zhang

Fan

et al . (2014) Application of wavelet transform with tunable Q-factor to analysis of non-stationary harmonics. Applied Mechanics and Materials 526: 182–186.