Title Page
 Table of Contents
 A review of time-frequency...
 Estimation of instantaneous frequency...
 Detection and localization of narrowband...
 Application of the wigner distribution...
 Biographical sketch

Title: Application of the Wigner distribution to problems in time-varying signal analysis
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00082240/00001
 Material Information
Title: Application of the Wigner distribution to problems in time-varying signal analysis
Physical Description: vii, 119 leaves : ill. ; 28 cm.
Language: English
Creator: Rao, Preeti, 1961-
Publication Date: 1990
Subject: Wigner distribution   ( lcsh )
Signal processing -- Mathematics   ( lcsh )
Electrical Engineering thesis Ph. D
Dissertations, Academic -- Electrical Engineering -- UF
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
Thesis: Thesis (Ph. D.)--University of Florida, 1990.
Bibliography: Includes bibliographical references (leaves 113-118).
Statement of Responsibility: by Preeti Rao.
General Note: Typescript.
General Note: Vita.
 Record Information
Bibliographic ID: UF00082240
Volume ID: VID00001
Source Institution: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: aleph - 001659139
oclc - 24529461
notis - AHX0866

Table of Contents
    Title Page
        Page i
        Page ii
    Table of Contents
        Page iii
        Page iv
        Page v
        Page vi
        Page 1
        Page 2
        Page 3
        Page 4
    A review of time-frequency distributions
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
    Estimation of instantaneous frequency and frequency-rate using the wigner distribution
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
    Detection and localization of narrowband transient signals
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
    Application of the wigner distribution to vibration signal processing
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
    Biographical sketch
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Copyright 2
        Copyright 1
Full Text








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-


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 -


ACKNOW LEDGMENTS ................. ................................... ii

ABSTRACT ......................................................... v


1 INTRODUCTION ................................................ 1


2.1 The Need for a Time-Frequency 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 ...


. . . . . 5

. . . . . . . . . 5
.................. 6
.................. 8
. . . . . . . . . 9
................. 11
.............. ... 11
................. 15
. . . . . . . . . 20

. . . . . . . . . 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 Frequency-Rate ..............
3.3.1 Derivation of a New Kernel ......... ...................
3.3.2 Computer Simulations ............................

- iii -

- iv -

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 Time-Frequency 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
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



Preeti Rao

December 1990

Chairman: Dr. Fred J. Taylor
Major Department: Electrical Engineering

The Wigner distribution (WD) is a bilinear time-frequency signal representa-
tion that possesses a number of characteristics desirable in time-varying 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
time-frequency distributions. The WD has been previously applied to the estimation
of instantaneous frequency of phase-modulated signals. Analytical expressions for
the performance of the discrete-time 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-

to-noise 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 frequency-rate of the

arbitrary phase-modulated signal. A computationally simple implementation of the

frequency-rate 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

time-frequency methods. The advantage provided by the WD over the traditional

short-time spectral magnitude is demonstrated by examples in the vibration monitor-

ing of machinery, where nonstationary signals arise due to speed variations and other

time-varying phenomena. Specifically, we consider the problems of envelope detect-
ing time-varying narrowband signals, and the detection of periodic transients of low

energy from the sidebands of a high frequency, time-varying carrier component.

- vi -


Signal analysis tools for nonstationary, or time-varying 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 time-frequen-

cy representations. Alternative signal features such as the time-extent of elementary

signal components that comprise the input signal may also be of interest, giving rise

to time-scale methods.

The short-time 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 time-scale in the case of the WT. When the scale parameter in the WT is chosen

inversely proportional to frequency, we get a time-frequency representation similar

to the STFT, but characterized by time-frequency 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 time-frequency



distributions to problems in signal processing. A number of properties that are desir-

able in a time-frequency 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 time-fre-

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 time-frequency 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 time-varying signal analy-


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

time-varying signals and systems. Some of the applications considered here have

traditionally benefitted from short-time spectral methods. The performance of WD-

based methods is investigated and compared with available techniques.

A review of time-frequency distributions is presented in Chapter 2 and an at-

tempt is made to place the WD in the perspective of time-frequency 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.


In Chapter 3 we address the problem of the joint estimation of frequency and

frequency-rate of arbitrary phase-modulated 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 discrete-time 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 frequency-rate of an arbitrary f.m. signal can be estimated. A

computationally simple implementation of the frequency-rate 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 time-frequency 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 smoothed-pseudo-WD

allows the independent control of time- and frequency-resolution in the signal repre-

sentation, and hence provides a more flexible tool in applications to nonstationary

situations than the fixed-window spectrogram. A detector based on the WD is pres-


ented that provides good time-frequency 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

short-time 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 time-varying nature of the

source of vibration. In Chapter 5 we demonstrate, by examples, the advantage pro-

vided by WD-based methods in the monitoring of nonstationary vibration signals.

Specifically, we investigate the application of the WD to the problems of envelope

detection of narrowband, time-varying signals, and the detection of periodic tran-

sients in time-varying signals.


2.1 The Need For A Time-Frequency 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 time-varying signals are speech, seismic

signals and biological signals. The phase-modulated signals of radar and sonar, as

well as the acoustic signals from rotating machinery, also belong in this class. The
analysis of such time-varying 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 time-varying spectral analy-

ses has been the short-time Fourier transform, which is implemented using either a

bank of time-windows or a bank of bandpass filters. The former is more common
and is based on decomposing the signal into short, and possibly overlapping, time



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 short-time window. Any short-time

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 short-time spec-

trum, however, the amount of time-frequency trade-off 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 "time-frequency 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 time-frequency 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 time-frequency 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 time-frequency 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 time-frequency signal energy is obtained.


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 two-dimensional probability distribution in time and frequency

of a quantum mechanical particle, defining what is now known as the Wigner-Ville

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 one-dimensional energy density functions in time and

frequency. Each of the distributions enumerated here has been derived independent-

ly from plausible time-frequency concepts and each has its own set of distinct prop-

erties that make it a good candidate for time-frequency analyses. It was realized by

Cohen [Coh66] that these time-frequency 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 time-frequency representations, a unifying framework for

several of the time-frequency 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 Choi-Williams distribution

[Cho89] and the cone-shaped kernel distribution proposed by Zhao, Atlas and

Marks [Zha90].


2.3 A Summary of Work to Date

Over the past decade there has been a renewed interest in the application of

time-frequency 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 time-frequency distributions. Much research has been di-

rected toward establishing distributions that satisfy specific properties deemed desir-

able in time-frequency analyses and in their interpretation. The quality of the esti-

mates of time-frequency 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 inter-relationships between the different

distributions [Coh89]. Efficient digital implementations have been developed
[Wil87]. Time-frequency 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 time-varying charac-

teristics of the signal or ii) the use of a particular property of the time-frequency

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

time-evolution 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, time-frequency representations
have been applied to the classification of muscle sounds [Bar90] and to event-related


potentials [Cho87]. These signals are inherently nonstationary and time-frequency
distributions provide consistent patterns rich in detail, making the classification easi-
er than with standard time-domain or frequency-domain 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
time-frequency 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 two-dimensional 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/
spatial-frequency resolution, the WD also has the advantage that it implicitly en-
codes phase in a single real-valued 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 time-frequency representations has been estab-
lished by Cohen [Coh66] which provides a unifying framework for a number of dis-
tributions that have been proposed for time-frequency 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 (2-1)

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 time-frequency shift invaria-

nce that characterizes the Cohen's class of distributions. That (2-1) provides a rea-

sonable definition of a general time-frequency 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]
Cf(t, a<; 1) = K(t, T)e-Ydr (2-2)
00 00
K(t, T) = J f f(u + r/2)f*(t T/2)P(4, r)eJ()t-u)dudQ (2-3)
-00 -00
We see that the generalized autocorrelation function is a time-averaged and lag-

weighted value of the time-dependent 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 time-frequency repre-

sentations are obtained. Based on the requirements of time-varying signal analysis

several properties desirable of a time-frequency 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 time-frequency 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-


- 11 -

2.5 The Wigner Distribution

The WD belongs to the Cohen's class of bilinear time-frequency distributions. It

is defined for a signal f(t) by,
Wf(t, w) = f (t + r/2)f*(t r/2) exp(- jor)dr (2-4)
and is obtained from (2-1) 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 frequency-shifts of the

signal are reflected as corresponding shifts in the WD time-frequency plane, so that

for g(t) = f(t -to), we get

Wg(t, ) = Wf(t- to, 0) (2-5)

For g(t) = f(t) exp(jot), we get

Wg(t, ) =Wf(t, o O) (2-6)

The WD is real-valued, 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,
Sf Wf(t, o)do = If(t)|2 (2-7)

- 12 -

fWf(t, )dt = F(o) 12 (2-8)
The WD preserves the time- and frequency-support of the signal, so that

f(t) = 0 for |tj > T => Wf(t, w) = 0 for jt]>T. (2-9)

F(O))=0 for |wo > Q=>Wf(t, o)=0 for |o|>Q0 (2-10)

However the WD is not necessarily zero at all time-instants and frequencies that the

signal goes to zero [Coh87].

The first-order frequency moment of the WD is given by,

,(t)12 JWf(t, w)do) = O(t) (2-11)

Qf)(t) = Im = Im-dln f(t) (2-12)
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. (2-12),

Qf(t) = '(t) (2-13)

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


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 (2-14)
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 (2-15)
--00 -CO -00
This property has been used to formulate the optimum detection of signals in the
time-frequency plane. [Fla88].

The WD preserves convolution and modulation. If g(t) = f(r)h(t-r)dr then

Wg(t, o) = J Wf(r, )Wh(t- -', w)dT (2-16)
For g(t)=f(t)h(t),

Wg(t, 0) = Wf(t, n)Wh(t, () nj)dri (2-17)
In practice the WD is implemented after applying a finite-length, moving window to
the data, centered at the time-instant of interest, giving rise to the pseudo-WD
(pWD). For a signal f(t) and real, symmetric window h(t) the pWD is computed as

pWf(t, ()= f(t+,/2)f*(t-r/2)h( )h*(--)exp(-jwor)dr (2-18)

The pWD is related to the WD by

pWf(t, 0) J Wf(t, q)Wh(0, ( 7)dj (2-19)


- 14-

Wh(0, ) = h2(r/2) exp(-jwr)dr (2-20)
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 slowly-varying, real

function, and 0(t) is a smooth, real function, so that the signal can be considered to

be phase-modulated with a slowly-varying envelope, the WD is concentrated about

the curve 0'(t) and has significantly lower spread compared to any other time-fre-

quency distribution. In the case of linear frequency-modulated signals, the WD at
any time-instant is a delta function at the instantaneous frequency. An enlightening

interpretation of this general behavior of the WD for phase-modulated 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*(t-T/2). The bandwidth of
the short-time 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 2-1.

- 15 -

Vj Vi

0 0
(a) (b)

Figure 2-1. Effective bandwidth (shown by shaded area) of windowed
frequency-modulated signal. (a) STFT (b) pWD

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 2-1 it is evident that for a given window, the

ideal situation is now only a linear frequency modulation approximation (or more

generally, a skew-symmetric 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 time-varying
frequencies. In any case, the instantaneous frequency can always be recovered ex-
actly from the first moment of the WD as given by (2-13).

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 time-frequency 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))] (2-21)

Wf,g(t, C) = ff(t + T/2)g*(t /2) exp(- jwr)dT (2-22)
is the cross-WD of f(t) and g(t). Thus the WD contains cross-components arising

from the interaction of the distinct time-frequency 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 (2-21) where

Wf(t, 0) = |aI2Wh(t tf, ) (Of)

Wg(t, o) = Ib2Wh(t tg, 0 cg) (2-23)

and the cross-term 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 ) (2-24)

Thus in addition to the auto-components given by the corresponding shifted versions

of Wh(t, 0), the resultant WD contains a cross-component 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 time-frequency 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 time-frequency plane. Figure 2-2 illustrates these results.



Figure 2-2. WD of a bicomponent signal [Hla84].

It has been found that the above simple geometrical laws on the nature of the
cross-term 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 time-frequency,
and hence is present even in monocomponent signals where the various time-shifted
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 short-time 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,

cross-terms arise from the interaction of each signal component with the remaining
window components not corresponding with it.

- 18 -

The cross-components 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 time-frequency plane.

Cross-terms 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 signal-to-noise 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 two-dimensional smoothing in

the time-frequency 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 (2-25)
-00 -00
where 1I(t, o) is the chosen smoothing function that is localized in the time-frequency

plane and is sufficiently regular, but otherwise arbitrary. Since the convolution in

(2-25) effects a lowpass filtering in the time-frequency plane, the highly oscillatory

cross-components are attenuated while the auto-components are generally smeared.

The optimal smoothing of the WD in terms of eliminating cross-terms while retain-

ing a highly concentrated distribution, would require adapting the smoothing func-

tion in (2-25) to match the local time-frequency characteristics of the signal. An

example of such an approach applied to unimodular signals has been presented

[And87]. A special case of the time-frequency smoothing represented by (2-25) 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 smoothed-pseudo WD (spWD) [Fla84a],

- 19 -

spWf(t, w) = f h(r) f g(t u)f(u + 2)f*(u )du e-Jdr (2-26)
-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 time-frequency analyses. The time-frequency separable

smoothing function, however, restricts the form of the smoothing region to be paral-

lel to the t- or f-axis.

In general, time-frequency smoothing of the WD leads to the sacrifice of the

marginal and local properties. The exponential kernel of the Choi-Williams' distri-

bution is an example of time-frequency smoothing where the two smoothing func-

tions are related in a specific manner that allows the marginal and local properties to

hold, while the cross-terms are suppressed by their spreading out in the time-fre-

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 cross-terms 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 non-overlapping fre-

quency-time components is free of cross-terms, 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 time-fre-
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 cross-compon-
ents but with concentrated auto-components. 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 cross-terms 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
short-time Fourier transform given for a signal f(t) and window h(t) by

Ft(w) = J f(r)h(r t) exp(- jor)dr (2-27)

The spectrogram is then

Sf(t, ) = IFt(()|2 (2-28)

The spectrogram is a member of the Cohen's class of bilinear time-frequency
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 (2-29)

00 OD
S(t,w ) =-- f (, r) exp(j [t-rw])drd (2-30)
-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 (2-1) by
the two-dimensional Fourier transform of (2-30). 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 (2-31)

We see that the spectrogram is a time-frequency 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 time-fre-
quency plane of the spectrogram is always slightly greater than that of the WD
[Nut88]. The time-frequency extent of the spectrogram is minimized when the win-
dow function is matched to the signal. Further, the time-frequency 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 (2-32)
-00 -00
00 00
J Sf(t,)dt j IF(r)121nH(w )2d (2-33)
-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 (2-26), the time-smoothing and frequency-smooth-

ing functions applied to the WD in (2-31) to get the spectrogram, are not indepen-

dent but rather, are related through the predetermined window h(t). This results in

the time-frequency trade-off that is characteristic of the spectrogram. However, the

amount of time-frequency smoothing effected by Wh(t, cw), the WD of the window

function, is sufficient to eliminate negative values and cross-terms 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 non-overlapping 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 time-width, 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 cross-term 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 frequency-averaging 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 time-frequency plane. Since the WD provides a maximally

concentrated distribution and the spectrogram is related to it via the time-frequency

averaging of (2-31) it seems obvious that a window function whose WD is most

concentrated in the time-frequency plane would achieve the minimal spread. The

Gaussian window is such a function with its WD being a double Gaussian function in

the time-frequency plane. In order to achieve optimal concentration, the parameters
of the window (namely its time-/frequency-width and orientation in the time-fre-

23 -

quency plane) must be matched to those of the signal components. Such an approach

is embodied in the data-adaptive transform of [Jon87] where the Gaussian window

parameters are adapted to minimize a concentration measure at each location in



3.1 The Concept of Instantaneous Frequency

The need for the analysis of frequency-modulated (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) (3-1)

where 9(t) is the Hilbert transform of s(t) :

1 ( s(T)
9(t) = (Hs)(t) = v.p. (t- dr (3-2)

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 (3-3)

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) (3-4)
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 (3-1), 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)) (3-5)

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 (3-4) have a meaningful interpretation as the IF of the signal.

3.2 Estimation of Instantaneous Frequency from the WD

The use of a time-frequency 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 time-frequency 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 (2-13). The discrete-time WD (DWD) is defined by

Ws(n, f) = 2 1 s(n+k)s'(n-k) exp(-j4;rkf) (3-6)

The average frequency at any time-instant n is computed as,
fi(n) =- arg[ f exp(j44rf)W,(n, f)df ] (3-7)

- 26 -

For a discrete-time complex signal given by,

s(n) = a(n) exp(je(n)) (3-8)

where a(n) and j(n) are real functions, (3-7) evaluates to [Cla80b]

fi(n) = [(n- 1)] mod 2 (3-9)

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
complex-valued 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. (3-9) holds exactly
for the pseudo-DWD computed using any (real-valued) 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 (3-8),
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) (3-10)

of a complex sinusoid in noise, where z(n) is complex white Gaussian noise (w.g.n.)

with zero-mean and variance r 2, it has been shown [Tre85] that for a high enough

input signal-to-noise ratio (SNR) A2/a (3-10) may be approximated by,

s(n) = A exp(j [2nrfn + + w(n)])


- 27 -

where {w(n)}, the phase noise sequence, is real w.g.n. Hence for such an input signal
the estimate of (3-9) 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 time-frequency 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 (3-6),
equivalent to the best least-squares fit of a complex sinusoid to the kernel sequence
{s(n+k)s* (n-k)} 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 phase-modulated
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)


- 28 -

where z(n) is complex w.g.n. with zero-mean and variance or z. The DWD of s(n)

is given by (3-6). In practice the DWD is evaluated over a finite data record of N

samples and is given by

Ws(n,f)=2 I s(n+k)s'(n-k)exp(-j4nkf) (3-13)

with N=2L+1. The DWD can be looked upon as the frequency-scaled discrete-time

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 least-square's sense, the kernel

sequence computed over the finite interval {-L+1,L-1}, at the time-instant n. For the

signal of (3-12) 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'(n-k)

+A exp(-j[ (n-k)])z(n+k)+z(n+k)z*(n-k) (3-14)

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 (3-14) yields the constant frequency

sinusoid A2exp(j8rajnk). The remaining three terms in (3-14) 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+ (3-15)

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

SNRDWD = 2A2oz + (3-16)

The problem is now reduced to that of estimating a constant frequency sinusoid in
white noise with SNR given by (3-16) using the peak of the DFT magnitude. It is well
known that the peak of the DFT is a maximum-likelihood estimate (MLE) of the
frequency of a sinusoid in w.g.n. At high enough input SNR this MLE meets the
Cramer-Rao 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 (3-17)

where z(n) is complex w.g.n. with zero-mean and variance a 2, the variance of the
MLE of frequency at high input SNR is given by

varDFT(f) = (3-18)
(2a)2(SNR) (N2 1)

where SNR in (3-18) stands for the SNR in the DFT of the noisy sequence of (3-17)

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 (3-16), into (3-18) 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) (3-19)

At high input SNR (A2 > o2) the above reduces to

varDwD(f) = (2x)2A2N(N2 1) (3-20)

which equals the variance of the MLE of the frequency of a stationary sinusoid in

w.g.n., as given by (3-18).

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 (3-16). Hence the threshold SNR for the DWD peak estimator is given by

SNRDWD = A4xN/2)= 15dB (3-21)
2A2o +o -

which at high input SNRs (A2 > oj) reduces to the condition,

.- = 21dB (3-22)

- 31 -

Since the DWD is always real-valued, 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 (fi-fi)2 (3-23)

where M is the number of trials, fi is the actual IF and ft the estimated value. A

64-sample data window was used and a large, zero-padded 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 3-1

compares the simulation results with the theoretical MSE given by (3-19). 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 skew-symmetry within

- 32 -

(MSE)-1 dB

0 2 4 6 8 10 12


Figure 3-1.

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 least-square'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 time-varying frequency-rate that varies continuously over a range

of values from zero to a maximum. In Figure 3-2 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 time-instants 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 mid-frequency, a point about

which the IF curve is skew-symmetric 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 trade-off between bias and variance of the estima-


- 34 -


1 3 5 7 9 11 13 15

Figure 3-2.

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 Frequency-Rate

The joint estimation of instantaneous frequency and frequency-rate of phase-
modulated signals is of importance in problems such as the orbit determination of
satellites, motion compensation in synthetic-aperture radar processing, and the high
resolution signature analysis of targets. In all of these situations, the received signal
is the return from a doppler-type radar system, and the frequency and its time-deri-
vative are important observation parameters. The joint estimation of doppler and

doppler-rate 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-

cy-rate have been proposed. These include using a linear least-squares fit estima-

tion of these parameters on a finite data set of frequency samples [Kno85], the maxi-
mum-likelihood estimation of the instantaneous frequency and frequency-rate of the
parameters of a linear chirp signal [Aba86], and the estimation of frequency-rate 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 frequency-rate 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 frequency-rate 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 frequency-rate 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) (3-24)

where q(t, r) is the WD kernel of the signal s(t) and is given by,

q(t, ) = s(t+-)s (t--) (3-25)

We show here that for the input signal given by s(t) = exp(j (t)) the two-dimen-

sional Fourier spectrum of r,(t, T1, T2) given by R,(t, 0i, 0)2) satisfies the following


Sf 0122Rs(t, C01,02)dowidCO2 = C x 0"(t) (3-26)
-00 -00
where C =-j(2n)2. For s(t) = exp(jo(t)), we get from (3-24),

r,(t, Ti,T2) = exp(j 1(t, Ti,T2)) where,

01(t, rI,r2) = ,(t + + ) (t + -)

(t ) + (t- 2 ) (3-27)


rs(t, Tr, T2) = R,(t, ow1,2) exp(j [w t1 + (2r2])d(old2 (3-28)
-00 -0
we have,

62-rt, T1, V2) - 1 [2
6t162 o 4o= f \ 2Rs(t, ox,w2)dwldw2 (3-29)
-00 -00

But from (3-27) we get,

- 37 -

62rs(t, r1, r2) 1r1=O,r2=o = jp"(t) (3-30)

The result (3-26) then follows from (3-29) and (3-30).
The implementation of (3-26) involves obtaining the two-dimensional 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 frequency-rate. We give the simpli-
fication here.

If one of the parameters T, or T say T1, is fixed at a value T, then the one-di-

mensional DFT of rs(t, T, r2) with respect to T2, denoted by R&,T(t, (02) can be

shown, by retracing the steps (3-27) to (3-30), to satisfy,
limr--o j0 2Rs,T(t, (02)do2 =
= limT-+o [0'(t +,-) '(t --)] = 0"(t) (3-31)

We can rewrite the above expression to obtain an estimate of the frequency-rate as,

T (2Rs,T(t, (02)dw2 p 0"(t) (3-32)
with the approximation becoming increasingly accurate as the interval T is made
smaller. For signals characterized by a purely quadratic frequency modulation,
(3-31) holds without the limit, and hence the expression (3-32) holds exactly, irre-
spective of the choice of the time-interval T.
We can simplify computation further, by basing the estimate of frequency-rate
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 discrete-time estimator based on the

- 38 -

frequency moment of Rs,T(t, (02) as in (3-32) would be more sensitive to noise than

one based on the peak of the same quantity. For f.m. signals for which the frequen-
cy-rate 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 frequency-rate equals 6at)

the spectrum Rs,T(t, W2) is a delta-function 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 frequency-rate. In

analogy with the behavior of the WD peak in the estimation of IF, the estimate of

frequency-rate 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 frequency-rate versus time curve. Nevertheless, we have now an
order-of-magnitude improvement in the tracking of time-varying frequency-rate
over methods that are based on the assumption of constant frequency-rate within the
data window.

3.3.2 Computer Simulations

In this section we investigate the performance of the frequency-rate estimator
based on the peak of the new transform for a discrete-time 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) (3-33)

which is a complex, linear f.m. signal in additive, complex w.g.n. of zero-mean and

variance ao 2, so that the SNR is equal to A2/o The actual frequency-rate 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 discrete-time by,

Rs,T=2(n, f) = > q(n+ 1,k)q'(n- 1,k)exp(-j4uakf) (3-34)

q(n,k) = s(n + k)s'(n k) kE (-L+1,L-1) (3-35)

It is seen that Rs,T(n, f) can be computed easily as the (frequency-scaled) DFT of the

product of DWD kernels separated in time by T. For the signal of (3-33) the kernel

of (3-34) is equal to exp(j4zaok) and hence the peak of (3-34) with respect to f

provides an unbiased estimate of the frequency-rate.

Figure 3-3 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-

quency-rate estimator of the signal (3-33). This lower bound is given by [Dju89]

var(d) = 2A2N( 1)(N 4) (3-36)

We see that the MSE of the frequency-rate 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 (3-34)

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 (3-34).

- 40 -




(MSE)-1 dB

7 8 9 10 11 12 13 14


Figure 3-3.

Performance of the frequency-rate estimator for a
linear f.m. signal in additive w.g.n. (N=64);
101ogio(1/MSE) versus input SNR.



-: C 3LB

x : si ulatior result


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 quasi-harmonic 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 time-frequency representation which

is closely related to the short-time 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,
y(t) = I Cm,n g(t- na) exp(j2mafit) (4-1)

- 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 time-scale variations of the signal. This

problem of the variability in detector performance with changes in time-scale of the

transient signal has been addressed through the use of time-scale 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

time-time-scale 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 time-scale. 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 short-time 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 trade-off determined by the window


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 time-frequency 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 non-overlapping 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 time-scale changes of the transient waveform, so that the detector

maintains a reasonable level of performance over a wide range of transient time


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

short-time power spectrum and the Wigner distribution are among the possible

time-frequency 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 time-frequency plane than does the spectrogram.

The extent of time-frequency 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 time-frequency plane, which for nonstationary input signals

can lead to a considerable simplification in detector structure. In this chapter, we
address the application of time-frequency 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 Time-Frequency 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 pseudo-WD, however, allows the increase of

data window size to achieve high enough frequency resolution without a sacrifice in

. time localization. The difference in combined time-frequency resolution between

- 47 -

the spectrogram and the pWD is illustrated by Figures 4-1 to 4-3. Figure 4-1 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 short-time 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 frequency-scaled (by a factor of two) version of

the actual signal DWD.) Spectrograms obtained with two different data window

lengths and spanning the time-instant of arrival of the transient are compared in

Figure 4-2. (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 time-averaging window used

to compute the spDWD. We use rectangular data and time-averaging 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 4-3 which shows the spDWD of the signal using the same two data

windows, does not display any perceptible trade-off 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 cross-terms in the WD on the

detection of the transient waveform when the received signal contains multiple,

closely spaced frequency components. The cross-terms 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 time-smoothing of the WD in order to improve the

frequency resolution by attenuating the cross-terms. This operation leads to a

corresponding loss of time resolution. However, since the time-smoothing window is

- 48 -



BEG .05120 SEC

MID .1280 SEC

ENO .2048 SEC

Figure 4-1. Gaussian transient in sinusoid.


wDn i



cfU I40 .

- 49 -





Figure 4-2. Effect of data window length on spectrogram of Gaussian
transient in sinusoid. (a) N=128 (b) N=512.

- 50 -





Figure 4-3. 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 trade-off between

time and frequency resolutions as in the short-time 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 cross-term interference,

choose the appropriate length for the time-averaging window. In most cases a
time-smoothing window that is much shorter in duration than the data window will

be adequate to attenuate the cross-terms 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 4-4 to 4-6 illustrate this point
well. Figure 4-4 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 4-5 and 4-6 show respectively the spectrogram and the

spDWD of the waveform in Figure 4-4, 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 4-6 shows the 512-sample spDWD of the signal with a

time-smoothing window length of 40 samples, in order to attenuate the cross-term
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 -




BEG .02560 SEC

1i ll

MID .1024 SEC

END .1792

Figure 4-4. Gaussian transient in sum of two closely spaced sinusoids.

VOr- W03i.

-edU4 I ---- -------- I ----- I ----- I _____ I I I I_________ | __



I I I i

- 53 -


__AFigure 4-5.

Figure 4-5.


Spectrogram of Gaussian transient in sum of two closely
spaced sinusoids (N=84).

- 54 -



Figure 4-6. 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 pulse-width 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
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 time-marginal property of the WD we have for a signal u(t),

I Wu(t, f)df = u(t)2 (4-2)
We can thus recover the instantaneous signal power at any time by integrating the
WD over all frequency, at that time. Equation (4-2) 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(r-t)j2dr (4-3)
-00 -00

Figure 4-8 compares the robustness of the time-marginals of the spectrogram and
the WD, computed using fixed data windows, to variations in the time duration of the

received transient signal. Figure 4-7 shows a time waveform containing three

- 56 -

BEG .05120 SEC MID .1280 SEC END .2048 SEC

Figure 4-7. Sum of Gaussian transients of same peak amplitude, but dif-
fering pulse durations.

- 57 -

60 -







Figure 4-8.

1000 1200 1400 1600 1800


Time-marginals of the spectrogram and the WD a compar-
ison; --- spectrogram; WD.


- 58 -

transient signals of the same center frequency and peak amplitude but different time
durations. Figure 4-8 compares the quantities (4-2) and (4-3) 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 (f-Af,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 Wu(t, f)df = |u(t) (4-4)
Integrating the WD in frequency provides the further advantage of smoothing
out cross-components between multiple transient signals that are closely spaced in
time. Signal components that are separated in time give rise to cross-terms in the
WD midway between the two interacting components. These cross-terms are
oscillatory in the frequency direction and, therefore, are suppressed by averaging the
WD in frequency, an effect achieved by (4-4). Figure 4-9 shows the time waveform
of two almost overlapping transient signals. Figure 4-10 shows the 200-sample
spectrogram and the 512-sample DWD, for this signal. Cross-terms are clearly seen
in the DWD and in this case, due to the overlap between signal components, in the
spectrogram as well. Figure 4-11 compares the quantities of (4-2) and (4-3) for this

signal, where the DWD is computed using a 512-sample data window and the

- 59 -

BEG .05120 SEC MID .1152 SEC END .1792 SEC

Sum of two closely spaced in time Gaussian transients.

Figure 4-9.

- 60 -





Figure 4-10. Cross-terms in the representation of overlapping Gaussian
transients. (a) Spectrogram (N=200) (b) pDWD (N=512).

- 61 -




3 30\


10 -" -,

900 950 1000 1050 1100 1150 1200 1250


Figure 4-11. Time-marginals of the spectrogram and the WD of overlap-
ping transients a comparison; --- spectrogram; WD.

- 62 -

spectrogram using a 128-sample 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 4-12 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

Window threshold 1p.
v(t) Detector

Figure 4-12. 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 finite-length window v(t) and then the WD is computed,
incorporating time-smoothing 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,
r(t, f)= J Wy(t, v)dv (4-5)
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 (4-5) 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 time-frequency resolution characteristic of the

WD-based detector via the example of the time-series of Figure 4-13. The signal is

comprised of two short-duration Gaussian pulses in an additive background created

by a closely-spaced (in frequency) narrowband component with slowly-varying
amplitude. Figure 4-14 shows the detection statistic of (4-5) 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 4-15 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 WD-based

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 -


1000 -





600 800 1000 1200 1400 1600 1800


Figure 4-13. Transient pulses in a narrowband background with slowly
varying amplitude.

- 65 -




S 5

0 4




600 800 1000 1200 1400 1600 1800


Figure 4-14. The WD detection-statistic (N=512, M=128) at the frequency
of the transient signal.

- 66 -






2 -


0 -------------

600 800 1000 1200 1400 1600 1800


Figure 4-15. 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 short-time 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) (4-6)

where u(t) is the transient signal waveform, and n(t) is zero-mean, stationary,
complex Gaussian noise that satisfies,

n(t)= 0 and n(tl)n(t2) = 0. (4-7)

The noise power spectral density is given by

Gn(f) = n(t + ) n*(t ) exp(- j2arfr) dr (4-8)

To compute the WD, we must window the received signal. The windowed waveform
is given by,

y(t) = v(t) [u(t) + n(t)] (4-9)

where v(t) is the chosen, finite-duration window function. The WD is then calculated

- 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 +-)

u*(t -2)

exp(- j2;rfr) dr

exp(- j2rfr) dr

n*(t -) exp(-j2nfr) dr

u 2

exp(-j2nfr) dr

Rw(t, r)=v(t+) v'(t- -)
2 2


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)]


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

var[Wy(t, f)] = 2W,(t, f) Ga(f) + 2 Gn(v) IB(t, f 12 dv
(NxN) (SxN)









- 69 -


Gnn) Gn(f) *G(f) (4-18)

B(t, f) Wv(t, f -) U(v) exp(j2n t) dv (4-19)

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 (4-17) 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

time-windowing of the data, then Wv(t, f) = 6(f), and the NxN term in (4-17)

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 4-17 compares the 512-sample DWD with the 256-sample spectrogram

of the noisy signal of Figure 4-16, 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 short-time Fourier transform for transients with time

durations that are short with respect to the observation window. For the sake of

- 70 -

4096 VALUES WO31.



4095 VALUES W032.


BEG 0.0 SEC MID .1280 SEC END .2560 SEC

Figure 4-16. Gaussian transient in white, Gaussian noise (Eo/NO = 7.6).

- 71 -



!M ijAy JA^^A-/-AM< W-A^^^^AA^~^lf^ ^f^ ~ ^


Figure 4-17. 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 (4-6) where,

u(t) = a exp(- t2/2oI) (4-20)

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 4-12, where

v(t)= exp(- t2/2p2) (4-21)

H(f) = exp(- f2/2B2) (4-22)

The filtered noise spectrum is therefore,

Gn(f) = NoIH(f)12 = No exp(- f2/B2) (4-23)

The quantities required to calculate the mean and the variance in (4-16) and (4-17)
for this example are given below.

U(f) = aoo 2 exp(- o2xr2f2) (4-24)

W(t, f)= 2a20ro I exp(- t2/o2 oO42f2) (4-25)

W(t, f) = 2p r exp(- t2/p2 p242f2) (4-26)

Gmn(f) = N2 B exp(- 2f2/B2) (4-27)

- 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 (4-29)
[var(Wy(t, f))] II

with the variance computed as in (4-17). The expression (4-29) 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
4-18 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 (4-30)

The mean signal output at t=0, f=0 is then,

Sy(0, 0) = U(f) V(f) df (4-31)

The variance of the output noise is,

var[Sy(0, 0)] = Gn(f) IV(f)2 df (4-32)

- 74 -



1 --7


Figure 4-18.

1 2 3 4

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 (4-24) and (4-23) respectively, and V(f) is the
spectrum of v(t),

V(f) = p exp(- p 222f2) (4-33)

Substituting the above expressions in (4-28) 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 (4-34)
[N0(ol + p2)]1/

From Figure 4-18 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 time-frequency 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 4-12 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 (f-Af,f+Af). At time t and frequency f, this quantity
may be written in the following equivalent form,
r(t, f) = JWy(t, v) dv = Wy(t, f) *f D(f) (4-35)

- 76 -


D(f) = 1 for f E (-Af,+Af)
=0 else (4-36)

But y(t) = v(t)x(t) and hence we have,

Wy(t, f) = Wv(t, f) *f Wx(t, f) (4-37)

Substituting in (4-35) we get,

r(t, f) = Wv(t, f) *f Wx(t, f) *f D(f) (4-38)


r(t, f) = Weq(t, f) *f Wx(t, f) (4-39)

Weq(t, f) = Wv(t, f) *f D(f) = JWy(t, v) dv (4-40)
From (4-39) we see that the effect of integrating the WD of the received signal over

the frequency interval (f-Af,f+Af) is equivalent to a modification of the window v(t)
so that its WD is now given by (4-40). For the window given by (4-21), this
expression evaluates to

1 t2
Weq(t, f) = 1exp(-) [erf(2np(f + Af)) erf(2arp(f- Af))] (4-41)

where erf refers to the error function and is given by

erf(y) = 2 exp(- t2)dt (4-42)

- 77 -

The mean detector output when the signal is present is then obtained from (4-16) by
replacing W,(t, f) with Weq(t, f) to get

El [r(t, f)] = Weq(t, f) *f [Wu(t, f) + Gn(f)] (4-43)

The mean detector output when the signal is absent is

Eo[r(t, f)] = Weq(t, f) *f Gn(f) (4-44)

The variance of the detection statistic when the signal is present is obtained from
(4-17) 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 (4-45)
(NxN) (SxN)

B(t, f) = Weq(t, f ) U(v) exp(j2nrvt) dv (4-46)

The variance of the detection statistic with only noise present is simply the first term
(NxN) of (4-45). 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 (4-21) and (4-22), the above expression was
evaluated at the time of arrival, t=0, and center frequency, f=0, of the transient signal
u(t) of (4-20), 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 4-19 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 4-19 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 4-19 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 (4-35), 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

(4-43), 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 (4-45) 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 (4-45) is constant since it does not depend on the pulse

- 79 -

SNR 4-

3- 4




Figure 4-19.

SNR of the WD-based 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 4-19. 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 high-speed 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 short-time 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

short-time 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 time-varying frequency components as might happen if the spindle

speed were to vary with time. In this case the short-time power spectrum would be

smeared due to the time-varying 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 high-speed milling machine during a cut. Figure 4-20 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 time-frame. The sampling rate

is 4000 Hz, and the total duration of the signal in Figure 4-20 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 4-21 shows the 256-sample spectrogram for a narrow time frame

of 0.2 sec spanning the occurrence of chatter, and the spDWD computed using a

512-sample data window and 50-sample time-smoothing 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 4-20. Milling machine data.

- 83 -





Figure 4-21. 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 4-22 and 4-23 we plot the

ki+Ak kitAk
quantities Z Sy(n, k) and Wy(n, k) respectively, versus time-instant
k=ki-Ak k=ki-Ak

n, where Wy(n, k) is the 512-sample spDWD with a 50-sample time-smoothing

window, and Sy(n, k) is the 256-sample 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


Figure 4-22. Estimation of power in chatter frequency region (continuous
line), and two other closely spaced frequency regions (bro-
ken lines) using the spectrogram.

15 -

10 -



- 86 -

) 15 -

5 -
5 -/ ------------------

1000 1050 1100 1150 1200 1250 1300 1350 1400 1450


Figure 4-23. Estimation of power in chatter frequency region (continuous
line), and two other closely spaced frequency regions (bro-
ken lines) using the spDWD.


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 short-time 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 time-vary-
ing nature of the source of vibration. Short-time 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 time-frequency 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, time-varying signals, and the detection of periodic transients in ti-
me-varying signals. It must be noted that the WD has been applied previously in

vibration signal studies [Fla89b] where time-frequency 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, time-va-

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:
2 Wf(t, ()d0) = If(t)|2 (5-1)
The above time-domain properties holds for the pWD irrespective of the data win-

dow applied. For signals of the form a(t)exp(j4(t)), (5-1) yields |a(t)12 which is the

square of the time-varying amplitude of the signal. Although (5-1) 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 cross-terms that arise from the application of the WD. By using

the analytic signal derived from the real signal, the cross-terms 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 cross-components 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 cross-terms 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 auto-terms, leading to a corresponding loss of concen-

- 90 -

tration in the time-frequency 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 trade-offs 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
high-power gas turbine. A brief overview of the application and the WD-based
algorithm for vibration-level monitoring are described. Simulations are used to eval-
uate the performance of the method. A detailed implementation of the real-time
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 5-1. 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 real-time 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 non-overlapping. 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.

tachometer 1 generator Power tachometer 2
SGG turbine

accelerometer 2
accelerometer 1

Figure 5-1. Gas turbine setup

5.2.2 Algorithm and Simulations

In Figure 5-2 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 WD-based

monitoring method, we simulate the time-varying 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 -


250 hz

Figure 5-2.

Magnitude spectrum of the turbine vibration velocity

- 93 -

We use the analytic signal derived from the samples of the real signal in order to
compute the discrete-time spWD (spDWD). This allows us to sample the input
signal at the Nyquist rate and also eliminates cross-terms 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
non-causal 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(n-N/2), ... ,f(n),..., f(n+N/2-1) } 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 time-averaging window g(m) of dura-
tion M, and data-window h(m) of length N. Take the DFT of the kernel, given by,

m=N/2-1 i=n+-l1
W(n,k) = 2 [ Ih(m)|2 e-j 2 k/N g(i) z(i+m)z*(i-m) ] (5-2)
m=-N/2 i=n-

where C =- is the analog frequency in radians/sec, T is the sampling period and
an N-point FFT is utilized. Figures 5-3 (a) and (b) show the DWD (N=512) and
spDWD (N=512,M=40) and h(n) is a Gaussian window of half-power 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 time-varying frequencies are clearly tracked by both the DWD

- 94 -

time .


time .


Figure 5-3. Effect of time-smoothing 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.

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs