# Flutter Test Data Processing Based on Improved Hilbert-Huang Transform

Authors: Hua Zheng; Junhao Liu; Shiqiang Duan
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Flutter tests are conducted primarily for the purpose of modal parameter estimation and flutter boundary prediction, the accuracy of which is severely affected by the acquired data quality, structural modal density, and nonstationary conditions. An improved Hilbert-Huang Transform (HHT) algorithm is presented in this paper which mitigates the typical mode mixing effect via modulation. The algorithm is validated by theory, by numerical simulation, and per actual flight flutter test data. The results show that the proposed method could extract the flutter model parameters and predict the flutter speed more accurately, which is feasible for the current flutter test data processing.


---

## Body

## 1. Introduction

“Flutter” is the self-excited vibration of an elastic structure under the coupling of aerodynamic force, elastic force, and inertial force; it is often accompanied by catastrophic structural damage [1]. Flutter analysis is a crucial aspect in the design of new or modified aircraft. Wind tunnel tests, aeroelastic models, and flight flutter tests are important components of flutter analysis. Such tests are risky, time-consumptive, and costly to conduct. To this effect, accurate and effective flutter test data processing techniques are in high demand. One of the most popular flutter boundary prediction (FBP) methods currently is the damping-based method [2, 3], which serves to extrapolate a curve fitted by the estimated structural modal damping factor against the airspeed to the abscissa while zero damping is considered the critical flutter criterion. The key to this type of FBP method is the accurate identification of modal damping factors.Many algorithms have been previously developed for identifying flutter test modal parameters, including fast Fourier transform- (FFT-) based methods [4, 5], Random Decrement Technique (RDT) [6], natural excitation technique combined with the eigensystem realization algorithm (NExT-ERA) [7], time series analysis based on the Autoregressive (AR) model [8, 9], and Stochastic Subspace Identification (SSI) [10, 11]. Although some of these algorithms are effective, none is ideal; for instance, nonstationary measured data and low signal to noise ratio (SNR) affect the Fourier-based methods. The AR model requires an appropriate model order and algorithm to function properly; RDT functionality is limited by the number of the main modals contained in the structure [12]. SSI and NExT-ERA are problematic in terms of their spurious mode [13, 14]. The Hilbert-Huang Transform (HHT) proposed in 1998 is an adaptive scheme well-suited to nonlinear, nonstationary time series analysis; however, the mode mixing effect [15] which emerges when dealing with signals over multifrequencies in each frequency band severely limits its application to flutter test data processing due to the inherent density of modal problems. This paper proposes an improved HHT which applies to flutter test modal parameter identification.


## 2. Hilbert-Huang Transform Theory

HHT is a relatively new nonstationary and nonlinear signal processing approach which is not limited under linear and stationary spectral analysis theory based on the Fourier transform. HHT has two basic steps. First, the original signal is decomposed into a series of Intrinsic Mode Function (IMF) components via Empirical Mode Decomposition (EMD); second, a Hilbert transform is performed on each IMF component to obtain instantaneous amplitude and phase information.According to classical Fourier theory, the local frequency only can be defined when there is at least one complete sine or cosine oscillation. This strict definition is not suitable for nonstationary signals which have ever-changing frequencies. A Hilbert transformation for the original signalx(t) is performed as follows to reveal the instantaneous frequency of the signal at any moment:(1)yt=Hxt=1πP∫-∞∞xτt-τdτwhere P denotes the Cauchy principal value; the transformation is available for all LP classes. According to this definition, when a complex conjugate is formed by x(t) and y(t), the analytical signal z(t) is(2)zt=xt+jyt=atejθtwhere a(t)=x2(t)+y2(t) denotes the instantaneous amplitude and θ(t)=artan⁡(y(t)/x(t)) denotes the instantaneous phase.In this context, the instantaneous frequency can be defined as follows:(3)ft=dθt2π·dtIn a sense, the local properties ofx(t) are emphasized because the Hilbert transform is defined as a convolution between x(t) and 1/t under (1). The local characteristics are further expressed as a polar coordinate in (2), which is the optimal approximation of the trigonometric functions with variable amplitude and phase.Although the definition of the instantaneous frequency and phase is given, some extra conditions are required to ensure the instantaneous frequencies have physical significance. The signals must be symmetrical with a local mean of zero and must have the same zero crossings and extreme points. A general signal can be described by the Hilbert transform in the frequency domain only after EMD.The EMD decomposition is based on the assumption that any signal is comprised of a series of different IMF components. The envelope of each IMF is defined by the local maximum and the minimum is symmetrical about the abscissa. The number of extreme points and zero crossing points should be equal or no more than one in the entire data sequence. Ideally, each IMF component contains a single modal and the number of IMFs equals the modal number of the original signal.The EMD steps to decompose any signalx(t) are as follows.(a)
Seek all the maximum and minimum points.(b)
Use cubic spline interpolation (CSI) to interpolate every maxima and minima point sequence to obtain the upper and lower envelopesU(t) and L(t).(c)
Calculate the average envelope linem(t) by m(t)=(U(t)-L(t))/2 and extract the details by c(t)=x(t)-m(t).(d)
Determine whetherc(t) meets the two IMF conditions given above. If so, c(t) is an IMF, as imfi. Otherwise, let x(t)=c(t) and repeat steps (a)-(d) until the conditions are met.(e)
Let the residualri+1=x(t)-imfi be the new signal to be decomposed, and repeat steps (a)-(d) for other IMFs. Repeat the steps above until the residual is a monotonic signal or falls below a certain threshold; then the decomposition is complete.The EMD procedure is also shown in Figure1. The IMF, per its theoretical definition, meets the conditions of the Hilbert transform as the instantaneous frequencies have physical significance. Thus, the HHT can be used for complex nonstationary signal analysis.Figure 1
EMD method.


## 3. Improved HHT and Modal Parameter Identification

### 3.1. Adjacent Mode Separation

Mode mixing problems are inherent to the EMD method. When the signal contains multiple modals with similar frequencies, the energy between the IMF components is widely divergent and the instantaneous frequencies of the IMFs may intermittently fluctuate, so the IMF components will not be decomposed correctly. As an example, consider the signalx(t)=A1cos⁡2πf1+A2cos⁡(2πf2). If the two frequencies are adjacent, the low-frequency component will be mixed into the high-frequency IMF component during the EMD process to create a mode mixing effect. After several simulation experiments, when 0.5<f1/f2<1, mode mixing is inevitable during the decomposition process and grows increasingly severely as the ratio approaches 1; this is a serious problem in regard to flutter test data analysis due to the dense modal and modal coupling.Without loss of generality, givenf1<f2, for any positive number f0 less than f1, (f1-f0)/(f2-f0)<f1/f2. The mode mixing effect can be reduced when the two close frequency components shift leftward simultaneously. Therefore, any signal X(t)=∑i=1nxi(t) can be processed via the following steps.(a) Construct a Complex Analytical Signal (4) z t = X t + j H X t = ∑ i = 1 n x i t + j H x i t = ∑ i = 1 n a i t e j θ i tin which H{•} is the Hilbert operator.(b) Modulation. Multiply the constructed complex analytical signal with the complex exponential signal w(t)=e-j2πf0·t.(5)z~t=zt·wt=Xt+jHXt·e-j2πf0·t=Xtcos⁡2πf0·t+HXtsin⁡2πf0·t+j-Xtsin⁡2πf0·t+HXtcos⁡2πf0·tEquation (5) shows that the modulated signal z~(t) is a complex signal. Its real and imaginary parts are a pair of Hilbert transform pairs. Therefore, if the original signal X(t) is processed as above, the corresponding modulated signal X~(t) can be obtained as follows:(6)X~t=realz~tand (5) can be rewritten as(7)z~t=zt·wt=X~t+jHX~t=∑i=1naitejθit-2πf0·tBy definition, the instantaneous frequency of each IMF component of the modulated signal is(8)f~it=dθit-2πf0·t2π·dt=fit-f0The instantaneous frequency is reduced by the above process. As long as the proper selection is made, the alias mode can be completely separated.


### 3.2. Modal Parameter Identification

The aeroelastic system is a typical multiple degree-of-freedom (DOF) dynamic system. The impulse response function can be expressed as follows:(9)Xt=∑i=1naie-ζiωitcos⁡ωdit+θiwhere ωi denotes modal natural frequencies, ζi denotes modal damping ratios, and ωdi denotes natural damping frequencies with ωdi=ωi1-ζi2.The proposed “improved” HHT method was designed to overcome the mode mixing problem caused by the dense modal and modal coupling for effective modal parameter identification in flutter tests. It works in the following stepwise process.(a) Select an appropriate complex exponential signal to modulate the flutter test signal to reduce the denseness of the adjacent modal by (5).(b) Obtain the IMF components by using EMD:(10)x~it=aie-ζiωitcos⁡ωdit+θi(c) Apply the Hilbert transform to each IMF component to calculate the modal parameters. The instantaneous amplitudeai(t) and instantaneous phase θi(t) can be expressed as(11)ait=aie-ζiωit(12)θit=ωdit+θiLogarithms and differential operators are performed for (11) and (12), respectively:(13)ln⁡ait=-ζiωit+ln⁡ai(14)ωit=dθitdt=ωdiThen the natural frequency and damping ratio are identified by straight line fitting with the instantaneous frequency ωi(t) and the logarithm of instantaneous amplitude ln⁡ai(t).(d) The IMF component is the decomposition ofX~(t), so the actual natural frequencies and damping ratios can be calculated from (15) and (16) by demodulation:(15)ln⁡ait=-ζiωi+ω0t+ln⁡ai(16)ωit=dθitdt-ω0=ωdi


## 3.1. Adjacent Mode Separation

Mode mixing problems are inherent to the EMD method. When the signal contains multiple modals with similar frequencies, the energy between the IMF components is widely divergent and the instantaneous frequencies of the IMFs may intermittently fluctuate, so the IMF components will not be decomposed correctly. As an example, consider the signalx(t)=A1cos⁡2πf1+A2cos⁡(2πf2). If the two frequencies are adjacent, the low-frequency component will be mixed into the high-frequency IMF component during the EMD process to create a mode mixing effect. After several simulation experiments, when 0.5<f1/f2<1, mode mixing is inevitable during the decomposition process and grows increasingly severely as the ratio approaches 1; this is a serious problem in regard to flutter test data analysis due to the dense modal and modal coupling.Without loss of generality, givenf1<f2, for any positive number f0 less than f1, (f1-f0)/(f2-f0)<f1/f2. The mode mixing effect can be reduced when the two close frequency components shift leftward simultaneously. Therefore, any signal X(t)=∑i=1nxi(t) can be processed via the following steps.(a) Construct a Complex Analytical Signal (4) z t = X t + j H X t = ∑ i = 1 n x i t + j H x i t = ∑ i = 1 n a i t e j θ i tin which H{•} is the Hilbert operator.(b) Modulation. Multiply the constructed complex analytical signal with the complex exponential signal w(t)=e-j2πf0·t.(5)z~t=zt·wt=Xt+jHXt·e-j2πf0·t=Xtcos⁡2πf0·t+HXtsin⁡2πf0·t+j-Xtsin⁡2πf0·t+HXtcos⁡2πf0·tEquation (5) shows that the modulated signal z~(t) is a complex signal. Its real and imaginary parts are a pair of Hilbert transform pairs. Therefore, if the original signal X(t) is processed as above, the corresponding modulated signal X~(t) can be obtained as follows:(6)X~t=realz~tand (5) can be rewritten as(7)z~t=zt·wt=X~t+jHX~t=∑i=1naitejθit-2πf0·tBy definition, the instantaneous frequency of each IMF component of the modulated signal is(8)f~it=dθit-2πf0·t2π·dt=fit-f0The instantaneous frequency is reduced by the above process. As long as the proper selection is made, the alias mode can be completely separated.


## 3.2. Modal Parameter Identification

The aeroelastic system is a typical multiple degree-of-freedom (DOF) dynamic system. The impulse response function can be expressed as follows:(9)Xt=∑i=1naie-ζiωitcos⁡ωdit+θiwhere ωi denotes modal natural frequencies, ζi denotes modal damping ratios, and ωdi denotes natural damping frequencies with ωdi=ωi1-ζi2.The proposed “improved” HHT method was designed to overcome the mode mixing problem caused by the dense modal and modal coupling for effective modal parameter identification in flutter tests. It works in the following stepwise process.(a) Select an appropriate complex exponential signal to modulate the flutter test signal to reduce the denseness of the adjacent modal by (5).(b) Obtain the IMF components by using EMD:(10)x~it=aie-ζiωitcos⁡ωdit+θi(c) Apply the Hilbert transform to each IMF component to calculate the modal parameters. The instantaneous amplitudeai(t) and instantaneous phase θi(t) can be expressed as(11)ait=aie-ζiωit(12)θit=ωdit+θiLogarithms and differential operators are performed for (11) and (12), respectively:(13)ln⁡ait=-ζiωit+ln⁡ai(14)ωit=dθitdt=ωdiThen the natural frequency and damping ratio are identified by straight line fitting with the instantaneous frequency ωi(t) and the logarithm of instantaneous amplitude ln⁡ai(t).(d) The IMF component is the decomposition ofX~(t), so the actual natural frequencies and damping ratios can be calculated from (15) and (16) by demodulation:(15)ln⁡ait=-ζiωi+ω0t+ln⁡ai(16)ωit=dθitdt-ω0=ωdi


## 4. Numerical Simulation and Application

### 4.1. Modal Parameter Identification

The effects of flutter boundary prediction are dependent on accurate modal parameter identification. The following simulation was primarily run to test the performance of HHT in terms of modal parameter estimation. The response of a 3DOF dynamic system was constructed according to (9) with the modal parameters shown in Table 1. The temporal history of the response waveform excited by the impulse signal is shown in Figure 2, where the duration is 1 s and sampling rate is 1000 Hz.Table 1
Modal parameters of 3DOF system.
The first modal The second modal The third modal Natural frequency (Hz)fi 50 70 150 Damping ratio (%)ζi 1.0 0.7 0.3Figure 2
Simulated signal temporal history.Since the first and second modal frequencies are adjacent,f1/f2=0.7, the traditional HHT method inevitably results in modal aliasing as shown in Figure 3, where each IMF component and its frequency spectrum are listed. Figure 3 also shows that the frequencies in the first and third IMFs are pure of any other frequency components, but the second contains both 70 Hz and 50 Hz frequency components.Figure 3
Traditional HHT decomposition results.In the modulation frequencyf0=30 Hz, all the frequencies contained are shifted by 30 Hz to the left; this reduces the ratio of f1/f2 to 0.5 and prevents mode mixing. The improved HHT was performed subsequently to obtain the decomposition results shown in Figure 4. Each IMF component was separated clearly with an individual modal. After the right decomposition, the instantaneous frequency and instantaneous amplitude curve of each IMF were calculated according to (13) and (14) as shown in Figure 5. Linear fitting was then conducted separately to obtain each modal natural frequency and damping ratio of the modulated signal. The final modal parameters of the original signal were obtained by demodulation according to (15) and (16), as shown in Table 2. Compared with the true values, the results in Table 2 indicate that the proposed method yields very accurate modal parameters. There is slight (tolerable) error attributable to endpoint effects and short sampling length due to the envelope fitting method.Table 2
Estimated modal parameters.
Natural frequency (Hz)fi Damping ratio (%)ζi First modal Second modal Third modal First modal Second modal Third modal True value 50 70 150 1.0 0.7 0.3  Original HHT 49.239 70.522 4.357 0.913 0.532 0.8 Error /% 1.52 0.75 14997 8.7 24 236  Improved method 49.305 69.185 149.340 0.949 0.681 0.301 Error /% 1.39 1.16 0.44 5.10 2.71 0.33Figure 4
Improved HHT decomposition results.Figure 5
Linear curve fitting of logarithmic IMF amplitude.


### 4.2. Flutter Test Application

We next applied the proposed method to actual flutter flight test data (China Flight Test Establishment) to determine its practical applicability. Four main acceleration sensors were installed on the auxiliary tank and wingtip of the aircraft symmetrically, and the structural responses were excited by a small rocket as pulse excitations. The flutter flight test was conducted at six dynamic airspeeds from 600 Km/h to 780 Km/h. The data digitized at sampling rate of 64 Hz and 1024 data points were supplied for the improved HHT method. The flutter boundary designed for the aircraft is 850 Km/h. The lowest three natural frequencies measured by the vibration test and the FEM analysis are given in Table3. The signals from the acceleration sensor located at the left wingtip were selected and the power spectral density of every airspeed was determined as shown in Figure 6. The data have strong noise in the lower frequency rage.Table 3
Wing modal frequencies.
Mode No. FEM (Hz) Vibration experiment (Hz) 1 0.73 0.81 2 4.93 5.15 3 22.12 21.79Figure 6
Spectrums of vibration signals at every airspeed.The third mode was not effectively activated due to limited excited energy. The first and second modal parameters including the frequency and damping at each test stage were estimated by the improved HHT method as shown in Table4; the second mode is the primary flutter mode with smaller damping and higher energy. We conducted extrapolation and curve fitting based on the second modal damping factors; the FBP result is 824 Km/h, as shown in Figure 7, which is consistent with the design value.Table 4
Changes in modal parameters with wind speed.
airspeed(km/h) Natural frequency (Hz) Damping ratio (%) f 1 f 2 ζ 1 ζ 2 600 0.415 3.315 7.521 2.239 650 0.347 4.235 53.51 4.175 700 0.375 4.552 21.72 3.721 730 0.647 4.731 12.37 3.109 760 0.524 4.679 32.73 2.829 780 0.672 5.044 23.15 2.222Figure 7
FBP results of flight flutter tests.


## 4.1. Modal Parameter Identification

The effects of flutter boundary prediction are dependent on accurate modal parameter identification. The following simulation was primarily run to test the performance of HHT in terms of modal parameter estimation. The response of a 3DOF dynamic system was constructed according to (9) with the modal parameters shown in Table 1. The temporal history of the response waveform excited by the impulse signal is shown in Figure 2, where the duration is 1 s and sampling rate is 1000 Hz.Table 1
Modal parameters of 3DOF system.
The first modal The second modal The third modal Natural frequency (Hz)fi 50 70 150 Damping ratio (%)ζi 1.0 0.7 0.3Figure 2
Simulated signal temporal history.Since the first and second modal frequencies are adjacent,f1/f2=0.7, the traditional HHT method inevitably results in modal aliasing as shown in Figure 3, where each IMF component and its frequency spectrum are listed. Figure 3 also shows that the frequencies in the first and third IMFs are pure of any other frequency components, but the second contains both 70 Hz and 50 Hz frequency components.Figure 3
Traditional HHT decomposition results.In the modulation frequencyf0=30 Hz, all the frequencies contained are shifted by 30 Hz to the left; this reduces the ratio of f1/f2 to 0.5 and prevents mode mixing. The improved HHT was performed subsequently to obtain the decomposition results shown in Figure 4. Each IMF component was separated clearly with an individual modal. After the right decomposition, the instantaneous frequency and instantaneous amplitude curve of each IMF were calculated according to (13) and (14) as shown in Figure 5. Linear fitting was then conducted separately to obtain each modal natural frequency and damping ratio of the modulated signal. The final modal parameters of the original signal were obtained by demodulation according to (15) and (16), as shown in Table 2. Compared with the true values, the results in Table 2 indicate that the proposed method yields very accurate modal parameters. There is slight (tolerable) error attributable to endpoint effects and short sampling length due to the envelope fitting method.Table 2
Estimated modal parameters.
Natural frequency (Hz)fi Damping ratio (%)ζi First modal Second modal Third modal First modal Second modal Third modal True value 50 70 150 1.0 0.7 0.3  Original HHT 49.239 70.522 4.357 0.913 0.532 0.8 Error /% 1.52 0.75 14997 8.7 24 236  Improved method 49.305 69.185 149.340 0.949 0.681 0.301 Error /% 1.39 1.16 0.44 5.10 2.71 0.33Figure 4
Improved HHT decomposition results.Figure 5
Linear curve fitting of logarithmic IMF amplitude.


## 4.2. Flutter Test Application

We next applied the proposed method to actual flutter flight test data (China Flight Test Establishment) to determine its practical applicability. Four main acceleration sensors were installed on the auxiliary tank and wingtip of the aircraft symmetrically, and the structural responses were excited by a small rocket as pulse excitations. The flutter flight test was conducted at six dynamic airspeeds from 600 Km/h to 780 Km/h. The data digitized at sampling rate of 64 Hz and 1024 data points were supplied for the improved HHT method. The flutter boundary designed for the aircraft is 850 Km/h. The lowest three natural frequencies measured by the vibration test and the FEM analysis are given in Table3. The signals from the acceleration sensor located at the left wingtip were selected and the power spectral density of every airspeed was determined as shown in Figure 6. The data have strong noise in the lower frequency rage.Table 3
Wing modal frequencies.
Mode No. FEM (Hz) Vibration experiment (Hz) 1 0.73 0.81 2 4.93 5.15 3 22.12 21.79Figure 6
Spectrums of vibration signals at every airspeed.The third mode was not effectively activated due to limited excited energy. The first and second modal parameters including the frequency and damping at each test stage were estimated by the improved HHT method as shown in Table4; the second mode is the primary flutter mode with smaller damping and higher energy. We conducted extrapolation and curve fitting based on the second modal damping factors; the FBP result is 824 Km/h, as shown in Figure 7, which is consistent with the design value.Table 4
Changes in modal parameters with wind speed.
airspeed(km/h) Natural frequency (Hz) Damping ratio (%) f 1 f 2 ζ 1 ζ 2 600 0.415 3.315 7.521 2.239 650 0.347 4.235 53.51 4.175 700 0.375 4.552 21.72 3.721 730 0.647 4.731 12.37 3.109 760 0.524 4.679 32.73 2.829 780 0.672 5.044 23.15 2.222Figure 7
FBP results of flight flutter tests.


## 5. Conclusion

An improved HHT method was proposed in this paper for identifying flight flutter modal parameters and predicting flutter boundaries. Typical mode mixing effect issues were resolved by applying a modulation algorithm that solely requires a structural response in the form of measured acceleration. Simulation and flight flutter test data results altogether indicated that this approach is simple, effective, and feasible for flutter test data processing.



---

*Source: 3496870-2018-08-12.xml*