Identification of modal parameters of a time invariant linear system by continuous wavelet transformation

Mechanical Systems

and

Signal Processing

Mechanical Systems and Signal Processing 21 (2007) 1642–1664

Identification of modal parameters of a time invariant linear

system by continuous wavelet transformation

C.S. Huang



, W.C. Su

Department of Civil Engineering, National Chiao Tung University, 1001 Ta-Hsueh Road, Hsinchu 30050, Taiwan Received 9 January 2006; received in revised form 17 July 2006; accepted 17 July 2006

Available online 22 September 2006

Abstract

The applications of wavelet transforms have received significant attentions in many fields. This work proposes a procedure for identifying modal parameters of a linear system using the continuous wavelet transform. The merits of the proposed procedure over the exiting schemes of applying the wavelet transform to system identification for a linear system are in use of the time invariance property and filtering ability of the transform to enhance the efficiency of identifying the modal parameters of a structure from its earthquake responses or free vibration responses. The effectiveness and accuracy of the proposed procedure are validated via numerical simulations. The effects of noise and wavelet function on identifying the modal parameters of the structure are also explored in processing the numerically simulated acceleration responses of a six-story shear building subjected to base excitation. The dynamic characteristics of close modes are accurately determined. Finally, the proposed procedure is adopted to obtain the modal parameters of a three-story non-symmetric steel frame from its measured acceleration responses in a shaking table test. A total of nine modes are identified, including modes with high frequencies and very small amplitude.

r2006 Elsevier Ltd. All rights reserved.

Keywords: System identification; Linear system; Wavelet transformation

1. Introduction

Health monitoring of a system is an important and interesting issue in both academia and industry. Structural health monitoring aims at determining the durability and safety of structures. A structural system may sustain damage when subjected to severe loading such as a strong earthquake, or when its material deteriorates. The dynamic characteristics of a system are often identified from its dynamic responses in a health monitoring program of the system for the following reasons. Damage to the system reduces its natural frequencies, raises its modal damping and alters its mode shapes. The identified dynamic features can be employed to ensure the accuracy of the numerical model of the system in the design stage, and to correct the model if necessary. The corrected numerical model can then be applied to future damage assessment of the system.

www.elsevier.com/locate/jnlabr/ymssp

0888-3270/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2006.07.011

Corresponding author. Tel.: +886 3 5712121x54962; fax: +886 3 5716257. E-mail address:[email protected] (C.S. Huang).

The dynamic characteristics of a system are often identified either in the frequency domain or the time domain. Identification techniques in the frequency domain (i.e. spectrum analysis[1]) are widely recognised as simple, but do not accurately estimate the spectra of measured responses, leading to the loss of accuracy of identified modal parameters, particularly for highly damped systems and systems with severe modal interference. Conversely, methods in the time domain (e.g. time series methods [2], the subspace approach

[3,4], and the ANN approach[5]) frequently provide accurate results if the measured responses are not severely

contaminated by noise. However, filtering noise from the measured responses in the time domain is not so convenient as that in the frequency domain.

A time-frequency analysis method, called wavelet transform, has been fully developed in the theoretical aspect over the past two decades, having been successfully applied in many fields. Various techniques for identifying modal parameters of a linear system have been developed based on discrete wavelet transform or continuous wavelet transform. Two approaches using the discrete wavelet transform have been proposed. Robertson et al.[6,7]proposed a procedure for extracting impulse response characteristics of a system from its forced responses by expanding the impulse response function and input force in the Duhamel integral equation in terms of wavelet basis functions. Robertson et al. utilised Daubechies wavelets[8]in their study. Then, an eigensystem realisation algorithm was adopted to identify the dynamic characteristics of the system from the estimated impulse response functions. Huang et al.[9]applied the discrete wavelet transform to the discrete equations of motion and determined the natural frequencies, modal damping ratios and mode shapes using either free decay responses or earthquake responses. They examined three different mother wavelet functions called ‘‘sym1’’, ‘‘sym4’’, and ‘‘sym10’’[8]in their numerical studies.

Some works have been published on identifying the dynamic characteristics of a linear system according to the continuous wavelet transform. Schoenwald[10]applied the continuous wavelet transform to the equation of motion of a single-degree-freedom system and identified the parameters in the equation of motion. Ruzzene et al. [11], Gouttebroze and Lardies[12], Lardies and Gouttebroze[13], and Le and Argoul[14]applied the continuous wavelet transform to the free decay responses of a system using various wavelet functions and proposed different procedures for determining the natural frequencies, modal damping ratios, or mode shapes of the system from the modulus or the phase of the wavelet transform. These procedures are valid for systems with weak damping. Furthermore, the success of these techniques depends on the ability of the frequency localisation of the wavelet transform to decompose the responses of multiple modes accurately into the responses of the single mode of interest. Mode decomposition is typically difficult for close modes with severe modal interference.

Previous works have some disadvantages. The schemes using the discrete wavelet transform cannot employ the time invariance property of continuous wavelet transform, which gives a simple relationship between the wavelet transform of f(t) and f(ttn), where tnis a constant. Consequently, these schemes consume much

computation time in evaluating the discrete transforms of functions with different time shifts involved in their formulations. Moreover, since the scale parameters in the discrete wavelet transform are often defined in powers of two, an efficient frequency filter cannot be achieved from the discrete wavelet transform. These drawbacks do not occur in the continuous wavelet transform. However, existing algorithms based on the continuous wavelet transform can only process the free decay responses of systems with weak damping. This study develops a novel scheme free from these shortcomings.

The proposed approach applies the continuous wavelet transform to the measured dynamic responses of a linear structure. The discrete linear equations of motion among the measured degrees of freedom are then established in the wavelet domain. A fixed value or several fixed values of the scale parameter in the transformation can be adopted to filter out the uninteresting responses. The coefficients in the discrete equations of motion, which determine the dynamic characteristics of the system, are found in the wavelet domain through the least squares approach. The modal parameters are then computed directly from these coefficients by solving an eigenvalue problem. The numerically simulated earthquake responses of a six-story shear building are first processed to verify the effectiveness and accuracy of the proposed approach. Numerical experiments are then undertaken to determine the effects of noise, wavelet function and scale parameter selection, and other factors, on the accuracy of determining dynamic characteristics. Several wavelet functions such as the Shannon, Meyer, Morlet and Haar functions are studied. The feasibility of the proposed approach to elucidate real structures is demonstrated by processing the measured responses of a three-story steel frame

in shaking table tests. The frame is a non-symmetric structure, and its horizontal and torsional responses of floor are coupled.

2. The continuous wavelet transform

The wavelet transform of a function of time, f(t), belonging to L2space, which meansR_11 f2ðtÞ dt finite, is defined as Wcf ða; bÞ ¼ f ðtÞ; ca;bðtÞ   ¼ 1_ffiffiffi a p Z 1 1 f ðtÞc t  b a   dt, (1)

where /,S and the superscript * denote inner product and the complex conjugate, respectively; a is a dilation or scale parameter, which is typically a positive real and plays the role of the inverse of frequency; b is a translation parameter, which indicates the locality of the transformation; c_a;bðtÞ ¼ ðpffiffiffiaÞ1cððt  bÞ=aÞ, and c(t) is a mother wavelet function. A mother wavelet function must satisfy the admissibility condition,

Cc¼ Z 1 1 ^ cðoÞ    2 o j j doo1, (2)

where ‘‘||’’ represents the magnitude, and ^cðoÞ is the Fourier transform of c(t)[15]. The inverse of the wavelet transform is given by f ðtÞ ¼ 1 Cc Z 1 1 Z 1 1 Wcf ða; bÞca;bðtÞ 1 a2da db. (3)

The wavelet transform is a linear transform with the following time invariance property[16]: Let gðtÞ ¼ f ðt  tÞ,

then

Wcgða; bÞ ¼ Wcf ða; b  tÞ. (4)

This property , which is very important in this work, can be easily proved by using the definition given in Eq. (1). This relationship given by Eq. (4) cannot usually be employed in the discrete wavelet transform because, in the discrete wavelet transform, a and b are often chosen as a ¼ am₀ and b ¼ nb0am0, where m and n belong to Z.

The wavelet transform decomposes an arbitrary function f(t) into the functions ca;bðtÞ with different a and b.

Hence, the wavelet transform alters the frequency contents of f(t) if the scale parameter is fixed, meaning that the transformation has a filtering effect. Applying Fourier transformation to Eq. (1) with a fixed value of the scale parameter a yields

j ^Wcf ða; oÞj ¼ ffiffiffia

p

j ^cðaoÞj j ^f ðoÞj, (5)

where ^Wcf ða; oÞ and ^f ðoÞ are the Fourier transforms of Wcf ða; bÞ with respect to b and f(t), respectively.

The wavelet transform gives the localised information of the spectrum of f(t) with a frequency window[15], o a  D_c^ a ; o a þ D_c^ a  , (6) where o_{=a and 2D} ^

c=a are the center and width of the window; respectively, and o

*

and D_c^ are defined as

o _¼ 1 ^ cðoÞ 2 Z 1 1 o ^_cðoÞ_2do and Dc^ ¼ 1 ^ cðoÞ Z 1 1 ðo  o_Þ2_cðoÞ^  2do  1=2 , (7) where ^ cðoÞ 2¼ Z 1 1 ^ cðoÞ    2do.

3. Methodology 3.1. Background

The equations of motion of a time invariant linear system with multiple degrees of freedom are

½Mf €xg þ ½C f _xg þ ½K fxg ¼ ff g, (8)

where [M], [C], and [K] are the mass, damping and stiffness matrices of the structure system, respectively; f €xg, f_xg, and {x} are the acceleration, velocity and displacement response vectors of the system; and {f} is the input force vector. Eq. (8) is not restricted to the proportional damping system and can be accurately discretised by the impulse invariant transformation or the method of mapping of differential[17]and expressed as[18–20]

xðtÞ  ¼X 2 i¼1 ½ ¯Ui xðt  iDtÞ  þ H¯ ₁ f ðt  DtÞ, (9a) _xðtÞ  ¼X 2 i¼1 ½ ¯U_i _xðt  i DtÞþX 1 j¼0 ^ H h i j f ðt  j DtÞ  , (9b) or €xðtÞ  ¼X 2 i¼1 ½ ¯U_i €xðt  i DtÞþX 2 j¼0 ~ H   j f ðt  j DtÞ  , (9c)

where ½ ¯U_i, ½ ¯H₁, ½ ^H_j, and ½ ~H_jare coefficient matrices related to [M], [C], [K], and Dt that represents the time increment.

Real applications in structural health monitoring often only measure acceleration or velocity responses at some (not all) degrees of freedom due to economic concern. The discrete equations of motion corresponding to the measured degrees of freedom are[18,20]

yðtÞ



¼X

I i¼1;2

½U_i yðt  iÞþX

J j¼0;1

½H_j f ðt  jÞ, (10)

where {y(ti)} and {f(ti)}are the measured responses, in terms of acceleration or velocity, and the forces at time ti Dt, respectively; I and J denote the lags of output and input, respectively; ½U_iand ½H_j are coefficient matrices. Notably, Eq. (10) closely resembles the time series model ARX with multiple variables. The problem of identifying a linear dynamic system is converted to a linear time series problem with one-step-ahead prediction. When all ½Hj(j ¼ 0,1,2, y, J) vanish, Eq. (10) describes free decayed vibration responses of the dynamic system.

The values of I and J in Eq. (10) can theoretically be determined in advance from the measured degrees of freedom and the true total degrees of freedom of the system under consideration if no noise is present[18,20]. However, the measured responses in reality always contain noise, and the total number degrees of freedom of the system are unknown because the system is generally continuous. Different values of I and J are often used in establishing Eq. (10). The identified modal parameters of a dynamic system usually converge to accurate results as I and J increase.

3.2. Solution in the wavelet domain

Treating {y(ti)} and {f(ti)} as vector functions, and applying the continuous wavelet transform to Eq. (10) yields

Wcyða; ¯bÞ ¼

XI i¼1

½U_iWcyða; ¯b  iÞ þ

XJ j¼0

½H_jWcf ða; ¯b  jÞ, (11)

where the translation parameter b is set to ¯b Dt, and ¯b is an integer because b must be a discrete number when the transformation is applied to discrete responses. Eq. (11) is derived using the time invariance property given

in Eq. (4), so the continuous wavelet transform on {y(t)} and {f(t)} needs to be performed only once. If the discrete wavelet transform is applied to Eq. (10) as in Huang et al.[9], then the wavelet transform has to be applied to all {y(ti)} and {f(ti)} for different i, which is time-consuming.

Constructing Eq. (11) for different ¯b and rearranging the resulting equations produces ½Yð0Þ ¼ ½ ^C  ½Y

½F " #

, (12)

where

½ ^C  ¼ ½U1 ½U2 . . . ½UI ½H0 ½H1 . . . ½HJ

h i

, (13a)

½Y ¼ ½Yh ð1ÞT ½Yð2ÞT . . . ½YðI ÞTiT, (13b)

½F ¼ ½Fh ð0ÞT ½Fð1ÞT . . . ½FðJÞTiT, (13c)

½YðiÞ ¼Wcyða; max iÞ Wcyða; max i þ 1Þ . . . Wcyða; max i þ mÞ, (13d)

½FðiÞ ¼ Wcf ða; max iÞ Wcf ða; max i þ 1Þ . . . Wcf ða; max i þ mÞ

 

(13e) and max is the larger value of I and J, so that ¯b does not become negative in Eq. (11). In Eqs. (13d) and (13e), m is a number sufficiently large to make Eq. (12) be an overdeterminate system of linear algebraic equations for unknown coefficients in ½ ^C. In the following analysis, m is set to the data number of responses under consideration.

If only one value of a utilised in building Eq. (12), then the responses and input forces in a certain range of frequency are considered. The frequency range is determined by the mother wavelet function and the chosen scale parameter a. To employ a wide frequency range in responses and input forces, several values of a can be adopted in establishing Eq. (12). In this case, the expressions of ½YðiÞand ½FðiÞare modified as

½YðiÞ ¼hYi₁ Yi₂ . . . Yi_Niand ½FðiÞ ¼hFi₁ Fi₂ . . . Fi_Ni, (14) where

Yi_j¼ Wcyðaj; max iÞ Wcyðaj; max i þ 1Þ . . . Wcyðaj; max i þ mÞ

h i

, (15a)

Fi_j¼ Wcf ðaj; max iÞ Wcf aj; max i þ 1Þ . . . Wcf ðaj; max i þ mÞ

h i

. (15b)

The coefficient matrix ½ ^C  in the overdeterminate system of Eq. (12) is calculated by the least squares or the generalised inverse approaches. Hence,

½ ^C  ¼ ½Yð0Þ ½Y ½F " #þ

, (16)

where the superscript ‘‘+’’ denotes the generalised inverse operation. 3.3. Estimation of the dynamic characteristics of structures

Apparently, Eq. (10) without the force terms describes free decay responses, explaining why the dynamic characteristics of the structure are determined from the coefficient matrices ½U_i. By adopting the concept behind the Ibrahim time domain system identification technique, Huang[21]and Yang et al.[22]proved that if

one constructs a matrix from these coefficient matrices as follows: G ½  ¼ 0 I 0 0 . . . 0 0 0 I 0 . . . 0 .. .

½UI ½UI 1 ½U1

2 6 6 6 6 4 3 7 7 7 7 5, (17)

where I is an l*l unit matrix, and l is the dimension of {y(t)}, then the dynamic characteristics of the system under consideration can be computed from the eigenvalues and eigenvectors of [G].

In Eq. (17), [G] is an ¯N  ¯N matrix, where ¯N ¼ I  l. Let lk and {fk} represent the kth eigenvalue and

eigenvector of [G], respectively. The eigenvalue lkis normally a complex number, and is set to ~akþi ~bk. The

frequency and damping ratio of the system are computed by ~bk¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a2 kþb 2 k q ; xk¼ ak= ~bk, (18) where b_k¼ 1 Dttan 1 b~k ~ak ! ; ak¼ 1 2Dtln ð ~a 2 kþ ~b 2 kÞ. (19)

In these equations, ~bk is the pseudo-undamped circular natural frequency, and xk is the modal damping

ratio.

The eigenvector {fk} can be written as ffkgT1 ffkgT2 . . . ffkgTI

 T

with each ffkgi having l components.

Since ½G fk  ¼lk fk  , (20)

substituting Eq. (17) into Eq. (20) results in

ff_kg_i¼lkffkgi1 ði ¼ 2; 3; . . . ; I Þ. (21)

The relation given in Eq. (21) yields

ff_kg ¼ ff_kgT₁; lkffkgT1; l2kffkgT1; . . . ; ln1k ffkgT1

 T

, (22)

where ff_kg1 corresponds to a mode shape of the measured degrees of freedom of the structural system.

Notably, spurious modes other than real structural modes occur when ¯N42n where n is the number of degrees of freedom of the structural system under consideration and is usually unknown. The real structural modes consistently arise as ¯N increases and when no spurious mode occurs with strong interference with the structural modes. The modal parameters of a structure can be determined from the corresponding stabilisation diagrams, which display the variations of the identified modal parameters with the values of I and J in Eq. (10). Sometimes, good engineering judgment based on the knowledge on the structural system must be exercised to confirm the structural modes.

4. Numerical verification

Numerically simulated acceleration responses of a six-story shear building subjected to base excitation were processed to demonstrate the accuracy and effectiveness of the proposed approach for identifying modal parameters. The theoretical natural frequencies of the system were 0.679, 2.18, 3.48, 6.00, 6.20 and 7.08 Hz, and the modal damping ratios were set to 5%. The numerical solutions were obtained by applying the Laplace transform to the equations of motion. Fig. 1 displays the input and response time histories of the second, fourth, and top floors and their corresponding Fourier spectra. The Fourier spectra reveal that the amplitude of the sixth mode is smaller than that of the first mode by two orders. The close modes, the 4th and 5th modes, are not easily separated in these spectra. The input base excitation and acceleration responses of the six degrees of freedom with t ¼ 0–60 s (see Fig. 1)

with a sampling rate of 250 Hz were used to estimate the corresponding wavelet transforms Wcyða; ¯bÞ and

Wcf ða; ¯bÞ where ¯b ¼ 1; 2; . . . ; 15; 000. Wcyða; ¯bÞ and Wcf ða; ¯bÞ with ¯b ¼ 200024000 were used in constructing

Eq. (12), to avoid the border problems (or edge effects) in the wavelet transforms[23].

The modal assurance criterion (MAC)[24]was applied to verify the agreement between the identified and theoretical mode shapes. The modal assurance criterion is defined as

MAC fj _iIg; fj_iAg¼ fjiIg T_fj iAg   2 fj_iIgTfj_iIgfj_iAgTfj_iAg, (23)

where fjiIgis the identified ith mode shape, and fjiAgis the corresponding analytical mode shape. The value

of MAC is between zero and unity. MAC is close to one when two mode shapes are similar, and is zero when two mode shapes are orthogonal to each other.

The following analysis considers the identified results as ‘‘accurate’’ when the identified frequencies and modal damping ratios are within 2% and 20% of the theoretical ones, respectively, and the MAC values exceed 0.9 as well. The identified natural frequencies are typically much more accurate than the identified modal damping ratios, and therefore smaller errors are required for the identified natural frequencies than for identified modal damping ratios when defining accurate results.

4.1. Shannon wavelet function

The mother wavelet function was not specified in Section 3. Eq. (5) indicates that the wavelet transform exhibits a frequency filtering effect. The wavelet transform should not significantly alter the frequency characteristics (the locations of the peaks) of the original responses in the frequency range of interest when Eq. (11) or (12) is constructed. The Shannon wavelet function [16,25], an ideal band-pass filter, seems to be a favorable choice.

The Shannon wavelet and its Fourier transform are, respectively (see Fig. 2)

c_SðtÞ ¼ sinð2ptÞ  sinðptÞð Þ=pt (24a)

and ^

c_SðoÞ ¼ 1 for po oj jo2p; 0 otherwise: 

(24b) The Fourier transform of the Shannon wavelet is compactly supported, and the wavelet decays slowly in time. Theoretically, WcSf ða; bÞ, the wavelet transform of f(t) obtained using the Shannon wavelet, completely

preserves the characteristics of f(t) in the frequency interval ½1=2a; 1=a (in terms of Hz) when only positive frequencies are considered. The frequency interval ½1=2a; 1=a is defined as the preserved frequency range, corresponding to the Shannon wavelet with a fixed scale parameter a.

4.2. Effects of noise

The measured data always contain some level of corrupted noise. To simulate the noise effect, 10% or 20% variance of the noise-to-signal ratio (NSR) was randomly added to the computed responses and input.Table 1

lists the identified modal parameters obtained from the acceleration responses of all six degrees of freedom and the base excitation input with various noise levels, using different pairs of (I, J) in Eq. (10). For simplicity and based on Eq. (9c), J is set equal to I in all the numerical results. The aim of this experiment was to identify the higher modes, that is, modes 4–6, so that a ¼ 0.092 was used for the results inTable 1.

When no noise is present, the results demonstrate that I ¼ J ¼ 2 yielded very accurate modal parameters. Raising I and J still yielded the results with very high accuracy. However, the processing of noisy responses reveals different observations. When NSR ¼ 10%, I ¼ JX10 yielded accurate results, while NSR ¼ 20% needed I ¼ JX15. Noise always causes discrepancies between the identified modal parameters and the analytical ones. Notably, spurious modes other than real structural modes occur when I42. The real structural modes are consistently identified as I increases.

Fig. 3illustrates the stabilisation diagrams of the results obtained using different values of a in processing

the responses with NSR ¼ 20%. The values of a were 1.0 and 0.25 for identifying mode 1 and modes 2–3, respectively. As before, the results of the 4th–6th modes were obtained using a ¼ 0.092. Using I ¼ JX10

yielded accurate results for the first to third modes. Modes 4–6 were more difficult to identify accurately than modes 1–3 and needed larger values of I and J.Fig. 3presents that the structural modes are stably identified by increasing I and J. Notably, the 6th mode shape obtained from I ¼ J ¼ 13 was worse than those obtained from different values of I and J because a spurious mode occurred at a frequency close to that of the 6th mode and interferes significantly with that mode.

-8 -6 -4 -2 0 2 4 6 8 -1 -0.5 0 0.5 1 ψS ψMe ψMo ψH 0 2 0 0.2 0.4 0.6 0.8 1 1.2 -8 -6 -4 -2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 0 0.1 0.2 0.3 0.4 0.5 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0 1 2 3 0 0.2 0.4 0.6 0.8 1 time(sec) -1 -0.5 0 0.5 1 0 10 frequency(Hz) time(sec) frequency(Hz) time(sec) frequency(Hz) time(sec) frequency(Hz) 0 0.2 0.4 0.6 0.8 |ψS | |ψMe | |ψMo | |ψH | 2 4 6 8 0 2 4 6 8 -8 -6 -4 -2 0 2 4 6 8 1 3 (a) (b) (c) (d)

Fig. 2. Mother wavelet functions and their Fourier moduli: (a) Shannon function, (b) Meyer function, (c) Morlet function, and (d) Haar function.

4.3. Effects of scale parameter selection

Selecting an appropriate value of the scale parameter a seems to be important. It is well known that the scale parameter affects the frequency resolution in wavelet transformation of responses. However, the time-frequency resolution affected by the scale parameter does not need to be specially considered in the present approach, which can be observed from the results shown in this section.

Table 2presents the results for the 2nd and 3rd modes obtained using various values of a and different pairs

of (I, J). The responses with NSR ¼ 20% were processed. The values of a in Table 2 yield the preserved frequency range ½1=2a; 1=a, covering the natural frequencies of the 2nd and 3rd modes. The lower bound of

Table 1

The effects of noise and mother wavelet functions on identifying the 4th–6th modes

Wavelet function Data condition I & J Natural frequency (Hz) Damping ratio (%) MAC

4 5 6 4 5 6 4 5 6 Shannon No noise 2 6.00 6.20 7.08 5.00 5.00 5.00 1.00 1.00 1.00 5 6.00 6.20 7.08 5.00 5.00 5.00 1.00 1.00 1.00 10 6.00 6.20 7.08 5.00 5.00 5.00 1.00 1.00 1.00 15 6.00 6.20 7.08 5.00 5.00 5.00 1.00 1.00 1.00 10% noise 5 6.07 6.17 7.08 4.97 6.61 5.62 0.95 0.94 1.00 10 5.99 6.20 7.06 4.85 5.11 4.62 1.00 0.99 1.00 15 5.97 6.22 7.10 4.76 4.89 4.54 0.99 0.98 0.99 20 5.98 6.22 7.11 4.63 4.97 4.51 0.99 0.97 0.99 20% noise 10 5.99 6.22 7.18 4.32 4.96 4.27 0.98 0.96 0.97 15 6.03 6.26 7.19 4.93 5.12 3.90 0.96 0.97 0.97 20 6.04 6.22 7.12 5.01 5.78 4.75 0.99 0.98 0.97 25 6.03 6.18 7.09 4.87 5.38 4.56 0.99 1.00 0.98 Meyer 10% noise 5 6.06 6.27 7.16 4.99 5.63 4.43 0.97 0.94 0.99 10 6.06 6.27 7.21 4.76 6.03 5.06 0.99 0.93 1.00 15 6.00 6.22 7.12 4.91 5.12 4.92 0.99 1.00 0.99 20 6.01 6.21 7.13 5.01 5.06 4.62 1.00 0.99 1.00 20% noise 10 6.11 6.48 7.44 5.03 8.11 6.05 0.91 0.89 0.94 15 6.02 6.26 7.14 4.70 5.60 5.11 0.98 0.99 0.96 20 6.02 6.25 7.17 4.69 5.24 4.52 0.98 1.00 0.97 25 6.02 6.23 7.19 5.04 5.41 4.13 0.99 0.95 0.95 Morlet 10% noise 5 6.00 6.23 7.12 4.31 5.22 4.63 0.99 0.96 0.97 10 6.00 6.24 7.12 4.87 5.02 4.54 0.98 1.00 1.00 15 6.02 6.23 7.13 4.93 5.20 4.38 0.99 0.99 1.00 20 6.01 6.19 7.11 5.04 5.33 4.55 0.99 0.99 0.99 20% noise 10 6.06 6.22 7.16 4.82 5.32 3.75 0.99 0.97 0.98 15 6.04 6.22 7.17 4.72 5.44 4.34 1.00 1.00 0.98 20 6.02 6.19 7.11 4.78 5.56 5.27 0.99 1.00 0.97 25 5.98 6.22 7.09 4.32 4.78 5.71 0.99 0.95 0.96 Haar 10% noise 5 6.09 5.57 0.92 10 6.09 6.28 7.31 5.39 7.53 6.08 0.96 0.69 0.99 15 6.04 6.21 7.12 5.02 5.62 5.14 1.00 0.98 1.00 20 6.02 6.21 7.11 4.88 5.26 4.63 1.00 1.00 1.00 20% noise 10 6.11 5.24 0.60 15 6.08 6.25 7.21 5.23 7.34 6.57 0.98 0.86 0.99 20 6.07 6.21 7.18 5.11 6.16 5.13 0.99 0.92 0.98 25 6.00 6.21 7.13 4.81 5.44 5.90 1.00 0.99 0.96

the frequency interval corresponding to a ¼ 0.23 only slightly exceeds the second natural frequency. Notably, when a ¼ 0.23, the modal parameters of the 3rd mode identified using I ¼ J ¼ 20 are not so good as those obtained using other I ¼ J larger than 10 because a spurious mode arose with strong interference with the 3rd mode.

Table 2reveals that although different values of a may yield markedly different results with small I and J

(I ¼ J ¼ 7, 8, or 9), all of the identified results are accurate when I and J are larger. Accordingly, as long as the natural frequencies are covered by the preserved frequency interval corresponding to the chosen a, these modes can always be accurately identified by increasing I and J.

0 2 4 6 8 frequenc y(Hz) 3 4 5 6 7 8 9 damping ratio(%) 1st mode 2nd mode 3rd mode 4th mode 5th mode 6th mode 6 8 10 12 14 16 18 20 22 24 I,J 0.6 0.8 1 MA C

As mentioned in Section 3.2, numerous values of a can be simultaneously employed to establish Eq. (12), so a wide range of frequencies can be covered and many modes can be identified at the same time.Fig. 4displays the stabilisation diagrams of the results obtained using a ¼ 1.0, 0.25, and 0.092 to process the responses with NSR ¼ 20%. Accurate results for all six modes were obtained with I ¼ J440 that is much larger than the minimum value of I and J needed for accurate results in Fig. 3.

4.4. Effects of wavelet function

The Meyer, Morlet, and Haar wavelet functions are applied to process the simulated responses with NSR ¼ 10% or 20% to explore the effects of using different mother wavelet functions on identifying the modal parameters.

The Meyer wavelet is defined in the frequency domain as [8]

^ c_Með Þ ¼o ð2pÞ1=2eio=2_sin p 2v 2p3j j o 1     _2p 3 p oj jp4p3 ; ð2pÞ1=2eio=2cosp₂v_4p3 j j o 1 4p₃ p oj jp8p₃ ; 0 j jeo 2p₃ ;8p 3   ; 8 > > < > > : (25)

where vðsÞ ¼ s4_{ð35  84s þ 70s}2_20s3_{Þand s 2 ½0; 1. No explicit analytical expression exists for the Meyer}

wavelet in the time domain. The Meyer wavelet and its Fourier transform modulus are also depicted inFig. 2. Apparently, the Meyer function is compactly supported in the frequency domain and smoothes out the discontinuity of the Shannon wavelet function. The Meyer wavelet is not compactly supported in the time domain and decays faster than the Shannon wavelet, as time goes to infinity.

The Morlet wavelet and its Fourier transform are, respectively,[14]

c_MoðtÞ ¼ et2=ð2d2Þþi ¯bt _(26a)

and ^

c_MoðoÞ ¼ dpffiffiffiffiffiffi2peðo ¯bÞ2d2=2, (26b)

where d and ¯b are constants. The Morlet wavelet does not have compact support in the time and frequency domains, but decays more quickly than the Shannon and Meyer wavelets as time approaches infinity. Since

^

c_Moð0Þa0, the Morlet wavelet does not fulfill the admissibility condition given in Eq. (2). However, the Morlet wavelet is numerically admissible when d ¯b is large enough.

A d ¯b value of larger than 5 is adopted in practice. In the literature on applying the continuous wavelet transform to system identification[10–14], d ¼ 1 was commonly used to provide a high frequency resolution. However, d ¼ 1 also causes a very sharp peak of ^_c_MoðoÞ_{ at o ¼ ¯b, which may yield the wavelet transform} significantly altering the frequency characteristics (locations of the peaks) of a function in the frequency range

Table 2

The effects of scale parameter selection on identifying 2nd and 3rd modes

I ¼ J Natural frequency (Hz) Damping ratio (%) MAC

0.23 0.25 0.27 0.23 0.25 0.27 0.23 0.25 0.27 7 3.47 2.17 3.49 2.17 3.50 9.00 7.86 7.51 6.93 8.92 0.99 1.00 1.00 1.00 0.99 8 2.09 3.49 2.18 3.50 2.17 3.48 7.61 6.64 6.56 6.59 6.51 7.92 1.00 1.00 1.00 1.00 1.00 0.99 9 2.15 3.49 2.17 3.49 2.17 3.48 5.40 5.54 6.18 5.58 5.37 6.30 1.00 1.00 1.00 1.00 1.00 1.00 10 2.18 3.49 2.16 3.49 2.18 3.49 4.62 5.55 5.35 5.32 4.85 5.59 1.00 1.00 1.00 0.99 1.00 1.00 15 2.19 3.49 2.17 3.48 2.18 3.49 5.26 5.15 4.65 5.06 4.86 5.14 0.99 0.98 0.99 0.99 1.00 0.99 20 2.19 3.44 2.18 3.48 2.18 3.49 4.95 5.74 4.65 4.92 4.87 5.10 0.99 0.92 1.00 0.98 1.00 1.00 25 2.19 3.49 2.18 3.48 2.18 3.48 4.63 5.04 4.65 4.85 4.82 5.03 1.00 0.97 1.00 0.98 1.00 0.99 30 2.18 3.50 2.19 3.50 2.18 3.48 4.82 5.50 4.86 5.14 4.83 4.86 1.00 0.99 1.00 0.96 1.00 0.98

of interest. In the following analysis, d ¼ 1=pffiffiffi2 and ¯b ¼ 2:26p were used, and the corresponding Morlet wavelet and its Fourier transform modulus are also shown inFig. 2.

The Haar wavelet function is the simplest available wavelet function. The Haar function and its Fourier transform are, respectively,

c_HðtÞ ¼ 1 0ptp1 2; 1 1 2ptp1; 0 other 8 > < > : (27a) 0 2 4 6 8 frequenc y(Hz) 4 5 6 7 8 9 damping ratio(%) 1st mode 2nd mode 3rd mode 4th mode 5th mode 6th mode 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 I,J 0.6 0.8 1 MA C

and ^

c_HðoÞ ¼ 1 ioð1 þ e

io_2eio=2_{Þ. (27b)}

The Haar wavelet and its Fourier transform modulus are also given inFig. 2. The Haar wavelet has the shortest support in time among all orthogonal wavelets.

The Meyer, Morlet, and Haar wavelets are not ideal band-pass filters. In particular, j ^cHðoÞj has numerous

side lobes, so that the Haar wavelet function is far away from an ideal band-pass filter. The preserved frequency range for these wavelets is defined simply as the frequency interval in which the Fourier transform modulus of the wavelet exceeds 80% of its maximum value. Consequently, the preserved frequency ranges for the Meyer, Morlet, and Haar wavelets are [0.5348/a, 0.9311/a], [0.9875/a, 1.275/a], and [0.4313/a, 1.0807/a], respectively.

The values of a were set to 0.1, 0.12, and 0.08 for the Meyer, Morlet, and Haar wavelets, respectively, to identify the modal parameters of the 4th–6th modes from the responses with NSR ¼ 10% and 20%. The results obtained with different pairs of I and J are also given inTable 1.Table 1reveals that the used wavelet mother function considerably affects the accuracy of identified modal parameters obtained using small I and J (smaller than 10). Nevertheless, using sufficiently high values of I and J yields accurate results, regardless of which wavelet function is employed.

Figs. 5–7 illustrate the stabilisation diagrams of the identified results obtained using different wavelet

functions in processing the responses with NSR ¼ 20%. To identify the modal parameters of mode one, a was set equal to 1.0 and 0.9 for the Meyer and Haar wavelets, respectively, while a was set equal to 0.26 and 0.22 for modes 2 and 3. When the Morlet wavelet function was used, a was set to 1.5, 0.50, and 0.32 for identifying modes 1–3, respectively, because a suitable value of a was not found to make the corresponding preserved frequency range cover the natural frequencies of modes 2 and 3 simultaneously. Notably, the occurrence of spurious modes near mode 6 worsens the accuracy of the identified modal damping of mode 6 obtained using I ¼ J ¼ 24 and 25 when the Meyer wavelet was used (see Fig. 5). These results ofFigs. 5–7further confirm that using sufficiently high values of I and J yields accurate results, regardless of which wavelet function is employed.

Figs. 3 and 5–7and Table 1indicate that the Haar function is not as efficient as the Shannon, Meyer or

Morlet function in identifying the modal parameters because the Haar function needs larger I and J to give accurate results. Using larger I and J increases computational time. Comparing the results obtained using the Shannon, Meyer, and Morlet functions reveals that the Morlet function can be utilised with smaller I and J values to yield accurate results. Notably, without showing results here, these four wavelet functions exhibited the same efficiency in processing the data without noise, and the modal parameters were identified precisely using I ¼ J ¼ 2.

Finally, the preserved frequency range is the frequency interval in which most of the frequency contents of f(t), representing a series of signals, are preserved after the wavelet transformation. Although the value of 80% chosen in defining the preserved frequency range may be somewhat empirical, the definition is simple, and such preserved frequency range is very easy to be determined and works well for identifying modal parameters. The value of 80% used is not crucial in identifying the modal parameters. Different values can be adopted to define the preserved frequency range if j ^Wcf ða; oÞj and j ^f ðoÞj have peaks at the possible natural frequencies in

the preserved frequency range. Then, the preserved frequency ranges defined by different percentages of the maximum Fourier transform modulus of the wavelet may yield different values of the scale parameter for identifying the modal parameters of certain modes. The effects of using different values of a on the accuracy of the identified results were studied in the preceding section.

5. Application

The measured responses of a three-story non-symmetric steel frame in shaking table tests (seeFig. 8) were processed to demonstrate the applicability of the present approach to real measured data. Shaking table tests are often carried out in a laboratory to examine the behavior of structures in earthquakes. The three-story steel frame under consideration was 3 m long, 2 m wide, and 9 m high. Lead blocks were piled on each floor,

such that the mass of each floor was approximately 3500 kg. The frame was subjected to base excitation of the El Centro earthquake with a reduced level. Twelve accelerometers were installed on the floors, as indicated in

Fig. 9(i.e. Y3b, Y2b, X3a, X2a, etc.), to measure the acceleration responses of the frame at different locations

and directions. The accelerometers were placed near the corners of floors. The data sampling rate was 200 Hz. As shown inFig. 8, the asymmetry of the frame was mainly due to the incorporation of three stiffening braces into each floor. The asymmetry led to the coupling between the responses in long axis and short axis directions, called the X and Y directions, respectively. The responses measured by the 12 accelerometers, representing 12 degrees of freedom, were used in the following analysis.

Fig. 10depicts the base excitation input in the X direction and some of the responses at different degrees of

freedom, as well as the corresponding Fourier spectra. Since the base excitation was in the X direction, the 0 2 4 6 8 frequenc y(Hz) 4 5 6 7 8 9 damping ratio(%) 1st mode 2nd mode 3rd mode 4th mode 5th mode 6th mode 6 8 10 12 14 16 18 20 22 24 I,J 0.6 0.8 1 MA C

responses in the X direction were larger than those in the Y direction. The Fourier spectra of the responses reveal peaks in the frequency ranges of 3–5, 10–14, 15–18, 21–23, 29–31 and 35–38 Hz, based on which appropriate values of a were chosen for identifying the modal parameters of the frame.

The responses at t ¼ 0–40 s (see Fig. 10) were used in the wavelet transform with the Morlet wavelet function, and only Wcyða; ¯bÞ and Wcf ða; ¯bÞ with ¯b ¼ 100022000 were used in constructing Eq. (12). Different

values of a were adopted to identify different modes. The values of a were set to 0.35, 0.25, 0.095, 0.064, 0.05, and 0.033 to identify the natural frequencies in [2.82,3.64], [3.95,5.10], [10.4,13.4], [15.4,19.9], [19.8,25.5], and [29.9,38.6] Hz, respectively.Fig. 11depicts the stabilisation diagrams of the identified natural frequencies and modal damping ratios. The modal parameters were stably identified by using large values of I and J.

0 2 4 6 8 frequenc y(Hz) 3 4 5 6 7 8 9 damping ratio(%) 1st mode 2nd mode 3rd mode 4th mode 5th mode 6th mode 6 8 10 12 14 16 18 20 22 24 I,J 0.6 0.8 1 MA C

Accurately identifying the parameters for higher modes (5th–9th modes) needs I and J values larger than those for identifying the lower modes.

Table 3lists the identified natural frequencies and modal damping ratios, which are the averages of the

results obtained using the largest five I and J values for each mode shown inFig. 11, while the mode shapes are depicted inFig. 12. The floor was assumed to be rigid when plotting these mode shapes. Nine modes were identified in total. InFig. 12, the black dotted, black solid, gray solid and gray dotted lines are corresponding to columns 1, 2, 3 and 4 inFig. 9, respectively.Fig. 12shows that the first and second modes were mainly transverse modes in the X and Y directions, respectively, while the third mode was a torsional mode. The coupling of the responses in different directions was very severe in the higher modes.

0 2 4 6 8 frequenc y (Hz) 4 5 6 7 8 9 10 damping ratio(%) 1st mode 2nd mode 3rd mode 4th mode 5th mode 6th mode 10 15 20 25 30 35 I,J 0.7 0.8 0.9 1 1.1 MA C

6. Concluding remarks

This work presented an approach based on the continuous wavelet transform for identifying modal parameters of a linear system from its seismic responses or free vibration responses. The continuous wavelet transform is applied to the discrete equations of motion among the measured degrees of freedom. The transforms of the functions in the equations are calculated from the measured responses and using the time invariance property of the continuous wavelet transform. By choosing appropriate scale parameters in the transformation, data filtering is also carried out in the transformation to enhance the efficiency of identifying the modal parameters (natural frequencies, modal damping ratios, and mode shapes), which are estimated directly from the coefficients of the discrete equations of motion.

The main theoretical advantage of the proposed approach over other existing approaches based on the continuous wavelet transform is that its measured responses processed are not limited to free decay responses. Unlike the existing approaches based on the discrete wavelet transform, the proposed approach employs the time invariance property and filtering ability of the continuous wavelet transform to improve data processing efficiency.

The proposed approach was validated first by successfully processing the numerically simulated acceleration responses of a six-story shear building subjected to base excitation. The modal parameters of all the degrees of

freedom were easily and accurately identified even when the responses and input included noise with NSR ¼ 20%. The Shannon, Meyer, Morlet, and Haar wavelet functions were applied to process these responses. The wavelet function and the chosen scale parameter can significantly affect the accuracy of the identified modal parameters obtained using small values of I and J in Eq. (10). Nevertheless, using large values of I and J always yields accurate results, regardless of which wavelet function and the scale parameter are employed if the corresponding preserved frequency range covers the frequencies of interest.

2 3 4 5 frequenc y(Hz) 6 8 10 12 14 16 18 20 I,J 3 4 5 6 7 damping ratio(%) frequenc y(Hz) damping ratio(%) frequenc y(Hz) damping ratio(%) frequenc y(Hz) damping ratio(%) 1st mode 2nd mode 9 10 11 12 13 0 1 2 3 4 3rd mode 4th mode 16 18 20 22 24 10 12 14 16 18 20 22 24 26 28 30 32 I,J 6 8 10 12 14 16 18 20 I,J 10 12 14 16 18 20 22 24 26 28 30 32 I,J 0 1 2 3 _{5th mode} 6th mode 26 28 30 32 34 36 38 40 0 2 4 6 8 10 7th mode 8th mode 9th mode

Fig. 11. Stabilisation diagrams of the identified modal parameters for a three-story steel frame.

Table 3

Identified frequencies and modal damping ratios

Natural frequency (Hz) Damping ratio (%)

3.04 4.01 4.23 5.57 10.50 3.31 12.97 0.61 17.38 0.89 21.97 0.61 30.35 1.42 31.84 1.33 37.28 1.25

The measured responses of a three-story non-symmetric steel frame with 3 m long, 2 m wide, and 9 m high from shaking table tests were analysed to demonstrate the applicability of the present approach in processing real measured data. A total of nine coupling modes were identified.

Although this work only demonstrates the feasibility of the proposed approach in processing the responses of a linear system subjected to base excitation, the proposed approach is certainly suitable for dealing with free decay responses or the responses due to different types of force input if the input forces are measured. Applying the proposed approach to determine the modal parameters of a system from its ambient

measurements requires a technique such as the random decrement technique to convert the random responses into free decay responses.

Acknowledgements

This work reported herein was supported by the National Science Council, ROC through research grant No. NSC92-2211-E-009-043. The support is gratefully acknowledged. The appreciation is also extended to the National Center for Research on Earthquake Engineering for providing shaking table test data.

References

[1] J.S. Bendat, A.G. Piersol, Engineering Applications of Correlation and Spectral Analysis, second ed., Wiley, New York, 1993. [2] E. Safak, M. Celebi, Seismic response of Transamerica building. II: system identification, Journal of Structural Engineering, ASCE

117 (8) (1991) 2405–2425.

[3] A. VanDerVeen, E.F. Deprettere, A.L. Swindlehurst, Subspace-based signal analysis using singular value decomposition, Proceedings of the IEEE 81 (9) (1993) 1277–1308.

[4] C.S. Huang, H.L. Lin, Modal identification of structures from ambient vibration, free vibration, and seismic response data via a subspace approach, Earthquake Engineering and Structural Dynamics 30 (12) (2001) 1857–1878.

[5] C.S. Huang, S.L. Hung, C.M. Wen, T.T. Tu, A neural network approach for structural identification and diagnosis of a building from seismic response dat, Earthquake Engineering and Structural Dynamics 32 (2003) 187–206.

[6] A.N. Robertson, K.C. Park, K.F. Alvin, Extraction of impulse response data via wavelet transform for structural system identification, Journal of Vibration and Acoustics, Transactions of the ASME 120 (1) (1998) 252–260.

[7] A.N. Robertson, K.C. Park, K.F. Alvin, Identification of structural dynamics models using wavelet-generated impulse response data, Journal of Vibration and Acoustics, Transactions of the ASME 120 (1) (1998) 261–266.

[8] I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia, 1992.

[9] C.S. Huang, S.L. Hung, C.I. Lin, W.C. Su, A wavelet-based approach to identifying structural modal parameters from seismic response and free vibration data, Computer-Aided Civil and Infrastructure Engineering 20 (2005) 408–423.

[10] D.A. Schoenwald, System identification using a wavelet-based approach, in: Proceedings of the 32nd IEEE Conference on Decision and Control, San Antionio, TX, 1993, pp. 3064–3065.

[11] M. Ruzzene, A. Fasana, L. Garibaldi, B. Piombo, Natural frequencies and dampings identification using wavelet transform: application to real data, Mechanical System and Signal Processing 11 (2) (1997) 207–218.

[12] S. Gouttebroze, J. Lardies, On using the wavelet transform in modal analysis, Mechanics Research Communications 28 (5) (2001) 561–569.

[13] J. Lardies, S. Gouttebroze, Identification of modal parameters using the wavelet transform, International Journal of Mechanical Sciences 44 (11) (2002) 2263–2283.

[14] T.-P. Le, P. Argoul, Continuous wavelet transform for modal identification using free decay response, Journal of Sound and Vibration 277 (2004) 73–100.

[15] C.K. Chui, An Introduction to Wavelets. Wavelet Analysis and its Application, vol. 1, Academic Press, Boston, 1992. [16] A. Teolis, Computational Signal Processing with Wavelets, Birkhauser, Boston, 1998.

[17] L.R. Rabiner, B. Gold, Theory and Application of Digital Signal Processing, Prentice-Hall, Inc, Englewood Cliffs, NJ, 1975. [18] J.M. Leuridan, Some direct parameters model identification Methods applicable for multiple input modal analysis, Ph.D.

Dissertation, Department of Mechanical and Industrial Engineering, University of Cincinnati, USA, 1984.

[19] V. Sotoudeh, Measurement and analysis of structural vibrations. Ph.D. Dissertation, Department of Civil Engineering, Stanford University, USA, 1986.

[20] C.S. Huang, Structural identification from ambient vibration measurement using the multivariate AR model, Journal of Sound and Vibration 241 (3) (2001) 337–359.

[21] C.S. Huang, A study on techniques for analyzing ambient vibration measurement (II)—time series methods, Report No. NCREE-99-018, National Center for Research on Earthquake Engineering, Taipei, Taiwan, 1999 (in Chinese).

[22] Q.J. Yang, P.Q. Zhang, C.Q. Liand, X.P.A. Wu, System theory approach to multi-input multi-output modal parameters identification method, Mechanical Systems and Signal Processing 8 (1994) 159–174.

[23] S. Mallat, A Wavelet Tour of Signal Processing, Academic Press, London, 1998.

[24] R.L. Allemang, D.L. Brown, A correlation coefficient for modal vector analysis, in: Proceeding of the First International Modal Analysis Conference, Bethel, CN, USA, 1983, pp. 110–116.