Transient Feature Extraction method based on adaptive TQWT sparse optimization

: Aiming at the problem of strong impact, short response period and wide resonance frequency bandwidth of transient vibration signals, a transient feature extraction method based on adaptive Tunable Q-factor Wavelet Transform (TQWT) was proposed. Firstly, the characteristic frequency band of the vibration signal was selected according to the time-frequency distribution. Based on the characteristic frequency band, the sub-band average energy weighted wavelet Shannon entropy was used to optimize the number of decomposition layers, quality factor and redundancy of TQWT, so as to achieve the adaptive optimal matching of the impact characteristic components in the vibration signal. Then, according to the characteristics of the transient impact of the telemetry vibration signal, the TQWT decomposition coefficients were sparse reconstructed to obtain more sparse impact characteristics, and the weighted power spectrum kurtosis was used as the impact characteristic index to select the optimal sub-band, Finally, the inverse transform of TQWT was used to reconstruct the optimal sub-band to enhance its weak impact features. The simulation and measured signal processing results verify the effectiveness of the algorithm. A transient feature extraction method based on adaptive TQWT sparse optimization is proposed to extract the impact characteristics of transient vibration signals. Based on the time-frequency distribution of vibration signal, an adaptive selection strategy for key parameters of TQWT is designed. In addition, the optimization idea was introduced into the TQWT decomposition process to further highlight the more sparse impact components in the signal according to the feature frequency band, which effectively improved the ability to extract the transient features of the signal. The results of simulation experiments show that the new method is superior to the WPT decomposition and EEMD method in terms of noise removal and transient feature extraction.


Introduction
The telemetry vibration signals are time series that include the system operating state, which are collected by sensors such as vibration acceleration or displacement, temperature, and pressure installed in the aircraft. Affected by the vibration of the aircraft itself, the flight environment, the electromagnetic environment and the transmission conditions, the vibrations from the various structural sections will affect, modulate and superimpose each other. The transmission path of the vibration is complex and variable, and the collected telemetry vibration signals are often mixed with a large amount of different frequencies (high frequency and low frequency), impact noise and harmonic components. Therefore, the spectral components are extremely complex and highly coupled, and show strong nonlinearity, non-stationary and transient impact [1].
The anomaly of aircraft system state is usually not obvious, and the signal level is very weak and sparse, under the influence of test environment, noise and system vibration complexity, the impact characteristics are very weak in the transient vibration signal, which is difficult to extract. At present, the research of feature enhancement focuses on linear filtering, adaptive decomposition combined with resonance demodulation and time-frequency analysis [2]. For example, Combet et al. [3] proposed a filtering method based on Spectral Kurtosis (SK) to realize fault feature extraction for gearbox fault diagnosis.
The idea is to find the appropriate center frequency and bandwidth with SK as the index, so as to retain the impact signal and remove the noise and interference signal. However, SK has a high sensitivity to noise (especially in-band noise), which seriously affects the accuracy of resonance band selection. In addition, for non-stationary signals, when there is interference signal with the same frequency as the fault signal component in the signal, SK value will also be large and the in-band interference is difficult to remove. Wang et al. [4] realized transient feature extraction of weak faults by modeling transient features and identifying transient signal model parameters with correlation filtering. However, this modeling method requires a lot of prior information, so it is difficult to establish an accurate transient signal model in practice. Georgoulas [5] et al. and He et al. [6] proposed hybrid intelligent fault diagnosis models respectively for rolling bearings based on EMD and EEMD resonance demodulation. The adaptive decomposition method was used to decompose the non-stationary signal into several Intrinsic Mode Functions (IMF) of instantaneous frequencies with physical significance, and a specific IMF was selected for resonance demodulation to extract the weak fault features, but the inherent problems of the adaptive decomposition method, such as mode aliasing, endpoint effect and pseudo-component, will restrict the effect of signal decomposition and affect the accuracy of feature extraction.
As a representative multi-scale analysis technique, Wavelet analysis can detect the signal feature information hidden in a specific scale by matching the basis function of the inner product transformation principle with fault feature similarity, For example, Fyfe et al. [7] took Morlet wavelet as the wavelet basis, proposed a wavelet threshold de-noising method based on maximum likelihood estimation, and realized the de-noising and fault feature extraction of bearing monitoring signals. Lin et al. [8] proposed a wavelet coefficient screening method for vibration signal analysis based on data analysis, and then completed the reconstruction and resonance demodulation of the optimized wavelet coefficient, and extracted the transient fault information of vibration signal. Nguyen et al. [9] proposed a feature extraction and enhancement method for bearing state monitoring signal based on improved wavelet soft threshold and Hilbert envelope demodulation analysis. The effectiveness of wavelet transform lies in that the wavelet base that matches the fault feature waveform best can best represent the fault feature information under multi-source noise interference, but the wavelet transform has its limitation due to its constant quality factor (which is defined as the ratio of center frequency to bandwidth). Once the basis function and decomposition layer number are determined, the quality factor is fixed. When the central frequencies of each component in the signal to be decomposed differ greatly, good decomposition effect will be achieved. However, when the center frequency of impulse signal is close to that of other signal components, good decomposition effect cannot be obtained. In view of the shortcoming of constant quality factor in traditional Wavelet Transform, the Tunable enactor Wavelet Transform (TQWT) proposed by Selesnick [10], by adjusting the quality factor wavelet, the optimal matching is realized for the characteristic signal components with specific oscillation behavior. Even if the central frequency of the signal component is close to each other, the signal can be separated effectively as long as it has different bandwidth. For example, Luo et al. [11] used kurtosis indicator to select fault feature frequency band, and used the reconstruction signal of optimal TQWT decomposition coefficient and demodulation analysis to realize the extraction of bearing fault features. Li et al. [12] took the normalized weighted value of kurtosis and smoothness index as the objective function and proposed the fault diagnosis method of rolling bearing based on the optimal quality factor signal resonance sparse decomposition. He et al. [13] proposed a feature enhancement method of adjacent coefficient bearing monitoring signal based on TQWT, according to the fault characteristics of Hilbert envelope demodulation, the wavelet coefficient of adjustable quality factor was optimized to realize the feature enhancement of weak bearing signals. Song et al. [14] proposed a new automated epilepsy diagnosis scheme based on TQWT, which can extract various spectral features and realize the intelligent fault diagnosis for roller bearings.
Although the above application has achieved a good practical effect, how to select the key parameters of TQWT to achieve the optimal matching of the characteristic signal components with specific oscillation behavior needs further research.
In view of the above problems, this paper proposes a transient feature extraction method based on adaptive TQWT, which takes time-frequency distribution as the basis, and selects the characteristic frequency band to constrain the number of decomposition layers. Then the sub-band average energy weighted wavelet Shannon entropy is used as the objective function to optimize the quality factor and redundancy. The wavelet basis function matching with the feature waveform is constructed adaptively to improve the feature extraction ability of TQWT. Finally, the weighted power spectrum kurtosis is used as the impact feature index to select the optimal characteristic sub-band, and the optimal sub-band was reconstructed by the inverse transform of TQWT. Hilbert envelope demodulation is used to extract and enhance the weak impact characteristics. Simulation and measured signal processing results show that the method can effectively extract the impact characteristics of transient vibration signals.
The rest of this paper is organized as follows. Section 2 gives the principle of TQWT and presents the major principle and details of the proposed method. Section 3 uses the simulation and measured transient vibration signal to verify the performance of the proposed method as compared to other traditional methods. Finally, Section 4 draws conclusions.

Principle of TQWT
TQWT is a new nonlinear signal decomposition method, which sets the appropriate quality factor wavelet basis function from the frequency domain and establishes the redundant dictionary library, and has the characteristics of perfect reconstruction and moderate completeness, and realizes the decomposition and reconstruction of signals iteratively by means of a double channel band pass filter bank [10]. As shown in Fig.1.
It can be seen that TQWT can be equivalent to the band-pass filter under multi-scale decomposition, and the central Bw of the equivalent filter under each scale are shown as follows Where J is the number of decomposition layers, for the given quality factor Q and redundancy parameter r . In order to ensure that the analysis time domain length of all scales does not exceed the time domain of the analyzed signal, There is a theoretical maximum of J , which is Where N is the signal length, .   represents rounding down. Although TQWT has the advantage of overcoming the traditional wavelet constant quality factor, However, its time-domain and frequency-domain characteristics are uniquely determined by combination of three key parameters{ , , } Q r J . Therefore, in order to achieve the optimal matching of the characteristic signal components with specific oscillation behavior and improve the decomposition performance of TQWT, it is necessary to optimize { , , } Q r J .

Adaptive optimization of TQWT parameters
For the optimization of wavelet parameters, Zhang et al. [15] optimized the waveform parameters, bandwidth and center frequency parameters of the Morlet wavelet basis function by minimizing the Shannon entropy of wavelet to achieve the optimal matching with the impact characteristics of the fault. Taking its ideas for reference, wavelet Shannon entropy can measure the sparsity of wavelet coefficients on all scales. Then, the sparsity of the wavelet coefficients can represent the matching degree between the wavelet basis function and the characteristic waveform. This paper proposes an adaptive optimization strategy for TQWT parameters.
Compared with traditional continuous wavelet transform, TQWT has the following characteristics: a) The sequence length of wavelet coefficients at all scales obtained by continuous wavelet transform is fixed, determined by the number of decomposition layers and the sampling frequency, however, TQWT realizes signal decomposition in an iterative manner by using a double-channel band pass filter, the length of the obtained wavelet coefficient sequence decreases with the increase of the number of decomposition layers as the parameter { , , } Q r J is adjusted.
b) Once the basis function and decomposition layer number are determined, the quality factor of the continuous wavelet transform is constant, constrained by the normalization factor, wavelets of different scales maintain the same energy, but the quality factor of TQWT is adjustable, and the wavelet energy varies with the scale.
Based on the above characteristics, Traditional wavelet Shannon entropy cannot directly measure the sparsity of the TQWT decomposition wavelet coefficients effectively. For this reason, this paper proposes to optimize parameter { , , } Q r J by using wavelet Shannon entropy weighted by sub-band average energy. The definition is as follows: represent the wavelet energy, the length of the wavelet coefficient sequence and the ith element in the jth layer wavelet coefficient sequence of TQWT, respectively. According to Eq. (7) and (8), the sub-band average energy-weighted wavelet Shannon entropy (SAEWSE) comprehensively takes into account the series length unfixed and the energy-scale correlation of TQWT wavelet coefficients, and can measure the sparsity of the wavelet coefficients.
Although there is a theoretical maximum decomposition layers of TQWT, moreover, the increase of decomposition layers mainly affects the decomposition performance density in the low-frequency region. Too many decomposition layers may not only lead to excessive decomposition of feature band information, but also cause redundant decomposition of irrelevant band information and increase the computational burden. Based on the analysis of the time-frequency distribution, the 3dB bandwidth attenuation method is introduced to determine the characteristic frequency band of the signal, so as to constrain the decomposition layer number of TQWT. For signal () n x to be analyzed, the time frequency distribution was obtained by using STFT Where k and The adaptive optimization method of TQWT parameters is summarized as follows:

Sparse optimization of TQWT wavelet coefficients
Although the parameter adaptive optimization method proposed in Chapter B of Section II can extract the characteristic frequency band of the signal and realize the optimal matching of TQWT wavelet to the intrinsic structure (characteristic signal component) of the signal, in essence, TQWT is over-complete redundancy decomposition, and the frequency responses overlap on the adjacent scales. In combination with the characteristics of transient impact of telemetry vibration signals (which meet the sparse decomposition condition), in order to remove these redundant information and improve the energy aggregation of the characteristic frequency band, this section proposes sparse reconstruction of the wavelet coefficients obtained from TQWT decomposition to further highlight its impact characteristics. The objective function of sparse optimization is 12 12 arg min || || S.T. || Where, w is the TQWT wavelet coefficient matrix, ŵ is the sparse expression matrix of the wavelet coefficient, and 1 TQWT  is the inverse transformation of TQWT. Lagrange multiplier method is used to transform equation (12) Where,  is the Lagrange multiplier, and the range of the selected decomposition layers is determined by the upper and lower limits of the characteristic frequency band, that is, the selection range of the wavelet layers is determined by combining equations (11) and (14) () For the optimization problem of equation (12), ADMM algorithm is adopted in this paper for quick solution. Here, the alternating iterative formula is directly given Where s is the auxiliary variable introduced by ADMM algorithm, and  is the penalty parameter.

Optimal selection of TQWT wavelet sub-band
Kurtosis contains time-domain kurtosis and frequency-domain kurtosis. Literatures [16]- [17] show that, as an index to measure the impulse property of the signal, transient impulse components can be effectively detected in the diagnosis of rotating machinery fault. Considering that the TQWT wavelet coefficient has the characteristics of energy-scale correlation, the power spectrum kurtosis is used to measure the overall sharpness of the sub-band signal spectrum, and then to conduct the optimal selection of the TQWT wavelet characteristic sub-band. The definition is as follows: After obtaining the sparse expression matrix ŵ of TQWT wavelet coefficients, the sub-band signal matrix   (18) Where ( , ) j  E P is the mean of ( , ) j  E P .

Steps of the proposed method
The impact characteristics of the vibration signals are very weak and sparse, which is difficult to extract and enhance. In order to solve the problem, the adaptive TQWT is used to decompose the vibration signals into multi-scale, and the specific sub-band are extracted to enhance the impact characteristics of the signals. The algorithm is as follows: a) Pre-processing: The vibration signals need to pre-processing: zero drift correction, trend item removal, and outlier elimination, etc.

b) Adaptive TQWT decomposition:
The method in Chapter 2.2 is adopted to adaptively select the optimal TQWT parameter ( , , ) opt opt sopt Q r J , and the signal () n x is decomposed to obtain the TQWT wavelet coefficient matrix w .
c) Sparse optimization of TQWT coefficients: the sparse expression matrix ŵ of TQWT wavelet coefficient matrix w was obtained by using Eq. (13) -Eq. (15). e) Envelope spectrum analysis：The sub-band signal envelope with the maximum power spectrum kurtosis is selected for spectrum analysis to extract its impact characteristic frequency.

Experiment 1：Simulation signal analysis
A simulation experiment is designed to compare the performance of the EEMD-based feature extraction method, WPT-based feature extraction method and the proposed method. The signal to noise ratio (SNR) of the signal and the joint time-frequency entropy (JTFE) [18]   bandwidth are 1386Hz and 673Hz, respectively. It can be clearly seen that the range of 3dB bandwidth can contain the main information of the pulse, which verifies the effectiveness of the proposed optimal feature band selection method. Fig.2(a) shows the schematic diagram of TQWT decomposition layers selected from the upper and lower limits of the characteristic frequency band. From right to left, the frequency response curve of TQWT equivalent filter from low level to high level is presented. The number of TQWT decomposition layers is determined by the lower frequency limit of the characteristic frequency band. Fig.2(a) shows that when Q=3.1 and r=3.5, the number of decomposition layers can be intuitively determined to be 11. the parameter { , , } Q r J was optimized and selected. As can be seen from Fig.3, when Q=3.1and r=3.5, the value of SAEWSE is the smallest.
The simulation signal is decomposed according to the optimal parameter ( , , ) opt opt sopt Q r J , and the results are shown in Fig.4. Then, the optimal characteristic sub-band is determined by the power spectrum kurtosis, as shown in Fig.5. Sub-band 8 is the selected optimal characteristic sub-band, which has the highest PSK value. The envelope-spectrum analysis is performed and the result is shown in Fig. 6. The impact characteristics with a period of 0.02s could be clearly seen. Moreover, the impact characteristic frequency and its frequency harmonic components of 1-5 can be clearly identified, so it can be seen that the proposed method can extract the impact characteristic under the condition of strong noise.  Fig.7 and Fig.9 show the time domain waveform and envelope spectrum of characteristic sub-band signal extracted by wavelet packet decomposition and EEMD decomposition. The wavelet packet basis function adopts db4 wavelet. According to the frequency range of the characteristic frequency band, the 3-layer wavelet packet is determined to be decomposed into 8 nodes. The single-branch reconstruction of the wavelet coefficients of 8 nodes was carried out successively. The single-branch reconstructed signal with the optimal characteristic sub-band of node (3,2) was determined according to the maximum power spectrum smoothness index, and its envelope spectrum was analyzed, as shown in Fig.8. The EEMD method obtained 9 IMFs, and the same maximum PSK value determined the optimal characteristic sub-band as the IMF2, which was analyzed by envelope spectrum.  More obvious impact characteristics can also be observed from Fig.7. However, the extraction of the impact characteristic frequency and its frequency doubling components is greatly affected by the noise, especially the basic annihilation of the frequency harmonic components 4 and 5 in the noise. The amplitude of other characteristic frequencies is also smaller than that in Fig.6, and there is still an obvious noise interference frequency component. It can be seen in Fig.9, the EEMD method reduces the degree of mode aliasing in the signal decomposition process by adding random noise. However, with the introduction of noise, there are more noise interference components, but the impact characteristic frequency and its frequency doubling components are mixed with the interference frequency, so it is difficult to identify. Table.1 presents the quantitative analysis results of the above three methods, by comparison, it can be seen that the adaptive TQWT sparse optimization method has the ability to extract transient features efficiently and accurately. The output SNR was 3.08db, and the envelope spectrum smooth index JTFE was 0.7816 at the minimum. It shows that the extracted characteristic sub-band energy distribution is more concentrated and has less in-band noise, which indicates that the proposed method has better capability of transient feature extraction with WPT decomposition and EEMD method.

Experiment 2：Verification of measured signals
In order to verify the effectiveness of the proposed method, the telemetry vibration signals collected by high-frequency vibration sensors in the test mission of a certain aircraft were used for processing and verification. The sampling frequency was 5.12KHz. The results are shown in Fig.10. As SNR calculation requires the prior information of noiseless signal, it can be used as a good quantitative index in simulation research, but it cannot be used for the measured signal. Therefore, for the measured and telemetry vibration signals, we only use PSK and JTFE to quantify the time-frequency energy distribution.  On the SSE curve, the driving frequency corresponding to the maximum value is 1.763KHz, and the upper and lower limits of the characteristic frequency band determined by the 3dB bandwidth are respectively 1973Hz and 1592Hz. It is clear that the 3dB bandwidth range can contain the main information of the pulse. According to the frequency response diagram of TQWT band-pass filter in Fig.10(a), the number of decomposition layers can be intuitively determined to be 10. , According to the weighted wavelet Shannon entropy of sub-band average energy, the parameter { , , } Q r J was optimized and selected. As can be seen from Fig.11, when Q=6.2 and r=8, the SAEWSE value is the smallest, thus determining the optimal parameter to decompose the measured signal. The decomposition results of TQWT are shown in Fig.12. According to Fig.10(a), select the number of layers (3 to 10 layers) of the sparse reconstruction of the wavelet coefficient, complete the sparse optimization of the wavelet coefficient, and further highlight the impact characteristics of the signal. Then, single branch reconstruction is carried out on the wavelet coefficients after sparse optimization. The power spectrum smoothness index was used to determine the sub-band 6 as the optimal characteristic sub-band, as shown in Fig.13, and the results were shown in Fig.14 for envelope spectrum analysis.  Similarly, the time domain waveform and envelope spectrum of characteristic sub-band signal are extracted by wavelet packet decomposition and EEMD decomposition in Fig.15 and Fig.16. The wavelet packet basis function also adopts db4 wavelet. According to the frequency range of the characteristic frequency band, the two-layer wavelet packet is determined to be decomposed into four nodes, and the wavelet coefficients of the four nodes are reconstructed by single branch. The single-branch reconstructed signal with the optimal characteristic sub-band of node (2,2) was determined according to the maximum power spectrum smoothness index, and its envelope spectrum was analyzed. The EEMD method obtained 10 IMFs, and the same maximum power spectrum smoothing index determined the optimal characteristic sub-band as the second decomposition component imf2, which was analyzed by envelope spectrum. Table 2 presents the quantitative analysis results of the above three methods, due to the small difference in the center frequency of each impact component in the measured transient vibration signal, the wavelet packet decomposition has its limitations under the influence of constant quality factor. Once the basis function and decomposition layer number are determined, the analysis frequency band of each layer is fixed and invariable, which is likely to span or not include the whole characteristic frequency band, so it is difficult to achieve good decomposition effect. Its JTFE is only 0.8395, and PSK is 13.03. The EEMD method reduces the degree of mode aliasing in the signal decomposition process by adding random noise, but when the modal frequency components are similar, it cannot completely eliminate the mode aliasing, in addition, with the introduction of noise, there are inevitably more interference frequency components, resulting in a larger JTFE of only 0.8742 and PSK of 15.26. The comparison shows that the adaptive TQWT sparse optimization method in this paper has efficient and accurate transient feature extraction capability. Envelope spectrum smooth index PSK is 20.71, JTFE is the minimum of 0.8028, it reflects that the extracted characteristic sub-bands have more concentrated energy distribution and less in-band noise. Therefore, it can be concluded that the proposed method has better extraction capability of transient features with wave packet decomposition and EEMD method.