Locating
Faulty
Rolling Element Bearing
Signal by Simulated Annealing
Mid Year Report
AMSC 663,
2012
Jing Tian
Course Advisor:
Professor Radu Balan
, Professor Kayo Ide
Research
Advisor: Dr. Carlos Morillo
University of Maryland, College Park
Abstract
Vibration
acceleration
signal
is
widely
used in
the health
monitoring of rolling element
bearing
s
.
A critical work of the bearing fault diagnosis is
l
ocating the
optimum
frequency
band
that
contains faulty bearing signal
, which is usually masked by noise
.
This project diagnoses the
bearing fault with a combination of methods. The optimum frequency band is located by
maximizing spectral kurtosis with simulated annealing. Then the fau
lt feature frequency
component is extracted from the frequency band by using envelope analysis.
Modules of the
program
were
validated by a combination of analytic
work and simulation work. The
whole
program will be validated by
an open database
of
both nor
mal and faulty bearing vibration
signals.
I.
Introduction
Bearing provides relative rotational freedom and transmits a load between two structures.
It
is widely used in electromechanical systems.
Rolling element bearing
is a major source of failure
in electromechanical systems. For
example,
bearing faults account for more than 40% of th
e induction motor’s failure [1], and
g
earbox bearing failure is the top contributor of the wind turbines downtime
[2, 3]. Bearings are
inexpensive devices, but the fa
ilure of bearing is costly.
A $5,000 wind turbine bearing
replacement can easily turn into a $250,000 project, not to mention the cost of downtime [4]
.
In
1987, LOT Polish Airlines Flight 5055 Il

62M crashed because of failed bearings in one engine,
killin
g all
the183 people on the plane [5]. In

situ health monitoring is used to improve the
condition

based maintenance, which reduces the frequency and the loss of the bearing failure.
In the bearing health monitoring, e
arly detection of the bearing fault is a
major concern for
the industry
.
Vibration acceleration signal is widely
used in this purpose because it
is sensitive to
the bearing fault
and it can
be monitored in

situ
.
The objective of
the vibration signal
bearing fault detection is to test if the vibration signal
x(t)
contains the faulty bearing signal
s(t)
Faulty bearing:
x(t) = s(t) + ν(t)
Normal bearing:
x(t) = ν(t)
where
x(t)
is the monitored vibration signal;
s(t)
is the faulty bearing signal;
v(t)
is the noise,
which is unknown
.
An industrial practice to
test the existence of
s(t)
is to test if a unique frequency component
of
s(t)

the fault feature frequency component
can be extracted from
x(t)
or not.
If the fault
feature frequency component is ex
tracted, the hypothesis that the bearing is faulty is true,
otherwise the hypothesis is false.
According to the research in [6], f
aulty bearing signal
s(t)
is a modulated signal
s(t) = d(t)c(t)
where
d(t)
is the modulating signal.
It is a result of the periodic impact between the bearing’s
rolling elements and the fault on the bearing’s contact surface.
Its frequency co
mponent is the
fault feature frequency, which is illustrated in a simulated faulty bearing signal in Fig.1.
The
fr
equency is provi
ded by the bearing manufacturer;
c(t)
is the carrier signal, which is a result of
the loading and vibration transfer function. This signal is usually unknown.
Fig. 1, Faulty bearing signal
s(t)
f
Fault
is the fault feature frequency
Met
hods like envelope analysis have been developed to extract the fault feature frequency.
The problem is that i
n the presence of noise the extraction may fail.
The solution is to
band

pass
filter the vibratio
n signal in the frequency domain, as shown in
a si
mulated vibration signal in
Fig. 2.
Fig. 2, Vibration signal in the frequency domain
The challenge to design the filter is that the optimum frequency band to band

pass filter the
faulty bearing signal is usually unknown. This project provides a
solution to find the optimum
frequency band.
II.
Methodology
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
0.5
0
0.5
1
Time(s)
Amplitude
0
2000
4000
6000
8000
10000
12000
14000
0
2
4
x 10
4
Frequency (Hz)
Amplitude
1/
f
Fault
0
500
1000
1500
2000
2500
0
0.005
0.01
Frequency(Hz)
Amplitude
Faulty bearing
signal
This project
locates
the optimum frequency band by optimizing the band

pass filter with
simulated annealing
(SA).
The idea is,
t
he frequency band
dominated
by
the faulty bearing signal
is non

Gaussi
an, and
therefore it
has
a high
spectral kurtosis value
[7].
In the presence of white Gaussian noise, by
maximizing the SK
, the optimum frequency band for the faulty bearing signal
can be
found. The
optimization problem is to maximize SK in terms of the ce
ntral frequency, bandwidth, and the
order of the finite impulse response (FIR) band

pass filter.
where
f
c
is the frequency band’s central frequency; Δ
f
is the width of the band; M is the order of
FIR filter;
f
Faul
is the fault feature frequency;
f
s
is
the sampling rate.
When the optimum frequency band is obtained, envelope analysis is applied to the filtered
signal to extract the bearing faulty feature frequency.
Fig. 3 shows the flow chart of the algorithm.
Fig. 3, Flow chart of the algorithm
III.
Implementation and Validation
(1)
Filter

bank
x(n)
is the sampled version of the vibration signal
x(t)
. It has N points.
At first, the vibration
signal x(n) is band

pass filtered by a FIR filter
h
to produce the filtered signal
y(n)
:
2
2
;
2
)
,
,
(
f
f
f
f
f
f
f
to
Subject
M
f
f
SK
Maximize
s
c
s
Fault
c
FIR filter
h
i
(
f
ci
,
Δ
f
i
, M
i
)
SK
SA
Maximize SK by
f
c
,
Δ
f
, M
Optimized
FIR filter
h(
f
co
,
Δ
f
o
, M
o
)
EA
x(n)
y
i
(n)
SK
i
x(n)
y
o
(n)
FFT
a(n)
Magnitude
A(
f
)
A(
f
)
Maximized SK
SK
o
f=
f
Fault
?
The bearing
is faulty
The bearing
is normal
Yes
No
x(n)
is the sampled vibration signal;
y
i
(n) is filtered output of the
ith
FIR filter
h
i
;
SK
i
is the SK of the
y
i
(n);
y
o
(n) is the output of the optimized FIR filter;
a(n)
is the envelope of
y
o
(n) ;
A(
f
)
is the FFT of
a(n)
h
d
(n)
is the
impulse response of the filter
w(n) is the window function
. In this project, Hamming window will be used:
This algorithm is implemented by Matlab’s build

in function “fir1”.
Before optimizing the filter, initial input is obtained by calculating SK fo
r the signal filtered
by an FIR filter

bank.
The filter

bank has a
structure
of binary tree as
shown
in Fig.
4
.
S
k,j
denotes the
j
th filter at the
k
th level.
When the signal
is
processed by the filter

bank,
sub

signals
corresponding
to
the filters are
obta
ined.
Fig.
4
Structure of the FIR filter

bank
To validate the program, a test signal having four frequency components is used. The
sampling rate is 1024Hz, and 10240 data points are used.
The signal is shown in Fig.
5
.
)
3
/
510
2
cos(
4
)
2
/
384
2
cos(
6
)
4
/
200
2
cos(
8
)
2
2
cos(
10
t
t
t
t
x
h
n
x
n
y
)
(
)
(
)
(
)
(
n
w
n
h
h
d
)
2
/
(
]
2
/
2
/
)
2
/
sin[(
]
2
/
2
/
)
2
/
sin[(
)
(
M
n
f
f
f
M
n
f
f
f
M
n
n
h
s
c
s
c
d
M
n
M
n
n
w
0
),
2
cos(
46
.
0
54
.
0
)
(
Level 0
Level 1
Level 2
Level 3
…
Level k
S
k,j
…
0
f
s
/2
…
Frequency
S
0,1
…
…
S
2,3
S
2,2
S
2,1
S
1,2
S
1,1
…
S
3,1
S
3,2
S
3,3
S
3,4
S
3,5
S
3,6
S
3,7
S
3,8
S
2,4
Fig.
5
. Test Signal (Left: time domain; right: frequency domain; dotted red line: edge of a
frequency band)
A two

level filter

bank is used, which has 4 filters. For sampling rate = 1024Hz, the filters
are [0, 128]Hz, [128, 256]Hz, [256,384]Hz, [384, 512]
Hz
.
Frequency components are observed in the frequency bands
. Results in the frequency domain
are illustrated in Fig. 6
.
All the frequency components locate on the right place of the frequency
axis. But the magnitude was changed.
Fig.
6
.
Frequency components obtained by applying the filter

bank
(Above left: pass

filtered
in
[0, 128]Hz
; above right: pass

filtered in
[
128
,
256
]Hz
;
below left: pass

filtered in
[
256
,
384
]Hz
;
below right: pass

filtered in
[
384
,
512
]Hz
)
Response of the magnitu
de is obtained by calculating the magnitude of the Fourier transform
of the filters
’
coefficients, as illustrated in Fig.
7.
W
e can observe that the
magnitude response
at
384Hz
is smaller than one
. That
explains why the
384Hz frequency component has
smaller
magnitude than its real value.
Fig. 7.
Magnitude response of the filter

bank
(2)
Spectral Kurtosis
Then the spectral kurtosis of the filtered signal
y(n)
is calculated. Spectral kurtosis is d
efined
as follows:
where
κ
r
is the
r
th order cumulant.
Y(m)
is the DFT of the signal
y(n)
:
Both y(n) and Y(m) are N points sequences. SK is a real number.
To estimate SK, the formula for joint cumulant is used:
According to [8],
DFT of a stationary signal is a circular complex random
variable, and
E[Y(m)
2
]=0, E[Y
*
(m)
2
]=0.
Therefore, we have
Program of SK is validated by a white Gaussi
an noise and a cosine function
, which are
plotted in Fig. 8.
The
sampling rate is
fs
=
1024Hz, and
the number of the
data points
is N=
1024.
Fig. 8. White Gaussian noise(Above left: time domain; above right: frequency domain), and
cosine function (Below left: time domain; below right: frequency domain)
2
*
2
*
*
4
)}]
(
),
(
{
[
)}
(
),
(
),
(
),
(
{
m
Y
m
Y
m
Y
m
Y
m
Y
m
Y
SK
1
,...,
1
,
0
,
)
(
)
(
1
0
2
N
m
e
n
y
m
Y
N
n
N
n
m
i
)]
(
)
(
[
2
)]
(
)
(
[
)]
(
)
(
[
)]
(
)
(
)
(
)
(
[
*
*
*
*
*
m
Y
m
Y
E
m
Y
m
Y
E
m
Y
m
Y
E
m
Y
m
Y
m
Y
m
Y
E
SK
2
}]

)
(
{
[
}

)
(
{
}]

)
(
{
[
}]

)
(
{
[
2
}

)
(
{
2
2
4
2
2
2
2
4
m
Y
E
m
Y
E
m
Y
E
m
Y
E
m
Y
E
SK
If we investigate two signals
analytically,
we can get the following results.
Fourier transf
orm
of the white Gaussian noise also follows Gaussian
distribution
, and therefore its spectral kurtosis
is 0. For a cosine function
y=cos(
2
π
f
t)
,
f=20,
its Fourier transform is
)]
(
)
(
[
2
1
)
(
f
m
f
m
m
Y
where
δ
(m) is the delta function
.
W
hen
m=f
and
m=

f
,
Y(m)=1/2
, otherwise
Y(m)=0
.
SK of the cosine function is estimated with a given number of data points N=1024 as follows:
510
2
/
)
5
.
0
5
.
0
[(
/
)
5
.
0
5
.
0
(
2
]
/
)

)
(

[(
/
)

)
(

(
2
}]

)
(
{
[
}

)
(
{
2
2
2
2
4
4
2
1
2
1
4
2
2
4
N
N
N
m
Y
N
m
Y
m
Y
E
m
Y
E
SK
N
k
k
N
k
k
The numerical result is obtained as follows:
SK for the white noise is

0.0406
; SK for the
cosine function
is 510
.0000
.
Both
values
match the analytical results.
(3)
Simulated annealing
T
h
e process of estimating SK as a function of the FIR filter is optimized by simulated
annealing
(SA) [9]
, which is
a
metaheuristic global optimization tool.
The flowchart of
implementing is illustrated in Fig.
6
. In reach iteration, there is a chance that a worse case would
be accepted and thus simulated annealing can avoid the searching being trapped in a local
extremum.
The algorithm is validated by a one
dimension
function and a three dimension function.
The one dimension function tested the algorithm
’
s ability to locate the global optimum when the
function has many local extremum.
The function is plotted in Fig. 9. The global minim
um is y=

100 when x
=
π
.
2
)
(
)
cos(
100
x
x
y
Fig. 9. One dimension function
Setting of the simulated annealing is illustrated in Fig.
10
. Initial temperature T=1000.
Every time when a new result is accepted, T is reduced as T=0.99T. The initial input is a random
number between [

10, 10]. In every iteration a step is added to the input. The step is a random
10
0
10
20
100
0
100
200
x
y
number, the step size is reduced
as S=0.9999S when T drops.
The algorithm stops at 1
,
000
iterations.
Fig.
10
,
Setting of the simulated annealing for one dimension function
The minimum found was y=

99.999
7. The corresponding
minimize
r is x=
3.1440
.
The
searching process is illustrate
d in Fig. 11.
Fig. 11. Searching process (Left: variable; right: function value)
The algorithm was further validated by a three dimension function. This function can be
regarded as the addition of three functions
(Fig. 12)
. Global minimum is achieved
when all the
three functions
’
global minimum are reached. The global minimum is y=

150, when
x
1
= π, x
2
=15,
and x
3
=

4π.
)
4
/(
)
4
sin(
50
)
15
(
)
(
)
cos(
100
3
3
2
2
2
1
1
x
x
x
x
x
y
Initialize the temperature
T
End a round of searching
Use the initial input vector
W
Compute function value
SK(W)
Generate a
step
S
Compute function value
SK
(W+S)
SK(
W
+S)
<
SK(W)
exp[(
SK(W)

SK(
W
+S)
)/T] >
rand ?
Termination
criteria reached?
Replace
W
with
W+S
, reduce
T
Keep x unchanged, reduce
T
Yes
No
Yes
Yes
No
T = 1000
W: A random number in [

10,10]
S: A random number in
[

10,10]
T = 0.99T
T = 0.99T
1,000 iterations
Fig. 12. Three dimension function
Fig.
13.
describes the setting of the algorithm.
Initial temperature T=1000. Every time when a
new result is accepted, T is reduced as T=0.99T. The initial input is a three

element vector. Each
element is a random number between [

10, 10]. In every iteration a step vector is added to the
input. Elements
of the step vector are random numbers. Amplitude of the vector is reduced as
S=0.9999S when T drops. The algorithm stops at 100,000 iterations.
Fig. 13. Setting of the simulated annealing for three dimension function
The minimum found was y=

149.7608
. The corresponding
minimize
r is x=[
3.179
6,
14.808
8,

12.60
60].
The searching process is illustrated in Fig. 14.
20
10
0
10
20
30
100
0
100
200
x
y
2
1
1
)
(
)
cos(
100
x
x
y
2
2
)
15
(
x
y
)
4
/(
)
4
sin(
50
3
3
x
x
y
Initialize the temperature
T
End a round of searching
Use the initial input vector
W
Compute function value
SK(W)
Generate a random step
S
Compute function value
SK
(W+S)
SK(
W
+S)
<
SK(W)
exp[(
SK(W)

SK(
W
+S)
)/T] >
rand ?
Termination
criteria reached?
Replace
W
with
W+S
, reduce
T
Keep x unchanged, reduce
T
Yes
No
Yes
Yes
No
T = 1000
W: A vector, each element is a
number in [

10,10]
S: A vector, each element is a
number in [

10,10]
T = 0.99T
T = 0.99T
100,000 iterations
Fig. 14. Searching process (Left: variable; right: function value)
(4)
Envelope analysis
When the optimized frequency band is found,
envelope analysis is applied to the filtered
signal.
Th
e enveloped signal is obtained
from the m
agnitude of the analytic signal
which is
constructed
via Hilbert transform
:
Analytic signal
The envelope is the magnitude of the analytic signal
Fig.
15
shows the effect of envelope analysis on a modulated signal
Fig.
15
, Effect of envelope analysis
The program is validated by a modulated signal
in Fig. 16
:
2
1
)
1
(
x
x
y
where
)
10
2
cos(
1
t
x
)
4
/
200
2
cos(
100
2
t
x
d
t
h
y
t
y
o
o
)
(
)
(
)
(
ˆ
t
t
h
1
)
(
)
(
ˆ
)
(
)
(
t
y
j
t
y
t
y
o
o
a

)
(

)
(
t
y
t
a
a
0.032
0.033
0.034
0.035
0.036
0.037
0.038
0.5
0
0.5
Time(s)
Amplitude
Original signal
Hilbert transform of
the original signal
Enveloped signal
Fig.
16.
Modulated signal (Left: time domain; right: frequency domain)
Frequency component of the modulating signal is expected to be extracted after the envelope
analysis and Fourier transform.
After envelope analysis,
de

modulated signal is obtained
(Fig. 17)
.
The de

modulated signal
has the same frequency as the modulating signal, and the same amplitude as the original
modulated signal.
Fig.
17.
De

modulated signal (Left: time domain; right: frequency domain)
IV.
Summary
All the modules of the program achieved their designated function. Some
topics can be
investigated to improve the method.
Firstly, the reduced magnitude after filtering needs to be
compensated; secondly
, the performance of the envelope analysis module in c
ase
of noise has
not been examined.
In 2013 term, all the modules will be assembled to achieve the task of bearing fault detection.
And parallel computing will be implemented to improve the program.
2012
•
October

Literature review; exact validation methods
; code writing
•
November

Middle: code writing

End: Validation for envelope analysis and spectral kurtosis
•
December

Semester project report and presentation
2013
•
February

Complete validation
•
March

Adapt the code for parallel computing
•
April

Validate the
parallel version
•
May

Final report and presentation
References
[1] L. M. Popa, B.

B. Jensen, E. Ritchie, and I. Boldea, “Condition monitoring of wind
generators,” in
Proc. IAS Annu. Meeting
, vol. 3, 2003, pp. 1839

1846.
[2] Wind Stats Newsletter, 2003
–
200
9, vol. 16, no. 1 to vol. 22, no. 4, Haymarket Business
Media, London, UK
[3] H. Link; W. LaCava, J. van Dam, B. McNiff, S. Sheng, R. Wallen, M. McDade, S. Lambert,
S. Butterfield, and F. Oyague,“Gearbox Reliability Collaborative Project Report: Findings
from
Phase 1 and Phase 2 Testing", NREL Report No. TP

5000

51885, 2011
[4] C. Hatch, “Improved wind turbine condition monitoring using acceleration enveloping,”
Orbit, pp. 58

61, 2004.
[5] Plane crash information
http://www.planecrashinfo.com/1987/1987

26.htm
[6]
P. D. Mcfadden, and J. D. Smith, “Model for the vibration produced by a single
point defect in a rolling element bearing,” Journal of Sound and Vibration, vol. 96,
pp. 69

82, 1984.
[7] J. Antoni, “The spectral kurtosis: a useful tool for characterising non

stationary signals”,
Mechanical Systems and Signal Processing
,
20,
pp.282

307, 2006
[8] P. O. Amblard, M. Gaeta, J. L. Lacoume, “Statistics for complex variables and signals

Par
t I:
Variables”, Signal Processing
53
, pp. 1

13, 1996
[9] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, "Optimization by Simulated Annealing".
Science
220
(4598), pp. 671
–
680, 1983
Comments 0
Log in to post a comment