Characterization of Non-linear Bearings Using the Hilbert-Huang Transform

Changes in the performance of bearings can significantly vary the distribution of internal forces and moments in a structure as a result of environmental or operational loads. The response of a bearing has been traditionally idealised using a linear model but a nonlinear representation produces a more accurate picture at the expense of modelling complexity and computational time. In this paper, a lead-rubber bearing is idealised using the hysteretic Bouc-Wen model. The Hilbert-Huang Transform is then employed to characterize features of the non-linear system from the instantaneous frequencies of the bearing response to a time-varying force. Instantaneous frequencies are also shown to be a useful tool in detecting sudden damage to the bearings simulated by a reduction in the effective stiffness of the force-deformation loop.


Introduction
Generally, it is reasonable to assume that the response of a structure to moderate dynamic loads will remain in the elastic region, however, large deformations or specific components with non-linear physical properties are represented best by non-linear models [1]. An example of a nonlinear structural component is an elastomeric bearing. Elastomeric bearings are the most commonly used type of seismic isolator [2] and their main purpose is to allow lateral movement due to shrinkage, temperature variation and earthquakes. In addition to this lateral flexibility, bearings can have a degree of compressibility in the vertical direction, typically modelled with translational springs. In some cases, rotational springs are also needed to accurately resemble field measurements from static and dynamic loading tests. For instance, there is experimental evidence that changes in joints and elastomeric bearings are one of the main causes of changes in bridge frequencies [3][4][5].
Of relevance to this investigation, is the non-linear nature of elastomeric bearings due their inherent damping properties [6]. The nonlinearity can be described by the restoring force and displacement behaviour, i.e., a hysteresis loop describing the recoverable and permanent deformations. A number of mathematical hysteresis models that predict the bearing response for different levels of stress are presented in the literature [7][8][9][10][11]. Hwang et al [7] experimentally validates a mathematical model that replicates the hysteretic behaviour of hard rubber bearings for a variety of external loading conditions. Matsagar and Jangid [12] investigate the influence of different isolator characteristics, such as the shape of the isolator force-deformation loop and yield displacement of the isolator, on the seismic response of the structure. They find that the seismic response is significantly influenced by the shape of the loop and that lower yield displacements increase the acceleration response of the structure. In this paper, elastomeric bearings are assumed to be lead rubber bearings modelled as vertical springs. Lead rubber bearings present several convenient features with fundamental drawbacks related to the negligible increase in damping and high deformability for low static loads. The lead plug increases the horizontal stiffness significantly compared to other bearings (e.g. laminated rubber bearings). With the placing of lead plugs ( Fig. 1(a) where ( ) is the external vertical load applied to the bearing), aside from energy dissipation for seismic response and stiffness for static loads, a single compact device is obtained, able to satisfy most of the requirements associated with a good bearing. For that reason rubber bearings are used extensively in practical applications in bridge structures [13]. Fig. 1(b) represents a bilinear hysteretic model for simulating the response a lead rubber bearing. In this figure, ( ) is the time-varying vertical displacement of the bearing and ( ) is the restoring force in the bearing. The stiffness of the bearing (slope ) is related to post-yield stiffness by . The lead responds with elastic perfectly plastic loops, therefore, after yielding (f s (t) > f sy in Fig. 1(b)), the stiffness is equal to the stiffness of the rubber bearing alone (K v ) and the loops are almost bilinear [14].
is the effective stiffness value given by the ratio of the maximum restoring force ( ) to the maximum displacement ( ) [15]. Deterioration in elastomeric bearings is caused by external restraint forces which are generated in the movement of concrete in response to temperature, creep and shrinkage. The restraint may consist of friction at the bearings, bonding to already hardened concrete or by attachment to other components of the structure. Cracks resulting from the actions of external restraint forces develop in a similar manner to those caused by external load. Concern about the reliability of elastomeric bearings is increasing and studies on experimental and in-service bearings reveal the occurrence of permanent changes in the engineering properties of the bearings [4,5]. Cyclic loading is found to cause progressive damage, analogous to fatigue in metals [16]. Most damage cases are in the form of progressive debonding between the rubber and steel shims. A degree of stiffness loss can be expected due to repeated cyclic loading causing breading of molecular bonds between polymers [17]. A method to characterize and monitor the structural parameters characterizing the bearing response is clearly needed.
Many methods have been proposed for parameter identification of Single Degree-Of-Freedom (SDOF) and Multiple Degree-Of-Freedom (MDOF) structural systems. Recently, wavelet ridges of Continuous Wavelet Transform (CWT) have been employed to identify the Instantaneous Frequency of a time-varying structure, where the frequency change is gradual representing a linear system with time-varying frequency [18]. Nagarajaiah and Basu [19] investigate modal identification of linear time-variant systems and they conclude that the wavelet analysis is more appropriate than the Hilbert-Huang Transform (HHT) in identifying changes in the frequency of linear systems due to the number of scales and coefficients that can be applied. Feldman [20] applies the Hilbert Transform (HT) to identify the nonlinear instantaneous modal parameters of SDOF systems in forced and free vibration. However, the application of HT is limited as it only defines a 'mono-component' signal (i.e., obtain only one frequency value using a single component) [21]. For 'multicomponent' signals, the concept of instantaneous frequency becomes meaningless, given that the HT is inherently limited to 'mono-component' signals and breakdown of the signal components is needed. The signal needs to be defined in terms of its single 'mono-components' each having a different instantaneous frequency [8] as addressed by Huang et al [22]. They decompose the signal using Empirical Mode Decomposition (EMD) to obtain the Intrinsic Mode Functions (IMFs) and then use the HT to obtain the instantaneous frequency of the signal. This methodology is called the HHT and it is one of the most popular methods for analysis of nonlinear and non-stationary data [22][23][24]. In the following sections, the HHT is applied to the characterization of the response of nonlinear systems with hysteresis that define lead rubber bearings. First, basics on the HHT are reviewed. Second, the mathematical model of the bearing (a SDOF) is introduced. Then, the response of the model to a time-varying force and to a sudden change in the forcedisplacement hysteresis loop is investigated via the HHT.

The Hilbert Huang Transform
There are two main steps in the HHT analysis process [22,25]: EMD and the HT. The HT characterizes a time-domain signal ( ) by extracting its envelope ( ) and its phase ( ) which can be expressed as ( ) = ( )cos θ( ). Huang et al [22] applies the EMD process first to overcome the limitation of the HT to only 'mono-component' signals.

First step: EMD
The signal is decomposed into a number of IMFs using the 'sifting' process. The IMFs can be defined as a representation of oscillatory modes of variable frequency and amplitude that are time dependent, as opposed to simple harmonic functions of constant amplitude and frequency. The IMFs must satisfy two conditions: the number of extrema and the number of zero crossings must differ by no more than one, (b) the mean value of the envelope defined by the maxima and the minima at any given time must be zero.
The sifting process to obtain the IMFs is as follows: i) identify the local maxima and minima, ii) find the mean m 1 (t) of the upper and lower envelopes, iii) obtain the difference ℎ 1 ( ) between the signal x(t) and the mean m 1 (t), i.e., ℎ 1 ( ) = ( ) − iv) Check if ℎ 1 satisfies the two conditions ((a) and (b) above) to be an IMF. If the two conditions are met, then h 1 will be the first IMF. However, this rarely occurs in the first iteration, and therefore ℎ 1 must be treated as the signal and the process i) to iii) repeated again (i.e., ℎ 11 ( ) = ℎ 1 ( ) − 11 ( )), until the sifted signal satisfies the two conditions for an IMF.
After the sifting process i) to iii) has been done a number of times i.e., k, the two conditions will be met and the first IMF,   1 ct, will be obtained (Equation (1)).
However, too many sifting cycles could reduce all the components to a constant-amplitude signal with frequency modulation only, i.e., only frequency variation is retained. To guarantee that the IMF components retain enough physical sense of amplitude and frequency modulations, the number of times, k, the sifting process repeats has to be limited. There are two main mathematical stoppage criteria proposed in the literature that are used to ensure that the sifting process stops when conditions (a) and (b) above are met [22,23]:  A first criterion based on limiting the normalised Square Difference (SD) between two successive sifting operations to a small value, i.e.: at time and for a total duration of the signal.
 A second criterion based on assuming the sifting process has ended after a predefined number of times. For instance, if after S sifting times there is no difference between the number of zero crossings and the extrema, or differ at most by one, then the process is stopped.
The first criterion is difficult to implement in practice due to the need to select an accurate value of to stop the sifting process. A small value of SD does not necessarily satisfy the number of zero crossings and extrema. For that reason, the second criteria is adopted in this paper, although there also exists difficulty in predefining the value of S which it is generally taken between 4 and 8 [23]. The value of S adopted here is 4 to save computational time given that choosing a higher value has not shown to significantly enhance the results of the analysed signals.
The process of sifting starts by separating the finest local mode   1 ct(Equation (1)) from the data x(t). This gives the residual The residual is then used as the input signal to obtain the second IMF 2 ( ( )) c t and so on. The decomposition into higher order IMFs is repeated m times until the value of () m rt becomes a 'monotonic' component from which no more IMFs are extracted.

Second step: The HT
Once the EMD process has been finalised, the second step of the HHT analysis is the HT of the IMFs [24,25]. A complex analytical signal is formed for each IMF consisting of a real part defined as the th i IMF ( ( )), and an imaginary part defined by the HT ( () i y t in Equation (3)).
The complex signal can be expressed in exponential format as in Equation (4).
The instantaneous frequency   i IF t associated to the i th IMF at time t can be obtained from: Using the HHT, it is possible to characterize the frequency of a nonlinear system such as a nonlinear support. The potential of the HHT in detecting changes in bearing stiffness is examined in the following sections.

Description of Mathematical Model of Bearing
An elastomeric bearing is idealised as a SDOF system consisting of a mass (M) supported on a spring ( () v k t ) and damper (C v ) system ( Fig. 2(a)). The response of a SDOF system with a nonlinear hysteretic stressstrain relationship has been analysed by Dobson et al [9], who provide a mathematical model based on Bouc-Wen capable of simulating the asymmetric hysteretic response found in bearings. The asymmetric hysteresis is composed of four different curves that describe the force-displacement response of the bearings. Here, the nonlinear stiffness in the spring, K v (t), is given by the slope of the curve in Fig. 2(b), which is a smoother and more realistic representation of a bearing response than the simplified bilinear mathematical model consisting of straight lines shown in Fig. 1(b) [26]. In Fig. 2(b), z is the hysteretic displacement which defines the permanent displacement in the hysteresis, i.e. at unloading the displacement does not return to the original value. given in Equation (7).
where   ut,   ut and () ut are the acceleration, velocity, and displacement respectively of the SDOF, which can be obtained by integration of Equation (7) using the Wilson-θ method. The response of the nonlinear SDOF system in Equation (7) can be re-formulated [27] as: where  is damping, () t  is the natural frequency given by ft is the restoring force for nonlinear system given by: In Equation (9),  is the ratio of post-yield to pre-yield stiffness and () zt is the hysteretic displacement (i.e., the difference between the displacement at loading and unloading (permanent displacement)). The post-yield to pre-yield stiffness ratio (  ) regulates the shape of the hysteresis. If  =1 the force-displacement relationship reverts to linear case, while for  =0.2 (0≤ ≤ 1) a hysteresis force-displacement relationship, as illustrated in Fig. 2 (b), is formed.
The Bouc-Wen hysteretic displacement ( () zt , which is given by the horizontal axis in Fig. 2(b)) is used to guarantee a smooth transition between pre-yield and post-yield stiffness that represents the real forcedeformation behaviour of the bearing. () zt is governed by Equations (10) and (11) where the influences of the characteristic strength, damping, post to pre-yield stiffness ratio and other parameters are taken into account via pre-established coefficients [27]: where   ztis the incremental hysteretic displacement,   utis the velocity of the SDOF system, are coefficients that define the shape of the hysteresis ( A ,  , influence the loop size and n influence its smoothness), t  is the incremental time step and () zt t   is the hysteresis displacement at the previous time step. In general, by using Equation (10) the coefficients are dependent on each other to produce a hysteretic loop that is representative of the hysteresis bearing response. In order to have more flexibility and better control over the coefficients that define the hysteresis shape, Equation (12) is used instead of Equation (11) to obtain the force-displacement loops presented in Fig. 2(b) as described by Dobson et al [9]: Equation (12) decouples the hysteretic loop in four loading regimes each of which being regulated by a different constant shape coefficient: Fig. 2(b)). The constants A and n define the shape and level of nonlinearity respectively in the restoring force-deformation behaviour, and the constants B, C, D and E shape the hysteresis loop. The constant shape parameter to employ will depend on the location of the force-displacement curve and the signs of the velocity ( () ut ) and the hysteretic displacement ( () zt ). Three of the nonlinear terms in Equation (12) Fig. 2(b).
When ( ) 0 ut  and ( ) 0 zt  the system is linear and there is no constants that need to be defined for the force-displacement curve.
The main advantages of the model in Equation (12) over the Bouc-Wen model in Equation (11) is the added flexibility to the hysteretic cycle as each curve is defined separately. Setting the nonlinear parameters B, C, D and E equal to zero results in a linear relationship with a slope equal to that of the initial tangent slope. Depending on the selected coefficients, choosing a positive value shall bow the curve inwards (or outwards) from the loop, with the negative value producing the opposite effect.

Characterization of Non-Linear SDOF Bearing Model Using the HHT
The response of the SDOF in Fig. 2 (   ut) is obtained for a sinusoidal load   ft given in N by Equation (13) for testing the ability of the HHT-based method in characterizing the structural system.
In this theoretical test, the SDOF bearing model has a mass of M =3135.8 kg and a viscous damping of 6% (ξ = 0.06) [5,28,29]. A vertical effective stiffness of 1.14375×10 9 N/m is assumed which leads to a 'effective' natural frequency of 30.7 Hz, The following values are assigned to the constants of the hysteresis model of Equation (12) to define the shape of Fig. 2 , and 500 E  . These values are selected to create a force-displacement curve that resembles previous experimental investigations on the response of hysteretic bearings [30][31][32].

Figs. 3(a) and (b) show the acceleration (  
ut) and displacement ( () ut ) respectively of the SDOF system under the sinusoidal load specified in Equation (13). The four stages (B, C, D, and E) of the hysteresis in Fig. 2(b) are shown by rounded dashed circles in Fig. 3 for the response between 1.1 s to 1.73 s (which covers one full hysteresis loop 0.63 s long, i.e., the time from start of loading at B to unloading at D and E to reloading at C and back to B in the curve). When the curve is at C, the stiffness is relatively high and the acceleration response (from 1.1 s to 1.225 s) has lower amplitude, which is similar to the stage when the curve is at D (from 1.425 s to 1.54 s). When the curve is at B (from 1.225 s to 1.425 s) the lower stiffness produces higher amplitude and a relatively linear force-displacement relationship ( Fig. 2(b)), similarly to stage E from 1.54 s to 1.73 s.
(a) (b) Figure 3. Response of non-linear SDOF system to a sinusoidal load: a) acceleration, b) displacement The HHT is applied to the acceleration signal from Fig. 3(a). Figs. 4(a) and b) show the first four IMFs resulting from the EMD process and the associated instantaneous frequencies via the HT respectively.
Although the total number of IMFs extracted from the acceleration is seven, only the first four IMFs are shown in the figure. The last three IMFs (IMF5-7) are residues of small magnitude that do not provide meaningful information on the system behaviour. Lower IMFs contain higher frequency components of the signal, and IMF1 appears to have higher coefficient values than other IMFs and to dominate the response.

Figure 4. HHT: (a) IMFs 1-4 using EMD. (b) associated instantaneous frequencies (IF 1 (t) -IF 4 (t)) in Hz using the HT
The frequency of the applied sinusoidal load (Equation (13)) is 1.59 Hz and it is captured by IMF4 in Fig. 4. The natural frequency of the SDOF system, NF(t), is time-varying and determined in Hz by: where the stiffness () v kt is obtained from the force-displacement curve in Fig. 2(b) at each point in time t . Fig. 5 shows the close correlation between the IF 1 (t) extracted from the IMF1 and the natural frequency of the system, NF(t). Discrepancies between IF 1 (t) and NF(t) are noticed for the lowest stiffness in the loop (stages B and E), i.e., IF 1 (t) predicts about 12 Hz for a minimum NF of 14 Hz. The latter can be partially attributed to mode mixing in the EMD process. Although the IF1 and using HHT does not exactly match that of the 'actual values' of the system, it closely approximates the exact frequency and stiffness of the system. Figs. 7 and 8 show the analysis of the acceleration in Fig. 3 using Fast Fourier and Continuous Wavelets Transforms (CWT) respectively to compare with the HHT. In Fig. 7, the Power Spectral Density (PSD) shows the range of possible frequencies of the system with the highest peaks taking place at 1.6 Hz and 4.8 Hz. These two large peaks are followed by others of smaller magnitude up to 42 Hz representing the continuous change in frequency of the system with time. It is evident that the PSD is very limited in characterizing a non-linear system due to the lack of time information. There is no way of knowing if these frequencies correspond to stages B, C, D, or E on the force-displacement hysteresis curve.  Fig. 3(a) The CWT uses inner products to measure the similarity between a signal and a function known as the mother wavelet. In the CWT, the signal is compared to shifted and compressed or stretched (scaled) versions of the mother wavelet. By comparing the signal to the wavelet at various scales and points in time, a frequency-time map of wavelet coefficients is obtained [33]. Fig. 8(a) shows the wavelet coefficients versus pseudo-frequency and time that result from applying the CWT to the acceleration response of Fig. 3(a) using the Mexican Hat mother wavelet. The pseudo-frequency is directly proportional to the centre frequency of PSD (m 2 /s the wavelet in Hz (For the Mexican Hat wavelet, this centre frequency is 0.25 Hz), and inversely proportional to the scale of the wavelet and to the sampling period. It is difficult to interpret the behaviour of the system from Fig. 8 alone, so sections at times 1.1, 1.225, 1.425, 1.54, and 1.73 s (A-A, B-B, C-C, D-D, and E-E respectively) are plotted versus frequency in Fig. 8(b) and sections at frequencies 14 and 41 Hz (F-F and G-G respectively) are plotted versus time in Fig. 8(c). In Fig. 8(b), section A-A corresponds to initial stage of C, while section B-B corresponds to end of stage C and beginning of stage B. Sections through C-C and D-D correspond to end of stage B and beginning of stage D. It can be seen that the sections C-C and D-D have the same magnitude than sections A-A and B-B but they are in the negative region. Section E-E corresponds to end of stage E and it is almost equal in magnitude to section A-A suggesting that the stiffness is returned to its original value and a full hysteresis cycle is complete. Fig. 8(c) shows the variation in content of the pseudo-frequencies 14 and 41 Hz (lower and upper values respectively of the range of variation of the frequency of the system) with time. However, the range of variation of the frequency of the system, the frequency at which the system is responding for each point in time and its duration is not as clearly defined by wavelets as in Fig. 5 simply using the IF 1 (t) by HHT. It has been shown that HHT can be used to characterize frequency and stiffness changes due to nonlinearity properties of a SDOF system, although one problem that may arise with the EMD process is over decomposition of the signal resulting into excessive IMFs and mode mixing, i.e., IMFs sharing similar frequencies. It is also important to note that this is a theoretical investigation of a SDOF system with a specific hysteretic loop and that experimental results may contain other complexities involving external factors such temperature and/or inconsistencies of the applied load and the specimen, and/or internal factors such as friction between steel shims and rubber (for steel reinforced rubber bearings).

Identification of Changes in the Bearing Behaviour Using the HHT
Although the HHT was initially develop to characterize non-linear systems, its application to structural damage detection has been limited mostly to sudden losses in stiffness using simple periodic functions, experimental beams with crack, or linear time variant systems [19,[34][35][36][37]. The reason for this possibly lays on the lack of data to characterize how damage affects the overall performance of non-linear structural systems. In this investigation, damage is theoretically simulated through the introduction of a reduction in the effective stiffness of a hysteretic loop as shown in Fig. 9(a) [38]. Acceleration response of the SDOF system is illustrated in Fig. 9(b) where a 10% reduction in effective stiffness is introduced after 2.5 s. This change in the effective stiffness of the hysteretic loop is visualized as an increase in the maximum amplitude of the acceleration signal of 0.282 m/s 2 with respect to the original acceleration (i.e., maximum acceleration is 28.8% higher after the stiffness reduction). This large increase in acceleration is attributed to the low stiffness values experienced by the bearing at loading and unloading. In comparison, at higher stiffness the change in acceleration is much lower (6.2%).  Figure 9. (a) Force-displacement relationship for SDOF system before and after 10% loss in effective stiffness, (b) acceleration response before and after 10% loss in effective stiffness Four IMFs of interest can be obtained from Fig. 9(b) using EMD and the IFs of these IMFs are shown in Fig. 10(a). Edge effects (before 0.5 s and after 4.5 s) appear in the IFs as a result of the decomposition process [22]. IF 1 (t) to IF 4 (t) show a sudden peak at the instant of change in effective stiffness. This peak can be clearly visualized in IF 3 (t) and IF 4 (t) of Fig. 10(a), due to the relative high frequency content at this location with respect to the rest of the IMF. The dominant IMF is IMF1 and Fig. 10(b) zooms in on the IF 1 (t) shown in Fig. 10(a). The average range of IF 1 (t) is 42. 35-10.64 Hz prior to the simulated damage and 38. 35-9.46 Hz after damage has occurred. Once edge effects are ignored, the reduction in frequency as a result of the stiffness change is 9.6% (-4 Hz) at the upper range (cases C and D in Fig. 2(b)) and 11.1 % (-1.18 Hz) at the lower range (cases B and E in Fig. 2(b)). The average drop in stiffness is (9.6%+11.1%)/2=10.35%, i.e., approximately the 10% reduction in effective stiffness.
(a) (b) Figure 10. Instantaneous frequencies in Hz before and after 10% effective stiffness loss of (a) IMFs 1-4 This section has demonstrated the capability of the HTT to extract important time-frequency information about the nonlinear behaviour of a SDOF system. IF 1 (t) has been able to capture a small change in effective stiffness of 10%. It is acknowledged that experimental setups and more complex systems will require IFs from more than one IMF to gather sufficient information to characterize the system.

Influence of Noise in the Performance of the HHT
Noise can be caused by a number of factors which include the natural intermittent instabilities in the system, the concurrent phenomena in the environment where investigations are conducted, and/or the sensors and recording systems. As a result, data is always an amalgamation of signal and noise, that here it is simulated using the additive model: where {̈( )} is the acceleration corrupted by noise, {̈( )} is the acceleration of the system at time , is the noise level, { } is the standard normal distribution vector with zero mean and unit standard deviation, and is the standard deviation of the noise-free signal. For testing purposes, the original acceleration data of Fig. 3(a) is corrupted with 1% noise level (i.e., = 0.01 in Equation (15)) which results into Fig. 11. Figure 11. Acceleration response of non-linear SDOF system to a sinusoidal load with 1% noise Following EMD of the corrupted acceleration signal, the corresponding IMFs 1 -8 are represented in Figs. 12(a) and (b). Compared to the noise-free signal (Fig. 4), there is a higher number of IMFs that need to be considered as a result of the presence of noise. It can be seen that even low noise levels can have a significant impact on the decomposition process. Figs. 12(c) and (d) show the instantaneous frequencies associated to each IMF. IMFs 1-3 reveal high variability of frequency content between 0-500Hz, 0-300Hz and 0-150Hz for IMF1, IMF2, and IMF3 respectively, which can be attributed to noise. IMFs 4-8 contain a smoother variation of frequencies that appears to be in the range of instantaneous frequencies obtained in Fig. 4(b) for noise-free acceleration. However, when plotting the instantaneous frequencies of IMF4 and IMF5 (Fig. 12), together with the timevarying natural frequency of the system in Fig. 13, the influence of noise is clear. The extraction of a clean instantaneous frequency that matches the frequency of the system is rendered unsatisfactory.
(a) (b) Figure 13. Comparison of instantaneous frequency of IMF using HT on a 1% noisy signal and actual instantaneous natural frequency (NF(t)) of the system: (a) Instantaneous frequency of IMF4 and NF(t).
For cases where noise has distinct time or frequency scales from those of the true signal, Fourier transform filters can be applied to separate the noise from the signal. Once the signal is filtered from noise, the HT has been shown to be capable of extracting phases and instantaneous frequencies for damage identification of linear scenarios [39,40]. Similarly, Sections 4 and 5 have demonstrated that the HHT is a valid tool to characterize instantaneous changes using the response of a non-linear SDOF bearing system where noise is non-existing or has been safely removed. However, traditional filtering methods may not be as reliable in separating mixed noise and harmonics in non-linear systems, and Section 6 has shown that the direct application of EMD and the HT to a corrupted signal can produce misleading results. In this regard, Ensemble Empirical Mode Decomposition (EEMD) [41] is an alternative to EMD that attenuates mode mixing and noise to some extent.

Conclusions
A hysteretic Bouc-Wen model has been implemented in a SDOF system to simulate the response of a lead rubber bearing. An instantaneous frequency analysis of the response of the SDOF via the HHT has provided information that allows characterizing the different stages of stiffness in this non-linear system (i.e., stages B , C , D , and E of the hysteretic loop). The HHT has overcome the inability of the PSD to distinguish the multiple frequencies associated to one mode of vibration in a non-linear system. It has also outperformed wavelet analysis in characterizing the system. A sudden change in effective stiffness has been successfully identified via the instantaneous frequency. Finally, it has been shown that the accuracy of the method can be severely limited by noise. Further research is needed to address the impact of noise and of integrating the bearing as part of a more sophisticated MDOF system on the performance of the HHT.