18.4 Curve Fitting of Transfer Functions

Back

Parameter estimation in vibrating systems can be interpreted as a technique of experimental modeling.

This process requires experimental data in a suitable form, preferably excitation – response data, and is

often represented as a set of transfer functions in the frequency domain. Parameter estimation using

18-10 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

measured response data is termed model identification or, simply, identification in the literature on

systems and control. We shall present a parameter estimation procedure that involves frequency transfer

functions, which is particularly useful in EMA.

18.4.1 Problem Identification

Transfer functions that are computed from measured time histories using digital Fourier analysis

(e.g., FFT) cannot be directly used in the modal analysis computations. The data must be available

as analytical transfer functions. Therefore, it is important to represent the computed transfer

functions with suitable analytical expressions. This is done, in practice, either by curve fitting

Box 18.1

MAIN STEPS OF EXPERIMENTAL

MODAL ANALYSIS

1. Measure an admissible set of excitation ðuÞ and response ðyÞ signals. (Cover all DoF; one

response measurement should be for the excitation location.)

2. Group the signals, assign windows, and filter the signals.

3. Compute transfer functions using FFT and the spectral density method:

Guy ¼ Fuy =Fuu

4. Compute ordinary coherence functions

g2

uy ¼

lFuy l2

FuuFyy

and choose the accurate transfer functions on this basis (guy close to 1 ) accept; guy close

to 0 ) reject).

5. Curve fit n admissible transfer functions to expressions:

Gik ¼

Xn

r¼1

􀀑 ðcickÞr

v2r

2 v2 þ 2jzrvrv

􀀜

Hence, extract:

Residues ðcickÞr ) mode shapes vectors cr which are M-normal;

Natural frequencies (undamped), vr ;

Modal damping ratios (viscous), zr :

6. Form the modal matrix C ¼ ½c1; c2; …; cn􀀉;

Compute C21:

7. Modal mass matrix M􀀊 ¼ I;

Modal stiffness matrix K􀀊 ¼ diag

􀀑

v21

; v22

; …; v2

n

􀀜

;

Modal damping matrix C􀀊 ¼ diag½2z1v1; 2z2v2;…; 2znvn􀀉;

8. Compute the system model:

Mass matrix M ¼ ðC21ÞC21;

Stiffness matrix K ¼ ðC21ÞTK􀀊C21;

Damping matrix C ¼ ðC21Þ C􀀊C21.

Experimental Modal Analysis 18-11

© 2005 by Taylor & Francis Group, LLC

a suitable transfer function model into the computed data or by simplified methods such as “peak

picking.” Accordingly, this conversion of data is an experimental modeling technique.

Identification of transfer function models from measured data is an essential step in EMA. Apart from

that, it has other important advantages. In particular, analytical transfer function plots clearly identify

system resonances and generate numerical values for the corresponding parameters (resonant

frequencies, damping, phase angles, and magnitudes) in a convenient manner. This form represents a

significant improvement over the crude transfer function plots, which are normally far less presentable

and rather difficult to interpret.

18.4.2 Single- and Multi-Degree-of-Freedom Techniques

Several single-DoF techniques exist for extracting analytical parameters from experimental transfer

functions. In particular, the methods of curve fitting (circle fitting) and peak picking are considered here.

In a single-DoF method, only one resonance is considered at a time.

Single-DoF curve fitting, or more correctly, single-resonance curve fitting is the term used to denote any

curve fitting procedure that fits a quadratic (second-order) transfer function into each resonance in the

measured transfer function, one at a time. In the case of closely spaced modes (or closely spaced

resonances), the associated error can be very large. The accuracy is improved if expressions of a higher

order than quadratic are used for this purpose, but unacceptable errors can still exist. In peak picking,

each resonance of experimental transfer function data is examined individually; the resonant frequency

and the damping constant corresponding to that resonance are determined by comparing with an

analytical single-DoF transfer function.

In multi-DoF curve fitting, or more appropriately, multiresonance curve fitting, all resonances (or

modes) of importance are considered simultaneously and fitted into an analytical transfer function of

suitable order. This method is generally more accurate but computationally more demanding than the

single-resonance method. In choosing between the single-resonance and multiresonance methods, the

required accuracy should be weighted against the cost and speed of computation.

18.4.3 Single-Degree-of-Freedom Parameter Extraction in

the Frequency Domain

The theory of curve fitting by a circle (i.e., circle fitting) for each resonance of an experimentally

determined transfer function is presented first. Next, the peak picking method will be described.

18.4.3.1 Circle-Fit Method

It can be shown that the mobility transfer function (velocity/force) of a single-DoF system with linear

viscous damping, when plotted on the Nyquist plane of real axis and imaginary axis for the frequency

transfer function, is a circle. Similarly, it can be shown that the receptance or dynamic flexibility or

compliance transfer function (displacement/force) of a single-DoF system with hysteretic damping, when

plotted on the Nyquist plane, is also a circle. Note that, for hysteretic damping, the damping constant (in

the time domain) is not actually a constant but is inversely proportional to the frequency of motion.

However, in this case, in the frequency domain, the damping term will be independent of frequency. The

fact that such circle representations are possible for transfer functions of a single-DoF system may be used

in fitting a circle to a transfer function that is computed from experimental data. This will lead to

determining the analytical parameters for the transfer function. This approach is illustrated now through

analytical development.

Case of Viscous Damping

Consider a single-DoF system with linear, viscous damping, as given by

my€ þ cy_ þ ky ¼ f ðtÞ ð18:39Þ

18-12 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

where m; k; and c are the mass, stiffness, and damping constants of the system elements, respectively, f ðtÞ

is the excitation force, and y is the displacement response. Equation 18.39 may be expressed in the

standard form:

y€ þ 2zvny_ þv2

ny ¼

1

m

f ðtÞ ð18:40Þ

Receptance

Y ðsÞ

FðsÞ ¼

1

m

􀀑

s2 þ 2zvns þ v2

n

􀀜 with s ¼ jv ð18:41Þ

Mobility

V ðsÞ

FðsÞ ¼

sY ðsÞ

FðsÞ ¼

s

m

􀀑

s2 þ 2zvns þ v2

n

􀀜 with s ¼ jv ð18:42Þ

Consider the mobility (velocity/force) transfer function given by

GðsÞ ¼

s

s2 þ 2zvns þ v2

n ð18:43Þ

where the constant parameter, m, in Equation 18.42 has been omitted, without loss of generality. In the

frequency domain ðs ¼ jvÞ, we have

GðjvÞ ¼

jv

v2

n 2 v2 þ 2jzvnv ð18:44Þ

Multiply the numerator and the denominator of GðjvÞ in Equation 18.44 by the complex conjugate of the

denominator (i.e., v2

n 2 v2 2 2jzvnv). Then, the denominator is converted to the square of its original

magnitude, as given by

D ¼ ðv2

n 2 v2Þ2 þ ð2zvnvÞ2 ð18:45Þ

and the frequency transfer function (Equation 18.44) is converted into the form

GðjvÞ ¼

jv

D

􀀑

v2

n 2 v2 2 2jzvnv

􀀜

ð18:46Þ

Gð jvÞ ¼ Re þ j Im where Re ¼

2zvnv2

D

and Im ¼

v

D ðv2

n 2 v2Þ ð18:47Þ

Now, we can write

Re 2

1

4zvn ¼

8z2v2

nv2 2 4z2v2

nv2 2 ðv2

n 2 v2Þ2

4zvnD ¼

4z2v2

nv2 2 ðv2

n 2 v2Þ2

4zvnD

Hence, in view of Equation 18.47 we have

Re 2

1

4zvn

􀀒 􀀓2

þIm2 ¼

4z2v2

nv2 2 ðv2

n 2 v2Þ2

4zvnD

" #2

þ

v2ðv2

n 2 v2Þ2

D2

¼

16z 4v4

nv4 2 8z2v2

nv2ðv2

n 2 v2Þ2 þ ðv2

n 2 v2Þ4 þ 16z2v2

nv2ðv2

n 2 v2Þ2

16z2v2

nD2

¼

􀀑

4z2v2

nv2 þ ðv2

n 2 v2Þ2􀀜2

16z2v2

nD2 ¼

D2

16z2v2

nD2 ¼

1

16z2v2

n ¼ R2

It follows that the transfer function, GðjvÞ; represents a circle in the real – imaginary plane, with the

following properties:

Circle radius R ¼

1

4zvn ð18:48Þ

Experimental Modal Analysis 18-13

© 2005 by Taylor & Francis Group, LLC

Circle center ¼

1

4zvn

; 0

􀀏 􀀐

ð18:49Þ

Now, we may reintroduce the constant parameter,

m; back into the transfer function, as in Equation

18.42. Then, we have

Circle radius R ¼

1

4zvnm

;

Center ¼

1

4zvnm

; 0

􀀏 􀀐 ð18:50Þ

A sketch of this circle is shown in Figure 18.3(a). As

mentioned before, the plane formed by the real and

imaginary parts of GðjvÞ as the Cartesian x and y

axes, respectively, is the Nyquist plane. The plot of

GðjvÞon this plane is the Nyquist diagram. It follows

that the Nyquist diagram of the mobility function

(Equation 18.42 or Equation 18.44) is a circle.

Case of Hysteretic Damping

Consider a single-DoF system with hysteretic

damping. The equation motion is given by

my€ þ

h

v

y_ þ ky ¼ f ðtÞ

with f ðtÞ ¼ f0 sin vt

ð18:51Þ

Note the frequency dependent damping constant, with the hysteretic damping parameter, h; in the time

domain. The receptance function, Gð jvÞ; is given by

Gð jvÞ ¼

1

ðk 2 mv2Þ þ jh

Note that the damping term, jh; is independent of frequency in the frequency domain for this case. As for

the case of viscous damping, we can easily show that the Nyquist plot of this transfer function is a circle with

Radius ¼

1

2h

and Center ¼ 0; 2

1

2h

􀀏 􀀐

ð18:52Þ

A sketch of the resulting circle is shown in Figure 18.3(b).

In general, for a multi-DoF viscous-damped system we have the “mobility” function

Gikð jvÞ ¼ jv

Xn

r¼1

􀀑 ðcickÞr

v2r

2 v2 þ 2jzrvrv

􀀜 ð18:53Þ

If the resonances are not closely spaced we can assume that, near each resonance ðrÞ

Gik ¼ constant offset ðcomplexÞ þ single DoF mobility function

¼ constant offset þ

jvðcickÞr

v2r

2 v2 þ 2jzrvrv

􀀒 􀀓

ð18:54Þ

We can curve fit each resonance r to a circle this way and thereby extract the ðcickÞr value (the residue)

from the radius of the circle fit.

Note: This method cannot be used if the resonances are closely spaced and consequently if significant

modal interactions are present.

1

4zwnm 2zwnm

1 Re

Im

0

G( jw) Plane

G( jw) Plane

Re

Im

2h

− 1

0

(a)

(b)

FIGURE 18.3 (a) Circle fit of a mobility function with

viscous damping; (b) circle fit of a receptance function

with hysteretic damping.

18-14 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

18.4.3.2 Peak-Picking Method

The peak-picking method is also a single-DoF method in view of the fact that each resonance of an

experimentally determined transfer function is considered separately. The approach is to compare the

resonance region with an analytical transfer function of a damped single-DoF system. One of three types of

transfer functions, receptance, mobility, and accelerance as listed in Table 18.1, can be used for this

purpose. Note that, when the level of damping is small, it can be assumed (approximately) that the

resonance is at the undamped natural frequency vn ¼

ffiffiffiffiffi

k=m p : Substituting this value for v in each of

the frequency transfer functions, we can determine the transfer function value at resonance, denoted by

Gpeak ðjvÞ: It is noted from Table 18.1 that this function value in general depends on the damping constant

and the natural frequency. Since vn is known directly from the peak location of the transfer function, it is

possible to compute c (or the damping ratio, z) by first determining the corresponding peak magnitude.

Specifically, from Table 18.1, it is clear that we should pick the imaginary part of the frequency transfer

function for receptance or accelerance data and the real part of the transfer function for mobility data.

Then, we pick the peak value of the chosen part of the transfer function and the frequency at the peak.

Table 18.2 gives normalized expressions for the three frequency transfer functions, receptance ¼

displacement/force; mobility ¼ velocity/force; accelerance ¼ acceleration/force, in the frequency

domain, in the case of a single-DoF mechanical system with (1) viscous damping and (2)

hysteretic damping. It may be verified that their Nyquist diagrams are circles (either exactly or

approximately), thereby enabling one to use the circle-fit method.

Note that r ¼ v=vn; where v is the excitation frequency and vn is the undamped natural frequency;

z is the damping ratio in the case of viscous damping; and d ¼ h=k where h is the hysteretic damping

parameter and k is the system stiffness.

Peak picking is good in cases where modes are well separated and lightly damped. It does not

work when the system is highly damped (or overdamped) or when the damping is zero (infinite

peak). It is a quick approach that is appropriate for initial evaluations and trouble shooting.

TABLE 18.1 Some Frequency Transfer Functions Used in Peak Picking

Single-DoF system my€ þ cy_ þ ky ¼ f ðtÞ

Receptance (dynamic flexibility, compliance)

Gr ðjvÞ ¼

Displacement y

Force f

1

ms2 þ cs þ k

with s ¼ jv

Mobility Gm ðjvÞ ¼

Velocity v

Force f

s

ms2 þ cs þ k

with s ¼ jv

Accelerance Ga ðjvÞ ¼

Acceleration a

Force f

s2

ms2 þ cs þ k

with s ¼ jv

Resonant peaks Gpeak ðjvÞ (occur approximately

at v ¼ vn for light damping)

Gr

peak ¼ 2

j

cvn

Gm

peak ¼

1

c

Ga

peak ¼

jvn

c

TABLE 18.2 Normalized Frequency Response Functions for Single-DoF Curve Fitting

Frequency Response Function With Viscous Damping With Hysteretic Damping

Receptance

1

1 2 r2 þ 2jzr

1

1 2 r2 þ jd

Mobility

jr

1 2 r2 þ 2jzr

jr

1 2 r2 þ jd

Accelerance

2r2

1 2 r2 þ 2jzr

2r2

1 2 r2 þ jd

Experimental Modal Analysis 18-15

© 2005 by Taylor & Francis Group, LLC

18.4.4 Multi-Degree of Freedom Curve Fitting

We shall now discuss a general multiresonance curve fitting method. The corresponding single-resonance

method should also be clear from this general procedure. Note that many different versions of problem

formulation and algorithm development are possible for least squares curve fitting, but the results should

be essentially the same. The method presented here is a frequency-domain method as we are dealing with

experimentally determined frequency transfer functions. In a comparable time-domain method,

a suitable analytical expression of the complex exponential form is fitted into the experimental impulse

response function obtained by the inverse Fourier transformation of measured transfer function.

That method acquires additional error due to truncation (leakage) and finite sampling rate (aliasing)

during the inverse FFT.

18.4.4.1 Formulation of the Method

The objective of the present multiresonance (multi-DoF) curve fitting procedure is to fit the computed

(measured) transfer function data into an analytical expression of the form

GðsÞ ¼

􀀑

b0 þ b1s þ · · · þ bmsm􀀜

􀀑

a0 þ a1s þ · · · þ ap21sp21 þ sp

􀀜 for m # p ð18:55Þ

The data for curve fitting are the N complex transfer function values ½G1; G2; …; GN 􀀉 computed at discrete

frequencies ½v1; v2; …; vN 􀀉: Typically, if 1024 samples of time history were used in the FFT

computations to determine the transfer function, we would have 512 valid spectral lines of transfer

function data. However, near the high-frequency end, these data values become excessively distorted due

to the aliasing error; only a part of the 512 spectral lines will be usable, typically the first 400 lines. In that

case, we have N ¼ 400: This value can be doubled by doubling the FFT block size (to 2048 words in the

buffer), thereby doubling the record length or the sampling rate. It is acceptable to leave out part of the

computed transfer function, not for poor accuracy but because that part falls outside the frequency band

of interest in that particular modal analysis problem. A less wasteful practice is to pick the sampling rate

of the measured time-history data to reflect the highest frequency of interest in the modal analysis.

The (complex) error in the estimated value at each frequency point (spectral line), vi; is given by

e~i ¼ GðviÞ 2 Gi ¼

􀀑

b0 þ b1si þ · · · þ bmsmi 􀀜

􀀑

a0 þ a1si þ · · · þ ap21s p21

i þ sp

􀀜 2 Gi ð18:56Þ

with s ¼ jv:

The characteristic equation of the dynamic model is given by

DðsÞ ¼ a0 þ a1s þ · · · þ ap21sp21 þ sp ¼ 0 ð18:57Þ

Its roots are the eigenvalues of the system. For damped systems, they occur in complex conjugates with

negative real parts (note: p ¼ 2 £ number of DoFs, in typical cases). For systems with rigid-body modes,

zero eigenvalues will also be present. However, since there is some damping in the system and since the

lowest frequency that is tested and analyzed is normally greater than zero even for systems with rigidbody

modes, we have

DðjvÞ – 0 ð18:58Þ

in the frequency range of interest. Hence, the estimation error given by Equation 18.56 can be expressed as

ei ¼

􀀑

b0 þ b1si þ · · · þ bmsmi

􀀜

2 Gi

􀀑

a0 þ a1si þ · · · þ ap21sp21

i sp

i

􀀜

ð18:59Þ

with s ¼ jv:

The quadratic error function is given by the sum of the squares of magnitude error for all discrete

frequency points used in modal analysis; thus

J ¼

XN

i¼1

leil2 ¼

XN

i¼1

ep

i ei ð18:60Þ

18-16 Vibration and Shock Handbook

© 2005 by Taylor & Francis Group, LLC

Note that [ ]p denotes the complex conjugate. Complex conjugation is achieved by simply replacing ðjvÞ

with ð2jvÞ in Equation 18.59. It follows that Equation 18.60 can be written as

J ¼

XN

i¼1

h

b0 þ b1ð2jviÞ þ · · · þ bmð2jviÞm 2 Gp

i {a0 þ a1ð2jviÞ þ · · · þ ap21ð2jviÞp21

D

þð2jviÞp}

ih

b0 þ b1ðjviÞ þ · · · þ bmðjviÞm 2 Gi{a0 þ a1ðjviÞ þ · · · þ ap21ðjviÞp21 þ ðjviÞp}

i E

ð18:61Þ

The basis of the least squares curve fitting method of parameter estimation is to pick the transfer function

parameters bi; i ¼ 0; 1; …; m and ai; i ¼ 0; …; p 2 1; such that the quadratic error function, J; is a

minimum. Analytically, this requires

›J

dbk ¼ 0 k ¼ 0; 1; …; m ð18:62Þ

›J

dak ¼ 0 k ¼ 0; 1; …; p 2 1 ð18:63Þ

Note that Equation 18.62 and Equation 18.63 correspond to m þ p þ 1 linear equations in the m þ p þ 1

unknowns bi; i ¼ 0; 1; …; m and ai; i ¼ 0; 1; …; p 2 1: A well-defined solution exists to this set of

nonhomogeneous equations provided that the equations are linearly independent, which is guaranteed if

the determinant of the coefficients of the unknown parameters does not vanish. It is a good practice to

check for linear independence of the set of m þ p þ 1 equations using this determinant condition prior

to performing further computations to solve the equations. The solution approach itself is primarily

computational in nature and is not presented here. Figure 18.4 shows a result of multi-DoF curve fitting

on an experimental frequency transfer function, as collected from a civil engineering structure. Note the

close match of the magnitude but not the phase angle. This analysis resulted in the resonant frequency

and damping ratio values that are given in Table 18.3.

18.4.5 A Comment on Static Modes and Rigid-Body Modes

Some test systems may possess static modes, and rigid-body modes under rare circumstances. Static

modes arise in analytical models if we fail to assign an inertia (mass) element for each DoF. Rigid-body

modes arise in analytical models if proper restraints are not provided for the inertia elements. In practice,

FIGURE 18.4 An example of multi-DoF curve fitting on experimental data.

Experimental Modal Analysis 18-17

© 2005 by Taylor & Francis Group, LLC

however, static modes arise if a coordinate is assigned to a DoF that actually does not exist, or if some

parts of the physical system are relatively light with stiff restraints (i.e., very high natural frequencies), and

rigid-body modes arise in the presence of relatively heavy components restrained by very flexible

elements (i.e., very low natural frequencies). Note that the assumed transfer function (Equation 18.55)

allows for both these extremes. Specifically, if static modes are present it is necessary that the transfer

function can be expressed as a sum of a constant term (static mode) and an ordinary transfer function

(without a static mode). Hence, it will approach a nonzero constant value as the frequency, v; increases.

This requires m ¼ p: If rigid-body modes are present, the characteristic polynomial DðsÞ of the model

should have a factor s2: This corresponds to a0 ¼ a1 ¼ 0:

18.4.6 Residue Extraction

The estimated transfer functions as given by Equation 18.55 is in the form of the ratio of two

polynomials; the rational fraction form. This has to be converted into the partial fraction form given by

Equation 18.17 in order to extract the residues ðcickÞr that are needed for determining the mode shapes.

For this, the natural frequencies, vr ; and the modal damping ratios, zr ; should be computed first. These

are given by the roots of the characteristic Equation 18.57 as

lr ; lp

r ¼ zrvr ^ j

ffiffiffiffiffiffiffiffi

1 2 z2r

q

vr ; r ¼ 1; 2; …; n ð18:64Þ

Once these eigenvalues are known, by solving Equation 18.57 using the estimated values for

a0; a1; …; ap21; it is a straightforward task to compute the quadratic factors

Dr ðsÞ ¼ s2 þ 2zrvr s þ v2r

r ¼ 1; 2; …; n ð18:65Þ

Note that, from Equation 18.17

ðcicrÞ ¼ ½GikðsÞDr ðsÞ􀀉s¼lr ð18:66Þ

assuming distinct eigenvalues. This is true because, when the partial fraction form is multiplied by

Dr ðsÞ; it will cancel out with the denominator of the partial fraction corresponding to the rth mode,

leaving its residue. Then, when s is set equal to Dr ; all the remaining partial fraction terms will vanish

due to the fact that Dr ðlrÞ ¼ 0; provided that the eigenvalues are distinct. Since GikðsÞ are known

from the estimated transfer functions, the residues can be computed using Equation 18.66. Finally,

the mode shapes are determined using the procedure outlined earlier. Some curve fitting approaches

are summarized in Box 18.2.