APPLICATION OF THE WIGNER DISTRIBUTION TO PROBLEMS IN TIME
VARYING SIGNAL ANALYSIS
By
PREETI RAO
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1990
ACKNOWLEDGMENTS
I would like to express my heartfelt thanks to my advisor, Dr. Fred Taylor, for
introducing me to a fine research topic and for his continuous encouragement, sup
port and kindness over the past three years.
I would also like to sincerely thank Dr. Childers for his helpful criticisms and
for generously allowing me access to his excellent library and other laboratory facili
ties.
I am very grateful to Dr. Principe for his advice and for maintaining a sincere
interest in my progress all along.
I further thank Dr. Childers, Dr. Principe, Dr. Latchman and Dr. Y.C. Chow for
serving on my committee.
I also wish to thank Dr. Tlusty and his graduate students of the Machine Tool
Lab for their time and efforts spent in collecting data for my use.
Finally, I take this opportunity to thank the many people I have come to know in
Gainesville as a student here, who have made this one of the most pleasant and
fruitful periods in my life.
 ii 
TABLE OF CONTENTS
ACKNOW LEDGMENTS ................. ................................... ii
ABSTRACT ......................................................... v
CHAPTERS
1 INTRODUCTION ................................................ 1
2 A REVIEW OF TIMEFREQUENCY DISTRIBUTIONS ......
2.1 The Need for a TimeFrequency Representation
2.2 A Brief Historical Review ...................
2.3 A Summary of Work to Date ................
2.4 Cohen's Class of Bilinear Representations .....
2.5 The Wigner Distribution ....................
2.5.1 Desirable Properties .................
2.5.2 Interference Terms ...................
2.6 The Spectrogram and Its Relation to the WD ...
3 ESTIMATION OF INSTANTANEOUS FREQUENCY
RATE USING THE WIGNER DISTRIBUTION ......
. . 5
. . . 5
.................. 6
.................. 8
. . . 9
................. 11
.............. ... 11
................. 15
. . . 20
AND FREQUENCY
. . . 24
3.1 The Concept of Instantaneous Frequency .......................
3.2 Estimation of Instantaneous Frequency from the WD .............
3.2.1 Performance of the DWD Peak for Linear F.M. Signals .....
3.2.2 Computer Simulations .................................
3.2.3 Estimation of Instantaneous Frequency of Nonlinear F.M.
Signals ..... ......................................
3.3 Joint Estimation of Frequency and FrequencyRate ..............
3.3.1 Derivation of a New Kernel ......... ...................
3.3.2 Computer Simulations ............................
 iii 
 iv 
4 DETECTION AND LOCALIZATION OF NARROWBAND TRANSIENT
SIGNALS ....................................................... 41
4.1 Detection of Transient Signals ............................... 41
4.2 Existing Approaches to the Detection and Localization of Narrowband
Transient Signals ...................... ...... ................. 42
4.3 Desirable Detector Characteristics ............................. 44
4.3.1 High TimeFrequency Resolution ........................ 46
4.3.2 Robustness to Transient Duration ........................ 51
4.4 A Detector based on Instantaneous Power ...................... 62
4.5 Noise Performance of the WD ............................... 67
4.6 Noise Performance of the WD Transient Detector ................ 75
4.7 Analysis of Experimental Milling Machine Data ................. 80
5 APPLICATION OF THE WIGNER DISTRIBUTION TO VIBRATION
SIGNAL PROCESSING .......................................... 87
5.1 Spectral Analysis in Vibration Signal Processing ................. 87
5.2 Envelope Detection ......................................... 88
5.2.1 Vibration Monitoring in the Gas Turbine ................ 90
5.2.2 Algorithm and Simulations ............................. ..91
5.3 Detection of Periodic Transients .............................. 99
6 SUMMARY ................................................... 104
APPENDIX ....................................................... 108
REFERENCES .................................................... 112
BIOGRAPHICAL SKETCH ......................................... 117
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
APPLICATION OF THE WIGNER DISTRIBUTION TO PROBLEMS IN TIME
VARYING SIGNAL ANALYSIS
By
Preeti Rao
December 1990
Chairman: Dr. Fred J. Taylor
Major Department: Electrical Engineering
The Wigner distribution (WD) is a bilinear timefrequency signal representa
tion that possesses a number of characteristics desirable in timevarying signal anal
ysis. This dissertation investigates and extends the application of the WD and its
derivatives to new areas.
The WD is placed in the general framework of nonstationary signal analysis and
timefrequency distributions. The WD has been previously applied to the estimation
of instantaneous frequency of phasemodulated signals. Analytical expressions for
the performance of the discretetime WD (DWD) in estimating the instantaneous
frequency of linear f.m. signals in additive white Gaussian noise are derived and are
verified using simulation. It is shown that the peak of the DWD provides an optimal
estimate of the instantaneous frequency for such signals at sufficiently high signal
tonoise ratios. The applicability of these results to the general case of nonlinear
f.m. signals is discussed. Based on a methodology suggested by the WD, a new
signal kernel is derived whose spectrum is concentrated at the frequencyrate of the
arbitrary phasemodulated signal. A computationally simple implementation of the
frequencyrate estimator is proposed and computer simulations are used to investi
gate the noise performance of this method.
The WD is applied to the detection of narrowband transient signals of unknown
waveform and arrival time in additive noise. A specific class of received signals is
considered, motivated by examples from underwater acoustics and vibration signal
processing. A detector based on the WD is proposed that provides good time local
ization of the.transient waveform and is not strongly sensitive to changes in the pulse
duration or frequency of the received transient signal. The performance of the detec
tor in additive Gaussian noise is analyzed.
The field of vibration signal processing is rich with potential applications for
timefrequency methods. The advantage provided by the WD over the traditional
shorttime spectral magnitude is demonstrated by examples in the vibration monitor
ing of machinery, where nonstationary signals arise due to speed variations and other
timevarying phenomena. Specifically, we consider the problems of envelope detect
ing timevarying narrowband signals, and the detection of periodic transients of low
energy from the sidebands of a high frequency, timevarying carrier component.
 vi 
CHAPTER 1
INTRODUCTION
Signal analysis tools for nonstationary, or timevarying signals, need to incorpo
rate the time dependence of signal characteristics in their formulation. Most often
the signal feature chosen to be represented is frequency, giving rise to timefrequen
cy representations. Alternative signal features such as the timeextent of elementary
signal components that comprise the input signal may also be of interest, giving rise
to timescale methods.
The shorttime Fourier transform (STFT) and the wavelet transform (WT) are
linear signal transforms that analyze the signal in terms of elementary waveforms
indexed by time and another parameter, which is frequency in the case of the STFT,
and timescale in the case of the WT. When the scale parameter in the WT is chosen
inversely proportional to frequency, we get a timefrequency representation similar
to the STFT, but characterized by timefrequency resolution that is frequency depen
dent, instead of fixed. When one is interested in the energy distribution of a signal in
time and frequency, it is natural to consider transforms that are characterized by a
bilinear dependence on the signal. If the further constraint of invariance with re
spect to shifts in both time and frequency is applied, the admissible transforms re
duce to those which belong to the "Cohen's class." Similarly if the transform is
constrained to be invariant with respect to shifts and dilations in time, we get the
affine distributions [Fla89a].
Although the question of a joint signal representation in time and frequency was
first addressed as early as 1946 by Gabor and soon after by Ville, it is only in the past
decade that there has been a surge of activity in the application of timefrequency
1
2
distributions to problems in signal processing. A number of properties that are desir
able in a timefrequency signal representation have been identified. These are gen
erally based on properties recognized in the spectral density functions of stationary
signal analysis with the addition of an explicit time dependence. It turns out that
there is no single distribution that satisfies all the requirements. Much research has
been addressed toward obtaining distributions that satisfy specific properties and
their interpretation. Various theoretical aspects of these distributions have been
studied and a number of practical applications proposed.
The Wigner distribution (WD), a member of Cohen's class of bilinear timefre
quency representations, possesses a large number of desirable properties. Of the
representations that satisfy the finite support property and preserve inner products,
the WD provides the most concentration in the timefrequency plane for a given
signal. The WD also satisfies the marginal and local properties in time and frequen
cy. These characteristics have made it an attractive tool in timevarying signal analy
sis.
A primary objective of this work is to extend the existing applications of time
frequency methods, based on the WD and its derivatives, to new areas involving
timevarying signals and systems. Some of the applications considered here have
traditionally benefitted from shorttime spectral methods. The performance of WD
based methods is investigated and compared with available techniques.
A review of timefrequency distributions is presented in Chapter 2 and an at
tempt is made to place the WD in the perspective of timefrequency signal represen
tations and spectral analysis methods. The properties of the WD that form the bases
for the applications considered in this work are discussed in depth. The WD of
multicomponent signals is characterized by interference terms that obscure its inter
pretation as a signal energy distribution. The theory of interference in the WD is
given and the methods available for its suppression are described.
3
In Chapter 3 we address the problem of the joint estimation of frequency and
frequencyrate of arbitrary phasemodulated signals, one that arises in situations
where the velocity of an object and the curvature of its path must be measured simul
taneously. The ability of the WD to provide a concentrated distribution at the instan
taneous frequency of a narrowband signal has been exploited widely in practice. We
examine analytically the performance of the peak of the discretetime signal WD in
the estimation of the instantaneous frequency of f.m. signals in additive noise. An
exact expression for the variance of the estimate for the case of a linear f.m. signal in,
additive white Gaussian noise is derived and verified using computer simulations.
The extension of this result to the more general case of nonlinear f.m. signals is
discussed. Based on a methodology suggested by the WD, a new signal transform is
obtained by which the frequencyrate of an arbitrary f.m. signal can be estimated. A
computationally simple implementation of the frequencyrate estimator is proposed
and computer simulations are used to investigate its performance in noise.
Chapter 4 investigates the application of the WD to the detection of unknown
transient signals in additive noise. Since the received signal is inherently nonstation
ary, it is expected to benefit from the application of timefrequency methods. Moti
vated by examples from underwater acoustics and vibration signal processing, we
consider transient signals that are monochromatic and of narrow time durations rela
tive to the observation interval. Further the transient waveform is assumed to be
embedded in a noisy background comprised of deterministic and random compo
nents. We seek to detect the transient signal and localize it in time. The viability of
the spectrogram (which is closely related to existing methods for such problems) and
the WD in achieving this are explored and compared. The smoothedpseudoWD
allows the independent control of time and frequencyresolution in the signal repre
sentation, and hence provides a more flexible tool in applications to nonstationary
situations than the fixedwindow spectrogram. A detector based on the WD is pres
4
ented that provides good timefrequency resolution and is relatively insensitive to the
time duration of the transient signal. The noise properties of the WD estimate are
discussed and the performance of the detector for transient signals in additive Gaus
sian noise is derived. The detection algorithm is applied to a sample of experimental
ly obtained milling machine data containing transient vibration due to tool chatter.
The spectral analysis of vibration signals is an integral component of machinery
monitoring and fault detection. It is typically carried out using the magnitude of the
shorttime Fourier transform. While this method provides satisfactory results when
the signal is stationary over the duration of the analysis window, there are situations
when the vibration signal is nonstationary due to the timevarying nature of the
source of vibration. In Chapter 5 we demonstrate, by examples, the advantage pro
vided by WDbased methods in the monitoring of nonstationary vibration signals.
Specifically, we investigate the application of the WD to the problems of envelope
detection of narrowband, timevarying signals, and the detection of periodic tran
sients in timevarying signals.
CHAPTER 2
A REVIEW OF TIMEFREQUENCY DISTRIBUTIONS
2.1 The Need For A TimeFrequency Representation
The Fourier representation of a signal enables the decomposition of the signal
into its individual frequency components and yields the energy density of each of
these components. Standard spectral estimation methods are based on the implicit
assumption that the signals being analyzed are characterized by a stationary spec
trum or a spectrum that does not vary with time. In a great number of signal process
ing applications such an assumption is not justified. The spectral characteristics of
signals arising in nature often evolve in time due to some changing characteristic of
the underlying process. Examples of such timevarying signals are speech, seismic
signals and biological signals. The phasemodulated signals of radar and sonar, as
well as the acoustic signals from rotating machinery, also belong in this class. The
analysis of such timevarying signals cannot be carried out entirely using standard
Fourier techniques or other spectral estimators based on the stationarity assumption.
For while these methods give the spectral components of the signal, the time localiza
tion of the frequency components is often concealed in the values of the different
phases, making it difficult to access this important element in the investigation of
such signals..
For the past several decades the common tool for timevarying spectral analy
ses has been the shorttime Fourier transform, which is implemented using either a
bank of timewindows or a bank of bandpass filters. The former is more common
and is based on decomposing the signal into short, and possibly overlapping, time
5
6
segments in each of which the signal is assumed to be approximately stationary, and
computing the Fourier spectrum of each segment to estimate the frequency content
of the signal at that time. There are many situations in which such an analysis is
neither adequate nor appropriate. This can occur when the signal spectrum is chang
ing so rapidly that it is difficult to find a suitable shorttime window. Any shorttime
spectral energy measurement requires a minimum effective bandwidth and observa
tion time, and frequency shifts within this bandwidth as well as power variations
within this observation time remain obscured. In the case of the shorttime spec
trum, however, the amount of timefrequency tradeoff is fixed by the arbitrary,
often predetermined, selection of windows/filters and the reciprocal relation between
spectral and temporal resolutions. These considerations have given rise to the need
to turn to the practical application of joint functions of time and frequency, historical
ly termed "timefrequency distributions" in order to describe the energy density of a
signal in time and frequency simultaneously.
In this chapter we give a historical review of timefrequency representations and
summarize briefly the research to date. We discuss in greater depth the Wigner
Distribution (WD), whose properties and applications are the subject of this thesis.
We also discuss the spectrogram and its relation to the WD.
2.2 A Brief Historical Review
The question of a joint timefrequency distribution was first addressed by Ga
bor [Gab46] and Ville [Vil48], inspired by similar developments in quantum me
chanics, where there is a partial mathematical resemblance to timefrequency analy
sis. Gabor proposed an expansion of a given waveform into the sum of elementary
signals of "minimum" spread in time and frequency each with different center fre
quency and epoch. Since the elementary signals are chosen to be concentrated in
time and frequency, a local measure of timefrequency signal energy is obtained.
7
Helstrom [Hel66] made Gabor's representation exact by using an integral expansion.
The elementary signals of Gabor's representation overlap, and hence the kernel does
not represent the energy distribution in time and frequency. Ville [Vil48] attempted
to derive an exact signal energy distribution, based on one presented by Wigner in
1932 to calculate the twodimensional probability distribution in time and frequency
of a quantum mechanical particle, defining what is now known as the WignerVille
distribution, or simply, the Wigner distribution. Its characteristic function, intro
duced in modified form to radar theory by Woodward, is well known as the ambiguity
function. Page [Pag52] defined an "instantaneous power spectrum" as the rate of
change of the energy density spectrum of the signal segment from oo to T, as T is
increased. Levin [Lev64] used the same definition for the segment (T, oo) and de
fined a new function as the average of both types of instantaneous power spectra. In
an insightful paper, Rihaczek [Rih68] derived a complex energy density function
based on the concepts of the onedimensional energy density functions in time and
frequency. Each of the distributions enumerated here has been derived independent
ly from plausible timefrequency concepts and each has its own set of distinct prop
erties that make it a good candidate for timefrequency analyses. It was realized by
Cohen [Coh66] that these timefrequency representations have more in common
than is readily apparent, and he proposed what has now become known as the Co
hen's class of bilinear timefrequency representations, a unifying framework for
several of the timefrequency distributions that have been proposed in the literature.
Some recent research has been directed toward deriving representations starting
from the definition of Cohen's class to satisfy specific constraints required by a par
ticular application. Examples of this approach are the ChoiWilliams distribution
[Cho89] and the coneshaped kernel distribution proposed by Zhao, Atlas and
Marks [Zha90].
8
2.3 A Summary of Work to Date
Over the past decade there has been a renewed interest in the application of
timefrequency distributions to problems in signal processing. Starting with the
work of Claasen and Mecklenbrauker [Cla80a,b,c] who presented a set of fine papers
on the WD, there has been much research on the theoretical aspects and on the
practical applications of timefrequency distributions. Much research has been di
rected toward establishing distributions that satisfy specific properties deemed desir
able in timefrequency analyses and in their interpretation. The quality of the esti
mates of timefrequency energy distribution and methods to improve these have also
been studied widely. Other theoretical aspects that have received attention include
the inversion of the representation and the interrelationships between the different
distributions [Coh89]. Efficient digital implementations have been developed
[Wil87]. Timefrequency distributions have been applied in various fields where
nonstationary signals arise. The applications have generally been based upon one of
the following i) use of the distribution to present graphically the timevarying charac
teristics of the signal or ii) the use of a particular property of the timefrequency
distribution to describe the evolution with time or frequency of an important signal
parameter. We give here some examples of the wide variety of applications of time
frequency distributions that have appeared in the literature over the past few years.
Imberger and Boashash exploited the instantaneous frequency property of the
WD in the analysis of geophysical and oceanographic signals [Imb86], where the
timeevolution of a significant parameter is represented by the variation of signal
frequency with time. Janse and Kaizer [Jan83] used the WD to graphically represent
the nonstationary signals encountered in loudspeaker design. In a similar manner,
the WD has been used as an aid in the design and analysis of ultrasonic transducers
[Mar86]. In the area of biological signal processing, timefrequency representations
have been applied to the classification of muscle sounds [Bar90] and to eventrelated
9
potentials [Cho87]. These signals are inherently nonstationary and timefrequency
distributions provide consistent patterns rich in detail, making the classification easi
er than with standard timedomain or frequencydomain techniques alone. Starting
with the work of Chester [Che84], several researchers have applied the WD in the
analysis and recognition of speech. The difficulties presented by artifacts in the
analysis of speech have been sought to be overcome by the use of smoothed variants
of the WD [Ril87]. A WD based detector has been developed for the detection of
known signals in noise [Kum84]. A general approach to the detection problem in the
timefrequency plane has been formulated by Flandrin [Fla88]. Breed and Posch
[Bre84] utilized the spatial WD to estimate the range of a target based on the quadrat
ic phase variation due to waveform curvature across the array. The WD has also
been applied in image analysis. The WD, defined for a twodimensional image, is
used to achieve the invariant recognition of objects [Jac84] and applied to the seg
mentation of textured images [Ree90]. In addition to its superiority in joint space/
spatialfrequency resolution, the WD also has the advantage that it implicitly en
codes phase in a single realvalued function. This is a valuable characteristic in
applications to vision where the phase is a critical component.
2.4 Cohen's Class of Bilinear Representations
A generalized class of bilinear timefrequency representations has been estab
lished by Cohen [Coh66] which provides a unifying framework for a number of dis
tributions that have been proposed for timefrequency signal analysis. Cohen's class
is defined by,
Cf(t, ct; 0) = f J f(u + r/2)f*(u r/2)0 (, r)eJ(tu)dudrd (21)
where 0 is an arbitrary kernel that defines the particular distribution, and all inte
grals go from oo to oo. The kernel is chosen to be independent of the signal and of
 10 
time and frequency, which results in the bilinearity and timefrequency shift invaria
nce that characterizes the Cohen's class of distributions. That (21) provides a rea
sonable definition of a general timefrequency representation can be understood by
rewriting it as the Fourier transform of a generalized autocorrelation function K(t,T)
of the signal f(t) [Chd89]
00
Cf(t, a<; 1) = K(t, T)eYdr (22)
00
where
00 00
K(t, T) = J f f(u + r/2)f*(t T/2)P(4, r)eJ()tu)dudQ (23)
00 00
We see that the generalized autocorrelation function is a timeaveraged and lag
weighted value of the timedependent autocorrelation function of the signal. In gen
eral, the two weighting functions represented by the kernel 4((,r) are not indepen
dent. Depending on the choice of the kernel different bilinear timefrequency repre
sentations are obtained. Based on the requirements of timevarying signal analysis
several properties desirable of a timefrequency distribution have been identified
[Cla80c]. These include shift invariance (which holds for all members of the Co
hen's class), marginal properties by which the instantaneous power and spectral den
sity of the signal are obtained, realness and positivity, obtention of instantaneous
frequency and group delay from the local moments of the distribution, preservation
of support in time and frequency, inner product conservation, and compatibility with
linear filtering and modulation. Furthermore, a valuable characteristic of any time
frequency distribution would be for it to be concentrated mainly in those regions of
the timefrequency plane where, according to intuition or experience, signal energy
is expected to be located. Each of the above properties imposes different constraints
on the kernel function, thus limiting the number of suitable representations obtain
able.
 11 
2.5 The Wigner Distribution
The WD belongs to the Cohen's class of bilinear timefrequency distributions. It
is defined for a signal f(t) by,
00
Wf(t, w) = f (t + r/2)f*(t r/2) exp( jor)dr (24)
00
and is obtained from (21) with the kernel /(, r) = 1 The WD satisfies a large
number of desirable properties. Since we refer to one or another of these throughout
this thesis, these properties are summarized here. Proofs can be found in [Cla80a]
and [DeB73].
2.5.1 Desirable Properties
Since the WD is a member of Cohen's class, time and frequencyshifts of the
signal are reflected as corresponding shifts in the WD timefrequency plane, so that
for g(t) = f(t to), we get
Wg(t, ) = Wf(t to, 0) (25)
For g(t) = f(t) exp(jot), we get
Wg(t, ) =Wf(t, o O) (26)
The WD is realvalued, being the Fourier transform of the conjugate symmetric
inner product of the signal. However it is not always positive making it difficult to
interpret as a true signal energy distribution.
The instantaneous signal power and the energy spectral density of the signal are
obtained from the corresponding marginal quantities of the WD as,
00
Sf Wf(t, o)do = If(t)2 (27)
00
 12 
00
fWf(t, )dt = F(o) 12 (28)
00
The WD preserves the time and frequencysupport of the signal, so that
f(t) = 0 for tj > T => Wf(t, w) = 0 for jt]>T. (29)
F(O))=0 for wo > Q=>Wf(t, o)=0 for o>Q0 (210)
However the WD is not necessarily zero at all timeinstants and frequencies that the
signal goes to zero [Coh87].
The firstorder frequency moment of the WD is given by,
00
,(t)12 JWf(t, w)do) = O(t) (211)
GO
where
Qf)(t) = Im = Imdln f(t) (212)
f(t) dt
For complex signals written in the form f(t) = v(t) exp(jo(t)) where v(t), the enve
lope, and O(t), the phase, are real functions, we get from eq. (212),
Qf(t) = '(t) (213)
Hence the average frequency of the WD at any time is the derivative of the phase.
This quantity can be interpreted as the instantaneous frequency of the signal for a
certain class of signals. In an analogous manner, the first moment in time of the WD
of the impulse response of a linear system turns out to be equal to the group delay
[Cla80a].
The WD is invertible in that the signal f(t) can be recovered up to the constant
factor f* (0), as given by
 13 
f(t) = 2 o 1 J W(t/2, o) exp(jt)do (214)
00
This ambiguity is represented in the fact that the signals f(t) and f(t)exp(ja),
where a is a real constant, have the same WD.
An interesting relation exists between the inner product of two signals and that of the
corresponding WDs, known as Moyal's formula and is given by,
I f(t)g*(t)dt 12=2 f f Wf(t, o)Wg(t, o)dtdo (215)
00 CO 00
This property has been used to formulate the optimum detection of signals in the
timefrequency plane. [Fla88].
The WD preserves convolution and modulation. If g(t) = f(r)h(tr)dr then
Wg(t, o) = J Wf(r, )Wh(t ', w)dT (216)
00
For g(t)=f(t)h(t),
Wg(t, 0) = Wf(t, n)Wh(t, () nj)dri (217)
00
In practice the WD is implemented after applying a finitelength, moving window to
the data, centered at the timeinstant of interest, giving rise to the pseudoWD
(pWD). For a signal f(t) and real, symmetric window h(t) the pWD is computed as
pWf(t, ()= f(t+,/2)f*(tr/2)h( )h*()exp(jwor)dr (218)
The pWD is related to the WD by
pWf(t, 0) J Wf(t, q)Wh(0, ( 7)dj (219)
Since
 14
00
Wh(0, ) = h2(r/2) exp(jwr)dr (220)
00
is always a lowpass type of function, the pWD is a frequency smoothed version of the
original WD. This implies that the marginal and local properties of the WD in time
are maintained in the pWD, while the corresponding frequency properties are re
placed by weighted frequency averages.
Further, it has been shown [Jan82] that of members of the Cohen's class that
satisfy the finite support property and the Moyal's formula, the WD has the least
amount of global spread about its center of gravity for any given signal. When sig
nals of the type v(t) exp(j0(t)) are considered, where v(t) is a slowlyvarying, real
function, and 0(t) is a smooth, real function, so that the signal can be considered to
be phasemodulated with a slowlyvarying envelope, the WD is concentrated about
the curve 0'(t) and has significantly lower spread compared to any other timefre
quency distribution. In the case of linear frequencymodulated signals, the WD at
any timeinstant is a delta function at the instantaneous frequency. An enlightening
interpretation of this general behavior of the WD for phasemodulated signals, pres
ented in [Fla84b], is given here. The WD of a signal f(t) computed over a finite data
window to obtain the pWD, can be regarded as an STFT performed on a modified
version of the signal given by the inner product f(t+T/2)f*(tT/2). The bandwidth of
the shorttime spectrum is related to the frequency deviation of the signal within the
data window. Let us assume, without loss of generality, that the data window is
centered at t=0. Then for a signal given by f((T) within the window, the instantaneous
frequency of the signal is represented by Vi(r). However the "instantaneous fre
quency" of the inner product f(T/2)f*(T/2) is given by the quanti
ty [Vj ( ) + Vi( )] This is illustrated by Figure 21.
 15 
Vj Vi
0 0
(a) (b)
Figure 21. Effective bandwidth (shown by shaded area) of windowed
frequencymodulated signal. (a) STFT (b) pWD
[Fla84b]
We see that the total frequency deviation within the window of the inner product is
less than or equal to that of the signal itself, leading to a reduction in bandwidth of
the resulting spectrum. From Figure 21 it is evident that for a given window, the
ideal situation is now only a linear frequency modulation approximation (or more
generally, a skewsymmetric instantaneous frequency about the window center
which gives rise to a constant frequency inner product function). This is a weaker
restriction than the constant frequency law of the signal within the window required
by the STFT. The WD therefore acts as a kind of autoadaptive frequency demodula
tor and provides an order of magnitude improvement in the tracking of timevarying
frequencies. In any case, the instantaneous frequency can always be recovered ex
actly from the first moment of the WD as given by (213).
2.5.2 Interference Terms
Despite the many desirable properties of the WD, enumerated in the previous
section, its application in practice has been limited by the occurrence of interference
 16
terms, or artifacts, in the representation of multicomponent signals. This phenome
non is attributable to the bilinear structure of timefrequency distributions belonging
to the Cohen's class. The WD of the sum of two signals f(t) and g(t) is given by
Wf+g(t, w) = Wf(t, 0w) + Wg(t, w) + 2 Re[Wf.g(t, o))] (221)
where
00
Wf,g(t, C) = ff(t + T/2)g*(t /2) exp( jwr)dT (222)
00
is the crossWD of f(t) and g(t). Thus the WD contains crosscomponents arising
from the interaction of the distinct timefrequency components of the signal. Some
general conclusions about the nature of the interference terms can be drawn from the
following example [Hla84]. Consider a signal h(t) concentrated at (t=0,f=0) and
shifted in time and frequency to get the signals f(t) = a h(t tf) exp(jwft) and
g(t) = b h(t tg) exp(jwgt). The WD of the signal f(t)+g(t) is given by (221) where
Wf(t, 0) = aI2Wh(t tf, ) (Of)
Wg(t, o) = Ib2Wh(t tg, 0 cg) (223)
and the crossterm denoted by If,g(t, 0) is
2 Re[Wf,g(t, w)] = 21abl cos[wUd(t tm) td(0) m) + 0]
X Wh(t tm, CO (m)
with (td, (Od) = (tg tf, COg (i)
Stf +tg Of + cog
(tm, (Om) = (. 9 ) (224)
Thus in addition to the autocomponents given by the corresponding shifted versions
of Wh(t, 0), the resultant WD contains a crosscomponent which is a modulated
 17 
version of Wh(t, w0) located at the midpoint of the interacting components. The mod
ulation is represented as an oscillation in the timefrequency plane in a direction that
is orthogonal to the line connecting the two signal components. The frequency of this
oscillation is proportional to the distance of separation between the two signal com
ponents in the timefrequency plane. Figure 22 illustrates these results.
WW
t
Figure 22. WD of a bicomponent signal [Hla84].
It has been found that the above simple geometrical laws on the nature of the
crossterm hold quite generally even when the interacting signal components are not
of the same form. It should be noted that an interference term always arises in the
WD from the interaction of every two signal components disjoint in timefrequency,
and hence is present even in monocomponent signals where the various timeshifted
components of the signal interact.
An interesting interpretation of the interference phenomenon is provided in
[Jon89] where the WD is viewed as the matched window shorttime Fourier spec
trum of the multicomponent signal. While highly concentrated signal terms are pro
duced in the WD due to the matched windowing of the input signal components,
crossterms arise from the interaction of each signal component with the remaining
window components not corresponding with it.
 18 
The crosscomponents are an integral part of the representation. In them the
relative phase relations of the different frequency components are preserved, and
they enable the marginal properties to hold. However, the presence of interference
terms makes the WD difficult to interpret since these have no real physical signifi
cance in terms of the energy distribution of the signal in the timefrequency plane.
Crossterms have been demonstrated to result in a loss of resolution in the WD of
multicomponent signals [Jon89]. In the case of signals contaminated by additive
noise, the effect of noise on the resulting distribution is severe, limiting the applica
tion of the WD to situations in which the signaltonoise ratios are high.
Much research has been directed toward the suppression of interference terms
in the WD. The most important approach is the use of twodimensional smoothing in
the timefrequency plane, the basis for this being the oscillatory nature of the cross
components [Fla84a]. This operation is given by
00 00
Wf(t, ) = f f(t u, w v)Wf(u, v)dudv (225)
00 00
where 1I(t, o) is the chosen smoothing function that is localized in the timefrequency
plane and is sufficiently regular, but otherwise arbitrary. Since the convolution in
(225) effects a lowpass filtering in the timefrequency plane, the highly oscillatory
crosscomponents are attenuated while the autocomponents are generally smeared.
The optimal smoothing of the WD in terms of eliminating crossterms while retain
ing a highly concentrated distribution, would require adapting the smoothing func
tion in (225) to match the local timefrequency characteristics of the signal. An
example of such an approach applied to unimodular signals has been presented
[And87]. A special case of the timefrequency smoothing represented by (225) is
when the function Il(t,wo) is chosen to be separable in time and frequency, of the form
fl(t,(o) = g(t)xH((o), to get the smoothedpseudo WD (spWD) [Fla84a],
 19 
spWf(t, w) = f h(r) f g(t u)f(u + 2)f*(u )du eJdr (226)
00 00
where H((w) is the spectrum of h(t). Now the frequency resolution is determined by
h(t) while the time resolution is set by g(t), making available a flexible tool with two
degrees of freedom for timefrequency analyses. The timefrequency separable
smoothing function, however, restricts the form of the smoothing region to be paral
lel to the t or faxis.
In general, timefrequency smoothing of the WD leads to the sacrifice of the
marginal and local properties. The exponential kernel of the ChoiWilliams' distri
bution is an example of timefrequency smoothing where the two smoothing func
tions are related in a specific manner that allows the marginal and local properties to
hold, while the crossterms are suppressed by their spreading out in the timefre
quency plane. The resulting distribution is no longer as concentrated as the WD but
provides improved interpretability in the case of multicomponent signals.
In the case of real signals, the WD contains artifacts at d.c. due to the interaction
*of the positive and negative frequency components of the signal. The analytic signal,
derived from the real signal by suppressing the negative frequency spectrum, is usu
ally utilized to overcome this particular problem [Boa87]. Another approach that has
been studied for the elimination of artifacts in the WD of a noisy signal is via the
truncation of an outer product expansion based on the SVD of the WD [Mar85].
Noise components that share the same space as the signal components are not elimi
nated, and some signal components may be truncated leading to signal distortion.
An approach to the elimination of crossterms at the "display" level, by using the
correlation that exists between the WD and the spectrogram of the signal, has been
proposed [Rao89]. The spectrogram of a signal comprised of nonoverlapping fre
quencytime components is free of crossterms, while providing a smoothed repre
sentation of the signal components. A comparison of the WD with the spectrogram
 20 
of the signal could provide an indication of the approximate regions in the timefre
quency plane where signal components exist and hence used to derive a "mask",
which when applied to the WD would yield a representation free of crosscompon
ents but with concentrated autocomponents. Simulation results show that in the
case of narrowband signals in noise and multicomponent signals with individual com
ponents that are not spaced very closely, excellent results can be obtained without
any a priori knowledge of the signal structure. However this method needs further
refinement before it can be successfully applied to signals of more complex nature in
which case it is possible that crossterms may overlap with signal components in the
WD plane, making them difficult to identify.
2.6 The Spectrogram and its Relation to the WD
The spectrogram of a signal is obtained by taking the magnitude squared of the
shorttime Fourier transform given for a signal f(t) and window h(t) by
Ft(w) = J f(r)h(r t) exp( jor)dr (227)
The spectrogram is then
Sf(t, ) = IFt(()2 (228)
The spectrogram is a member of the Cohen's class of bilinear timefrequency
representations. The Cohen's class can also be expressed in the following alternative
form, in terms of the WD, according to,
Cf(t,0;0) = f JO(t r, o )Wf(r, )drd (229)
where
00 OD
S(t,w ) = f (, r) exp(j [trw])drd (230)
0 00
 21
Hence any member of the Cohen's class can be obtained by a linear transformation
of the WD through the kernel c?(t,o) which is related to the original kernel of (21) by
the twodimensional Fourier transform of (230). It is possible therefore to examine
the properties of other members of the Cohen's class based on the known properties
of the WD.
The following expression relates the spectrogram of a signal to its WD [Cla80Oc],
Sf(t, 2)) Jf Wf(r, )Wh(r t, () 4)drdQ (231)
We see that the spectrogram is a timefrequency smoothed version of the WD, with
the smoothing function given by the WD of the window h(t). A direct consequence of
this is the fact that for a given signal waveform, the effective area in the timefre
quency plane of the spectrogram is always slightly greater than that of the WD
[Nut88]. The timefrequency extent of the spectrogram is minimized when the win
dow function is matched to the signal. Further, the timefrequency properties that
hold exactly in the case of the WD, hold only in the sense of averages over the
duration of the window in the spectrogram. That is, the marginals of the spectrogram
in time and frequency are given by
J S(t;, o)do = If(r)121h( t)12dr (232)
00 00
00 00
J Sf(t,)dt j IF(r)121nH(w )2d (233)
00 00
where H(o) is the window spectrum. Similar properties hold for the average frequen
cy and the average time, where the corresponding desired local quantities are aver
aged over the duration of the window with the square of the window, or its spectrum,
as the weighting function. Hence the spectrogram can provide useful information on
 22 
these quantities only for signals for which they do not vary appreciably over the
extent of the. window.
In contrast to the spWD of (226), the timesmoothing and frequencysmooth
ing functions applied to the WD in (231) to get the spectrogram, are not indepen
dent but rather, are related through the predetermined window h(t). This results in
the timefrequency tradeoff that is characteristic of the spectrogram. However, the
amount of timefrequency smoothing effected by Wh(t, cw), the WD of the window
function, is sufficient to eliminate negative values and crossterms between disjoint
components in the spectrogram. The reason for this can be appreciated through the
following example. Assume a signal composed of two components separated in
time. If these components must be nonoverlapping in the spectrogram, we need to
choose a time window that is equal to or shorter than the time interval separating
them. Hence the frequency spread of the WD of the window, being inversely propor
tional to its timewidth, will be equal to or greater than a quantity inversely propor
tional to the time interval between the components. Since the period of oscillation in
the frequency direction of the crossterm in the WD of such a signal is inversely
proportional to the time interval between the two components, it will effectively be
suppressed by the frequencyaveraging wrought by the WD of the window function.
The nature of the spectrogram, by which it is relatively free of interference
terms, has prompted research on determining windows that would maximize concen
tration of signals in the timefrequency plane. Since the WD provides a maximally
concentrated distribution and the spectrogram is related to it via the timefrequency
averaging of (231) it seems obvious that a window function whose WD is most
concentrated in the timefrequency plane would achieve the minimal spread. The
Gaussian window is such a function with its WD being a double Gaussian function in
the timefrequency plane. In order to achieve optimal concentration, the parameters
of the window (namely its time/frequencywidth and orientation in the timefre
23 
quency plane) must be matched to those of the signal components. Such an approach
is embodied in the dataadaptive transform of [Jon87] where the Gaussian window
parameters are adapted to minimize a concentration measure at each location in
timefrequency.
CHAPTER 3
THE ESTIMATION OF FREQUENCY AND FREQUENCYRATE USING THE
WIGNER DISTRIBUTION
3.1 The Concept of Instantaneous Frequency
The need for the analysis of frequencymodulated (f.m.) signals gave rise to the
mathematical description of the concept of instantaneous frequency (referred hence
forth to as IF). An intuitively satisfying and useful description was first provided by
Ville [Vil48]. Ville, based on earlier work by Gabor, related the IF of a real signal to
the derivative of the phase of the corresponding analytic signal. For a real signal s(t),
the analytic signal is given by
a(t) = s(t) + ji(t) (31)
where 9(t) is the Hilbert transform of s(t) :
1 ( s(T)
9(t) = (Hs)(t) = v.p. (t dr (32)
with v.p. representing the Cauchy principle value. The analytic signal has a spectrum
given by,
Sa(f) = 2S(f) f > 0
= S(0) f = 0
=0 f < 0 (33)
If the analytic signal is expressed in the form Sa(t) = a(t) exp(j (t)) then the IF of
s(t) is defined as
 24 
 25 
f(t) 1 d(t) (34)
2an dt
In order for the above definition to be meaningful, it is necessary that the analytic
form of a real signal s(t), given by (31), correspond to the complex form given by
the signal plus its imaginary quadrature component. For a real signal of the form
a(t) cos(o(t)) the following holds [Bed63] :
a(t)cos(o(t)) + jH[a(t)cos(0(t))] = a(t)exp(ja(t)) (35)
if and only if A(f), the spectrum of a(t), lies entirely in the region If < fo and 9{cos
4(t)} exists only outside this region. If the spectra of a(t) and cos(4(t)) are not
separated in frequency the analytic form of their product will be a distorted version of
a(t) exp(j (t)) due to the spectral foldover from overlapping spectra. Only when the
above condition is met, or equivalently when the signal is narrowband, does the de
rivative of the phase in (34) have a meaningful interpretation as the IF of the signal.
3.2 Estimation of Instantaneous Frequency from the WD
The use of a timefrequency distribution to estimate the IF of a signal has been
considered widely, and applied in practice [Boa89]. As stated in Chapter 2 a proper
ty desirable in any timefrequency distribution is for its average frequency at any
time to be equal to the IF of the complex signal. The WD satisfies this property as
given by (213). The discretetime WD (DWD) is defined by
Ws(n, f) = 2 1 s(n+k)s'(nk) exp(j4;rkf) (36)
The average frequency at any timeinstant n is computed as,
1/4
fi(n) = arg[ f exp(j44rf)W,(n, f)df ] (37)
1/4
 26 
For a discretetime complex signal given by,
s(n) = a(n) exp(je(n)) (38)
where a(n) and j(n) are real functions, (37) evaluates to [Cla80b]
fi(n) = [(n 1)] mod 2 (39)
Hence the average frequency of the DWD is the arithmetic mean of the forward and
backward differences of the phase sequence, taken modulo IT (since the DWD is
periodic in frequency with a period equal to ir). Although in the case of a general
complexvalued signal, this reflects an inability to discriminate between frequencies
in the interval (0,ir) and (ir,2ir), it provides an acceptable definition of the IF of an
analytic signal, whose IF is constrained to an interval of length 1r. (39) holds exactly
for the pseudoDWD computed using any (realvalued) window function. Thus the
DWD can be used to estimate the IF since its first moment with respect to frequency
provides an unbiased estimate of the IF for a signal of the general form of (38),
independent of the actual amplitude modulation a(n), or phase modulation <(n).
The presence of noise, however, leads to a serious degradation of the IF estimated
from the first moment of the DWD. The reason for this is made evident by the
following example. For a sequence given by,
s(n) = A exp(j [2=fn + 0]) + z(n) (310)
of a complex sinusoid in noise, where z(n) is complex white Gaussian noise (w.g.n.)
with zeromean and variance r 2, it has been shown [Tre85] that for a high enough
input signaltonoise ratio (SNR) A2/a (310) may be approximated by,
s(n) = A exp(j [2nrfn + + w(n)])
(311)
 27 
where {w(n)}, the phase noise sequence, is real w.g.n. Hence for such an input signal
the estimate of (39) which consists of simply the difference in successive samples of
the noisy phase sequence with no averaging, would exhibit a high statistical variance
even at high input SNR. Since the WD provides a highly concentrated distribution of
signal energy in the timefrequency plane for narrowband signals, a natural alterna
tive is the use of peak detection on the WD at any time to estimate the IF. Taking the
peak of the DWD or maximizing it with respect to frequency is, as seen from (36),
equivalent to the best leastsquares fit of a complex sinusoid to the kernel sequence
{s(n+k)s* (nk)} as a function of k. When the kernel is indeed perfectly sinusoidal, as
is the case when the input signal is of constant amplitude and quadratic phase (linear
f.m.) the DWD peak estimate of IF is unbiased. For signals with nonlinear frequency
modulations and also for any arbitrary amplitude modulation, the DWD peak esti
mate is biased to the extent that the kernel diverges from a perfect sinusoid.
The DWD peak has been applied to the estimation of the IF of phasemodulated
signals in practice [Imb86, Boa89]. Experimental results comparing the perform
ance of the DWD estimate to other conventional methods for the estimation of the IF
and amplitude of f.m. signals in additive, white noise have been reported [Har88]. It
is of interest, therefore, to derive an analytical expression for the performance in
noise of the DWD peak estimator of IF. We begin with an analysis of the DWD peak
in the estimation of IF of linear f.m. signals in additive w.g.n. The theoretical results
are then verified using computer simulations. Finally we discuss, qualitatively, the
behavior of the DWD peak estimate for nonlinear f.m. signals.
3.2.1 Performance of the DWD Peak for Linear F.M. Signals
The input signal is assumed to be given in the complex form by
s(n) = Aexp(j (n)) + z(n)
(312)
 28 
where z(n) is complex w.g.n. with zeromean and variance or z. The DWD of s(n)
is given by (36). In practice the DWD is evaluated over a finite data record of N
samples and is given by
k=L1
Ws(n,f)=2 I s(n+k)s'(nk)exp(j4nkf) (313)
k=L+l
with N=2L+1. The DWD can be looked upon as the frequencyscaled discretetime
Fourier transform (DFT) of the kernel sequence with a scaling factor of 2 and hence
computed using a standard DFT routine. Taking the peak of the DWD is equivalent
to finding the complex sinusoid that fits best, in the leastsquare's sense, the kernel
sequence computed over the finite interval {L+1,L1}, at the timeinstant n. For the
signal of (312) the kernel is given by
s(n + k)s*(n k) = A2 exp(j[O(n + k) (n k)])
+A exp(j[O(n+k)])z'(nk)
+A exp(j[ (nk)])z(n+k)+z(n+k)z*(nk) (314)
The first term is a SxS term, the next two are SxN terms and the last a NxN term,
where S refers to signal, and N to noise. For signals with a quadratic phase function,
given by A exp(j27[ain2 + ao]) the first term in (314) yields the constant frequency
sinusoid A2exp(j8rajnk). The remaining three terms in (314) represent the noise
in the kernel sequence. Since the input noise sequence z(n) is assumed to be zero
mean, white and Gaussian, the distinct noise terms in the kernel are uncorrelated
with each other and with the signal term. The resultant noise power is the sum of the
variances of the three uncorrelated noise components, and is white since the noise
components at different lags are uncorrelated. The noise power in the kernel is thus
given by
 29 
NPkernel = 2A2O+ (315)
It is not Gaussian however due to the introduction of the product term. However, at
high input SNR, the product term is negligible and hence the noise is close to being
Gaussian. Due to the fact that there are only N/2 independent samples in the (conju
gate symmetric) kernel, the variance of the noise in the DWD spectrum is scaled by
the factor N/2 rather than the expected factor N over the input SNR. Hence the SNR
in the DWD spectrum is
A4x(N/2)
SNRDWD = 2A2oz + (316)
The problem is now reduced to that of estimating a constant frequency sinusoid in
white noise with SNR given by (316) using the peak of the DFT magnitude. It is well
known that the peak of the DFT is a maximumlikelihood estimate (MLE) of the
frequency of a sinusoid in w.g.n. At high enough input SNR this MLE meets the
CramerRao lower bound (CRLB) on the variance of the frequency estimate [Rif74].
For an input sequence given by
x(n) =Aexp(j[2rfn + ])+z(n) forn = 0,1, ...N 1 (317)
where z(n) is complex w.g.n. with zeromean and variance a 2, the variance of the
MLE of frequency at high input SNR is given by
6
varDFT(f) = (318)
(2a)2(SNR) (N2 1)
where SNR in (318) stands for the SNR in the DFT of the noisy sequence of (317)
and is given by NA2/ z2. We use this result to derive the variance of the frequency
estimated from the peak of the DWD for a linear f.m. signal in w.g.n.
Due to the frequency scaling by a factor of 2. inherent in the DWD the frequency
estimate obtained from its peak must be corrected by the same factor. This leads to a
 30 
reduction in variance by a factor of 4 for the DWD Substituting the SNR of the
DWD spectrum, as given by (316), into (318) and scaling by 4, we obtain the
variance of the DWD peak estimate of IF for a linear f.m. signal in white Gaussian
noise, which is
6(2A2o2 + o )
varD(f) = (2)22A4N(N2 1) (319)
At high input SNR (A2 > o2) the above reduces to
6o
varDwD(f) = (2x)2A2N(N2 1) (320)
which equals the variance of the MLE of the frequency of a stationary sinusoid in
w.g.n., as given by (318).
We see therefore that the DWD peak is an optimal estimator of IF for linear
f.m. signals in w.g.n. at high input SNRs. As the SNR decreases a threshold effect is
expected to occur as does in the case of the DFT estimate of a stationary sinusoid in
noise. The threshold effect is due the phenomenon (known as outliers) of frequency
estimates falling outside the main lobe in the DWD spectrum and on one of the minor
maxima leading to a sharp increase in the mean squared error (MSE). The probabili
ty of occurrence of an outlier is related to the SNR in the spectrum, and it is known
from the DFT case that the threshold SNR, below which the MSE starts to increase
abruptly, is about 15 dB [Rif74]. The SNR in the DWD is related to the input SNR
through (316). Hence the threshold SNR for the DWD peak estimator is given by
SNRDWD = A4xN/2)= 15dB (321)
2A2o +o 
which at high input SNRs (A2 > oj) reduces to the condition,
A2N
. = 21dB (322)
 31 
Since the DWD is always realvalued, we can detect the peak simply from the
function itself, instead of computing its absolute value. This should lead to a further
reduction in the probability of an outlier since the minor maxima in the DWD spec
trum oscillate in sign. However, computer simulations indicate that this reduction is
slight (between 0.5 to 0.75) resulting in a nearly insignificant effect on the MSE.
3.2.2 Computer Simulations
Computer simulations were used to verify the performance of IF estimation
using the DWD peak. The MSE at various input SNRs for a linear f.m. signal in
w.g.n. was estimated by taking the average of the squared errors between the actual
and estimated frequencies, obtained over between 300 to 3000 trials. The expression
for the MSE is thus,
MSE M (fifi)2 (323)
where M is the number of trials, fi is the actual IF and ft the estimated value. A
64sample data window was used and a large, zeropadded FFT (of length 215)
computed in order to minimize errors in the IF estimate due to frequency quantiza
tion (in practice an efficient interpolation scheme can be substituted). Figure 31
compares the simulation results with the theoretical MSE given by (319). It is seen
that at SNRs above the threshold, the experimentally obtained values closely
match the theoretical values of MSE. The threshold SNR is correctly predicted at 3
dB for the given value of data window size N=64.
3.2.3 Estimation of Instantaneous Frequency of Nonlinear F.M. Signals
In the case of nonlinear f.m. signals, the DWD peak estimator is generally
biased. The magnitude of the bias depends on the extent of the nonlinearity of the IF
curve (or more generally, on the deviation of the IF law from skewsymmetry within
 32 
(MSE)1 dB
0 2 4 6 8 10 12
SNR dB
Figure 31.
Performance of the DWD peak estimator of IF
for a linear f.m. signal in additive w.g.n. (N=64);
101oglo(1/MSE) versus input SNR.
 33 
the data window, as discussed in Section 2.5.1). The bias arises in the attempt to
apply the leastsquare's fit of a constant frequency sinusoid to the kernel sequence
which is no longer narrowband, within the data window. The dependence of the bias
on the nature of the IF curve is illustrated well by the sinusoidal f.m. signal which is
characterized by a timevarying frequencyrate that varies continuously over a range
of values from zero to a maximum. In Figure 32 we plot the IF and the squared
error (SE), between the actual IF and the DWD peak estimate, versus time for a
sinusoidal f.m. signal with a data window of length approximately 7 percent of the
modulation period. It is seen that the bias of the DWD peak estimate, represented by
the squared error, reaches its maximum at the peaks of the IF versus time curve.
These timeinstants represent the maximum deviation of frequency of the kernel, as
a function of lag, from a constant, within the data window. The squared error is zero,
indicating an unbiased estimate of frequency, at the midfrequency, a point about
which the IF curve is skewsymmetric so that the kernel is reduced to a constant
frequency sinusoid within the data window.
For an arbitrary nonlinear f.m. signal, the bias can be controlled by suitably
windowing the input signal before generating the DWD. For a slowly varying f.m.
signal, the window length can be chosen small enough so that the IF law is practically
linear within the window. In this case the performance results of the DWD peak
estimator derived in the previous section can be applied directly. Due to the in
creased main lobe width in the DWD created by the use of a shorter data window, this
approach reduces the bias at the expense of an increase in the variance of the esti
mate in the presence of noise. Alternative window types such as, for example the
Gaussian window which weights the data nearer the window center more, can be
employed in order to reduce the tradeoff between bias and variance of the estima
tor.
 34 
0.5
1.5
1 3 5 7 9 11 13 15
time
Figure 32.
Bias in the DWD peak estimate of IF of a sinusoi
dal f.m. signal as seen by the squared error (SE).
 35 
3.3 Joint Estimation of Frequency and FrequencyRate
The joint estimation of instantaneous frequency and frequencyrate of phase
modulated signals is of importance in problems such as the orbit determination of
satellites, motion compensation in syntheticaperture radar processing, and the high
resolution signature analysis of targets. In all of these situations, the received signal
is the return from a dopplertype radar system, and the frequency and its timederi
vative are important observation parameters. The joint estimation of doppler and
dopplerrate allows the simultaneous measurement of object velocity and the curva
ture of the object path.
Several different approaches to the joint estimation of frequency and frequen
cyrate have been proposed. These include using a linear leastsquares fit estima
tion of these parameters on a finite data set of frequency samples [Kno85], the maxi
mumlikelihood estimation of the instantaneous frequency and frequencyrate of the
parameters of a linear chirp signal [Aba86], and the estimation of frequencyrate of
a sequence of linear f.m. samples after carrying out phase data differencing twice on
the input sequence so that the problem is reduced to that of estimating a constant
phase signal in colored noise [Dju89]. In all the above approaches the form of the
phase function of the input sequence is assumed to be a quadratic polynomial (linear
f.m.) with unknown parameters which are then estimated from the data.
We present here a new approach to the estimation of frequencyrate based on a
methodology suggested by the Wigner distribution. We derive a new kernel from the
signal and show that its spectrum is concentrated at the frequencyrate of the (arbi
trary f.m.) input signal.
3.3.1 Derivation of a New Kernel
The first moment property of the WD by which the average frequency of the WD
at any time equals the IF of the signal, suggests that it is possible to develop a similar
 36
formulation for the frequencyrate based on the difference of phase of WD kernels.
We define a new kernel given by,
rs(t, rj,2) = q(t+T, r2)q'(t, 2 2) (324)
where q(t, r) is the WD kernel of the signal s(t) and is given by,
q(t, ) = s(t+)s (t) (325)
We show here that for the input signal given by s(t) = exp(j (t)) the twodimen
sional Fourier spectrum of r,(t, T1, T2) given by R,(t, 0i, 0)2) satisfies the following
property,
Sf 0122Rs(t, C01,02)dowidCO2 = C x 0"(t) (326)
00 00
where C =j(2n)2. For s(t) = exp(jo(t)), we get from (324),
r,(t, Ti,T2) = exp(j 1(t, Ti,T2)) where,
01(t, rI,r2) = ,(t + + ) (t + )
(t ) + (t 2 ) (327)
Since,
rs(t, Tr, T2) = R,(t, ow1,2) exp(j [w t1 + (2r2])d(old2 (328)
00 0
we have,
62rt, T1, V2) 1 [2
6t162 o 4o= f \ 2Rs(t, ox,w2)dwldw2 (329)
00 00
But from (327) we get,
 37 
62rs(t, r1, r2) 1r1=O,r2=o = jp"(t) (330)
The result (326) then follows from (329) and (330).
The implementation of (326) involves obtaining the twodimensional kernel
function and its DFT prior to taking the joint moment, a computationally expensive
procedure. It is possible to simplify the computation by an approximation at the cost
of a slight loss in the accuracy of the estimate of frequencyrate. We give the simpli
fication here.
If one of the parameters T, or T say T1, is fixed at a value T, then the onedi
mensional DFT of rs(t, T, r2) with respect to T2, denoted by R&,T(t, (02) can be
shown, by retracing the steps (327) to (330), to satisfy,
00
limro j0 2Rs,T(t, (02)do2 =
T T
= limT+o [0'(t +,) '(t )] = 0"(t) (331)
We can rewrite the above expression to obtain an estimate of the frequencyrate as,
T (2Rs,T(t, (02)dw2 p 0"(t) (332)
00
with the approximation becoming increasingly accurate as the interval T is made
smaller. For signals characterized by a purely quadratic frequency modulation,
(331) holds without the limit, and hence the expression (332) holds exactly, irre
spective of the choice of the timeinterval T.
We can simplify computation further, by basing the estimate of frequencyrate
on simply the peak of Rs,T(t, 02) with respect to 0)2, instead of the first moment.
Another advantage of this modification is the anticipated reduction in noise sensitiv
ity of the estimator since, as in the case of the instantaneous frequency estimator
discussed in Section 3.2, it is expected that the discretetime estimator based on the
 38 
frequency moment of Rs,T(t, (02) as in (332) would be more sensitive to noise than
one based on the peak of the same quantity. For f.m. signals for which the frequen
cyrate versus time curve is linear, Rs,T(t, 0)2) is concentrated on this curve. For
example, for the signal given by s(t) = exp(jnat3) (whose frequencyrate equals 6at)
the spectrum Rs,T(t, W2) is a deltafunction at 0"(t) T = 6;ratT and hence the peak
of Rs,T(t, 0)2) at any time t provides an unbiased estimate of the frequencyrate. In
analogy with the behavior of the WD peak in the estimation of IF, the estimate of
frequencyrate from the peak of Rs,T(t, 02) is expected to be biased in the case of
arbitrary f.m. signals with the amount of bias depending on the deviation from skew
symmetry of the frequencyrate versus time curve. Nevertheless, we have now an
orderofmagnitude improvement in the tracking of timevarying frequencyrate
over methods that are based on the assumption of constant frequencyrate within the
data window.
3.3.2 Computer Simulations
In this section we investigate the performance of the frequencyrate estimator
based on the peak of the new transform for a discretetime signal in additive w.g.n.
Since an analytical formulation for the variance of the estimate is far too complex,
we content ourselves with computer simulations to compare the MSE of the estimator
(using an input signal for which it is unbiased) with the corresponding CRLB. The
input signal is given by,
s(n) = A exp(j raon2) + z(n) (333)
which is a complex, linear f.m. signal in additive, complex w.g.n. of zeromean and
variance ao 2, so that the SNR is equal to A2/o The actual frequencyrate is
 39 
equal to ao at all time. The MSE at various input SNRs is determined by averaging
the squared estimation error over a large number of trials.
The estimator is implemented in discretetime by,
k=L1
Rs,T=2(n, f) = > q(n+ 1,k)q'(n 1,k)exp(j4uakf) (334)
k=L+l
where,
q(n,k) = s(n + k)s'(n k) kE (L+1,L1) (335)
It is seen that Rs,T(n, f) can be computed easily as the (frequencyscaled) DFT of the
product of DWD kernels separated in time by T. For the signal of (333) the kernel
of (334) is equal to exp(j4zaok) and hence the peak of (334) with respect to f
provides an unbiased estimate of the frequencyrate.
Figure 33 shows the simulation results using a data window of length 64 sam
ples, and compares the MSE obtained, with the CRLB for the variance of any fre
quencyrate estimator of the signal (333). This lower bound is given by [Dju89]
90oo
var(d) = 2A2N( 1)(N 4) (336)
We see that the MSE of the frequencyrate estimate almost meets the CRLB at high
input SNR. As the SNR decreases a threshold is reached at which the MSE rises
sharply due to the occurrence of outliers in the spectrum R5,T(n, f). The relatively
high value of threshold SNR is attributable to the noisy nature of the kernel of (334)
which is a quadratic function of the signal, and hence contains a number of signal x
noise terms of various powers. Above the threshold SNR, the MSE of the frequency
rate estimate is low. For a given value of SNR, it is lower than that of the IF estimate
from the DWD peak. The reason for this is the narrower main lobe width in the
spectrum Rs,T(n, f) which is the convolution of DWD's as seen from (334).
 40 
100
90
80
(MSE)1 dB
7 8 9 10 11 12 13 14
SNR dB
Figure 33.
Performance of the frequencyrate estimator for a
linear f.m. signal in additive w.g.n. (N=64);
101ogio(1/MSE) versus input SNR.
x
x
: C 3LB
x : si ulatior result
CHAPTER 4
DETECTION AND LOCALIZATION OF NARROWBAND TRANSIENT
SIGNALS
4.1 Detection of Transient Signals
A transient signal can be defined as a waveform of arbitrary shape with a time
duration that is short compared with the observation interval. The detection and
analysis of transient signals is a problem of importance in fields such as biomedical
signal processing, seismic signal processing and acoustics. The detection of a signal
requires processing that is determined by what is known of the signal characteristics,
and when noise is present, what is known about the noise. Depending on the
available knowledge, various methods have been developed for the particular signal
detection problem. A commonly occurring case is that of a signal of known form in
additive Gaussian noise. The appropriate processing is then based on a matched
filter. A classic problem in seismic signal processing, underwater surveillance and
radar is the detection of transient signals with unknown waveforms. When the signal
can be described in terms of unknown parameters it is sometimes appropriate to
consider the signal as a sample function of a random process. When the signal
statistics are known, this knowledge can be used to design a suitable detector
[Hel68]. When the signal is deterministic but of unknown form, a multitude of
different detection methods is available ranging from simple energy detection
[Urk67] to adaptive signal processing techniques that are based on assumed signal
models. The unknown signal parameters are estimated from the observed data in
order to derive the detection statistic [Por86, Liv88]. Another type of unknown signal
 41 
 42 
situation is that of a composite signal of possibly overlapping, unknown, multiple
wavelets which is effectively treated by cepstrum techinques [Kem72].
Here we address the problem of the detection and localization of narrowband
transient signals of unknown waveform and short duration that are embedded in a
noisy background comprised of quasiharmonic and random components. Such
signal conditions are found in passive sonar and in vibration signal monitoring,
where the transients are typically of low energy and only an approximate knowledge
of their frequency location is available. We begin with a review of some recently
proposed methods for the detection of such signals. We then give some practical
examples of the problem we seek to solve and describe desirable detector
characteristics. We discuss and compare systematically, the application of the
Wigner distribution and the spectrogram (which is closely related to some of the
existing approaches to the similar problem). We propose a detection scheme based
on the Wigner distribution and derive analytical results for the performance of the
detector for transient signals in additive Gaussian noise.
4.2 Existing Approaches To The Detection And Localization Of Narrowband
Transient Signals
A recently proposed approach to the detection and localization of transient
signals is based on the Gabor representation, a timefrequency representation which
is closely related to the shorttime Fourier transform [Fri89]. The Gabor
representation is inherently localized and therefore well suited for the modeling of
the transients that are assumed to be concentrated in time and frequency. The Gabor
coefficients { Cm,n } of the given signal, that contains possibly the transient
waveform and additive noise, are computed from the following expression,
00
y(t) = I Cm,n g(t na) exp(j2mafit) (41)
m,n=oo
 43 
where y(t) is the received signal, g(t) is a chosen window function and a and 13 are
detector parameters, with op = 1. The detector, based on a generalized likelihood
ratio test, examines sets of coefficients { Cm,n } spanning relatively short time
intervals in order to detect the transient waveform and localize it in time. The
performance of the detector is dependent on how closely the window function, g(t),
matches the transient waveform. Thus the Gabor detector is not robust to parameter
variations in the transient waveform. The authors remark, in particular, the need to
increase the robustness of the detector to timescale variations of the signal. This
problem of the variability in detector performance with changes in timescale of the
transient signal has been addressed through the use of timescale analysis [Gar90].
The wavelet transform is proposed for the detection of transient signals with
unknown arrival times and waveforms that are known to within a scale factor. When
the analyzing wavelet is chosen to be of the form of the transient signal to be
detected, the wavelet transform provides good localization of the transient in
timetimescale space Since the parameters that are represented by the wavelet
transform are the time of arrival and scale, the frequency content of the transient
cannot be arbitrary but rather must be a strict function of the timescale. The author
restricts himself to the class of transient signals that are monochromatic and have a
frequency that is inversely proportional to the time duration, an example of which is
the class of vibration transients with a damping factor proportional to the frequency.
The shorttime power spectrum is widely applied in passive sonar signal
processing to detect transient signals of short duration by the temporary increase in
energy in portions of the received signal spectrum, caused by the arrival of a
transient [Kni81, Wol83]. This procedure is equivalent to performing narrowband
filtering within energy detection. The detector parameters are the filter bandwidths
and averaging times which determine its performance for given signal conditions.
 44 
All of the approaches to transient signal detection and localization discussed
above are based on the analysis of the received signal using fixed window functions.
Since the signal waveform is not known a priori the window. function cannot be
optimally chosen, and hence these methods suffer from performance constraints
imposed by the particular window and the consequent lack of robustness to
parameter variations of the transient waveform. Further the time and frequency
resolutions obtainable in the representation of the received signal are not
independent but are characterized by a fixed tradeoff determined by the window
function.
4.3 Desirable Detector Characteristics
To motivate our approach to transient detection, we give some practical
examples of the type of problem we have attempted here. One such is the detection
and classification of mechanical parameters of ships based on underwater acoustic
sounds [Lou87]. Ship propeller noise can be modeled as a burst of pressure pulses
caused by collapsing cavitation bubbles once every rotation for each blade. The
received signal is further contaminated by noise from other sources like gearboxes,
turbines and electric motors. Detection of the cavitation and the power in the
cavitation frequency is of interest and needs to be achieved without a precise
knowledge of the cavitation frequency or bandwidth. A similar problem arises in the
detection of transient vibration signals in the field of machine monitoring. The
transient vibration signal is typically of low energy and is buried in a background
comprised of a number of harmonic components of relatively high energy,
contributed by other sources interacting with the system in which the transient
vibration originates. Additionally the vibration signal may contain random noise.
We concern ourselves here with the class of transient signals that are
monochromatic and of short, but otherwise arbitrary, time durations. The transient
 45 
signal is possibly in background noise comprised of deterministic and random
components. Since the transient signals of interest are localized in time and in
frequency, it is natural to consider the application of a timefrequency representation
to the detection problem. Our goal is to develop a detector which:
(i) requires no a priori information on the actual time waveform or the frequency of
the transient signal;
(ii) provides high enough frequency resolution to discriminate the transient signal
from the background signal components, so that the detector output is not influenced
by energy changes in the background, that are nonoverlapping in frequency with the
transient signal. Further, we would like the detector to provide high time resolution,
so that the transient signal is well localized in time;
(iii) is robust to timescale changes of the transient waveform, so that the detector
maintains a reasonable level of performance over a wide range of transient time
durations.
Since we do not intend to make any assumption on the time waveform of the
transient signal but intend to detect the transient by recognizing its spectral signature,
or the temporary increase in energy in a region of the spectrum, we are restricted to
the choice of a nonparametric method for processing the received signal. The
shorttime power spectrum and the Wigner distribution are among the possible
timefrequency representations that suggest themselves for the detection problem.
As discussed in Section 2.6, for a given waveform, the WD provides a more
concentrated representation in the timefrequency plane than does the spectrogram.
The extent of timefrequency spread of the spectrogram is minimized when it is
computed using a window matched to the signal waveform, which is not possible
when the transient is of unknown form.
The application of the WD to the detection of signals in noise has been
considered before [Kum84, Fla88]. Flandrin [Fla88] has shown how the optimum
 46 
detection of a signal of known structure (known deterministic signal, or random
signal of known statistics) in Gaussian noise can be achieved using the WD. The
inner product invariance property of the WD permits the implementation of the
optimum detector in the timefrequency plane, which for nonstationary input signals
can lead to a considerable simplification in detector structure. In this chapter, we
address the application of timefrequency methods to the detection of signals of
unknown form in additive noise.
In the following subsections, we discuss the application of the WD and the
spectrogram to the detection of transient signals in noise, where the received signal
satisfies the assumptions outlined in this section, and compare their performance.
4.3.1 High TimeFrequency Resolution
Since the received signal contains the transient waveform embedded in a noisy
background comprised of spectral components originating in other interacting
phenomena and which could be varying slowly in amplitude, we need high frequency
resolution in our signal representation, in order to isolate the energy contribution in
the narrow spectral band where the transient is expected to be concentrated. For
relatively closely spaced components, this would require us to choose a greater
transform length or data window size. In the case of the spectrogram, an increase in
data window length has a direct bearing on the time resolution of the transient signal.
Since a longer data window would lead to the spreading of the transient signal energy
over a greater time duration, it is expected that the time localization of the transient
waveform would degrade as the window length is increased. Further, increasing the
data window length also contributes to a decrease in the SNR and hence reduces the
detectability of the signal in noise. The pseudoWD, however, allows the increase of
data window size to achieve high enough frequency resolution without a sacrifice in
. time localization. The difference in combined timefrequency resolution between
 47 
the spectrogram and the pWD is illustrated by Figures 41 to 43. Figure 41 shows
a time signal, comprised of the sum of a sinusoid of constant amplitude at frequency
0.09 and a transient waveform with a Gaussian envelope at frequency 0.12. In this
example and all other examples in this chapter, we use the complex form of the
signal when computing both the DWD and shorttime power spectrum. The
simulated complex signals have spectra confined to the frequency region [0,'rr] to
avoid aliasing in the DWD [Cla80]. (The DWD shown here, is computed using a
standard FFT routine and hence is a frequencyscaled (by a factor of two) version of
the actual signal DWD.) Spectrograms obtained with two different data window
lengths and spanning the timeinstant of arrival of the transient are compared in
Figure 42. (In this and all other figures in this chapter, N refers to the data window
length in samples and M to the length in samples of the timeaveraging window used
to compute the spDWD. We use rectangular data and timeaveraging windows
throughout.) It is seen that the longer data window provides the superior frequency
resolution of the distinct signal components, but the energy of the transient is spread
in time, so that the actual time of arrival of the transient signal is not obvious. In
contrast, Figure 43 which shows the spDWD of the signal using the same two data
windows, does not display any perceptible tradeoff in the time localization of the
transient with the improved frequency resolution of the longer data window.
One may be concerned about the effect of crossterms in the WD on the
detection of the transient waveform when the received signal contains multiple,
closely spaced frequency components. The crossterms contributed by these
components degrade the resolution of the WD and may lie in the spectral region of
the transient signal and hence influence the detector output. In such cases it is
necessary to employ some timesmoothing of the WD in order to improve the
frequency resolution by attenuating the crossterms. This operation leads to a
corresponding loss of time resolution. However, since the timesmoothing window is
 48 
. VALUES
I I I I I I I I I
BEG .05120 SEC
MID .1280 SEC
ENO .2048 SEC
Figure 41. Gaussian transient in sinusoid.
1024k
wDn i
0
1024
cfU I40 .
 49 
time
time
frequency
(a)
frequency
(b)
Figure 42. Effect of data window length on spectrogram of Gaussian
transient in sinusoid. (a) N=128 (b) N=512.
 50 
time
frequency
(a)
time
frequency
(b)
Figure 43. Effect of data window length on spDWD of Gaussian tran
sient in sinusoid. (a) N=128, M=40 (b) N=512, M=40.
 51 
chosen independently of the data window, there is no immediate tradeoff between
time and frequency resolutions as in the shorttime power spectrum. For instance,
we can set the data window duration in time, long enough to achieve the frequency
resolution we need, and then depending on the extent of crossterm interference,
choose the appropriate length for the timeaveraging window. In most cases a
timesmoothing window that is much shorter in duration than the data window will
be adequate to attenuate the crossterms sufficiently to prevent them from
influencing the output of the detector. Hence using the spWD allows the independent
control of time and frequency resolutions. Figures 44 to 46 illustrate this point
well. Figure 44 shows the time waveform of a Gaussian transient of frequency 0.12,
embedded in a signal comprised of two harmonics closely spaced in frequency with
the transient signal. Figures 45 and 46 show respectively the spectrogram and the
spDWD of the waveform in Figure 44, spanning the time of arrival of the transient.
The spectrogram was computed using a data window length of 84 samples, to match
the duration of the transient waveform, so as to achieve good time localization. We
see that the distinct frequency components are just about well resolved in the
spectrogram. Figure 46 shows the 512sample spDWD of the signal with a
timesmoothing window length of 40 samples, in order to attenuate the crossterm
generated by the sinusoids. It is seen that we achieve better time and frequency
resolutions in the spDWD over the spectrogram.
4.3.2 Robustness to Transient Duration
In order to maintain a good detection performance even with relatively short
signal durations, it is necessary for the detector not to be strongly sensitive to changes
in the time duration of the received transient waveform. For this condition to hold,
we must base our detection statistic on a parameter or feature of the signal that is
invariant to such changes, but represents only the transient signal amplitude. A
 52 
. VALUES
a
1024
BEG .02560 SEC
1i ll
MID .1024 SEC
END .1792
Figure 44. Gaussian transient in sum of two closely spaced sinusoids.
VOr W03i.
edU4 I   I  I  I _____ I I I I_________  __
,,
J
I I I i
 53 
time
__AFigure 45.
Figure 45.
frequency
Spectrogram of Gaussian transient in sum of two closely
spaced sinusoids (N=84).
 54 
time
frequency
Figure 46. spDWD of Gaussian transient in sum of two closely spaced
sinusoids (N=512, M=40).
 55 
change in the time duration of the transient signal corresponds to a change in the total
signal energy, and therefore will affect the magnitude of the peak of the transient
signal spectrum. In other words, a decrease in pulsewidth will lead to a
corresponding decrease in magnitude of the peaks of both, the spectrogram and the
WD, at the frequency and time of arrival of the transient signal and hence a
degradation in performance of any detector based on these quantities. The
instantaneous signal power, on the other hand, is a quantity that reflects the
amplitude of the transient signal and does not depend on its time duration. The WD
of the received signal can be used to obtain an estimate of the instantaneous signal
power.
The WD itself is simply the Fourier spectrum of the signal kernel and hence the
magnitude of the peak of the WD is a function of the total signal energy. However,
from the timemarginal property of the WD we have for a signal u(t),
I Wu(t, f)df = u(t)2 (42)
00
We can thus recover the instantaneous signal power at any time by integrating the
WD over all frequency, at that time. Equation (42) holds also for the pWD
irrespective of the time window used to compute it. The spectrogram, on the other
hand, gives the average signal energy over the duration of the window h(t) with the
square of the window as a weighting function,
Su(t, f)df= J u(r)2h(rt)j2dr (43)
00 00
Figure 48 compares the robustness of the timemarginals of the spectrogram and
the WD, computed using fixed data windows, to variations in the time duration of the
received transient signal. Figure 47 shows a time waveform containing three
 56 
BEG .05120 SEC MID .1280 SEC END .2048 SEC
Figure 47. Sum of Gaussian transients of same peak amplitude, but dif
fering pulse durations.
 57 
60 
50
40
30
20
10
0
600
Figure 48.
1000 1200 1400 1600 1800
time
Timemarginals of the spectrogram and the WD a compar
ison;  spectrogram; WD.
800
 58 
transient signals of the same center frequency and peak amplitude but different time
durations. Figure 48 compares the quantities (42) and (43) as a function of time
for this waveform. The magnitude of the integral of the spectrogram is a strong
function of the time duration of the transient, since it reflects the total energy of the
signal over the duration of the window. The integral of the WD on the other hand is
proportional to the instantaneous signal power and hence reaches the same peak
value for all three transient signals.
Since the transient signals of interest are assumed to be concentrated in
frequency, the above property of the WD by which the instantaneous signal power
can be recovered, suggests a detector based on integrating the WD of the received
signal over a narrow bandwidth corresponding to the region of significant spectral
support of the transient signal. If the WD of the signal u(t) is contained entirely in the
region (fAf,f+Af) at any time and is zero outside, then the instantaneous signal
power may be recovered by integrating the WD only over this frequency band thus,
f+Af
f Wu(t, f)df = u(t) (44)
fAf
Integrating the WD in frequency provides the further advantage of smoothing
out crosscomponents between multiple transient signals that are closely spaced in
time. Signal components that are separated in time give rise to crossterms in the
WD midway between the two interacting components. These crossterms are
oscillatory in the frequency direction and, therefore, are suppressed by averaging the
WD in frequency, an effect achieved by (44). Figure 49 shows the time waveform
of two almost overlapping transient signals. Figure 410 shows the 200sample
spectrogram and the 512sample DWD, for this signal. Crossterms are clearly seen
in the DWD and in this case, due to the overlap between signal components, in the
spectrogram as well. Figure 411 compares the quantities of (42) and (43) for this
signal, where the DWD is computed using a 512sample data window and the
 59 
BEG .05120 SEC MID .1152 SEC END .1792 SEC
Sum of two closely spaced in time Gaussian transients.
Figure 49.
 60 
time
frequency
(a)
time
frequency
(b)
Figure 410. Crossterms in the representation of overlapping Gaussian
transients. (a) Spectrogram (N=200) (b) pDWD (N=512).
 61 
60
50
40
3
3 30\
20
10 " ,
900 950 1000 1050 1100 1150 1200 1250
time
Figure 411. Timemarginals of the spectrogram and the WD of overlap
ping transients a comparison;  spectrogram; WD.
 62 
spectrogram using a 128sample window. The transient signals are well resolved in
time in the DWD, but clearly are not as well resolved in the spectrogram inspite of
having matched the spectrogram window duration to the duration of each of the
transient signals.
4.4 A Detector Based on Instantaneous Power
Based on the considerations of the earlier sections we propose here a detector
for transient signals in additive noise based on monitoring the received signal power
in localized regions of the spectrum. The received signal is assumed to meet the
assumptions stated in Section 4.3. Figure 412 shows the configuration of the
detector. The received signal is filtered so as to bandlimit the noise to the
u(t)+n(t) Filter Compute Wy(t, v) maxf ,
SH(f) WD A
Decisior
Window threshold 1p.
v(t) Detector
Figure 412. Narrowband transient signal detector
approximately known broad spectral region that the transient signal is expected to lie
within, while the signal of interest is passed unattenuated. This operation is
necessitated by the fact that the WD estimate for a signal in additive, white noise is
characterized by infinite variance. (Detailed analytical results on the noise
properties of the WD are presented in the next section.) The input signal is next
windowed by the finitelength window v(t) and then the WD is computed,
incorporating timesmoothing if necessary. The WD is integrated over a narrow
 63 
frequency band surrounding the center frequency of the transient waveform to obtain
an estimate of the instantaneous signal power as,
f+Af
r(t, f)= J Wy(t, v)dv (45)
fAf
This quantity is then compared with a threshold to determine whether a transient
signal is present. When the exact frequency of the transient signal is not known, we
obtain a maximum likelihood estimate of the frequency by computing (45) for
several (overlapping) frequency bands throughout the approximately known spectral
region of the transient signal, and choosing the maximum output to detect the
presence of the transient signal. The quantity Af is a parameter of the detector and is
fixed at a predetermined value. It represents the largest bandwidth or equivalently,
the narrowest duration transient pulse, that can be effectively detected.
We illustrate the high combined timefrequency resolution characteristic of the
WDbased detector via the example of the timeseries of Figure 413. The signal is
comprised of two shortduration Gaussian pulses in an additive background created
by a closelyspaced (in frequency) narrowband component with slowlyvarying
amplitude. Figure 414 shows the detection statistic of (45) with the region of
integration being the narrow frequency band surrounding the center frequency of the
transient pulse. It is seen that the detection statistic is influenced mostly by the
arrival of the transient signals. Figure 415 shows the similar detection statistic
computed using the spectrogram with two different data windows. The spectrogram
succeeds in resolving the transient signal energy from the background spectral
energy changes only at the expense of a loss in the time localization of the transient.
In the following sections, we analyze the performance of the WDbased
transient signal detector in the presence of additive random noise. We begin by
studying the SNR in the WD estimated from a noisy signal and then derive the
performance of the detector.
 64 
1500
1000 
500
0
500
1000
1500
600 800 1000 1200 1400 1600 1800
time
Figure 413. Transient pulses in a narrowband background with slowly
varying amplitude.
 65 
8
7
6
S 5
0 4
3
2
1
600 800 1000 1200 1400 1600 1800
time
Figure 414. The WD detectionstatistic (N=512, M=128) at the frequency
of the transient signal.
 66 
7
6
5
4
3
2 
I II I
0 
600 800 1000 1200 1400 1600 1800
time
Figure 415. The detection statistic based on the spectrogram with differ
ent data window lengths : N=128;  N=256, at the
frequency of the transient signal.
 67 
4.5 Noise Performance of the WD
In the previous sections we have considered the effect of deterministic,
background signal components on the detection of transient signals using the WD.
We now examine the properties of the WD estimated from a signal contaminated
with additive random noise, in order to investigate the noise performance of the
detector. Nuttall [Nut88] has derived the mean and the variance of the WD estimated
from a such a noisy waveform, and we begin by summarizing his results. We then
discuss the bearing of these results on the choice of detector parameters, and
compare the performance in noise of the WD estimate with the performance of the
magnitude of the shorttime Fourier transform under similar signal conditions.
(Unless otherwise indicated, all integrals go from oo to oo.)
The received signal is given by,
x(t) = u(t) + n(t) (46)
where u(t) is the transient signal waveform, and n(t) is zeromean, stationary,
complex Gaussian noise that satisfies,
n(t)= 0 and n(tl)n(t2) = 0. (47)
The noise power spectral density is given by
Gn(f) = n(t + ) n*(t ) exp( j2arfr) dr (48)
To compute the WD, we must window the received signal. The windowed waveform
is given by,
y(t) = v(t) [u(t) + n(t)] (49)
where v(t) is the chosen, finiteduration window function. The WD is then calculated
as,
 68 
Wy(t,f) = y(t +2) y'(t 2) exp(j2'fr) dr
= a+b+c+d
a = Rwv(t, r)
b= fRvv(tTr)
c =f Rwvv(t, r)
d = fRw(t, T)
d= JRw(t,7)
u(t + )
n(t + )
u(t+ )
n(t +)
2
u*(t 2)
exp( j2;rfr) dr
exp( j2rfr) dr
n*(t ) exp(j2nfr) dr
2j
u 2
exp(j2nfr) dr
Rw(t, r)=v(t+) v'(t )
2 2
(415)
The quantities a and b, are respectively SxS and NxN terms, while c and d are SxN
terms, where S refers to signal and N to noise. Based on the assumptions stated
earlier in this section, Nuttall [Nut88] has derived the following expressions for the
mean and the variance of the WD estimate of the noisy signal. The mean is given by,
Wy(t, f) = Wv(t, f) *f [Wu(t, f) + Gn(f)]
(416)
where W,(t, f) is the WD of the window v(t) and Wu(t, f) that of u(t), and *f refers
to convolution in the frequency variable. The variance of the WD estimate is given
by
var[Wy(t, f)] = 2W,(t, f) Ga(f) + 2 Gn(v) IB(t, f 12 dv
(NxN) (SxN)
(417)
where
(410)
(411)
where
(412)
(413)
(414)
 69 
where
Gnn) Gn(f) *G(f) (418)
B(t, f) Wv(t, f ) U(v) exp(j2n t) dv (419)
and U(f) is the Fourier transform of u(t). The variance of the WD estimate is
comprised of two terms, one each contributed by the NxN portion and the SxN
portion of the total noise. It can be seen from the NxN term in (417) that the
variance of the WD estimate depends on both the window function v(t) and on the
noise power spectral density Gn(f). If we set v(t) = 1 for all t, or do not use any
timewindowing of the data, then Wv(t, f) = 6(f), and the NxN term in (417)
becomes infinite. Furthermore, irrespective of the time window, if the input noise is
white, Gn(f) = No for all f, making Gnn(f) and hence the resultant variance infinite.
Hence, in order for the variance to be finite it is necessary for the time window v(t) to
be of finite length and further, also necessary for the noise to be bandlimited.
Figure 417 compares the 512sample DWD with the 256sample spectrogram
of the noisy signal of Figure 416, consisting of a Gaussian transient in w.g.n. of
power spectral density No and total transient signal energy Eo so that Eo/No =
7.6. The DWD is very noisy compared to the spectrogram. This is due in part to
the longer data window used in computing the DWD but, more significantly, due to
the noisy nature of the DWD attributable to the existence of NxN as well as SxN
terms, and the fact that the noise is white.
We next study, quantitatively, the SNR of the WD estimated from a noisy
waveform for an example case. Our aim is to compare the SNR of the WD with that
of the magnitude of the shorttime Fourier transform for transients with time
durations that are short with respect to the observation window. For the sake of
 70 
4096 VALUES WO31.
2048
4096
4095 VALUES W032.
2048
4096
BEG 0.0 SEC MID .1280 SEC END .2560 SEC
Figure 416. Gaussian transient in white, Gaussian noise (Eo/NO = 7.6).
 71 
time
time
frequency
(a)
!M ijAy JA^^A/AM< WA^^^^AA^~^lf^ ^f^ ~ ^
frequency
(b)
Figure 417. Effect of additive noise on (a) Spectrogram (N=256) (b)
pDWD (N=512).
 72 
analytical tractability we choose the signal, window and filter to be Gaussian in form.
The received signal is given by (46) where,
u(t) = a exp( t2/2oI) (420)
and n(t) is white Gaussian noise of power spectral density No. The received signal is
windowed in time by a window v(t) of effective length p, and filtered by a transfer
function H(f) of effective bandwidth B so that the WD estimate has finite variance.
The window length p and the filter bandwidth B are assumed to be large enough so
that the transient signal u(t) is passed essentially unaltered, while the noise is
attenuated. The configuration is as shown in Figure 412, where
v(t)= exp( t2/2p2) (421)
H(f) = exp( f2/2B2) (422)
The filtered noise spectrum is therefore,
Gn(f) = NoIH(f)12 = No exp( f2/B2) (423)
The quantities required to calculate the mean and the variance in (416) and (417)
for this example are given below.
U(f) = aoo 2 exp( o2xr2f2) (424)
W(t, f)= 2a20ro I exp( t2/o2 oO42f2) (425)
W(t, f) = 2p r exp( t2/p2 p242f2) (426)
Gmn(f) = N2 B exp( 2f2/B2) (427)
 73 
We are interested in the SNR of the WD estimate given by
= difference of mean outputs
standard deviation of output
where the numerator refers to the difference in mean output with signal present from
that with signal absent. We compute the SNR of the WD estimate at the time of
arrival t=0, and the center frequency f=0, of the transient signal,
Wy(t, f) *' Wu(t, f)
SNRWD = W(tf) W(tf) at t = 0, f = 0 (429)
[var(Wy(t, f))] II
with the variance computed as in (417). The expression (429) was evaluated
numerically over a range of values of pulse width oo, keeping the parameters a, p, B
and No fixed. The values chosen are p = 40.0, B = 4.0, No = 0.2, and a = 2.0. Figure
418 shows a plot of the SNR of the WD estimate versus pulse width. In order to get a
better understanding of what the numbers represent, we plot on the same graph the
SNR of the, more familiar, magnitude of the STFT under identical conditions, i.e.
signal, noise, window and filter. The analytical expression for the SNR of the
magnitude of the STFT estimated at the time of arrival and center frequency of the
transient is derived below.
The STFT is computed at t=0, by applying the window v(t) to the received signal
x(t) to get,
Sy(0, f) = x(r) v(r) exp(j2orfr) dr (430)
The mean signal output at t=0, f=0 is then,
Sy(0, 0) = U(f) V(f) df (431)
The variance of the output noise is,
var[Sy(0, 0)] = Gn(f) IV(f)2 df (432)
 74 
SNR
3
2
1 7
0
0
Figure 418.
1 2 3 4
00
Output SNR of the WD and STFTI estimates of a Gaus
sian pulse in w.g.n. versus pulse duration.
 75 
where U(f) and Gn(f) are given by (424) and (423) respectively, and V(f) is the
spectrum of v(t),
V(f) = p exp( p 222f2) (433)
Substituting the above expressions in (428) for the output SNR and simplifying, we
get the SNR in the magnitude of the STFT estimate as,
SNRsTFT = a oo x [I 1 241/4 (434)
[N0(ol + p2)]1/
From Figure 418 we see that for the parameters of this example, the SNR in the WD
estimate is comparable to, and in fact slightly better than, that of the magnitude of
the STFT. The WD estimate however, degrades rapidly with an increase in filter
bandwidth B. The STFT magnitude is not as sensitive to changes in B. In summary,
with the proper choice of filter and window parameters it is possible to achieve in the
estimation of the WD a noise performance on the same order as that of the
magnitude of the STFT, while obtaining the superior timefrequency resolution of
the WD.
4.6 Noise Performance of the WD Transient Detector
A transient signal detector based on estimating the instantaneous signal power
using the WD, was presented in Section 4.4. We use the results presented in the
previous section on the noise performance of the WD to analyze the performance of
the detector of Figure 412 for a transient signal in additive noise. The detection is
based on the integral of the WD estimated from the noisy signal, computed over a
fixed frequency region given by (fAf,f+Af). At time t and frequency f, this quantity
may be written in the following equivalent form,
f+Af
r(t, f) = JWy(t, v) dv = Wy(t, f) *f D(f) (435)
fAf
 76 
where,
D(f) = 1 for f E (Af,+Af)
=0 else (436)
But y(t) = v(t)x(t) and hence we have,
Wy(t, f) = Wv(t, f) *f Wx(t, f) (437)
Substituting in (435) we get,
r(t, f) = Wv(t, f) *f Wx(t, f) *f D(f) (438)
or,
r(t, f) = Weq(t, f) *f Wx(t, f) (439)
where,
f+At
Weq(t, f) = Wv(t, f) *f D(f) = JWy(t, v) dv (440)
fAf
From (439) we see that the effect of integrating the WD of the received signal over
the frequency interval (fAf,f+Af) is equivalent to a modification of the window v(t)
so that its WD is now given by (440). For the window given by (421), this
expression evaluates to
1 t2
Weq(t, f) = 1exp() [erf(2np(f + Af)) erf(2arp(f Af))] (441)
where erf refers to the error function and is given by
erf(y) = 2 exp( t2)dt (442)
0
 77 
The mean detector output when the signal is present is then obtained from (416) by
replacing W,(t, f) with Weq(t, f) to get
El [r(t, f)] = Weq(t, f) *f [Wu(t, f) + Gn(f)] (443)
The mean detector output when the signal is absent is
Eo[r(t, f)] = Weq(t, f) *f Gn(f) (444)
The variance of the detection statistic when the signal is present is obtained from
(417) and is given by,
var [r(t, f)] = 2Wiq(t, f) *f G(f) + 2 f Gn(v) IB(t, f )2 dv
J 2 (445)
(NxN) (SxN)
where
B(t, f) = Weq(t, f ) U(v) exp(j2nrvt) dv (446)
The variance of the detection statistic with only noise present is simply the first term
(NxN) of (445). Since an exact derivation of the probability distribution of the
detection statistic is not mathematically tractable due to the presence of nonlinear
noise terms, we do not attempt to determine the probability of detection. We
characterize the performance of the detector instead by a SNR measure given by,
SNRr El[r(t, f)] Eo[r(t, f)]
I2 (var [r(t, f)] + varo[r(t, f)])
Using the window and filter of (421) and (422), the above expression was
evaluated at the time of arrival, t=0, and center frequency, f=0, of the transient signal
u(t) of (420), using error function tables and numerical integration over a range of
pulse widths. The constant parameters were: p=40.0, B=4.0, a=2.0, and No = 0.2.
The results are plotted in Figure 419 which shows the output SNR as a function of
 78 
the transient pulse width co relative to the window length p, for different values of
Af. The energy contained in the transient signal is proportional to the pulse width oo
and hence Figure 419 also represents the variation of output SNR with input SNR
(given by signal energy /noise power spectral density), with peak signal amplitude
kept constant. The SNR obtained in the WD peak, which is the limiting case with Af
> 0, is also plotted for comparison. From Figure 419 we see that the output SNR
for Af > 0, rises rapidly with transient pulse duration, and reaches a level that is then
maintained as the transient duration increases further. The insensitivity of the output
SNR to the energy of the received transient signal beyond a threshold pulse duration,
is due to the detection statistic being based on an estimate of instantaneous signal
power, or equivalently on the signal amplitude, rather than on the total received
energy. For a given Af, we obtain from (435), an estimate of the signal power in the
bandwidth (Af,+Af) at the instant t. As the pulse width oo is increased starting
from a small value, the bandwidth of the WD spectrum of the transient pulse
decreases correspondingly. This leads to an increasingly greater fraction 6f the
transient signal's WD spectrum falling within the region of integration (Af,+Af),
and hence the rapid rise in output SNR. Above a threshold pulse width 00, the entire
transient signal WD spectrum is concentrated within the frequency region (Af,+Af),
so that now any further increase in GO does not add to the mean signal output in
(443), but leaves it unchanged. The variance of the output is also a function of the
pulse width, when all other parameters are fixed as is the case we are considering
here. The power in the SxN term of (445) is dependent on both the received signal
power and the noise power. Since only a signal parameter, i.e. the pulse width, is
being varied here, the received noise power is expected to remain constant. Hence
the power in the SxN term increases with the pulse duration until the threshold pulse
width and then maintains a constant level for any further increases of pulse width.
The power in the NxN term of (445) is constant since it does not depend on the pulse
 79 
SNR 4
3 4
2
1
0
0
Figure 419.
SNR of the WDbased detection statistic versus pulse
width, at various Af.
 80 
width. The above discussion on the dependence of the mean and variance of the
detector output on the received transient pulse width, with all other parameters kept
constant, explains the behavior of the plot of output SNR versus pulse width. When
the detector parameter Af is varied we get a different curve in Figure 419. As Af is
increased we are integrating over an increasing region of the noisy WD spectrum,
which leads to the increase in the total noise power and hence the decrease in the
output SNR of the estimate for a given received signal pulse width. On the other
hand, for higher Af the output SNR attains a constant level more rapidly and at a
lower threshold pulse width since now the wider region of integration includes the
entire WD spectrum of the transient waveform at t=0 for even narrower pulse widths.
The WD peak is a limiting case (obtained as Af>0). As ao is increased, the WD
spectrum at t=0 of the transient pulse becomes increasingly concentrated and hence
the mean detector output, which is simply the value of the WD at f=0 increases. The
power of the NxN term remains constant while that of the SxN term increases with
oo. The net effect is a rapid rise of output SNR with increasing pulse width.
4.7 Analysis of Experimental Milling Machine Data
Detection of tool chatter during highspeed milling operations is a problem that
has benefitted from spectral analysis methods [Smi89]. In stable milling conditions
the vibration signal from the cutting tool is characterized by a spectrum containing
primarily the tooth frequency and its harmonics. As the cut becomes unstable due to
changes in operating conditions, such as the depth of cut or spindle speed, a new
component, the chatter frequency, begins to dominate the vibration signal spectrum.
The exact frequency of the chatter induced vibration is generally not known a priori.
In [Smi90] the magnitude of the shorttime Fourier spectrum of the vibration signal
is used to detect the occurrence of chatter by observing the peaks in the spectrum.
This method has proven adequate for the commonly occurring condition of fixed
 81 
spindle speed (characterized by a stationary signal spectrum) and a long duration of
chatter. Experiments on vibration data from a cut obtained while continuously
increasing the axial depth of cut, so that chatter set in and continued to increase after
the threshold depth was crossed, indicate that the spectrogram and the WD perform
equivalently in the detection of the onset of chatter. It might be expected however,
from the discussions of the earlier sections, that if the chatter were to be intermittent,
lasting for time durations that are short compared to the observation interval, then
the WD would perform better in the detection and localization of chatter compared to
shorttime spectral techniques. Another case that would benefit greatly from the
application of the WD is when the vibration signal is nonstationary, so that it is
characterized by timevarying frequency components as might happen if the spindle
speed were to vary with time. In this case the shorttime power spectrum would be
smeared due to the timevarying tooth frequency and harmonics, possibly obscuring
the chatter frequency component and making it difficult to detect.
We present here the analysis of a data sample of the acoustic signal collected
from a highspeed milling machine during a cut. Figure 420 shows a portion of the
time signal in which the axial depth of cut was gradually increased and then
decreased so that the tool entered the chatter mode for a short time duration, marked
by the increase in signal power near the center of the timeframe. The sampling rate
is 4000 Hz, and the total duration of the signal in Figure 420 is 0.38 sec. The spindle
speed was 4000 rpm, the number of teeth equal to 6 and the feed rate was set at 3000
mm/min. Figure 421 shows the 256sample spectrogram for a narrow time frame
of 0.2 sec spanning the occurrence of chatter, and the spDWD computed using a
512sample data window and 50sample timesmoothing window for the same
signal region. The development and the decay of the chatter frequency component
are clearly perceptible in both the spectrogram and the spDWD. The tooth frequency
and its harmonics are also visible. While the spectrogram and the spDWD show
 82 
2048 VALUES WD52.
1024 
1. 1L I Jill I I ,, l ,l i l ll 1 l I l U l
I I [II I 
BEG .1024 SEC MID .1536 SEC END .2048 SEC
Figure 420. Milling machine data.
 83 
time
frequency
(a)
time
frequency
(b)
Figure 421. Milling machine data (a) Spectrogram (N=256) (b) spDWD
(N=512, M=50).
 84 
similar time resolution, the latter displays better frequency resolution due to the
longer data window used. In Figures 422 and 423 we plot the
ki+Ak kitAk
quantities Z Sy(n, k) and Wy(n, k) respectively, versus timeinstant
k=kiAk k=kiAk
n, where Wy(n, k) is the 512sample spDWD with a 50sample timesmoothing
window, and Sy(n, k) is the 256sample spectrogram. The region of the frequency
summation (ki Ak, ki + Ak) is chosen to be a small fixed interval centered at each of
three distinct closely spaced frequencies including the chatter frequency. Since we
do not have a definite knowledge of the true signal spectrum, it is difficult to interpret
the performance of the spectrogram and spWD as represented by these plots. We,
however, make the following general observations. The power in the chatter
frequency component is seen to rise rapidly with time in both, but is more
concentrated in time in the spDWD. The curves at the other closely spaced
frequencies arise from the "leakage" of the chatter frequency spectrum into the
neighboring frequency bands and hence could be representative of the frequency
resolution. With this interpretation it is seen that the frequency resolution of the
spectrogram is poorer than that of the spDWD.
 85 
1050 1100 1150 1200 1250 1300 1350 1400 1450
time
Figure 422. Estimation of power in chatter frequency region (continuous
line), and two other closely spaced frequency regions (bro
ken lines) using the spectrogram.
15 
10 
5
0
1000
 86 
) 15 
5 
5 / 
1000 1050 1100 1150 1200 1250 1300 1350 1400 1450
time
Figure 423. Estimation of power in chatter frequency region (continuous
line), and two other closely spaced frequency regions (bro
ken lines) using the spDWD.
CHAPTER 5
APPLICATION OF THE WIGNER DISTRIBUTION TO VIBRATION SIGNAL
MONITORING
5.1 Spectral Analysis in Vibration Signal Processing
The vibration monitoring of machines, structures and other systems is an impor
tant process in the effective maintenance and operation of these systems. Vibration
signals emitted by machinery carry valuable information about the general condition
of the various system components. Vibration signals may be periodic or random,
stationary or nonstationary and may further contain transient phenomena and me
chanical shocks. Each component of the system contributes its peculiar spectral
signature. For example, in the case of rotating machinery the vibration signals are
characterized by high frequencies originating in rolling element bearings, gear teeth,
turbomachinery blading, etc., lower frequencies from shaft misalignments and the
essentially random contributions associated with fluid flow, jet noise and cavitation.
Spectral analysis is necessary to reveal the individual frequency components that
make up the wideband vibration signal.
.Vibration signal monitoring is typically carried out using the shorttime spectral
magnitude computed in an FFT analyzer. This method provides satisfactory results
when the signal is stationary over the duration of the analysis window. There are
situations however, when the vibration signal is nonstationary due to the timevary
ing nature of the source of vibration. Shorttime spectral techniques are then not
always adequate for the monitoring of such signals. The purpose of this chapter is to
show, by examples, what can be gained by the application of timefrequency meth
 87 
 88 
ods to the problem of vibration monitoring involving nonstationary signals. Specifi
cally, we investigate the application of the WD to the problems of envelope detection
of narrowband, timevarying signals, and the detection of periodic transients in ti
mevarying signals. It must be noted that the WD has been applied previously in
vibration signal studies [Fla89b] where timefrequency signatures of the vibration
signal were used to detect and identify abnormal machine operating conditions. Here
we propose to apply the WD in a more quantitative formulation.
5.2 Envelope Detection
The estimation of the envelope of a specific component of the vibration signal
related to a phenomenon of interest is a common operation in machinery fault detec
tion and monitoring. The observation and analysis of the vibration level at specific
frequencies, such as for example the structural resonances of the system, can indi
cate the onset or development of a fault. Typically, envelope detection is achieved
via a narrowband filter centered at the frequency of interest in order to isolate the
particular component. The filter output is rectified and lowpass filtered to obtain the
envelope, which may further be subjected to spectrum analysis, if necessary. Very
often the component of interest is nonstationary due, for example, to changes in
speed, and the above method must be modified to allow for the tracking of the time
varying frequency component. A commonly used method to achieve this is by the
use of an analog tracking filter that is adaptively tuned to the frequency of interest,
and determines the level of the signal in a small bandwidth around this frequency.
Since the WD provides a highly concentrated representation of narrowband, timeva
rying signals, its application to the envelope detection of such signals is of interest.
We review here the properties of the WD relevant to such an application.
For complex signals of the form a(t)exp(ja(t)), with dI(t) a real, smooth func
tion and a(t) slowly varying, the WD spectrum is concentrated about the instanta
 89 
neous frequency d4/dt. The average frequency of the WD, at any given time instant,
exactly equals the instantaneous frequency of the signal. Further, the instantaneous
signal power can be recovered from the WD by integrating the WD spectrum over all
frequency at the time instant of interest. That is:
00
2 Wf(t, ()d0) = If(t)2 (51)
00
The above timedomain properties holds for the pWD irrespective of the data win
dow applied. For signals of the form a(t)exp(j4(t)), (51) yields a(t)12 which is the
square of the timevarying amplitude of the signal. Although (51) indicates that the
integration of the WD must be performed over all frequency, in the case of narrow
band signals whose WD spectrum is concentrated in a narrow frequency band locally
defined around the instantaneous frequency at any time instant, the instantaneous
signal power can be obtained by integrating over this narrow subrange. This property
will be applied to the problem of the amplitude tracking of nonstationary signal com
ponents. Vibration signals are typically expected to be multicomponent and we need
to contend with the crossterms that arise from the application of the WD. By using
the analytic signal derived from the real signal, the crossterms that arise from the
interaction between positive and negative spectral components can be eliminated.
Since the intent is to recover the amplitudes of the individual signal components that
comprise the input signal by integrating the WD over narrow bandwidths in frequen
cy, care must be taken to ensure that the effect of the other signal components and
the crosscomponents on the integral is minimal. This can be achieved through the
spWD. The extent of smoothing necessary for the effective suppression of cross
terms depends on the separation of the interacting primary signal components, with
those crossterms that arise from the interaction of closely spaced signal components
needing a relatively greater amount of smoothing. However, smoothing the WD
results in the smearing of the autoterms, leading to a corresponding loss of concen
 90 
tration in the timefrequency plane. The integral over frequency of the spWD at any
time, would yield the average signal power over the duration of the smoothing time
window which would imply a corresponding loss of accuracy in the tracking of time
varying amplitudes. Hence there are tradeoffs to be considered when choosing the
smoothing window.
The above considerations suggest that the spWD can provide a powerful and
flexible tool in the amplitude tracking of nonstationary, narrowband signals. We
discuss here a specific application, that of the monitoring of vibration levels in a
highpower gas turbine. A brief overview of the application and the WDbased
algorithm for vibrationlevel monitoring are described. Simulations are used to eval
uate the performance of the method. A detailed implementation of the realtime
digital vibration monitor is presented in [Rao90b].
5.2.1 Vibration Monitoring in the Gas Turbine
The gas turbine consists of two separate rotating units coupled by air pressure.
These are the Gas Generator (GG) and the Power Turbine (PT) as illustrated in
Figure 51. Since there is no mechanical coupling between the two units their rota
tional frequencies are not the same. The available vibration signal is the output of an
accelerometer mounted on each turbine unit. Based on the analyses of these signals
the levels of vibration of each of the two units must be monitored in realtime in
order to facilitate decisions on speed, loading, and emergency shutdown of the tur
bine. The output of each accelerometer contains two dominant components one each
at the rotational frequency of the GG and the PT. These are likely to vary in time due
to changes in rotation speed. However, the frequency bands in which vibration anal
ysis is performed are nonoverlapping. The vibration signal also contains a relatively
low level of background noise contributed by other sources coupled to the turbine.
 91 
Tachometer signals are available that provides data continually on the rotational
velocities of each of the two units.
Gas
tachometer 1 generator Power tachometer 2
SGG turbine
PT
accelerometer 2
accelerometer 1
Figure 51. Gas turbine setup
5.2.2 Algorithm and Simulations
In Figure 52 the magnitude spectrum of the simulated (stationary) turbine vi
bration velocity signal from the accelerometer connected to the PT is shown. The
signal is comprised of the two dominant frequency components contributed by each
of the GG and PT, and a background made up of harmonics from various other
sources coupled to the system. The low noise floor in the spectrum is generated by
the lowpass filtering of a series of impulses. In order to demonstrate the WDbased
monitoring method, we simulate the timevarying nature of the vibration signal by
frequency modulating the PT and GG frequency components over frequency bands
specified by the manufacturer. The vibration components are specified to lie be
tween 15 hz and 70 hz for the PT, and 84 hz to 200 hz for the GG. The WD is used to
track the simulated amplitude changes of each of the vibration components in the
velocity signal in conjunction with information provided by the tachometer signal on
the actual, or instantaneous, vibration frequencies at any time. The sampling fre
quency is set at 500 hz.
 92 
magnitude
250 hz
frequency
Figure 52.
Magnitude spectrum of the turbine vibration velocity
signal.
 93 
We use the analytic signal derived from the samples of the real signal in order to
compute the discretetime spWD (spDWD). This allows us to sample the input
signal at the Nyquist rate and also eliminates crossterms in the DWD that arise from
the interaction of positive and negative spectral components of the real signal. The
analytic signal can be obtained from the input signal by using an FIR filter implemen
tation of the Hilbert transform. The Hilbert tranform is well approximated by a
noncausal and symmetric FIR filter of order 79 [Boa87]. So in order to compute the
analytic signal at each input sample we need to "look ahead" about 39 time samples.
In order to estimate the vibration amplitude in the velocity signal at a time in
stant n, the sampled velocity signal denoted by f(n) is processed as follows :
1. Select a data set { f(nN/2), ... ,f(n),..., f(n+N/21) } and compute the correspond
ing analytic signal z(n) using the FIR filtering implementation. Since the FIR filter is
symmetric and of order 79 samples, we discard 39 samples at each end of the output
data set and replace the discarded samples with zeros.
2. Generate the spDWD kernel by using the timeaveraging window g(m) of dura
tion M, and datawindow h(m) of length N. Take the DFT of the kernel, given by,
m=N/21 i=n+l1
W(n,k) = 2 [ Ih(m)2 ej 2 k/N g(i) z(i+m)z*(im) ] (52)
m=N/2 i=n
where C = is the analog frequency in radians/sec, T is the sampling period and
NT
an Npoint FFT is utilized. Figures 53 (a) and (b) show the DWD (N=512) and
spDWD (N=512,M=40) and h(n) is a Gaussian window of halfpower width 200 sam
ples) computed from the analytic signal derived from the simulated velocity signal.
We have taken g(m) as simply a rectangular window, but another smoothing function
could have been chosen so as to weight those products nearer the time instant n
more, relative to those that are far removed. The time resolution between plots is 50
samples (0.1s). The timevarying frequencies are clearly tracked by both the DWD
 94 
time .
frequency
(a)
time .
frequency
(b)
Figure 53. Effect of timesmoothing on the pDWD of the simulated
turbine vibration signal. (a) pDWD (N=512) of vibration sig
nal. (b) spDWD (N=512, M=40) of vibration signal.
