Instantaneous Frequency Identification Using Adaptive Linear Chirplet Transform and Matching Pursuit

An instantaneous frequency identification method of vibration signal based on linear chirplet transform and Wigner-Ville distribution is presented. This method has an obvious advantage in identifying closely spaced and time-varying frequencies. The matching pursuit algorithm is employed to select optimal chirplets, and a modified version of chirplet transform is presented to estimate nonlinear varying frequencies. Because of the high time resolution, the modified chirplet transform is superior to the original method. The proposed method is applied to time-varying systems with both linear and nonlinear varying stiffness and systems with closely spaced modes. A wavelet-based identification method is simulated to compare with the proposed method, with the comparison results showing that the chirplet-based method is effective and accurate in identifying both time-varying and closely spaced frequencies. A bat echolocation signal is used to verify the effectiveness of themodified chirplet transform.The result shows that it will significantly increase the accuracy of nonlinear frequency trajectory identification.


Introduction
Natural frequency is a crucial parameter in dynamic systems, and most of the frequency identification methods are based on the Fourier transform.However, these methods are effective only when the frequency contents of the vibration signals are time-invariant.Linear time-varying (LTV) systems are often used in engineering, such as large flexible structures for outer space exploration, spacecraft, launch vehicles, and vehicle-bridge systems [1,2].Classical linear time-invariant models are not capable of describing the nonlinear and/or time-varying characteristics of such systems.Parameter identification methods for LTV systems have been widely studied in the last two decades, but most of them are based on inaccurate and computationally sectional time invariance [3,4].
Because responses of LTV systems are nonstationary, time-frequency analysis becomes a necessary and unique approach for identification of time-varying parameters of LTV systems.Shi et al. [5] utilized the Hilbert-Huang transform [6] to identify time-varying system parameters by using free and forced responses and assuming intrinsic mode functions from the empirical mode decomposition (EMD) to be orthogonal to each other.Unfortunately, EMD is not capable of decomposing a signal having multiple components of close frequencies.Chen and Wang [7] proposed a signal decomposition theorem called analytical mode decomposition to address the challenges in closely spaced modes.
Due to its multiresolution analysis ability, wavelet transform has been widely studied for linear and nonlinear system identification.Ruzzene et al. [8] demonstrated the extraction of modal parameters from the continuous wavelet transform (CWT) of a dynamical system's impulse response.A comparative study of the Morlet wavelet, Cauchy wavelet, and harmonic wavelet for system identification was performed by Le and Argoul [9], and a criterion for choosing mother wavelets was also introduced, with this technique being extended to the response from ambient excitation tests [10].Based on the orthogonality of wavelet scaling functions, Ghanem and Romeo [11] proposed a method to transform 2 Shock and Vibration the differential equation of an LTV system into simple linear equations and hence the parameter identification problem became a simple problem of solving linear equations.Xu et al. [12] proposed the use of state-space equations instead of the original high-order differential equation in order to decrease calculation costs.Recently an exploratory study of wavelet-based frequency response functions for timevariant systems was presented by Staszewski and Wallace [13].It is a new way to describe time-variant systems in the time-frequency domain using CWT.However, according to the Heisenberg uncertainty principle, wavelet transform is incapable of simultaneously providing high resolutions in both time and frequency domains.For this reason, great difficulties are encountered in using wavelet-based methods for identification of dynamical systems having close and/or fastvarying modal frequencies.Unfortunately, modal frequencies of a highly flexible structure are often very close to each other and they vary with time and the structure's deformed geometry.
Wigner-Ville distribution (WVD) is a quadratic timefrequency representation of a dynamical signal, and its time-frequency resolution is higher than that of CWT [14].However, cross-terms occur in WVD if the processed signal contains multiple components.Many studies have been done to reduce these cross-terms [15].Chen et al. [16] presented a method called chirplet Wigner-Ville distribution that was capable of representing nonlinear and nonstationary signals effectively.The matching pursuit (MP) algorithm [17] was employed in the chirplet signal decomposition, and the vibration signal is decomposed adaptively into a series of linear chirplets.The method has been applied in fault diagnosis [18] and earthquake signal analysis [19].A frequency banding operator was introduced by Angrisani and D' Arco [20] in order to represent the nonlinear frequency trajectories.Wang and Jiang [21] introduced a parameter called curvature to extend the original chirplet atom to cubic form, and results show that the accuracy of time-frequency representation is improved.Yin et al. [22] presented a fast refinement algorithm in matching pursuit, which not only reduced the computation cost but also improved the accuracy of adaptive chirplet decomposition.
In this paper, an instantaneous frequency identification method based on chirplet transform (CT) and MP is proposed.The effectiveness of the method in identifying close modes and time-varying systems is investigated in simulated mass-stiffness-damping systems.A modified CT is presented to reduce the errors occurring in identifying nonlinear varying frequencies.An echolocation chirp signal is utilized to assess the method for practical signals.The results are compared with the CWT-based method.For completeness, the basic theory of CWT and CT is introduced in Section 2; full identification procedures of CWT and CT are established in Section 3; Section 4 presents the application of the methods to simulated LTV and close-mode systems and practical signals.Conclusions are drawn in Section 5.

Time-Frequency Representation of
Free Vibration Signal where Φ() is the Fourier transform of ().This kind of function is fast attenuated and supports a short duration in the time domain.A whole family of these functions can be generated by shifts in the time domain and scaling in the scale domain: where  is the scale parameter which is related to frequency and  is the translation parameter related to time.
The local characteristic and time-frequency representation of a signal () can be illustrated by decomposing () into wavelet coefficients using the basis of wavelet functions  , .This can be expressed by the following equations: where  , () is the complex conjugate of  , and (, ) are the wavelet coefficients.The resolution of the wavelet transform in time and frequency domains is defined by the following equations: where Δ  and Δ  are the duration and bandwidth of the wavelet function, respectively.The work presented in this paper utilizes the Morlet wavelet, defined as The Morlet wavelet is a complex wavelet as shown in Figure 1.Each wavelet basis supports a time-frequency window restricted by It is noted that an increase of time resolution will lead to a decrease of frequency resolution, which meets the requirements of the Heisenberg uncertainty principle.
Most of the real signals are asymptotic and can be given in the form of the sum of single components:  where   () are the amplitudes and   () are the phases.It is assumed that, in an asymptotic signal, the variation of phase is much faster than the variation of amplitude.The analytical signal of a single component can be given by The ridge and skeleton of the CWT contain most of the information of the original signal, with the ridge being defined as where  is a constant depending on the mother wavelet.
Several algorithms can be found to extract the ridge.In this paper, the local maximum algorithm is performed.Once the ridge is obtained, the instantaneous angle frequency () can be obtained by where  0 is the center angle frequency of the wavelet.

WVD and Linear Chirplet Decomposition.
For an arbitrarily signal (), the WVD is defined as WVD has the highest time-frequency resolution, and it is an ideal approach to represent mono component nonstationary signals.As shown in (10), the WVD is not linear, and the WVD of a sum of multiple signal components is not equal to the sum of the WVD of these signal components.For example, the WVD of signal where is the cross-term.Due to the cross-term interference issue, researchers seldom use it directly to identify time-frequency structures of signals; however, adaptive linear chirplet decomposition combined with WVD is a promising approach to avoid cross-term and obtain high resolution time-frequency representation.
The linear chirplet transform is defined as the inner product of the signal () and the chirplet basis function () [23]: in this paper, Gauss chirplet is used and its analytical expression is The Gauss chirplet is generated by time shifting, frequency shifting, scaling, and rotating of the original Gauss function.A typical Gauss mother chirplet and its WVD are shown in Figure 2.

Matching Pursuit.
Since the Gauss chirplets do not form an orthogonal basis, matching pursuit is a promising solution to realize adaptive decomposition.The vital idea of the method is to select optimal chirplets () to approach the original signal () and decompose () step by step using the MP algorithm.() can be written as where the residual signals are In the first step,  = 0,  0 () = (), and the weighting factor   is the canonical inner product of the residual signal   and chirplet function (); that is, The object of MP is to find optimal parameters of chirplets and generate () adaptively to minimize the energy of the residual signal: This is equivalent to the following equation:  This optimization problem is solved using the method presented by O'Neill and Flandrin [24].In order to identify the instantaneous frequency, the chirplet spectrum is utilized: where AD represents adaptive decomposition and WVD   is the WVD of the optimal chirplets (): The cross-term does not exist in the chirplet spectrogram, and the resolution of time-frequency depends on the durations of the optimal chirplets.Once the chirplet spectrum is obtained, a local maximum algorithm [25] similar to the ridge extraction method is utilized to estimate the instantaneous frequency.

Continuous Wavelet Transform-Based Method.
In the continuous wavelet transform-based (CWT) method, the vibration signal is decomposed using wavelet functions.The resolution is dependent on the time and frequency windows of each wavelet.Instantaneous frequencies can be estimated in different sizes of frequency windows by controlling translation and dilatation.However, as mentioned before, an increase of time resolution will lead to a decrease of frequency resolution, and hence different resolutions will be used in different modes.In the time-varying systems, the frequency resolution varies with time.The identification procedure using CWT is as follows: (1) Calculate the CWT of the vibration signal using the Morlet wavelet.
(2) Extract the ridges of the CWT using local maximum algorithm.

Linear Chirplet Transform-Based Method.
In the CTbased method, the cross-term of the classical WVD is eliminated.The resolution is dependent on the time-frequency property of optimal chirplets.The MP algorithm ensures the sparsity of the signal decomposition.The noise is departed in the residual signal while most of the time-frequency characteristics of the original response are retained in the selected chirplets.It is also noticed that CWT is a special case of CT where only time shifting and scaling parameters are available and other parameters are equal to zero.The identification procedure using CT is as follows: (1) Obtain the parameters of optimal chirplets via MP algorithm.
(2) Calculate the WVD of each optimal chirplet basis.The chirplet spectrum is obtained by the linear weighted sum of the WVDs.
(3) Estimate the instantaneous frequencies using local maximum algorithm.

A Modified Version of CT.
Since the basis functions of CT are linear chirplet functions, the method is effective in identifying systems with linear varying frequencies; however, in case of nonlinear varying frequencies, the nonlinear trajectories are approached by several linear segments, and thus significant errors occur.In order to solve this problem, a  modified version of CT is utilized.In the matching pursuit process, the key step is to solve the optimization problem presented in (19).In the original CT method, the constraints of parameters in the optimal chirplet   () are not strict.In order to obtain higher time resolution, the duration of the optimal chirplet is restricted in the modified CT, and hence the number of optimal chirplets becomes larger and the time resolution is improved.
In the modified CT, the duration of optimal chirplets is restricted to 1/ of the signal length.Here  is defined as the restriction parameter.The time support of selected optimal chirplets depends on .A large value for  will lead to high time resolution and vice versa.Thus, the optimization process only chooses the chirplets with short time support.The vibration signal is decomposed using selected short chirplets.
The procedure of the modified CT based method is shown in Figure 3; the difference between the modified version and the classic one is only in the first step matching pursuit where the parameter  is introduced.With this parameter, chirplets with wide time duration will not be considered as optimal chirplets to represent the vibration signals; in other words, a constraint condition is added in the original optimization problem.This is significant in analyzing vibration signals with nonlinear varying frequencies.
It should also be noted that the bandwidth of each chirplet gets wider because of the restriction of duration, and the frequency resolution is reduced.Nevertheless, it is still much better than that of the CWT-based method.The effectiveness of the modified CT is verified in the following section.

Numerical and Practical Signal Analysis
In this section, several examples of time-varying systems and systems with closely spaced modes are investigated to illustrate the effectiveness and accuracy of the algorithms presented in Figure 3.A two-degree-of-freedom system with closely spaced modes and two single-degree-of-freedom time-varying systems are studied.The identification procedures are performed using the free decay responses of each system.The responses are assumed to be contained with zero mean Gaussian noises and the signal-to-noise ratio is 10 dB.The sampling frequency is 10 Hz and the number of samples is equal to 1024.A practical signal is also used to verify the effectiveness of the modified CT-based method.

Linear Time-Varying
Frequency.The single-degree-offreedom linear time-varying system is governed by the general equation: where (), (), and () are time-varying mass, damping, and stiffness, respectively. is assumed equal to 0 since free vibration is considered.In this case, () and () are assumed constant, and () varies with time linearly:  = 1,  = 0.02, () = 300 + 5.The free decay response is calculated by the 4th Runge-Kutta method and is shown in Figure 4. Gaussian white noise is added to the response.It can be seen from Figure 5 that the system property can hardly be identified from the time domain response because of the effect of noise.
The FFT spectrum of the response is shown in Figure 6.The signal is analyzed by both CWT-and CT-based methods.The CWT scalogram and the CT spectrogram are shown in Figures 7 and 8, respectively.
From the comparison, it can be seen that both of the methods present the true natural frequency of the system.However, a difference is clear: the frequency resolution of CT is much higher than that of CWT.In the CWT scalogram, higher frequency has lower resolution, while the CT spectrogram has uniform resolution for all frequencies.The identification results obtained by local maximum algorithm of both CWT and CT are shown in Figure 9.It is noted that the effect of noise is obvious in the CWT scalogram; however, this phenomenon is not observed in the CT spectrogram.This is because the noise is automatically ignored in the matching pursuit process and only a few  The CWT scalogram and the CT spectrogram are shown in Figures 10 and 11, respectively.As can be seen in Figure 11, several linear segments are used to approach the sinusoidal wave.The variation of frequency is not well tracked.Although the frequency resolution of CT is much higher than CWT, CWT is more suitable for tracking a nonlinear trajectory (see Figure 10).Similar to the previous case, the results also show that CWT is more sensitive to noise than CT.In this study, a restriction of the duration of chirplets is introduced to overcome the problem.The duration of each chirplet in the time domain is restricted to 5% of the original signal length or less in the iteration process.As a result, more optimal chirplets are required; nevertheless, the decomposition of the signal is still sparse.The improved result is shown in Figure 12, where it can be seen that the frequency trajectory is well tracked, and the high frequency resolution and good antinoise performance are retained in the improved method.The frequency trajectories extracted from both the CWT scalogram and the modified CT spectrogram are shown in Figure 13, which shows that the modified CT is effective and accurate.
In order to evaluate the accuracy of the modified CT, the residual sum of squares (RSS) is introduced as follows: where  is the exact value of the parameter,   is the identified value,  is the total sample number, and  is the sample  number for each time instant.In order to obtain accurate results, small RSS is needed.
Table 1 shows the RSS of the identification results obtained by CWT, CT, and modified CT.It can be seen that the CT-based method is more accurate than the CWT-based method, and the modified CT is more accurate than the original CT.Thus, the modified CT is a necessary approach when high accuracy is required.
In order to find the optimal value of the restriction parameters , a series of simulations is performed.The RSS of each identification result with different  is shown in Figure 14.It can be seen from Figure 14 that, with an increase of , the RSS decrease as expected.After  reaches 50, the decrease rate becomes relatively small.Considering that the increase of  will increase the computation cost, 50 is the optimal value of  in this case.It is also worth noting that the RSS will increase after  exceeds 70.This is because the frequency resolution becomes lower, and the influence of noise increases.
Both CWT-and CT-based methods are performed to estimate the natural frequencies of the two modes.Figure 16 shows the CWT scalogram and Figure 17 shows the timefrequency distribution of CT.From Figure 17, two closely spaced modes are completely decoupled, and the correct frequencies are identified.However, in Figure 16, two modes are separated only in the very first period; most of the time, they are overlapped.Thus, the CT-based method is more powerful in separating closely spaced modes than the CWTbased method.It is also worth noting that only four linear chirplet basis functions are sufficient to represent the original signal.In other words, only four iterations are required in  the matching pursuit, which means the CT-based method can also reduce computational complexity.

Practical Signal Analysis.
In order to verify the effectiveness of the proposed method in the application to real data, the Eptesicus fuscus bat echolocation signal [26] is used.This signal is a typical multicomponent and nonstationary signal.It contains three different modes which are overlapped in both time and frequency domains.practical signals contain nonlinear varying frequencies, so the modified CT is a necessary approach to identify frequency trajectories in real systems accurately.
4.5.Discussions.Time-frequency resolution is a key factor in instantaneous frequency identification.As mentioned in Section 2, CWT is not able to obtain high time and frequency resolution simultaneously, and its resolution varies with scale.This disadvantage is obvious in Section 4.3 where closely spaced time-varying frequency is considered, due to the variation of resolution, CWT-based method is not able to track two frequency trajectories accurately.On the other hand, the identification accuracy of CWT can be easily affected by noise.Since the CT method presented in this paper is based on matching pursuit, it is less sensitive to noise.Each selected chirplet possesses independent timefrequency support and their WVDs represent special areas in spectrogram with highest time-frequency resolution.The accuracy of CT and WVD based method mostly depends on the time-frequency support of selected optimal chirplets.It is verified in this section that the accuracy is improved significantly by using modified CT.
In cases of nonlinear varying frequencies, computation cost remains a challenge.Due to the restriction in time support, the number of optimal chirplets increases in the proposed modified approach; in other words, more iterations are required in MP and the computation cost is significantly increased.Introducing a nonlinear basis chirplet is also a promising approach to improve accuracy; see reference  [ 20,21]; in these methods, only one pursuit is needed for each frequency component; this is similar to "curve fitting" while the proposed method in this paper is similar to "piecewise linear approximation" of nonlinear frequency trajectories.In practical application, the balance between accuracy and efficiency needs to be considered.

Conclusions
Instantaneous frequency identification based on CT and WVD is proposed.The adaptive chirplet transform using MP is employed to represent the vibration signal in the timefrequency plane.Complete frequency identification procedure is established, and local maximum algorithm is utilized to extract the instantaneous frequency from the CT spectrogram.
A system with closely spaced modes and two kinds of time-varying systems are investigated to validate the proposed method.The results are compared with an instantaneous frequency identification method based on CWT.Results show that both the CWT-and CT-based methods present the true frequencies.However, the frequency resolution of the CT-based method is much higher than that of the CWT-based method; CT clearly separated the close modes while CWT failed.The antinoise performances of the two methods are compared by adding a noise of 10 dB signalto-noise ratio.Results show that the CT-based identification method has an excellent antinoise performance.
A modified version of CT is proposed to identify nonlinear varying frequencies, and its effectiveness is verified using both simulated and practical signals.Results show that the modified CT successfully improves the accuracy in identifying nonlinear varying frequencies.

4. 2 .
Nonlinear Time-Varying Frequency.The methods are applied to another kind of time-varying system: the parameters vary nonlinearly.In this case, () varies with time nonlinearly:  = 1,  = 0.02, () = 300 + 200 sin 0.01.The frequency of vibration is modulated by a sinusoidal wave.

Figure 14 :
Figure 14: RSS with the increase of .

Figure 18 Figure 17 :
Figure 17: CT spectrogram of displacement response of case 3.

Table 1 :
RSS of different identification results of case 2.