Mechanical Systems and Signal Processing (1998) 12(3), 415–426

MODELING OF LOW SHAFT SPEED BEARING

FAULTS FOR CONDITION MONITORING
Y.-F. W*  P. J. K
CRASys, Department of Systems Engineering, ANU, Canberra, ACT 0200, Australia

(Received February 1997, accepted after revisions November 1997)

A general model of faulty rolling element bearing vibration signals is established. The
efficacy of the envelope-autocorrelation technique for condition monitoring of such
bearings leads us to derive the envelope-autocorrelation function of the model in the case
of very low shaft speed. For application to the bearing condition monitoring, a simplified
model is given. Finally the comparison of simulated data using the model with in situ data
confirms the satisfactory performance of the model and the effectiveness of
envelope-autocorrelation as a fault detection technique for low shaft speed bearings.
7 1998 Academic Press Limited

1. INTRODUCTION
The condition monitoring of rotating machinery is important in terms of system
maintenance and process automation, especially for the condition monitoring of rolling
element bearings which are the most common components in industrial rotating
machinery. They are used in industries from agriculture to aerospace, in equipment as
diverse as descaler pinch rolls to the space shuttle main engine turbopumps. For this
reason, a variety of bearing fault detection techniques have been proposed; but as far as
we know very few articles have reported the theoretical modeling of bearing faults with
mathematical descriptions.
Braun [1] gave a theoretical model of the rolling element bearing with multipoint bearing
defects. In his model, the repetitious impulse responses induced by bearing defects and the
modulation by the bearing structural vibration were introduced.
McFadden and Smith [2] improved the model by taking into account the impulse series,
the modulation of the periodic signal caused by non-uniform load distribution and of the
vibration transmission of rolling element bearing, as well as the exponential decay of
vibration.
Wang and Harrap [3] improved the model further by considering the impulse series, the
modulation of the periodic signal from non-uniform load distribution and of the first
bearing vibration mode with mathematical descriptions. They introduce the envelope
autocorrelation analysis for both simulated signals and experimental data.
Both [2] and [3] consider the modulation of cage frequency caused by a non-uniform
load distribution. If the fault occurs in a rolling element, such as a ball or roller, then the
cage frequency is what needs to be considered. If the fault occurs elsewhere, then due to
the non-uniform load distribution, the period of the signal will be different [4, 5].
* On leave from the Department of Mechanical Engineering, Zhejiang University of Technology,
Hangzhou, 310014, P. R. China. Present address: Department of Aerospace and Mechanical Engineering,
University College, UNSW, Australian Defence Force Academy, Canberra, ACT 2600, Australia.

0888–3270/98/030415 + 12 $25.00/0/pg970149 7 1998 Academic Press Limited

416 .-.   . . 
This paper is motivated by the poor fit of previous signal models to an industrial case
examined by the authors: condition monitoring of descaler pinch rolls in a steel mill [6].
Several attempts at standard (and not so standard) analysis of accelerometer signals had
been unable to predict failure in the rollers—failure which occurred rapidly and, in some
cases, catastrophically.
The remainder of this paper is organised as follows. The next section presents our
modified model for rolling element bearing fault signals. Section 3 derives the
envelope-autocorrelation of this signal model, and also give a model simplification for
application to the bearing condition monitoring. Section 4 compares signals generated
using our proposed model with signals obtained in an industrial condition monitoring
environment. It is concluded that the signal model presented accurately mimics the
industrial example.

2. A MODEL FOR ROLLING ELEMENT BEARING FAULTS

Suppose x(t) is an original signal from a rolling element bearing with a single point fault,
it can be expressed as
x(t) = xf (t) · xq (t) · xbs (t) + xs (t) + n(t), (1)
where: xf (t) is the basic impulse series produced by the fault which impacts repetitiously
with another surface in the bearing (the repetition period is Tf = 1/ff and ff is the
characteristic frequency of the bearing fault); xq (t) is the modulation effect caused by
non-uniform load distribution of the bearings and the cyclic variation of transmission path
between the fault impact site and the transducer; and xbs (t) is the bearing induced
vibration. This modulating effect is caused partly by the bearing geometry, but more
significantly by the position of the vibration sensor relative to the fault location. It is this
element which, if not taken into account, significantly effects the efficacy of the signal
model in real-life situations where optimal sensor position is not possible or perhaps
undesirable. xs (t) is the machinery induced vibration. The signal xs (t) is an unwanted,
structured, and predictable damped harmonic signal which is also considered to be a source
of noise. n(t) is a Gaussian white noise sequence with variance sn2 . This is unpredictable
measurement noise present in any practical measurement system.

2.1.    

The basic impulse series due to the bearing fault is assumed to be [1–3]
a Nf

xf (t) = d0 s d(t − jTf ) 1 A0 s cos (2plff t + ffl ), (2a)

j=1 l=1

where A0 = d0 ff is the amplitude, ffl is the initial phase of lth harmonic, and Nf is the
number of harmonics induced by impulse series. The characteristic frequency of a bearing
fault can be the characteristic frequencies of an inner race fault (fbpfi ), an outer race fault
(fbpfo ), and a rolling element fault (fbsf ), depending on where the fault occurs.

2.2. - 

The modulation effect due to non-uniform loading is assumed to be a cosine pulse series
[2, 3]; therefore, one has
Nq

xq (t) = s Aqk cos (2pkfq t + fqk ), (2b)

k=1
        417
T 1
Periodic characteristics of bearings with faults under various loadings and transmission path
effects
Outer race Inner race Rolling element
Causes of periodicity faults faults faults
Stationary loading No effect fr fc
Shaft unbalance fr No effect fr − fc
Rolling element fc fr − fc No effect
diameter variations
Transmission path No effect fr fc
variations fbsf /2
fr , Shaft rotational frequency; fc , cage frequency; fbsf , rolling element fault of frequency.

which is a periodic function with the period of Tq (fq = 1/Tq ). The amplitude and initial
phase of the kth harmonic are Aqk and fqk , respectively, and Nq is the number of harmonics.
The periodicity fq depends on the elements of bearing fault as shown in Table 1 [4, 5].

2.3. - 

If the vibration induced by the bearing is assumed to be the superposition of all
vibrational modes, then
Nbs

xbs (t) = s Absj e−abs t cos (2pfbsj t + fbsj ),

j
(2c)
j=1

where fbsj , Absj , fbsj , and absj are the resonance frequency, amplitude, initial phase, and
damping factor of jth vibrational mode of the bearing, respectively, and Nbs is the modal
order of bearing vibration signal.

2.4. - 

If the vibration due to machinery other than the bearing in question has the same form
as equation (2c), then
Ns

xs (t) = s Asn e−as t cos (2pfsn t + fsn ),

n
(2d)
n=1

where fsn , Asn , fsn , and asn are the resonance frequency, amplitude, initial phase, and
damping factor of nth vibrational mode of the machinery system, respectively, and Ns is
the modal order of machinery.
Theoretically, Nq , Nf , Nbs , and Ns are infinite. The band-limited nature of most
measurements means that it is reasonable to assume finite values for these in practice.

3. THE ENVELOPE-AUTOCORRELATION OF THE PROPOSED MODEL

In the very low shaft speed case, the following assumptions are made. The modulation
frequencies due to non-uniform load distribution and fault frequency are much lower than
the resonance frequencies from the bearing and machinery:
min {fbsj , fsn }max {ff , fq }, (j = 1, 2, . . . , Nbs ; n = 1, 2, . . . , Ns ,) (3)
418 .-.   . . 
and the basic impulse series are dominant in the vibrational signal from the faulty bearing:

A0max {Absj , Aqk , Asn }. (j = 1, 2, . . . , Nbs ; k = 1, 2, . . . , Nq ; n = 1, 2, . . . , Ns ) (4)

Because of inequality (3), xs (t) in equation (1) can be filtered out easily.
Furthermore it is assumed that damping and initial phase may be neglected in the
derivation of the envelope-autocorrelation, because they only affect the magnitudes and
do not affect the frequency components.
The concepts of envelope and autocorrelation are defined as follows.
The signal envelope of a real-valued signal x(t) is =z(t)=, where

z(t),x(t) + ix̃(t), (5)

and z(t) is the analytic signal [7] associated with x(t), i is the complex operator, x̃(t) is
the Hilbert transform [8] of x(t)

g
a
1 x(t)
x̃(t) = dt. (6)
p −a
t −t

The envelope-autocorrelation of an envelope signal =z(t)= is defined by

Rzz (t),E[=z(t)= · =z(t + t)=], (7a)

if =z(t)= can be considered ergodic, then the envelope-autocorrelation can be estimated as

g
T
1
Rzz (t) = =z(t)= · =z(t + t)= dt, (7b)
T 0

where T is some finite time interval. For estimating convenience, here another function
Fzz (t) which is related to Rzz (t) is introduced:

g
T
1
Fzz (t), =z(t)=2 · =z(t + t)=2 dt. (8)
T 0

Under the assumptions stated above and using the above definitions, one has

Nbs ,Nq ,Nf 4

A0
x(t) 1 xbs (t)xq (t)xf (t) + n(t) = s Absj Aqk s Bp (t) + n(t), (9a)
4 j,k,l = 1 p=1

where

B1 (t) = B1j,k,l (t) = cos [2p(fbsj + kfq + lff )t],

B2 (t) = B2j,k,l (t) = cos [2p(fbsj − kfq − lff )t],

B3 (t) = B3j,k,l (t) = cos [2p(fbsj + kfq − lff )t],

B4 (t) = B4j,k,l (t) = cos [2p(fbsj − kfq + lff )t]. (9b)

        419
Therefore, the function Fzz (t) can be expressed as (see Appendix A)
Nbs ,Nq ,Nf
1
Fzz (t) = C + s (Dj,k )2{2 cos [2p(2lff )t] + 2 cos [2p(2kfq )t]
2 j,k,l = 1

+ cos [2p(2lff + 2kfq )t] + cos [2p(2lff − 2kfq )t]}

Nq Nf 7
+ s s (Ek,k 1 )2 s Hm 1 (t), (except k = k1 + l = l1 ) (10)
k,k 1 = 1 l,l 1 = 1 m1 = 0

where
H0,1,...,7 (t) = cos {2p[(l 2 l1 )ff 2 (k 2 k1 )fq ]t}, (11a)
and
C,C12 + 2snb
2
,
Nbs ,Nq ,Nf
(A0 )2
C1, s (Absj )2(Aqk )2,
4 j,k,l = 1

(A0 )2(Absj )2(Aqk )2

Dj,k, ,
8
bs N
(A )2
Ek,k 1, 0 Aqk Aqk1 s (Absj )2, (11b)
8 j=1

where C, C1 , Dj,k , Ek,k 1 are constants. In the rolling element bearings, usually it is true for
fq Q ff . By comparing Rzz (t) with Fzz (t), the main frequencies in the envelope-autocorre-
lation function Rzz (t) are the fault characteristic frequency ff and its harmonics
lff (l = 2, 3, . . . , Nf ). kfq (k = 1, 2, . . . , Nq ) are the side frequencies of main frequencies, but
their amplitudes are smaller than those of main frequencies. There are also some
components of small amplitude with frequencies of 12 (l 2 l1 )ff 2 12 (k + k1 )fq (here (l 2 l1 ) is
odd, l, l1 = 1, 2, . . . , Nf ; k, k1 = 1, 2, . . . , Nq , except l = l1 + k = k1 ) in Rzz (t).

3.1.        

If the bearings are running under the condition of a light load, and the clearances of
bearings are not sufficient, the periodicity of non-uniform load distribution imposes hardly
any modulation on the impulse series caused by the bearing fault. In this case, it is
supposed
xq (t) = 1, (12)
then the theoretical model of this special case is
x'(t) 1 xbs (t)xf (t) + n(t), (13)
the function F'zz (t) which is related to envelope-autocorrelation function will be as follows,
N ,N
1 bs f
F'zz (t) = C' + s (D' )2 cos [2p(2lff )t]
2 j,l = 1 j

Nbs ,Nf ,Nf

+ s (D'j )2{cos [2p(l + l1 )ff t] + cos [2p(l − l1 )ff t]}, (l $ l1 ) (14)

j,l,l 1 = 1
420 .-.   . . 
where

C' = (C'1 )2 + 2(s'nb )2,

N ,N
(A0 )2 bs f j 2
C'1 = s (A ) ,
2 j,l = 1 bs

(A0 )2(Absj )2
D'j = , (15)
2

where C', C'1 , and D'j are constants. In this special case, the main frequencies in the
envelope-autocorrelation function are the fault characteristic frequency ff and its
harmonics lff (l = 2, 3, . . . , Nf ). They are dominant in the spectrum. There are still some
small-amplitude components with the frequencies of 12 (l 2 l1 )ff (here l 2 l1 is odd,
l, l1 = 1, 2, . . . , Nf , l $ l1 ). Only if the bearings run under the heavy load and their
clearances are big enough, will the general model be effective.

Figure 1. Simulated signal of a bearing with an inner race fault (fsamp = 5000 Hz): (a) simulated signal generated
by model; (b) power spectrum of (a); (c) envelope of (a); (d) envelope-autocorrelation of (a); (e)
envelope-autocorrelation power spectrum of (a).
        421

Figure 2. Simulated signal of a bearing in good condition (fsamp = 5000 Hz): (a) simulated signal generated by
model; (b) power spectrum of (a); (c) envelope of (a); (d) envelope-autocorrelation of (a); (e)
envelope-autocorrelation power spectrum of (a).

4. COMPARISON OF MODEL WITH IN SITU DATA

4.1.     
The in situ vibration data are the signals from a bearing housing which supports the
descaler pinch rolls. The mean rotational speed of the pinch rolls is 36.06 rpm (0.601 Hz),
and the shaft speed varies randomly in a small range about the mean. In our case, the
clearance of bearing is small. The pinch rolls are used to move a 25-mm diameter bar, and
the bearings run under a light load. A simulated signal from bearing with an inner race
fault corresponding to the in situ condition is generated with the model [shown in Fig. 1(a)]
with damping neglected and the sampling frequency fsamp = 5000 Hz. The modal orders of
vibrational signals from the bearing, machinery system, and impulse series are Nbs = 5,
Ns = 5 and Nf = 7, respectively, with a Gaussian white noise sequence (standard deviation
sn = 0.0122). The amplitude of the basic impulse series is dominant. The fault frequency
and its harmonics are submerged completely in the white noise in power spectrum [shown
in Fig. 1(b)]. The envelope contains the white noise. However, the period Tf and the fault
characteristic frequency ff are obvious in the envelope-autocorrelation and its power
spectrum [shown in Fig. 1(d) and (e)].
422 .-.   . . 
4.2.     
In a simulation of a bearing in good condition, no impulse series are present in the signal
Fig. 2. Without other excitations, the vibrational signals are only of small amplitude. In
this case, the white noise sequence is large in comparison to the vibrational signal. Here
the simulated signal is generated with xf (t) = 1, with other parameters as those in Section
4.1.

4.3.   :  

The in situ vibration data are from a bearing housing of the descaler pinch rolls. The
rolls have a very low shaft speed with a small range of random variation about the mean.
The bearings are SKF23226, double-row, spherical roller bearings. The signal in Fig. 3 is
from a bearing with inner race fault. The mean shaft speed of rolls is 36.06 RPM
(0.601 Hz). In the power spectrum the fault frequency and its harmonics are completely
submerged in the background noise [Fig. 3(b)]. The peak frequencies f1 = 8.5 Hz,
f2 = 50 Hz and their harmonics are from the background noise. These components are also

Figure 3. In situ data from a bearing with an inner race fault (fsamp = 5000 Hz): (a) in situ data from bar 1365;
(b) power spectrum of (a); (c) envelope of (a); (d) envelope-autocorrelation of (a); (e) envelope-autocorrelation
power spectrum of (a).
        423

Figure 4. In situ data from a bearing in good condition (fsamp = 5000 Hz): (a) in situ data from bar 1370; (b)
power spectrum of (a); (c) envelope of (a); (d) envelope-autocorrelation of (a); (e) envelope-autocorrelation power
spectrum of (a).

present in the signals from healthy bearings [Fig. 4(b)]. From Fig. 3(d), there is an obvious
periodicity in the envelope-autocorrelation: its period is Tf = 0.154 s, hence ff = 1/
Tf = 6.5 Hz, which coincides with the characteristic frequency from Fig. 3(e). Comparing
simulated signal and in situ data from the bearing with an inner race fault, they have similar
frequency components in the envelope-autocorrelation power spectra [Figs 1(e) and 3(e)],
and their envelope-autocorrelations [Figs 1(d) and 3(d)] have the same basic periodicity.

4.4.   : - 

Comparing the simulated signal shown in Fig. 2 and in situ data shown in Fig. 4 from
the good condition bearing, each of the original signals [Figs 2(a) and 4(a)],
envelope-autocorrelation [Figs 2(d) and 4(d)] and envelope-autocorrelation power spectra
[Figs 2(e) and 4(e)], is very close, which proves the hypothesis. For the healthy bearing,
the other excitations are small, and the white noise sequence is relatively dominant. In the
power spectra of in situ data, there are some peak frequency components from background
noise [Figs 3(b) and 4(b)].
424 .-.   . . 
5. CONCLUSIONS
In this paper, a theoretical model of bearing fault is established. The function Fzz (t),
which is related to the envelope-autocorrelation of signal in the case of very low shaft
speed, is derived. According to the in situ condition, a special case of the model is given
for condition monitoring. By comparing the simulated signals and in situ data, the
following conclusions can be obtained.
, The simulated data from the model confirms the in situ data, and the model has a
satisfactory performance.
, It is proved theoretically that the main frequencies in the envelope-autocorrelation of
signals from faulty bearings running under a very low shaft speed are the fault
characteristic frequency ff and its harmonics lff (l = 1, 2, . . . . , Nf ) [equation (14)].
, Envelope-autocorrelation [Figs 1(c) and 3(c)] from faulty bearing in the very low shaft
speed case is periodical, and the period is Tf = 1/ff . However, for healthy bearings, there
is no periodicity in the envelope-autocorrelation [Figs 2(d) and 4(d)].
, From the envelope-autocorrelation power spectra, the fault frequency and its harmonics
are dominant, but in the power spectra, they are completely submerged in the back-
ground noise, and there are only some peak frequency components from the background
noise [Figs 1(b) and 4(b)].

ACKNOWLEDGEMENTS
The authors acknowledge the funding of the activities of the Cooperative Research
Centre for Robust and Adaptive Systems by the Australian Commonwealth Government
under the Cooperative Research Centres Programme.

REFERENCES
1. S. G. B 1980 IEEE Transactions on Sonics and Ultrasonics SU-27, 317–328. The signature
analysis of sonic bearing vibrations.
2. P. D. MF and J. D. S 1984 Journal of Sound and Vibration 96, 69–82. Model for
the vibration produced by a single point defect in a rolling element bearing.
3. W. Y. W and M. J. H 1996 Machine Vibration 5, 34–44. Condition monitoring of ball
bearings using envelope autocorrelation technique.
4. Y.-T. S, S.-J. L 1992 Journal of Sound and Vibration 155, 75–84. On initial fault detection
on a tapered roller bearing: frequency domain analysis.
5. I. M. H 1994 Noise and Vibration Conference, Perth. A review of factors affecting bearing
vibration monitoring.
6. J. M. S, P. J. S, A. J. C, F. B and L. B. W 1995
Proceedings of Control ’95, Melbourne, Australia, 117–121. Spectral estimation of variable speed
rotating machinery at a hot strip mill.
7. J. V 1948 Cables et Transmissions 2, 61–74. Theorie et applications de la notion de signal
analytique.
8. A. V. O and R. W. S 1989 Discrete-Time Signal Processing. Englewood Cliffs,
NJ: Prentice-Hall.

APPENDIX A: DERIVATION OF ENVELOPE-AUTOCORRELATION

FOR THE PROPOSED MODEL
For a given signal shown in equation (9), with equation (6), one has
Nbs ,Nq ,Nf 4
A0
x̃(t) 1 s Absj Aqk s B p (t) + ñ(t), (A1a)
4 j,k,l = 1 p=1
        425
where

B 1 (t) = B 1j,k,l (t) = sin [2p(fbsj + kfq + lff )t],

B 2 (t) = B 2j,k,l (t) = sin [2p(fbsj − kfq − lff )t],

B 3 (t) = B 3j,k,l (t) = sin [2p(fbsj + kfq − lff )t],

B 4 (t) = B 4j,k,l (t) = sin [2p(fbsj − kfq + lff )t], (A1b)

and ñ(t) is the Hilbert transform of n(t). Then the envelope =z(t)= of the original signal x(t)
can be calculated as follows

=z(t)=2 = x 2(t) + x̃ 2(t)

6
Nbs ,Nq ,Nf 4
(A0 )2
= s (Absj )2(Aqk )2 s [Bp2 (t) + B p2 (t)]
16 j,k,l = 1 p=1

7
3 4
+2 s s [Bp (t)Bp 1 (t) + B p (t)B p 1 (t)]
p = 1 p1 = p + 1

Nbs ,Nq ,Nf Nbs ,Nq ,Nf

(A0 )2
+ s s Absj Aqk Absj1 Aqk1
8 j,k,l = 1 j 1 ,k 1 ,l 1 = 1

6 7
4 4 4 4
× s Bp (t) s B'p (t) + s B p (t) s B p' (t)
p=1 p=1 p=1 p=1

Nbs ,Nq ,Nf 4 Nbs ,Nq ,Nf 4

A0 A0
+ n(t) s Absj Aqk s Bp (t) + ñ(t) s Absj Aqk s B p (t)
2 j,k,l = 1 p=1
2 j,k,l = 1 p=1

+ n 2(t) + ñ 2(t), (except j = j1 + k = k1 + l = l1 ) (A2)

where

B'p (t) = Bpj1 ,k1 ,l1 (t), (p = 1, 2, 3, 4) (A3a)

and

B p' (t) = B pj1 ,k1 ,l1 (t). (p = 1, 2, 3, 4) (A3b)

426 .-.   . . 
In the presence of a bearing fault, n(t) is very small compared with the signal
xbs (t)xq (t)xf (t), and n 2(t) + ñ 2(t) is rather smaller and can be ignored. Expanding and
simplifying equation (A2), one has
Nbs ,Nq ,Nf
(A0 )2
=z(t)=2 = s (Absj )2(Aqk )2{2 + 2 cos [2p(2lff )t] + 2 cos [2p(2kfq )t]
8 j,k,l = 1

+ cos [2p(2lff + 2kfq )t] + cos [2p(2lff − 2kfq )t]}

Nbs ,Nq ,Nf Nbs ,Nq ,Nf 15
(A0 )2
+ s s Absj Aqk Absj1 Aqk1 s H'm 2 (t)
8 j,k,l = 1 j 1 ,k 1 ,l 1 = 1 m2 = 0

Nbs ,Nq ,Nf 4 Nbs ,Nq ,Nf 4

A0 A0
+ n(t) s Absj Aqk s Bp (t) + ñ(t) s Absj Aqk s B p (t),
2 j,k,l = 1 p=1
2 j,k,l = 1 p=1

(except j = j1 + k = k1 + l = l1 ) (A4)
where
H'0,1,...,15 (t) = cos {2p[(fbsj − fbsj1 ) 2 (l 2 l1 )ff 2 (k 2 k1 )fq ]t}. (A5)
If j $ j1 , with inequality (3), the frequencies [(fbsj − fbsj1 ) 2 (l 2 l1 )ff 2 (k 2 k1 )fq ](l, l1 = 1,
2, . . . , Nf ; k, k1 = 1, 2, . . . , Nq ) are in high-frequency, so these components can be ignored
in the very low shaft speed case. In the second term of equation (A4), only the components
with j = j1 are remained, and one denotes
Nbs ,Nq ,Nf 4
A0
nb(t), n(t) s Absj Aqk s Bp (t),
2 j,k,l = 1 p=1

0
Nbs ,Nq ,Nf 4
A
s Absj Aqk s B p (t).
0
nb (t), ñ(t) (A6)
2 j,k,l = 1 p=1

If the autocorrelation is calculated for the envelope signal as shown in equation (8),
because of
E[nb(t)nb(t + t)] = snb
2
, (A7)
E[0
nb (t)0
nb (t + t)] = snb
2
, (A8)
where s 2
nb is the variance of the signals nb(t) and 0
nb (t), with equation (11b), the function
Fzz (t) can be expressed as
Nbs ,Nq ,Nf
1
Fzz (t) = C + s (Dj,k )2{2 cos [2p(2lff )t] + 2 cos [2p(2kfq )t]
2 j,k,l = 1

+ cos [2p(2lff + 2kfq )t] + cos [2p(2lff − 2kfq )t]}

Nq Nf 7
+ s s (Ek,k 1 )2 s Hm 1 (t), (except k = k1 + l = l1 ) (A9)
k,k 1 = 1 l,l 1 = 1 m1 = 0

where
H0,1,...,7 (t) = cos {2p[(l 2 l1 )ff 2 (k 2 k1 )fq ]t}, (A10)
and Fzz (t) is a function which is related to Rzz (t).