- Email: [email protected]
Adaptive Vold–Kalman filtering order tracking
ARTICLE IN PRESS Mechanical Systems and Signal Processing Mechanical Systems and Signal Processing 21 (2007) 2957–2969 www.elsevier.com/locate/jnlabr/ymssp
Adaptive Vold–Kalman filtering order tracking Min-Chun Pan, Cheng-Xue Wu Department of Mechanical Engineering, National Central University, Jhongli City, Taoyuan County 320, Taiwan, ROC Received 26 February 2007; received in revised form 7 June 2007; accepted 11 June 2007 Available online 14 June 2007
Abstract This study proposes and implements an adaptive Vold–Kalman filtering order tracking (VKF_OT) approach to overcome the drawbacks of the original VKF_OT scheme for condition monitoring and diagnosis of rotary machinery. The paper comprises theoretical derivation and numerical implementation. Comparisons of the adaptive VKF_OT scheme to the original are accomplished through processing two synthetic signals composed of close orders and crossing orders, respectively. Parameters such as the weighting factor and the correlation matrix of process noise, which influence tracking performance, are investigated. The adaptive scheme coping with end effects of computation can simultaneously extract multiple order/spectral components, and effectively decouple close and/or crossing orders associated with multi-axial reference rotating speeds. Furthermore, the adaptive OT scheme is realized through Kalman filtering based upon adapted one-step prediction scheme, where a parameter, the weighting factor, is newly introduced into the computation. Thus the technique can be computed on-line and implemented as a real-time processing application. r 2007 Elsevier Ltd. All rights reserved. Keywords: Adaptive order tracking; Vold–Kalman filtering; Acoustic and vibration signals; Rotary machine monitoring
1. Introduction Mechanical systems are the devices to transform various energy, e.g. hydraulic, thermal and electric energy into mechanical energy, or mechanical energy to other forms. Mechanical systems under periodic loading due to rotary operation usually respond in measurements with a superposition of sinusoids whose frequencies are integer (or fractional integer) multiples of the reference shaft speed, usually the revolution speed of the loading. The fundamental frequency corresponding to the reference speed is called a basic order. Most rotary machines are generally composed of various rotating elements, such as gear, chain, belt and shaft for power transmission. The frequencies of noise or vibration signals arising from these parts relate to the revolution speed of the shaft connected with the rotating parts. Order tracking (OT) is one of effective tools for the diagnosis of rotary machinery. This technique is to analyze and extract the sinusoidal contents of a dynamic signal acquired from an acousto-mechanical system under periodic loading. Conventional OT approaches such as the windowed Fourier transform (WFT) and the resampling methods are primarily based on Fourier analysis, and have limited resolution in some situations and suffer from Corresponding author. Tel./fax: +886 3 4267312.
E-mail address: [email protected] (M.-C. Pan). 0888-3270/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2007.06.002
ARTICLE IN PRESS 2958
M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
a number of shortcomings [1,2]. Vold and Leuridan proposed a novel algorithm, Vold–Kalman filtering order tracking (VKF_OT), for the estimation of a single order component, and further simultaneous estimation of multiple orders [3–6]. Pan and Lin [7,8] further explored the theoretical details of VKF_OT and compared the differences of two VKF_OT schemes: the angular-velocity and angular-displacement VKF_OTs. But it is noted that these two VKF_OT schemes must be calculated off-line and implemented as post-processing techniques. They are not real-time applications yet. The computations for both the angular-velocity and angular-displacement VKF_OTs are rather time-consuming since all the data of an acquired signal to be tracked for specific order/spectral components are included in huge matrix manipulation [9]. To cope with the computation complexity and consider for on-line applications, an adaptive real-time processing algorithm based on the Kalman filter is derived and realized in this paper. This implementation makes concurrent computation possible while the measured data streams are being acquired. 2. Theoretical basis The kth-order component arising from the operation of a rotary machine can be shown as [7] xk ðtÞ ¼ ak ðtÞyk ðtÞ þ ak ðtÞyk ðtÞ,
(1)
where ak(t) denotes the complex envelope, and ak(t) is the complex conjugate of ak(t) to make xk(t) a real waveform. It is noted that yk(t) is a carrier wave in the formality of [5] Z t yk ðtÞ ¼ exp ki oðuÞ du , (2) 0
where o(u) represents the speed of the reference axle, and the discrete form of Eq. (2) can be expressed as ! n X yk ðnÞ ¼ exp ki oðmÞDT .
Rt 0
oðuÞ du is the elapsed angular displacement. Thus
(3)
m¼0
The basic idea of the angular-displacement VKF_OT is to define local constraints, the so-called structural equation and data equation [5], which make the unknown complex envelopes smooth and relate the tracked orders to the measured signal. These two equations are solved by using a least-squares technique through a large number of matrix computations. It resembles the Wiener filter in mathematical operation. Furthermore, the Kalman filter is essentially a recursive algorithm whose solution converges to the optimum Wiener solution in some statistical sense. The mathematical model of the Kalman filter can be embodied in a pair of equations, i.e., the process and measurement equation [10], as below: Process equation : xðn þ 1Þ ¼ Fðn þ 1; nÞxðnÞ þ m1 ðnÞ,
(4)
Measurement equation : yðnÞ ¼ CðnÞxðnÞ þ m2 ðnÞ.
(5)
The proposed adaptive angular-displacement VKF_OT here provides us a recursive algorithm without having to solve the complex matrix problem and evaluate a huge inverse matrix with all observed time points. Thus the original (off-line) angular-displacement VKF_OT [5–7] with its computational complexity in solving for the inverse matrix can be improved. 2.1. Structural equation As the order component xk(t) to be tracked is formulated by Eq. (1), the envelope ak(t) needs to be evaluated for extracting the order components. Generally, ak(t) is a relatively smooth polynomial with a low degree due to the intrinsic inertia of a rotary mechanical system. It fulfills [5] d s ak ðtÞ ¼ ck ðtÞ, dts
(6)
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2959
where ck(t) denotes a higher-degree term in ak(t) and physically delivers the information of variation intensity in amplitude for order component k; and the index s is the differentiation order. Thus using a difference operator their discrete form can be expressed as rs ak ðnÞ ¼ ck ðnÞ,
(7)
where X denotes the difference operator. As we choose s ¼ 2 to alter Eq. (7) in a form of the process equation like Eq. (4), then Eq. (7) can be shown as ak ðn þ 1Þ 2ak ðnÞ þ ak ðn 1Þ ¼ ck ðnÞ. Rewriting and arranging Eq. (8) leads to matrix form like # " # # " " 0 ak ðnÞ 0 1 ak ðn 1Þ . ¼ þ ck ðnÞ ak ðn þ 1Þ ak ðnÞ 1 2
(8)
(9)
Then we symbolize Eq. (9) as ak ðn þ 1Þ ¼ Mak ðnÞ þ Ck ðnÞ.
(10)
To simultaneously track multiple order/spectral components, Eq. (10) needs to be augmented to all K-order/spectral components for data extraction. For multiple-OT and decoupling, Eq. (10) is further extended as 3 2 M 0 0 0 0 0 32 3 2 3 2 a1 ðn þ 1Þ a1 ðnÞ C1 ðnÞ 7 6 0 M 0 0 0 0 7 7 6 7 6 76 a ðnÞ 7 6 C2 ðnÞ 7 6 a2 ðn þ 1Þ 7 6 76 7 6 7 6 6 2 7 6 7 6 0 0 M 0 0 0 (11) 7¼6 6 .. . 7 þ 6 . 7. 76 7 6 . 6 6 .. 7 6 .. 7 . . . . . . . 7 5 4 .. 5 4 5 4 . . .. .. 54 .. .. .. CK ðnÞ aK ðn þ 1Þ aK ðnÞ 0 0 0 0 0 M Likewise, the symbolized form of Eq. (11) is Aðn þ 1Þ ¼ Fðn þ 1; nÞAðnÞ þ ZðnÞ,
(12)
where A(n) is a 2K 1 matrix (‘‘accumulated’’ K column vectors ak(n), which represents the kth order component at time instant n), Z(n) is a 2K 1 matrix (‘‘accumulated’’ K column vectors Ck(n), in fact), and F(n+1, n) is the transition matrix with a dimension 2K 2K. Eq. (12) is the structural equation of the adaptive angular-displacement VKF_OT, whose form is identical to the process equation (Eq. (4)) of the Kalman filter. The structural equation, A(n+1), is like an IIR lowpass filter, and characterizes the smoothness condition of unknown complex envelopes. A(n+1) can be regarded as the output of the lowpass filter excited by white noise Z(n).
2.2. Data equation A measured signal y(n) can be considered as a combination of several order/spectral components xk(n) expressed by Eq. (1), and measurement noise [5], i.e. X yðnÞ ¼ ak ðnÞyk ðnÞ þ xðnÞ, (13) k2j
where the integral j( ¼ 71, 72, 73, y, and/or 7K) denotes the order of spectral components to be tracked, and x(n) comprises unwanted spectral components and measurement errors, i.e., a difference of the measured signal from all tracked order components of interest. It is known y(n) is the measured signal and x(n) the unwanted part at time instant n. Let Bk(n) comprise a carrier signal with the form Bk ðnÞ ¼ 0 yk ðnÞ .
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2960
Therefore, Eq. (13) can be rewritten as h
yðnÞ B1 ðnÞ B2 ðnÞ B3 ðnÞ
2
a1 ðnÞ
3
7 6 i6 a2 ðnÞ 7 6 BK ðnÞ 6 . 7 7 ¼ xðnÞ 6 .. 7 5 4
(14)
aK ðnÞ
or symbolized as yðnÞ ¼ BðnÞAðnÞ þ xðnÞ,
h
where the carrier term BðnÞ ¼ B1 ðnÞ and 3 2 a1 ðnÞ 7 6 6 a2 ðnÞ 7 7 6 AðnÞ ¼ 6 . 7 6 .. 7 5 4 aK ðnÞ
B2 ðnÞ
B3 ðnÞ
(15) i BK ðnÞ is a known 1 2K measurement matrix,
is the complex envelope that modulates B(n). Eq. (15) is the data equation of the adaptive angulardisplacement VKF_OT, and has the same form as the measurement equation, Eq. (5), of a Kalman filter. 2.3. Kalman filter based on one-step prediction As derived before, we deduce Eqs. (12) and (15) of the adaptive angular-displacement VKF_OT through the structural and data equations (Eqs. (8) and (13)) of the original scheme that conducts order-tracking computation in an off-line manner. These two equations are formed identical to the process and measurement equation of Kalman filtering. Consequently, the order-tracking problem can be considered as a stateestimation problem like the procedure of Kalman filtering. In the process equation of a Kalman filter, the 2K 1 vector Z(n) represents process noise, modeled as a zero-mean, white-noise process whose correlation matrix is defined by ( Q1 ðnÞ; n ¼ m; H E½ZðnÞZ ðmÞ ¼ (16) 0; nam: Accordingly, we can define the correlation matrix of the process noise as ( Qk;l n ¼ m \ k ¼ l; 1 ðnÞ; n E½ck ðnÞcl ðmÞ ¼ 0; nam [ kal;
(17)
where ck(n) or c1(m) is the element in the process noise vector Z(n), k, lAj( ¼ 1,2,3,y, and/or K) represents the index of the order components to be tracked, Qk;l 1 ðnÞ is the element at the kth row and lth column of Q1(n), and the asterisk denotes complex conjugation. Likewise, the variable x(n) in the measurement equation of Kalman filtering is the measurement noise, modeled as a zero-mean, white-noise process whose correlation is defined by ( Q2 ðnÞ; n ¼ m; n E½xðnÞx ðmÞ ¼ (18) 0; nam: It is worth mentioning that the process noise ck(n) regards the higher order terms of the envelope (ak(n)) of each order component xk(n), and the measurement noise x(n) relates to the artifacts nothing to do with order components of interest. Thus the process noise and the measurement noise are statistically independent, so that E½ck ðnÞxn ðmÞ ¼ 0; for all n and m:
(19)
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2961
To design the parameters Q1(n) and Q2(n), we employ the weighting factor ever introduced in the angularvelocity VKF_OT [3] to our proposed Kalman filtering process that enables to carry out the tracking process in an adaptive way along a whole measurement. In the angular-velocity VKF_OT, the measure Se(n) and SZ(n) are defined as the standard deviations of the heterogeneity e(n) and nuisance component Z(n) of the structural and data equations, respectively. The ratio of the two standard deviation functions written as rðnÞ S Z ðnÞ=S ðnÞ is the weighting factor of the angular-velocity VKF_OT. Correspondingly, the introduced weighting factor r(n) can also satisfy Eq. (19) to be defined as Eq. (17) for all n6¼m. For a real scalar weighting factor corresponding to constant bandwidth [11], Q1(n) is a diagonal matrix taken to be a 2K 2K identity matrix multiplied by a designated variance value. The relation between Q1(n) and Q2(n) for the adaptive VKF_OT scheme can be acquired as S2Z ðnÞ S2 ðnÞ
¼ r2 ðnÞ ¼
Q2 ðnÞ Qk;l 1 ðnÞ
(20)
and then Q2 ðnÞ ¼ r2 ðnÞQk;l 1 ðnÞ.
(21)
As a result, the tracked order envelopes, A(n), can be recursively solved by using Kalman filtering based on one-step prediction. The following summarized computation steps for each measurement instant n (n ¼ 1, 2, 3, y, N; and N is the number of a measurement string) are adapted from [10] to accommodate the weighting factor, r(n), and parameters of the proposed adaptive OT here, i.e., for n ¼ 1, 2, 3, y, N, GðnÞ ¼ rðnÞFðn þ 1; nÞKðn; n 1ÞBH ðnÞ½BðnÞKðn; n 1ÞBH ðnÞ þ Q2 ðnÞ1 ,
(22)
^ aðnÞ ¼ yðnÞ BðnÞAðnjy n1 Þ,
(23)
^ þ 1jy Þ ¼ rðnÞFðn þ 1; nÞAðnjy ^ rðnÞAðn n n1 Þ þ GðnÞaðnÞ,
(24)
KðnÞ ¼ Kðn; n 1Þ r1 ðnÞFðn; n þ 1ÞGðnÞBðnÞKðn; n 1Þ,
(25)
r2 ðnÞKðn þ 1; nÞ ¼ rðnÞFðn þ 1; nÞKðnÞrðnÞFH ðn þ 1; nÞ þ r2 ðnÞQ1 ðnÞ.
(26)
Additionally, Table 1 presents a summary of the variables and parameters used to solve the Kalman-filter problem for the adaptive OT scheme; and Fig. 1 illustrates a block diagram of the Kalman-filtering process. Table 1 Summary of the Kalman variables and parameters used in the adaptive OT scheme (adapted from [10]) Variable
Definition
Dimension
A(n) y(n) F(n+1, n) B(n) Q1(n) Q2(n) ^ þ 1jy Þ Aðn n ^ Aðnjy nÞ G(n) a(n) K(n+1, n) K(n) r(n)
State at time n Observation at time n Transition matrix from time n to time n+1 Measurement matrix at time n Correlation matrix of process noise Z(n) Correlation matrix of measurement noise x(n) Predicted estimate of the state at time n+1 given the observations y(1), y(2), y, y(n) Filtered estimate of the state at time n given the observations y(1), y(2), y, y(n) Kalman gain at time n Innovations vector at time n ^ þ 1jy Þ Correlation matrix of the error in Aðn n ^ þ 1jy Þ Correlation matrix of the error in Aðn n Weighting factor at time n
2K 1K 11 2K 2K 1K 2K 2K 2K 11 2K 1K 2K 1K 2K 1K 11 2K 2K 2K 2K 11
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2962
Initial condition ˆ (1| y ) A 0
Observation vector y(n)
K(n, n-1)
Kalman gain computer
K(n, n-1)
One-step predictor
ˆ (n+1| y ) A n F(n, n+1)
Filtered state estimate ˆ (n | y ) A n
G(n)
Riccati equation solver
K(n+1, n)
initial condition K(1,0) z-1I Fig. 1. Block diagram of Kalman filtering based on one-step prediction [10] for proposed adaptive OT. Table 2 Order components composed in synthetic signal 1 Order components
1
4
4.2
9
Amplitude
Linearly increasing from 0 to 10
Linearly increasing from 3 to 13
Linearly increasing from 5 to 15
Linearly 10 (fixed)
3. Numerical implementation Two sets of synthetic signals are designated to validate the proposed adaptive OT scheme, and to compare with the original angular-displacement method as well. These synthetic signals comprise specific order components and resonances, where synthetic signal 1 is to justify the scheme for the effectiveness of distinguishing close orders, and synthetic signal 2 is to simulate multiple crossing orders arising from two different axle systems. Data sampling frequency, 2000 Hz, is employed for these data. Computation parameters that influence tracking performance, such as the weighting factor and the correlation matrix of process noise, are investigated here. 3.1. Two close-orders distinguishing—synthetic signal 1 Synthetic signal 1 comprises four order components including orders 1, 4, 4.2, and 9. The assumed reference shaft speed linearly increases from stationary to 3000 rpm in 5 s. Table 2 illustrates the amplitudes of these orders. This signal is designated to validate the proposed OT scheme in discriminating close orders. Fig. 2 displays the rpm-frequency spectrum evaluated by using the conventional WFT that characterizes two close orders, orders 4 and 4.2. Fig. 3 presents the computed amplitudes versus reference shaft speed tracked by using the original and adaptive scheme, respectively. Additionally, the waveforms of order 4 and 4.2 are extracted as shown in Fig. 4. Both techniques can discriminate two close order components effectively. It is, however, noted that the amplitudes computed by the angular-displacement scheme decline to
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2963
Rpm-frequency spectrum 1000 900
7
800 6
Frequency (Hz)
700 5 600 4
500 400
3
300 2 200 1
100 0
0
500
1000
1500
2000
2500
3000
Reference shaft revolution speed (rpm) Fig. 2. Illustration of rpm–frequency spectrum (synthetic signal 1).
zero at two ends of measurement duration due to the constraint of structural equation. More explanatorily, the reason can be given that the structural equation is locally represented by a low polynomial, but its matrix form for the computation of envelope is not completely defined at initial and final data points for all measured data length N. The weighting factor also affects the results of order components to be extracted. It is known that as the weighting factor is large, the tracking time response is long; and vice versa. That is, in the case of a large weighting factor, the estimated amplitudes are severely influenced by the structural equation [7,11]. For the use of the adaptive scheme, the computed amplitudes have slight indentation due to the weighting factor tuning even though the structural equation has constrained the estimates of the complex envelopes with smoothness conditions. As to computation end effects, no end effect at the latter portion of data although there still exist various initial time responses caused by the choice of weighing factors. To compare these two OT schemes, the original one yields smoother amplitudes than the adaptive one, but shows end effects at the initial and latter portion of measurement data. While applied in rotary machines, one significant feature for the adaptive scheme superior to the original is that the former is a recursive algorithm and can be calculated in a real-time manner. Moreover, computation efficiency can be a crucial issue for reallife applications. Table 3 shows the elapsed CPU time, and exhibits that the adaptive VKF_OT is around ten times efficient to the original in this case. The reason is that the time for evaluating the original scheme is mainly spent in the computation of inverse matrices using the conjugate gradient method. 3.2. Multi-axle crossing order decoupling—synthetic signal 2 Two sets of reference speed with order crossing are designated, where one shaft speed is the same as before, and the other shaft speeds up from stationary to 9000 rpm at 2 and 5 s. Accordingly, this synthetic signal comprises two sets of order components, one with orders 1, 4, and 9, and the other with orders 2 and 5 in accordance with respective reference speeds. The details of this synthetic signal are given in Table 4. The signal is designed to justify the performance of the proposed VKF_OT scheme in discriminating crossing orders arising from different rotating systems.
ARTICLE IN PRESS 2964
M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
Fig. 3. Amplitudes of orders 4 and 4.2 tracked by using (synthetic signal 1): (a) the angular-displacement VKF_OT, and (b) the adaptive scheme, respectively.
Fig. 5 shows the time–frequency spectrum of synthetic Signal 2 via the computation of the WFT, which illustrates two sets of orders and characterizes order-crossing occurrences. The computed amplitudes of orders 1, 4, and 9 using the original and adaptive VKF_OTs are illustrated in Fig. 6. Both techniques are capable of tracking specific orders crossing with other spectral components. The estimated amplitudes using the adaptive scheme show indentations although they have fewer end effects than the original scheme. The computed Order 1 shows unevenness at around 2 s due to the mixing of three run-up orders crossing concurrently. Fig. 7 displays the waveform estimations of specific orders extracted by the adaptive method. As to the computation efficiency for these two OT methods, Table 5 lists CPU time spent in extracting orders 9, 4 and 1, respectively. It is noted that computation consumption using the original scheme is 3.5 times as much as the adaptive scheme for tracking Order 9, which crosses once with other spectral components, as shown in Fig. 5. The CPU consumption dramatically increases up to around 17 times when extracting Orders 4 and 1, both of which cross twice with the other set of orders. In general, the computation time required for using VKF_OT methods increases as more observed time points and number of tracked order/spectral components increase are acquired. Another factor dominating computation efficiency is the number of crossing occurrences for an order to be tracked with other spectral components. While more crossing occurrences exist in an order, this increases computational complexity in solving inverse matrix. 3.3. Influence of the designed parameters To investigate the design parameters of adaptive angular-displacement VKF_OT on tracking performance, amplitude estimation of order 9 presented in the synthetic signals is conducted. The relation of main design parameters in the adaptive scheme formulates as shown in Eq. (21). As the weighting factor r and the correlation matrix of process noise Q1 are assigned, the correlation of measurement noise Q2 can be decided.
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2965
Fig. 4. Waveforms of orders 4 and 4.2 tracked by using (synthetic signal 1): (a) the angular-displacement VKF_OT, and (b) the adaptive scheme, respectively.
Table 3 Comparison of computation efficiency on synthetic signal 1 for the adaptive and original VKF_OT schemes Computation schemes
Angular-displacement VKF_OT
Adaptive angular-displacement VKF_OT
CPU time (s)
205.1
21.7
Table 4 Spectral components composed in synthetic signal 2 Set 1
Order numbers Amplitude
1 Linearly increasing from 0 to 10
4 Linearly increasing from 3 to 13
Set 2
Order numbers Amplitude
2 15 (fixed)
5 20 (fixed)
9 10 (fixed)
Fig. 8 illustrates tracked amplitudes using different weighting factors and a fixed Q1 ¼ 102 I2 2, where the designed amplitude of the synthesized order 9 is shown in the top. The results demonstrate that a large weighting factor may slow down convergence rate of the tracked order component, but obtain smooth steadystate value of the tracked amplitude. It signifies that the weighting-factor tuning will affect both the convergence rate and the smoothness of the complex envelops defined by the structural equation. Moreover, it reveals that a tradeoff exists between the rate of convergence and tracking accuracy. Fig. 9 illustrates a comparison of amplitude estimation with different Q1 and a fixed weighting factor r ¼ 1 106, where Q1 defined in the process equation of Kalman filter is designated to be a 2 2 identity matrix multiplied by different variance values q1 ¼ 102, 105, 108, and 1015, respectively. It is noted that as the variance value q1
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2966
Time-frequency spectrum 1000
9
900 8 800 7 700 Frequency (Hz)
6 600 5 500 4 400 3
300 200
2
100
1
0
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
Time (sec) Fig. 5. Spectrogram of two sets of synthesized order components (synthetic signal 2).
Fig. 6. Amplitudes of orders 1, 4, and 9 tracked by using (synthetic signal 2): (a) the angular-displacement VKF_OT, and (b) the adaptive scheme, respectively.
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2967
Fig. 7. Waveform of (synthetic signal 2): (a) order 9, (b) order 4, and (c) order 1 tracked by using the adaptive scheme, respectively. Table 5 Comparison of computation efficiency to synthetic signal 2 for using two VKF_OT schemes (CPU time in seconds) Computation schemes
Angular-displacement VKF_OT
Adaptive angular-displacement VKF_OT
72.6 473.8 511
20.5 30.4 32.3
Order components Order 9 Order 4 Order 1
is small, amplitude estimation converges fast to reach its designated value. Furthermore, it is observed, for the value of q1 ¼ 1 1015, the fast convergence to the amplitude, 10, accompanies acute initial transient. Observing both Figs. 8 and 9, it is worth mentioning that a large weighting factor r associated with a small variance matrix Q1 can yield both fast convergence and extract tracking performance, but a side-effect of acute initial transient. 4. Conclusions The study proposes, derives and implements an advanced VKF_OT approach, i.e. the adaptive angulardisplacement VKF_OT technique, to overcome the deficiencies of the original VKF_OT scheme. Comparisons and discussions of the developed scheme to the original are accomplished through processing two synthetic signals. The adapted scheme can simultaneously extract multiple order/spectral components, and effectively decouple close and/or crossing orders. As mentioned before, the adaptive scheme extracts
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2968
Fig. 8. Amplitudes estimation of order-9 computed by using the adaptive scheme with various weighting factor r.
Table 6 Comparisons on different aspects between three VKF_OT techniques Analysis schemes
Angular-velocity VKF_OT [3,4,7]
Angular-displacement VKF_OT [5–7]
Adaptive angulardisplacement VKF_OT
Fair No Good Fair (one order per computation) Not available
Poor No Good Good
Good Yes Fair Good
Good
Good
TF resolution tuning by weighting factor Not affected Necessary
TF resolution tuning by weighting factor
TF resolution tuning by weighting factor
Not affected Necessary
Not affected Necessary
Comparison items Computation efficiency Real time processing Waveform reconstruction Discrimination of two close orders Discrimination of crossing orders (spectral components) Time and frequency (TF) resolution Influences of high slew-rate Prior knowledge on system specification
order components with slight indentation, but its computation is extremely efficient. This merit benefits real-time processing, especially for the purpose of industrial applications. As the OT technique is not exact science, it may well be that the decrease in computation time can justify the reduced accuracy. As a summary, Table 6 compares characteristics of the adaptive VKF_OT scheme with other two techniques [7] in different aspects.
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2969
Fig. 9. Amplitude estimation of order 9 computed by using the adaptive scheme with various correlation matrix of process noise Q1, where q1 ¼ 102, 105, 108 and 1015, respectively.
Acknowledgements This work was supported by the National Science Council of Republic of China under the Grant NSC-952622-E-008-008-CC3. This funding is gratefully acknowledged. References [1] D.K. Bandhopadhyay, D. Griffiths, Methods for analyzing order spectra, SAE paper 951273, 1995. [2] M.Ch. Pan, J.X. Chen, Transmission noise identification using two-dimensional dynamic signal analysis, Journal of Sound and Vibration 262 (1) (2003) 117–140. [3] H. Vold, J. Leuridan, High resolution order tracking at extreme slew rates, using Kalman tracking filters, SAE Paper No. 931288, 1993. [4] H. Vold, J. Leuridan, High resolution order tracking using Kalman tracking filters-theory and applications, SAE Paper No. 951332, 1995. [5] H. Vold, M. Mains, J. Blough, Theoretical foundation for high performance order tracking with the Vold–Kalman tracking filter, SAE Paper No. 972007, 1997. [6] H. Vold, H. Herlufsen, Multi axle order tracking with the Vold–Kalman tracking filter, Sound and Vibration 31 (5) (1997) 30–34. [7] M.Ch. Pan, Y.F. Lin, Further exploration of Vold–Kalman filtering order tracking with shaft-speed information—(I) theoretical part, numerical implementation and parameter investigations, Mechanical Systems and Signal Processing 20 (5) (2006) 1134–1154. [8] M.Ch. Pan, Y.F. Lin, Further exploration of Vold–Kalman filtering order tracking with shaft-speed information—(II) engineering applications, Mechanical Systems and Signal Processing 20 (6) (2006) 1410–1428. [9] M.Ch. Pan, S.W. Liao, C.C. Chiu, Improvement on gabor order tracking and objective comparison with Vold–Kalman-filtering order tracking, Mechanical Systems and Signal Processing 21 (2) (2007) 653–667. [10] S. Haykin, Adaptive Filter Theory, Prentice-Hall, Upper Saddle River, NJ, 07458, 2002. [11] C. Feldbauer, R. Holdrich, Realization of a Vold–Kalman tracking filter—A least squares problem, Proceedings of the COST G-6 Conference on Digital Audio Effects, DAFX 1-4, 2000.
Adaptive Vold–Kalman filtering order tracking Min-Chun Pan, Cheng-Xue Wu Department of Mechanical Engineering, National Central University, Jhongli City, Taoyuan County 320, Taiwan, ROC Received 26 February 2007; received in revised form 7 June 2007; accepted 11 June 2007 Available online 14 June 2007
Abstract This study proposes and implements an adaptive Vold–Kalman filtering order tracking (VKF_OT) approach to overcome the drawbacks of the original VKF_OT scheme for condition monitoring and diagnosis of rotary machinery. The paper comprises theoretical derivation and numerical implementation. Comparisons of the adaptive VKF_OT scheme to the original are accomplished through processing two synthetic signals composed of close orders and crossing orders, respectively. Parameters such as the weighting factor and the correlation matrix of process noise, which influence tracking performance, are investigated. The adaptive scheme coping with end effects of computation can simultaneously extract multiple order/spectral components, and effectively decouple close and/or crossing orders associated with multi-axial reference rotating speeds. Furthermore, the adaptive OT scheme is realized through Kalman filtering based upon adapted one-step prediction scheme, where a parameter, the weighting factor, is newly introduced into the computation. Thus the technique can be computed on-line and implemented as a real-time processing application. r 2007 Elsevier Ltd. All rights reserved. Keywords: Adaptive order tracking; Vold–Kalman filtering; Acoustic and vibration signals; Rotary machine monitoring
1. Introduction Mechanical systems are the devices to transform various energy, e.g. hydraulic, thermal and electric energy into mechanical energy, or mechanical energy to other forms. Mechanical systems under periodic loading due to rotary operation usually respond in measurements with a superposition of sinusoids whose frequencies are integer (or fractional integer) multiples of the reference shaft speed, usually the revolution speed of the loading. The fundamental frequency corresponding to the reference speed is called a basic order. Most rotary machines are generally composed of various rotating elements, such as gear, chain, belt and shaft for power transmission. The frequencies of noise or vibration signals arising from these parts relate to the revolution speed of the shaft connected with the rotating parts. Order tracking (OT) is one of effective tools for the diagnosis of rotary machinery. This technique is to analyze and extract the sinusoidal contents of a dynamic signal acquired from an acousto-mechanical system under periodic loading. Conventional OT approaches such as the windowed Fourier transform (WFT) and the resampling methods are primarily based on Fourier analysis, and have limited resolution in some situations and suffer from Corresponding author. Tel./fax: +886 3 4267312.
E-mail address: [email protected] (M.-C. Pan). 0888-3270/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2007.06.002
ARTICLE IN PRESS 2958
M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
a number of shortcomings [1,2]. Vold and Leuridan proposed a novel algorithm, Vold–Kalman filtering order tracking (VKF_OT), for the estimation of a single order component, and further simultaneous estimation of multiple orders [3–6]. Pan and Lin [7,8] further explored the theoretical details of VKF_OT and compared the differences of two VKF_OT schemes: the angular-velocity and angular-displacement VKF_OTs. But it is noted that these two VKF_OT schemes must be calculated off-line and implemented as post-processing techniques. They are not real-time applications yet. The computations for both the angular-velocity and angular-displacement VKF_OTs are rather time-consuming since all the data of an acquired signal to be tracked for specific order/spectral components are included in huge matrix manipulation [9]. To cope with the computation complexity and consider for on-line applications, an adaptive real-time processing algorithm based on the Kalman filter is derived and realized in this paper. This implementation makes concurrent computation possible while the measured data streams are being acquired. 2. Theoretical basis The kth-order component arising from the operation of a rotary machine can be shown as [7] xk ðtÞ ¼ ak ðtÞyk ðtÞ þ ak ðtÞyk ðtÞ,
(1)
where ak(t) denotes the complex envelope, and ak(t) is the complex conjugate of ak(t) to make xk(t) a real waveform. It is noted that yk(t) is a carrier wave in the formality of [5] Z t yk ðtÞ ¼ exp ki oðuÞ du , (2) 0
where o(u) represents the speed of the reference axle, and the discrete form of Eq. (2) can be expressed as ! n X yk ðnÞ ¼ exp ki oðmÞDT .
Rt 0
oðuÞ du is the elapsed angular displacement. Thus
(3)
m¼0
The basic idea of the angular-displacement VKF_OT is to define local constraints, the so-called structural equation and data equation [5], which make the unknown complex envelopes smooth and relate the tracked orders to the measured signal. These two equations are solved by using a least-squares technique through a large number of matrix computations. It resembles the Wiener filter in mathematical operation. Furthermore, the Kalman filter is essentially a recursive algorithm whose solution converges to the optimum Wiener solution in some statistical sense. The mathematical model of the Kalman filter can be embodied in a pair of equations, i.e., the process and measurement equation [10], as below: Process equation : xðn þ 1Þ ¼ Fðn þ 1; nÞxðnÞ þ m1 ðnÞ,
(4)
Measurement equation : yðnÞ ¼ CðnÞxðnÞ þ m2 ðnÞ.
(5)
The proposed adaptive angular-displacement VKF_OT here provides us a recursive algorithm without having to solve the complex matrix problem and evaluate a huge inverse matrix with all observed time points. Thus the original (off-line) angular-displacement VKF_OT [5–7] with its computational complexity in solving for the inverse matrix can be improved. 2.1. Structural equation As the order component xk(t) to be tracked is formulated by Eq. (1), the envelope ak(t) needs to be evaluated for extracting the order components. Generally, ak(t) is a relatively smooth polynomial with a low degree due to the intrinsic inertia of a rotary mechanical system. It fulfills [5] d s ak ðtÞ ¼ ck ðtÞ, dts
(6)
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2959
where ck(t) denotes a higher-degree term in ak(t) and physically delivers the information of variation intensity in amplitude for order component k; and the index s is the differentiation order. Thus using a difference operator their discrete form can be expressed as rs ak ðnÞ ¼ ck ðnÞ,
(7)
where X denotes the difference operator. As we choose s ¼ 2 to alter Eq. (7) in a form of the process equation like Eq. (4), then Eq. (7) can be shown as ak ðn þ 1Þ 2ak ðnÞ þ ak ðn 1Þ ¼ ck ðnÞ. Rewriting and arranging Eq. (8) leads to matrix form like # " # # " " 0 ak ðnÞ 0 1 ak ðn 1Þ . ¼ þ ck ðnÞ ak ðn þ 1Þ ak ðnÞ 1 2
(8)
(9)
Then we symbolize Eq. (9) as ak ðn þ 1Þ ¼ Mak ðnÞ þ Ck ðnÞ.
(10)
To simultaneously track multiple order/spectral components, Eq. (10) needs to be augmented to all K-order/spectral components for data extraction. For multiple-OT and decoupling, Eq. (10) is further extended as 3 2 M 0 0 0 0 0 32 3 2 3 2 a1 ðn þ 1Þ a1 ðnÞ C1 ðnÞ 7 6 0 M 0 0 0 0 7 7 6 7 6 76 a ðnÞ 7 6 C2 ðnÞ 7 6 a2 ðn þ 1Þ 7 6 76 7 6 7 6 6 2 7 6 7 6 0 0 M 0 0 0 (11) 7¼6 6 .. . 7 þ 6 . 7. 76 7 6 . 6 6 .. 7 6 .. 7 . . . . . . . 7 5 4 .. 5 4 5 4 . . .. .. 54 .. .. .. CK ðnÞ aK ðn þ 1Þ aK ðnÞ 0 0 0 0 0 M Likewise, the symbolized form of Eq. (11) is Aðn þ 1Þ ¼ Fðn þ 1; nÞAðnÞ þ ZðnÞ,
(12)
where A(n) is a 2K 1 matrix (‘‘accumulated’’ K column vectors ak(n), which represents the kth order component at time instant n), Z(n) is a 2K 1 matrix (‘‘accumulated’’ K column vectors Ck(n), in fact), and F(n+1, n) is the transition matrix with a dimension 2K 2K. Eq. (12) is the structural equation of the adaptive angular-displacement VKF_OT, whose form is identical to the process equation (Eq. (4)) of the Kalman filter. The structural equation, A(n+1), is like an IIR lowpass filter, and characterizes the smoothness condition of unknown complex envelopes. A(n+1) can be regarded as the output of the lowpass filter excited by white noise Z(n).
2.2. Data equation A measured signal y(n) can be considered as a combination of several order/spectral components xk(n) expressed by Eq. (1), and measurement noise [5], i.e. X yðnÞ ¼ ak ðnÞyk ðnÞ þ xðnÞ, (13) k2j
where the integral j( ¼ 71, 72, 73, y, and/or 7K) denotes the order of spectral components to be tracked, and x(n) comprises unwanted spectral components and measurement errors, i.e., a difference of the measured signal from all tracked order components of interest. It is known y(n) is the measured signal and x(n) the unwanted part at time instant n. Let Bk(n) comprise a carrier signal with the form Bk ðnÞ ¼ 0 yk ðnÞ .
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2960
Therefore, Eq. (13) can be rewritten as h
yðnÞ B1 ðnÞ B2 ðnÞ B3 ðnÞ
2
a1 ðnÞ
3
7 6 i6 a2 ðnÞ 7 6 BK ðnÞ 6 . 7 7 ¼ xðnÞ 6 .. 7 5 4
(14)
aK ðnÞ
or symbolized as yðnÞ ¼ BðnÞAðnÞ þ xðnÞ,
h
where the carrier term BðnÞ ¼ B1 ðnÞ and 3 2 a1 ðnÞ 7 6 6 a2 ðnÞ 7 7 6 AðnÞ ¼ 6 . 7 6 .. 7 5 4 aK ðnÞ
B2 ðnÞ
B3 ðnÞ
(15) i BK ðnÞ is a known 1 2K measurement matrix,
is the complex envelope that modulates B(n). Eq. (15) is the data equation of the adaptive angulardisplacement VKF_OT, and has the same form as the measurement equation, Eq. (5), of a Kalman filter. 2.3. Kalman filter based on one-step prediction As derived before, we deduce Eqs. (12) and (15) of the adaptive angular-displacement VKF_OT through the structural and data equations (Eqs. (8) and (13)) of the original scheme that conducts order-tracking computation in an off-line manner. These two equations are formed identical to the process and measurement equation of Kalman filtering. Consequently, the order-tracking problem can be considered as a stateestimation problem like the procedure of Kalman filtering. In the process equation of a Kalman filter, the 2K 1 vector Z(n) represents process noise, modeled as a zero-mean, white-noise process whose correlation matrix is defined by ( Q1 ðnÞ; n ¼ m; H E½ZðnÞZ ðmÞ ¼ (16) 0; nam: Accordingly, we can define the correlation matrix of the process noise as ( Qk;l n ¼ m \ k ¼ l; 1 ðnÞ; n E½ck ðnÞcl ðmÞ ¼ 0; nam [ kal;
(17)
where ck(n) or c1(m) is the element in the process noise vector Z(n), k, lAj( ¼ 1,2,3,y, and/or K) represents the index of the order components to be tracked, Qk;l 1 ðnÞ is the element at the kth row and lth column of Q1(n), and the asterisk denotes complex conjugation. Likewise, the variable x(n) in the measurement equation of Kalman filtering is the measurement noise, modeled as a zero-mean, white-noise process whose correlation is defined by ( Q2 ðnÞ; n ¼ m; n E½xðnÞx ðmÞ ¼ (18) 0; nam: It is worth mentioning that the process noise ck(n) regards the higher order terms of the envelope (ak(n)) of each order component xk(n), and the measurement noise x(n) relates to the artifacts nothing to do with order components of interest. Thus the process noise and the measurement noise are statistically independent, so that E½ck ðnÞxn ðmÞ ¼ 0; for all n and m:
(19)
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2961
To design the parameters Q1(n) and Q2(n), we employ the weighting factor ever introduced in the angularvelocity VKF_OT [3] to our proposed Kalman filtering process that enables to carry out the tracking process in an adaptive way along a whole measurement. In the angular-velocity VKF_OT, the measure Se(n) and SZ(n) are defined as the standard deviations of the heterogeneity e(n) and nuisance component Z(n) of the structural and data equations, respectively. The ratio of the two standard deviation functions written as rðnÞ S Z ðnÞ=S ðnÞ is the weighting factor of the angular-velocity VKF_OT. Correspondingly, the introduced weighting factor r(n) can also satisfy Eq. (19) to be defined as Eq. (17) for all n6¼m. For a real scalar weighting factor corresponding to constant bandwidth [11], Q1(n) is a diagonal matrix taken to be a 2K 2K identity matrix multiplied by a designated variance value. The relation between Q1(n) and Q2(n) for the adaptive VKF_OT scheme can be acquired as S2Z ðnÞ S2 ðnÞ
¼ r2 ðnÞ ¼
Q2 ðnÞ Qk;l 1 ðnÞ
(20)
and then Q2 ðnÞ ¼ r2 ðnÞQk;l 1 ðnÞ.
(21)
As a result, the tracked order envelopes, A(n), can be recursively solved by using Kalman filtering based on one-step prediction. The following summarized computation steps for each measurement instant n (n ¼ 1, 2, 3, y, N; and N is the number of a measurement string) are adapted from [10] to accommodate the weighting factor, r(n), and parameters of the proposed adaptive OT here, i.e., for n ¼ 1, 2, 3, y, N, GðnÞ ¼ rðnÞFðn þ 1; nÞKðn; n 1ÞBH ðnÞ½BðnÞKðn; n 1ÞBH ðnÞ þ Q2 ðnÞ1 ,
(22)
^ aðnÞ ¼ yðnÞ BðnÞAðnjy n1 Þ,
(23)
^ þ 1jy Þ ¼ rðnÞFðn þ 1; nÞAðnjy ^ rðnÞAðn n n1 Þ þ GðnÞaðnÞ,
(24)
KðnÞ ¼ Kðn; n 1Þ r1 ðnÞFðn; n þ 1ÞGðnÞBðnÞKðn; n 1Þ,
(25)
r2 ðnÞKðn þ 1; nÞ ¼ rðnÞFðn þ 1; nÞKðnÞrðnÞFH ðn þ 1; nÞ þ r2 ðnÞQ1 ðnÞ.
(26)
Additionally, Table 1 presents a summary of the variables and parameters used to solve the Kalman-filter problem for the adaptive OT scheme; and Fig. 1 illustrates a block diagram of the Kalman-filtering process. Table 1 Summary of the Kalman variables and parameters used in the adaptive OT scheme (adapted from [10]) Variable
Definition
Dimension
A(n) y(n) F(n+1, n) B(n) Q1(n) Q2(n) ^ þ 1jy Þ Aðn n ^ Aðnjy nÞ G(n) a(n) K(n+1, n) K(n) r(n)
State at time n Observation at time n Transition matrix from time n to time n+1 Measurement matrix at time n Correlation matrix of process noise Z(n) Correlation matrix of measurement noise x(n) Predicted estimate of the state at time n+1 given the observations y(1), y(2), y, y(n) Filtered estimate of the state at time n given the observations y(1), y(2), y, y(n) Kalman gain at time n Innovations vector at time n ^ þ 1jy Þ Correlation matrix of the error in Aðn n ^ þ 1jy Þ Correlation matrix of the error in Aðn n Weighting factor at time n
2K 1K 11 2K 2K 1K 2K 2K 2K 11 2K 1K 2K 1K 2K 1K 11 2K 2K 2K 2K 11
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2962
Initial condition ˆ (1| y ) A 0
Observation vector y(n)
K(n, n-1)
Kalman gain computer
K(n, n-1)
One-step predictor
ˆ (n+1| y ) A n F(n, n+1)
Filtered state estimate ˆ (n | y ) A n
G(n)
Riccati equation solver
K(n+1, n)
initial condition K(1,0) z-1I Fig. 1. Block diagram of Kalman filtering based on one-step prediction [10] for proposed adaptive OT. Table 2 Order components composed in synthetic signal 1 Order components
1
4
4.2
9
Amplitude
Linearly increasing from 0 to 10
Linearly increasing from 3 to 13
Linearly increasing from 5 to 15
Linearly 10 (fixed)
3. Numerical implementation Two sets of synthetic signals are designated to validate the proposed adaptive OT scheme, and to compare with the original angular-displacement method as well. These synthetic signals comprise specific order components and resonances, where synthetic signal 1 is to justify the scheme for the effectiveness of distinguishing close orders, and synthetic signal 2 is to simulate multiple crossing orders arising from two different axle systems. Data sampling frequency, 2000 Hz, is employed for these data. Computation parameters that influence tracking performance, such as the weighting factor and the correlation matrix of process noise, are investigated here. 3.1. Two close-orders distinguishing—synthetic signal 1 Synthetic signal 1 comprises four order components including orders 1, 4, 4.2, and 9. The assumed reference shaft speed linearly increases from stationary to 3000 rpm in 5 s. Table 2 illustrates the amplitudes of these orders. This signal is designated to validate the proposed OT scheme in discriminating close orders. Fig. 2 displays the rpm-frequency spectrum evaluated by using the conventional WFT that characterizes two close orders, orders 4 and 4.2. Fig. 3 presents the computed amplitudes versus reference shaft speed tracked by using the original and adaptive scheme, respectively. Additionally, the waveforms of order 4 and 4.2 are extracted as shown in Fig. 4. Both techniques can discriminate two close order components effectively. It is, however, noted that the amplitudes computed by the angular-displacement scheme decline to
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2963
Rpm-frequency spectrum 1000 900
7
800 6
Frequency (Hz)
700 5 600 4
500 400
3
300 2 200 1
100 0
0
500
1000
1500
2000
2500
3000
Reference shaft revolution speed (rpm) Fig. 2. Illustration of rpm–frequency spectrum (synthetic signal 1).
zero at two ends of measurement duration due to the constraint of structural equation. More explanatorily, the reason can be given that the structural equation is locally represented by a low polynomial, but its matrix form for the computation of envelope is not completely defined at initial and final data points for all measured data length N. The weighting factor also affects the results of order components to be extracted. It is known that as the weighting factor is large, the tracking time response is long; and vice versa. That is, in the case of a large weighting factor, the estimated amplitudes are severely influenced by the structural equation [7,11]. For the use of the adaptive scheme, the computed amplitudes have slight indentation due to the weighting factor tuning even though the structural equation has constrained the estimates of the complex envelopes with smoothness conditions. As to computation end effects, no end effect at the latter portion of data although there still exist various initial time responses caused by the choice of weighing factors. To compare these two OT schemes, the original one yields smoother amplitudes than the adaptive one, but shows end effects at the initial and latter portion of measurement data. While applied in rotary machines, one significant feature for the adaptive scheme superior to the original is that the former is a recursive algorithm and can be calculated in a real-time manner. Moreover, computation efficiency can be a crucial issue for reallife applications. Table 3 shows the elapsed CPU time, and exhibits that the adaptive VKF_OT is around ten times efficient to the original in this case. The reason is that the time for evaluating the original scheme is mainly spent in the computation of inverse matrices using the conjugate gradient method. 3.2. Multi-axle crossing order decoupling—synthetic signal 2 Two sets of reference speed with order crossing are designated, where one shaft speed is the same as before, and the other shaft speeds up from stationary to 9000 rpm at 2 and 5 s. Accordingly, this synthetic signal comprises two sets of order components, one with orders 1, 4, and 9, and the other with orders 2 and 5 in accordance with respective reference speeds. The details of this synthetic signal are given in Table 4. The signal is designed to justify the performance of the proposed VKF_OT scheme in discriminating crossing orders arising from different rotating systems.
ARTICLE IN PRESS 2964
M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
Fig. 3. Amplitudes of orders 4 and 4.2 tracked by using (synthetic signal 1): (a) the angular-displacement VKF_OT, and (b) the adaptive scheme, respectively.
Fig. 5 shows the time–frequency spectrum of synthetic Signal 2 via the computation of the WFT, which illustrates two sets of orders and characterizes order-crossing occurrences. The computed amplitudes of orders 1, 4, and 9 using the original and adaptive VKF_OTs are illustrated in Fig. 6. Both techniques are capable of tracking specific orders crossing with other spectral components. The estimated amplitudes using the adaptive scheme show indentations although they have fewer end effects than the original scheme. The computed Order 1 shows unevenness at around 2 s due to the mixing of three run-up orders crossing concurrently. Fig. 7 displays the waveform estimations of specific orders extracted by the adaptive method. As to the computation efficiency for these two OT methods, Table 5 lists CPU time spent in extracting orders 9, 4 and 1, respectively. It is noted that computation consumption using the original scheme is 3.5 times as much as the adaptive scheme for tracking Order 9, which crosses once with other spectral components, as shown in Fig. 5. The CPU consumption dramatically increases up to around 17 times when extracting Orders 4 and 1, both of which cross twice with the other set of orders. In general, the computation time required for using VKF_OT methods increases as more observed time points and number of tracked order/spectral components increase are acquired. Another factor dominating computation efficiency is the number of crossing occurrences for an order to be tracked with other spectral components. While more crossing occurrences exist in an order, this increases computational complexity in solving inverse matrix. 3.3. Influence of the designed parameters To investigate the design parameters of adaptive angular-displacement VKF_OT on tracking performance, amplitude estimation of order 9 presented in the synthetic signals is conducted. The relation of main design parameters in the adaptive scheme formulates as shown in Eq. (21). As the weighting factor r and the correlation matrix of process noise Q1 are assigned, the correlation of measurement noise Q2 can be decided.
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2965
Fig. 4. Waveforms of orders 4 and 4.2 tracked by using (synthetic signal 1): (a) the angular-displacement VKF_OT, and (b) the adaptive scheme, respectively.
Table 3 Comparison of computation efficiency on synthetic signal 1 for the adaptive and original VKF_OT schemes Computation schemes
Angular-displacement VKF_OT
Adaptive angular-displacement VKF_OT
CPU time (s)
205.1
21.7
Table 4 Spectral components composed in synthetic signal 2 Set 1
Order numbers Amplitude
1 Linearly increasing from 0 to 10
4 Linearly increasing from 3 to 13
Set 2
Order numbers Amplitude
2 15 (fixed)
5 20 (fixed)
9 10 (fixed)
Fig. 8 illustrates tracked amplitudes using different weighting factors and a fixed Q1 ¼ 102 I2 2, where the designed amplitude of the synthesized order 9 is shown in the top. The results demonstrate that a large weighting factor may slow down convergence rate of the tracked order component, but obtain smooth steadystate value of the tracked amplitude. It signifies that the weighting-factor tuning will affect both the convergence rate and the smoothness of the complex envelops defined by the structural equation. Moreover, it reveals that a tradeoff exists between the rate of convergence and tracking accuracy. Fig. 9 illustrates a comparison of amplitude estimation with different Q1 and a fixed weighting factor r ¼ 1 106, where Q1 defined in the process equation of Kalman filter is designated to be a 2 2 identity matrix multiplied by different variance values q1 ¼ 102, 105, 108, and 1015, respectively. It is noted that as the variance value q1
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2966
Time-frequency spectrum 1000
9
900 8 800 7 700 Frequency (Hz)
6 600 5 500 4 400 3
300 200
2
100
1
0
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
Time (sec) Fig. 5. Spectrogram of two sets of synthesized order components (synthetic signal 2).
Fig. 6. Amplitudes of orders 1, 4, and 9 tracked by using (synthetic signal 2): (a) the angular-displacement VKF_OT, and (b) the adaptive scheme, respectively.
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2967
Fig. 7. Waveform of (synthetic signal 2): (a) order 9, (b) order 4, and (c) order 1 tracked by using the adaptive scheme, respectively. Table 5 Comparison of computation efficiency to synthetic signal 2 for using two VKF_OT schemes (CPU time in seconds) Computation schemes
Angular-displacement VKF_OT
Adaptive angular-displacement VKF_OT
72.6 473.8 511
20.5 30.4 32.3
Order components Order 9 Order 4 Order 1
is small, amplitude estimation converges fast to reach its designated value. Furthermore, it is observed, for the value of q1 ¼ 1 1015, the fast convergence to the amplitude, 10, accompanies acute initial transient. Observing both Figs. 8 and 9, it is worth mentioning that a large weighting factor r associated with a small variance matrix Q1 can yield both fast convergence and extract tracking performance, but a side-effect of acute initial transient. 4. Conclusions The study proposes, derives and implements an advanced VKF_OT approach, i.e. the adaptive angulardisplacement VKF_OT technique, to overcome the deficiencies of the original VKF_OT scheme. Comparisons and discussions of the developed scheme to the original are accomplished through processing two synthetic signals. The adapted scheme can simultaneously extract multiple order/spectral components, and effectively decouple close and/or crossing orders. As mentioned before, the adaptive scheme extracts
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2968
Fig. 8. Amplitudes estimation of order-9 computed by using the adaptive scheme with various weighting factor r.
Table 6 Comparisons on different aspects between three VKF_OT techniques Analysis schemes
Angular-velocity VKF_OT [3,4,7]
Angular-displacement VKF_OT [5–7]
Adaptive angulardisplacement VKF_OT
Fair No Good Fair (one order per computation) Not available
Poor No Good Good
Good Yes Fair Good
Good
Good
TF resolution tuning by weighting factor Not affected Necessary
TF resolution tuning by weighting factor
TF resolution tuning by weighting factor
Not affected Necessary
Not affected Necessary
Comparison items Computation efficiency Real time processing Waveform reconstruction Discrimination of two close orders Discrimination of crossing orders (spectral components) Time and frequency (TF) resolution Influences of high slew-rate Prior knowledge on system specification
order components with slight indentation, but its computation is extremely efficient. This merit benefits real-time processing, especially for the purpose of industrial applications. As the OT technique is not exact science, it may well be that the decrease in computation time can justify the reduced accuracy. As a summary, Table 6 compares characteristics of the adaptive VKF_OT scheme with other two techniques [7] in different aspects.
ARTICLE IN PRESS M.-C. Pan, C.-X. Wu / Mechanical Systems and Signal Processing 21 (2007) 2957–2969
2969
Fig. 9. Amplitude estimation of order 9 computed by using the adaptive scheme with various correlation matrix of process noise Q1, where q1 ¼ 102, 105, 108 and 1015, respectively.
Acknowledgements This work was supported by the National Science Council of Republic of China under the Grant NSC-952622-E-008-008-CC3. This funding is gratefully acknowledged. References [1] D.K. Bandhopadhyay, D. Griffiths, Methods for analyzing order spectra, SAE paper 951273, 1995. [2] M.Ch. Pan, J.X. Chen, Transmission noise identification using two-dimensional dynamic signal analysis, Journal of Sound and Vibration 262 (1) (2003) 117–140. [3] H. Vold, J. Leuridan, High resolution order tracking at extreme slew rates, using Kalman tracking filters, SAE Paper No. 931288, 1993. [4] H. Vold, J. Leuridan, High resolution order tracking using Kalman tracking filters-theory and applications, SAE Paper No. 951332, 1995. [5] H. Vold, M. Mains, J. Blough, Theoretical foundation for high performance order tracking with the Vold–Kalman tracking filter, SAE Paper No. 972007, 1997. [6] H. Vold, H. Herlufsen, Multi axle order tracking with the Vold–Kalman tracking filter, Sound and Vibration 31 (5) (1997) 30–34. [7] M.Ch. Pan, Y.F. Lin, Further exploration of Vold–Kalman filtering order tracking with shaft-speed information—(I) theoretical part, numerical implementation and parameter investigations, Mechanical Systems and Signal Processing 20 (5) (2006) 1134–1154. [8] M.Ch. Pan, Y.F. Lin, Further exploration of Vold–Kalman filtering order tracking with shaft-speed information—(II) engineering applications, Mechanical Systems and Signal Processing 20 (6) (2006) 1410–1428. [9] M.Ch. Pan, S.W. Liao, C.C. Chiu, Improvement on gabor order tracking and objective comparison with Vold–Kalman-filtering order tracking, Mechanical Systems and Signal Processing 21 (2) (2007) 653–667. [10] S. Haykin, Adaptive Filter Theory, Prentice-Hall, Upper Saddle River, NJ, 07458, 2002. [11] C. Feldbauer, R. Holdrich, Realization of a Vold–Kalman tracking filter—A least squares problem, Proceedings of the COST G-6 Conference on Digital Audio Effects, DAFX 1-4, 2000.