Time–frequency analysis for parametric and non-parametric identification of nonlinear dynamical systems

Mechanical Systems and Signal Processing 36 (2013) 332–353

Contents lists available at SciVerse ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Time–frequency analysis for parametric and non-parametric identification of nonlinear dynamical systems P. Frank Pai n Department of Mechanical and Aerospace Engineering, University of Missouri, Columbia, MO 65211, United States

a r t i c l e in f o

abstract

Article history: Received 16 February 2012 Received in revised form 27 November 2012 Accepted 6 December 2012 Available online 18 January 2013

This paper points out the differences between linear and nonlinear system identification tasks, shows that time–frequency analysis is most appropriate for nonlinearity identification, and presents advanced signal processing techniques that combine time–frequency decomposition and perturbation methods for parametric and non-parametric identification of nonlinear dynamical systems. Hilbert–Huang transform (HHT) is a recent datadriven adaptive time–frequency analysis technique that combines the use of empirical mode decomposition (EMD) and Hilbert transform (HT). Because EMD does not use predetermined basis functions and function orthogonality for component extraction, HHT provides more concise component decomposition and more accurate time– frequency analysis than the short-time Fourier transform and wavelet transform for extraction of system characteristics and nonlinearities. However, HHT’s accuracy seriously suffers from the end effect caused by the discontinuity-induced Gibbs’ phenomenon. Moreover, because HHT requires a long set of data obtained by highfrequency sampling, it is not appropriate for online frequency tracking. This paper presents a conjugate-pair decomposition (CPD) method that requires only a few recent data points sampled at a low-frequency for sliding-window point-by-point adaptive time–frequency analysis and can be used for online frequency tracking. To improve adaptive time–frequency analysis, a methodology is developed by combining EMD and CPD for noise filtering in the time domain, reducing the end effect, and dissolving other mathematical and numerical problems in time–frequency analysis. For parametric identification of a nonlinear system, the methodology processes one steady-state response and/or one free damped transient response and uses amplitude-dependent dynamic characteristics derived from perturbation analysis to determine the type and order of nonlinearity and system parameters. For non-parametric identification, the methodology uses the maximum displacement states to determine the displacement– stiffness curve and the maximum velocity states to determine the velocity–damping curve. Numerical simulations and experimental verifications of several nonlinear discrete and continuous systems show that the proposed methodology can provide accurate parametric and non-parametric identifications of different nonlinear dynamical systems. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Time–frequency analysis Parametric and non-parametric nonlinearity identification Conjugate-pair decomposition Perturbation analysis Hilbert–Huang transform Signal processing techniques

1. Introduction Mechanical systems are often built with unwanted or wanted nonlinearities [1–10]. Large elastic deformations introduce geometric nonlinearity, deformation-dependent material properties result in material nonlinearity, and other

n

Tel.: þ 1 573 884 1474; fax: þ 1 573 884 5090. E-mail address: [email protected]

0888-3270/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ymssp.2012.12.002

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

333

nonlinearities can be caused by backlash, clearances between mounting brackets, geometric constraints on deformation, misalignment of substructures, dry friction, and many types of nonlinear hysteretic damping, including aerodynamic damping, damping of shock absorbers for vehicles, and material damping of shape memory alloys and other materials [1–6]. For example, highly flexible deployable/inflatable structures are often used in space structural systems, advanced vehicles, many mechanical and civil systems, and medical devices and implants in order to reduce weight/space, increase speed/efficiency, and/or provide special mechanisms [2,4,6]. Such structures are built with wanted geometric nonlinearities. Damage introduces extra nonlinearities, and aging aggravates nonlinearities. Hence, dynamics of actual mechanical systems is often nonlinear [7–10]. For example, dry friction between fractured surfaces is a nonlinear effect, and the breathing of a crack causes nonlinear intermittent transient response. Also, demountable/retractable civil structures are built with nonlinearities caused by loose joints, and the usage/aging aggravates such nonlinearities. Because parametric models of such nonlinearities are often difficult to determine, non-parametric system identification is more feasible for such systems. Nonlinearity identification and damage inspection are very challenging reverse engineering tasks. For several decades structural engineers have been developing dynamics-based methods for fast and reliable system identification and damage inspection of large structures [7–10]. For reverse engineering, it is important to design experiments to limit the number of system parameters involved in each set of experimental data in order to increase the possibility of unique identification results. But, it is even more important to have a signal processing and data mining method that can extract from each set of experimental data as many system parameters as possible in order to reduce the time and cost of experiments. Hence, signal processing plays the key role in nonlinearity identification, but the challenge is how to obtain accurate extraction/ tracking of dynamic characteristics (e.g., natural frequencies, mode shapes, mass, damping, stiffness, and external loading) from nonlinear noise-contaminated signals. Most signal processing methods for dynamics characterization require frequency- and/or time-domain signal processing techniques. Transforming time-domain data into the frequency domain introduces errors and numerical noise, and some methods even require the frequency domain data to be transformed back into the time domain, which introduces even more errors. For transformation from the time domain to the frequency domain, the first step is to choose a complete set of basis functions with each one having a pre-determined frequency (e.g., Fourier transform) or a narrow frequency band (e.g., wavelet transform), and the second step is to perform convolution computation using the chosen basis functions to extract components similar to the basis functions. Unfortunately, this approach is not suitable for a nonlinear and/or nonstationary signal because its frequency changes with time. Dynamics with a changing frequency is sometimes called high dynamics. For high dynamics, adaptive basis functions need to be used, and they can only be derived from the signal itself. Current literature review indicates that Hilbert–Huang transform (HHT) is better than shorttime Fourier transform, Wigner–Ville distribution, and wavelet transform for time–frequency analysis of nonlinear nonstationary signals because HHT does not use predetermined basis functions and function conformality for component extraction [11,12]. Moreover, HHT provides concise decomposition of a signal into just a few intrinsic mode functions (IMFs) by using the empirical mode decomposition (EMD), and each extracted IMF is often physically meaningful [11,12]. Furthermore, HHT provides more accurate time-varying amplitudes and frequencies of each IMF by using Hilbert transform (HT). Unfortunately, the accuracy and capability of HHT suffers from the end effect and several mathematical and numerical problems [13]. Moreover, how to identify and estimate nonlinearities from time-varying amplitudes and frequencies is an important but challenging task [14]. Hence, new developments in time–frequency analysis and nonlinearity identification are needed. Because dynamic characteristics of nonlinear systems change with time, this work shows that time–frequency analysis is most appropriate for their characterization and identification. Moreover, this work develops and experimentally verifies time–frequency signal processing methods and identification techniques for parametric and non-parametric identification of nonlinear dynamical systems. We show that the type, order, and magnitude of nonlinearity of a nonlinear system can be determined by performing time–frequency analysis of just a few transient responses and/or steady-state responses to harmonic excitations. 2. Identification of nonlinear systems Nonlinear systems have dynamic characteristics very different from those of linear systems in many aspects and in time and frequency domains, and hence approaches for their identification are different. For example, the steady-state response of a linear system to a harmonic excitation is not affected by initial conditions and it vibrates at the excitation frequency with a constant amplitude. Hence, the start-up transient part due to non-zero initial conditions can be well separated from the steady-state part, which enables the use of just one forced transient response for time-domain system identification [15]. On the other hand, a harmonically excited response of a nonlinear system cannot be decomposed into a constant amplitude steady-state part and a transient part only due to non-zero initial conditions. Hence, for a nonlinear system, a transient response and a steady-state response need to be separately obtained in order to perform accurate system identification, as shown later in Section 4. For a nonlinear system, because the principle of superposition is invalid and its steady-state response frequency can be different from the excitation frequency, system identification is difficult to be performed in the frequency domain. Next we

334

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

show that time–frequency analysis provides a unique tool for parametric and non-parametric identification of nonlinear systems. Parametric identification of nonlinearity is to determine the types, order, and functional form of nonlinearity and parameter values. First we show that parametric identification of a nonlinear dynamical system can be performed with the help of perturbation solutions. 2.1. Perturbation solutions for parametric identification To demonstrate the use of perturbation analysis for nonlinearity identification, we consider the following Duffing oscillator subjected to a harmonic excitation at a frequency O close to its linear natural frequency o and the corresponding steady-state second-order perturbation solution [1]: u€ þ 2Bou_ þ o2 u þ au3 ¼ FcosðOtÞ, F ¼ F 0 =m, O  o uðt Þ ¼ a1 cosðOtfÞ þ a3 cosð3Ot3fÞ ¼ AcosðOtf þ YÞ,a3  AðtÞ 

B¼

32O2

5 a1

a3 sinð2Ot2fÞ a3  sinð2Ot2fÞ a1 þ a3 cosð2Ot2fÞ a1

^ ðt Þ  O þ Y _ ¼ Oþ O 3a 8o

aa31

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a21 þ a23 þ 2a1 a3 cosð2Ot2fÞ  a1 þ a3 cosð2Ot2fÞ

YðtÞ  tan1



ð1aÞ

2

2Oa23 þ 2Oa1 a3 cosð2Ot2fÞ 2Oa3  Oþ cosð2Ot2fÞ a1 a21 þ a23 þ 2a1 a3 cosð2Ot2fÞ

 2  3as 4  2 F a61  a1 þ s þ B2 o2 a21  ¼ 0, s  Oo 2o 4o

  3aa21 s  tan f 8o2 o

^ ¼ oþ o

3a 2 a 8o 1

ð1bÞ

ð1cÞ ð1dÞ

ð1eÞ

ð1fÞ

ð1gÞ

ð1hÞ

where f is the phase difference between the forcing function and the response, and B, a, F, and f are constants. The amplitude a1 is a nonlinear function of F,B,o, and a , as shown in Eq. (1f) from first-order perturbation analysis [1]. Eq. (1b) shows that u(t) consists of two synchronous harmonics (i.e., reaching maxima at the same time), but, because ^ varying at a frequency 2 O. a1 ba3 , u(t) would appear as one distorted harmonic with an amplitude A and a frequency O This amplitude- and frequency-modulation phenomenon can be used to determine the order (cubic or other) of ^ and A are at their maxima when u(t) is at its maxima or minima nonlinearity. Moreover, if a 40, a3/a1 40, O ^ and A are at their minima when uðtÞ ^ _ (i.e.,u_ ¼ 0), and O is close to its maxima or minima (i.e., u¼0). If a o0, a3/a1 o0, O ^ and A are at their maxima when uðtÞ _ and A are at their minima when u(t) is at its maxima or minima (i.e., u_ ¼ 0), and O is close to its maxima or minima (i.e., u ¼0). This phenomenon can be used to determine the type (hardening or softening) of nonlinearity, and the magnitude of nonlinearity can be estimated using Eq. (1b) and the time-domain data a1(t) and a3(t) as a ¼ 32O2 a3 =a31 . Moreover, one can use Eq. (1g) to estimate the damping ratio B, use Eq. (1h) to estimate the free ^ , and use Eq. (1f) to obtain F and then estimate the mass as m¼ F0/F, undamped amplitude-dependent natural frequency o where F0 is the known excitation amplitude. Numerical results show that all these characteristics under a harmonic excitation also exist in transient vibrations [16]. For systems with different orders of nonlinearity, similar perturbation solutions can be derived for identification of different nonlinearities [1,5]. For a system with quadratic nonlinearity, for example, its steady-state response can be shown to be a distorted harmonic with its amplitude and frequency varying at a frequency O. 2.2. Non-parametric identification If a system’s order and/or type of nonlinearity is unknown, non-parametric identification is needed. For a dynamical system, when the displacement u is at its maximum or minimum, u_ ¼ 0 and it is equivalent to a constant–displacement state that is used in experiments to determine the spring force. Similarly, when the velocity u_ is at its maximum or minimum, u€ ¼ 0 and it is equivalent to a constant–velocity state that is used in experiments to determine the damping force. Hence, for non-parametric identification of nonlinearities one can use the maximum displacement states to determine the displacement–stiffness curve and the maximum velocity states to determine the velocity–damping curve. For a free damped vibration of the following nonlinear system with unknown orders p and q:   ^ 2 uq ¼ 0 ^ u_ p þ o m u€ þ 2Bo ð2aÞ

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

335

the specific spring force (i.e., force per unit mass) Fs and the specific damping force Fd can be obtained as ^ 2 A^ Fs ¼ o

ð2bÞ

^ B^ Fd ¼ 2Bo

ð2cÞ

^ in Eq. (2b) are the instantaneous amplitude and frequency of u when u is at its extrema (i.e., peaks or where the A^ and o ^ and B^ in Eq. (2c) are the instantaneous damping ratio, frequency, and amplitude of u_ when u_ is at its valleys), and the B, o extrema. However, difficulties exist because different material points of a nonlinear continuous or discrete system undergoing arbitrary vibration may arrive at their maximum displacements and/or velocities at different times (i.e., the so-called complex operation deflection shape (ODS)). Such a complex ODS needs to be decomposed into normal ODSs in order for this approach to work, but this task is challenging, as shown next. 2.3. Challenging issues It is almost impossible to perform nonlinearity identification in the frequency domain alone. To show that we consider the following amplitude-modulated signal:     uðtÞ ¼ ð1 þ ecos ot Þcos Ot ¼ cos Ot þ e=2 cosðO þ oÞt þ e=2 cosðOoÞt ð3aÞ where O is the main vibration frequency, oð 5 OÞ a low modulation frequency, and e a small parameter. Eq. (3a) shows that u(t) consists of three different regular harmonics. For the following frequency-modulated signal, we have     uðtÞ ¼ cosðOt þ ecos ot Þ  1e2 =4 cos Ot e=2e3 =16 ½sinðO þ oÞt þsinðOoÞt     ð3bÞ  e2 =8 ½cosðO þ 2oÞt þ cosðO2oÞt þ e3 =48 ½sinðO þ 3oÞt þsinðO3oÞt where the Taylor expansion is obtained by treating ecos ot as a small parameter. Eq. (3b) shows that the Fourier spectrum of u(t) would consist of several small but uniformly spaced discrete spectral lines around the main frequency O, but the actual frequency of u(t) varies within a continuous frequency bandwidth from O  e o to O þ e o. Eq. (3a,b) show that a distorted harmonic can be decomposed into many regular harmonics and hence multiple discrete spectral lines exist on the corresponding spectrum. Hence, it is difficult to perform quantitative nonlinearity identification using frequencydomain data. Moreover, Eq. (1b) shows that a distorted harmonic in the time domain is more meaningful and convenient than several regular harmonics for characterization of a nonlinear system. Fig. 1 shows the nonlinear frequency response function (FRF) obtained from the first-order perturbation solution (1f) and the direct numerical integration using the Runge–Kutta method. Note that the perturbation solution agrees well with the numerical one, especially around the linear resonance frequency. Direct numerical integrations show that u(t) is a regular harmonic when A is small, but u(t) becomes a distorted harmonic consisting of more and more regular harmonics when A increases, as the perturbation solution predicts. As shown in Section 2.1, perturbation solutions can be used for estimation of nonlinearity. However, because multiple solutions may coexist for one value of O, nonlinear FRFs cannot be experimentally obtained using impulse excitations, random excitations, periodic random excitations, pseudo-random excitations, burst random excitations, chirp excitations, or other broad-band excitations [3]. A broad-band excitation tends to simultaneously excite several modes. If each modal vibration is a distorted harmonic due to nonlinearity, the Fourier 1.0

0.8

A

0.6

0.4

0.2

0 0.4

0.6

0.8

1.0

1.2

1.4

1.6

Ω /ω Fig. 1. Nonlinear FRF of Eq. (1a) with B ¼0.05, o ¼ 1, F¼ 0.1, and a ¼1, where red dots are from direct numerical integrations and black solid and broken lines represent stable and unstable perturbation solutions, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

336

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

spectrum loses its physical meaning because each distorted harmonic is represented by several spectral lines and hence couplings exist and peaks due to higher harmonics (see Eq. (1b)) may be erroneously identified as extra modes. More seriously, it is impossible to determine the order and magnitude of nonlinearity from Fourier spectra because the unstable branch in Fig. 1 cannot be experimentally obtained and the backbone curve (from Eq. (1h)) cannot be determined. Furthermore, Fourier spectra cannot reveal transient nonlinearities because it is an averaged representation of the whole sampled signal duration. To identify damping in the time domain requires free damped vibrations or forced responses to a non-stationary excitation. Then, how to extract accurate time-varying dynamic characteristics from a nonlinear non-stationary signal is very challenging because, except the distorted harmonic response to a harmonic excitation, there are many other nonlinear phenomena that may coexist. For example, a nonlinear system has amplitude-dependent natural frequencies (see Eq. (1h)), and it may vibrate at a frequency much lower than an externally applied harmonic excitation [17]. Other nonlinear phenomena include multiple-mode vibrations caused by external resonances and/or internal resonances (i.e., interwave modulation), intermittent nonlinear response, bifurcation, and chaotic vibrations [1,2,5]. Hence, the key for successful nonlinearity identification is the ability to extract accurate time-varying dynamic characteristics from nonstationary multiple-component responses. Fourier transform (FT) is often used to transform time-domain data into frequency-domain data to identify system characteristics from the obtained spectrum, but a spectrum cannot show time-varying characteristics of the processed signal [3,18]. To overcome this problem, short-time Fourier transform (STFT, spectrogram), wavelet transform (WT), and Hilbert–Huang transform (HHT) have been developed for time–frequency analysis. STFT uses windowed regular harmonics and function orthogonality to simultaneously extract time-localized regular harmonics, and WT uses scaled and shifted wavelets and function orthogonality to simultaneously extract time-localized components similar to the wavelets [19,20]. Unfortunately, STFT and WT cannot perform adaptive time–frequency analysis and cannot provide accurate time-varying frequencies and amplitudes because of the use of pre-determined basis functions and function orthogonality for simultaneous extraction of components [19–22]. Although the matching pursuit method has no cross-term problems and is capable of adaptive time–frequency analysis, its ‘‘best matching’’ search for a linear combination of time–frequency atoms from an over-complete dictionary in order to have most of the combination coefficients being zero is a computationally expensive nonlinear process [23]. Moreover, choosing an appropriate over-complete dictionary is a difficult, subjective task, and its accuracy in signal decomposition is still an open question. For example, the perturbation solution shown in Eq. (1b) can be a good time–frequency atom for use in matching pursuit analysis of signals from dynamic systems with cubic nonlinearity, but it is difficult to build a dictionary by simple dilation, translation, and modulation of this atom. On the other hand, HHT uses the apparent time scales revealed by the signal’s local maxima and minima to sequentially sift components of different time scales and can provide accurate and fast adaptive time–frequency analysis [11,12]. Because HHT does not use predetermined basis functions and function orthogonality for component extraction, each of its extracted intrinsic mode functions (IMFs) often has physical meaning (as shown next in Section 3.2), and it provides accurate instantaneous amplitudes and frequencies of extracted components for accurate estimation of system characteristics and nonlinearities [11,24]. However, the accuracy of HHT analysis suffers from several mathematical and numerical effects that require further study and improvement [11,13]. 3. Time–frequency analysis techniques First we show the physical meaning of the time-varying instantaneous frequency (IF) and instantaneous amplitude because they play the key role in time–frequency analysis. 3.1. Instantaneous frequency and amplitude To show the calculation of instantaneous frequency (IF) of a general signal, we consider a time function u(t) defined over 0 rt rT and extended by repeating itself to become a periodic function. The extended function can be expressed in terms of regular harmonics with constant amplitudes and frequencies by using the discrete Fourier transform (DFT) as [18] uðtÞ ¼

N=2 X

ðai cosoi t þbi sinoi t Þ ¼ a0 þ 2

i ¼ N=2

ai 

N=2 X

i¼1

ðai cosoi t þ bi sinoi t Þ ¼ a0 þ

N=2 X

Ai cosyi

i¼1

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N N 1 X 1 X b 2 uk cosoi t k , bi  uk sinoi t k , Ai  2 a2i þ bi , yi  oi ttan1 i N k¼1 N k¼1 ai

ð4Þ

where oi  2pi/T ,uk u(tk), tk  (k  1)Dt, k¼1,...,N, N is the total number of samples (assumed to be even here), Dt is the sampling interval, 1/Dt is the sampling frequency, T ( ¼NDt) is the sampled period (or signal duration), Df ( ¼1/T) is the frequency resolution, 1/(2Dt) ( ¼N/(2T)) is the Nyquist (or maximum) frequency, and ai and bi are constant amplitudes of the regular harmonics. Eq. (4) shows that the conformality (or orthogonality) between u(t) and cosoi t and sinoi t is used to calculate the timeaveraged amplitude of each regular harmonic component. This is a brute-force method because regular harmonics of uniformly spaced frequencies are chosen as basis functions and u(t) is expressed as a summation of all N/2 regular

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

337

harmonics of constant amplitudes and frequencies. Eq. (4) also shows that u(t) is assumed to be a periodic function having a period T although the original function may not be periodic. Because the lowest-frequency harmonic in Eq. (4) has a frequency o1 ¼2p/T, the T needs to cover at least one period of the signal’s lowest-frequency component. Because uN þ 1 ¼u1 is implicitly assumed in DFT, if sudden change (i.e., discontinuity) exists between uN and u1 due to sampling and/ or transient response, the o1 harmonic is not a real component. Moreover, in order to account for this discontinuity, many high-frequency harmonics exist and result in the so-called leakage problem in the frequency domain and the end effect (also called Gibbs’ pffiffiffiffiffiffiffi effect) in the time domain at the two data ends. The leakage makes the Fourier spectrum U i ð  ai jbi ,j  1Þ difficult to understand. Gibbs’ effect happens because of the use of continuous functions (i.e., cosoi t and sinoi t in Eq. (4)) to fit a discontinuous function, and it results in overshoots around the discontinuities [18,25].   Some researchers used the Fourier spectrum Ui to compute the power spectral density (PSD) Si  U i U ni , bi-spectrum Bij ð  U i U j U ni þ j Þ, and tri-spectrum T ijk ð  U ni U nj U nk U i þ j þ k Þ, and use them for nonlinearity identification [26,27]. Here, U ni denotes the complex conjugate of Ui and Ui þ j þ k  U(oi þ oj þ ok). The Si provides information on the second-order properties (i.e., energy) of the signal. Bij and Tijk can reveal the existence of super-harmonic, sub-harmonic, and combination resonances, and hence they can be used for detection and rough estimation of quadratic and cubic nonlinearities, respectively [27]. However, because of the leakage problem and a distorted harmonic being presented by several spectral lines in the Fourier spectrum (see Sections 2.1 and 2.3), these high-order spectra cannot provide accurate and well separated quantities for nonlinearity identification, especially if the system is not weakly nonlinear [3,26,27]. In order to reduce the leakage and Gibbs’ effect, the sampled signal duration T needs to cover at least three cycles of the lowest-frequency component to be extracted [28]. Moreover, Shannon’s sampling theory requires more than two samples per period of the component to be extracted in order to capture the component’s frequency. Hence, more than six samples covering three periods are needed in order to extract the correct but averaged frequency. More seriously, a small T results in a lowfrequency resolution Df(¼1/T), and a small N results in a small maximum frequency (¼N/(2T)) that can be revealed. Hence, T and N are often numerically increased by zero-filling in using the short-time Fourier transform (STFT) to capture time-varying amplitudes and frequencies. Unfortunately, different lengths of zero-filling result in different values for the sampled signal duration T and hence different oi(¼2pi/T) for the basis functions in Eq. (4). To match with the zero-valued section of the modified function and the discontinuity due to zero-filling, additional artificial harmonics are often introduced into Eq. (4). Hence, this frequency domain method by STFT is not able to provide accurate time-varying frequency and amplitude at all. Another way to eliminate the leakage problem and Gibbs’ effect is to extend the two data ends with smooth curves to _ and u€ continuous at t ¼0 & T. Of course, this process can only be done case by case, and it is impossible to make u, u, program such a signal processing process for arbitrary signals. However, it is theoretically possible to exactly present an arbitrary continuous signal in terms of regular harmonics of constant frequencies and amplitudes by the Fourier transform with appropriate sampling and data extension. After Eq. (4) is obtained, because Ai and y_ i are constant, the exact conjugate part of u(t) can be obtained using Hilbert Transform (HT) as [11,29] ! N=2 N=2 X X Ai cosyi ¼ Ai sinyi ð5Þ vðtÞ ¼ HT a0 þ i¼1

i¼1

Moreover, it follows from Eqs. (4) and (5) that _ ¼ uðtÞ

N=2 X

_ ¼ oi Ai sinyi , vðtÞ

i¼1

N=2 X

oi Ai cosyi

ð6Þ

i¼1

Then, the instantaneous frequency (IF), o, can be defined using the phase angle y on the u–v plane as   _ uv _ dy d tan1 v=u uv ¼ 2 ¼ o¼ dt dt u þ v2

ð7Þ

Even without Gibbs’ effect the IF calculated from Eq. (7) is inaccurate or even erroneous. To show that we consider the following signal u(t) consisting of two regular harmonics having constant amplitudes a1 and a2: uðtÞ ¼ a1 cosy1 þ a2 cosy2

ð8aÞ

where y1  o1t, y2  o2t, and o1 and o2 are constant. It follows from Fig. 2a and the law of cosines that qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a1 siny1 þ a2 siny2 a12 ¼ a21 þ a22 þ2a1 a2 cosðy1 y2 Þ, y12 ¼ tan1 a1 cosy1 þ a2 cosy2 1 2

1 2

o  y_ 12 ¼ ðo1 þ o2 Þ þ ðo1 o2 Þ

a21 a22 a212

ð8bÞ

ð8cÞ

if a1 ¼a2 and o2 b o1  0, uðtÞ  u þ a2 coso2 t with u  a1 cosy1 being a constant around t¼tn. The IF should be o E o2 because the instantaneous trajectory circles around Point O in Fig. 2a with a frequency close to o2. On the other hand, it follows from Eq. (8c) that o ¼(o1 þ o2)/2E o2/2, which is an IF with respect to Point O1. Because o1 E0, the period 2p/o1-N and any finite-length sampling will present the signal as uðtÞ  u þ a2 coso2 t with u being a constant. Hence, vðtÞ ¼ HT ðuðt ÞÞ  a2 sino2 t, and the actual frequency should be close to o2. Obtaining v ¼ HTðuÞ ¼ a1 siny1 þ a2 siny2 requires an infinite length of sampling in order to present the signal as uðtÞ ¼ a1 cosy1 þa2 cosy2 , which is almost impossible to be

338

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

tn+1 tn tn−1

a2

v a12

O

θ2 O2

a1

O1

θ12

θ1

u

u

v

(u,v) a ˜˜ (u,v)

+

θ

a3

θ3 a2

a1 a0

θ1

θ2 u

Fig. 2. Trajectory of a signal on the phase plane: (a) two-harmonic signal, and (b) three-harmonic signal.

numerically realized. The erroneous prediction of Eq. (8c) results from the assumptions uðtÞ ¼ a1 cosy1 þ a2 cosy2 and vðtÞ ¼ a1 siny1 þ a2 siny2 and the use of Point O1 as the instantaneous trajectory center. If the slow component a1 coso1 t is known and temporarily fixed or separated, the modified trajectory will center at point O2 and the frequency will be o2. All the mistakes and paradoxical theories about IF in the literature are mainly caused by the co-existence of these three centers O, O1, and O2 in Fig. 2a and the difficulty in tracking them on the time and frequency domains. Fig. 2a shows that, without removing the slow moving component, the instantaneous trajectory circles around Point O with a frequency close to o2. Hence, for a general signal consisting of multiple harmonic components, one can define its instantaneous frequency to be the circular instantaneous frequency (cIF) of the curvature radius of the signal’s trajectory on the phase plane [24]. It follows from Fig. 2b that the curvature, cIF, and instantaneous center of the trajectory on the u–v (phase) plane can be obtained as Curvature :

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi € u€ v_ u_ v u_ v_ v_ u_ 1 dit ¼ dia ¼ , V  u_ 2 þ v_ 2 , it  i þ j, ia  i j 3 V V V V a Vdt V

Circular IF : o ¼

€ u€ v_ uv _ uv _ dy u_ v u_ v atan1 ¼ a 2 , y  tan1 v_ u dt u þ v2 V2

~ v~ Þ ¼ ðuacosy,vasinyÞ Instantaneous Center : ðu,

ð9Þ

where a is the radius of curvature, i is the unit vector along the u-axis, j is the unit vector along the v-axis, it is the unit vector tangent to the trajectory, and ia is the outward unit vector normal to the trajectory. Because Eq. (9) shows that € u€ v_ ), the curvature and cIF have the same sign, if the is assumed to be always positive (i.e., a  V 3 = u_ v radius of curvature 2 € u€ v_ =V ). This concept of circular instantaneous frequency enables the the cIF will also be always positive (i.e., o ¼ u_ v development of time-domain-only techniques for online frequency tracking [24]. However, because a general signal of a dynamical system of multiple degrees of freedom contains multiple modal vibrations, its cIF varies dramatically and is not useful for system identification if modal vibrations are not separated [24]. If a general nonlinear signal can be decomposed into sort-of modal vibration components without moving average, each component has no local extrema within each fundamental period and no local loops on the phase plane, each component’s referred instantaneous frequency (rIF) with respect to the origin on the phase plane may be non-circular but is always non-negative, and the time-varying rIF and referred instantaneous amplitude (rIA) are convenient for combining the perturbation analysis shown in Section 2.1 for system identification. The empirical mode decomposition (EMD) of Hilbert– Huang transform is a valuable method for decomposing a general nonlinear nonstationary signal into zero-mean intrinsic mode functions (IMFs), and Hilbert transform (HT) enables accurate calculation of rIF and rIA of each IMF, as shown next.

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

339

3.2. Hilbert–Huang transform Hilbert–Huang transform (HHT) combines EMD and HT for time–frequency analysis of nonlinear non-stationary signals [11,12]. EMD is a data-driven signal decomposition technique that sequentially extracts zero-mean regular/distorted harmonics from a signal, starting from high- to low-frequency components, and it is a dyadic filter that can perform adaptive time–frequency analysis. The first step of HHT is to use EMD to sequentially decompose a signal u(t) into n intrinsic mode functions (IMFs) ci(t) and a residual rn as [11,12] uðtÞ ¼

n X

ci ðtÞ þ r n ðtÞ,

ð10Þ

i¼1

where c1 has the shortest characteristic time scale and is the first extracted IMF. The characteristic time scale of c1 is defined by the time lapse between the extrema of u. Once the extrema are identified, compute the upper envelope by connecting all the local maxima using a natural cubic spline, compute the lower envelope by connecting all the local minima using another natural cubic spline, subtract the mean of the upper and lower envelopes, m11, from the signal, and then treat the residuary signal as a new signal. Repeat these steps for K times until the left signal has a pair of symmetric envelopes (i.e., m1K E0), and then define c1 as c1  um11    m1K

ð11Þ

This sifting process eliminates low-frequency riding waves, makes the wave profile symmetric, and separates the highest-frequency IMF from the current residuary signal. During the sifting process for each IMF, a deviation Dv is computed from two consecutive sifting results as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N N uX X ð12Þ Dv  t ½c1k ðt i Þc1k1 ðt i Þ2 = c21k1 ðt i Þ, i¼1

i¼1

where, for example, c1k u m11...  m1k, ti ¼(i  1)Dt, and T( ¼NDt) is the sampled period. A systematic method to end the iteration is to limit Dv to be a small number (e.g., Dv ¼0.001) and/or to limit the maximum number of iterations. However, the dyadic property of EMD is valid only if one kept the iteration number is fixed at about 10 because the dyadic property would break down if the iteration number is too high or too low [12]. After c1 is obtained, define the residual r1 (  u  c1), treat r1 as the new data, and repeat the steps shown in Eq. (11) to obtain other ci (i ¼2,...,n) as cn ¼ r n1 mn1    mnK ,

r n1  uðtÞc1    cn1

ð13Þ

The whole sifting process can be stopped when the residual rn ( ¼rn  1  cn) becomes a monotonic function from which no more IMFs can be extracted. For data with a trend, the last IMF rn should be the trend and has at most one local extremum. As shown later in Section 4 by numerical examples, such IMFs often have physical meanings, especially for nonlinear mechanical systems. The second step of HHT is to perform Hilbert transform and compute the time-varying frequency oi and amplitude Ai of each ci. After all ci(t) are extracted, one can perform Hilbert transform to obtain di(t) from each ci(t) and then define zi(t) and use Eq. (10) to obtain qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 z ðtÞ  c ðtÞ þ jd ðtÞ ¼ A ejyi , A  c2 þ d , y  tan1 d =c , o  y_ i

i

uðtÞ ¼ r n þRe

i

i

n X i¼1

i

i

i

! ½ci ðtÞ þjdi ðtÞ

¼ r n þ Re

i

n X

i

i

i

i

! Ai ðtÞejyi ðtÞ

ð14Þ

i¼1

Comparison of the time-varying Ai(t) in Eq. (14) with the constant Ai



qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2  2 a2i þbi in Eq. (4) reveals that HHT is an

adaptive signal decomposition method and it allows the use of distorted harmonics to decompose a complex signal into just a few IMFs. Because the moving average of each ci(t) is removed by EMD, the accuracy of instantaneous frequency oi from Eq. (14) is higher than those from Eqs. (8c) and (9). Because the so-obtained frequency is referred to the origin O1 on the phase plane (Fig. 2a), it is a referred IF (i.e., rIF) based on Eq. (7). Eq. (7) and (9) show that rIF and cIF are the same only if the signal is a regular harmonic. Moreover, u€ and v€ are involved in the calculation of cIF but not in the calculation of rIF. It follows from Eq. (14) that oi and other time derivatives can be computed using c c_ þ di d_ i , A_ i ¼ i i Ai

2 2 c_ 2 þ ci c€ i þ d_ i þ di d€ i A_ i A€ i ¼ i , Ai

oi ¼

ci d_ i c_ i di A2i

,

o_ i ¼

ci d€ i c€ i di 2oi ðci c_ i þ di d_ i Þ A2i

ð15Þ

The c_ i , d_ i , c€ i , and d€ i needed in Eq. (15) can be computed using the inverse Fourier transform shown in Eq. (6) without numerical differentiation in time domain. Some nonlinearity identification methods based on the use of the instantaneous frequency and amplitude and their derivatives shown in Eq. (15) have been proposed [30–32]. Unfortunately, our numerical

340

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

simulations show that this method is not able to provide accurate estimation of nonlinear dynamic characteristics because of the Gibbs’ effect and other numerical problems [33]. _ 1k and m € 1k are Because the m1k in Eq. (11) is the average of the upper and lower cubic spline envelopes, m1k, m _ and u€ remain in c1(t), and all other ci, c_ i , and c€ i are continuous time functions. In other words, the discontinuities of u, u, continuous everywhere. Hence the o1 and A1 of c1 become excellent indicators for pinpointing time instants of intermittent transient effects (e.g., impacting loads on structures and breathing of cracks). Moreover, if a highfrequency wave is added to the signal, the added wave and the signal’s discontinuities will be extracted as the first IMF. Then the discontinuities can be accurately extracted by subtracting the added wave from the first IMF, and the extracted discontinuities are more accurate than those from wavelet analysis. This property of HHT can be also used to filter noise [33]. Hence, HHT is more accurate and versatile than short-time Fourier transform and wavelet transform for time–frequency analysis of nonlinear non-stationary signals. Unfortunately, the accuracy of HHT analysis suffers from the end effect and several other mathematical and numerical problems. (#1) Because of the use of natural boundary conditions, EMD cannot provide accurate decomposition at the two data ends and hence introduces the first end effect that is later worsened by HT [34,35]. (#2) If ci(0)aci(T), HT aggravates the end effect introduced by EMD and decreases the accuracy of the extracted frequencies and amplitudes at data ends. (#3) Even if the two ends of an extracted IMF are smoothly extended to begin and end at zero, the calculated frequencies and amplitudes are inaccurate if the envelopes of the two extended segments are asymmetric. (#4) To accurately estimate, for example, the nonlinearity coefficient a in Eq. (1a), a3(t) needs to be separated from A(t) , which HHT cannot provide because a3 is usually small. (#5) HHT cannot accurately decompose a signal consisting of two non-synchronous regular harmonics with a frequency ratio 2/3o o1/o2 o1.5 when they reach their extrema around the same time [33,36], which causes the serious intermittency problem (i.e., a mode only appears in a segment of time) and decreases the credibility of HHT. (#6) If a signal contains two harmonics of close frequencies and amplitudes, the computed time-varying frequency and amplitude are inaccurate or erroneous because its envelope is asymmetric or its amplitude changes sign. (#7) If a transient signal contains a fast decaying (or increasing) moving average, the moving average alters the signal’s peaks and reduces the accuracy of extracted components because the computed lower and upper envelopes are inaccurate. (#8) If the sampling rate is not high enough (e.g., less than 10 samples per cycle), the extracted IMFs are inaccurate because local extrema (i.e., peaks and valleys) are not accurately located during EMD. Problems #1 to #7 can be reduced by pre-processing the signal before HHT analysis and post-processing the extracted IMFs by further decomposition [33]. Problem #1 can be reduced by using slope-based methods to extend the two data ends [34,35]. However, our simulations indicate that slope-based methods reduce the end effect for some signals but actually worsen the end effect for some other signals. Problems #2 and #3 can be relieved by adding a characteristic wave with a symmetric envelope to the two data ends and making the data to begin and end at zero, adding a mirror image, and/ or using a smooth near-square window. Problem #7 can be relieved by curve-fitting using low-order polynomials to remove the moving average before EMD analysis. Problems #4 to #6 can also be solved by extending the synchronous detection method (SDM) used in radio-signal demodulation and lock-in amplifiers to further decompose an IMF [33,37]. For time–frequency analysis without these complicated retrofitting techniques, we propose in the next section a conjugatepair decomposition (CPD) method that is more appropriate for signals from low-frequency sampling and can be used for online frequency tracking. 3.3. Conjugate-pair decomposition We propose an enhanced conjugate-pair decomposition (CPD) method that improves a three-point method for tracking of instantaneous frequencies of an arbitrary signal [13]. To enable frequency extraction of a signal u(t) containing a moving average C0, we assume the signal to have the following form: uðtÞ ¼ C 0 þ

p X

p X       ai cos oi t þ fi ¼ C 0 þ ½C i cos oi t Di sin oi t 

i¼1

ð16aÞ

i¼1

    C i  ai cos oi t n þ fi , Di  ai sin oi t n þ fi

ð16bÞ

where oi are assumed to be the frequencies estimated at the previous time step tn  1, t ð  tt n Þ is a shifted time, and tn is the current time instant. The unknown functions C0, Ci, and Di for the current data point at t ¼ 0ði:e:,t ¼ t n Þ can be determined by minimizing the following squared error: Serror 

ðm1 XÞ=2

 2 U n þ k un þ k

ð17Þ

k ¼ ðm1Þ=2

where m (Z1 þ2p) is the total number of processed data points and is assumed to be an odd number here, un þ k denotes the function form shown in Eq. (16a) with t  t k ¼ kDt, and Un þ k denotes the actual u(t) at t¼tn þkDt. After Ci and Di are determined, the instantaneous amplitude Ai, phase angle yi, and frequency oi can be obtained as qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   ð18aÞ Ai ¼ C 2i þD2i , yi  oi t n þ fi ¼ tan1 Di =C i

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

oi  y_ i 

yi ðt n Þyi ðtn qDtÞ qDt

341

ð18bÞ

The backward finite difference is used to average the frequency over qDt in order to reduce the influence of noise in frequency tracking. For post-processing instead of online tracking, a central finite difference scheme should be used to have better accuracy. One can see from Eqs. (16b) and (18a) that Di is the Hilbert transform of Ci. This frequency tracking method is essentially based on circle-fitting of data on the u–v plane to determine C0, Ai, and yi [13,24]. Hence, the so-obtained frequency is a circular IF (i.e., cIF) w.r.t. the center O in Fig. 2a and is based on Eq. (9). To initiate this frequency tracking process starting from t ¼t1 þ 2p, one can use the Teager–Kaiser algorithm (TKA) to provide the initial guess of oi for Eq. (16a) [13,38] or use the frequency spectrum of a small section of the beginning time trace. Note that, because the oi are assumed to be known constants in Eq. (16a), the algebraic equations for least-squares fitting are linear ones, instead of nonlinear ones when the oi are treated as unknowns. Because there are only 1þ2p unknowns in Eq. (16a), only 1 þ2p local data points (i.e., m ¼1þ2p) need to be processed in Eq. (17). Noise filtering is possible if more than 1þ2p points are used in Eq. (17). However, although more data points can be used to further reduce the influence of measurement noise and/or signal processing errors, the estimated frequency is more a locally averaged value. For any point other than the starting point at t ¼t1 þ 2p , the previous point’s frequency can be used as an initial guess of frequency without using TKA. Because the initial guess of oi is updated using Eq. (18b), this method is able to adapt and capture the actual variation of frequency [13]. Moreover, the accuracy of oi for each time instant can be improved by running several iterations of the circle-fitting process and updating the value of oi. Instead of using just C0, one may consider a linear polynomial or even a high-order polynomial in Eq. (16a) for separating any possible linear or high-order transient parts. However, because CPD depends on local least-squares curve-fitting using just a few data points, the orthogonality between the functions used in Eq. (16a) cannot be accurately enforced. When more unknowns are used in Eq. (16a), the estimations become less accurate. Hence, Eq. (16a) is recommended. Moreover, if an IMF from EMD is to be further analyzed by using Eq. (16a), one can set C0 ¼0 to accelerate the computation and increase the accuracy because any IMF from EMD analysis should have no moving average. Accurate frequency tracking is the key for online identification and health monitoring of dynamical systems. For a long stationary signal, one can easily figure out its average frequency and amplitude over one period. But, if only four (or three) consecutive data points covering less than one period of the signal are given, how to estimate the instantaneous frequency and amplitude is challenging and is the so-called frequency tracking problem. Because HHT requires a certain length of data and the use of HT, it cannot be used for frequency tracking because a set of three or four data points is too short to have maxima and minima and to perform HT. To the author’s knowledge, TKA is the only method in the literature that can process only four data points for online frequency tracking without any other help [38]. Because the instantaneous frequency from TKA is also based on Eq. (9), it is a cIF [13]. Unfortunately, the accuracy of TKA is easily destroyed by measurement noise due to its use of finite difference [13]. Moreover, because a signal is assumed to be a pure harmonic in TKA, any moving average in the signal can destroy the accuracy of TKA. On the other hand, CPD uses a constant and one (or more) pair of adaptive harmonics to fit three (or more) recent data points and estimate the instantaneous frequency and amplitude, and it dramatically reduces the influence of any moving average. Moreover, noise filtering is an implicit capability of CPD because of its use of curve-fitting, especially when the number of processed points increases. Furthermore, different from HHT, CPD requires the processed m data points to cover at least ¼ period of the signal and hence is more appropriate for high-frequency signals sampled at low frequencies. Hence, CPD and HHT can actually complement each other for more accurate adaptive time–frequency analysis for system identification. To use the CPD shown in Eqs. (16a)–(18b) for signal decomposition, data from both sides of the point under evaluation need to be used in Eq. (17) and (18b). For online frequency tracking, data from the left side (previous time instants) need to be used. For a highly nonlinear system, the recommended approach is to use p ¼1 in Eq. (16a) and then find out the order of nonlinearity from the obtained time-varying amplitude and frequency. After that, use p¼ 2 in Eq. (16a) with o2 ¼3o1 (if the nonlinearity is identified to be cubic) or o2 ¼2o1 (if the nonlinearity is quadratic). Then, use Eq. (1a)–(1h) to estimate system parameters. For a signal consisting of several regular/distorted harmonics, the recommended approach is to use EMD to decompose it into IMFs and then use CPD to perform time–frequency analysis on each IMF. 4. Numerical and experimental results 4.1. Extraction of instantaneous frequency and amplitude To demonstrate CPD and HHT for extracting instantaneous frequencies and amplitudes we consider the following signal: uðtÞ ¼ e0:03t ð1 þ 0:2sinoa t Þsinðot þ sinop tÞ þ 0:01  noise,

o ¼ 2p, oa ¼ 0:1o, op ¼ 0:2o

ð19Þ

where Dt¼1/16 s is used and noise is a vector of normally distributed zero-mean random numbers with a standard deviation of one. One may examine the signal’s power spectral density (PSD, Fig. 3b) and follow Eq. (4) to explain that the signal consists of many regular harmonics, but it is almost impossible to extract quantitative or even just qualitative information about

342

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

Fig. 3. HHT and CPD analyses of Eq. (19): (a) signal, (b) PSD, (c and d) from HHT analysis, and (e and f) from CPD analysis. Blue lines are exact solutions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

nonlinearity. Fig. 3c compares the original time-varying amplitude and frequency (thin blue lines) with those extracted by HHT and CPD with p¼1 in Eq. (16a) and m¼5 in Eq. (17). They show that results from CPD are slightly better than those from HHT. Moreover, HHT suffers from the end effect, but CPD does not. However, HHT and CPD have different types of errors, and the accuracy of HHT can be improved by using a smaller sampling interval Dt because HHT depends on accurate capture of the signal’s local peaks and valleys. On the other hand, the accuracy of CPD can be improved by using a larger Dt to increase the range covered by 4Dt (since m¼ 5 is used in Eq. (17)). In other words, when Dt is fixed, HHT is good for processing low-frequency signals and CPD is good for high-frequency signals. Hence, HHT and CPD can complement each other. One can also improve the time–frequency analysis by extracting the noise from the signal before HHT or CPD analysis. As explained right after Eq. (15), one can add a high-frequency harmonic (e.g., 2sinð8ptÞ) and then extract the added harmonic with noise as the first IMF to make the signal free from noise. Because HHT processes the whole set of data to obtain the time–frequency decomposition and because the end effect exists, HHT is not capable of online frequency tracking. On the other hand, CPD processes only few (three or more) recent data points and hence is very appropriate for online frequency tracking [13]. However, to use CPD for tracking a lowfrequency signal, down-sampling may be necessary in order to have at least ¼ period of the component to be tracked covered by the few recent data points.

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

343

4.2. Identification of transient events To show HHT and CPD for capturing multiple transient events we consider the following nonlinear non-stationary system:    _ ¼0 m u€ þ 2Bou_ þ r o2 u þ 0:05u3 ¼ F 0 , uð0 Þ ¼ 2, uð0Þ 8 t5 > > for 0 rt r 7:5 0:980:02tanh > > > tr > > > > t10 > > > < 0:940:02tanh t for 7:5o t r12:5 r ð20Þ r¼ t15 > > for 12:5 ot r 17:5 > 0:900:02tanh > > tr > > > > t20 > > > : 0:860:02tanh t for 17:5 ot r

o ¼ 2p, B ¼ 0:01, m ¼ 2, tr ¼ 0:05 where tanh½ðtt a Þ=t r  is a function that changes from  1 to þ1 around t¼ ta with a transition period determined by tr. In addition, there is an impact force F0 ¼5 N applied between 25 rtr25.05 s. Dt ¼0.1 s is used and a vector of normally distributed zero-mean noise with a small standard deviation of 0.003 is added to the signal obtained from direct numerical integration using the Runge–Kutta method. The PSD in Fig. 4b shows a wide frequency bandwidth, which is caused by the multiple transient events. Fig. 4c shows that the frequency tracked by CPD clearly shows the impact instant at t¼25 s and the four time instants when the stiffness suddenly decreases by 4% and the natural frequency drops by 2% because the pffiffiffi vibration frequency is proportional to r . Fig. 4c also shows that the free vibration frequency decreases when the amplitude decreases, indicating a hardening–type nonlinearity. Fig. 4d shows that the frequency extracted by HHT does not clearly show the five transient events (especially the impact event) because the low sampling rate, noise, and the end effect introduced by HT decrease the accuracy. On the other hand, a low sampling rate is favorable to CPD, CPD averages out noise because it is based on curve-fitting, and the end effect does not exist in CPD because HT is not used. This example shows that CPD is better than HHT for capturing transient events. The accuracy of that shown in Fig. 4d from HHT can be improved by using a smaller Dt, but, for this example, the results can never be as good as those from CPD. 4.3. Parametric identification First we use the following nonlinear oscillator to show that at least one free transient response and one forced steadystate response are needed for parametric identification:   m u€ þ 2Bou_ þ o2 u þ au3 ¼ F 0 cos Ot ð21aÞ m ¼ 2,

_ o ¼ 2p, B ¼ 0:01, a ¼ 10,F 0 ¼ 30, O ¼ p, uð0Þ ¼ 0, uð0Þ ¼ 20

ð21bÞ

Dt ¼0.1s is used here. Fig. 5b shows the result from EMD analysis. It is obvious that c1 is the damped natural harmonic and c2 is the distorted harmonic caused by the harmonic excitation with O ¼0.5 Hz and the cubic nonlinearity. For a linear system, the particular solution under a harmonic excitation should be independent of the transient response and is a steady-state response with a constant amplitude starting from t ¼0. The nonlinear term couples the forced vibration with the transient vibration, and the input excitation energy is partially used to keep the transient response from decaying too fast. This explains why it often takes a long time for a nonlinear system to achieve a steady-state forced vibration. This coupling also changes the amplitude of c2 from a constant in a linear case to a gradually increasing value in this nonlinear case. Hence, one cannot use the separated transient response c1 to estimate the damping ratio B because c1 contains some input energy from the applied force. If it is used, the estimated damping ratio from c1 would be too low at the beginning. Moreover, the separated c2 cannot be used to estimate the mass m as discussed in Section 2.1 because it is not a steady-state vibration. These are also the main reasons why some researchers (e.g., Refs. [31,32]) failed to obtain accurate nonlinearity identification by using the time-varying frequencies and amplitudes and Eq. (15). Hence, to obtain accurate nonlinearity identification, it is necessary to separately process one transient response and one steady-state response, instead of one response containing both transient and steady-state responses. Another important point here is that it is more accurate to estimate nonlinearities using a large-amplitude transient response than a steady-state response because the latter depends on the use of second-order asymptotic solutions. However, it is easier to detect nonlinearities using a small-amplitude steady-state response than using a small-amplitude transient response. This concept is very useful for early detection of nonlinear structural vibrations (e.g., flutter of aircraft, and chattering of machining) before the vibration is too big to be controlled. To show the use of time–frequency analysis for nonlinearity identification, we consider the nonlinear oscillator shown in Eq. (21a) with m ¼ 2,

o ¼ 2p, B ¼ 0:01, a ¼ 10, F 0 ð ¼ FmÞ ¼ 50, O ¼ o

ð22Þ

344

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

Fig. 4. Time-varying frequency of the oscillator shown in Eq. (20): (a) u(t), (b) PSD, (c) from 5-point CPD, and (d) from HHT.

c1

c2

Fig. 5. EMD analysis of Eq. (21a,b): (a) u(t), and (b) ci(t).

Direct numerical integration is performed with Dt¼1/30 s for 150 s to damp out any transient effect. Fig. 6b shows that the steady-state response consists of several harmonics equally separated by 2 O(¼ 2 Hz) , indicating a frequency modulated at 2 O. Fig. 6c and d show that the time-varying frequency and amplitude from HHT analysis modulate at 2 O, indicating a cubic nonlinearity (see Eqs. (1c) and (1e)). Moreover, because the frequency and amplitude are at their maxima when u is at its local maxima or minima and are around their minima when u ¼0, we know that a 4 0. If only one pair of conjugate functions are used in Eq. (16a), Fig. 6c and d show that the time-varying amplitude and frequency from CPD do modulate at 2 O and the amplitude is at its maxima when u is at its local extrema, indicating a hardening cubic nonlinearity. However, the time-varying frequency does not clearly indicate a hardening nonlinearity and the variations of amplitude and frequency are not the same as those from HHT analysis because they are different by definition (see Eq. (7) and (9)). Unfortunately, EMD is not able to decompose the time-varying amplitude shown in Fig. 6c into a1 and a3 that are needed for using perturbation solutions to estimate system parameters. On the other hand, CPD is able to. Fig. 6e and f show that the amplitudes and frequencies from CPD using two pairs of conjugate functions with o1 ¼ O and o2 ¼3 O in Eq. (16a) are constant and a1 ¼1.4834 and a3 ¼0.0286. Note that a1 and a3 are equal to the average and varying amplitudes in Fig. 6c from HHT. Because there are five unknowns, at least five data points are needed for frequency extraction and

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

CPD

HHT

345

CPD

HHT

a1

ω3 = 3ω1 ω1

a3

Fig. 6. Analysis of the steady-state response of Eqs. (21a) and (22): (a) u(t), (b) PSD, (c) amplitude, (d) frequency, and (e and f) from CPD with p ¼ 2 in Eq. (16a).

seven points are used here. The transient effect at the beginning of Fig. 6e and f is caused by the use of TKA to estimate the frequency at ts ¼3Dt. One can use Eq. (1b) to estimate a ¼ 32O2 a3 =a31 ¼ 11:07. Furthermore, we obtain f ¼0.0440 from Fig. 6a and Eqs. (1a) and (1b) and B ¼ ½3aa21 =ð8oÞstan f=o ¼0.0102 from Eq. (1g) with s( ¼ O  o)¼0. Then Eq. (1f) with s ¼0 is used to obtain F ¼27.127 ( ¼F0/m) and hence m¼F0/F ¼1.843 is estimated. The key point here is that all a, B and m can be estimated with an error less than 10% from just one steady-state response to a harmonic excitation at O E o. ^ ¼ 1:21o, which It follows from Eq. (1h) that the undamped natural frequency under a vibration amplitude a1 ¼1.4834 is o indicates that this is not really a weakly nonlinear problem. Because the perturbation solutions are derived based on the assumption of weak nonlinearity, the estimation accuracy increases when the nonlinearity coefficient a and the excitation force F0 decreases. Also, the accuracy of estimations can be further improved when more sets of responses to different F0 and/or O are available. Fig. 6e and f confirm again that the use of CPD for further decomposition of an IMF is better than the use of the synchronous detection method used in radio-signal demodulation and lock-in amplifiers [33]. Next we consider a transient response of the nonlinear oscillator shown in (21a) with m ¼ 2,

_ o ¼ 2p, B ¼ 0:01, a ¼ 10, F 0 ¼ 0, uð0Þ ¼ 3, uð0Þ ¼0

ð23Þ

because Eqs. (1b and 1e) and Fig. 6d show that the referred instantaneous frequency (rIF) from HHT can be directly used for nonlinearity identification but the circular instantaneous frequency (cIF) from CPD using p ¼1 in Eq. (16a) cannot,

346

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

200ς from  u(t)

200ς from u(t)

min

max

Fig. 7. HHT analysis of Eqs. (21a) and (23): (a) o1, (b) A1&B, (c) o1 , and (d) o1  A1. (For interpretation of the references to color in this figure, the reader is referred to the web version of this article.)

we will use HHT for this numerical case. The only IMF of this transient vibration is c1( ¼u), i.e., the damped natural harmonic. Fig. 7a and b show that, when the vibration amplitude A1 is large, o1 exactly modulates at 2o1 and decreases with A1, indicating the existence of hardening cubic nonlinearity. The large variations around the two ends in Fig. 7a and b are caused by the discontinuity c1(0)ac1(30) and Gibbs’ effect, but the computed o1 and A1 corresponding to peaks (black circles in Fig. 7c) and valleys (blue crosses in Fig. 7c) of c1 are still accurate, as indicated by the two broken tangent lines in Fig. 7a and b. The local mimina in Fig. 7c indicate that the computed o1 and A1 corresponding to peaks and valleys of c_ 1 are _ inaccurate because the o1 and A1 here are computed from u(t) but u(t) E0 at these time instants. If the velocity signal uðtÞ is processed, the computed o1 and A1 corresponding to peaks and valleys of c_ 1 are accurate, but the computed o1 and A1 corresponding to peaks and valleys of c1 will be inaccurate. The black broken line in Fig. 7b is the damping ratio calculated using finite difference and the o1 and A1 of u(t), which is inaccurate because the correct answer should be 200B ¼2. On the _ other hand, the damping ratio calculated using the o1 and A1 of uðtÞ is more accurate, as shown by the blue solid line in _ for signal processing, and it is better to use uðtÞ _ Fig. 7b. Hence, it is better to have both u(t) and uðtÞ for the estimation of the damping ratio and damping force and u(t) for the estimation of the spring force and backbone curve. If only the velocity _ uðtÞ is directly measured (e.g., using a laser vibrometer), u(t) can be easily obtained by numerical integration and the integration constant due to the unknown initial displacement will be automatically removed by the EMD process in HHT _ analysis. If only the displacement u(t) is directly measured, it is better to use just u(t) because uðtÞ from numerical differentiation is often inaccurate, especially if measurement noise is significant. In Fig. 7d, when A1 is large, the backbone curve from the perturbation solution (o1 ¼ o þ3aA21 =8o, Eq. (1h)) deviates from the actual backbone curve (the dots) because A1 o1 is assumed in the perturbation analysis. Because Eq. (1h) does not include the influence of a3, the match can be improved by using CPD to separate a1 from a3 like that shown in Fig. 6e. Results shown in Figs. 6 and 7 clearly demonstrate that HHT and CPD complement each other for signal processing and system identification.

4.4. Non-parametric identification of nonlinearity The time-varying frequency and amplitude shown in Fig. 7a and b can be used with Eq. (2b) to obtain the nonlinear specific spring force Fs shown Fig. 8a, which shows a small deviation from the exact one (Fs ¼ o2u þ a3u3) because the amplitude and frequency of c1 modulate at 2o1. However, when the signal is further decomposed into two components

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

347

Fig. 8. Non-parametric identification of Eqs. (21a) and (23): (a) u Fs curve, and (b) u Fd curve. The solid lines in (a) and (b) represent the exact spring and damping forces, respectively.

by using CPD with p ¼2 in Eq. (16a) and frequencies o1 and 3o1 (o1 ¼the sliding-window average of the o1 in Fig. 7a), the estimated Fs agrees well with the exact curve. One can use this u Fs curve to estimate the type and order of nonlinearity. The specific damping force Fd shown Fig. 8b is obtained using the time-varying frequency and amplitude from HHT _ analysis of uðtÞ and the damping ratio calculated and presented in Fig. 7b (blue solid line). It deviates from the exact one ðF d ¼ 2Bou_ Þ, but it shows the correct, linear damping. The deviation is mainly caused by the use of finite difference and the linear formula to calculate B (see Fig. 7b). Again, the match can be further improved by using CPD with p ¼2 in Eq. (16a) and frequencies o1 and 3o1 . Fortunately, such further decomposition process is only needed for highly nonlinear _ d shown in Fig. 8a and b should be accurate enough for problems. For common weakly nonlinear problems, the u-Fs and u-F non-parametric identification. 4.5. Identification of nonlinear continuous systems Here we consider a 1000  25.4  6.4 mm3( ¼L  b  h along the x, y, and z axes) 2024-T851 aluminum beam hinged at x ¼0.2L & 0.8L and subjected to two harmonic base excitations wð0,t Þ ¼ wðL,t Þ ¼ 3:0sinðOtÞ mm along the z direction. The material properties are Young’s modulus E ¼72.4 GPa, Poisson’s ratio n ¼0.33, and mass density r ¼2780 kg/m3. The first four linear natural frequencies oi are obtained to be 32.91, 72.15, 110.77, and 221.10 Hz. Transient analysis by direct numerical integration using the Newmark-beta method is performed using an in-house fully nonlinear finite-element code GESA (Geometrically Exact Structural Analysis) [6]. GESA has been implemented with geometrically exact 1D and 2D structural theories for modeling and analysis of highly flexible structures undergoing large elastic deformations, and its accuracy has been verified through many numerical and experimental studies [6]. The beam here is modeled using twenty fully nonlinear beam elements (eBeam29n in Ref. [6]), and the damping matrix is computed using a modal damping ratio of 0.03 for each mode. If the excitation frequency is O ¼32.91 Hz ( ¼ o1 ), Fig. 9 shows the steady-state response w(0.5L,t) at Node 11 and its time-varying frequency and amplitude from HHT analysis. The only IMF extracted by EMD is the signal itself. Because Fig. 9a, c and d show that o1 and A1 modulate exactly at 2 O and o1 and A1 are at their maxima when w(0.5L,t) is at its maxima or minima, it clearly reveals that the beam has a hardening cubic nonlinearity. Fig. 9d shows that the continuous frequency band of w(0.5L,t) is 30.3–35.6 Hz but that of w(0.1L,t) is 32.5–33.3 Hz. Fig. 9c and d also show that w(0.1L,t) and w(0.35L,t) have about the same amplitude but very different frequency bandwidths. On the other hand, w(0.35L,t) and w(0.5L,t) have about the same frequency bandwidth but very different amplitudes. The reason is that, because the beam is hinged at x ¼0.2L & 0.8L, w(0.35L,t) and w(0.5L,t) are under the stretching effect caused by geometric nonlinearity but w(0.1L,t) is not. Hence, the frequency variation of a signal is valuable for detection of nonlinearity even under small-amplitude vibration. For regular on-earth structures, a common understanding is that the amplitude of a linear vibration needs to be less the structure’s thickness. Fig. 9d shows that, although the vibration amplitude of w(0.5L,t) is less than h/2, the frequency variation is large ( 415%). Hence, it may be better to define linear vibration by giving a limit on the continuous frequency bandwidth, instead of the vibration amplitude. When CPD with p ¼2 in Eq. (16a) is used to decompose the w(0.5L,t) shown in Fig. 9a into two components, we obtain a1 ¼2.153 mm and a3 ¼0.0878 mm. Using Eq. (1b) we obtain a ¼ 32O2 a3 =a31 ¼ 0:2815O2 ¼ 12037:3 (with w(x,t) in mm), which is not a small number. One can use this nonlinear coefficient in the study of idealized single-mode nonlinear vibration. Experimental results presented next show that the proposed time–frequency analysis is accurate for delicate identification of nonlinearities. 4.6. Experimental results Here we present three real examples to demonstrate applications of the proposed time–frequency analysis for system identification. They are the head motion of a small electromagnetic shaker, the free vibration of a clamped-free beam, and the free vibration of a clamped–clamped beam. We do not use hinged–hinged beams (see Fig. 9b) because accurate hinge

348

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

0.5L

0.1L

0.5L

0.35L 0.1L

0.35L

Fig. 9. HHT analysis of response of a hinged-hinged beam: (a) w(0.5L,t), (b) 25 ODSs during one period, and (c,d) A(t) and o(t) of w(0.1L,t), w(0.35L,t) and w(0.5L,t).

supports for large dynamic loading under nonlinear vibration are difficult to design and manufacture. Moreover, misalignment of hinges and contact friction cause dynamic disturbances. One may consider the use of ball bearings or roller bearings for end supports, but the non-smooth contact and clearance in such bearings cause serious dynamic disturbances. Furthermore, if shakers are used to obtain forced vibrations, the dynamical shaker-structure interaction can distort measured vibration signals. Also, the shaker itself can distort the input signal from a function generator, as shown next in Example #1. Hence, it is better to use free vibrations in order to have accurate system identification results.

4.6.1. Example #1: shaker head For structural vibration testing using electromagnetic shakers, it is important for the shaker not to distort the input signal in order to have high-fidelity experimental results. If the shaker head’s motion is too much distorted from the needed regular harmonic, it is difficult to correct the motion even by feedback control using an accelerometer for monitoring. Here we consider the no-load head motion of a K2007E0 miniature shaker (Fig. 10a) from the Modal Shop, Inc. When the shaker is fed with a regular harmonic signal from a Wavetek FG2A functional generator, its head’s velocity is measured using a Polytec PSV-200 scanning laser vibrometer with a sampling frequency of 2560 Hz (i.e.,Dt ¼1/2560 s). When the shaker is driven by a 10-Hz regular harmonic function from the function generator, Fig. 10b shows that the shaker head’s velocity is a distorted harmonic, and Fig. 10c shows that the continuous frequency bandwidth is about 6.0 Hz (i.e., 6.5–12.5 Hz), which is 60% of the input frequency. If the acceleration (obtained by finite difference) is analyzed, the frequency bandwidth is even larger. However, if the displacement (obtained by integrating velocity) is analyzed, the frequency bandwidth is reduced to 3.2 Hz (8.5–11.7 Hz). This signal distortion phenomenon often happens to electromagnetic shakers under low-frequency excitation. The obtained time-varying amplitude and frequency are clean because noise is filtered out by adding a high-frequency signal 0:5sin 100pt to the signal and then extracting it with the noise contained in the signal as the first IMF (see explanation right after Eq. (15)). The frequency is at its maximum (points 1 & 5 in Fig. 10c) when the velocity is at its maximum (i.e., zero displacement). The frequency is at its minimum (point 2) when the velocity is at its zero (i.e., max. displacement). Moreover, when the velocity decreases from zero at point 2 to point 3 in Fig. 10b, the frequency increases. Hence, we conclude that the signal contains softening cubic nonlinearity. However, when the velocity is right after its local minimum (point 3), the frequency tends to be constant for a while and then increases, indicating a hardening nonlinearity. But, when the frequency increases from point 4 to point 5, the displacement decreases (because the velocity increases), indicating a softening nonlinearity. Hence, the time–frequency analysis reveals that the shaker head mechanism contains hardening (point 3 to point 4) and softening (all other steps) nonlinearities. The peak-topeak displacement amplitude can be calculated to be 3.8 mm, and the shaker’s stroke length is 13 mm. The observed asymmetric nonlinearity is probably caused by the armature’s supporting structure.

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

4

1

3 5,1

349

2

5 5

1 4 3

2 4

2

3

5

1

2 4 2 4

1 3

5

3

Fig. 10. Motion of a shaker head: (a) shaker, (b and c) measured velocity and extracted A(t) and o(t) when O ¼ 10 Hz, and (d and e) measured velocity and extracted A(t) and (t) when O ¼94 Hz.

When the shaker is driven by a 94-Hz regular harmonic function, Fig. 10d and e show that the amplitude is about constant and the continuous frequency bandwidth is about 7 Hz, which is 7.4% of the input frequency. The extrema and zeros of the signal shown in Fig. 10d and the corresponding extrema of the frequency shown in Fig. 10e indicate the existence of a hardening cubic nonlinearity because the amplitude and frequency modulate at 2 O and the frequency is at its maximum (points 2 & 4) when the displacement is at its maximum (i.e., zero velocity). The peak-to-peak displacement amplitude can be calculated to be 0.8 mm, which is much smaller than that of the 10-Hz vibration. This is probably the main reason why the frequency bandwidth and hence signal distortion decrease when the vibration frequency increases, as shown by Fig. 10b and d.

4.6.2. Example #2: clamped-free beam We consider a horizontal clamped-free 784.2  31.75  3.175 mm3 steel beam subject to an initial static transverse tip displacement about 11.5 cm, as shown in Fig. 11a (without the stopper). Fig. 11b shows the transverse vibration velocity of a point at x ¼158.7 mm measured using a Polytec PSV-200 scanning laser vibrometer with a sampling frequency of 2560 Hz (i.e., Dt ¼1/2560 s). The first four natural frequencies oi are identified from the signal’s power spectral density (PSD, Fig. 12c) to be 3.82, 24.1, 67.5, and 132.0 Hz, and the mass density r is measured to be r ¼6951 kg/m3. Hence, Young’s modulus E is estimated by matching the first four experimental natural frequencies with the theoretical linear frequencies (3.85, 24.12, 67.53, 132.34 Hz) from the Euler–Bernoulli beam theory to be E¼148 GPa.

350

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

y

x

z

784.2mm×31.75mm×3.175mm

158.7mm

w(x,t)

stopper

x = 611.7 mm

PSV-200 laser vibrometer

_ Fig. 11. Transient vibration of a clamped-free beam: (a) setup (no stopper), (b) w(0.1587 m,t), and (c) IMFs from EMD.

Fig. 11c shows the three IMFs sequentially extracted using EMD and a signal conditioning method. We add 0:6sinðo3 t Þ to the signal and then extract the third and other high-frequency modes and noise into c1 and then obtain the second and first modal vibrations as c2 and c3, respectively. Fig. 12a and b show the time-varying amplitudes and frequencies of IMFs from HHT analysis. Because the first IMF c1 (mainly the third modal vibration) contains other high-frequency modes and noise of the original signal, its calculated frequency o1 ð  o3 Þ shows significant variation in Fig. 12b. Moreover, Fig. 12b shows that, when an IMF’s amplitude is small, its calculated frequency is inaccurate because a small ci in Eq. (14) makes   tan1 di =ci difficult to be accurate. For HHT, the inaccurate amplitudes and frequencies at the two data ends in Fig. 12a and b are caused by the discontinuity-induced Gibbs’ effect. The spectrum of c1 (not included here) clearly reveals that the 3rd and other high-frequency modes are all extracted into c1. Fig. 12a and d show the time-varying amplitude (broken black line) and frequency of the first-mode vibration obtained by using CPD to process the c3 in Fig. 11c after downsampling using Dt~ ¼ 40Dt and five points (4Dt~ cover about one period of c3). Fig. 12d clearly shows that the frequency increases about 0.02 Hz when the amplitude decreases during 6.4 s, and the frequency varies at a frequency of 2o1 . These two phenomena reveal that the first-mode vibration c3 contains a softening cubic nonlinearity, which agrees with previous numerical results from multiple shooting analysis [6]. Moreover, Figs. 11c and 12d show that the frequency is at its maximum when c3 is at its maximum or minimum (i.e., locations of zero displacement), and the frequency is at its minimum when c3 is at its zeros (i.e., locations of maximum or minimum displacement). This also indicates the existence of softening nonlinearity. This example demonstrates how delicate nonlinearity identification can be performed by using the proposed time–frequency analysis that combines HHT and CPD analyses. 4.6.3. Example #3: clamped–clamped beam We consider the beam shown in Fig. 11a (without the stopper) with clamped–clamped boundary conditions and a length L¼859.2 mm. The first four natural frequencies oi are identified from the signal’s power spectral density (PSD) to be 19.5, 54.5, 107.7, and 179.9 Hz, and the theoretical linear frequencies from the Euler-Bernoulli beam theory are 20.4, 56.2, 110.2, and 182.3 Hz. The mid-span point of the beam is given an initial static displacement about 2 mm, and a sampling frequency of 5120 Hz (i.e., Dt ¼1/5120 s) is used. Fig. 13a and b show the transverse vibration velocity measured at x ¼L/2 and the two IMFs extracted from the signal using EMD with the signal conditioning method by adding 0:6sinðo3 t Þ to the signal before

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

351

ω1

A3

ω2

A2 A1

ω3

ω3

ω2 ω1

ω3

_ _ Fig. 12. Time–frequency analysis of w(0.1587 m,t) : (a and b) amplitudes and frequencies of IMFs from HHT, (c) PSD of w(0.1587 m,t), and (d) frequency of c3(t) from CPD. The black broken line in Fig. 12a is from CPD analysis.

_ Fig. 13. Transient vibration of a clamped–clamped beam: (a) w(0.5L,t), and (b) IMFs from EMD.

extracting the first IMF. Because it is a symmetric loading, the second mode is not excited and the c1 and c2 are the third and first modal vibrations, respectively. Figs. 13b and 14a and b show that the amplitude and frequency of c2 (i.e., the first mode) modulate at 2o1 and they are at their maximum when c2 is at its zero (i.e., locations of maximum or minimum

352

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

Fig. 14. Time–frequency analysis of c2(t) from Fig. 13b: (a and b) amplitude and frequency from HHT, and (c and d) amplitude and frequency from CPD.

displacement) and at their minimum when c2 is at its maximum or minimum (i.e., locations of zero displacement). These clearly reveal the existence of a hardening cubic nonlinearity. Because the frequency from HHT analysis suffers from the end effect, it cannot show a clear trend of frequency decrease when the amplitude decays. However, Fig. 14c and d show the time-varying amplitude and frequency of the first-mode obtained by using CPD to process c2 in Fig. 13b after downsampling using Dt~ ¼ 16Dt and five points to cover about one period of c2. Fig. 14d clearly shows that the frequency decreases about 0.1 Hz when the amplitude decreases during 3.2 s, and the frequency varies at a frequency 2o1 . These two phenomena again reveal that the first-mode vibration contains a hardening cubic nonlinearity. The maximum vibration amplitude can be calculated from Fig. 14c and d to be as small as 1.8 mm ( oh¼ 3.175 mm), but the proposed method can still identify the hardening cubic nonlinearity. Because the proposed system identification method is based on time–frequency analysis of pointwise time responses, it works for any signals measured from 1D or 2D structures or 3D solids. Of course, more experimental verifications on different mechanical systems with different nonlinearities are needed in order to validate this technology. 5. Concluding remarks This work presents a methodology for nonlinear system identification using perturbation solutions and two algorithms that are complement to each other for adaptive time–frequency analysis. Hilbert–Huang transform (HHT) can be used for offline signal decomposition using the empirical mode decomposition (EMD) and extraction of the referred instantaneous frequency and amplitude of each extracted intrinsic mode function (IMF) using Hilbert transform (HT). The conjugate-pair decomposition (CPD) algorithm can be used for online tracking and offline post-processing of the circular instantaneous frequency and amplitude of an arbitrary signal and for further decomposition of IMFs in order to use perturbation solutions for accurate parameter estimation. Perturbation analysis provides guidance and approximate analytical solutions for system identification using each IMF’s time-varying frequency and amplitude. CPD can also use EMD with signal conditioning techniques to avoid the intermittence problem and then use one or more pairs of conjugate functions to extract the time-varying frequency and amplitude of each IMF. CPD has an implicit noise filtering capability, and this capability can be enhanced by processing more neighboring data points for each estimation. Furthermore, CPD is better than HHT for capturing transient events because CPD is not affected by the end effect caused by signal discontinuity. However, CPD needs the Teager–Kaiser algorithm to estimate the frequency of the first point in order to start the tracking process, and it is accurate only if the used local data points (three or more points) can cover an enough portion (¼ or more) of the period of the component to be extracted. Simulated dynamic responses of discrete and continuous systems and

P. Frank Pai / Mechanical Systems and Signal Processing 36 (2013) 332–353

353

experimental nonlinear vibrations of a shaker head and clamped–free and clamped–clamped beams are used to verify the proposed methodology. Numerical and experimental results show that the proposed methodology is capable of accurate parametric and non-parametric identification of nonlinear dynamical systems. This work also proposes to define the so-called linear vibration by giving a limit on the continuous frequency bandwidth of a system’s response to a harmonic excitation, instead of a limit on the vibration amplitude.

Acknowledgment This work is partially supported by the Office of Naval Research under Grant no. 0014-11-1-0334 and the National Science Foundation under Grants CMMI-1039433 and IIP-1026883. The support is gratefully acknowledged. References [1] A.H. Nayfeh, D.T. Mook, Nonlinear Oscillations, Wiley-Interscience, New York, 1979. [2] F.C. Moon, Chaotic Vibrations, An Introduction for Applied Scientists and Engineers, Wiley-Interscience, New York, 1987. [3] K. Worden, G.R. Tomlinson, Nonlinearity in Structural Dynamics: Detection, Identification and Modelling, Institute of Physics Publishing, Bristol and Philadelphia, 2001. [4] C.H.M. Jenkins, Gossamer Spacecraft: Membrane and Inflatable Structures Technology for Space Applications, AIAA, Reston, Virginia, 2001. [5] A.H. Nayfeh, P.F. Pai, Linear and Nonlinear Structural Mechanics, Wiley-Interscience, New York, 2004. [6] P.F. Pai, Highly Flexible Structures: Modeling, Computation and Experimentation, AIAA, Reston, Virginia, 2007. [7] S.W. Doebling, C.R. Farrar, M.B. Prime, D.W. Shevitz, Damage Identification and Health Monitoring of Structural and Mechanical Systems from Changes in Their Vibration Characteristics: A Literature Review, Report no. LA-13070-MS, Los Alamos National Laboratory, 1996. [8] H. Sohn, C.R. Farrar, F.M. Hemez, D.D. Shunk, D.W. Stinemates, B.R. Nadler, A Review of Structural Health Monitoring Literature: 1996-2001, Report no. LA-13976-MS, Los Alamos National Laboratory, 2003. ¨ [9] D. Balageas, C.P. Fritzen, A. Guemes, Structural Health Monitoring, International Society for Technology in Education (ISTE), London, UK, 2006. [10] D.E. Adams, Health Monitoring of Structural Materials and Components: Methods with Applications, Wiley-Interscience, New York, 2007. [11] N.E. Huang, N.O. Attoh-Okine (Eds.), The Hilbert–Huang Transform in Engineering, CRC Press, Boca Raton, FL, 2005. [12] N.E. Huang, S.S.P. Shen, Hilbert–Huang Transform and its Applications, World Scientific Publishing, Hackensack, New Jersey, 2005. [13] P.F. Pai, Three-point frequency tracking method, Struct. Health Monit. 8 (6) (2009) 425–442. [14] P.F. Pai, Time–frequency characterization of nonlinear normal modes and challenges in nonlinearity identification of dynamical systems, Mech. Syst. Signal Process. 25 (2011) 2358–2374. [15] J.N. Juang, R.S. Pappa, An eigensystem realization algorithm for modal parameter identification and model reduction, J. Guidance Control Dyn. 8 (5) (1985) 620–627. [16] P.F. Pai, Long transient phenomenon in nonlinear structural vibration, J. Sound Vib. 321 (2009) 305–318. [17] S.A. Nayfeh, A.H. Nayfeh, Energy transfer from high- to low-frequency modes in a flexible structure via modulation, J. Vib. Acoust. 116 (1994) 203–207. [18] E.O. Brigham, The Fast Fourier Transform, Prentice-Hall, Englewood Cliffs, New Jersey, 1974. [19] Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia, 1992, CBMS-NSF Lecture Notes 61. [20] G. Strang, T. Nguyen, Wavelets and Filter Banks, Wellesley-Cambridge Press, Wellesley, MA, 1997. [21] S.A. Neild, P.D. McFadden, M.S. Williams, A review of time–frequency methods for structural vibration analysis, Eng. Struct. 25 (2003) 713–728. [22] D.F. Shi, L.S. Qu, N.N. Gindy, General interpolated fast fourier transform: a new tool for diagnosing large rotating machinery, J. Vib. Acoust. 127 (2005) 351–361. [23] S.G. Mallat, Z. Zhang, Matching pursuits with time–frequency dictionaries, IEEE Trans. Signal Process. 41 (12) (1993) 3397–3415. [24] P.F. Pai, Circular Instantaneous Frequency, Adv. Adaptive Data Anal. 2 (1) (2010) 39–64. [25] S.L. Hahn, Hilbert Transforms in Signal Processing, Artech House, Boston, 1996. [26] M.R. Hajj, R.W. Miksad, E.J. Powers, Perspective: measurements and analyses of nonlinear wave interactions with higher-order spectral moments, J. Fluids Eng. 119 (1) (1997) 3–13. [27] J.M. Nichols, A. Milanese, P. Marzocca, The analytical trispectrum for multiple degree-of-freedom systems possessing cubic nonlinearity, in: Proceedings of the SPIE 15th International Symposium on Smart Structures and Materials & Nondestructive Evaluation and Health Monitoring, San Diego, California, March 9–13, 2008. [28] S. Rossignol, P. Desain, , State-of-the-art in fundamental frequency tracking, in: Proceedings of the Workshop on Current Research Directions in Computer Music, Barcelona, Spain, 2001, pp. 244–254. [29] P. Flandrin, Time–Frequency/Time–Scale Analysis, Academic Press, 1999. [30] G. Kerschen, K. Worden, A.F. Vakakis, J.C. Golinval, Past, present and future of nonlinear system identification in structural dynamics, Mechanical Systems and Signal Processing 20 (2006) 505–592. [31] M. Feldman, Nonlinear system vibration analysis using the Hilbert transform—I. Free vibration analysis method ‘FREEVIB’, Mech. Syst. Signal Process. 8 (1994) 119–127. [32] M. Feldman, Nonlinear system vibration analysis using the Hilbert transform—II. Forced vibration analysis method ‘FORCEVIB’, Mech. Syst. Signal Process. 8 (1994) 309–318. [33] P.F. Pai, A.N. Palazotto, HHT-based nonlinear signal processing method for parametric and non-parametric identification of dynamical systems, Int. J. Mech. Sci. 50 (12) (2008) 1619–1635. ¨ [34] M. Datig, T. Schlurmann, Performance and limitations of the Hilbert–Huang transformation (HHT) with an application to irregular water waves, Ocean Eng. 31 (2004) 1783–1834. [35] F. Wu, L. Qu, An improved method for restraining the end effect in empirical mode decomposition and its applications to the fault diagnosis of large rotating machinery, J. Sound Vib. 314 (2008) 586–602. [36] G. Rilling, P. Flandrin, One or two frequencies? The empirical mode decomposition answers, IEEE Trans. Signal Process. 56 (1) (2008) 85–95. [37] R.M. Mersereau, M.J.T. Smith, Digital Filtering, Wiley-Interscience, New York, 1981. [38] P. Maragos, J.F. Kaiser, T.F. Quatieri, Energy separation in signal modulations with application to speech analysis, IEEE Trans. Signal Process. 41 (10) (1993) 3024–3051.