Transcript 
AUTOMATED TREND EXTRACTION OF SENSOR
SIGNALS FOR PATTERN BASED DATA
ANALYSIS
By
SRINIVAS S. GANTT
Bachelor of Engineering
Karnataka Regional Engineering College
Surathkal, India
1992
Submitted to the Faculty of the
Graduate College of the
Oklahoma State University
in partial fulfillment of
the requirements for
the Degree of
MASTER OF SCIENCE
December, 1996
AUTOMATED TREND EXTRACTION OF SENSOR
SIGNALS FOR PATTERN BASED DATA
ANALYSIS
Thesis Approved:
/{.
~eciP~
Dean of Graduate College
ii
ACKNOWLEDGMENTS
I'd like to take this opportunity to thank the School of Chemical Engineering,
Oklahoma State University for providing me with financial support.
I'd like to thank Dr. James R. Whiteley for being my advisor, and wish to
acknowledge his constant support and guidance throughout my stay at OSU. Also, a
large note of appreciation goes to Dr. Arland H. Johannes and Dr. Karen A. High for
serving on my thesis committee.
I'd like to thank Bruce Colgate, Jack DeVeaux, Julia A. Rumsey, and Loy at the
R&D center of the Phillips Petroleum Co., Bartlesville.
A personal note of gratitude goes to my friends in Surathkal, Visakhapatnam and
Bombay, most of all Sarin and Raghu. Without their support, I would never have
achieved all that I did. Another personal note of appreciation goes to my friends at OSU
wlth whom 1 spent many a sleepless night at school, especially during the first two
semesters. They made life at grad school real fun despite its trials and tribulations.
The deepest appreciation and thanks must go to my mother and brothers, Krishna
and lagan. for all that they had done for me. I'd like to thank my sister Bindu for her love
and affection. A special note ofthanks to Padma for bringing the sunshine to my life.
r wish to dedicate this work to the fond memory of my father.
]JJ
TABLE OF CONTENTS
CHAPTER PAGE
I. INTRODUCTION.......................................................... I
Thesis outline........................................................... 6
2. CONVENTIONAL SIGNAL PROCESSING TECHNIQUES.. ..... 8
Digital signal representation and sampling......................... 8
Motivation for signal filtering........................................ I I
Conventional signal processing techniques....... ... II
The direct methods............................................ 12
Box car method. . .. .. . ..... . . . .. . .. .... .. . .. .. .. 12
Backward slope method...................... 13
Combination box car and backward slope
method...................... 19
Digital filters.................................................... 23
Single exponential filter... .. .. .. 23
Double exponential filter... .. .. 26
Moving average filter.. .. .. .. .. ... 27
Rate of change filter.. .. .. .. .. .. .. .. 29
Performance of wavelets............................................... 31
Summary 33
3. THE FOURIER AND WAVELET TRANSFORMS.................. 35
Digital signal analysis and processing............................... 35
Signal Decomposition................................................. 37
The Fourier transform................................................. 37
Continuoustime Fourier series............................. 38
Discretetime Fourier series.. . .. . .. . . . .. .. .... . .. . .. . . ... .. . 39
Timefrequency localization.......................................... 47
The wavelet transform. . . ... ... ... .. . .. . ... .. .. . . . . .. .. .. .. . . .. . . . . .. .. 51
Discrete wavelet transform................................... 51
Scaling and wavelet functions....... .. 52
The Daubechies family... .. .. . .. . 53
Construction of first order Daubechies
wavelet................................................ 55
iv
Signal filtering using scaling and wavelet functions........... 57
Multi  resolution analysis..... ....... .. .. . ... .... .... .. . .. .. .... .. 62
4. AUTOMATED TREND EXTRACTION.......... 69
Need for automated trend extraction. .. ... . .. .. . .. . . . . . . ... .... . ... 69
Determination of the degree of smoothing from the
user input window length....... .. 70
Determination of the minimum degree of smoothing
from raw signal characteristics......... .. .. .. ... .. ..... ..... .... .... 73
Dominant frequency identification........................ 74
The cumulative power spectrum........................... 75
Determination of the minimum number of levels of
decomposition from the cutoff frequency index........ 78
Determination of the minimum window length
required for process monitoring. .. .. . . .. . .. .. . .. . 81
Automated Trend Extraction algorithm. . . .. .. .. .. 83
Chapter Summary 88
5. CASE STUDIES........... 89
Introduction............................................................ 90
Individual Case Studies...... 91
Case 1 91
Case 2 91
Case 3 91
Case 3a 95
Case 4.... 95
Case 4a 99
Case 4b 101
Case 4c ,. 10·\
Discussion 106
Case 5 107
Chapter Summary , . . .. . . . ... . .. .. ... . . ... . .. . .. . 112
6. CONCLUSIONS.......................................................... 114
Recommendations and future work................................. 1\5
v
LIST OF TABLES
Table Page
1. Number of levels of wavelet decomposition for various
cutoff frequency indices..................................... 80
2. Summary of the results obtained by the ATE code............. 90
vi
Figure
LIST OF FIGURES
Page
1. I. Typical sensor signal and its smoothed approximation. . .... .... 2
1.2. Sampled representation of a signal... .. 4
1.3. Sensor signal and its wavelet smoothed approximations
at different levels '" .. . . .. . . .. ... 5
2.1. Periodic sampling of an analog signal... .. .. 9
2.2. Boxcar algorithm. . . ... . ..... .. .. ... .. . .. . .. .. . ... ..... . .. . .. .. .. .. . ... 14
2.3. Performance of the Boxcar algorithm for a sensor signal........ . 15
2.4. Backward slope algorithm.................. .. 16
2.5. Performance of the Backward slope algorithm for a
sensor signal. . . .. . .. . . . . .. .. ...... .. . .. . . . .. . ... .. . . . . .. .. .. .. .. .. . .. .. .. 18
2.6. Combination boxcar and backward slope algorithm... ... ... .. .. . . 21
2.7. Performance of the Combination boxcar and backward
slope algorithm for a sensor signal........... .. .... ........ .... .. ..... 22
2.8. Performance of the Single Exponential filter........................ 25
2.9. Performance of the Double Exponential filter ..... '" .. ... . . . .. . . ... 28
2.10. Performance of the Moving average filter............................ 30
2.) I. Performance of the Rate of change filter.. . .. . . . . .. .. .... .... .. . . . . .. 32
2.12. Performance of sixth order Daubechies wavelets.................... 34
3.1. Example problem on Fourier coefficients............ .. 42
3.2. Cosine and sine coefficients of a Fourier decomposed signal..... 46
VII
3.3. Timefrequency localization................. 49
3.4a. Box function.......................................................... 55
3.4b. Haar wavelet................. 57
3.5. Representation of basic decomposition and reconstruction
of a signal 59
3.6. Representation of multiresolution algorithm... .. 63
3.7. Detail and blurred coefficients of a signal......... ... ..... ... .... 64
3.8. Smoothed signal at various levels of decomposition............ 66
4.1. Span of the scaling functions in the process monitoring
window 7\
4.2. Cumulative power spect.rum of a signal............................ 77
4.3. Multi resolution tree analysis................. 78
4.4. Representation of the dominant frequencies.. . ... .... . . . .. . ..... . 82
5. J. Case study I 92
5.2. Case study 2 93
5.3. Case study 3 94
5.4. Case study 3<1 .. ... .. . . .. 96
5.5. Case study 4 98
5.6. Case study 4a lOO
5.7. Case study 4b ]02
5.8. Case study 4c 104
5.9. Case study 5 108
5.\ O. Cumulati ve power spect.ra for different data windows.. . .. . . ... 110
VIII
b(t)
respectively, and on
(d) Higher value of recording limit, 0.05. High frequency noise is retained.
25
(a) Low value of alpha, 0.01. Smoothing desirable, fundamental
trend of the original signal is not preserved.
(c) Low value of alpha, 0.01. Fundamental trend of original signal
preserved with slight delay.
(b) Relatively higher value of alpha, 0.05. Ability to catch up with
abrupt changes poor.
Figure 2.8. Performance of the Single Exponential filter for typical sensor signals.
0.60
0.50
0.55
0.55
0.50
0.60
0.50
0.60
0.50
0.70
0.60
0.70
26
(14)
( 13)
(12)
(11 )
y;' = (y)(a)xlI + y(l a) YIJI +(I  y) Y;'_I
Y' = a 2 /I"x + 2( I  a)yI'II  (I  a)2 Y1''12
Writing the filter equation (II) for the previous sampling instant yields
The advantage of the double exponential filter over the single exponential filter is that it
Solving for YnJ from ( 13), substituting in (12), rearranging, and selecting y = a gives
The Matlab code for this filter is dubexpfil.m and is included in the Appendix.
provides better filtering of high frequency noise, especially if y =ex.
removing high frequency noise. This is a second order filter which is essentially
from the exponential filter in (10). The second filter can be expressed a'i
equivalent to two first order filters in series where the second filter treats the output signal
because of the averaging effect of a. The of performance of this filter is not acceptable
2.3.2.2 Double exponential filter
The double exponential filter is an improvement over the single exponential filter for
the trend in the original signal are recorded at a much later time in the smoothed signal
for realtime patternbased process monitoring.
the algorithm to record abrupt changes in the smoothed signal is limited. The changes in
retains some high frequency noise and is not smooth enough. Further, the capability of
the trend of the smoothed signal is closer to the trend. However, the smoothed signal till
27
Figures 2.9a, 2.9b, 2.9c and 2.9d show two signals with different characteristics
and their filtered signals, for different values of a. Depending on the value of a u ed, the
filtered signal may be smooth at the expense of preserving the fundamental trend of the
original signal, or it may retain noise. For low values of a the smoothing obtained may
be desirable, but the initial response of the filter is in the opposite direction and the
fundamental trend of the filtered signal is totally distorted. This could be explained by
equation (14). The processed value at any instant of time is dependent on the weighted
sum of the value of the original signal at that time and the processed values at the
previous two instants of time. This inverse response performance could lead to a wrong
diagnosis of the plant situation if the filtered signal were used for process monitoring
purposes. For a relatively higher value of a, the filtered signal retains some high
frequency noise and is unacceptable for patternbased process monitoring purposes.
2.3.2.3 Moving average fitter
This filter averages a specified number of past data points, by giving equal weight to each
data point. The movingaverage fHter is usually less effective than the exponential filter,
which gives more weight to the most recent data.
The movingaverage filter can be expressed as
1 n
y = Lx
• 1/ J i=nJ+J I
( 15)
where J is the number of past data points that are being averaged. The previous fiitered
value YflJ can be expressed as
28
(a) Low value of alpha, 0.01. Smoothing desirable but fundamental
trend of original signal totally distorted.
Figure 2.9. Performance of the Double Exponential filter for various values of alpha.
0.50
(d) Low value of alpha, 0.01. Fundamental trend of filtered signal
preserved but delayed in time.
0.55
(c) Higher value of alpha, 0.25. High frequency noise retained.
0.60
(b) Relatiively higher value of alpha, 0.1. Initial trend of filtered signal in
0.70 opposite direction and ability to catch up with abrupt changes is poor.
0.50
0.60
0.50
0.70
0.40
0.80
0.60
signal. These regions are indicated by the dotted circle in Figure 2.10. For higher values
low values of J, the high frequency noise of the original signal is retained. The response
above reasons make the filtered signals obtained using this filter unacceptable for pattern
29
( 17)
(16)
InI
y = ~ x
. III J L.i ;
;=nJ
I
)'" = )'nI +j(X" XnJ)
The Matlab code for this filter is mov_ave...fil.m, a copy of which is included in
If a noisy measurement changes suddenly by a large amount and then returns to the
preserving the fundamental trend of the original signal, or it may retain noise. Further,
original value at the next sampling instant, or close to it, a noise spike is said to have
2.3.2.4 Rate of change (Noise spike) filter
occurred. This filter is designed to eliminate such noise spikes. These filters are used to
the first J values of the filtered signal are the same as that of the original signal. The
Depending on the value of J used, the filtered signal may be smooth at the expense of
based monitoring purposes.
the Appendix. The performance of this filter is shown for various values of J in Figure
signals with different characteristics, and their filtered signals for different values of 1.
2.10. The first J data points of the original signal are initialized as filtered values. For
of J, the fundamental trend gets distorted. Figures 2.1 Oa, 2. lOb and 2. IOc show two
of this filter is too slow to accurately capture abrupt changes in the trend of the original
Subtracting (16) from (IS) gives the recursive form of the movingaverage filter:
0.50 ..L _
30
I
/
/
I
Fi gure 2.10. Performance of the Moving Average filter.
(a) Past data points used for averaging, J=25. Filtered signal not very smooth.
,
\
(c) Past data points used for averaging, J=100. High frequency noise retained.
0.50 .1.. _
(b) Past data points used for averaging, J =35. Filtered signal smooth,
fundamental trend of original signal not preserved.
0.60
" 0.70 /
"
I
\
I
0.60
I
0.50
0.60
0.70


31
limit how much the filtered output is pennitted to change from one sampling instant to
the next. If L1x is the maximum allowable change, the noise spike filter can be written as
{ x.
if Ix"  YII_I! ~ ~.x
Y/, = Ynl  tu: if Y"l  xlI)tu: (18)
YII_I +tu: if YII
_ 1  XII (tu:
If a large change in the measurement occurs, the filter replaces the measurement by the
previous filter output plus (or minus) the maximum aIJowable change.
The Matlab code for this filter is rate_of_chJil.m, a copy of which is included in
the Appendix.
The performance of this filter for two signals with different characteristics is
shown in Figure 2.11. The first signal and its overlapped filtered output are shown in
Figures 2.lla and 2.11 b for ~x equal to 0.0005 and 0.00075 respectively. The second
signal and its overlapped filtered output are shown in Figures 2.1 Ic and 2.1 Id for ~x
equal to 0.00009 and 0.002 respectively. Depending on the value of ~x used, the filtered
signal may be smooth at the expense of preserving the fundamental trend of the original
signal, or it may retain noise. A smoothed signal preserving the fundamental trend of the
original signal is not obtained using this filter. Consequently, this filter is not acceptable
for smoothing signals for patternbased data analysis.
2.4 Performance using wavelets
The basic problem with conventional methods described previously is the inability to
reject high frequency noise while simultaneously retaining the ability to capture abrupt
32
0.70
0.60
0.50 L _
(a) Low value of Delta x =0.0005. Fundamental trend of original signal
0.70 not preserved.
0.60
/
Low value of Delta x =0.00009. Fundamental trend of original signal
not preserved.
0.50
(b) Delta x =0.00075, ability to catch up with abrupt changes poor.
0.60
0.55
0.50
(c)
0.60
0.55
0.50
(d) Delta x =0.002, high frequency noise retained.
Figure 2.11. Performance of the Rate of Change filter.

33
changes, which are genuine. Wavelets can achieve this. For demonstration purposes
Figure 2.12 shows the same signals considered in the above techniques smoothed using
wavelets. The technical discussion on wavelets is presented in the foHowing chapter.
The intent here is to simply show the reader that a superior technique exists and to offer
visual proof.
Summary of this chapter
This chapter presented an overview of the conventional methods of signal processing and
explains why they are not suitable for accurate patternbased process monitoring.
Moreover, the conventional techniques involve the use of parameters that need to be
hardcoded by the user. Further experimentation is needed to find the optimum values of
these parameters which give as close to the desired performance as possible.
UIlfortunately, even when optimum values of these parameters are selected, the
conventional techniques do not perform as well as wavelet smoothing. Signal processing
using wavelets is highly desirable because of the high quality of resolution obtained.
Moreover, this method does not involve hard coded filter constants or recording limits.
The degree of smoothing can be adjusted by simply altering the number of levels of
decomposition of the signal. The next chapter describes the mechanics of wavelet
smoothing, with a brief introduction to the Fourier transform.
Figure 2.12. Performance of sixth order Daubechies wavelets.
Original signal and smoothed approximation obtained using sixth order
Daubechies wavelets at the fifth level of decomposition.
34
Original signal and its smoothed approximation obtained using sixth order
Daubechies wavelets at the seventh level of decomposition.
0.50
0.55
0.60
0.70
0.50
0.60
......
Chapter 3
THE FOURIER AND WAVELET TRANSFORMS
This chapter describes the fundamentals associated with two different signal processing
techniques which we propose to use at OSU for automated trend extraction. The first
technique employs a Fourier or frequencydomain decomposition of a timedomain
signal. The material covered serves two purposes. The first is to describe the analytic
approach used to characterize a timedomain signal. The second is to provide an
introduction to the primary signal processing technique that we use for automated trend
extraction, wavelet signal decomposition. This chapter lays the foundation for the
automated trend extraction procedure proposed in the following chapter.
3.1 Digital Signal analysis and processing
The study of a signal as a function of time is called time analysis. Such an analysis can
provide information about the source and medium of propagation. However, if other
properties of the signal are to be studied, it would be better to study the signal using a
different representation. One of the most powerful representations is the spectral
representation, which is generated using a frequencydomain analysis. For stationary
signals, i.e., signals whose properties are statistically invariant over time, the most
common method to perform a frequency analysis is the Fourier transform. The Fourier
transform coefficients give a clear picture of the frequency content of the signal. This
35
36
method works well when the signal is fundamentally composed of a few periodic
components. However, if the signal is nonstationary, i.e., exhibits abrupt changes due to
unexpected transient events, a Fourier decomposition yields a wide range of frequency
components as analytically required to reproduce abrupt changes of the signal in time.
Many of the frequency components are not representative of the process which is
generating the signal but are simply analytical artifacts required to approximate sharp
transitions. Therefore an analysis of nonstationary signals requires more than the Fourier
transform.
Two techniques commonly employed for analysis of nonstationary signals
include wavelets decomposition and the ShortTime Fourier transform. The latter
analyzes a signal in (data) windows of fixed length, using the Fourier transform, to
determine the dominant frequencies prevalent in the signaL. A recently developed
technique is the wavelet transform [Rioul, 1991]. The ShortTime Fourier transform uses
a single analysis window whereas the wavelet transform uses short windows (contraction)
at high frequencies and long windows (dilation) at low frequencies. The contractions and
dilations are called scalings. The notion of scale in wavelet transform could be envisaged
as an alternative to frequency in ShortTime Fourier transform. This effectively means
that a signal is mapped into a timescale plane in wavelet transform as compared to the
timefrequency plane in the ShortTime Fourier transform. The most attractive feature of
the wavelet transform is that a signal can be analyzed at various scales and resolutions.
Thus, an explicit representation of the temporal features of a signal in a joint timefrequency
plane can be obtained.
.......
37
3.2 The Fourier transform
respectively.
JCt) = L,gkUk(t)+ L,hk~k(t) (19)
k k
the decomposition coefficients associated with the basis functions aCt) and ~(t)
sine functions, the Fourier series representation is obtained. The Fourier transform, when
magnitude for each of the frequencies covered. The original signal can be perfectly
applied to any signal, decomposes the signal into sine and cosine coefficients of various
When the basis functions of equation (19) are represented by the traditional cosine and
Here, k is the index defining the length of the signal. Various basis functions could be
adopted depending on the nature of the problem at hand. In equation (19), gk and hk are
weighted responses to each of the basis functions. Generally, a signal can be decomposed
into the sum of a pair of basis functions aCt) and pet).
way of solving many engineering problems. Fourier and wavelet decompositions both
utilize this approach. When a decomposed signal is input to a linear, time invariant
system, then superposition theory can be applied and the output obtained by adding up the
3.1.1 Signal Decomposition
Decomposition of a given signal into a weighted sum of basis functions i a conveni.ent

reconstructed by using all the sine and cosine coefficients in equation (19).
In this section, the Fourier series expansion for continuoustime periodic signals is
reviewed, followed by a discussion of the Fourier series expansion for discretetime
signals. Next, the definition and properties of the discrete Fourier transform (DFT) using
.......
38
using the DFf to perform linear convolution, which is the basic filtering operation.
(20a)
for all t (20)
k=oo
f(t) = I Fkejklrlo
k=~
by the following equations:
representation of a continuoustime periodic signalf(t). Here, the complex term ejkrOu
implicitly represents the sine and cosine basis functions by virtue of the definition
where j is the complex variable with the property f = 1. In equation (20), Fk are the
coefficients of the expansion, and Qo, the fundamental frequency. They are determined
way [Ludeman, 1986].
The expansion given by equation (20) is called the exponential Fourier series
A periodic, continuous time signal f( t) with period T can be expressed as a weighted sum
of a countable number of complex exponential continuous time functions in the following
interpreting DFf results in terms of frequency content of discretetime signals and by
transform, are discussed. Linear filtering and spectral analysis are addressed by
the fast Fourier transform, an efficient method for computing the discrete Fourier
3.2.1 Continuoustime Fourier series
(21 )
Qo =21rJT (2Ia)
39
A convenient trigonometric representation of equation (21) in terms of a weighted sum of
sine and cosine basis functions is
..
f(t) = a o + L{aA cosknot +bk sinkQot)
k=1
where the constants ao, ({k, and bl< can be determined as
(22)
1 T
an = Jf(t)dt
To
(23a)
2 T
ak = Jf(t)cosknotdt k =1,2, 00
To
(23b)
2 T
hk =  Jf(t)sin knot dt k = 1,2, 00
To
(23c)
of complex exponential sequences. By virtue of the fact that sinusoidal sequences are
A real, periodic discretetime signal x(n) of period N can be expressed as a weighted sum
(24c)
(24b)
(24a)
k=l, 2, .
Clo=Fo
The constants ak, and bk• of the trigonometric fonn, and Fk, of the complex exponential
form are related as follows:
3.2.2 Discretetime Fourier series
unique only for digital frequencies from 0 to 27[, the expansion contains only a finite
number of complex exponentials as follows:
where the coefficients of the expansion X(k) and the fundamental digital frequency ~ are
INI
x(n) =LX(k)eikWoJ'
N k=O
given by
(25)
40
NI
X(k) = Lx(n)eikWnn
11=0
for all k (26a)
An alternative form of the discrete Fourier series, for a periodic discretetime
Equations (25) and (26a) are called the discrete Fourier transform series (DFS).
(26d)
(26b)
(26c)
NI
X(k) = Lx(n)W~N
11=0
~ = 2rrIN
equations, for odd and even N.
This is called the trigonometric form of the Fourier series and is represented by two
signal in terms of the si ne and cosine basis functions can be written as shown below.
for a periodic discretetime signal. Equation (26) is sometimes written in the following
equivalent form:
The expression given in (25) is referred to as the exponential form of the Fourier series
EvenN
N/21 21t NI2I 2n: N
x(n) = A(O) + L A(k)cos(kn) + L B(k)sin(kn) + A()cosll1t (27)
k=1 N k=l N 2
The 1ast term contains cos nn:, which is just (1t, and is associated with the highest
frequency possible.
......
41
The coefficientA(O), is the most dominant Fourier coefficient and represents the
The relationships between the A(k), B(k) of the trigonometric form, and the X(k) of
(29c)
(28)
(29d)
(29a)
(29b)
k=/,2, ,(NI2)1
k=l,2, ,(NI2)l
INI
A(N /2) = L.x(n)cosn1t
N 11=0
1 Nl
A(O) =L.X(II)
N 11=0
2 NI 21t
A(k) =L.x(n)cos(kn),
N n=O N
2 NI 21t
B(k) = N L.x(n)sin(kn),
,,=() N
(Nl)12 21t (NI)/2 21t
x(ll) = A(O)+ L. A(k)cos(kn) + L. B(k)sin(kn)
k=1 N k=1 N
the discrete Fourier coefficients for a real x(n) with an even N are given by
average value of the original signal. The coefficientA(O) can also be considered to be the
OddN
The constants A(O), A(k), and B(k) for even N can be shown to be
If N is odd, equations (29a c) apply but A(NI2) =0.
coefficient at k =0 is zero and hence is not considered.
cosine coefficient at frequency index zero (i.e., atk =0). Note that B(O), the sine
A(O) =X(O)/N (30a)
A(k) = [X(k) + X(N  k)]/N,
B(k) = j[X(k)  X(N  k)]/N,
k = 1,2, ....... , (N/2)1
k = 1,2, ....... , (NI2)l
(30b)
(30c)
A(NI2) =X(N/2)/N (30d)
......
42
Again, if N is odd, equations (30a  c) apply and A(NI2) =O. The following example
illustrates the complex exponentials and the trigonometric forms of the discrete Fourier
representation. From Figure 3.1, x(n) is [J 2 3 4J
x(n)
4  r
 r3
 r
I ....,r Periodic
4 3 2 I 0 I 2 3 n
Figure 3.1. Example problem
The exponential form of the Fourier series is specified once the X(k) are calculated from
equation (26a). Before carrying out these calculations using equation (26a), the powers of
WN must be determined for N =4. Here, N =4 because the period associated with this
signal is 4. WI is given by
W4 =e j (2rr/4) = cos(7t /2)  j sin(n /2) =j
The other powers of W.; are easily seen to be
W4
2 =W4
1
•W4
1 =( j)( j) =I
W;' =W4
2
•W4
1 = (1)( j) = j
W4
4 = W4
3
•W4
1 =(j)(  j) = I
Using these powers of W.;, the X(k) for k =0, I, 2, 3 are calcluated as follows:
43
41
X (0) = L x(n)W40
11 = X(O) + x(l) + x(2) + x(3)
=1+2+3+4=10
41
X(1) =L x(n)W4
11l =X(O)W40 + X(l)W4
1 + X(2)W4
2 + X(3)W4
3
1/=0
=1( I) + 2( j) + 3(1) + 4(j) = 2 + 2 j
41
X (2) =L x(n)W/ol/ = x(O)W40 + x(l)W4
2 + X(2)W4
4 + X(3)W4
6
11=0
=1(1) + 2( 1) + 3(1) + 4(1) =2
41
X (3) =L x(n)W4
3o
/l =x(O)W40 + x(1)W4
3 + X(2)W4
6 + x(3)W4l)
1/=0
= 1(1) + 2(j) + 3(1) + 4( j) =2  2j
Now, the trigonometric coefficients are calculated from the exponential coefficients as
shown below, using equation (30). As N = 4, the number of cosine coefficients would be
(N/2 +/) =3 and the number of sine coefficients would be (N/2 /) =1.
10 5
A(O) = X(O)/ N ==4
2
A(l) = [X(l)+ X(41)]/4= [(2+2)+(22)]/4=1
B(l) = j[X (I)  X (4 1)]/ 4 = i[(2 + 2j)  (2  2j)]/ 4 = I
4
A(2) = X() = X (2) /4 =1/ 2
2
Now, x(n) can be determined using these trigonometric coefficients as shown, using
equation. (20) for even N.
......
(412)1 21t
x(O) = A(O) + L A(k) cos(k0)
Ie=l 4
(412)1 2n
+ L B(k)sin(kO)+ A(4/2)cos(07t)
1e=1 4
= A(O) + A(1) cos(07t) + BO) sin(On)
=(512) + (1)1 + 0 + (112)
=1
Similarly,
(412)1 2n:
xCI) = A(O) + L A(k) cos(k1)
Ie=! 4
(412)1 2n:
+ ~ B(k)sin(k 4 I) + A(4/2)cos(1n:)
= A(O) + A(l)cos(n /2) +B(l) sin(n /2) + A(2)cos(ln)
= (512) + (1)0 + (1)1 + (112)(1)
=2
(412)1 2n:
x(2) = A(O) + L A(k) cos(k 2)
1e=1 4
(4/2)1 2n:
+ L B(k)sin(k2)+A(4/2)cos(2n)
k=! 4
=A(O)+ A(l)cosn +B(l) sinn +A(2)cos2n:
= (512) + (1)(1) + (1)0 + (112)/
=3
(412)1 2n:
x(3) = A(O) + :L A(k) cos(k  3)
k=1 4
(4/2)1 2n:
+ :L B(k)sin(k3) + A(4/ 2)cos(3n:)
Ie=l 4
44
......
45
= A(O) + A(1) cos(37t /2) + B(1) sin(ll: /2) + A(2) cos(ll:)
=(512) + (1)0 + (1)(1) + (112)(1)
=4
From the above calculations it is apparent that the X(k)'s can be written in matrix form a
follows:
X(O) WO WO WO WO x(O) 4 4 4 4
X(l) WO Wi W2 W3 x(l) 4 4 4 4 = WO (31 )
X(2) W2 WO W2 x(2) 4 4 4 4
X(3) WO W3 W2 WI x(3) 4 4 4 4
The entries in the coefficient matrix have been reduced in powers of W4 by using the fact
that W4
4 =1. Fast ways of making the calculations above can be shown to be equivalent
to decomposing the square matrix into products of relatively sparse matrices. This idea is
exploited further in the fast Fourier transform.
The Fourier transform is an excellent tool to determine the frequency content of a
signal. Figure 3.2 shows a raw sensor signal and its Fourier decomposed sine and cosine
coefficients. Figure 3.2u shows the sensor signal in the time domain i.e., the magnitude
of the parameter measured at various instants of time. Figures 3.2b and 3.2c show the
same signal transformed to the frequency domain by the Fourier transform. The cosine
and sine coefficients obtained by performing the Fourier transform operation on the time
domain representation of the sensor signal are a measure of the contribution from the
different frequencies which make up the original signal.
The magnitude of the average value of the signal (cosine coefficient at frequency
index zero), A(O), is 0.5977. The sine coefficient, B(O), at frequency index zero is zero.
The cosine and sine coefficients have not been drawn to the same scale so that the reader
.....

46
0.70
0.60
0.50
50 100 150 200 250 300 350 400 450 500
(b) Cosine coefficients against frequency index. Most dominant
Fourier coefficient, A(O) =0.5977. is not shown in plot.
101 151 201 251
Frequency index
51
(a). Original signal.
0.006
0.005
~c 0.004
eLl 'u 0.003
E
g 0.002
u
~ 0.001
V> 8 0.000
0.001
0.002
0.06
0.05
CIl i: 0.04
.~
u 0.03 lEeLl
0 0.02 u
eLl c
i;) 0.01
0.00
0.01
51 101 151 201
Frequency index
251
(c) Sine coefficients against frequency index.
Figure 3.2. Original signal decomposed into cosine and sine components.
.....
47
can appreciate the difference in their magnitude. From Figures 3.2b and 3.2c, it can be
observed that the first few coefficients are relatively large in magnitude. These
coefficients capture the amplitude and frequency of the main body of the signal. The
other relatively small coefficients represent the amplitude and frequency of sensor noise
introduced in the signal. Thus, the Fourier transfonn gi ves a complete picture of the
frequency content of the signal.
The Fourier transform is an excellent tool that could be used to alternate between
the time and frequency domains of a signal representation. However, it fails to provide
infonnation about the time and frequency domains of a signal simultaneously i.e., the
Fourier transfonn gives infonnation of the different frequencies that existed for the total
duration of the signal but not the frequencies which exist at a particular time.
Simultaneous information on the time and frequency domains of a sensor signal
representation is known as timefrequency localization and is explained in the following
section.
3.3 Timefrequency localization
The need for a simultaneous timefrequency analysis is discussed below by considering a
progression of examples. If a simple sine wave that lasts forever is considered, then to
describe it completely it would suffice to say that it is a sine wave of fixed magnitude for
a particular frequency for all time. The timefrequency description would show only one
particular frequency at all times. Now, if the sum of three sine waves that last forever,
with different frequencies are considered, the timefrequency spectrum would show three
different frequencies present for all times. Thus time has played a passive role since
.......

48
whatever is occurring, occurs for all time, and the standard amplitudetime spectrum
provides a satisfactory description [Cohen, 1995].
Figure 3.3a illustrates the case where the signal is formed by concatenating three
different sinusoidal components. To fully describe such a situation, the frequencies for
each time ought to be given. The timefrequency plot for this illustration (top left) shows
three different frequencies and their duration. The power spectrum shown just below the
timefrequency plot shows the frequencies 1, 2, and 3 but does not tell when they existed.
Figure 3.3b is an extension of the previous illustration. This timefrequency plot shows
that there are times when all the three frequencies coexist (time index 4 to 6.5), only two
of the frequencies exist (time index 2 to 4 and 6.5 to 8.5), and only one of the three
frequencies exists (time index 0 t02 and 8.5 to 10.5). This information cannot be
extracted from the power spectrum only. Figures 3.3c and 3.3d illustrate this idea further.
Note that the power spectrum is roughly the same in all four parts of Figure 3.3 and
indicates the different frequencies exist.ing but does not give any idea when the different
frequency components were present.
A technique that gives a good timefrequency description is needed. The Fourier
transform is only capable of identifying the frequencies present over the total duration of
the signal and not the frequencies that exist at a particular time. This is where the beauty
of the wavelet transform comes into play. Wavelet transforms provide for goodtime
frequency localization. The next section describes briefly the mathematical basis of the
wavelet transform, construction of simple wavelets, multiresolution analysis, and signal
decomposition and reconstruction.
49
~
10
:?: 9
~ 8 <.
7
Cl.I 6
.5 5 I<
4
3
2
1
0
1 2 3
(a) Frequency
Power spectrum
10
9
8
7
4> 6 .I5< 5 i
4
3
2
1
0
1 2 3
(b) Frequency
Power spectrum
Figure 3.3. The schematic timefrequency localization plots of finite duration sine
waves. From the power spectrum, it is not possible to estimate the duration of
waves. A timefrequency plot clearly shows the frequencies existing for each time.
4
50
~
10
9
8
j~ 7
6
5
.S 4
3
2
1
0
1 2 3
(c) Frequency
Power spectrum
10
~
~ 9
~ 8
7 ~~ 6
5
4 ~ 3
<; 2
~ 1
0
1 2 3
(d) Frequency
Power spectrum
Figure 3.3 (continued)
51
3.4 The wavelet transform
decomposition, the basi' functions are the sine and cosine functions. As indicated
The wavelet transform is analogous to the Fourier transform. In the Fourier signal
previously, a discrete signalf(t) can be represented as foHows:
f (t) =I gk cos(rk) +I 17 k sin(tk)
k k
(32)
difference between the basis functions of the Fourier transform and the wavelet tran form
is that the scaling and wavelet basis functions are complex functions which are derived
scaling function. An introduction to the discrete wavelet transfolln is given below in
f(t)= Ig'(Xk(t)+ Ih,~,(t) (33)
k Ir.
Here, a and ~ are the scaling and wavelet basis functions, k is the index which defines the
length of the signal, and Rk and hk are the decomposition coefficients. An important
and do not occur naturally. Furthermore, the wavelet basis function is derived from the
function and the wavelet function.
In a similar fashion, the Wavelet series representation has two basis functions, the scaling
which the scaling and wavelet functions are discussed in more detail.
3.4.1 Discrete wavelet transforms
The discrete wavelet transform, like the fast Fourier transform, is a fast, linear operation
that operates on a data vector whose length is an integer power of two, transforming it
into a numerically different vector of the same length. The wavelet transform, like the
fast Fourier transform is invertible and, in fact, orthogonal [Press, 1992]. The inverse
transform, when applied as a matrix, is simply the tran pose of the transform matrix

52
[Strang, 1993, 1989]. Unlike the Fourier transform which is uniquely defined by ine and
cosine basis functions, there is not one single unique set of basis functions for the wavelet
transform; in fact, there an infinite number of possible basis functions. Each of these
scaling and wavelet basis functions has a unique representation over a definite interval
and vanish outside this interval. In wavelet mathematics, different scaling and wavelet
functions are obtained by specifying the family and order of the wavelet. A brief
explanation of the scaling and wavelet functions is given below.
3.4.1.1 Scaling functions and wavelet functions
Scaling functions and wavelet functions are represented by special types of equations
called dilation equations. The general form of a dilation equation is as follows:
(34)
Here,l/( is a vector of dilation function coefficients. The factor two provides for the
dilation, i.e., expansion or contraction whereas the factor k provides for translation.
Thus, frequency and time localization can both be achieved. The general equation of a
scaling function is given by
(35)
Here, c/( is the scaling function coefficient. The factor j provides for dilation and k
provides for translation. A very important point to note here is that the scaling function is
expressed as a sum of dilations and translations of the scaling function itself.
The wavelet equati~m is derived from the scaling function by taking "differences"
as shown
Daubechies wavelets are a distinct family of wavelets. These wavelets are orthogonal and
Motard, 1994]. A very brief introduction to the Daubechies family of wavelet is given in
[Daubechies 1992, 1993].
53
W(x) =Ld (36) k (2xk)
k
There are many different wavelet families; a discussion of which is beyond the
the following section. For more information, the interested reader may refer to
are compactly supported. Support. is defined as the span of the scaling or wavelet
Here, the wavelet coefficients are given by dk• Wavelets are defined as functions of the
function over which the function has a nonzero value. They are constructed in the
purview of this work. The interested reader may refer to [Kaiser, 1994; Walter, 1994;
scaling functions because the decomposition basis functions have to be interdependent.
3.4.2 The Daubechies Family

frequency domain where it is easy to handle the dilation and translation parameters.
Certain constraints are imposed on the construction of the scaling and wavelet functions
depending on the properlies of the wavelet desired. These constraints are ultimately
reflected in the respective sets of coefficients, Ck and d/;. For the Daubechies family, the
constraints are orthonormality and orthogonality as explained below. Note that the
following is a summary overview of the material presented in [Daubechies 1992, 1993].
A function (t) is said to be orthogonal if the following relation is satisfied:
(37)

54
unity.
where c is a constant. The function is said to be orthonormal when c take the value of
For the Daubechies family, the following two conditions are imposed. The
J(x) = I (38)
scaling function, must be orthonormal.
The wavelet functions must satisfy the condition
JW(x)dx = 0
Integrating equation (35) and applying the orthonormality condition equation (38)
Similarly, integrating equation (36) with the condition described by equation (39)
(39)
(40)
(41 )
(42)
(43)
The conditions for a wavelet to be both orthogonal and orthonormal are that the
sum of it's scaling function coefficients should be two and the sum of its wavelet function
coefficients should be zero. When the scaling function and the wavelet function are
quadrature mirror filters of each other [Akansu, 1993; Cohen, 1992a; Daubechies, 1988]
the relationship between their coefficients is given by
(44)

55
This implies that the wavelet coefficients are obtained directly from the scaling function
coefficients by reversing the order and changing the sign on every alternate coefficient.
Thus the tedious process of calculating them separately is avoided.
A brief description of the construction of first order Daubechies wavelets is given
in the following section. A very detailed treatment on this topic is available in [Akansu,
1993; Chui, 1992a; Daubechies, 1992; Haykin, 1991; Strang, 1989].
3.4.2.1 Construction of first order Daubechies wavelet
The first step in the construction of wavelets involves the computation of the scaling
function and wavelet function coefficients. Once these coefficients are computed, the
actual scaling and wavelet functions are then determined. The following describes the
least complicated of the Daubechies family of wavelets. The first order wavelet,
commonly referred to as the Haar wavelet, is constructed below.
The scaling function is a simple box function given by
{
I O~x~1
(x) = (2xk)
k
W(x) =do (2x)  d} q>(2x J)
The wavelet equation is given by
As k has two values 0 and 1, W(x) reduces to
by a unit value.
Equations (35) and (39) are both satisfied when do =I and d} =I. So,
as x lies between 0.5 and I. Substituting the orthogonality property in the above
between zero and one i.e., x lies between zero and 0.5. Similarly, (2xl) exists as long
This is where the beauty of orthogonality comes into play. The translates of the dilation
equation, the values of both Co and c} are obtained as I. This shows that the box function
equation are orthogonal. (2x) and $(2x1) are orthogonal. q>(2x) exists as long as 2x lies
W(x) =(2x) (2x1) (48)
The wavelet function takes the following values
{
I, 0 S; x S; I /2
W(x)= I, 1/2$x$1
0, otherwise

The
Tl:
I'
57
wavelet function is shown below.
I~ L 
0
0 0 () .
I 112 112 I I,
J
""
I  
coefficients in two different vectors a and b are considered, then the coefficients of the
is the multiplication of the coefficients of polynomials. If two polynomials with their
polynomial obtained by convolving these two polynomials is represented by a vector c,
length( a) + lellJ{th( b ) 1 (49)
Figure 3.4b. The Haar wavelet and its dilated and translated versions
the length of which is given by
Mathematically, signal filtering is performed by the convolution operation. Convolution
3.5 Signal filtering using scaling and wavelet functions
The kth element of the convolution product of c is given by
Ck =La(j)b(k + 1 j) (50)
The sum is over all the values ofj, which is from 1 to (length( a) + length( b ) 1). When
both vectors a and b are of the same length, n
c( 1) = a(l )b(1)
c(2) = a(l )b(2) + a(2)b( 1)
e(3) = a(/ )b(3) + a(2)b(2) + a(3)b(l)
58
c(k) =a(l )b(k) + a(2)b(kl) + + a(n)b(l)
c(2n2) = a(nJ )b(n) + a(n) b(nl)
c(2nJ) = a(n)h(n)
Wavelet decomposition is performed by convolving the signal with the scaling
and wavelet function coefficients. The scaling and wavelet function coefficients are also
known as filter coefficients. The resulting signal after convolution is dyadically sampled
i.e., every alternate sample is taken into consideration.
Iff is the signal, Hand G are the low pass and high pass filters respectively, then
the signal is decomposed as follows
A high pass filter retains the low frequency part and allows the high frequency parl of the
signal to pass through. Similarly, the low pass filter retains only the high frequency part
of the signal and allows low frequency part of the signal to pass through. Here, H is the
filter containing the scaling function coefficients and G is the filter containing the wavelet
function coefficients. The blurred signal contains the low frequency part of the original
signal and the detail signal retains the high frequency part of the signal. The notation "*,,
represents the convolution operation followed by downsampling in which every alternate
value is retained. The information content of the blurred and the detail signals is
mutually exclusive, and the original signal at any level of decomposition can be obtained
(51 )
(52)
Blurred signal =H*f
Detail signal =G*f
59
principle of orthogonality.
The original signal is reconstructed by reversing the decomposition operation.
HT"(Blurred signal) + CT"(Detail signal) =original signal (53)
by the combination of the blurred and detail signals at that level. This results from the
HT and CT are the transpose of Hand G. The notation "",, represents the upsampling
representation of the decomposition and reconstruction of a signal.
signal by inserting zeros between each value. Figure 3.5 illustrates the basic
Reconstructed
Signal
G
H
Original Signal
operation and convolution with the signal. Upsampling is doubling the length of the
Figure 3.5. Basic decomposition and reconstruction representation
Consider a signalfhaving n samples at its original resolution which is represented
in the time domain as the vector aD. At the original resolution, the elements of the vector
aO (the digitally sampled signal) are the values of the signal itself. At the first level of
decomposition, the decomposition coefficients are given by
j = 1, n/2; k = I, n (54)
_ ............
equally into blurred and detail signals). At any level of signal decomposition the values
constant Y2 is the normalization constant (this takes into account that the signal is divided
60
j =1 1112; k =1, n (55)
Here, Ck and dk are the scaling and wavelet function coefficients respectively. The
of a and h are computed by recursion from their values at the previous level. It is to be
noted that a signal can only be decomposed nI2 levels because at the nl2tll level, the
blurred and detail signal would contain only one coefficient each.
The decomposition procedure is illustrated below using the Haar wavelet and a
short signal so as to give the reader a good feel of this subject. The Haar scaling
coefficients [co CJ] are [J I] and the Haar wavelet coefficients [do dl ] are [1 1]. The
signal is given by [1 234]. The decomposition coefficients are given by
The blurred signal at the first level is al =[3/2712] and the detail signal at the first level
is hi =[1/2 112]. Reconstruction is performed as follows:
af = L..a~C2j_k + L..b;d2i k
j j
j = J, nJ2. (56)
For this example,
 
61
a~ = [1 x 3/2 + () x 7 / 2] + [1 x 1/2 + 0 x I] = I
a~ = [1 x 3/ 2 + 0 x 7 / 2] +[l x 11 2 +0 x I/2] = 2
a~ = [Ox3/2+ Ix7/2]+[Oxl+ 1I2xI]= 3
{l~ =[Ox3/2+lx7/2]+[OXl+lxI/2]=4
This recreates the original signal aO =[J 2 3 4].
The above example illustrates perfect reconstruction. For smoothing purposes,
pelfect reconstruction is not used. Instead, the high frequency part of the signal i'
removed. If the high frequency part is to be removed, the detail signal coefficients are
made zero and the signal reconstructed from the blurred signal coefficients and the zeroed
detail coefficients. In the above example, the level of decomposition is one. However, a
signal could be decomposed to more levels than one and viewed at each level of
decomposition (resolution). The process can be repeated until a single blurred nonzero
coefficient is obtained.
For smoothing purposes, the key is determining how many levels of
decomposition should be performed. The most straightforward approach is as follows.
First the signal is decomposed one level, setting the detail signal coefficients to zero (this
eliminates the high frequency part), and reconstructing the signal and checking the degree
of smoothing. If the reconstructed signal still contains some high frequency noise, the
blurred signal at the first level is decomposed again. The detail coefficients at the second
level are again made zero and the signal reconstructed from the blurred and detail signals
at the second level. This process is repeated until the desired degree of smoothing is
62
achieved. This process of viewing the signal at multiple resolutions is called multiresolution
analysis. A brief description of multiresolution analysis is given below.
3.6 Multiresolution analysis
As explained above, the process of decomposing a signal to various levels is known as
multiresolution analysis. The detail and blurred signal coefficients at any level are
computed by recursion from the results at the previous level. At any level, the blurred or
the detail signals are considered to be an averaged version of the signal (blurred or detail)
at the previous level. The frequency of the signals is twice that of the frequency at the
previous level (also known as a scale twice as large). If the blurred signal at any level is
considered the averaged version of the blurred signal at the previous level, then the detail
signal is obtained from the difference between the signal at the previous level and the
blurred signal generated at the current level, i.e., the infonnation present in the original
signal but filtered out in the averaged version. In other words, the blurred signal at any
level j is obtained by the combination of the blurred and detail signals at the next lower
level j+ / and the detail signal at level j from the difference between the blurred signals at
level j and level j J.
The basic idea behind multiresolution analysis can be summarized by Figure 3.6.
.....
Detail signal at level 3
Blurred signal at level 3
G
H
Detail signal at level 2
H
G
G
Original
Signal(level 0)
Figure 3.6. Basic representation of the multiresolution algorithm
Application of multiresolution analysis to an actual plant signal is illustarted in
63
H
Detail signal at level I
level is the sum of the detail and blurred signals at the previous or next higher level.
Figure 3.6 shows a signal decomposed three levels. It is evident that the signal at any

Figure 3.7. The original signal is decomposed into blurred (low frequency) and detail
(high frequency) signals at the first level. It can be seen that the low frequency or blurred
signal basically follows the fundamental trend of the original signal with magnitude
almost the same as that of the original signal. However, the high frequency part of the
original signal does not follow any particular trend relative to the trend of the original
signal and is of magnitude substantially lesser than that of the original signal. Subsequent
results after two, three, four, five and seven levels of decomposition are presented in
Figure 3.7.
The reconstructed signal with the detail coefficient zeroed is superimposed on the
original signal is presented in Figure 3.8. From Figure 3.8, it can be seen that the
......
0.70
0.65
0.60
0.55
0.50
64
Original signaL.
0.02
0.00
0.02
Detail coefficients at first leveL.
0.80
0.60
0.40
Blurred coefficients at first level.
O.Ol
0.00
0.01
0.80
0.60
0.40
Detai I coefficients at second level. Blurred coefficients at second level.
0.02
0.00
0.02
0.80
0.60
0.40
Detail coefficients at third level. Blurred coefficients at third level.
Figure 3.7. Detail and blurred coefficients for the ffrst four levels of decomposition
using sixth order Daubechies wavelets.
....

0.02
0.00
0.02
0.70
0.60
0.50
65
. ,
Detail coefficients at fourh level.
0.02
0.00
0.02
Detail coefficients at fifth level.
Blmred coefficients at fourth level.
0.70
0.60
0.50
Blurred coefficients at fifth level.
0.02
0.00
0.02
0.70
0.60
0.50
Detail coefficients at seventh level. Blurred coefficients at seventh level.
Figure 3.7. (contd.) Detail and blurred coefficients for fourth, fifth, and seventh
levels of decomposition using sixth order Daubechies wavelets.
....
F
0.70
0.60
0.50
0.70
0.60
0.50
0.70
0,60
0.50
First level of decomposition.
Second level of decomposition.
Third level of decomposition.
Figure 3.8. Smoothed signal obtained (with detail coefficients set to zero)
using sixth order Daubechies wavelets, superimposed on original signal for
various levels of decomposition.
66

0.70
0.60
0.50
0.70
0.60
0.50
0.70
0.60
0.50
Fourth level of decomposition.
Fifth level of decomposition.
Seventh level of decomposition.
Figure 3.8. (contd.) Smoothed signal obtained (with detail coefficients set to zero)
using sixth order Daubechies wavelets, superimposed on original signal for various
levels of decomposition.
67
...

68
smoothed signal obtained from reconstruction at the fourth level of decomposition is an
excellent representation of the raw sensor signal. Further, it can be seen that as the level
of decomposition increases beyond the optimum level of decomposition, the capacity of
the smoothed signal to retain the fundamental trend of the original raw signal is reduced.
The smoothed signals reconstructed from beyond the fifth level of decomposition prove
this beyond doubt.
Multi  resolution analysis [Cohen, 1992b; Daubechies, 199]; Mallat, 1989a;
Mallat, 1989b; Mallat, 1989cl is a very powerful tool for trend extraction and pattern
recognition. Accurate representation of the fundamental trend of the original sensor
signal can be obtained from decomposition of the signal followed by reconstruction
without the detail coefficients. However, increasing the number of levels of
decomposition also increases the degree of smoothing; so an optimum level must be
selected so that a smoothed representation is obtained without sacrificing the fundamental
trend of the original signal. The optimum level of decomposition was previously
determined by trial and error, i.e., by experimenting the signal with various levels of
decomposition and finding the level of decomposition which best smoothes the raw
sensor signal.
The optimum level of decomposition has previously been determined by
experimentation and depends on the pattern recognition problem being addressed. This
can be avoided and the optimum level or the level very near to the optimum level of
decomposition can be determined automatically. The next chapter presents our proposed
method for automated trend extraction.
..

Chapter 4
AUTOMATED TREND EXTRACTION
4.1 Need for automated trend extraction
In Chapter 3 the desired level of smoothing of signals is achieved empirically by trial and
error. Experimentally determining the optimum level of smoothing of a raw sensor signal
is not suitable for realtime applications.
An automated approach is needed if we expect to provide a stand alone system
which allows operators to create their own process monitoring applications. Furthermore,
different signals may require different degrees of smoothing depending on the
characteristics (e.g., noise content) of individual sensors. Ideally, we would like to have a
system which analyzes the characteristics of the monitoring application as well as the
sensor signals and recommends the appropriate degree of smoothing for each signal. The
system would graphically illustrate the recommended degree of smoothing and allow the
user (operator) to either accept the recommendation or specify more or less smoothing.
There is no a priori correct level of smoothing, as the desired degree of smoothing
is application dependent. We have empirically determined that the degree of smoothing
for patternbased analysis is usually determined by the length of the data window used for
observing and monitoring a process. Typically, it suffices to provide sufficient smoothi ng
so that approximately four blurred coefficients are used over the span of the window
length after decomposition and before reconstruction. The reconstructed signal can then
69
...
70
be sampled as desired to obtain a compact repre entation, as de cribed in Chapter 1. In
other words, the "correct" number of levels of decomposition normally depends on the
number of blurred coefficients after decomposition that are fixed over the window length
used for process monitoring.
4.2 Determination of the degree of smoothing from the user input window length
The smoothing obtained for a typical sensor signal fixing four blurred coefficients
after decomposition, over a 45 minute window of observation is shown in Figure 4.1.
The number of levels of decomposition of the raw signal is back calculated from the
window length used for process monitoring, the sampling period of the raw data, and the
desired number of blurred coefficients over the window of observation, after
decomposition.
The window length for process monitoring is problem specific and left to the
discretion of the operator (user). Nevertheless, a minimum window length required for
monitoring purposes can be helpful in providing insight when specifying the actual
window length. Determination of the minimum window length for process monitoring is
described later in the chapter. The following discussion describes how the recommended
degree of smoothing is determined from the window length.
The user input window length, winx (units of time) is first converted to the
associated number of data points with the aid of the signal sampling period, Ts and this
value is checked to see if it is less than half the length of the original signal, N. If the
user input window length is less than NI2, then winx is used to calculate the number of
levels of wavelet decomposition, otherwise the user is prompted to enter a window of
.....
71
I
I
I
I
.I 
I
I
I
I
~I 
I
I
I
I
I
I
I I I
~~~~. . I I I
I I I
I I
I I
I I
I
I ~.
Monitoring Window
Time ._..;:,.
I
I
I
iI 
I
I I
I I
 LI LI I ~""AII~I ~_____ _ ~ __
I I I
I I I
I I I
I I I
 IL IL ~I _
I I I
I I I
I I I
I I I
I I I
I
I
I
I
jI 
Original sampling frequency
~
a • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •
Decomposition
level
1 • • • • • • • • • • • • • • • •
•
•
• •
•
• •
.~
Desired frequency for smoothed signal
•
•
•
•
•
•
4
3
2
5 •
Figure 4.1 Span of scaling functions in the process monitoring window and
the desired frequency level for the smoothed signal.
...
...
72
length less than N*T.J2 units of time. Figure 4. I illustrates how the desired number of
levels of decomposition are determined so that there are approximately four points
spanning the window length. Analytically, the number of levels of decomposition,
aceLevels can be calculated according to the following formula,
act_levels =Log2(winx/4*T.rJ (57)
It is to be noted here that, for a specific problem, winx is a constant and acelevels
depends only on the sampling frequency. The number of levels of wavelet decomposition
obtained, acelevels, gives the best estimate of the optimum degree of smoothing.
Nevertheless, the user must verify this result and be given the option of specifying more
or Jess smoothing if deemed necessary.
The algorithm is built around the Matlab smoothing module first proposed by Mr.
V. Raghavan [Raghavan, 1994]. This smoothing module uses sixth order Daubechies
wavelets. The following discussion documents the Matlab mfiles which are utilized.
The scaling and wavelet function coefficients are determined by the function
daub.m. All smoothing techniques involve convolution which tend to distort the trend of
the smoothed signal towards the ends. Distortion of the smoothed signal towards the end
is unacceptable for patternbased process monitoring purposes. To avoid trend distortion,
the original signal is padded on either side by half its length by the function net.m
[Raghavan, 1994]. The padded signal is then decomposed to aceLeveLs number of levels
by the functionJwt.m and reconstructed with the blurred coefficients at this level. All
detail coefficients are set to zero. Signal reconstruction is then performed by the function
ifwt.m. The extensions to the signal are then removed to obtain the smoothed
representation of the original signal. The user is then prompted if the smoothed
..

73
representation is satisfactory, too smooth, or still needs to be smoothed. Depending on
the input of the user, appropriate action is taken.
4.3 Determination of the minimum degree of smoothing from raw signal
characteristics
The number of levels of decomposition obtained from the user input window length
generally produces the desired level of smoothing. However, in special cases where the
signal contains a steady, low frequency oscillation produced from control system
interactions, more smoothing may be required. In this situation, it is necessary to identify
the dominant low frequency component and ensure that the signal is decomposed beyond
this level.
The frequency spectrum of the original signal is useful for identifying the
dominant low frequency component. In this work, the term "cutoff' frequency is
introduced to differentiate the first few low frequencies of large magnitude that are of
interest, and other frequencies of smaller magnitude that are not of interest. For our
work, we assume that the cutoff frequency is the dominant frequency and the two terms
can be used interchangeably.
In this work, the term frequency index is often used. This is just another way of
representing frequency [Ludemann, 1992]. Referring to equation (29b) in Chapter 3, we
note that the frequency index k corresponds to the frequency "k2n IN" in a Fourier
decomposition.
..
74
4.3.1 Dominant frequency identification
This section describes how we quantify the cutoff frequency and use this information to
determine the minimum amount of decomposition required to smooth through low level,
steady oscillations.
As described in Chapter 3, the discrete Fourier transform (OFf) can be applied to
split a signal to obtain the frequencies associated with the signal in terms of cosine and
sine coefficients of various amplitudes and frequencies. We utilize the properties of the
discrete Fourier transform to determine the cutoff frequency index.
Figure 3.2 in the previous chapter shows the cosine and sine coefficients plotted
against their respective frequency indices for a representative time domain signal. By
visual inspection of Figure 3.2, it can be seen that the first few cosine coefficients are of
relatively large magnitude. The cosine coefficients around the frequency index forty and
above are of smaller magnitude and can be categorized as noise for purposes of
determining the cutoff frequency. The maximum sine coefficient occurs at frequency
index one. Sine coefficients from two to four are comparable to the maximum value but
coefficients above frequency index four can be ignored. This is an illustrative example
where the sine and cosine coefficients do not drop off at the same frequency index.
The question that arises now is the basis on which the cutoff frequency index is
to be determined, i.e., to select the sine or the cosine cutoff frequency index. Also, there
may be a few coefficients well within the cutoff frequency that are of very small
magnitude. This questions the validity of the cutoff frequency index itself. An
alternative approach, the cumulative power spectrum, gives a fairly good idea of the cutoff
frequency index. The next subsection addresses this topic.

75
4.3.2 The CumuJative power spectrum
The power spectra of a signal is determined from the output of a Fourier decomposition.
The value of the power of a signal (consisting of real values) at any frequency is
determined by taking the sum of the squares of the real and imaginary parts of the Fourier
coefficient obtained at that frequency and dividing it by the length of the signal. The
cumulative power is determined by adding the values of the power at each ucceeding
frequency. This is demonstrated by the following example.
The Fourier coefficients of the signal shown in Figure 3.1 of the previous chapter,
x(n) =[1 23 4J. are given by
X(O) = 10, X( 1) =2+2j, X(2) = 2, X(3) = 22j
The power at each of these points is determined as shown:
pro) =(101/4 =25
P( /) = [(21 + (2/1/4 = (4 + 4)/4 =2
P(2) = (21/4 = 4/4 =1
p(J) = [(21 + (211/4 = (4 + 4)/4 = 2
The cumulative power at each of these points is determined as shown below:
prO) =25
P(I) = prO) + P(I) = 25 + 2 = 27
P(2) = prO) + P(l) + P(2) = 25 + 2 +1 = 28
P(3) = prO) + P( 1) +P(2) +P(3) = 25 + 2 +1 +2 = 30
The cumulative power spectrum for the above signal is {25 2728 30}.
In a similar manner, the cumulative power spectrum of any signal can be
determined. However, the cumulative power spectrum obtained by this method is not
..
76
very smooth but is characterized by small changes in magnitude at some points. A
smoother estimate of the power could be obtained by using Welch's averaged
periodogram method. The signal is divided into overlapping sections, each of which is
detrended and windowed by a Hanning window [Oppenheim, 1975]. The adjacent
records of the power are averaged so as to obtain more reliable spectral estimates. All the
values thus obtained are in a confidence interval of 95%. The cumulative power
spectrum obtained by this method, and the source signal are shown in Figure 4.2.
As shown in Figure 4.2, the cumulative power spectrum increases steadily and
then flattens. The reason for the initial steep increase is that the first few values of the
power are relatively large due to the large magnitude of the first few low frequency cosine
and sine coefficients. Subsequent coefficients are of relatively smaller magnitude and
contribute less power. This effect is reflected by the fI attening of the cumulative power
spectrum after the first few low frequency coefficients. For our purposes, we define the
frequency index where the cumulative power spectrum "bends" as the cutoff or dominant
frequency index. The cutoff frequency index is a key parameter because it is used to
detelmine two important parameters,
(J) the minimum number of levels of decomposition of the signal so as to retain
the dominant frequencies and
(2) the minimum window length required for effective process monitoring, which
can be used as a basis for selecting the actual window length for process monitoring.
The cutoff frequency is determined by the examination of the second derivative
of the CPS. The interested reader may refer to the documentation in the function
delecl.m, a copy of which is included in the Appendix.
77
0.55
0.55
0.54
0.54
(.J
1~~l
s"
0.53
C...t~,
g'~
10 20 30 40 50 60
:,. , ~
I ,
'0"
(a) Signal. f...~J
:~;.
t;H ( __ I
.::;,
'1~_'i ii:l
I f .... "  0.00007 ..~
'.£
0.00006
0.00006
0.00005
0.00005
0.00004
0.00004
5 10 15 20 25 30
(b) The cumulative power spectrum obtained for the above signal.
Figure 4.2. Source signal and its cumulative power spectrum.
( ....

78
4.3.3 Determination of the minimum number of levels of signal decomposition from
the cutoff frequency index
Figure 4.3 lays the foundation for understanding how the cutoff frequency index
determines the minimum number of levels of wavelet signal decomposition.
Original Signal 211
Blurred Coefficients
2n_1
FirST level
Detail Coefficients
2n!
Detail Coefficients
2n2
Detai t Coefficients
211 ;
Blurred Coefficients 2',2
Second level
Blurred Coefficients
2n;
ilh level
nih level
Blurred Coefficients
J
Derail Coefficients
I
Figure 4.3. Multi resolution tree analysis for a signal of length N = 2/j
Consider a signal of some finite length N, where N is equal to 2/1. After
petforrning the fast Fourier transform on this signal, the number of cosine and sine
coefficients obtained are (NI2J+ I and (N/2)1 respectively. Now if the original signal is
decomposed fully to n levels of decomposition using the wavelet transform, one blurred
and a total of Nl detail coefficients would be obtained. The blurred and detail
coefficients at the nth level approximate the original signal with scaling and wavelet
79
functions which can be characterized by the lowest nonzero Fourier frequency. From the
definitions in Chapter 3, the lowest nonzero Fourier frequency corresponds to (1(21t1N))
and is associated with the a] coefficient. Therefore, a frequency index of I (equal to i)
corresponds to a wavelet signal decomposition of n levels.
If the wavelet signal decomposition is performed n1 levels. two blurred and a
total of N2 detail coefficients are obtained. The coefficients at the (nl) level of
decomposition are associated with scaling and wavelet functions which span half the
original signal and can be thought of as having the same frequency as thea2 and b2
Fourier coefficients with frequency equal to (2(2rrJN)). In other words, a Fourier
frequency index of 2 (equal to i) corresponds to a wavelet signal decomposition of nI
levels.
Likewise, if the wavelet signal decomposition is performed n2 levels, four
blurred and a total of (N4) detail coefficients would be obtained. The coefficients at
level (n2) can be visualized as having the same frequencies as the a4 and b4 Fourier
coefficients corresponding to frequency equal to (4(2rrJN)). In other words, a cutoff
frequency of 4 (equal to 22
) corresponds to a wavelet signal decomposition of n2 levels.
Similarly. for a wavelet signal decomposition of ni levels, i blurred and a total of
(N2i) detail coefficients would be obtained. The coefficients at the (ni) level can be
visualized as having the same frequencies as the Fourier coefficients a2i and b/ In other
words, a Fourier frequency index of i corresponds to a wavelet signal decomposition of
ni levels. The following table gives the relation between the cutoff frequency index and
the corresponding levels of wavelet decomposition.

80
Table 4.1.
The number of levels of wavelet decomposition for various cutoff indices
Frequency Corresponding level of
index decomposition
2°=1 n
2'=2 nl
22=4 n2
23=8. n3.
2; nl
2"'=NI2 0
From Table 4.1, we deduce that the minimum number of levels of decomposition,
min_levels, can be calculated from the dominant frequency index by the following
formula:
min_levels =logicutofffrequency index) (58)
From the above discussion, it is evident that for the number of levels of decomposition to
be an integer value, the cutoff frequency index needs to be a power of two. However,
the cutoff frequency obtained in practice mayor may not be a power of two. When not

81
an integer, the calculated minimum number of levels of wavelet decomposition is
incremented to the next integer. In doing so, a few of the low frequencies are di carded.
The number of levels of decomposition thus obtained min_levels, gives the
minimum number of levels of decomposition which retain all the dominant frequencies.
The smoothed signal thus obtained generally requires additional smoothing as discussed
previously.
The concept of a dominant frequency is further illustrated in Figure 4.4. A signal
with a strong periodic component is shown in Figure 4.4a. The cumulative power
spectrum obtained using the first 64 samples of the signal is presented in Figure 4.4b.
The cutoff frequency index is determined to be 8. The minimum number of levels of
decomposition is therefore 3. The original signal smoothed 3 levels using sixth order
Daubechies wavelets is shown in Figure 4.4c. From Figure 4.4c, it can be confirmed that
the smoothed signal still retains the periodic swings in the original signal. Figure 4.4d
verifies that decomposing the signal one level beyond this level yields a smoothed
approximation without the dominant low frequency oscillation.
4.3.4 Determination of the minimum window length required for process
monitoring
A minimum suggested window length for patternbased data analysis can be determined
from the cutoff frequency index by an empirical algorithm a described below. We
recommend that the minimum window length should be long enough to capture four
cycles of all the contain four cycles of all the frequencies below the cutoff, inclusive of
the cutoff.
82
200 400 600 800 1000
Original signal characterized with highly periodic swings.
..... .'
200 400 600 800 I000
Smoothed approximation obtained by decomposing the original
signal three levels retains the dominant frequencies as shown.
Signal needs to be smoothed beyond the dominant frequencies.
3 5 8 10 13 15 I8 20 23 25 28 30
Cumulative power spectrum obtained from the first 64 data points
of the original signal.
(b)
0.70
0.65
0.60
0.55
0.50
(c)
0.70
0.65
0.60
0.55
0.50
0.70
0.65
0.60
0.55
0.50
(a)
0.00125
0.00100
0.00075
0.00050
0.00025
200 400 600 800 I000
Cd) Smoothed approximation of the original signal obtained by
decomposing the original signal four levels using sixth order
Daubechies wavelets, and the original signal.
Figure 4.4: Representation of the dominant frequencies inherent in
a signal and the need to smooth beyond these frequencies.
83
This is an empirical recommendation and is intended to ensure that a sufficient portion of
the signal is used to determine the fundamental trend. The cutoff frequency index,
jieq_index, is first converted to its value in radians per unit time, using the sampling
interval of the data, T., and the length (number of data points) of the signal considered for
fast Fourier transform computation, N1, by the relation
This value of the minimum window size in time is displayed and the user prompted to
input the length of window of observation in time.
W =(freq_index) *(2rt)/(N I *T.~)
The value of W in cycles per unit time is given by
W =(freq_index)/N1*T.~
Then, the length of the minimum window in units of time is given by
min_win =4*(NI *T~)/(jieq_index)
(59)
(60)
(61 )
4.4 Automated Trend Extraction algorithm
Our proposed automated trend extraction (ATE) algorithm is described in this ection. A
now sheet that shows the various functions involved in the code, and the key variables is
included. The algorithm incorporates all the concepts described in this chapter.
FLOW OF CONTROL IN THE AUTOMATED TREND
EXTRACTION ALGORITHM
The raw sensor signal of length N, a
power of 2, that needs to be smoothed
is read from a data file and the raw
signal displayed. Here, N = 2n.
+
The raw sensor signal is input to the function di!)play.m where a steady part of
the raw sensor signal is chosen by the user for performing the fast Fourier
transform calculations. The length of this data window, winx needs to be a
power of 2. The length of winx is usually sleceted to be 64. For very noisy'
signals, a larger window size of 256 would be helpful. Once the length of winx
is selected, a window of length of winx is moved through the entire raw signal
so as to facilitate the user in selecting a steady part of the signal.
The data contained in winx is input to the function !asf.m, where the
Fourier coefficients of all the data points in winx are determined, the
I power at each data point calculated, and the cumulative power at each
point is obtained.
The cumulative power spectrum is used to determine the cutoff
frequency index. The cutoff frequency is a measure of the fundamental
frequency and other dominant frequencies. Frequencies beyond the
cutoff are not considered. The cutoff frequency index is detected by
the function detect.m.
2
84
d
85
The cutoff frequency index is input to the function window.m
to determine the minimum data window size required for
effective pattern recognition. The decision on the minimum
window size is based on the period of occurrence of all
frequencies till the cutoff. Four times the period is taken as a
judicious estimate of the length of the minimum window size
required for "effective pattern recognition".
The user is prompted to input a data window size
greater than the minimum window size, in the function
window.m. The size of this window, say x, is problem
specific and is left to the discretion of the user.
The user input window size is used to determine the
number of levels of wavelet signal decomposition, based
on an empirical algorithm given by
no levels2 = log2(N) log2(4+N*Tlx)
Minimum levels of
decomposition 1S
determined from the
cutoff frequency
index by the relation
no levels] =
log2(length (winx))log2 0)
no levels
Level ofsmoothing decided
by the Cumulative power
spectrum
Number of levels of wavelet signal decomposition
used for smoothing is given by
This is done in the function decide.m
no levels = no_levels1
..
Number of levels of wavelet signal decomposition
4 used for smoothing is given by
no levels = no levels2
This is done b the function decide.m.
86
The raw sensor signal is smoothed with the
3 number of levels obtained. Smoothing by wavelets
is done by the function power.m.
Smoothing less. Need to
get beyond the dominant
frequencies
Smooth one
more level
STOP
t Smoothing more than
desirable
Desmooth
one level
...
87
AUTOMATED TREND EXTRACTION CODE
START
di!)play.m
eruneh.m
no levels '>tI..
= no levelsl
Smoothing More
decide.m
no levels
deleet.m
j,no_levels]
window.m
min_win,no_levels2
no levels
no levels+
..

Chapter Summary
This chapter described the bases for an Automated Trend Extraction algorithm. The
algorithm utilizes information concerning the monitoring application (user specified
window size) and the characteristics of the signal to be smoothed to generate a
recommendation for the appropriate degree of wavelet smoothing. Application of the
ATE algorithm is illustrated in the next chapter.
88
..

Chapter 5
CASE STUDIES
5.1 Introduction
To evaluate the performance of the ATE algorithm, a set of signals with different
characteristics and noise content were considered. Plant data from temperature, pressure
and flow sensors was used. In every case, a window length of 45 minutes was selected
for process monitoring. Table 5.1 summarises the results obtained by the automated trend
extraction algorithm for all case studies. The foJ lowing section presents an analysis of all
test cases considered.
89
Table 5. t
Summary of the results obtained by the ATE code
90
Case Figure Signal Tf in freq. min_levels Window I acClevels
length m1OS. index length in mins.
. .
mInImum actual
I 5.1 1024 I 8 3 32 45 4
2 5.2 1024 1 ] 1 2.54 24 45 4
3 5.3 1024 I 9 2.83 29 45 4
3a 5.4 512 ) 8 3 32 45 4 I
4 5.5 8192 1/6 5 3.67 9 45 7
4a 5.6 2048 2/3 11 2.54 ]6 45 5 I
I 4b 5.7 4096 ]/6 6 3.42 8 45 7
4c 5.8 2048 1/6 6 3.68 9 45 7
5 5.9 8192 1/6 15 2.09 3 45 7
5a 5.10 4096 1/6 6 3.42 8 45 7
5b 5.11 1024 1/6 6 3.42 9 45 7
•
..
91
5.2 Individual Case Studies
5.2.1 Case IFigure 5.1: The signal considered in this case is characterized by periodic
swings of high frequency. The minimum number of levels of wavelet signal
decomposition obtained from the cutoff frequency index, min_levels, is 3. The number
of levels of wavelet signal decomposition determined from the user input process
monitoring window length, act_levels, is 3.49. As the number of levels of decomposition
needs to be an integer, the actual number of levels of decomposition, acelevels is
incremented to the next integer value 4. The smoothed signal representation obtained by
smoothing four levels using sixth order Daubechies wavelets is shown in Figure 5.1 d.
5.2.2 Case 2Figure 5.2: A signal with different characteristics is considered in this
case. Although this signal also appears to contain a highly periodic component, the
amplitude is much smaller and the frequency BlUCh higher than Case ~. The minimum
number of levels of wavelet signal decomposition obtained from the cutoff frequency
index, min_levels, is 2.54. The number of levels of wavelet signal decomposition
determined from the user input window size, acelevels. is 3.49, and is incremented to 4.
The smoothed signal is shown in Figure 5.2d. The smoothed signal approximation
clearly captures the true trend of the signal.
5.2.3 Case 3Figure 5.3: In this case, a signal characterized by moderately large, abrupt
changes in its fundamental trend is considered. The number of levels of wavelet signal
decomposition obtained from the cutoff frequency index, min_levels, is 2.83. The
0.70
0.65
0.60
0.55
0.50
92
10 20
0.64
0.63
0.62
0.61
0.60
0.00125
0.00100
0.00075
0.00050
0.00025
200
T
I
I
I
400
30
(a)
40
(b)
600
50
800
60 70
1000
80
0.70
0.65
0.60
0.55
0.50
3 5 8 10 13 15 18 20 23 25 28 30
(c)
200 400 600 800 1000
(d)
Figure 5.1. (a) original signal with periodic swings (b) First 64 data points of the
original signal for obtaining the CPS (c) Cumulative power spectrum
(d) Smoothed representation of the original signal using sixth order Daubechies
and decomposing four levels.

93
0.40
0.35
0.30
0.25
0.20
200 400 600 800 1000
(a) ..)
0.36 !~~
0.34 )I;'"
:[~
0.32
...0..;.) !:,.
!,"
0.30 I."~J
0.28 .;.~i
10 20 30
J~
40 50 60 .of
(b) .,
~.~ ~
"i ,~ ::f
0.0006 ./ ::1 '0:
.~
0.0004 d~
0.0002
3 5 8 10 13 15 I8 20 23 25 28 30
(c)
0.40
0.35
0.30
0.25
0.20
200 400 600 800 1000
(d)
Figure 5.2. (a) original signal (b) First 64 data points of the original signal for obtaining the
cumulative power spectrum (c) cumulative power spectrum (d) smoothed representation of
the original signal obtained by smoothing four levels, using sixth order Daubechies.
0.60
0.50
0040
0.30
0.20
94
0048
0046
0044
0042
0040
0.00100
0.00075
0.00050
0.00025
10
200
(
'i
20
400
(a)
30
(b)
600
40
800
50
1000
60
3 5 8 10 13 15 I8 20 23 25 28 30
(c)
0.60
0.50
0040
0.30
0.20
200 400
Cd)
600 800 1000
Figure 5.3. (a) original signal (b) First 64 data points of the original signal for for
obtaining the cumulative power spectrum (c) cumulative power spectrum
(d) Smoothed representation of the original signal, smoothed four levels using sixth
order Daubechies wavelets.

95
number of levels of wavelet signal decomposition determined from the user input window
size, acelevels, is 3.49, and is incremented to 4. The smoothed signal representation
obtained by smoothing four levels using sixth order Daubechies wavelets is shown in
Figure S.3d. Excellent performance is again noted.
5.2.3a Case3aFigure 5.4: The second half of the signal in Figure S.3a is considered
separately to verify that the smoothing algorithm is independent of signal length. The
algorithm is expected to provide the same degree of smoothing as in Figure S.3d as the
noise content in both these signals is the same. The calculated minimum number of
levels of wavelet signal decomposition obtained from the cutoff frequency index,
min_levels, is 3.0. The number of levels of wavelet signal decomposition detennined
from the user input window size, acelevels is 3.49, and is incremented to 4. The
smoothed signal representation is shown in Figure 5.4d. To compare the smoothing
obtained in Figure S.4d with that of Figure 5.3d, the second half of the original and
smoothed :ignals in S.3d is shown in Figure SAe. From figures SAd and 5.4 e, it is
observed that the degree of smoothing obtained is the same in both the cases, and the
smoothed patterns are identical.
5.2.4 Case 4Figure 5.5: The signal considered in this case is sampled at a rate (Ts =
0.167 min.) which is six times higher than the sampling rate of the signals considered to
this point (Ts = 1 min.). As a consequence, the high frequency content of the signal is
observed to be very prominent. Further, the signal is also characterized by sharp changes
.,e:
..'

96
0.60
0.50
0040
0.30
0.20
100 200 300 400 500
(a) .• )
0.38 ~~ ~
0.38 ...
.i~
0.37 :fA.. ,"
0.37 ~~
0.36 ~i
20 40 60 ~.
(b) , :; ~
~ ~
11:
o.cx:mo 'F; • ;~i 0.00020 .~
a.crow
0.00000
} 5 8 10 13 15 18 20 23 25 28 30
(c)
0.60
0.50
0040
0.30
0.20
100 200 300 400 500
(d)
Figure 504
...

0.60
0.50
0.40
0.30
0.20
97
513 613 713
(e)
813 913 1013
Figure 5.4. (a) original signal, last 512 samples of the signal in Figure 5.3 (a)
(b) First 64 data points of tthe signal in (a) for obtaining the cumulative
power spectrum (c) cumulative power spectrum (d) smoothed representation
of the signal in (a), smoothed four levels, using sixth order Daubechies wavelets
(e) last 512 samples of the original and smoothed signal in Figure 5.3(d).
•
98
0.60
0.58
0.55
0.53
0.50
2000 4000 6000 8000
(a)
0.55
)
~S
0.55 .1'0
.~
0.54
!...,
:'"
0.54
, ..
,)
0.53
:1
§I
10 20 30 40 50 60 :..
~
(b) ~,
:~
0.000065 ~~ 0.000060 'I J
0.000055 I
,3
I
0.000050 I
0.000045 I
I
0.000040 I
I
0.000035 I
3 5 8 10 13 15 18 20 23 25 28 30
(c)
0.60
0.58
0.55
0.53
0.50
2000 4000 6000 8000
(d)
Figure 5.5. (a) original signal (b) First 64 data points of the original signal taken
for obtaining the cumulative power spectrum (c) cumulative power spectrum
(d) Smoothed representation of the original signal, smoothed seven levels
using sixth order Daubechies wavelets.
•
99
in its fundamental trend. As the sampling interval is much smaller compared to the
previous cases, the number of levels of decomposition in this ca e should be relatively
high. The number of levels of wavelet signal decomposition obtained from the cutoff
frequency index, min_levels, is 3.67. The number of levels of wavelet signal
decomposition determined from the user input window size, acClevels, is 6.08, and is
incremented to 7. The smoothed signal representation obtained by smoothing seven
levels using sixth order Daubechies wavelets is shown in Figure 5.5d. Excellent
performance is again noted.
5.2.4a Case 4aFigure 5.6: For this case study, the previous signal was sampled at a
rate four times slower (To = 0.667 min. as compared to 0.167 min.). The purpose of this
case study is to demonstrate that smoothing performance is dependent on the sampling
rate. As compared to the signal in Figure 5.5a, the high frequency content of signal in
Figure 5.6a is reduced due to the longer sampling interval (0.667 min. against 0.167
min.). Consequently, the number of levels of decomposition based on the process
monitoring window should be two less than in Case 4.
The minimum number of levels of decomposition obtained from the cutoff
frequency index is 2.54. The number of levels of decomposition obtained from the user
input window length is 4.08 as expected. Figure 5.5d shows the original and smoothed
signals. The degree of smoothing is the same as produced in Case 4. This case study
demonstrates that the automated trend extraction code is adaptive and achieves a degree
of smoothing that is dependent only on the length of the pattern recognition window or
;'"
[~
~::>
I"I"
,"
, :>
.;.~ ;J ::,
Figure 5.6. (a) original signal, signal sampled every 40 seconds (b) First 64
data points of the signal in (a) for obtaining the cumulative power spectrum
(c) cumulative power spectrum (d) Smoothed representation of the original
signal, smoothed five levels using sixth order Daubechies wavelets.
101
the frequency content of the signal, irrespective of the length of the signal.
5.2.4b Case 4bFigure 5.7: For this case study, the second half of the signal in Figure
5.5a is considered. The purpose of this case is to again demonstrate that smoothing
performance is independent of the length of the signal. The noise content of this signal is
the same as that in the signal shown in Figure 5.5a; therefore the same degree of
smoothing is anticipated. The number of levels of decomposition obtained from the cutoff
frequency index, min_levels, is 3.42. The number of levels of decomposition from the
process monitoring window size, acelevels, is 6.08, and is incremented to 7 levels. The
signal in Figure 5.7a smoothed seven levels using sixth order Daubechies wavelets is
shown in Figure 5.7d. To compare the smoothing obtained in Figure 5.7d with that of
Figure 5.5d, the second half of the original and smoothed signals in 5.5d are shown in
Figure 5.7e. From figures 5.7d and 5.7e, it is observed that the degree of smoothing
obtained is the same in both the cases, and the smoothed patterns are identical.
5.2.4c Case 4cFigure 5.8: A signal consisting of the last 2048 samples of the signal in
Figure 5.5a is considered to demonstrate again that the smoothing algorithm is
independent of the length of the signal. The number of levels of decomposition obtained
from the cutoff frequency index, min_levels, is 3.68. The number of levels of
decomposition obtained from the user input process monitoring window size is 6.08, and
is incremented to 7. The signal in Figure 5a smoothed seven levels using sixth order
Daubechies wavelets is shown in Figure 5.8d. The last 2048 samples of the original and
smoothed signals in Figure 5.5d are shown in Figure 5.8e for comparing the degree of
·, : ~ ·.·
; •·.~ ••
·." ·;3. ; l
; f· ·.··:
c
0.60
0.58
0.56
0.54
102
0.52 l. , , . ...
.,
'04
:~
...,.. ,...
J'"
.~
••
•
103
0.60
0.58
0.56
0.54
0.52
4097 5097 6097
(e)
7097 8097
.:
Figure 5.7. (a) original signal, 4096 samples, second half of the signal in Figure 5.5(a)
(b) First 64 data points of the signal in (a) for obtaining the cumulative power spectrum
(c) cumulative power spectrum Cd) Smoothed representation of the signal in (a), smoothed
seven levels using sixth order Daubechies wavelets (e) second half of the original and
smoothed signal in Figure 5.5(d).
,
1
,I
I
104
0.60
0.58
0.56
0.54
0.52
500 1000 1500 2000
(a)
0.55
;J '.~
0.55
,.... •04
0.54 ;..
I
0.54 }
:f
0.54 .~
1
20 40 60 ,~
(b) ,1
:l
0.00003
:i
:1
:'J.
0.00003
0.00002
0.00002
3 5 8 10 13 15 18 20 23 25 28 30
(c)
0.60
0.58
0.56
0.54
0.52
500 1000 1500 2000
(d)
Figure 5.8
s
0.60
0.58
0.56
0.54
0.52
6145 6645 7145
(e)
7645 8145
105
Figure 5.8. (a) original signal, last 2048 samples of signal in Figure 5.5 (a)
(b) First 64 samples of signal in (a) for obtaining the cumulative power spectrum
(c) cumulative power spectrum (d) smoothed representation of the signal in (a),
smoothed seven levels using sixth order Daubechies wavelet (e) last 2048
samples of the original and smoothed signal in Figure 5.5(d).
,..... ....
;...
I~
.}
:f
,1
'I
",I
'~
:'J1
...J..
106
smoothing obtained in Figure 5.8d. From both these figures it can be observed that the
degree of smoothing obtained is the same in both cases, and the smoothed patterns are
identical in nature. The above two case studies demonstrate that the degree of smoothing
obtained is the same for signals of the same frequency content, even though they may be
of different lengths, when the process monitoring window length employed is the same.
5.3 Discussion
In all the test cases, the smoothing obtained is as desired and the smoothed representation
preserves the fundamental trend of the original signal. All the above cases show clearly
that the degree of smoothing obtained is dependent on the pattern recognition window
length and is independent of the length of the original signal, for a fixed sampling rate.
We note, however, that the length of the signal should be sufficiently large so as to
provide more available levels of decomposition than required to smooth the signal.
The accuracy of the cutoff frequency index is critical in determining the
minimum level of smoothing, as the actual window length used for process monitoring is
indirectly dependent on the value of the cutoff frequency obtained. In all the above
cases, the cumulative power spectrum obtained is well defined and the cutoff frequency
index obtained is accurate. However, if the cumulative power spectrum is not well
defined, determination of the cutoff frequency is difficult and the degree of smoothing
obtained would depend entirely dependent on the discretion of the user in selecting the
window length for process monitoring. The following case presents such an example.
,"..J
.\
107
5.3.1 Case 5: The signal considered in this case is characterized by a large amount of
high frequency content distributed unevenly, as can be seen from Figure 5.9a. The
cumulative power spectrum obtained is shown in Figure 5.9c. It can be observed that the
cumulative power spectrum does not increase steadily and level off. It increases in smaH
steps. Consequently, a single striking bend is not observed as in all the previous cases.
In such cases where the cumulative power spectrum continues to increase steadily without
leveling off, determining the cutoff frequency index accurately becomes extremely
difficult. In this case, the cutoff frequency index is determined to be 15. The number of
levels of wavelet signal decomposition obtained from the cutoff frequency index
min_levels, is 2.09.
The accuracy of the minimum window length required for process monitoring
depends on how accurately the cutoff frequency index is determined. In such cases, the
minimum window size required may be lesser or greater than the actual window size the
user expects to input based on his experience. If the minimum window size required is
less than the window size anticipated by the user, and if the user input window size is
judiciously selected, appropriate smoothing is expected. However, if the minimum
window size required is greater than the window size anticipated by the user, the user will
have to compromise and input a larger window size than really desired. This may bring
about greater smoothing than desired (refer equation 57).
In this specific case, a user input window size of 45 minutes is greater than the
minimum obtained from the ATE code, so the degree of smoothing obtained is as desired.
The number of levels of wavelet signal decomposition determined from the user input
window size, ace/eve/s, is 6.08, which is incremented to 7. The smoothed signal
108
0.60
0.50
0.40
0.30
0.20
2000 4000 6000 8000
(a)
0.50
0.48
0.45
0.43
0.40
20 40 60
(b)
0.004
0.003 "
1/
I 0.002 I
I
0.001 I
I
I
5 lO 15 20 25 30
(c)
0.60
0.50
0.40
0.30
0.20
2000 4000 6000 8000
(d)
Figure 5.9. (a) original signal (b) First 64 data points of the original signal taken for obtaining
the cumulative power spectrum (c) cumulative power specctrum (d) Smoothed representation
of the original signal, smoothed seven levels using sixth order Daubechies wavelets.
109
representation obtained by smoothing seven levels using sixth order Daubechies wavelets
is shown in Figure 5.9d. From this figure, it can be observed that the degree of
smoothing is as desired. However the desirable degree of smoothing obtained is entirely
dependent on the process monitoring window size as the cutoff frequency index obtained
is not accurate.
The presence of large quantity of noise distributed unevenly does not yield a well
behaved cumulative power spectrum. In such cases, using a larger window size for the
fast Fourier transform may yield a cumulative power spectrum that increases steeply
initially and then levels off.
Figure 5. lOa shows the first 128 data points of the signal in Figure 5.9a. The
associated cumulative power spectrum is shown in Figure 5.lOb. Figure 5.1 Dc and 5.1 Od
show the first 256 data points of the original signal in Figure 5.9a, and the resulting
cumulative power spectrum obtained. Tn both cases, the cumulative power spectrum
obtained is not as desired. However, when a data window of the first 512 points of the
signal in Figure 5.9a is used for obtaining the cumulative power spectrum, the cumulative
power spectrum obtained is defined much better as seen in Figure 5.1 Of. The cutoff
frequency index obtained is 15.
The second set of 512 data points of the original signal is considered and the
resulting cumulative power spectrum is shown in Figure 5.1 Oh. The cutoff frequency
obtained from this cumulative power spectrum is 19. All these results show that for a
signal with nonuniform noise, selecting a relatively larger window size for performing
the fast Fourier transform results in a well behaved cumulative power spectrum.
0.50
0.45
0.40
0.35
110
0.30 L._~'_' r ~
25 50
(a)
75 100 125
0.05
0.04
0.03
0.02
0.01
0.50
0.45
0.40
0.35
0.30
5 10 15 20 25 30 35 40 45 50 55 60
(b)
50 100
(c)
ISO 200 250
0.03
0.02
0.02
0.01
0.01
10 20 30 40 50 60 70 80 90 100 110 120
Cd)
Figure 5.10. (a) First 128 data points of the signal in Figure 6.9a
(b) cumulative power spectrum for the data in (a) (c) First 256 data
points of the signal in Figure 6.9a (d) cumulative power spectrum
for the data in (c).
1I 1
0.45
0.40
0.35
100 200
(e)
300 400 500
0.10
0.08
0.05
0.03
20 40 60 80 100 120 140 ]60 ]80 200 220 240
(f)
0.50
0.45
0.40
0.35
0.30
100 200 300 400 500
(g)
0.150
0.125
0.100
0.075
0.050
0.025
20 40 60 80 100 120 140 ]60 180 200 220 240
(h)
Figure 5.10. (contd.) (e) First 5]2 data points of signal in Figure 6.9a.
(f) cumulative power spectrum for data in (e) (g) Second set of 512 data
points of the signal in Figure 6.9a (h) cumulative power spectrum for data in (g).
112
5.4 Chapter Summary
Various signals with different characteristics have been studied to evaluate the
performance of the automated trend extraction algorithm. The important observations are
listed below:
(I) The smoothed representation obtained is normally determined by the user
input window length for process monitoring. The number of levels of decomposition
obtained from the cutoff frequency, min_levels, only gives the minimum number of
levels to be decomposed so as to retain the dominant frequencies present in the original
signal. To obtain a smoothed representation of the original signal, smoothing should be
performed beyond the dominant frequencies.
(2) The window length for process monitoring is determined by the user based on
experience with the problem at hand. The user input window length determines the
number of levels of decomposition beyond the dominant frequencies and results in the
desired degree of smoothing. The minimum window length required for process
monitoring is directly dependent on how accurately the cutoff frequency index is
determi ned.
(3) Accurate determination of the cutoff frequency index depends on the
behavior of the cumulative power spectrum. If the cumulative power spectrum is not well
defined, the cutoff frequency index determined may not be accurate. In such a situation,
the smoothed representation obtained depends totally on the discretion of the user in
selecting the window length for process monitoring.
(4) If the cumulati ve power spectrum is well defined, the required degree of
smoothing is the same for signals with the same frequency content, irrespective of the
length of the signal, for the same process monitoring window length.
L13

Chapter 6
CONCLUSIONS
Wavelets provide a much better alternative for trend extraction than conventional
methods of signal processing such as the direct method, digital filters, and the Fourier
transform. Wavelets possess the ability to extract essential trends from process ignals
thus helping to provide compact representation which is imperative for efficient realtime
patternbased monitoring and control.
This work presented an automated approach to obtain the desired degree of
smoothing as required for realtime patternbased monitoring applications. The
properties of the wavelet and fast Fourier transforms are exploited to achieve this
purpose. Our automated trend extraction system can be used as a standalone system
which helps operators create their own process monitoring applications.
The automatic trend extraction algorithm automatically recommends an
appropriate degree of smoothing by utilizing information concerning the monitoring
application (user specified window size) and the characteristics of the signal to be
smoothed. The minimum degree of smoothing is determined from the characteristics of
the signal and is dependent on how accurately the cutoff frequency index i determined.
The cutoff frequency is determined from the cumulative power spectrum, which is
determined by performing the fast Fourier transform operation on a selected part of the
114
115
original signal. In some cases, when the cumulative power spectrum is not be well
defined, the minimum window size determined would be inaccurate may even be greater
than the actual window size anticipated by the user. In such cases, the user wiJI have to
compromise and input a larger window size resulting in less smoothing than desired.
Performance of the automated trend extraction algorithm is independent of signal
length for a fixed sampljng rate. Consequently, the algorithm is suitable for widespread
application without constraint.
Recommendations and future work
The approach adopted for most of the work done is basically empirical in nature. This
work needs to be consolidated with a theoretical background. Most of the relations and
parameters used in this work are based purely on experience and knowledge of sensor
signal behavior. This empirical work needs to be supported by a more generalized
mathematical basis. Following are some of the recommendations
• The relations used to convert the cutoff frequency index and the user input
window size to the number of levels of wavelet signal decomposition are empirical.
These relations need to be given a strong theo 