Model order determination and noise removal for modal parameter estimation

ARTICLE IN PRESS Mechanical Systems and Signal Processing 24 (2010) 1605–1620

Contents lists available at ScienceDirect

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

Model order determination and noise removal for modal parameter estimation Sau-Lon James Hu a,, Xingxian Bao b, Huajun Li b a b

Department of Ocean Engineering, University of Rhode Island, Narragansett, RI 02882-1197, USA Department of Ocean Engineering, Ocean University of China, Qingdao 266100, PR China

a r t i c l e in fo

abstract

Article history: Received 16 September 2009 Received in revised form 21 October 2009 Accepted 13 January 2010 Available online 22 January 2010

Given a noisy impulsive response function (IRF) that has been contributed by an unknown number of modes, this article proposes a different approach from the traditional methods for estimating modal parameters from this noisy IRF. The major difference lies in the way of handling noise and choosing the computational model order. Whereas the traditional approach accommodates noise by purposely increasing the computational model order, the proposed approach uses the actual system order as the computational model order and rejects noise prior to performing the modal parameter estimation. The proposed approach includes three steps: (1) model order (or number of modes) determination from the measured IRF—by finding the rank of a Hankel matrix constructed from the measured IRF, (2) noise removal from the measured IRF to obtain a filtered IRF—by implementing Cadzow’s algorithm for the structured low rank approximation (SLRA) on the Hankel matrix, and (3) modal parameters estimation from the filtered IRF—by using the complex exponential method (Prony’s method). Numerical studies include both synthesized and experimental data. While measured IRFs with mild and strong noise levels are simulated for a 5 degree-of-freedom massspring-dashpot system, the modal parameter estimations based on the filtered IRFs are very good for both noise levels. While experimental data are measured from two accelerometers mounted at a cantilever beam, the modal parameters estimated from the filtered IRFs of the two accelerometers are in excellent agreement. & 2010 Elsevier Ltd. All rights reserved.

Keywords: Model order determination Noise removal Modal parameter estimation

1. Introduction Modal identification involves estimating the modal parameters of a structural system from measured input–output data [1,2]. Over the past thirty years, many algorithms have been developed to estimate modal parameters based upon the measured data being the frequency response function (FRF) or the equivalent impulse response function (IRF). Because modal identification algorithms require the prior knowledge of the model order (or number of modes), much of the work concerned with modal parameter estimation since late seventies has involved methodology for determining a right model order for the modal parameter model. Today, an appropriate estimate of the order of the model has become the single most important problem in modal parameter estimation [3]. Furthermore, measured signals are inevitably contaminated with noise when a data acquisition system is used for an experimental measurement, and the measured FRF and IRF may not be clean enough for estimating the modal parameters with proper accuracy.

 Corresponding author. Tel.: + 1 401 874 6688; fax: + 1 401 874 6837.

E-mail addresses: [email protected] (S.-L. James Hu), [email protected] (X. Bao), [email protected] (H. Li). 0888-3270/$ - see front matter & 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2010.01.005

ARTICLE IN PRESS 1606

S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

Traditional ways of handling the uncertainty of the model order are basically a trial-and-error procedure by repeating the same analysis many times with a varying model order until a stable result emerging. Most of the techniques are to establish a maximum model order first, then data are acquired based upon an assumption that the model order is equal to this maximum. In a sequential fashion, the data is evaluated to determine if a model order less than the maximum will describe the data sufficiently. The commonly used techniques include: measurement synthesis and comparison (curve-fit), error chart, stability diagram, etc. [3]. As the measured data are contaminated by noise, the estimated modal parameters must contain error even if the model order has been chosen correctly. In the modern modal analysis, it does not seem to have a particular procedure intending to eliminating the noise. In fact, the computational model order is often set higher than the true system order, so the socalled ‘‘noise modes’’ could be absorbed [1–6]. The inclusion of a larger number than the critical number of model order will cause the creation of ‘‘computational’’ modes in addition to the genuine ‘‘physical’’ modes which are of interest. Among many researchers who advocated the usage of an over-determined system (i.e., one whose order is higher than the physical system’s order) for improving the accuracy of parameter identification in noisy situations, Braun and Ram [4] described the strategy to determine the necessary extent of overdetermination and the procedure for distinguishing true system modes; in particular, an efficient perturbation method based on the singular value decomposition was presented. They also demonstrated that such an overdetermination is both possible and useful when working in the frequency domain [5,6]. In contrast to the classical ways – choosing a higher model order intentionally for absorbing noise – in the modal parameter estimation process, our proposed procedure is to firmly determine the model order and eliminate much noise from the measured IRF before applying a time domain technique for modal parameter estimation. In determining the model order, the first step is to construct a Hankel matrix based on the measured IRF. We then apply singular value decomposition (SVD), together with an observed noise threshold, to determine the rank of the Hankel matrix [7]. This rank in theory is equal to two times the number of modes embedded in the measured IRF. For removing noise from the measured IRF to yield a filtered IRF with a known model order, a structured low rank approximation (SLRA) method for the Hankel matrix must be carried out. Mathematically, the SLRA concerns the construction of the nearest approximation to a given matrix by a matrix with a specific rank and a specific linear structure [8]. In our case, the resulting matrix must maintain the Hankel structure, and be lowered in rank as well. The data matrix being Hankel structured is equivalent to the existence of a linear time-invariant system that fits the data and the rank constraint is related to a bound on the model order [9,10]. Once the filtered IRF is obtained, we can apply the complex exponential method (Prony’s method) to extract the modal parameters [1,2,4]. Several SLRA algorithms have been published, mainly in the journals related to signal processing and numerical linear algebra [9–15]. Most SLRA algorithms were treated as an optimization problem and emphasized on meeting a normoptimality criterion. For example, De Moor [11] treated the problem in the context of the so-called structured total least squares problem (STLS), and his approach was based on Lagrange multipliers. Another algorithm for solving the structured total least squares problem, called structured total least norm, has been proposed by Rosen et al. [12]. One major restriction of the solution algorithms cited above is that they aimed at rank reduction of the data matrix by one. Due to the limitation of the rank reduction by one, those algorithms are not attractive to our application. One simple algorithm does not have the reduction limitation in rank number is Cadzow’s algorithm [16], which is to be employed in this study. Cadzow’s algorithm is an iterative method that alternates between unstructured low rank approximation (truncated singular value decomposition) and structure enforcement via anti-subdiagonal averaging, thereby only requiring computation related to singular value decomposition and manipulation of the matrix entries. Synthesized measurements of a 5 degree-of-freedom (DOF) mass-spring-dashpot system will be used to demonstrate the performance, and illustrate the procedure as well, of the proposed scheme in the numerical studies. Corrupted IRFs with both mild and strong noise levels are to be investigated. The performance of the proposed method will also be investigated using real measurements from a cantilever beam. 2. Preliminaries The focus of this study is on employing the complex exponential method (Prony’s method) for estimating modal parameters from the measured IRFs of multi-degree-of-freedom (MDOF) systems [1,2]. The basic principle of the complex exponential method is that any free vibration response function, including IRF, can be expressed by a series of complex exponential components, and the properties of each of which contain the eigenvalue and eigenvector properties of one mode. For facilitating our later presentation, a brief description of the IRF of an N-degree-of-freedom (N-DOF) system and the essential material of Prony’s method are given below. 2.1. Impulse response function (IRF) An N-DOF dynamic system can be represented by a linear second-order differential equation as Mx€ þ Cx_ þ Kx ¼ 0

ð1Þ

ARTICLE IN PRESS S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

1607

where M, C and K 2 RNN are mass, damping and stiffness matrices, respectively; and x, x_ and x€ 2 RN the displacement, velocity, and acceleration vectors, respectively. A general expression for any frequency response function (FRF) of the system may be written as [1,2] 1 0 FðoÞ ¼

N B X B @

C Ak Ak   qffiffiffiffiffiffiffiffiffiffiffiffi þ qffiffiffiffiffiffiffiffiffiffiffiffiC A 2 2 ok xk þ i ook 1xk ok xk þ i o þ ok 1xk

k¼1

ð2Þ

where ok and xk are the natural frequency and damping ratio, respectively, of the kth mode, Ak a complex amplitude, and the superscript ‘‘*’’ denotes the complex conjugate operator. Eq. (2) can also be written as FðoÞ ¼

2N X

Ak

k¼1

ok xk þiðook0 Þ

ð3Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 where ok0 ¼ ok ð1xk Þ is the damped 0frequency, and ok þ N ¼ ok0 and Ak þ N ¼ Ak for k= 1, y, N. By taking the inverse Fourier transform of Eq. (3), the corresponding impulse response function (IRF) is obtained as hðtÞ ¼

2N X

Ak esk t

ð4Þ

k¼1

where sk ¼ ok xk þiok0 .

2.2. Discrete IRF If the original FRF has been obtained in a discrete form, and is thus described at each of a number of equally spaced frequencies, the resulting IRF (found via the inverse discrete Fourier transform of the FRF) will similarly be described at a corresponding number of equally spaced time intervals. We may conveniently express this data set as h‘  hð‘DtÞ;

‘ ¼ 0; 1; 2; . . .

ð5Þ

where Dt is the time interval. From Eqs. (4) and (5), we have h‘ ¼

2N X

Ak Vk‘

ð6Þ

k¼1

where Vk ¼ esk Dt . The above equation is a non-linear function of the unknown modal parameters. However, once Vk (i.e. modal frequencies and damping ratios) are known, Eq. (6) becomes linear with respect to the remaining unknown modal parameters Ak. For this reason, the numerical approach in many algorithms involves identifying first the modal frequencies and damping ratios, which are global properties of the system and share the same values in different IRFs, regardless of the input and output locations. Because an N-DOF dynamic system (with N simultaneous second-order linear differential equations) is mathematically equivalent to a linear differential equation of order 2N, the IRF sequence {h‘} must satisfy a linear difference equation [1,2,17]: h2N þ n þ

2N1 X

bm hm þ n ¼ 0;

n ¼ 0; 1; . . .

ð7Þ

m¼0

where bm, m=0, y, 2N 1 are real constant coefficients.

2.3. Prony’s method Similar to the discrete Fourier transform, the original Prony’s method extracts valuable information from a uniformly sampled signal and builds a series of complex exponentials. The key to Prony’s method is that the coefficients in Eq. (7) are related to the following polynomial [4,18]: 2N X m¼0

bm V m ¼

2N Y

ðVVk Þ

ð8Þ

k¼1

where b2N has been set equal to 1, and the roots of the polynomial are V1, V2, y, V2N. Thus, finding values of the bm coefficients can determine values of Vk, and hence the system natural frequencies and damping ratios.

ARTICLE IN PRESS 1608

S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

In solving bm from the IRF data, the original Prony’s method has been based on the successive applications of Eq. (7) for 2N times. Taking n = 0, y, 2N 1 with Eq. (7) results in a full set of 2N equations 8 9 9 38 2 h1 h2    h2N1 > h0 h2N > > > > b0 > > > > > > > > > 7> 6 < h2N þ 1 > = = h2 h3    h2N 7< b1 6 h1 7 6 ¼  ð9Þ 7> ^ 6^ ^ ^    ^ > > > ^ > > > 5> 4 > > > > > > > > :h ; ; h2N1 h2N h2N þ 1    h4N2 : b 2N1

4N1

or Hb ¼ h~

ð10Þ

where H 2 R2N2N is a square matrix with Hankel structure, which is referred to as a matrix with constant (positive sloping) skew-diagonals. From Eq. (10), we can obtain the unknown polynominal coefficients

b ¼ H1 h~

ð11Þ

If more than 2N equations have been included in Eq. (9), then a least-squares solution of b for an over-determined case could be obtained.

3. Hankel matrix and structured low rank approximation Let {h‘}, ‘ = 0, y, s, represent an IRF (or a segment of an IRF) sequence of an N-DOF system, and be sequentially filled to form a rectangular Hankel matrix, Hmn 2 Rmn ,with m,n Z2N as 3 2 h1    hn1 h0 7 6 h2    hn 7 6 h1 7 ð12Þ Hmn ¼ 6 7 6^ ^ & ^ 5 4 hm1 hm    hs where s= m+ n  2. The original sequence {h‘} is readily retrieved from the first row and last column, or the first column and last row, of matrix H. If {h‘} is an exact (noise-free) impulse response function (IRF) of an N-DOF dynamic system, then the rank of H should be equal to 2N, as long as m and n are both greater than or equal to 2N and all N modes contribute to the IRF. The above statement could be illuminated based on Eq. (7) as we rewrite it into a linear recursive form 2N1 X

h2N þ n ¼ 

bm hm þ n ;

n ¼ 0; 1; . . .

ð13Þ

m¼0

This recursive equation suggests that each term of the sequence is defined as a linear combination of the preceding 2N terms. Thus a discrete IRF of an N-DOF dynamic system can be fully described by any of its 2N consecutive values. Because only 2N independent values are possible in a discrete IRF, the rank of H constructed from {h‘} cannot exceed 2N.

3.1. Rank and singular value decomposition One common way of knowing the rank of a matrix is based on the singular value decomposition (SVD) of the matrix. Let the SVD of a matrix A 2 Rmm be expressed as [19] A ¼ URVT

ð14Þ

mm

T

nn

mn

where U 2 R and V 2 R are orthonormal matrices, and R 2 R a rectangular diagonal matrix whose diagonal elements are singular values, arranged in order of deceasing size. Because some of the singular values may be zero, we therefore partition R into a submatrix Rr of r nonzero singular values and several zero matrices as 

R¼

Rr

0

0

0

 ð15Þ

where Rr is a r  r diagonal matrix. The decomposition then becomes A ¼ URVT ¼ Ur Rr VTr , where Ur and Vr consist of the first r columns of U and V. This decomposition indicates the rank of A is r. For theoretical data, the singular values should go to zero when the rank of the matrix is exceeded. For measured data, however, due to random errors and small inconsistencies in the data, the singular values will not become zero but will become very small.

ARTICLE IN PRESS S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

1609

3.2. Corrupted IRF and its Hankel matrix When a measured IRF sequence {h‘} is contaminated by random noise, we can express it as h‘ ¼ h ‘ þ e‘ ¼

2M X

Ak esk ‘Dt þ e‘

‘ ¼ 0; 1; 2;   

ð16Þ

k¼1

where h ‘ and e‘ represent the true signal and noise, respectively, and M is the number of degrees of freedom that is necessary to represent the measured data. Because some modes might not contribute to this particular IRF, M is possible less than N, the number of degrees of freedom of the system. Theoretically, the Hankel matrix H corresponding to the noisy signal {h‘} in Eq. (16) can be partitioned into two parts Hmn ¼ H mn þ Emn

ð17Þ

where H and E represent Hankel-structured matrices associated with the uncontaminated signal and noise parts, respectively. Because the signal fh ‘ g is contributed by M modes, its corresponding Hankel matrix H must be with a rank 2M. A noise removal method for estimating fh ‘ g from {h‘} is thus by estimating H based on H, where the estimate of H, ^ must be a Hankel matrix with rank 2M. While approximating H to its ‘‘nearest’’ matrix H, ^ an often used denoted H, ^ ^ 2 [8]. criterion is based on minimizing the Frobenius norm (L2-norm) of the difference between H and H, denoted JHHJ 3.3. Unstructured low rank approximation The (unstructured) low rank approximation problem is to approximate a data matrix A with another matrix A^ that has a ^ 2 , and the unstructured matrix A^ can be obtained via the truncation of the SVD of specific lower rank by minimizing JAAJ the data matrix A. This approach has often been referred to as based on the truncated singular value decomposition (TSVD) technique or Eckart–Young theorem [20], which is summarized below: Eckart–Young theorem. Given an A matrix of dimensions m  n of rank R where Rrmin(m,n) and its singular value decomposition A ¼ URVT

ð18Þ

^ 2 when then there exists an m  n matrix A^ of rank r oR which minimizes JAAJ A^ ¼ URr VT

ð19Þ

where Rr is the same matrix as R except that it contains only the r largest singular values, i.e., the other singular values are replaced by zeros. 3.4. Cadzow’s algorithm The optimal L2-norm lower rank approximation to a given Hankel matrix H can be obtained from the truncated singular value decomposition (TSVD) of H, but the resulting matrix will not possess a Hankel structure. In mathematical terms, the noise reduction problem here is a structured low rank approximation (SLRA) problem. In addition to lower its rank, the resulting matrix must maintain a Hankel structure. A simple engineering algorithm for the restoration of a Hankel matrix is the iterative Cadzow’s algorithm, which has been also called the anti-diagonal averaging method [16,21]. It can be summarized by the following steps: (1) Use the TSVD technique with an appropriate value of rank r estimated from the data to obtain a low rank approximation to the Hankel data matrix. Note that the resulting matrix A^ will not be a Hankel matrix. ^ from A^ by replacing all elements of each anti-subdiagonal by the arithmetic average of the (2) Rebuild a Hankel matrix H ^ will not be rank r. elements along the anti-subdiagonal. Note that the resulting Hankel matrix H In Cadzow’s algorithm, the low rank approximation and anti-subdiagonal averaging are alternated iteratively until a convergence test has been met. It has been proven that such iterative usage always converges, but Cadzow’s algorithm has been criticized being suboptimal in L2-norm criterion [9,21]. Because there is no practical reason to suggest that H in Eq.

x1 c1 k1

x2 c2

m1

k2

x3 c3

m2

k3

x4 c4

m3

k4

Fig. 1. A 5-DOF mass-spring-dashpot system.

x5 c5

m4

k5

m5

ARTICLE IN PRESS 1610

S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

Fig. 2. Receptance FRFs of the 5-DOF system.

−5

x 10

10 IRF noise−free

FRF with 5% noise FRF noise−free

Magnitude

Impulse Response

−5

2 1 0 −1

−10

10

−2 0

0.05

0.1

0.15 0.2 0.25 Time [s]

0.3

0.35

0

0.4

50

100 150 Frequency [Hz]

200

250

x 10

200 5% noise

Phase [°]

Noise magnitude

−7

2 1 0 −1 −2 −3

FRF with 5% noise FRF noise−free

100 0 −100 −200

0

0.05

0.1

0.15

0.2 0.25 Time [s]

0.3

0.35

0.4

0

50

100 150 Frequency [Hz]

200

250

Fig. 3. (a) Time-domain signals of an exact IRF and a 5% white noise simulation and (b) the corresponding corrupted and noise-free FRFs in frequency domain.

(17) must be the nearest matrix to H in L2-norm sense, the sub-optimality in L2-norm is not a critical flaw in engineering applications. 4. Simulation study A 5-DOF mass-spring-dashpot system shown in Fig. 1 is chosen for the numerical example to demonstrate the performance, and to illustrate the procedure as well, of the proposed scheme. Uniform mass, stiffness and damping coefficients are taken to be mn = 50 kg, kn = 2.9  107 N/m, and cn = 1000 N s/m, respectively, for n =1, y, 4. The coordinates of the 5-DOF model are denoted by xn, with x1 at the fixed end and x5 at the free-end. Performing the eigen analysis of the 5DOF system, we obtain its five modal frequencies: 34.499, 100.700, 158.730, 203.880, 232.520 Hz, and the corresponding damping ratios: 0.0037374, 0.010909, 0.017197, 0.022092, 0.025198. We have employed Matlab function impulse to generate discrete IRFs in this study [22]. Each IRF contains 1024 time steps with a sampling rate 500 Hz, i.e. time resolution Dt = 0.002 s. In turn, FRFs obtained from the Fourier transform of IRFs have the Nyquist frequency 250 Hz, and frequency resolution 0.4883 Hz. Displayed in Fig. 2 are 5 receptance (the ratio between displacement response and force) FRFs of the above 5-DOF system, in which Fij denotes the receptance FRF with input at coordinate xj and output at xi. While each FRF curve always has 5 local maxima at the modal frequencies of the system, some peaks might become small. For instance, F15 has its fifth peak being very small relative to its other peaks, an indication that the contribution of the fifth mode to F15 is very small, thus the modal property of mode 5 might not be easy to be identified from a noisy F15. Because F11 contains 5 clear peaks with relatively comparable magnitude and we intend to demonstrate the identification of all 5 modal frequencies and damping ratios, corrupted IRFs associated with i =1 and j =1 are simulated for the numerical studies below. Corrupted IRFs are referred to as ‘‘measured’’ IRFs in this section. They are generated by adding a Gaussian white noise to the exact (noise-free) IRFs. The level of the additive white noise is quantified by a stated percentage, defined as the ratio of the standard deviation of the white noise to that of the exact IRF. Shown in Fig. 3(a) are a segment of the exact IRF and that of 5% white noise signal. If we overlapped the exact and corrupted IRFs in one figure, we would not be able to see the difference between them because the scale of the noise is much smaller than that of the exact IRF. However, the difference

ARTICLE IN PRESS S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

1611

100

Normalized singular values

5% noise 20% noise exact 10−1

10−2

10−3 2

4

6

8 10 12 14 Index of singular values

16

18

20

Fig. 4. Normalized singular values of Hankel matrices.

between the corresponding corrupted and exact FRFs could be clearly observed when we overlapped them in one figure, as demonstrated in Fig. 3(b). Thus, the evaluation of the performance of the noise elimination is better done by comparing FRFs, instead of IRFs. 4.1. Model order determination Determining the model order of an IRF is by estimating the rank of its related Hankel matrix. For an IRF with 1024 steps, the corresponding largest (nearly square) Hankel matrix Hm  n is either a 512  513 or a 513  512 matrix, noting that n +m  1 is equal to 1024, the length of the IRF. Without losing generality, we assume mZn for a Hankel matrix and the number of columns n is used to indicate the size of the Hankel matrix. Estimating the rank of a matrix (based on a chosen singular value threshold) could be easily done if the singular values of the matrix have been ordered sequentially from the largest to the smallest. The ordered singular values associated with the Hankel matrix of size n =512 for the exact IRF, as well as for the corrupted IRFs with noise levels 5% and 20%, are plotted in Fig. 4, where each singular value has been normalized by the first (largest) singular value. For a 5-DOF system, the rank of the Hankel matrix H513  512 associated with the exact IRF must be equal to 10, which is twice the number of modes of the system. Indeed, this theoretical prediction has been verified in Fig. 4 as the singular value sk drops to zero when the index k is greater than 10. For a corrupted IRF, due to the presence of noise, sk for k 410 would not drop to zero, but approach to an asymptote. Thus the asymptotic singular value could be chosen as the noise threshold in the rank determination for the measured IRF. Since the asymptotes associated with 5% and 20% noise levels begin at k= 11 and 9, respectively, the model order of these two measured IRFs should be taken equal to 10 and 8, respectively. The model order of the measured IRF with 20% noise is 8 because the signal of the fifth (weakest) mode is very weak and has been buried in the strong noise. As a result, the fifth mode will not be identifiable from this 20%-noise IRF. 4.2. Corrupted IRF with 5% noise The model order of the 5%-noise IRF, i.e. the rank of a Hankel matrix associated with this IRF, has been determined equal to 10 (see Fig. 4). In this section, this 5%-noise IRF is chosen to demonstrate the performance of the noise removal below. 4.2.1. Noise removal An important aspect of the proposed noise elimination scheme is related to the chosen size of a Hankel matrix produced from the measured IRF. To implement any SLRA algorithm for obtaining a rank-10 Hankelmatrix requires a minimal size of the Hankel matrix with n = 11. Although a larger n could be chosen, traditional wisdom suggests that a smaller n is preferred because operating the SVD on a smaller Hankel matrix could save computational time. In addition, the problem formulations and solution algorithms in most SLRA articles have been restricted to the rank reduction of a Hankel matrix by one only [10–12]. Fig. 5 shows the evolution of noise elimination after 100, 1000, and 10000 iterations, respectively, by using Cadzow’s algorithm for the Hankel matrix with n= 11. Although noise could be quickly and effectively eliminated at the low frequency region, noticeable noise at the high frequency region remained even after 10000 iterations. This numerical observation basically agrees with De Moor’s conclusion that applying Cadzow’s algorithm would converge slowly [9].

ARTICLE IN PRESS 1612

S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

−5

−5

10

10

exact filtered

−6

−6

10 Magnitude

Magnitude

10

−7

10

−8

−7

10

−8

10

10

−9

10

exact filtered

−9

0

50

100 150 Frequency [Hz]

200

250

10

0

50

100 150 Frequency [Hz]

200

250

−5

10 Magnitude

exact filtered

−10

10

0

50

100 150 Frequency [Hz]

Phase [°]

50 0 −50 −100 −150 −200

200

250

exact filtered

0

50

100 150 Frequency [Hz]

200

250

Fig. 5. Performance of noise elimination for the Hankel matrix with n= 11: (a) after 100 iterations; (b) after 1000 iterations and (c) after 10000 iterations

−5

10 Magnitude

exact filtered

−10

10

Phase [°]

0 50 0 −50 −100 −150 −200

50

100 150 Frequency [Hz]

200

250

exact filtered

0

50

100 150 Frequency [Hz]

200

250

Fig. 6. Performance of noise elimination for the Hankel matrix with n =512 after 2 iterations.

ARTICLE IN PRESS S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

1613

0.05 n = 512 n = 100 n = 30 n = 11

0.045 0.04 0.035 α

0.03 0.025 0.02 0.015 0.01 0.005 0 10

1

10

2

3

10 10 Number of iterations

4

10

Fig. 7. a against the number of iterations for Hankel matrices with n =11, 30, 100, and 512.

Table 1 The estimated modal frequencies from the clean, filtered and measured (5% corrupted) IRFs. Mode

1 2 3 4 5

Clean

Filtered

Measured

Relative difference (%)

fc

ff

fm

(ff  fc)/fc

(fm  fc)/fc

34.499 100.700 158.730 203.890 232.490

34.499 100.690 158.750 203.720 232.750

34.469 100.610 158.610 201.690 –

0  0.010 0.013  0.083 0.112

 0.087  0.089  0.076  1.079 –

Cadzow’s algorithm alternates iteratively between two steps: (unstructured) low rank approximation (via TSVD technique) and anti-subdiagonal averaging. Although it is more computationally efficient with a small size of Hankel matrix during the TSVD calculation, the convergence rate would improve if the anti-subdiagonal averaging has been taken from a larger size of Hankel matrix. Working on the Hankel matrix with n= 512, which is the maximal size for an IRF with 1024 samples, we obtain an excellent result of the noise elimination even after just 2 iterations (see Fig. 6). A way to quantify the performance of the noise reduction is based on the L2-norm of the difference between the true and filtered ^ and vectors F and F^ are the true and filtered discrete FRFs, respectively. Introducing a FRFs, denoted JDFJ2 , where DF ¼ FF as a normalized JDFJ2

a¼

JDFJ2 JFJ2

ð20Þ

we plot a against the number of iterations for Hankel matrices with n = 11, 30, 100, and 512, respectively (see Fig. 7). Note that a is numerically equivalent to the percentage noise level defined earlier, and their identity could be proved based on Parseval’s identity that the sum of the squares of the Fourier coefficients of a function is equal to the square integral of the function [23]. Thus, the a value between the true FRF and the measured 5%-noise FRF should be close to 0.05. It takes only 2 iterations for the Hankel matrix with n= 512 to reach an asymptotic value a = 0.007. In contrast, the a value associated with n =11 after 10000 iterations is about 0.01, which is still larger (worse) than that associated with n = 512 after just 1 iteration. According to our computer code, the computational time required in a single iteration for a Hankel matrix of n= 512 is about 30 times more than that of n = 11. However, the number of iterations required to converge for n =11 is about 10000 times more than that for n = 512. Overall, a more efficient choice for the size of the Hankel matrix should be n = 512. 4.2.2. System identification using complex exponential method After carrying out the noise removal to obtain the filtered IRF, we apply Prony’s method for estimating modal parameters. Since the determined model order is 10, the corresponding Hankel matrix with n = 10, together with the righthand-side vector, shown in Eq. (9) must be formed for solving 10 polynomial coefficients bm, m =0, y, 9. Because the filtered IRF has already been processed based on a model order 10, and would fit a tenth order polynomial exactly, it makes

ARTICLE IN PRESS 1614

S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

Table 2 The estimated damping ratios from the clean, filtered and measured (5% corrupted) IRFs. Mode

1 2 3 4 5

Clean

Filtered

Measured

Relative difference (%)

xc

xf

xm

(xf  xc)/xc

(xm  xc)/xc

0.0037368 0.010909 0.017195 0.022077 0.025392

0.0037213 0.010939 0.017508 0.022694 0.022477

0.0051915 0.011947 0.018685 0.021776 –

 0.415 0.275 1.820 2.795  11.480

38.929 9.515 8.665  1.363 –

−5

10

20% noise exact

−6

10 Magnitude

−7

10

−8

10

−9

10

−10

10

0

50

100 150 Frequency [Hz]

200

250

Fig. 8. A comparison of the corrupted and clean FRFs.

no difference whether the solution formulation for Eq. (9) has been just-determined or over-determined. Substituting the obtained polynomial coefficients bm, m= 0, y, 9, together with b10 = 1, into Eq. (8), we then numerically solve the roots of the tenth order polynomial. The 10 roots must be 5 complex conjugate pairs, which then are used to determine the frequency and damping ratio of each mode. Table 1 lists the resulting 5 modal frequencies estimated from the filtered IRF, as well as those computed from the clean (noise-free) and measured (5% corrupted) IRFs with model order 10. As expected, the modal frequencies computed from the noise-free IRF are virtually identical to their theoretical counterparts. All five modal frequencies estimated from the filtered IRF are very close to those calculated from the exact IRF, with the relative difference ranging from 0% to 0.112% (see Table 1). If the measured IRF was used, we could not obtain an estimate for the fifth modal frequency, because solving for the roots of the related tenth order polynomial resulted in only 4 complex conjugate pairs. The result of the estimated damping ratios is given in Table 2. The first four damping ratios estimated from the filtered IRF agree well with those calculated from the exact IRF, all being with less than 3% relative error. The relative error of the fifth damping ratio is 11.48%, which is not bad, especially recognizing that the fifth mode contributes only a small portion to the measured IRF and could not be identified if the measured IRF was used in the identification process. In brief, the proposed procedure can estimate the damping ratios properly, and the improvement of estimating damping ratios from using filtered IRF (instead of measured IRF) is more significant than that of modal frequencies.

4.3. Corrupted IRF with 20% noise We repeated the same analysis for the corrupted IRF with a 20% noise level. A comparison between the corrupted and exact FRFs is shown in Fig. 8. The model order of the corrupted IRF, or the rank of a Hankel matrix associated with this IRF, has been determined equal to 8 (see Fig. 4). Fig. 9 is the filtered FRF based on a model order 8 in comparison with the exact FRF, where the matching of the first four modes is excellent. If the model order has been set equal to 10, the resulting FRF would have a computational mode (a small spike) at frequency near 25 Hz, instead of the true mode at frequency near 232 Hz (see Fig. 10). This incapability of matching the fifth mode could be explained by observing Fig. 8, in which the peak of the true mode near 232 Hz has been completely buried in the noise.

ARTICLE IN PRESS S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

1615

−5

10

exact filtered

−6

Magnitude

10

−7

10

−8

10

−9

10

0

50

100 150 Frequency [Hz]

200

250

Fig. 9. The performance of noise elimination based on the model order 8.

−5

10

exact filtered

−6

Magnitude

10

−7

10

−8

10

−9

10

0

50

100 150 Frequency [Hz]

200

250

Fig. 10. The performance of noise elimination based on the model order 10.

Table 3 The estimated modal frequencies from the clean, filtered and measured (20% corrupted) IRFs. Mode

1 2 3 4 5

Clean

Filtered

Measured

Relative difference (%)

fc

ff

fm

(ff  fc)/fc

(fm  fc)/fc

34.499 100.700 158.730 203.890 232.490

34.499 100.690 158.840 203.150 –

34.650 101.180 161.820 198.660 –

0  0.010 0.069  0.363 –

0.438 0.477 1.947  2.565 –

Applying Prony’s method to the filtered IRF with model order 8 results in 4 estimated modal frequencies and damping ratios (see Tables 3 and 4). They are very close to their true values, with relative errors from 0% to 0.363% for modal frequencies, and 1.109% to 8.778% for damping ratios. In contrast, the relative errors would be from 0.438% to 2.565% for modal frequencies, and 11.905% to 207.275% for damping ratios, if the measured IRF has been utilized instead.

ARTICLE IN PRESS 1616

S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

Table 4 The estimated damping ratios from the clean, filtered and measured (20% corrupted) IRFs. Mode

1 2 3 4 5

Clean

Filtered

Measured

Relative difference (%)

xc

xf

xm

(xf  xc)/xc

(xm  xc)/xc

0.0037368 0.010909 0.017195 0.022077 0.025392

0.0036717 0.01103 0.018478 0.024015 –

0.0020967 0.0088486 0.019242 0.067837 –

 1.742 1.109 7.461 8.778 –

 43.890  18.887 11.905 207.275 –

Fig. 11. A sketch of the experimental set-up.

2

2 1.5

1

1

0.5

Acceleration (g)

Acceleration (g)

free−end 1.5

0 −0.5

0.5 0 −0.5

−1

−1

−1.5

−1.5

−2 0

5

10

15 20 Time (s)

25

30

35

midpoint

−2 0

5

10

15 20 Time (s)

25

30

35

Fig. 12. Measured IRFs at locations: free-end and midpoint.

5. Experimental study The performance of the proposed method is also investigated using real measurements. A sketch of the experimental set-up is shown in Fig. 11 where the test structure is a cantilever beam. Two accelerometers, termed Sensors 1 and 2, are mounted at the locations near the top and midpoint of the beam, respectively. Exciting the beam by an impulsive load, we

ARTICLE IN PRESS S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

1617

0

10

Normalized singular values

Sensor 1 Sensor 2

−1

10

−2

10

−3

10

0

5

10 15 Index of singular values

20

25

Fig. 13. Normalized singular values of the square Hankel matrices based on measured IRFs at two locations.

1

200

10

measured filtered

measured filtered

150 100

0

10

Phase [°]

Magnitude

50

−1

10

0 −50 −100 −150

−2

−200

10

0

20

40 60 Frequency [Hz]

80

100

0

20

40 60 Frequency [Hz]

80

100

80

100

Fig. 14. Comparison of the measured and filtered FRFs of Sensor 1 (free-end sensor).

1

10

200

measured filtered 0

10

100 Phase [°]

Magnitude

measured filtered

150

−1

10

50 0 −50

−2

10

−100 −150

−3

10

−200

0

10

20

30

40 50 60 70 Frequency [Hz]

80

90 100

0

20

40

60

Frequency [Hz]

Fig. 15. Comparison of the measured and filtered FRFs of Sensor 2 (midpoint sensor).

ARTICLE IN PRESS 1618

S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

recorded the corresponding acceleration signals (see Fig. 12) for about 33 s with sampling rate 200 Hz. However, only the segment between 1 and 2.5 s with 301 sample points would be taken for the later analysis. In this section, these 301-point signals are referred to as the measured IRFs, from which the corresponding 151  151 square Hankel matrices could be produced. The normalized singular values of the square Hankel matrices are shown in Fig. 13. Because the seventh normalized singular value drops to a much smaller value for both Sensors 1 and 2, the noise threshold is set to have the model order equal to 6. The measured FRF and the corresponding filtered FRF (with model order 6) associated with Sensor 1 are shown in Fig. 14. The filtered FRF has 3 sharp peaks within the frequency range 0–100 Hz, and the visible level of noise at the measured FRF has been successfully removed. Similarly, Fig. 15 shows the measured and filtered FRFs associated with Sensor 2, and the noise removal from the measured FRF appears to be very good also. Observing the FRFs in Figs. 14 and 15, we noticed that the magnitude of the first peak (the first mode) is much smaller than that of the second or third peak. This indicates that the first mode is a much weaker contributor to its FRF, and so estimating the modal parameters of the first mode, from either measured or filtered IRFs, would be influenced more by the noise. As the second mode is the strongest one, the modal parameter estimation for the second mode is expected to be more accurate, and less affected by the noise. The first three modal frequencies estimated from the measured and filtered signals are listed in Table 5. Their counterparts computed from a finite element (FE) model are also included for comparison. The three modal frequencies obtained from the experimental data basically agree with those of the FE model, recognizing that discrepancy always exists between the true

Table 5 The first three modal frequencies estimated from the measured and filtered signals. Mode

FEM

1 2 3

Sensor 1

5.8495 36.650 102.580

Sensor 2

Relative difference (%)

fm1

ff1

fm2

ff2

(fm1  fm2)/fm2

(ff1  ff2)/ff2

5.2934 34.691 97.112

5.5146 34.691 97.162

5.3470 34.691 97.119

5.5157 34.691 97.162

 1.002 0  0.007

 0.020 0 0

Table 6 The first three damping ratios estimated from the measured and filtered signals. Mode

Sensor 1

Relative difference (%)

xm1

xf1

xm2

xf2

(xm1  xm2)/xm2

(xf1  xf2)/xf2

 0.0038926 0.0016358 0.0023735

0.0018143 0.0018111 0.0012603

0.0040007 0.0017985 0.0022811

0.0016805 0.0018115 0.0012556

 197.298  9.046 4.051

7.962  0.022 0.374

35

35

30

30

25

25 Model order

Model order

1 2 3

Sensor 2

20 15

20 15

10

10

5

5 0

20

40 60 Frequency (Hz)

80

100

0

10

20

30

40 50 60 70 Frequency (Hz)

80

90 100

Fig. 16. Stability diagrams corresponding to measured and filtered signals of sensor 1. The used symbols are: ‘‘o’’, not stable; ‘‘*’’, stable frequency and damping.

ARTICLE IN PRESS S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

1619

physical structure and its FE model. Modal frequencies are global parameters and theoretically they should remain the same when they are estimated from different IRFs. Due to a lack of the true values for the modal frequencies of the test structure, examining the agreement among the modal frequencies estimated from different signals seems to be the most logical way to evaluate the estimation performance. As exhibited in Table 5, the values for the second modal frequency (associated with the strongest modal component) estimated from the measured and filtered IRFs of Sensors 1 and 2 are all equal to 34.691 Hz. This estimated value has not been influenced by the measurement noise, and is believed to be accurate because of the consistency among different signals. The modal frequencies estimated from the two filtered IRFs associated with Sensors 1 and 2 are in excellent agreement—with the relative difference 0.02%, 0%, and 0% for the first three modal frequencies, respectively. Shown in Table 6 are the values for the damping ratios of the first three modes estimated from the measured and filtered IRFs. The damping ratios estimated from the measured signals are poor; in particular at Sensor 1 the estimation for the first mode is a negative quantity which is physically impossible. The damping ratios estimated from the filtered signals seem reasonable as the relative differences between Sensors 1 and 2 of the first three damping ratios are 7.96%, 0.02%, and 0.37%, respectively. As explained previously, the estimated first modal parameters (associated with the weakest modal component) would be influenced most by the measurement noise, and improved most by the noise removal. In performing experimental modal analysis, over-determined models have been often used, and stability diagrams have been routinely plotted to eliminate spurious numerical poles (modes). When producing a stability diagram, the poles corresponding to a certain model order are compared to the poles of a one order-lower model [24]. If the modal frequency and the damping ratio differences are within preset limits, the pole is labeled as a stable one. Adopting the preset limits with 1% difference for both frequencies and damping ratios, we plotted the stability diagrams based on the measured and filtered signals of Sensor 1 (see Fig. 16). Although the stability diagram associated with the measured signal has consistent frequency estimates near the three system frequencies extending from higher order to lower order models, many poles have been labeled unstable because they did not meet the preset 1% stability standard. In contrast, the stability diagram associated with the filtered signal has consistent and stable poles all the way from models with higher order to the model with order 3. 6. Concluding remarks Due to the noise and the uncertainty of the number of modes (or model order) in a measured FRF/IRF, accurately estimating modal parameters has been a challenging task. We proposed a modal identification scheme that is fundamentally different from the established schemes, which purposely increased the model order of their identification models to absorb noise. Our proposed method included three steps: (1) determining the rank of a Hankel matrix associated with the measured IRF based on an observed ‘‘noise threshold,’’ where the rank is theoretically twice the number of modes contributing to the measured IRF, (2) implementing Cadzow’s algorithm for the structured low rank approximation (SLRA) on the Hankel matrix to achieve the noise removal from the measured IRF to obtain a filtered IRF, and (3) using the complex exponential method (Prony’s method) to estimate the modal parameters from the filtered IRF with a known model order. In implementing the iterative Cadzow’s algorithm, selecting a proper size of the Hankel matrix is important. We found that choosing a nearly square Hankel matrix would be more computationally efficient. Both synthesized and experimental data were included in the numerical studies. For a synthesized 5-DOF mass-springdashpot system, measured IRFs with two different levels of noise were investigated. When the noise level was low, we could accurately estimate the modal parameters (frequencies and damping ratios) of all 5 modes by using the proposed scheme. When the noise level was high, the modal parameters of the first 4 modes could be identified nicely, but those of the fifth mode was unidentifiable because the signal associated with the fifth mode was too weak relative to the noise. For both low and high noise cases, we observed significant improvement in modal parameter estimation when the filtered, rather than measured, IRF was used. Our experimental data were measured from two accelerometers mounted near the free-end and midpoint, respectively, of a cantilever beam. Although the true values of the modal parameters of the test structure were not known, we could still judge the performance of the proposed scheme based on examining the agreement among the modal parameters estimated from IRFs measured from those two accelerometers, because modal frequencies and damping ratios are global parameters and they should be identical when estimated from different IRFs. Carrying out the proposed method, we found that the modal parameters estimated from the filtered IRFs of the two accelerometers were in excellent agreement.

Acknowledgements The authors acknowledge the financial support by the China National Science Key Fund under the grant number 50739004, and by the 863-project of China (Program number 2008AA092701). References [1] D.J. Ewins, Modal Testing: Theory, Practice and Applications, 2nd edition, Research Studies Press, Baldock, Hertfordshire, England, 2000. [2] N. Maia, J. Silva, J. He, N. Lieven, R.-M. Lin, G. Skingle, W. To, A. Urgueira, Theoretical and Experimental Modal Analysis, Research Studies Press, Taunton, Somerset, England, 1997.

ARTICLE IN PRESS 1620

S.-L. James Hu et al. / Mechanical Systems and Signal Processing 24 (2010) 1605–1620

[3] R.J. Allemang, D.L. Brown, A unified matrix polynomial approach to modal identification, Journal of Sound and Vibration 211 (3) (1998) 301–322. [4] S. Braun, Y.M. Ram, Determination of structural modes via the Prony method: system order and noise induced poles, Journal of the Acoustical Society of America 81 (1987) 1447–1459. [5] S. Braun, Y.M. Ram, Time and frequency identification methods in over-determined systems, Mechanical Systems and Signal Processing 1 (1987) 245–257. [6] S. Braun, Y.M. Ram, Structural parameter identification in the frequency domain: the use of over-determined systems, ASME Transactions, Journal of Dynamic Systems, Measurement and Control 109 (1987) 120–123. [7] K.Y. Sanliturk, O. Cakar, Noise elimination from measured frequency response functions, Mechanical Systems and Signal Processing 19 (2005) 615–631. [8] M.T. Chu, R.E. Funderlic, R.J. Plemmons, Structured low rank approximation, Linear Algebra and its Applications 366 (2003) 157–172. [9] B. De Moor, Total least squares for affinely structured matrices and the noisy realization problem, IEEE Transactions on Signal Processing 42 (11) (1994) 3104–3113. [10] I. Markovsky, S.V. Huffel, Overview of total least squares methods, Signal Processing 87 (10) (2007) 2283–2302. [11] B. De Moor, Structured total least squares and L2 approximation problems, Linear Algebra and its Applications 188–189 (1993) 163–207. [12] J. Rosen, H. Park, J. Glick, Total least norm formulation and solution of structured problems, SIAM Journal on Matrix Analysis and Applications 17 (1996) 110–126. [13] H. Park, L. Zhang, J.B. Rosen, Low rank approximation of a Hankel matrix by structured total least norm, Bit Numerical Mathematics 39 (1999) 757–779. [14] M. Schuermans, P. Lemmerling, S.V. Huffel, Block-row Hankel weighted low rank approximation, Numerical Linear Algebra with Applications 13 (2006) 293–302. [15] I. Markovsky, Structured low-rank approximation and its applications, Automatica 44 (2008) 891–909. [16] J.A. Cadzow, Signal enhancement—a composite property mapping algorithm, IEEE Transactions on Acoustics, Speech, and Signal Processing 36 (1) (1988) 49–62. [17] A.V. Oppenheim, R.W. Schafer, Discrete-Time Signal Processing, Prentice-Hall Press, Upper Saddle River, New Jersey, 1999. [18] L. Weiss, R.N. McDonough, Prony’s method, Z-transforms, and pade approximation, SIAM Review 5 (2) (1963) 145–149. [19] G. Golub, C.V. Loan, Matrix Computations, 3rd edition, Johns Hopkins University Press, Baltimore, MD, 1996. [20] G. Eckart, G. Young, The approximation of one matrix by another of lower rank, Psychometrika 1 (1936) 211–218. [21] D. Tufts, A. Shah, Estimation of a signal waveform from noisy data using low-rank approximation to a data matrix, IEEE Transactions on Signal Processing 41 (4) (1993) 1716–1721. [22] The MathWorks, Control system toolbox user’s guide, The MathWorks Inc., Natic, MA, 2004. [23] L.W. Johnson, R.D. Riess, Numerical Analysis, 2nd edition, Addison-Wesley Press, Reading, Massachusetts, 1982. [24] B. Peeters, System identification and damage detection in civil engineering, Ph.D. thesis, Katholieke Universiteit Leuven, 2000.