- Email: [email protected]
Accurate Estimation of Multivariable Frequency Response Functions
Copyright © 1996 IFAt' 13Lh Triennial World Cong-ress, Sun Fran.:iscu. USA
3a-J3 6
ACCUR.ATE ESTIMATION OF MULTIVARIABLE FREQUENCY RESPONSE FUNCTIONS 1 P . G uillaunle, R. Pinteloll and J. Schoukens·
• D"pMimenl EL EC, Vcije Universi!'eil Brussel ( VUB), 8-1050 Brussei, Belgium e-lIlail: [email protected] - fa x: +J2-2-6:!9 2850
Ahstract. Th€' non parametric estimation of Jnultivariable fre.quency response functions is studied within an "errors-in-variables" 8tochastk framework. It is shown that , eve.n when t.he !;e(".ond-ord ~r momentl:! of the dist urbances are unknow n, it is possible by Ilsing periodic t:>xcitations to obtain consistf'.nt: and asymptotically effideut point es t,imat.es of t.h(' multivariable frequency re~p OJlse functions. This is not tht~ ca...:;e with ra.ndorlL oE'. x~~it. a ti ons. Besides , attention is paid to thE': desi gn of c)ptima i mu lti-input exdt.atioHs U1ax.imizing the Fisher information ma.tri x. Eventually, the th e, )retical results arc illust.rated by means of real measurements in the field of vihratio n analy sis.
K eywords. Frequency response, multi variable systems, nonparametric ident.ificat io n , maxim um likelihood, estimation theory.
I INTROD UCTION The nonparamet,rk estiruat.io ll or freq uency response func.t.iolls (or im pulse respo nse fUHd io ns ) is of primary importance ill many scientific inv e~Ligations . This is the reason why mu ch atte.nt.io n ha...:; l.)f~e n paid t o nonpa ra.Il1etric esti mators in th e litt>rahtrf' (Bendat and Pie:rso l, 1880; llrillinger, 1fill I ; All work was s up port~ b:), t he Belgia n National Fun d fo r Scient ific: R esearch, th e F lemish government (GOA-IMMI) a.nd t he Re.lgian government 11..0.; ;~ part of the Belgian programme 011 111teruniversity Pules of Attraction ( I UA P-~O) initiated by 'he Bdgian State. P"illlt! j\ 1iui ster's Office, S(;iellce Policy Progr.t.l1l ming.
The first es timat.ors which were developed consid ereti on ly the errors on t.he o utput measurem ents. In man y pra.ctica.1 s it.uations ~ however. the disturbancE':s on th e. input. mea....·m rem cnts can be as large as the ones on the outputs. To cope wit.h t he presence of input as well ~ o u tput. disturba.n ce::; "errors-in-vari ables" estim ators were develop ed. A possible approach is the instrumeut nl vnr'iables (IV) method (see for instance Soderstrom and Sto;"" ( 1989, Chapters 8 and 10) and Cooper and Emmett (1995 )). The IV met. J,od usually requi res an additional data se(luence (i.e. the im;lrumental variables) which can be- considered as a d ra.wback . A positive result is t1w fact that, unde.r some mild condit.ions, the IV est imat.e G m be prove.n to be stro ngly (;onsistent. Another possibili ty is to lU:;e the Hv f'!Stimator which rely on t he t.ot.a.l leas1. sq llar f'~o;; met hod. T he Ht) est.imates arc cOllsistent when the input/outpul (I/O) errors are all equall y large all J uncorrelated . If this is not the case, consistent estimates ('an st.ill be obtain~d by introduci ng a. scaling mat.ri x resulting in a gene ralized tot.alleast squares
4351
pro blem (Van Hu ffel and Vandewalle, 199 1). To derive t hls sealing matri x, howev~r . the "true" covariance m atr ix of th(" input. a nd output. errors s hould be a priori kn own. In th e presen t pa.per t.he non parame tric estimatio n of frequenc.y rt'ponse function wi ll be revisited ta king into accou nt th e specifk aspects of periodic excit aLions , i.e.. the availa.bility of replicat.ed obsf'rvations of the "true" I/ O da.ta. This wi ll result. in an a lmost t ri vial solution; a trivial solu tions whidl is , however , more perrOfmant. than t he. com mo nl y used estimators. Morp.-over, consistent est imates are obtai n even if t he "t rue" cava fiance matrix of t he I/O f:'rrors is unknown . Furt her , it will be s hown that. t he pro poS<'d estimator is a.'.iy rnpto t ically efficient. Confidence intervals for the estimated fr('qu ency res ponse functi ons will be given. Next, a.ttent ion will be s penl 1.0 the design of o ptimal mul t i-inpu t cxcitations.
u, £1/. 1
Eu . 2
U,
G=
UO , 2
[V"
G,~ ]
V,
(J~",
G~ l
E y •p
E .. , ",
y.
u~
Uo""
Fig.
I . Multivariable frequen cy-domain vari ables" set up.
"errors-in-
C pxM J relation (I) can h e expressed fo r al1 t he data together as follows Yo
= GUo
(2)
C learly, if t.he Po urier vectors were observabl e- witho ut errors, th e p robl em would be m f!Tf~ly a mathematical one of so lving thf' SPot. of equatio ns (:2)
(.r., = Yo{;t0 wh~re. (.) t stands for t he pse udo-inverse. Notice that t.he matrixUo has t o be offull row-rank I i.e. m of I.. he input,s Uo(i) have to be linear indep(mdent.
2. PRORLEM fORMU LATION Co nsider a rnultivariable system wi•.h m inputs a nd p outputs. If the "true" inpu t and ou tput signals I uo(t) E R m and yo{t) E R P, were ob.c;crvable at eq uidista.nt time instants t71 HT, 11 U,. " N - I , with T t he sampling period and T = N T the uuse rvat io n pf'r iod, and if the Saml)led signals ('.ould be discrete Fourier t ra nsformed without. introd ucin g errors ( whi ch is t he case for periodic signals) ~ t.h en, fo r a[\ a ngul a. r frequencies in n {Wk = 27fk/T : k 0 , .. . , N / 2 - I} . t.h e following linear mapping wo uld ho lel r:xact.)y
In pra.ct.i ce, however ; errors usually affect the measured Fourier ved.ors , a nd consequ('nt ly, what is obsp, rved a re DOt t.he "t.rue" I/O Fourier data., Zo (UJ' , YOH] H, b ut hO together with some random perturbatio ns l Ez [ E~f. E:]H (see f igure 1). The problem , thus, beco mes statistical, re:;ult ing in random observations I Z = [U Il , yH]H, defined by t.he fo ll owin g Iler1'Or.s.in-variobles" (EV) stochast ic model with M t he number of measurement.s done
=
=
=
Y, YU,I
UO.I
The a im of t his paper is t,o show that. wit h periodic excitati ons (i. t'~. m ult,isilles) and a ppropri at e e.stimators it is poss ible t o ob tain accurate estima.tes of a mullivaria ble f1Y'fJuCll cy response function (M PR F). Experim enta l r l"'sults obtained in the fiel d of vib mtioll analysis dearly demonstrate that. th e new method o logy outperforms t he daosical estimators ba.<:;cd OH spec.t ral analysis and randOIll excitations , especiall y for very ligh t ly damp ed m echan ical str uctures .
=
e ll ,\
E" . ,
=
=
Z(i) Zo(i) + E, (i ) [C;, -I,J Zo(i) 0
(I )
=
1. J
i~ I ,
=
. .. , M :::: m. (4)
U$ii ng th e a bove matrix notati ons, (4) becomes
=
where G G(w.d E C pxm stands for t he fUu ltivnriabl p !ny/u e'lcy response IUllciimt (MFR.P) a.t. angular frequenc:y Wk , while UQ == (la (::';1.) E c m r-tud Yo = YO(Wk ) E CP are t he correspo nding disc.rele:- F'o urie r transfo rms of t.he il true" itlput- 01llpu l (I/ O) signals.
wi t h Z, 20, fz E c (m+p )x M.
Consider now th e sd. ((Uo(i) , \'o(i)) : i = I , .. . , M :::: m} consisting: of M o bsr-rvat.ions of the I/O Fourier vectors at an gu la r freq ueIlcy Wk . Hy grouping the M input a nd output fi'ou rier vectors 1111..0 IwO matrict'S U o [Uo(J), ... ,Uo( M )] E e n" M and Yo = [Yo(I) , .. . , Yo(M)] E
An extens ive a mo unt of statist.ical public.:ations has been devoted to t.he closely related topic of EV regression an alysis (Kcndall and Stuart, 1979 ; Gleser, 1981 ; Fuller, 1987). Dependin g "po n the nature of the vector ZOI random or determini stic, (4) is called in the stat.isticallit-
=
{
4352
Z =Zo + r.. lp] Zo = 0
re;. -
(;,)
erature a st rudural or a funl~tional errors-in-va riables r~lationshir. , respectively. For tht:" s tructura.l relationship , Anderson (1985) pointed out, that when all measured variables, Z, in<.~lud E': ("TrOTS ) E:;, and the secondorder mOlJlents of tht'5e ('.rrors are unknown , it is lIut possible to uniquely identify the f'RF' , (.'. Ind eed, the mnximum likelihood t'.stimule (M LE :l) depends on the. second-order mornent.s of t.he crror.s . Consequently, only a dass of ('andidates for t.he' FRF call be found . Usually, wh e:n an BV model is considered, t he covaria.nce matrix of thE': pE'rturbations, E::r(m) , il:' a.-ssumed to be (I priori known (up to a. constant. multiplicative factor). An EV modd together with this additional assumption is c.ornmonly referred to as the classIcal BV model. Without t.his artditi ollal a:;..<)umpti on (or "prejudice"), the problem is not solved yet (Ljung, 1987, p. ~O:l). In the next section , an nlternative a.'?Sumpt.ioll will be proposed based on repeated observations of tht" signals. This l:all nowadays easily be realized in pra(';tic(" by applying periodic (broadband) excitations (e. g. multisines (Sc.houkens et al., 1993)) and by mea;uring, say. I' > I peri ods of the signals. To show t, hi~, the (unction a l model will be s f.udied in some more ct~t. a il s in the rIP-xI. section.
is assum ed to be (second-order) stationa.ry and unc.:orrelated for different measurements (i "# j) . Moreover, notice t.hat C z accounts for possible correlations am ong the I/O F'ourier <.~oefficients . For a wide class of probability density functions of the hme-domain noise, the Faurier coefficien ts a.re asymptotic.ally (as the number of tim e sam ples IV ----+ 00) complex normally distributed and independent over the frequ encies (Brillinger, (981). This motivates t.he use of a ,,:om plex Gaussian probability density functIOn (POP) t,o l'onstruct t.he MLE . The (asymp tot.ical) iudependencE' over the frequencies allow liS to write the M L equations for a single frequency witho ut loss of (statistical) efficie nc.~y. This r e~~tl lts in the fol· lowing (negative) log- likelihood function
In Append ix A it is shown that {j can be eliminated from (7), which p;i ves
{(Cn =
tr
((BCzilHr ' (BMzBH])
=
:1. ML ESTIMATION OF MFRF'S In the first place, a general functional model is considered without synchronization. This means that all the ~u(i) in (4) can be different. TllC' mnximum likelihood (M L) estimator is conside-red a.'>sumin g the second-order morncnt:i of th e disturban ces Zo (i) t. e-I be a pr1.m'i known, Nf'xt , it will be a'3sumed that the fIlt'.asurements c.an be sy nchronized (i. e., Uo(i) Uo(j) , 'ti , j) which is often true in practice. If thi s is t.h(' case, it turns out that it. is possible to obtain consist€'lnt ML estimat.p-s of t he MFRF' wit.hout, f(~qlliring any info rma.tio n ahont. t.hf': secondorder mOflleo"'s of the. erro rs.
=
=
lVith Mz ZZH and B [6, -Ip] . Minimizing (8) with respect to iJ is equivalf':nt to solving a generalh::ed eigenva.lue prohl~rn (t.he pro\)f for positive-definite Hermit-ian matrices is similar to the one for positive-definit.e symmetric matrices given in Fuller, (Fuller , 1987, C hap . 4))
IMz -
~ Cz l
=0
(9)
The ML estimates are obtained by l~hoosing the columns of BH to be linear combina.t io nt> or the p characterist.ie vectors of M z in th e rne Lric (?z that are associated with the p smallest root.s) Am+ 1 .:; . .. . 2: ).m+p ) of (9). Eventually, {;ML i!:l give.n by
(:~~.)
= -BMd: , 1 : m)IiMd: , m + I : m + p)- ' (10) =
3. 1 A synch 1'Otl ized M e(JSUll!merl ts Although we co nsider th e iuputs lfll(i) to be deterministic: (function a.l rf'la.tiol1ship), the results of this section ("an also be app lif'd to t.ht' ::;trudural model. It will be .,.;umoo that the rando m v",-tors {C., (i) : i = I , ... , M} ar(=' c:omplf'x normally di!:;t,ribukd with
cov(E,(i) , l,,(j))
(8)
=(:ZO' j
while (((lAIc) ~;= J ),m+'. The gen emlized singular valUf~ de(:omposition (GSVD) offern a numerir.a.1 (more) sta.ble way of comp uting (:ML (Golub ancl Van Loan , 1990). Notice that the spectra l density matrix of the I/O no i::;e is required to get eo nsistent estimates of the frequency response matrix. Replacing (.'z by the identity mat rix results in t.h e HT! E"stimator which is not consis.tent ill general (L euridan et al., 1986).
(U)
and Cz E c (m +v»«IU+p) all (J IJfwn known Hermitiansymmetric: c()variance matrix. Hence, the dist.urbing noise 2 Depelldillf,. on l ite c ontext, MLE can s tallci for mar.imum fik t' Jih uud estimat!' as well a s muxim1U11 likd ih :H)d t"stimat or
3.2 Synchronize.d Measuremn!ts When a synchroni ze<.:l measurement set,up is used , the above r ~ ults can be simpliJied and improved by taking thi s informati o n into account . For a multi-excitation
4353
set.up, a. set. of m linear independent stimuli , Uu = (Uo(J). .. . , (lo(m)] E cm "u , aw required . For the dilTerent s timuli , the res ulting I/O cia(,a are measured P tim es (i.e. M mP), which leads to l.hE" followin g EV modE")
=
Z(i),= Zo + f,(i) } [(, , -I,]Zo=O
;=1 , ... ,1'.
( I I)
In Appendix B it. is show n that. t.he MLE red uces in this
case to ( 12)
Notice t.hat. the c.ovariance matrix of I,he disturbances is not requirE'.d any more to obtain consistent estimatE'.s. This estima.tor belongs 1.0 the (:.lass of the rnaxiIlmUJ likelihood estimators ( Kendall and Stll.rt, 1979) for which it has heen prove<.i that they are consistent and asympt.oti(.'a.lly efficient.
4. [)ESJ(;N OF OPTIM AL PERIO[)]C EXCITATIONS Th e quality of t.h e FRF est.imat.f!s do not depends on the e-stirnator on ly~ but. t.he excitation plaYA an important. role too. In many syst.ems (e.q. l~ontrnl systems 1 biological syst.ems) the c.hoice of the excitation is very limit.c:d. However , ill an import.ant. da..<;.<; of problems, only the maximal values of the input and out.put. signals are rest.ricted Lo maintain the linear behavior of the d e.vice under t/:'.I;t and/or to avoid overflow of t he measureme nt eq uipme nL. This freedom can be used to design optimnl t>xpel'i m e nt~ . The dt."Sign of opt.imal ('xc-itations re ly on the minimization of th e C rarner-R.ao lower bound (or , ~llIiva lent l y, the maximizat.ion of t.hf' Fisher information matrix 1 wh ich €':(Iuals the in vers(' nft.he C ramer-R.ao matrix ) .
It can b e shown that t.he Cram er- Ra.o lower bound of the MFRF (':$timat.p. is given by
=
where vec((;) (Ul1, ... ,Gpl,GI2, ... ,Gpmf, C BZ = lICzB H alld B [(;, - Ip]. This reslIlt i, a gene.ralizat.ion of Theorem 8.2.0 in (Brillinger , 198 1) where only o utput. d is~url>a llces are considered. To tler-rease t.he 10wE'.r bound , Olle ha.. t.o minimize det «(:cR(vec(G))) (i. e. Dovtimal design (F'edorov, 197 2)) . whid, equals
=
._ . det«(.c R(vec«(,»)
=
I p mp
det(CBZ )m det(lf,fU;)P
(14)
Fig . 2. Optirnized multisine with peak value normalized to one.
Thus, the optimal inpu~ design is thE'i one which maximizes Idet(Uo}l . Assume that we have d esigned , for a g i v(~ 1I frequency band , a mul tisine with opt.im ized np.st.factor (GuiliallHlc cl 01'1 1991) . Its Fourier coefficient at a.ngular frequency Wk will be d enoted by, sa.y, X. vVe applied this signal (wit.h normalized peak value) to the different inputs taking into aeeollnt. the maximum allowed pp.ak valuE"$ Hi" 1 :::: 1, .. . , m, of th~ input signals tiQ,i,(t). Put /IoU) = (R,XQ'j ,R2X£J2" .. . , RmXQm,f. The qut"~tion t.hat w ill be addres~ed now is how to choose the matrix Q E [-I, l]mxm such that Idet(Uo}l att.ains its ext.remmn. It is readily Vt~rified that the ext.remum equals 2'171 and t.hat this extrcnl1ltn is attained when Q is a I1adarnard matrix of ord er m (Brigg~ and Godfrey, 1966). For iustance, for m
= 21 Q
=
[! ~ll
is
optimal.
5. EXPERIMEN'IAL RESULTS To illus t rat.e the t heoreti cal rE-suits, a lightly damped mechanl r.al structure (,.e. aD a.luminium beam) has been exc.ited at two locations by nleans of mini-shakers and the force and accelNation at. t.hese two points have been measured using pi~zoelect.ri (': impedance heads. An optimized rnultisine (see Figure 'J) containing 12R fr equcncif!s uniformly distributed bd-ween DC and 1800 Hz with a flat. amplitude spectrum has bc:cn used to excite the !3tructure. FirstlYl we: applied the multisine, x(t), to both powpr amplifiers of t he mini·shakers. Th e resulting measurements are given in Pigurereffigure:measpp. After .,ha1. 1 we applied -;1: (t) t.•" the first power amplifier 1 a nd left. the second o ne unchan ged. From these two sets of measurement it is now possible to estimate th e MFRF b.t,lVeen ne and 1800 Hz using (12) with P = I. The resulting nonparamctri(' e:;timat€'s are given ill Figure 1. RepeallIlg this process P t imes and averaging the mea.-
4354
., .
~r,.;.ou£r.l: v 1>1»
I"':R!"01C l ~~
~'.
,"<
>1: -• •
Fig. 6 . H I
li'ig . :l. Measnfe
('st imat~
( Han ning windo w, M -= 10).
Ii. CONC LUSIONS
, ! ~
NiR T.lHCl'i 122)
f. ~
•..•
~
, ~
.~,
,,~
-= .
''''(IIJlN(;V (>~)
Fig. 4. MLr. of the MfRF usillg ( 12) wit,h P
,,,
[n the present papf': r th e maximum likelihood estimation of no n parametric multivariable frequency response fun ctions from noisy measurements of periodic input/output (I /O) signals has b eell considered. T he propose<.i maximum likeli hood ('stimator is consisten t even when the I/O errors are correlated. BE'_"ide.s, no prior noise informatiuJI is required. Moreover , attentiou has been paid to tb e optilllal design of a n expN iment fo r mult i- input systems. Experim ental result:; in the fi eld of vibration analysis illustrat.es th e possib ili ties of the meth od.
L
Fig. 5. lil f'St imat.es (redangula r window 1 lv/ = 40). surements ilC{'ordin g to (12) will further improve the estimates. To be able t u CO lJl pare, fhf:' same measurements have been repeated wi t h uncorrelated noise sources. After 40 averages usin g t.he- das:<;ical (mu lt ivariable) HI estimator one obtains t he res uh.s g ivf'1l in Figure !J . Thebad quality of t he estinH..t.e~ LS m a inl y du e to impo rtant. leakage er rors. Usin g a Hallnin~ window improvH) th e results (see Figure 6).
7. REFEREN CES
Anderson , B. D. O. ( /g85). Id ent.ifi cation of scalar errors-in-variabl es modeb with c1ynamics. Automaticn 21(6), 709- 711), Bendat, .r. S. and A. G. Piersol ( 1980). Engin eering Applir.ation ... of Con'dallou and S~dml A,wlysis . .Jo hn Wiley & Sons. New York . Briggs, P. A. N. and K. R. Godfrey (1966 ). Pscudorandarn :-;ignals for the dynamic analysis of multivariable system'. Proceeding,1 of th e l EE 113 (7), 12591267 Brillinger, D. R. ( 1981), Time Series: Data Annlysis 01"/ 'rlJe,of·Y. expand ed ed .. McG raw-Hill. ~ew York. Cooper , .J. E. and P. R. 8mm ett (1995). A nonpa ra metric approach to instrumentation variables frequency re-sponse fu ndio n estimation. Modal Analysis: 111e In te..maliotlui Journal of Analytical nnr! E>7X'rimental Modal Aunlysis 10(2), 84- 94. Fedorov, V . V. (1 972). Theor!1 of Oplimnl Experim erlls, Acadf>:m i(' Pr('s8. Fuller, W . A. (1 987). Measurrment En'ol' Models. Joh n Wiley & SOilS. New York GI",",r , L ..1. ( 1981) . Estimation in a rnulti variable 'error&-in-variabl es' r egrC!~~ioD model: Largp sample results. The Am.. /s of ,%,tistics 9(1}, 24-44.
4355
Golub, (]. H. and C. F. Van Loan (1990). Matrix Computations. second ed .. The Johns Hopkins University Press. Ra.J t.imo re.. Cuillaume, P., .I . Sehoukens, I-l. Pinl-elo n and I. l
40(6),982- 989. Kend all , M. and A. Stuart (1979). Inference and Relallotl.s/lip. Val. 2 of The. Advcmced Theor,'.J of Statistics. fourth ed .. Charles Griffin & Company Limited . London. Leurid an, .J., D. De Vis, H.
Van
(A .:J) while t he derivative with respect to Re(L) and yield t hat
[mill (A.4)
After elimination of i and L from (A.I) using (A.2)(A.4), one obtains (8).
Der Auweraer
and F Lembregt.s (198f$). A cOlllpari::;on uf
::Iorn~
frt'quen ey respons(' me.
AppcIIlOlJS MEASUREMENT In t.his appendix it is s hown that t.he c:ost functi on~ (7) aud (8) have the- sa me stationary point~.... Thereto tht' met hod of t.he Lagrangian multjpli~rs is appJi{Xi to (7) , r ~s ulting ill
Appendix B. SYNCHRO NOUS MEASUREMENTS CO lls idf'r p
((; , i . L)
= I> (Cz'[Z(i) -
.=,
i ][z(i) - i]H) (B. I)
+ lte(tr(LHBi)) where i E C tnxm . Setting the derivative with respect to i, (; and L t.o zero gives respectively
_2C;;1
(~Z(i) -
Afi) + BH t = 0 =0 Bi = 0
WH
(B.2) (B.3) (B4)
As (j E cm )(m is a regular m at rix by a.ssum pt ion , (8.:3) can only be :latisfied when L r: c mxm equals zero. Consequently, (13 .2) reduces to p
i = +. 'LZ(i)
(R5)
i= 1
and, by way of (B.4), the t raw,fer matrix estimate {;ML bt:.co m es
= rr (C;;I[Z - Z][Z - i]H) +Rc (lr
Similarly, t.he derivative with respect to Re(<")) and Jm(C!J gtve
(A.I)
(LH Hi))
(B.6)
with B = [D, - Fp] and where L ii:i a Lagrangian multiplier matrix for th€': constraint. Hi = 0 (use of thE'! real Notic~ that t.he a-priori knowl edge of Cz is not required part. of th e trace of LH Ri is just a convenient way of a nym ore. to obtaiu tht' ML estim ate of C. summing R.e(L.j) Re(Lk B/,Z'j) and Im(i"j) Jm(L. B,.i. j ) ove r all i,i). In its ext.~em,UlIJ , ({Cl' , z , i) must i1f' ~tat ionary with resp.... t to G . Z and L. Setting iJfjoR.e(Z) + jilf/Ol m( Z) equal to zero gives (notice that is not an a na lytic function , thus . al/ai dON< !lot exist).
e
- 2('i ' [Z - il +
lifl
L= 0
(A .2)
4356
3a-J3 6
ACCUR.ATE ESTIMATION OF MULTIVARIABLE FREQUENCY RESPONSE FUNCTIONS 1 P . G uillaunle, R. Pinteloll and J. Schoukens·
• D"pMimenl EL EC, Vcije Universi!'eil Brussel ( VUB), 8-1050 Brussei, Belgium e-lIlail: [email protected] - fa x: +J2-2-6:!9 2850
Ahstract. Th€' non parametric estimation of Jnultivariable fre.quency response functions is studied within an "errors-in-variables" 8tochastk framework. It is shown that , eve.n when t.he !;e(".ond-ord ~r momentl:! of the dist urbances are unknow n, it is possible by Ilsing periodic t:>xcitations to obtain consistf'.nt: and asymptotically effideut point es t,imat.es of t.h(' multivariable frequency re~p OJlse functions. This is not tht~ ca...:;e with ra.ndorlL oE'. x~~it. a ti ons. Besides , attention is paid to thE': desi gn of c)ptima i mu lti-input exdt.atioHs U1ax.imizing the Fisher information ma.tri x. Eventually, the th e, )retical results arc illust.rated by means of real measurements in the field of vihratio n analy sis.
K eywords. Frequency response, multi variable systems, nonparametric ident.ificat io n , maxim um likelihood, estimation theory.
I INTROD UCTION The nonparamet,rk estiruat.io ll or freq uency response func.t.iolls (or im pulse respo nse fUHd io ns ) is of primary importance ill many scientific inv e~Ligations . This is the reason why mu ch atte.nt.io n ha...:; l.)f~e n paid t o nonpa ra.Il1etric esti mators in th e litt>rahtrf' (Bendat and Pie:rso l, 1880; llrillinger, 1fill I ; All
The first es timat.ors which were developed consid ereti on ly the errors on t.he o utput measurem ents. In man y pra.ctica.1 s it.uations ~ however. the disturbancE':s on th e. input. mea....·m rem cnts can be as large as the ones on the outputs. To cope wit.h t he presence of input as well ~ o u tput. disturba.n ce::; "errors-in-vari ables" estim ators were develop ed. A possible approach is the instrumeut nl vnr'iables (IV) method (see for instance Soderstrom and Sto;"" ( 1989, Chapters 8 and 10) and Cooper and Emmett (1995 )). The IV met. J,od usually requi res an additional data se(luence (i.e. the im;lrumental variables) which can be- considered as a d ra.wback . A positive result is t1w fact that, unde.r some mild condit.ions, the IV est imat.e G m be prove.n to be stro ngly (;onsistent. Another possibili ty is to lU:;e the Hv f'!Stimator which rely on t he t.ot.a.l leas1. sq llar f'~o;; met hod. T he Ht) est.imates arc cOllsistent when the input/outpul (I/O) errors are all equall y large all J uncorrelated . If this is not the case, consistent estimates ('an st.ill be obtain~d by introduci ng a. scaling mat.ri x resulting in a gene ralized tot.alleast squares
4351
pro blem (Van Hu ffel and Vandewalle, 199 1). To derive t hls sealing matri x, howev~r . the "true" covariance m atr ix of th(" input. a nd output. errors s hould be a priori kn own. In th e presen t pa.per t.he non parame tric estimatio n of frequenc.y rt'ponse function wi ll be revisited ta king into accou nt th e specifk aspects of periodic excit aLions , i.e.. the availa.bility of replicat.ed obsf'rvations of the "true" I/ O da.ta. This wi ll result. in an a lmost t ri vial solution; a trivial solu tions whidl is , however , more perrOfmant. than t he. com mo nl y used estimators. Morp.-over, consistent est imates are obtai n even if t he "t rue" cava fiance matrix of t he I/O f:'rrors is unknown . Furt her , it will be s hown that. t he pro poS<'d estimator is a.'.iy rnpto t ically efficient. Confidence intervals for the estimated fr('qu ency res ponse functi ons will be given. Next, a.ttent ion will be s penl 1.0 the design of o ptimal mul t i-inpu t cxcitations.
u, £1/. 1
Eu . 2
U,
G=
UO , 2
[V"
G,~ ]
V,
(J~",
G~ l
E y •p
E .. , ",
y.
u~
Uo""
Fig.
I . Multivariable frequen cy-domain vari ables" set up.
"errors-in-
C pxM J relation (I) can h e expressed fo r al1 t he data together as follows Yo
= GUo
(2)
C learly, if t.he Po urier vectors were observabl e- witho ut errors, th e p robl em would be m f!Tf~ly a mathematical one of so lving thf' SPot. of equatio ns (:2)
(.r., = Yo{;t0 wh~re. (.) t stands for t he pse udo-inverse. Notice that t.he matrixUo has t o be offull row-rank I i.e. m of I.. he input,s Uo(i) have to be linear indep(mdent.
2. PRORLEM fORMU LATION Co nsider a rnultivariable system wi•.h m inputs a nd p outputs. If the "true" inpu t and ou tput signals I uo(t) E R m and yo{t) E R P, were ob.c;crvable at eq uidista.nt time instants t71 HT, 11 U,. " N - I , with T t he sampling period and T = N T the uuse rvat io n pf'r iod, and if the Saml)led signals ('.ould be discrete Fourier t ra nsformed without. introd ucin g errors ( whi ch is t he case for periodic signals) ~ t.h en, fo r a[\ a ngul a. r frequencies in n {Wk = 27fk/T : k 0 , .. . , N / 2 - I} . t.h e following linear mapping wo uld ho lel r:xact.)y
In pra.ct.i ce, however ; errors usually affect the measured Fourier ved.ors , a nd consequ('nt ly, what is obsp, rved a re DOt t.he "t.rue" I/O Fourier data., Zo (UJ' , YOH] H, b ut hO together with some random perturbatio ns l Ez [ E~f. E:]H (see f igure 1). The problem , thus, beco mes statistical, re:;ult ing in random observations I Z = [U Il , yH]H, defined by t.he fo ll owin g Iler1'Or.s.in-variobles" (EV) stochast ic model with M t he number of measurement.s done
=
=
=
Y, YU,I
UO.I
The a im of t his paper is t,o show that. wit h periodic excitati ons (i. t'~. m ult,isilles) and a ppropri at e e.stimators it is poss ible t o ob tain accurate estima.tes of a mullivaria ble f1Y'fJuCll cy response function (M PR F). Experim enta l r l"'sults obtained in the fiel d of vib mtioll analysis dearly demonstrate that. th e new method o logy outperforms t he daosical estimators ba.<:;cd OH spec.t ral analysis and randOIll excitations , especiall y for very ligh t ly damp ed m echan ical str uctures .
=
e ll ,\
E" . ,
=
=
Z(i) Zo(i) + E, (i ) [C;, -I,J Zo(i) 0
(I )
=
1. J
i~ I ,
=
. .. , M :::: m. (4)
U$ii ng th e a bove matrix notati ons, (4) becomes
=
where G G(w.d E C pxm stands for t he fUu ltivnriabl p !ny/u e'lcy response IUllciimt (MFR.P) a.t. angular frequenc:y Wk , while UQ == (la (::';1.) E c m r-tud Yo = YO(Wk ) E CP are t he correspo nding disc.rele:- F'o urie r transfo rms of t.he il true" itlput- 01llpu l (I/ O) signals.
wi t h Z, 20, fz E c (m+p )x M.
Consider now th e sd. ((Uo(i) , \'o(i)) : i = I , .. . , M :::: m} consisting: of M o bsr-rvat.ions of the I/O Fourier vectors at an gu la r freq ueIlcy Wk . Hy grouping the M input a nd output fi'ou rier vectors 1111..0 IwO matrict'S U o [Uo(J), ... ,Uo( M )] E e n" M and Yo = [Yo(I) , .. . , Yo(M)] E
An extens ive a mo unt of statist.ical public.:ations has been devoted to t.he closely related topic of EV regression an alysis (Kcndall and Stuart, 1979 ; Gleser, 1981 ; Fuller, 1987). Dependin g "po n the nature of the vector ZOI random or determini stic, (4) is called in the stat.isticallit-
=
{
4352
Z =Zo + r.. lp] Zo = 0
re;. -
(;,)
erature a st rudural or a funl~tional errors-in-va riables r~lationshir. , respectively. For tht:" s tructura.l relationship , Anderson (1985) pointed out, that when all measured variables, Z, in<.~lud E': ("TrOTS ) E:;, and the secondorder mOlJlents of tht'5e ('.rrors are unknown , it is lIut possible to uniquely identify the f'RF' , (.'. Ind eed, the mnximum likelihood t'.stimule (M LE :l) depends on the. second-order mornent.s of t.he crror.s . Consequently, only a dass of ('andidates for t.he' FRF call be found . Usually, wh e:n an BV model is considered, t he covaria.nce matrix of thE': pE'rturbations, E::r(m) , il:' a.-ssumed to be (I priori known (up to a. constant. multiplicative factor). An EV modd together with this additional assumption is c.ornmonly referred to as the classIcal BV model. Without t.his artditi ollal a:;..<)umpti on (or "prejudice"), the problem is not solved yet (Ljung, 1987, p. ~O:l). In the next section , an nlternative a.'?Sumpt.ioll will be proposed based on repeated observations of tht" signals. This l:all nowadays easily be realized in pra(';tic(" by applying periodic (broadband) excitations (e. g. multisines (Sc.houkens et al., 1993)) and by mea;uring, say. I' > I peri ods of the signals. To show t, hi~, the (unction a l model will be s f.udied in some more ct~t. a il s in the rIP-xI. section.
is assum ed to be (second-order) stationa.ry and unc.:orrelated for different measurements (i "# j) . Moreover, notice t.hat C z accounts for possible correlations am ong the I/O F'ourier <.~oefficients . For a wide class of probability density functions of the hme-domain noise, the Faurier coefficien ts a.re asymptotic.ally (as the number of tim e sam ples IV ----+ 00) complex normally distributed and independent over the frequ encies (Brillinger, (981). This motivates t.he use of a ,,:om plex Gaussian probability density functIOn (POP) t,o l'onstruct t.he MLE . The (asymp tot.ical) iudependencE' over the frequencies allow liS to write the M L equations for a single frequency witho ut loss of (statistical) efficie nc.~y. This r e~~tl lts in the fol· lowing (negative) log- likelihood function
In Append ix A it is shown that {j can be eliminated from (7), which p;i ves
{(Cn =
tr
((BCzilHr ' (BMzBH])
=
:1. ML ESTIMATION OF MFRF'S In the first place, a general functional model is considered without synchronization. This means that all the ~u(i) in (4) can be different. TllC' mnximum likelihood (M L) estimator is conside-red a.'>sumin g the second-order morncnt:i of th e disturban ces Zo (i) t. e-I be a pr1.m'i known, Nf'xt , it will be a'3sumed that the fIlt'.asurements c.an be sy nchronized (i. e., Uo(i) Uo(j) , 'ti , j) which is often true in practice. If thi s is t.h(' case, it turns out that it. is possible to obtain consist€'lnt ML estimat.p-s of t he MFRF' wit.hout, f(~qlliring any info rma.tio n ahont. t.hf': secondorder mOflleo"'s of the. erro rs.
=
=
lVith Mz ZZH and B [6, -Ip] . Minimizing (8) with respect to iJ is equivalf':nt to solving a generalh::ed eigenva.lue prohl~rn (t.he pro\)f for positive-definite Hermit-ian matrices is similar to the one for positive-definit.e symmetric matrices given in Fuller, (Fuller , 1987, C hap . 4))
IMz -
~ Cz l
=0
(9)
The ML estimates are obtained by l~hoosing the columns of BH to be linear combina.t io nt> or the p characterist.ie vectors of M z in th e rne Lric (?z that are associated with the p smallest root.s) Am+ 1 .:; . .. . 2: ).m+p ) of (9). Eventually, {;ML i!:l give.n by
(:~~.)
= -BMd: , 1 : m)IiMd: , m + I : m + p)- ' (10) =
3. 1 A synch 1'Otl ized M e(JSUll!merl ts Although we co nsider th e iuputs lfll(i) to be deterministic: (function a.l rf'la.tiol1ship), the results of this section ("an also be app lif'd to t.ht' ::;trudural model. It will be .,.;umoo that the rando m v",-tors {C., (i) : i = I , ... , M} ar(=' c:omplf'x normally di!:;t,ribukd with
cov(E,(i) , l,,(j))
(8)
=(:ZO' j
while (((lAIc) ~;= J ),m+'. The gen emlized singular valUf~ de(:omposition (GSVD) offern a numerir.a.1 (more) sta.ble way of comp uting (:ML (Golub ancl Van Loan , 1990). Notice that the spectra l density matrix of the I/O no i::;e is required to get eo nsistent estimates of the frequency response matrix. Replacing (.'z by the identity mat rix results in t.h e HT! E"stimator which is not consis.tent ill general (L euridan et al., 1986).
(U)
and Cz E c (m +v»«IU+p) all (J IJfwn known Hermitiansymmetric: c()variance matrix. Hence, the dist.urbing noise 2 Depelldillf,. on l ite c ontext, MLE can s tallci for mar.imum fik t' Jih uud estimat!' as well a s muxim1U11 likd ih :H)d t"stimat or
3.2 Synchronize.d Measuremn!ts When a synchroni ze<.:l measurement set,up is used , the above r ~ ults can be simpliJied and improved by taking thi s informati o n into account . For a multi-excitation
4353
set.up, a. set. of m linear independent stimuli , Uu = (Uo(J). .. . , (lo(m)] E cm "u , aw required . For the dilTerent s timuli , the res ulting I/O cia(,a are measured P tim es (i.e. M mP), which leads to l.hE" followin g EV modE")
=
Z(i),= Zo + f,(i) } [(, , -I,]Zo=O
;=1 , ... ,1'.
( I I)
In Appendix B it. is show n that. t.he MLE red uces in this
case to ( 12)
Notice t.hat. the c.ovariance matrix of I,he disturbances is not requirE'.d any more to obtain consistent estimatE'.s. This estima.tor belongs 1.0 the (:.lass of the rnaxiIlmUJ likelihood estimators ( Kendall and Stll.rt, 1979) for which it has heen prove<.i that they are consistent and asympt.oti(.'a.lly efficient.
4. [)ESJ(;N OF OPTIM AL PERIO[)]C EXCITATIONS Th e quality of t.h e FRF est.imat.f!s do not depends on the e-stirnator on ly~ but. t.he excitation plaYA an important. role too. In many syst.ems (e.q. l~ontrnl systems 1 biological syst.ems) the c.hoice of the excitation is very limit.c:d. However , ill an import.ant. da..<;.<; of problems, only the maximal values of the input and out.put. signals are rest.ricted Lo maintain the linear behavior of the d e.vice under t/:'.I;t and/or to avoid overflow of t he measureme nt eq uipme nL. This freedom can be used to design optimnl t>xpel'i m e nt~ . The dt."Sign of opt.imal ('xc-itations re ly on the minimization of th e C rarner-R.ao lower bound (or , ~llIiva lent l y, the maximizat.ion of t.hf' Fisher information matrix 1 wh ich €':(Iuals the in vers(' nft.he C ramer-R.ao matrix ) .
It can b e shown that t.he Cram er- Ra.o lower bound of the MFRF (':$timat.p. is given by
=
where vec((;) (Ul1, ... ,Gpl,GI2, ... ,Gpmf, C BZ = lICzB H alld B [(;, - Ip]. This reslIlt i, a gene.ralizat.ion of Theorem 8.2.0 in (Brillinger , 198 1) where only o utput. d is~url>a llces are considered. To tler-rease t.he 10wE'.r bound , Olle ha.. t.o minimize det «(:cR(vec(G))) (i. e. Dovtimal design (F'edorov, 197 2)) . whid, equals
=
._ . det«(.c R(vec«(,»)
=
I p mp
det(CBZ )m det(lf,fU;)P
(14)
Fig . 2. Optirnized multisine with peak value normalized to one.
Thus, the optimal inpu~ design is thE'i one which maximizes Idet(Uo}l . Assume that we have d esigned , for a g i v(~ 1I frequency band , a mul tisine with opt.im ized np.st.factor (GuiliallHlc cl 01'1 1991) . Its Fourier coefficient at a.ngular frequency Wk will be d enoted by, sa.y, X. vVe applied this signal (wit.h normalized peak value) to the different inputs taking into aeeollnt. the maximum allowed pp.ak valuE"$ Hi" 1 :::: 1, .. . , m, of th~ input signals tiQ,i,(t). Put /IoU) = (R,XQ'j ,R2X£J2" .. . , RmXQm,f. The qut"~tion t.hat w ill be addres~ed now is how to choose the matrix Q E [-I, l]mxm such that Idet(Uo}l att.ains its ext.remmn. It is readily Vt~rified that the ext.remum equals 2'171 and t.hat this extrcnl1ltn is attained when Q is a I1adarnard matrix of ord er m (Brigg~ and Godfrey, 1966). For iustance, for m
= 21 Q
=
[! ~ll
is
optimal.
5. EXPERIMEN'IAL RESULTS To illus t rat.e the t heoreti cal rE-suits, a lightly damped mechanl r.al structure (,.e. aD a.luminium beam) has been exc.ited at two locations by nleans of mini-shakers and the force and accelNation at. t.hese two points have been measured using pi~zoelect.ri (': impedance heads. An optimized rnultisine (see Figure 'J) containing 12R fr equcncif!s uniformly distributed bd-ween DC and 1800 Hz with a flat. amplitude spectrum has bc:cn used to excite the !3tructure. FirstlYl we: applied the multisine, x(t), to both powpr amplifiers of t he mini·shakers. Th e resulting measurements are given in Pigurereffigure:measpp. After .,ha1. 1 we applied -;1: (t) t.•" the first power amplifier 1 a nd left. the second o ne unchan ged. From these two sets of measurement it is now possible to estimate th e MFRF b.t,lVeen ne and 1800 Hz using (12) with P = I. The resulting nonparamctri(' e:;timat€'s are given ill Figure 1. RepeallIlg this process P t imes and averaging the mea.-
4354
., .
~r,.;.ou£r.l: v 1>1»
I"':R!"01C l ~~
~'.
,"<
>1: -• •
Fig. 6 . H I
li'ig . :l. Measnfe
('st imat~
( Han ning windo w, M -= 10).
Ii. CONC LUSIONS
, ! ~
NiR T.lHCl'i 122)
f. ~
•..•
~
, ~
.~,
,,~
-= .
''''(IIJlN(;V (>~)
Fig. 4. MLr. of the MfRF usillg ( 12) wit,h P
,,,
[n the present papf': r th e maximum likelihood estimation of no n parametric multivariable frequency response fun ctions from noisy measurements of periodic input/output (I /O) signals has b eell considered. T he propose<.i maximum likeli hood ('stimator is consisten t even when the I/O errors are correlated. BE'_"ide.s, no prior noise informatiuJI is required. Moreover , attentiou has been paid to tb e optilllal design of a n expN iment fo r mult i- input systems. Experim ental result:; in the fi eld of vibration analysis illustrat.es th e possib ili ties of the meth od.
L
Fig. 5. lil f'St imat.es (redangula r window 1 lv/ = 40). surements ilC{'ordin g to (12) will further improve the estimates. To be able t u CO lJl pare, fhf:' same measurements have been repeated wi t h uncorrelated noise sources. After 40 averages usin g t.he- das:<;ical (mu lt ivariable) HI estimator one obtains t he res uh.s g ivf'1l in Figure !J . Thebad quality of t he estinH..t.e~ LS m a inl y du e to impo rtant. leakage er rors. Usin g a Hallnin~ window improvH) th e results (see Figure 6).
7. REFEREN CES
Anderson , B. D. O. ( /g85). Id ent.ifi cation of scalar errors-in-variabl es modeb with c1ynamics. Automaticn 21(6), 709- 711), Bendat, .r. S. and A. G. Piersol ( 1980). Engin eering Applir.ation ... of Con'dallou and S~dml A,wlysis . .Jo hn Wiley & Sons. New York . Briggs, P. A. N. and K. R. Godfrey (1966 ). Pscudorandarn :-;ignals for the dynamic analysis of multivariable system'. Proceeding,1 of th e l EE 113 (7), 12591267 Brillinger, D. R. ( 1981), Time Series: Data Annlysis 01"/ 'rlJe,of·Y. expand ed ed .. McG raw-Hill. ~ew York. Cooper , .J. E. and P. R. 8mm ett (1995). A nonpa ra metric approach to instrumentation variables frequency re-sponse fu ndio n estimation. Modal Analysis: 111e In te..maliotlui Journal of Analytical nnr! E>7X'rimental Modal Aunlysis 10(2), 84- 94. Fedorov, V . V. (1 972). Theor!1 of Oplimnl Experim erlls, Acadf>:m i(' Pr('s8. Fuller, W . A. (1 987). Measurrment En'ol' Models. Joh n Wiley & SOilS. New York GI",",r , L ..1. ( 1981) . Estimation in a rnulti variable 'error&-in-variabl es' r egrC!~~ioD model: Largp sample results. The Am.. /s of ,%,tistics 9(1}, 24-44.
4355
Golub, (]. H. and C. F. Van Loan (1990). Matrix Computations. second ed .. The Johns Hopkins University Press. Ra.J t.imo re.. Cuillaume, P., .I . Sehoukens, I-l. Pinl-elo n and I. l
40(6),982- 989. Kend all , M. and A. Stuart (1979). Inference and Relallotl.s/lip. Val. 2 of The. Advcmced Theor,'.J of Statistics. fourth ed .. Charles Griffin & Company Limited . London. Leurid an, .J., D. De Vis, H.
Van
(A .:J) while t he derivative with respect to Re(L) and yield t hat
[mill (A.4)
After elimination of i and L from (A.I) using (A.2)(A.4), one obtains (8).
Der Auweraer
and F Lembregt.s (198f$). A cOlllpari::;on uf
::Iorn~
frt'quen ey respons(' me.
AppcII
Appendix B. SYNCHRO NOUS MEASUREMENTS CO lls idf'r p
((; , i . L)
= I> (Cz'[Z(i) -
.=,
i ][z(i) - i]H) (B. I)
+ lte(tr(LHBi)) where i E C tnxm . Setting the derivative with respect to i, (; and L t.o zero gives respectively
_2C;;1
(~Z(i) -
Afi) + BH t = 0 =0 Bi = 0
WH
(B.2) (B.3) (B4)
As (j E cm )(m is a regular m at rix by a.ssum pt ion , (8.:3) can only be :latisfied when L r: c mxm equals zero. Consequently, (13 .2) reduces to p
i = +. 'LZ(i)
(R5)
i= 1
and, by way of (B.4), the t raw,fer matrix estimate {;ML bt:.co m es
= rr (C;;I[Z - Z][Z - i]H) +Rc (lr
Similarly, t.he derivative with respect to Re(<")) and Jm(C!J gtve
(A.I)
(LH Hi))
(B.6)
with B = [D, - Fp] and where L ii:i a Lagrangian multiplier matrix for th€': constraint. Hi = 0 (use of thE'! real Notic~ that t.he a-priori knowl edge of Cz is not required part. of th e trace of LH Ri is just a convenient way of a nym ore. to obtaiu tht' ML estim ate of C. summing R.e(L.j) Re(Lk B/,Z'j) and Im(i"j) Jm(L. B,.i. j ) ove r all i,i). In its ext.~em,UlIJ , ({Cl' , z , i) must i1f' ~tat ionary with resp.... t to G . Z and L. Setting iJfjoR.e(Z) + jilf/Ol m( Z) equal to zero gives (notice that is not an a na lytic function , thus . al/ai dON< !lot exist).
e
- 2('i ' [Z - il +
lifl
L= 0
(A .2)
4356