Hansen, Lars Kai

Pedersen, Thorkild Find

Online Research Database In Technology

BAYESIAN MULTICHANNEL TRACKING OF A PERIODIC SIGNALThorkild Find Pedersen∗Bru¨el & Kjær Sound & Vibration A/SSkodsborgvej 307DK-2850 Nærum, Denmarktfpedersen@bksv.comLars Kai HansenInformatics and Mathematical ModellingTechnical University of DenmarkRichard Petersens Plads, Building 321DK-2800 Kongens Lyngby, Denmarklkhansen@imm.dtu.dkABSTRACTThe problem of estimating and tracking the fundamentalfrequency of a periodic signal is formulated in Bayesianterms and a solution is proposed. For multichannel mea-surements both simulated and experimental data show thatthe estimate improves when more data channels are used forthe estimation. The method is tested on a multichannel vi-bration measurement from the interior of a passenger car.From the measurements the running speed of the engine issuccessfully estimated during a run-up/down experiment.1. INTRODUCTIONWhen analyzing the behaviour of rotating machines, e.g.the vibrations casued by the engine in a passenger car, itis important to know the running speed, or fundamentalfrequency, of the engine. The fundamental frequency istypically measured using dedicated sensors like proximityprobes or photosensors which require direct access to the ro-tating part of the machine. Rotating parts are often not easyaccessible. It is therefore of interest to investigate the pos-siblity of extracting the fundamental frequency from othersources, e.g. the vibration signals being analyzed.In section 2 Bayesian spectral estimation theory [1, 2] isreviewed in the light of estimating the posterior probabilityof the fundamental frequency, ωF , conditioned on the mea-sured data, i.e. p(ωF |Di, Ik) where Ik symbolizes the priorknowledge in the estimation model, and Di represent one ormore of the measured signals. Only a fixed model is consid-ered here, but within the Bayesian framework it expandableto multiple models.By segmenting the source signals into overlapping timerecords and introducing a prior on ωF conditioned on theprevious two records, the estimator is extended to includetracking, i.e p(ωF (n)|ωF (n−1), ωF (n−2)), where n is thenumber of the current record.∗This work is funded by the Danish Academi of Technical Sciencesthrough the industrial Ph.d fellowship EF-812The fundamental frequency estimator and tracker is firsttested on a simulated mechanical system with one degreeof freedom, desribed in section 3. Then the experimentaldata described in section 4 are used for validation. The re-sults from both simulated and experimental data are foundin section 5 and 6.2. BAYESIAN ESTIMATION AND TRACKINGIn the litterature [3, 4] the problem of detecting and trackingperiodic signals is often defined as estimating the amplitudeand phase of limited number of harmonic components, i.ed(t) =K∑k=1(ak cos(ωkt) + bk sin(ωkt)) + e(t) (1)e(t) ∼ N (0, σ2) (2)The frequencies in the harmonic sequence are defined as or-ders αk of the fundamental frequency, ω0, i.e. ωk = αkω0.In [4] the asymptotic Crame´r Rao Lower bound is derivedfor this type of harmonic signal:Var(ωˆF ) ≥ 12SNRω2eff N3(3)where SNR =∑nk=1 c2k/(2σ2N ), c2k = a2k + b2k, and ω2eff =∑nk=1 k2c2k/∑nk=1 c2kWriting (1) in vectorized form, it is seen to be a linearmodeld = G · b+ e (4)where d is the observed data, G is a matrix containing basisvectors as columns, b is the parameter vector and e is anoise vector:dN×1 = [d0, d1, . . . , dN−1]>eN×1 = [e0, e1, . . . , eN−1]>b2K×1 = [a1, b1, . . . , ak, bK ]>GN×2K = [cos(Ω1), sin(Ω1), . . . , cos(ΩK), sin(ΩK)]Ωi = [ωi0, ωi1, . . . , ωi(N − 1)]>(5)The likelihood of the observed data given all parameters andthe characteristics of the noise isp(d|b,G, IK)= (2piσ2)−N/2 exp{− (d−G · b)> · (d−G · b)2σ2}and the Maximum Likelihood estimate of the parameters,b, isbˆ = (G>G)−1G> · dSince the likelihood expresses the probability of the datagiven the parameters, Bayes Rule is used to express theprobability of the parameters, θ = [ω0, σ,b], given the data,d, and the model, IK :p(θ|d, IK) = p(d|θ, IK)p(θ|IK)p(d|IK) (6)Of the parameters in θ only ωF is of interest, therefore un-informative prior distributions are assigned to the nuisanceparameters, σ and b:p(b|IK) = kb p(σ|IK) = kσσ(7)Integrating out the nuisance parameters, the posterior distri-bution of ωF is [1]p(ωF |d, IK) ∝∫ ∞0∫<2Kp(d|ω0, σ,b, IK) 1σdb dσ(8)∝(d>d− f>f)−(N−2K)2√det(G>G)(9)where, f is the estimated signal based on the maximum like-lihood estimate, i.e.f = G · bˆ (10)The posterior is recognized as a student t-distribution.2.1. Multichannel MeasurementsWhen a measurement contains multiple data channels,D1, D2, . . . , DN , then if the data channels are assumed tobe independent, i.e. p(D1, D2, . . . , DN ) =∏Ni p(Di), themeasurements can be combined through the joint proba-bility for the fundamental frequency ωF given all the datachannels, I.e.p(ωF |D1, D2, . . . , DN ) ∝N∏i=1p(ωF |Di) (11)2.2. TrackingThe posteriori in (8) is valid for stationary signals only. Inorder to track ωF the measurment data is segmented intooverlapping records of suitable lengths where ωF is close tostationary.The posteriori p(ωF (n)|Dx(n), Ik) is estimated for eachdata segment (denoted by the segment number, n) and theproblem of tracking ωF is now solved by adding a trackingprior that conditions ωF (n) on e.g. the previous two record,i.e p(ωF (n)|ωF (n−1), ωF (n−2)). Assuming the first orderderivative of ωF is constant, i.e.ωF (n)− ωF (n− 1) =ωF (n− 1)− ωF (n− 2) +N (0, σT 2) (12)the tracking prior becomesp(ωF (n)|ωF (n− 1), ωF (n− 2)) =1√2piσT 2e−(ωF (n)−2ωF (n−1)+ωF (n−2))2/2σT2 (13)The combined estimation and tracking algorithm is listen intable 2.3. SIMULATION DATAThe vibration signals are here modelled as the forced re-sponse of a mechanical system with a single degree of free-dom (SDOF). The equation of motion is governed bymx¨(t) + c x˙(t) + k x(t) = F (t) (14)where m is the mass, c is the damping, k is the spring forc-ing and F is the external forcing on the mass. (˙) and (¨ ) in-dicates single and double differentiation in time. x(t) is thedisplacement of the mass. The frequency response for dis-placement, velocity and acceleration of the system is shownin figure 1 together with the system. The mechanical sys-tem is simulated with a digital filter, which can be derivedfrom the Laplace transform of (14) using the bilinear trans-for. The driving force is a harmonic signal composed ofthe first three harmonics of the fundamental frequency (15).The amplitude of the harmonic components are inverselyproportional to the order of the component.x(t) = A3∑n=11ncosnΩ(t), Ω(t) =∫ t0ωF (τ)dτ (15)As is shown in figure 2 the relative amplitude of the orderschange as ωF changes and the orders pass through the reso-nance of the system. The top graph shows the profile of theharmonic orders of ωF . The lower graph shows the RMSacceleration amplitudes of the harmonic orders1.1The units used have no physical interpretation.White Gaussian noise with a constant variance, σN 2, isadded to the simulated response of the mechanical system.The simulated response thus becomes with ? denoting con-volution and hr(t) being the impulse response of the systemmodelled . The noise variance is determined from an overallSNR level defined as SNRdb = 10 log(Var(y)/σ2N ). wherey is the response. Since the effect of y(t) varies over time,so will the actual SNR level which is also observed in prac-tice. In the simulations the three responses, acceleration,velocity and displacement are used as independent sourcesignals.4. EXPERIMENTAL DATAThe experimental data consists of three vibration channelsand one tacho channel for reference. The vibration signalsare measured with a tri-axial accelerometer placed on themounting bolt of the front seat in a passenger car. The sig-nals are labeled Dx, Dy and Dz . The tacho channel wasconnected to an photo-sensitive probe emitting one pulseper revolution on the crankshaft, being equivalent to halfthe fundamental frequency (it was a 4 stroke engine). Alldata was recorded with a B&K Pulse multi-annalyzer andthe digitized data was processed using Matlab. For the pur-pose of finding ωF the exact scaling of the signals are not ofimportance, therefore the digitized values of the conditionedanalog signals have been used unscaled.During the measurement the car is in neutral and therunning speed is increased in steps; for the first 15 secondsthe engine is running approximately 1000 RPM, for the next25 seconds approx. 2000 RPM, then for 10 seconds 3500RPM whereafter it drops to idle for the remaining 10 sec-onds. Figure 5 shows a typical spectrogram of one of thevibration signals. The signal is seen to contain many higherorder harmonics of the RPM profile, but no specific orderstands-out clearly. In general the signal is very noisy. Thevibration signals were sampled at 1024 samples per second,while the tacho signal was sampled at 65536 samples persecond and linear interpolation2 was used to improve theaccuracy of the ωF reference.5. SIMULATION RESULTS100 Monte Carlo simulations are made of the simulated vi-bration signals. Each signal is segmented into records of 1/4second with 75% overlap. In each simulation, the posteriordistribution of ωF is computed for the harmonic sequenceK = [1, 2, 3].The posteriors of ωF for each signal are computed, andthe joint posterior formed. The tracking prior is initialized2Linear interpolation does not give the best accuracy, but sufficient forthese comparisons. Better precision is obtained by opsampling the tachosignal and using higher order interpolation schemes [5].withp(ωF (0)) = N (10, 1), σT 2 = 1.0for each posterior distribution. The results are summa-rized in figure 4, where the standard deviation of the ωF es-timates are plotted with symbols; with lines are shown theCrame´r Rao lower bounds. The simulated results are veryclose to the CRLB which indicates that the Bayesian trackeris optimal. The tracker was initialized with6. EXPERIMENTAL RESULTSEach of vibration signals, Dx, Dy, Dz were segmented intorecords of 1/4 second, overlapping 75%. For each seg-ment the posterior of ωF was computed for the harmonicsequence K = [1, 2, 2.5]. The three posteriors were thencombined. The initial value for the tracker was set top(ωF (0)) = N (12, 3), σT 2 = 0.8Figure 7 shows the result of tracking on the individual pos-terior distributions p(ωF |Dx), p(ωF |Dy), p(ωF |Dz) and thejoint posterior, p(ωF |Dx, Dy, Dz). It is seen that each of theindividual tracks fail at different times: Dx fails initially totrack at 1000 RPM the first 10 seconds; Dy cannot keep upwith the jump from 2000 to 3500 RPM and looses track; Dzlooses track de-accelerating from 3500 RPM and does notreturn to idle. The tracking error is defined as the differencebetween the obtained track and the measured tacho refer-ence. Table 1 summarizes the result of tracking on all pos-terior combinations of the source signals. The table showsthe mean error in column 2, the standard deviation σ in col-umn 3 and the maximum absolute error in column 4. Thetable shows convincingly that the estimate improves as thesignals are combined.7. CONCLUSIONHere the Bayesian paradigm has proven succesful in de-vising a method to estimate and track the fundamental fre-quency of a rotating machine from a multichannel vibrationmeasurement.Starting with the conditioned posterior distribution ofωF given measured data, the model was first extended froma single channel estimation to multichannel estimation viathe formulation of the joint probability of the data. Then byintroduction of prior distribution for ωF (t) given previousestimates, tracking was readily at hand.The method was tested on simulated and experimentaldata with great success. Simulations showed that the esti-mate meets the Crame´r Rao lower bound, and the experi-mental data showed the methods usefulnes in solving thepractical problem of determining the running speed of theengine in a passenger car. Not only was the method capableof providing accurate estimates, but also able to track rapidspeed changes.For the method to be of practical use however, there assome issues to be addressed. Choosing the right trackingvariance can be difficult, perhaps this could be solved usingthe Baum-Welch algorithm. Choosing the most appropriateharmonic sequence for the spectrum estimator is important,the Bayesian paradigm allows the method e.g. to be ex-panded to operate with multiple harmonic sequences. Fur-thermore reversible jump MCMC procedures [6] could beinvestigated for selecting the best harmonic model.8. REFERENCES[1] Joseph J.K. ´O Ruanaidh and William J. Fitzgerald, Nu-merical Bayesian Methods Applied to Signal Process-ing, Statistics and Computing. Springer, 1998.[2] G. L. Bretthorst, Bayhesian Spectrum Analysis andParameter Estimation, vol. 48 of Lecture Notes inStatistics, Springer-Velag, 1988, Available online:www.bayes.wustl.edu.[3] P.J. Parker and B.D.O. Anderson, “Frequency trackingof nonsinusoidal periodic signals in noise,” Signal Pro-cessing, vol. 20, no. 2, pp. 127–52, 1990.[4] A. Nehorai and B. Porat, “Adaptive comb filtering forharmonic signal enhancement,” IEEE Transactions onAcoustics, Speech and Signal Processing, vol. ASSP-34, no. 5, pp. 1124–38, 1986.[5] K. R. Fyfe and E. D. S. Munck, “Analysis of computedorder tracking,” Mechanical Systems & Signal Process-ing, vol. 11, no. 2, pp. 187–205, 1997.[6] C. Andrieu and A. Doucet, “Joint bayesian model se-lection and estimation of noisy sinusoids via reversiblejump mcmc,” Signal Processing, IEEE Transactionson, vol. 47, no. 10, pp. 2667 –2676, 1999.Fig. 1. Frequency responses for the second order filter sim-ulating the mechanical system with a single degree of free-dom (insert): Acc - acceleration respones; Vel - velocityresponse and Disp - displacement response.Data Mean Error σ Max. ErrorDx −1.3 3.1 9.1Dy 3.3 18 64Dz 0.33 1.7 9.9Dxy 0.057 0.6 5.1Dxz 0.17 1.6 9.9Dyz 0.099 0.77 5.8Dxyz 0.059 0.57 5.1Table 1. Tracking Results. From the three source signals atotal of 7 posterior distributions can be formed. The result oftracking on each of the combinations are summarized here.For each track the mean tracking error, its standard devi-ation (σ) and the maximum absolute error is shown. Sig-nificant improvements are observed when the sources arecombined.For each data segment n do:• Compute p(ωF |Di, Ik) for all channels i = 1 . . . N• Compute joint posterior, p(ωF |D1, . . . , DN , Ik)• ωˆF (n) = maxωFp(ωF |D1, . . . , DN , Ik)p(ωF |ωˆF (n− 1), ωˆF (n− 2))EndTable 2. Tracking AlgorithmFig. 2. Harmonic Orders of periodic signal. Top graph:the three curves show frequency profile of the threee har-monic components in the driving force of the mechanicalsystem. The harmonic orders decrease 3dB per octave. Bot-tom graph: RMS amplitudes of the acceleration responseof the driving force. It shows that the relative amplitudeof the harmonic components changes with the fundamentalfrequency; its start at the low frequency with the 3rd orderdominating, then as the frequency increases first the 2nd or-der and then the 1st order dominates.Fig. 3. Variation in SNR. This figure shows the non-stationary behaviour of the signal to noise ratio when thefundamental frequency and harmonic orders changes asshown in figure 2. The noise variance is set for each ofthe three responses , such that an overall SNR of 10 dB isachieved.Fig. 4. Crame´r Rao Lower Bounds. The estimated standarddeviation, σ, of 100 Monte Carlo simulations of trackingfundamental frequency in the three different response sig-nals are plotted here using symbols. With lines are shownthe theoretical limits computed according to [4]. Thereis strong agreement between simulations and theory. It isnoted that tracking on the joined posterior improves the es-timate significantly.Fig. 5. Experimental Data. The spectrogram of the mea-sured vibrations in the x direction is shown here. In the twoother directions the spectrograms look similar (no shown).Fig. 6. Posterior distribtution of combined signals. For eachof the measured signals the posterior distribution for ωF iscomputed. Combining the three distribution results in thejoined posterior distribution shown here conditioned on allthree signals.Fig. 7. Tracking Result. The four tracks shown here corre-sponds to tracking on the posterior conditioned on each ofthe source signals. The fourth track is from the combinedposterior shown in figure 6. The deviations from the truetrack are summarized in table 1.

Bayesian multichannel tracking of periodic signals: A new way to detemine the running speed of mechanical systems

http://orbit.dtu.dk/ws/files/2789187/imm2888.pdf

Bayesian multichannel tracking of periodic signals: A new way to detemine the running speed of mechanical systems

Abstract

Similar works

Full text

Available Versions

Online Research Database In Technology