THE ANALYSIS OF NON-STATIONARY SIGNALS USING TIME-FREQUENCY METHODS

Journal of Sound and Vibration (1996) 190(3), 419–447

THE ANALYSIS OF NON-STATIONARY SIGNALS USING TIME–FREQUENCY METHODS J. K. H  P. R. W Institute of Sound and Vibration Research, University of Southampton, Southampton SO17 1BJ, England (Received 1 November 1995) This paper deals with the time–frequency analysis of deterministic and stochastic non-stationary signals. It includes the following: a brief review of the fundamentals of time–frequency analysis and various time–frequency distributions; a summary of the inter-relations between time–frequency distributions, emphasizing links between evolutionary spectra, the Wigner–Ville distribution, and Cohen-class distributions; demonstration of the effects of using different analysis methods on sample time histories; the analytical modelling of physical phenomena so as to predict theoretical time–frequency distributions, with a view to offering insight into the signal generation mechanisms. 7 1996 Academic Press Limited

PRELIMINARY

The first author’s interest in time–frequency methods was initiated when a research student in 1966. Brian Clarkson had obtained a preview of a paper on the spectral analysis of non-stationary signals by Maurice Priestley to be published in the Journal of Sound and Vibration in 1967 and said, ‘‘See what you can make of that!’’. It is somewhat sobering to realize that this author is still responding to this suggestion. Time–frequency analysis is a pre-occupation that is hard to resist. This paper is an attempt to review the many approaches to time–frequency analysis of non-stationary processes and is inevitably coloured by the ISVR experience over the past three decades, with some emphasis on results from Ph.D. theses of students from 1980 to the present. The second author is also an ISVR graduate—although of a much more recent vintage— who is able to present a more contemporary perspective on a topic that, in the past 15 years or so, has erupted with new theory, computational techniques and widespread applications. 1. INTRODUCTION

A non-stationary signal is one the ‘‘structure’’ of which changes as a function of time. We will define ‘‘structure’’ in due course. For now we will indicate what is meant by listing a few examples—birdsong, speech, noise and vibration signals from accelerating traffic, Doppler shifted sound, chirp signals, vibration response of machinery undergoing operational changes, the impulse response of a damped non-linear Duffing oscillator, trajectories of a Lorenz system, velocity of an impacting oscillator—the list is endless. In each case, it is quite likely that most people would agree that the phenomena exhibit ‘‘non-stationarity’’ even though we might not be too specific about what form this takes. 419 0022–460X/96/080419+29 $12.00/0

7 1996 Academic Press Limited

420

. .   . . 

Frequency analysis of time histories now occupies a central role in the vast subject of (digital) signal processing and, together with the apparent natural ability of human beings to carry out ‘‘local’’ frequency analysis, this motivates a more formal study of the characterization of non-stationary processes in time and frequency. Analytical methods that address non-stationarity by decomposing signals over two dimensions have been developed. Such decompositions may lead to ‘‘double-frequency’’ spectra or the more physically meaningful time–frequency spectra. The latter allow the ‘‘internal structure’’ of a time history to be examined in a way that the individual time and frequency domains do not. A variety of such characterizations exists, broadly categorized as linear representations (Gabor forms, evolutionary spectra, and also including the time– scale forms of the wavelet transform) and quadratic representations (centred on the Wigner–Ville distribution and the associated Cohen-class of distributions). In essence, the advantages are that the time–frequency decompositions permit projections of the time history on to a space that allows separation of components of the signal so as to facilitate enhancement, detection (of events), filtering, classification (of phenomena) and resynthesis. This paper will include the following: a brief review of the fundamentals of time–frequency distributions, a summary of the interrelations between time–frequency distributions; emphasizing links between evolutionary spectra, the Wigner–Ville distribution, and Cohen-class distributions; demonstration of the effects of using different analysis methods on sample time histories; the analytical modelling of physical phenomena so as to predict theoretical time–frequency distributions, with a view to offering insight into the signal generation mechanisms, and also to support the (empirical) data analysis of time series.

2. TIME–FREQUENCY METHODS FOR NON-STATIONARY SIGNALS

2.1.  A non-stationary stochastic signal is a signal the statistical structure of which changes as a function of time, e.g., the mean, variance, correlation (covariance), etc., may change as time evolves. Specifically, we may say that the joint probability density relating the values of x(t) at times t1 , t2 , . . . , varies under a shift in time. In view of the success that spectral analysis of stationary time series has enjoyed there is, naturally, considerable interest in extending this to the non-stationary case to consider spectral descriptions that are functions of both frequency and time. Signals that are not random may display a form that can also be described as non-stationary: e.g., a chirp signal, for which we could refer to a ‘‘changing frequency’’. In the past years there has been considerable activity in discussing time–frequency distributions for such signals. In this paper, the word ‘‘nonstationary’’ is taken to include both stochastic and non-stochastic signals. The aim is to demonstrate the use of several time–frequency decompositions on a variety of signals. Non-stationary random processes having a ‘‘frequency modulated’’ form often arise in the physical world. Examples include descriptions of the motion of a vehicle over rough terrain as it accelerates over the surface and the acoustical case of the sound perceived by an observer as a source travels past the observer. The second example is more complex than the first as, in addition to the Doppler effect, it also includes range and directionality effects. Some of the material presented here is a summary of aspects of studies carried out by the authors and their co-workers over a number of years on developing time–frequency methods for this class of processes. The first author’s interest in the field was initiated by Priestly [1, 2], which led to applying the method to problems in acoustics and vibrations [3].

– 

421

Time–frequency descriptions are very attractive to the engineering world and empirical analysis of non-stationary processes using ‘‘moving window spectral analysis’’ has been commonplace for a long time. There has been an upsurge of interest in a variety of descriptions other than the conventional spectrogram and the evolutionary spectrum. Accordingly, this presentation will begin with a very brief review of aspects of signal analysis fundamental to time–frequency analysis. 2.2.   2.2.1. The Uncertainty Principle Fourier analysis provides a decomposition of a time history in the frequency domain, i.e., the pair†

g

x(t)= ejvtX(v) dv,

and

X(v)=

1 2p

g

x(t) e−jvt dt

(1a, b)

are interpreted as follows. In equation (1a), the signal x(t) is considered to be made up of sines and cosines (which last forever) having (complex) amplitudes 0X(v) dv, the amplitude density X(v) being computed from equation (1b). Introducing both time and frequency into the picture simultaneously to address a question such as ‘‘When does frequency component v make its dominant contribution to the signal x(t)?’’ does not have an immediate answer. Indeed, this question seems at best to be very vague when one considers the ‘‘Uncertainty Principle’’ [4, 5] briefly summarized here. Define

g

>x >2= x 2(t) dtQa,

= t

1 >x>2

g

tx 2(t) dt,

(Dt)2=

1 >x>2

g

(t−t)x 2(t) dt.

>x >2 is the energy of the signal, t is the centre of gravity of the area under x 2(t) and Dt is a measure of the dispersion in time of the signal. Similarly,

g

>X >2= =X(v)=2 dv,

g

v¯ = v =X(v)=2 dv,

(Dv)2=

1 >X >2

g

(v−v¯ )2=X(v)=2 dv.

Starting from these definitions, one can show that Dv.Dte12 :

(2)

i.e., the bandwidth Dv (or B), time width Dt (or T) product, i.e., BT, is never less than 1/2. In essence, the Uncertainty Principle (equation (2)) says that the more precisely one specifies ‘‘when’’ for a signal (i.e., T reducing) the less specifically can one specify what frequencies are involved. It is important to note that often x(t) refers to the analytical signal x˜ (t) (see the next section) and is complex valued, so the absolute value =x˜ (t)= is used in the definitions leading up to the Uncertainty Principle. Accordingly, the bandwidth–time product of a signal will differ depending on whether it is the real valued signal x(t) or its analytic counterpart that 2 is used. For example, if x(t) is a real Gaussian pulse, say e−at , then the BT product takes 2 the minimum value of 1/2. But, if x(t) is e−at cos v0 t, the BT product increases. However, 2 if the analytic version is used, i.e., x˜ (t)=e−at ejv0 t , then the BT product reverts to 1/2. † Throughout this work, integrals without limits will be taken to be between −a and a.

422

. .   . . 

2.2.2. Instantaneous frequency and group delay The group delay of a signal x(t) is obtained from the Fourier transform X(v) as follows. Express X(v) as =X(v)= ejargX(v); then the group delay, tg (v), is tg (v)=−(d/dv) arg X(v).

(3)

This may be interpreted as the time at which frequency component v makes its maximum contribution to a signal [6]. This is, therefore, an attempt to link the concepts of time and frequency based completely on the Fourier transform. An alternative approach uses the concept of the Hilbert transform of a signal as follows. The instantaneous frequency is defined for x(t) by first forming the analytic signal, i.e., x˜ (t)=x(t)+jxˆ (t)=A(t) ejf(t) ,

(4)

where xˆ (t) is the Hilbert transform of x(t). (The Hilbert transform of x(t) is formed by convolving x(t) with 1/pt.) The ‘‘envelope’’ of the analytic signal is A(t) and the ‘‘instantaneous phase’’ is f(t). The instantaneous frequency, vi (t), is vi (t)=f (t).

(5)

The two apparently very different relationships, equations (3) and (5), link time and frequency for the signal x(t) and one might ask if they are essentially the same—specifically whether they are inverse functions of each other. In general they are not—and the papers [7–9] describe conditions that must be met by a signal x(t) for these two time–frequency laws to be equivalent. These conditions are that the signal should be ‘‘asymptotic’’, i.e., that the BT product should be ‘‘large’’ and that the signal should be ‘‘monocomponent’’, i.e., vi=f should be invertible. The concept of monocomponent is not rigorously defined but is intuitively clear. Two chirps each with different instantaneous frequency added together lead to a signal that is not monocomponent. Computation of an overall f would result in some form of average. So, if a signal is monocomponent and asymptotic, then the instantaneous frequency and the group delay are equivalent. We note that it is the analytic signal that is commonly used in most time–frequency analyses. The ‘‘mathematical’’ reason for this is because the Fourier transform of x˜ (t) satisfies 2X(v), vq0 F[x˜ (t)]= X(v), v=0 . 0, vQ0

8

9

The fact that it is zero for negative frequencies accounts for the BT product being 1/2 for 2 the example above (i.e., x˜ (t)=e−at ejv0 t ). It also turns out to be particularly helpful when computing the Wigner distribution (cross-terms are reduced). Another reason for the importance of the analytic signal is that it leads to an unambiguous definition of instantaneous frequency which often has immediate physical significance. In reference [8], the analytical example of a Gaussian pulse modulated chirp signal is given showing how the group delay and instantaneous frequency forms are the same when the BT product is large and how different they are when it is not. A verbal description of what a large BT product means [8] is that the (relative) amplitude variation is small with respect to the (relative) frequency variation, or (refer to equations (4) and (5)). A /A(d/dt)(f )/f ,

– 

423

2.3.    - In this section we shall give a brief overview of the various approaches to the definition of time–frequency spectra. We emphasize that we are concerned with fundamental definitions here and omit any significant discussion of estimation methods. A common starting point is the concept of time-varying periodogram, obtained by sliding a window across a time record, then performing a Fourier transform to obtain the Short Time Fourier Transform (STFT)

g

S(t, v)= x(t')w(t−t') e−jvt' ) dt'.

(6)

A method of displaying this is to form =S(t, v)=2 and plotting this as a function of v and t. The squared magnitude of this quantity is called the periodogram or spectrogram and this three-dimensional plot (crudely) shows how the energy of the process is decomposed over time and frequency. This is an intuitively appealing procedure and is often very helpful. It must be stressed, however, that it is empirical in the sense that it will change for each window w used and so if it is regarded as an ‘‘estimate’’ of something one must ask the question ‘‘An estimate of what?’’, which implies the existence of some underlying fundamental concept. It is to this problem that the majority of this section is addressed.

2.3.1. Linear representations This section is based on the idea of representing a signal in terms of a set of other (simpler) components which when added up yield the original signal. In this sense the Fourier series might be regarded as the simplest form (relating to periodic signals). We shall include discussions relating to both deterministic and stochastic signals, for which the concept of non-stationary means that we must modify the Fourier form in some way and yet we still wish to retain the concept of frequency.

2.3.1.1. Harmonizable processes Let us begin with stationary random signals. The Fourier approach means that we can write (see references [1, 2])

g

x(t)= ejvt dZx (v).

(7)

This is a linear representation of x(t) as a weighted sum of sines and cosines. If the process is stationary, then the increments dZx (v) are orthogonal, i.e., E[dZ* x (v1 ) dZx (v2 )]=0,

v1$v2 ,

(8)

which leads to the definition of the power spectral density, S(v), where

g

E[x 2(t)]= S(v) dv,

and

S(v) dv=E[= dZx (v)=2 ].

(9, 10)

424

. .   . . 

From the linear representation, it follows that the autocorrelation function R(t) is related to the spectral density S(v) by

g

R(t)=E[x(t)x(t+t)]= S(v) ejvt dt,

(11)

with the inverse S(v)=

1 2p

g

R(t) e−jvt dt.

(12)

These well-known results are basic to what follows. The linear representation of equation (7) leads to the squaring of dZx (v) to yield a power spectral density. However, that same power spectral density could have been obtained via a Fourier transform of the autocorrelation function. We shall refer to this latter route to the definition of a power decomposition in section 2.3.2, under the heading ‘‘Quadratic forms’’. What happens when x(t) is non-stationary? For a class of non-stationary processes representation (7) is still valid (so-called harmonizable processes) but now dZx (v) is no longer orthogonal, so that we obtain a double-frequency spectrum, i.e., S(v1 , v2 ) dv1 dv2=E[dZ* x (v1 ) dZx (v2 )],

(13)

with the corresponding relationship E[x(t1 )x(t2 )]=R(t1 , t2 )=

gg

S(v1 , v2 ) ej(v2 t2−v1 t1 ) dv1 dv2 .

(14)

The quantity S(v1 , v2 ) is a generalization of the spectral density for stationary processes but lacks physical appeal. That is not to say that it is not useful or widely used (see, for example, references [10, 11]). In contrast, however, a physically meaningful two-dimensional spectrum whose dimensions are time and frequency was proposed by Priestley [1, 2] and is called the evolutionary spectral density. 2.3.1.2. The evolutionary spectral density The essence of Priestley’s idea is that the basic building blocks in the representation, namely, sines and cosines ejvt are replaced by amplitude modulated sines and cosines A(t, v) ejvt , so that

g

x(t)= A(t, v) ejvt dZx (v).

(15)

For the class of processes x(t) admitting this representation, dZx (v) is again orthogonal, and this leads to the definition of a time-varying spectral density termed the evolutionary spectral density, Se (t, v)==A(t, v)=2S(v),

(16)

where S(v) is the ‘‘usual’’ spectral density arising from dZx (v). The class of processes for which this is valid was termed ‘‘oscillatory processes’’ by Priestley. The essential requirement is that A(t, v) be ‘‘slowly varying’’ (as a function of t for each v) compared to ejvt .

– 

425

Note that this representation was presented as applying to stochastic signals. Extensions to cross-spectra and estimation methods have been described [12, 13], but no explicit formalism for deterministic signals based on this has been given.

2.3.1.3. The Gabor representation In a classic paper [4], Gabor considered a representation of signals (not only stochastic) in terms of basic building blocks that were (localized) Gaussian pulse modulated sines and cosines of the form w(t−nD) e jvk t .

(17)

These basic elements have the lowest BT product achievable. The w(t) are the Gaussian pulses (reference [14] permitted w(t) to have a more general form). The signal x(t) was considered to be a summation of these: i.e., x(t)=s s cnk w(t−nD) e jvk t n

(18)

k

The index n positions the pulse in time and the index k relates to the frequency. The coefficients cnk are distributed over the time (n)–frequency (k) plane. The non-orthogonal nature of the basic components makes the calculation of the cnk ’s difficult [14, 15]. The representation equation (18) is discrete in ‘‘cells’’ over time and frequency and continuous (integral) representations were given, for example in [16, 17]. In reference [18] it was attempted to link these two-dimensional decompositions with Priestley’s representation.

2.3.1.4. The short-time Fourier transform It is interesting to note that the short-time Fourier transform is closely related to the above. If one considers the STFT of x(t), as defined by equation (6), then it can be shown that for any finite energy function v(t), where

g

v(t)w(t) dt=Cs$0

(19)

then

x(t)=

1 Cs

gg

S(t', v)v(t−t') ejvt' dv dt'

(20)

Indeed, in this context the Gabor representation can be considered as a sampled version of the STFT, where the sampling is two-dimensional and occurs in the time–frequency plane. It is important to note that it is not possible to create a Gabor representation in which the basis functions, defined in equation (17), are both orthogonal and have good localizing properties in the time and frequency domains [19]. However, wavelet transforms, discussed in the next section, can be constructed which do possess both of these desirable properties [19].

426

. .   . . 

2.3.1.5. Wavelets [19] Recently, there has been considerable interest in wavelet transforms. Here we discuss wavelets within the context of an alternative to time–frequency methods. The continuous wavelet transform is defined as

g

W(t, a)= x(t)w((t−t)/a)* dt;

(21)

Like the STFT, the wavelet transform creates a function of two variables from a time history. In the case of wavelets the time history is decomposed as a function of time and scale, a. From (21) we see that the wavelet transform is calculated by forming inner products of the time history with versions of a second function w(t). This function is dubbed the mother wavelet. In fact the conditions on w(t) under which the representation (21) is valid are only mild: for well-behaved functions one only has to ensure that the mother wavelet has zero mean. This ensures that the mother wavelet must oscillate in some manner, and since it should also have finite energy then it will also decay towards 2a. For the purposes of this paper we restrict our choice of mother wavelet to functions of the form w(t)=h(t) ejv0 t .

(22)

Many of the wavelet functions found in the continuous wavelet literature are of this form, largely because (as we shall see) they most easily relate to the concept of frequency. We further assume that h(t) has the same properties as a spectral windowing function. The complex form of equation (22) means that for sufficiently large v0 then w(t) is approximately analytic. The choice of using an analytic wavelet is made for the same reasons that we use the analytic form of x(t) in other areas of time–frequency analysis. It is unimportant which function in equation (21) is analytic, but it is advantageous to ensure that at least one of them is. The effect of the scale parameter, a, is similar to the inverse of frequency. Consider the scaling operation alone, i.e., consider w(t/a), then w(t/a)=h(t/a) ejv0 t/a=h(t/a) ej(v0 /a)t .

(23)

From this we see that the effect of a is two-fold. First, the centre frequency has become v0 /a, allowing us to assign a specific frequency to each value of a. Second, the windowing function is modified to h(t/a); such an operation serves to stretch (or compress) the analyzing window. Thus for large values of a, w(t/a) has a short duration and high frequency, while for small a it has long duration and low frequency. Usually, only the squared modulus of the wavelet transform is considered; this is called the scalogram. This relationship mimics that of the STFT to the spectrogram. From the wavelet transform one can recover the original time history by using x(t)=

1 Cg

gg

1 za

W(t, a)w((t−t)/a)

da dt . a2

(24)

where the constant Cg is defined by Cg =

g

=W(v)=2 dv. =v=

(25)

– 

427

The requirement that Cg be finite is the only constraint on the mother wavelet function. The reconstruction formula (24) shows that the wavelet transform can be thought of as decomposing the signal into elements of the form w((t−t)/a). The primary difference between the wavelet transform and the STFT is the nature of the set of elementary functions. In the case of the STFT each function has a constant envelope so that as the frequency changes the duration of the elementary function is fixed. This should be contrasted with the wavelet transform, where the elementary functions vary their duration in accordance with the centre frequency. In both cases the BT products of all the elementary functions are the same, but in the wavelet case the bandwidth reduces at lower frequency so one maintains a constant proportional bandwidth or, equivalently, all the functions have a constant Q factor. The advantages of the wavelet transform lie in that for certain classes of signal such a constant Q analysis is more appropriate than the fixed bandwidth analysis offered by the STFT. Indeed, the human ear performs an analysis which is approximately constant Q. Wavelets also allow the construction of orthonormal bases with good time frequency resolution. Such representations seek to represent a signal in as few points as possible. However, in signal analysis problems it is commonly the case that over-representing the data is beneficial in which case this property of the wavelet transform is only of secondary importance. 2.3.2. Quadratic forms This section is concerned with what is sometimes referred to as quadratic representations or the ‘‘energetic’’ point of view on time–frequency methods. The STFT and the wavelet transform permit the energy of the signal to be displayed in either the (t, v) or the (t–a) plane: i.e., for the spectrogram, 1 2p

gg

g

=S(t, v)=2 dt dv=Ex= =x(t)=2 dt,

(26)

and for the scalogram

gg

=W(t, a)=2

dt da =Ex . a2

(27)

Both admit linear representations of the signal; however, two-dimensional energy decompositions do not in general require this. 2.3.2.1. The autocorrelation form As indicated in section 2.3.1.1, the spectral density of a signal may be obtained from the autocorrelation function, i.e., S(v)=

1 2p

g

R(t) e−jvt dt,

(28)

so that

g

R(0)=E[x 2(t)]= S(v) dv.

(29)

For a non-stationary signal the concept of a ‘‘local’’ autocorrelation function is logical: i.e., one that quantifies the behaviour of the average of products x(t1 ) x(t2 ) for t1 and t2 ‘‘in the vicinity’’ of time t.

. .   . . 

428

We defer (briefly) the exact definition of the local autocorrelation and express the local dependence as R(t, t); then a corresponding time-dependent spectral density might be defined as

g

SQ (t, v)= R(t, t) e−jvt dt

(30)

(the subscript referring to quadratic). This will have the property of an ‘‘energetic’’ description if SQ (t, v) describes the energy (or, in the stochastic case, power) of the process over the t–v plane such that Ex =

gg

SQ (t, v) dt dv.

(31)

A variety of approaches to the definition of the local autocorrelation may be conceived. Let us approach this intuitively. For a deterministic signal the localization at time t of the average of the product x(t1 ) x(t2 ) may be considered as follows. Let t2−t1=t, the usual ‘‘lag’’ variable, and form the product x*(u−t/2)x(u+t/2).

(32)

Now ‘‘concentrate’’ the values of this product (for fixed t) near time t by weighting it with a function that ‘‘peaks’’ near t, i.e., g(u−t), and then add up all such products. Specifically, one forms

g

g(u−t)x*(u−t/2)x(u+t/2) du.

(33)

This is a candidate for the local autocorrelation R(t, t) and one can envisage a variety of functions g that lead to greater or lesser concentration ‘‘near’’ time t. One can introduce greater generality by letting g() also vary with the lag t, so that we write

g

R(t, t)= g(u−t, t)x*(u−t/2)x(u+t/2) du.

(34)

For the stochastic case we must average this, so that

g

R(t, t)= g (u−t, t)E[x*(u−t/2)x(u+t/2)] du.

(35)

In the following sections we will consider particular cases of this. 2.3.2.2. The Wigner distribution Let us take the very special case of the ‘‘ultimate’’ concentrating function g above, namely a delta function, such that g(t)=d(t). Then for the deterministic case

g

R(t, t)= d(u−t)x*(u−t/2)x(u+t/2) du=x*(t−t/2)x(t+t/2)

(36)

– 

429

and

g

Sw (t, v)= x*(t−t/2)x(t+t/2) e−jvt dt.

(37)

For the stochastic case we use

g

Sw (t, v)= E[x*(t−t/2)(t+t/2)] e−jvt dt.

(38)

This particular time–frequency distribution has come to play a central role and is known as the Wigner–Ville distribution, or simply the Wigner distribution. Its origins are in quantum mechanics in 1932 (see the review [20]). However, its importance in signal analysis was not recognized until much later; see, for example, references [21–23]. It is possible to show that the double integral of the Wigner distribution over time and frequency yields the energy and in this sense it is an energy density function. However, while it has some characteristics that might appear disconcerting, let us begin with a very attractive property. Consider the analytic form of the chirp signal, i.e., 2

x(t)=ej(at+bt /2)

(39)

f =a+bt.

(40)

having instantaneous frequency

The Wigner distribution for this can easily be shown to be Sw (t, v)=d(v−(a+bt)):

(41)

i.e., it picks out the instantaneous frequency law ‘‘perfectly’’. It is tempting to think that this will be mirrored for other signals—but it turns out that it is only for a single (monocomponent) linear chirp that it ‘‘works’’ so well. Departures from this lose the ‘‘sharpness’’ in the t–v plane. Its more disconcerting characteristics, however, relate to two other properties. Despite being an energy density distribution it is negative for some values. In fact, but for the example of the linear chirp, the Wigner distribution always has some negative values. Interestingly this property (which is a natural concomitant of its definition) has not inhibited its acceptance and use. It is perhaps important to note that Priestley’s evolutionary spectral density is never negative and so differs fundamentally in this respect. The second disconcerting attribute arises from the existence of ‘‘cross-terms’’. This is a consequence of the bilinear nature of its definition. For example, if x(t) is composed of two components, i.e., x(t)=x1 (t)+x2 (t), then the Wigner distribution of x, Sw (t, v), is Sw (t, v)=S1 (t, v)+S2 (t, v)+S12 (t, v)+S21 (t, v).

(42)

S1 (t, v) and S2 (t, v) are the auto-Wigner distributions, while S12 (t, v) and S21 (t, v) are cross-Wigner distributions. This results in unwelcome non-zero values in the t–v plane. These terms tend to be located mid-way between the two auto terms and are oscillatory in nature. The result is a time–frequency distribution in which there appear components which are not directly related to a single term in the input signal. There is always the possibility that as a result of these cross-terms incorrect inferences about the input signal can be drawn.

. .   . . 

430

Words such as ‘‘confusing’’ and ‘‘spurious’’ to describe the cross terms, together with negative energy distribution, seem to throw doubt on whether it is worth using the Wigner distribution. However, its very central (pivotal) role in time–frequency analysis negates many of its disadvantages. This is expanded upon in what follows. In calculation of the Wigner distribution one is always faced with a limited length of data and so the infinite integrals in equation (37), (38) have to be replaced by finite integrals or, equivalently, a windowing function employed. Thus in practice one (nearly) always calculates a so called pseudo-Wigner–Ville distribution (PWD), defined as Sp (t, v)=

g

T/2

h(t)x*(t−t/2)x(t+t/2) e−jvt dt,

(43)

−T/2

which is equivalent to the choice of g(t, t)=d(t)h(t) in equation (34). The function h(t) is defined on [−T/2, T/2]. The importance and ‘‘popularity’’ of the Wigner distribution has been increasingly recognized over the past decade, but it is worth noting that substantial work, see e.g., references [24, 25], has been reported on the so-called ‘‘instantaneous spectrum’’ for non-stationary stochastic processes many years ago. The instantaneous spectrum is essentially equation (38). 2.3.2.3. The Cohen class of distributions In a series of papers, Cohen generalized the definition of time–frequency distributions in such a way as to include a wide variety of different distributions (see the review [20]). These different distributions can be represented in several ways; here we elect to define Cohen’s class as the Fourier transform, with respect to t, of the generalized local correlation function: i.e.

g

Sg (t, v)= R(t, t) e−jvt dt=

gg

g(u−t, t)x*(u−t/2)x(u+t/2) e−jvt dt du.

(44)

and the function g(t, t) is referred to as the kernel function, the choice of which defines the different distributions. Some of the more common choices for kernel functions are illustrated in Table 1

T 1 Definitions of kernels for various time–frequency distributions Distribution

g(t,t) d(t)

Wigner–Ville

d(t)h(t)

Pseudo-Wigner distribution (PWD) Spectrogram Choi–Williams [26]

Cone-kernel [27]

Reduced interference distribution [28]

h(t−t/2)h(t+t/2) 1 z4pt/s

6

exp

$ % 7 −st 2 4t 2

h(t),

=t =Q=t=/2

0,

elsewhere h(t/t) =t =

– 

431

Several of these distributions have a further degree of ‘‘user choice’’; namely, one is free to select a windowing function h(t). For the constraints on the choice of h(t) the reader is referred to the individual references. Also for the Choi–Williams distribution there is a parameter s which allows control of the degree of smoothing. It should be noted that the spectrogram is in fact a member of Cohen’s class, even though the STFT is not (since the STFT is a linear, rather than quadratic, representation). The spectrogram also has a unique role in Cohen’s class since its is the only distribution which is guaranteed to be non-negative for all signals. Note that for specific signals other distributions can be non-negative; for example, we have seen that the Wigner distribution for a linear chirp is non-negative. 2.3.2.4. Signal recovery/reconstruction/synthesis The linear representations of section 2.3.1 obviously allow reconstruction of the signal if the time–frequency decomposition is known, but it is less obvious whether energetic decompositions allow recovery of the signal. For the Wigner distribution the signal can be recovered uniquely to within a constant phase factor [29] by using x(t)=

1 2p

g

Sw (t/2, v) ejvt dv.

(45)

(It is appropriate to note here that the volume [29] is a very useful and informative collection of papers relating to time–frequency analysis in general.) A problem of perhaps greater physical significance is that of synthesizing a signal from a prescribed time–frequency distribution. One might envisage filtering in the t–v plane, say, by forming [30, 31] S (t, v)=G(t, v)Sw (t, v),

(46)

where G(t, v) is some ‘‘mask’’ attempting to remove multicomponents, interference, etc. The question is whether S (t, v) is valid. The papers [30] and [31] address this by using least squares methods based on time–frequency eigendecompositions. The essence is to find a signal y(t) with time–frequency decomposition S(t, v) that best approximates the modified time–frequency distribution as that which minimizes min y

gg

=S(t, v)−S (t, v)=2 dt dv.

(47)

It is important to realize that simple use of the inverse formula (20) with S (t, v) or some chosen construction is an ad hoc way of synthesizing a signal. However, it is noted that practical requirements often result in very effective procedures for generating signals which are useful. Methods for simulating non-stationary processes for earthquake type and other phenomena having required time–frequency spectra is discussed in reference [32]. 2.3.3. Interlinking of distributions The apparent diversity of approaches to time–frequency analysis may leave practitioners with a somewhat insecure feeling. This situation has been significantly improved in the recent past [20, 33–36] with some very illuminating publications, some results of which are summarized here. As noted earlier, the Wigner distribution plays a pivotal role.

. .   . . 

432

Specifically, if x(t) has a Wigner distribution Sw (t, v) and if Sg (t, v) denotes a (general) distribution belonging to the Cohen class, then Sg (t, v)=

gg

Sw (t, n)P(t−t, n−v) dt dn,

(48)

where P(t, v) is an arbitrary function of time and frequency. The interpretation of this remarkable result is clearer if P(t, v) is ‘‘low pass’’ in the time–frequency plane, so that the Cohen class may be considered as a smoothed Wigner distribution. This process is equivalent to treating the original Wigner distribution as an image and performing image processing to try to suppress the cross-terms. Such a strategy relies on the fact that the interference terms oscillate and so tend to occupy the higher spatial frequencies. As a special case, consider the spectrogram using an analyzing window, of, say h. Then it turns out that =S(t, v)=2=

gg

Sw (t, n)Hw (t−t, n−v) dt dn:

(49)

i.e., where Hw (t, v) is the Wigner distribution of the windowing function and plays the role of the P(t, v) in this case. Thus the spectrogram is the two-dimensional convolution of the Wigner distributions of the signal and the windowing function. Smoothing of this nature is referred to as ‘‘regular’’ smoothing. An associated remarkable result uses the concept of affine (linearly scaled) smoothing. This result states that the scalogram is related to the Wigner–Ville distribution by affine smoothing: i.e., =W(t, a)=2=

gg

Sw (t, n)Hw

0

1

t−t , an dt dn, a

(50)

where Hw (t, v) now represents the Wigner distribution of the mother wavelet function. So, like the spectrogram, the scalogram can be obtained by combining (affine smoothing) two Wigner distributions, in this case those of the signal and the mother wavelet. Thus, with the Wigner distribution as the core distribution along with the concepts of regular and affine smoothing, we can obtain a unified perspective on a wide class of distributions. However, the relationship between the members of the Cohen class and the evolutionary spectral density (for stochastic processes) does not conform to this, but is amenable to an alternative analysis. It was shown in reference [37] that if Sw (t, v) is the Wigner distribution, and Se (t, v) is the evolutionary spectral density of a process with modulating function A(t, v), then W(t, n)=

1 2p

g

V(t, t, v)Se (t, v) dv,

(51)

where

g

V(t, n, v)= A*(t−t/2, v)A(t+t/2, v) e−jt(n−v) dt/=A(t, v)=2:

(52)

i.e., V(t, n, v) is (for each v) a normalized Wigner distribution of the modulating function.

– 

433

A more general version of this was given in references [38, 39], where it was shown that if Sg (t, v) denotes a member of the Cohen class, then

g

Sg (t, v)= Se (t, v')

VA (t, v−v') dv', =A(t, v')=2

(53)

where VA is the Cohen class distribution for A(t, v). Note that these results are a smoothing type expression over v for fixed t, in contrast to the two-dimensional regular and affine smoothing forms. 2.3.4. Cyclostationary processes and evolutionary spectra We note briefly the current interest in so-called cyclostationary processes. There is an extensive bibliography on the subject (see, for example, references [11, 40]). An example of an cyclostationary process is a non-stationary stochastic process which is periodically correlated; i.e., E[x*(t1 )x(t2 )]=R(t1 , t2 )=R(t1+T, t2+T).

(54)

If t1=t−t/2 and t2=t+t/2, then we can write R as Rt (t)=Rt (t+T); i.e., for fixed t, Rt is periodic in t. Both double-frequency spectra and time–frequency spectra are relevant to this class of processes. In references [11, 41], helicopter noise is modelled in this form, and also described is a general class of processes called Correlation Autoregression (CAR) processes which include cyclostationary processes. Specifically such processes have covariance functions that satisfy a linear relationship of the type N

R(t1 , t2 )= s aj R(t1+tj , t2+tj ).

(55)

j=1

In reference [38] it was shown that if x(t) is a CAR process described by equation (55), then the corresponding evolutionary spectral density is also auto-regressive related: i.e., N

Se (t, v)= s aj Se (t+tj , v).

(56)

j=1

Similar relationships have been developed [39] for the Cohen class of distributions. 2.3.5. Desirable properties of time-frequency distributions We have examined various distributions and considered some of their attributes and interrelations. In many presentations what might be considered as desirable properties of time–frequency distributions are often listed at the outset—we shall present a shortened form now (see Table 2) and indicate whether some distributions satisfy these or not. Additionally, some time–frequency distributions and their properties are summarized in Table 3, and three distributions are compared in Table 4. We note that various publications, e.g., the volume [29], and reference [41], provide tabulations of this nature.

3. APPLICATIONS OF TIME–FREQUENCY METHODS

This section and the next are focused on examples of problems and data occurring in acoustics and vibration. The aim is to show how the methods described earlier are relevant to these. We include aspects of estimation and theoretical considerations.

. .   . . 

434

T 2 ‘‘Desirable’’ attributes of time–frequency distributions No.

Property attribute description

Mathematical description

P1 P2

Shift in time and frequency

x(t−t0 ) \ Sg (t−t0 , v) x(t) ejv0 t \ Sg (t, v−v0 )

P3

Reality

Sg (t, v)=Sg (t, v)*

P4 P5

Finite support in time and frequency

If x(t)=0 for =t=qT, then Sg (t, v)=0 for =t=qT If X(v)=0 for =v=qV, then Sg (t, v)=0 for =v=qV

P6

Marginal over time

1 2p

P7

and frequency

P8

Group delay

P9

Instantaneous frequency

P10

Positivity

g g g

g

Sg (t, v) dv==x(t)=2

Sg (t, v) dt==X(v)=2 tSg (t, v) dt

>g >g

vSg (t, v) dv

Sg (t, v) dt=tg (v) Sg (t, v) dv=vi (t)

Sg (t, v)e0 for all t and v

T 3 Some time–frequency distributions and their properties Properties 1. 2. 3. 4. 5. 6. 7. 8. 9. 10.

Shift in time Shift in frequency Reality Finite support in time Finite support in frequency Integration over time=energy density spectrum Integration over frequency=instantaneous power Average frequency at time t=instantaneous frequency Average time at frequency v=group delay Positivity

Wigner–Ville Spectrogram Pseudo-Wigner–Ville Choi–Williams Z Z Z Z Z

Z Z Z × ×

Z Z Z Z ×

Z Z Z × ×

Z

×

×

Z

Z

×

Z

Z

Z

×

Z

Z

Z ×

× Z

× ×

Z ×

3.1.   Time–frequency distributions are used as a tool for signal analysis: i.e., they are used to gain insight into the structure of a signal. The results of such an analysis are often (although by no means exclusively) interpreted by a human operator, and hence judging what defines a good time–frequency method may be highly subjective. To illustrate the application of time–frequency analysis to measured signals we present an example based on structural response. This signal would generally be classed as deterministic. In

– 

435

T 4 Comparison between spectrogram, evolutionary spectrum and pseudo-Wigner–Ville distribution Property Types of signal Associated ‘‘frequency’’ Reality Positivity Uniqueness No interference terms if q one component Fourier spectrum if stationary input Time–frequency resolution Relationship to Wigner–Ville Average ‘‘frequency’’ Average time at given ‘‘frequency’’

Spectrogram

Evolutionary spectrum

Deterministic/random Random Fourier Related to Fourier Z Z Z Z Z × Z Z Z

Z

Limited by Uncertainty Principle

Limited by Uncertainty Principle

Spread version in time and frequency Averaged frequency over window Averaged group delay over window

Theoretical integral relationship ? ? ?

PWD Deterministic/random Instantaneous Z × Z × Z If frequency resolution increased, time resolution unchanged Spread version in frequency

Averaged group delay over window

time–frequency analysis, signals are generally classed as deterministic, not if they can be completely described via mathematical expressions, but merely if they essentially consist of a relatively small number of narrow-band (time varying) components. Each of these narrow band components may be adequately modelled as a random process, e.g., a randomly excited, time-varying, narrow-band filter. Even so, they are often treated as deterministic signals. Of course, for a sufficiently narrow-band process, of which only one realization is available, distinguishing been a random and a deterministic process is a difficult and, in general, only an academic, task. 3.1.1. Structural response Structural responses are dispersive in nature; i.e., different frequency components travel through the system at different speeds. In acoustic problems dispersion is commonly caused by geometry of the propagating medium: e.g., in a duct different modes propagate at different rates. This results in system impulse responses which can be gainfully examined with the aid of time–frequency analysis. The paper [42] is an interesting application using the Wigner distribution for the identification of structure-borne noise components. Our aim here is to use a similar (but simpler) situation to illustrate some of the general aspects of time–frequency analysis that we have discussed above. The example presented here consists of a beam with one end freely suspended and an anechoic termination at the other. The beam is excited by an impulse, and acceleration is measured at a single point on the beam. The signal at that point consists of two components: the first is the direct path from the impulse to the measurement point, and the second is due to the reflection of the travelling wave from the free end returning to

436

. .   . . 

Figure 1. The impulse resonse of the transverse vibration of a simple beam

the measurement point. A time history from one such experiment can be seen in Figure 1, where these components are clearly in evidence. To explore the internal structure of this signal a spectrogram, using a Gaussian window, is calculated initially. Care has to be taken to choose the duration of the window. While one would like to retain a reasonable window length to ensure good frequency resolution, if the window is too long then at their closest point the two signal components will be merged: a manifestation of the loss of temporal resolution. In Figure 2 is shown a spectrogram in which the window length has been selected to produce a visually appealing result. The amplitude of the spectrogram is displayed by using a linear grey scale, dark reflecting higher values.

Figure 2. A spectrogram of the acceleration using a Gaussian window.

– 

437

The structure of the signal is now clear. Each component consists of a frequency modulated narrow-band process. Closer examination reveals that this curve closely fits the v−1/2 law predicted by simple Euler–Bernoulli theory [43]. The second component has travelled a greater distance and so has dispersed more, resulting in greater curvature of the time-frequency law. There is, however, considerable width (spreading) to each of the curves in the time–frequency plane; consequently it is difficult to accurately define the curve which describes the time–frequency laws. The goal of most time–frequency methods is to ‘‘tighten’’ distributions such as the one shown in Figure 2, so that the curves in the t–v plane become more compact. The underlying assumption here is, of course, that the distribution consists of narrow-band processes, i.e., is deterministic in a time–frequency sense. We discuss two approaches to this: first, the application of other distributions drawn from Cohen’s class, essentially trading time–frequency resolution for the appearance of unwanted cross terms, the second approach (which we have not yet detailed) is that of modifying the spectrogram. 3.1.2. The application of quadratic (bilinear) time–frequency methods In Figure 3 is shown a pseudo-Wigner distribution (employing a Gaussian window) calculated from the data shown in Figure 1. In this case a Gaussian window is used and the results are displayed in a fashion similar to those shown for the spectrogram. The most striking feature of this plot compared to the spectrogram is the presence of the interference terms. These are largest between the two main features, where they appear as chevrons, and inside the concave region of the second component. It is readily apparent how such features can confuse anyone trying to interpret such a distribution. It should be noted that while one generally considers the cross-terms as products of the interaction of two different signal components, a counter-example is illustrated in Figure 3. Specifically, the cross-terms within the concave region of the second component are present regardless of the first term (they are due to the curvature of the time–frequency law), the

Figure 3. The pseudo-Wigner distribution using a Gaussian window of acceleration.

438

. .   . . 

Figure 4. The cone-kernel distribution of acceleration.

(pseudo) Wigner distribution is optimal for linear time–frequency laws, and any departure generates cross-terms.† The second effect one should notice when comparing Figures 2 and 3 is that the time–frequency resolution of the latter is greatly improved; the lines defining the time–frequency laws for the two components are finer and are thus better defined. As we have noted, with Cohen’s class of time–frequency methods one aims to try to retain the good resolution of the Wigner distribution while also endeavouring to reduce the cross-terms. As one example of the behaviour of such a distribution we present the results of using a cone-kernel distribution [27]. The kernel function of these distributions has been specifically designed to retain the support properties (P4 and P5 in Table 2), and it is this requirement which results in the cone shape for the support of the kernel. Several other properties are also retained. The results from this distribution are shown in Figure 4. The results from the cone-kernel distribution have poorer resolution than the Wigner distribution but suffer less from the effects of cross-terms. This result illustrates the compromise that is inherent in this form of time–frequency method. 3.1.3. Modified time–frequency methods An alternative approach to improving the performance of time–frequency methods is to use a so-called modification of the distribution. In two papers in 1976 and 1978, Kodera et al. [7, 8] introduced a modification to the spectrogram (which they called the modified moving window method) which uses the phase information, normally discarded, to reposition the energy =S(t0 , v0 )=2 from (t0 , v0 ) to (t¯ g , v¯ i ). The new positions t¯ g and v¯ i are given by the group delay and instantaneous frequency respectively, estimates of which can be produced by the calculation of t¯ g=t0−d[arg S(t0 , v0 )]/dv0 ,

v¯ i=d[arg S(t0 , v0 )]/dt0 .

(57, 58)

Kodera et al. gave two methods for deriving these new positions. This modification was also applied to the pseudo-Wigner distribution [44]. For both distributions, modification † The phrase ‘‘cross-term’’ is now misleading, since they may arise because of an interaction of one component with itself.

– 

439

Figure 5. A modified spectrogram of acceleration.

relocates the energy to the input time–frequency law. Relocation for the spectrogram takes place in both the time and the frequency directions, whereas relocation for the PWD takes place in frequency only. This was discussed fully in reference [44]. Recently, a detailed analysis and generalization of the Kodera et al. approach has been given [45]. The modification is applied to the spectrogram of the structural response data and the results are shown in Figure 5. It should be realized that plotting these functions requires considerably more effort than for a conventional time–frequency distribution, since the data consists of points which are now scattered in a haphazard fashion over the time–frequency plane, rather than being regularly spaced. This, in part, accounts for the rather odd nature of the plot in Figure 5. However, examination of this plot reveals that there has been considerable sharpening of the plot. The two time–frequency laws are visible and are arguably finer than those generated by the Wigner distribution. To conclude this section, it is noted that wavelet transforms may be used to advantage. Reference [46] describes how wavelet analysis may be used to deduce the dispersion law for the example above. However, care (see reference [46]) is needed to ensure that appropriate post-processing of the wavelet transform is required to extract the relevant information. These details are inappropriate for inclusion here. 4. EVOLUTIONARY TIME–FREQUENCY SPECTRA

The aim of this final section is to develop the analytical basis for calculating evolutionary spectra for a class of frequency-modulated process in which the frequency modulation is due to motion of a source (in the acoustic case) or a vehicle traversing rough terrain. A series of papers and Ph.D. theses has been directed at this, and we shall summarize the essential ideas here.

440

. .   . . 

4.1.   The fundamental question that needs to be addressed is how representation (15), which says that a signal is the sum of amplitude-modulated components, may describe frequency-modulated processes. The first step is to recall an interpretation of the evolutionary spectral density as noted by Priestley [1, 2], namely, that an oscillatory process may be regarded as the response of a time-varying filter to a stationary excitation: i.e., if z(t) is a stationary process which is operated upon by a variable filter such that

g

x(t)= h(t, u)z(t−u) du,

(59)

then we may regard x(t) as an oscillatory process with representation (15), where

g

A(t, v)= h(t, u) e−jvu du.

(60)

Such interpretations of non-stationary processes in terms of time-variable ‘‘shaping filters’’ are of course widely used: e.g., in speech modelling. Our concern is how we might develop analytic evolutionary spectral density models for frequency modulated processes and the approach is as follows. We shall start with a stationary process and then distort the independent variable so that the process is ‘‘stretched out’’ or ‘‘compressed’’ so as to create a non-stationary process. For the purposes of this presentation, which deals with random processes, we begin by defining a process that is ‘‘stationary’’ when viewed as a function of some independent variable, say s (s is not, in general, time). We shall call this y(s). First, obtain a shaping filter representation for y(s) in terms of filtered white noise w(s). Now let s become a function of time, i.e., s=s(t), such that s(t) is not, in general, constant. This means that y(s) may now be regarded as a function of time, i.e., y[s(t)]=y˜ (t) say, and the ‘‘dilation’’ of the s-axis non-linearly results in a process that is frequency-modulated in form. The question is whether the s-domain constant shaping filter representation transforms to a form analogous to equation (15). We shall show that this is indeed so, but in order to achieve it we introduce the concept of covariance-equivalence, using the following method. Full details may be found in the references cited: e.g., [39, 47–52]. The process y(s) (for when s may be a space variable, for example), is assumed stationary (i.e., homogeneous) in the s-domain, with zero mean, variance sy2 , and autocorrelation function E[ y(s1 )y(s2 )]=Ryy (s2−s1 ). Now, regard s as a (deterministic) function of another variable t (time), which creates y˜ (t)=y[s(t)]. The functional dependence of s on t will be described by s˙ (t). If s˙ is constant, then y˜ (t) is a stationary process, but if s˙ is not constant y˜ (t) is obviously a non-stationary process but with the properties E[ y˜ (t)]=E[ y(s)]=constant (assumed zero),

E[ y˜ 2(t)]=E[ y 2(s)]=sy2 (constant).

Note that we shall assume s˙ q0. Even though the mean and variance are constant, the ‘‘frequency structure’’ varies in accord with our notions of frequency modulation. This is, in turn, reflected in the autocorrelation function for y˜ (t), i.e., Ry˜y˜ (t1 , t2 )=E[ y˜ (t1 ) y˜ (t2 )]=Ryy (s(t2 )−s(t1 )), which is a function of t1 and t2 and not simply (t2−t1 ) only (unless s˙ is constant).

(61)

– 

441

Figure 6. The shaping filter form for y(s).

4.1.1. Shaping filter models In order to be able to develop evolutionary spectral forms for such processes we shall now assume that y(s) can be described as the output of a shaping filter (see Figure 6) that is driven by white noise. The arguments in reference [39], etc., use state-space forms for the filters but we will omit going into detail and concentrate on conveying the basis of the ideas. It is indicated in Figure 6 that process y(s) may be formally expressed as

g

y(s)= h(s−s1 )w(s1 ) ds1 .

(62)

If s=s(t), then y˜ (t)=y[s(t)]=

g

t

h(s(t)−s(t1 ))w(s(t1 )) ds(t1 ),

(63)

−a

or y˜ (t)=

g

t

h(s(t)−s(t1 ))w(s(t1 ))s˙ (t1 ) dt1 .

(64)

−a

This shows y˜ (t) to be the output of a time variable filter driven by process w[s(t)] which is an independent variable dilated white process which must be replaced by a function of time in order to be able to proceed. The treatment of this problem may be approached formally as follows, upon noting a property of the delta function: namely, if g(t) is a function with simple zeros at t=ti , then d(g(t)) is equivalent to d(t−ti )/=g˙ (ti )=. Generalizing this slightly and applying it to the covariance of the white noise w[s(t)] results in E{w[s(t1 )]w[s(t)]}=d(t1−t)/=s˙ (t)=.

(65)

Since we have assumed s˙ q0, we can dispense with the modulus sign in equation (65). It is equations (64) and (65) that we now use. An ‘‘equivalent’’ covariance function will arise if we conceive of another white noise process, written as w1 (t)/zs˙ (t), where w1 (t) is stationary with E[w1 (t1 )w1 (t2 )]=d(t1−t2 ),

(66)

E[w1 (t1 )/zs˙ (t1 )w1 (t)/zs˙ (t)]=d(t1−t)/s˙ (t).

(67)

so that

The process w1 (t) [s˙ (t)]−1/2 is non-stationary in that it is a uniformly modulated white process, having an autocorrelation function which is indistinguishable from the required

442

. .   . . 

form in equation (65). Accordingly, we shall use w1 (t)[s˙ (t)]−1/2 in place of w[s(t)] in equation (65) and so produce a process which we shall call y1 (t), satisfying y1 (t)=

g

t

h(s(t)−s(t1 ))zs˙ (t1 )w(t1 ) dt1 .

(68)

−a

We use the notation y1 rather than y˜ since it is apparent that y1 and y˜ must differ in some respects. However, in view of the fact that equations (64) and (68) are both driven by excitations that are ‘‘covariance-equivalent’’ (namely, w[s(t)] and w1 (t) [s˙ (t)]−1/2 ), then it is reasonable to expect that y˜ (t) and y1 (t) are also covariance-equivalent: i.e., Ry˜y˜ (t1 , t2 )=Ry 1 y 1 (t1 , t2 ). That this is indeed so can easily be demonstrated. We remark it can also be shown that y˜ (t) and y˜1 (t) are ‘‘higher distribution equivalent’’, but we are only concerned with second order properties here. Furthermore, we note that the equivalence of w(s(t)) and w1 (t)[s˙ (t)]−1/2 is a manifestation of ‘‘self-similarity’’; i.e., temporal scaling reveals a similar structure. 4.2.        Evolutionary spectral forms for frequency modulated processes follow directly from the results of the previous section. The important point is that we shall use y1 (t) in place of y˜ (t) and so will use equations (59), "60) and (68). Let us formally express the stationary process w1 (t) as

g

w1 (t)= ejvt dW1 (v),

(69)

with power spectral density for w 1 (t) written Sw 1 w 1 (v)=1; then the solution of equation (68) is

g

y1 (t)= ejvtA(t, v) dW1 (v),

(70)

g

(71)

where A(t, v)=

a

h(s(t)−s(t−t))zs˙ (t−t) e−jvt dt.

0

The evolutionary spectral density for y1 (t), and hence (by covariance-equivalence) for y˜ (t) is Se (t, v)==A(t, v)=2.

(72)

4.3.      Examples of non-stationary random processes having a frequency-modulated form have appeared in the literature cited and full descriptions have been given in references [47–52]. We also note that these concepts have been used in the context of control [53]. In this section we will briefly describe three examples. 4.3.1. Vehicle motion over rough terrain [48, 51] Let us consider the motion of the mass of a vehicle accelerating over rough ground, modelled here in simplest form as in Figure 7.

– 

443

Figure 7. A simple vehicle model; the equation of motion is my=−c(y˙−h )−k(y−h ).

If the ground is modelled as h(s) having a covariance structure Rhh (j)=E[h(s)h(s+j)]=s 2 e−a=j=

(73)

then the spatial shaping filter for the ground is dh/ds+ah=sz2a w(s),

(74)

where w(s) is white. If s is a function of time, then h(s):h(s(t))=h (t). Combining this with simple linear dynamics (with v02=k/m, 2zv0=c/m) yields

&' &

y 0 1 d y˙ = −v02 −2zv0 dt h 0 0

0 (v02−2azv0 s˙ ) −as˙

'& ' & '

y 0 y˙ + 2zv0 z2a ss˙ w[s(t)]. h 1

(75)

Now using the result w1 (s(t))=w1 (t)/zs˙ and the other results given earlier in this section, yields the evolutionary spectral density shown in Figure 8. Figure 8 provides a clear indication of how the spectral density of the excitation broadens with time as the vehicle accelerates to excite the resonant behaviour of the response. Note that the frequency variable v does not run from zero. This is because at t=0 (up to which time the vehicle is at rest) the ‘‘frequencies’’ perceived by the vehicle are zero and the spectral density of the mass is concentrated at v=0 represented as a delta function, the integral of which is s 2. Extensions to the above that have been described include general velocity variations, multi-wheel vehicles and non-linear dynamics [48, 51].

Figure 8. The evolutionary spectral density of mass displacement for an accelerating vehicle.

444

. .   . . 

4.3.2. Propagating acoustic sources [47, 49, 50, 52] We shall briefly describe the formulation required for the determination of the evolutionary spectral density of the acoustic signal perceived by a fixed observer when a moving acoustic source emitting a stationary random signal passes by. The non-stationarity in this situation arises owing to range, Doppler and directivity effects. We restrict discussions to the case of a monopole moving at constant speed, and a simple ‘‘flyover’’ geometry is depicted in Figure 9. The signal received by the observer at time t is due to that generated by the source some time earlier. If q(t) is the monopole source strength, which is assumed to be a stationary random process, then it is possible to show that the (far) free field pressure at the receiver may be expressed in the form p(t)=m(t)q(t(t)),

(76)

where t=t−r(t)/c0 is the so-called retarded time (c0 is the speed of sound). The flyover geometry defines m(t) and r(t). If we assume q(t) (stationary in t, i.e., the reference frame of the source) has a shaping filter representation, we can then conceive of a process p1 (t) which is covariance-equivalent to p(t) for which we can obtain the evolutionary spectral density. For this simple geometry, i.e., a source moving straight and level over an observer, a contour plot for the theoretical evolutionary spectral density is given in Figure 10. The source is modelled as a process with a single spectral peak which is apparent from the figure as a ‘‘high’’ frequency as the source approaches and is ‘‘low’’ as it recedes. The ‘‘flyover’’ point is apparent, where the spectral density ‘‘broadens’’ when rates of change are greatest. It is noted that the Wigner distribution has also been computed and, interestingly and unusually, for this particular case the differences are visually indistinguishable. 4.3.3. Directionality patterns of moving sources Equation (76) was generalized in references [49, 50] to include explicitly directionality effects and may be written as p(t)=m(t)D[c(t)] f [t(t)].

(77)

The additional term D[c(t)] accommodates the directionality of the source and c(t) is the radiation angle relative to the observer. Lee [49] considered the problem of estimating the source directionality pattern from sound measurements which are both amplitude- and frequency-modulated. The application related to underwater sources and so reflected

Figure 9. The ‘‘flyover’’ geometry; the signal received by the observer at time t is due to the source when at position A at time t−r(t)/c0.

– 

445

Figure 10. The evolutionary spectral density/Wigner distribution for a simple flyover case.

signals (such as off the water–air interface) were also included. Approaches based upon evolutionary spectra and Wigner distributions were used with success. 5. CONCLUDING REMARKS

Time–frequency analysis has long been used as an effective tool for practitioners. This has been based on the spectrogram. The past 20 years has seen major developments in theory and understanding which now offer the practical signal analyst considerably more scope, though it is fair to say that many of the distributions described are still far from being ‘‘standard’’ approaches. However, the increasing number of applications being reported and the ready availability of new software should help ensure rapid expansion of methods that truly offer great insight into a multitude of phenomena. REFERENCES 1. M. B. P 1965 Journal of the Royal Statistical Society B27, 204–237. Evolutionary spectra and nonstationary processes. 2. M. B. P 1967 Journal of Sound and Vibration 6, 86–97. Power spectral analysis of nonstationary processes. 3. J. K. H 1968 Journal of Sound and Vibration 7, 393–416. On the response of single and multidegree of freedom systems to nonstationary excitations. 4. D. G 1946 Journal of the IEEE, London 93(III), 429–457. Theory of communication. 5. H. P. H 1970 Fourier Analysis. New York: Simon and Schuster. 6. L. Z and C. D 1963 Linear System Theory. New York: McGraw-Hill. 7. K. K, C. de V and R. G 1976 Physics of the Earth Planetary Interiors 12, 142–150. A new method for the numerical analysis of nonstationary signals. 8. K. K, R. G and C. de V 1978 IEEE Transactions on Acoustics, Speech and Signal Processing ASSP-26(1), 64–76. Analysis of time-varying signals with small BT values. 9. B. B 1987 Proceedings of the ASTED International Symposium on Signal Processing and its Applications, Brisbane, Australia, 261–268. Nonstationary signals and the Wigner–Ville distributions: a review. 10. J. B. R 1965 Journal of Sound and Vibration 2, 375–390. On the harmonic analysis of evolutionary random vibration.

446

. .   . . 

11. J. C. H and A. G. M 1990 Journal of Sound and Vibration 142, 191–202. Correlation autoregressive processes with application to helicopter noise. 12. M. B. P 1966 Journal of the Royal Statistical Society B28, 228. Design relations for nonstationary processes. 13. M. B. P and H. T 1973 Journal of the Royal Statistical Society B35, 153–188. On the analysis of bivariate non-stationary processes. 14. G. J. D 1969. Proceedings of the Symposium on Applications and Methods of Random Data Analysis, ISVR, University of Southampton. The representation of signals in frequency and time. 15. R. Y. J 1991 Ph.D. Thesis, University of Southampton. Adaptive filtering and the identification of tones in broadband noise. 16. C. W. H 1966 IEEE Transactions on Information Theory IT12, 81. An expansion of a signal in Gaussian elementary signals. 17. L. K. M and I. S. R 1967 IEEE Transactions on Information Theory IT13, 344–345. A generalisation of the Gabor–Helstrom transform. 18. J. K. H 1971 Ph.D. Thesis, University of Southampton. Frequency time methods in vibrations. 19. I. D 1991 in Advances in Spectrum Analysis and Array Processing, Vol. 1 (S. Haykin, ed.) 366–417. Englewood Cliffs: Prentice-Hall, Chapter 8: The wavelet transform: A method for time-frequency localisation. 20. L. C 1989 Proceedings of the IEEE 77(7), 941–981. Time–frequency distributions—A review. 21. T. A. C. M. C and W. F. G. M 1980 Philips Journal of Research 35. The Wigner distribution—a tool for time-frequency analysis, parts I, II, III. 22. T. A. C. M. C and W. F. G. M 1984 Proceedings of the ICASSP, 41B7. 1–4. On the time–frequency discriminations: can they look sharper than Heisenberg? 23. C. P. J and A. J. M. K 1983 Journal of the Audio Engineering Society 31, 198–223. Time–frequency distributions of loudspeakers: the application of the Wigner distribution. 24. J. S. B and A. G. P 1966 Measurement and Analysis of Random Data. New York: John Wiley. 25. W. D. M 1970 Journal of Sound and Vibration 11, 19–64. Spectral analysis of the convolution and filtering of nonstationary stochastic processes. 26. H.-I. C and W. J. W 1989 IEEE Transactions on Acoustics, Speech and Signal Processing 37, 862–871. Improved time–frequency representation of multicomponent signals using exponential kernels. 27. Y. Z, L. E. A and R. M 1990 IEEE Transactions on Acoustics, Speech and Signal Processing 38(7), 1984–1091. The use of cone-shaped kernels for generalised time–frequency representations of nonstationary signals. 28. J. J and W. J. W 1992 IEEE Transactions on Acoustics, Speech and Signal Processing 40(2), 402–412. Kernel design for reduced interference distributions. 29. L. C 1992 in Time Frequency Signal Analysis Methods and Applications (B. Boashash, editor), 3–42. Longman Chesire Ltd. Introduction: a primer on time–frequency analysis. 30. J. J and W. J. W 1992 in Time Frequency Signal Analysis Methods and Applications (B. Boashash, editor), 389–405. Melbourne: Longman Chesire Ltd. Time-varying filtering and signal synthesis. 31. F. H and W. K 1994 IEEE Transactions on Signal Processing 42(12), 3321–3334. Time–frequency projection filters and time–frequency signal expansions. 32. J. C, M. B and J. B 1988 Random Processes: Measurement, Analysis and Simulation. Fundamental Studies in Engineering, 8. Amsterdam: Elsevier. 33. P. F and O. R 1990 Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing, Albuquerque, U.S.A., 2455–2458. Wavelets and affine smoothing of the Wigner–Ville distribution. 34. O. R and P. F 1992 IEEE Transactions on Signal Processing 40(7), 1746–1757. Time-scale energy distributions: a general class extending wavelet transforms. 35. O. R and M. V 1991 IEEE Signal Processing Magazine 14–38. Wavelets and signal processing. 36. F. H and G. F. B-B 1992 IEEE Signal Processing Magazine 21–55. Linear and quadratic time–frequency signal representations. 37. J. K. H, J. S. L and R. F. H 1986 Proceedings of the ICASSP, Tokyo. The relationship between Wigner–Ville and evolutionary spectra for frequency modulated random processes.

– 

447

38. J. K. H 1992 Proceedings of the 14th ICA Beijing, L2-1. Analytic time–frequency spectra for acoustic and vibration signals. 39. J. K. H, R. F. H, Y. H, T and J. S. L 1993 in Developments in Time Series Analysis (in Honour of M. B. Priestley) (T. Subba Rao, editor) London: Chapman and Hall, 355–373. The prediction of time–frequency spectra using covariance equivalent models. 40. W. A G 1987 Statistical Spectral Analysis: A Non-probabilistic Theory. Englewood Cliffs, New Jersey: Prentice-Hall. 41. P. J. L, J. W. P and L. E. A 1993 IEEE Transactions on Signal Processing 41(2), 750–767. Bilinear time-frequency representation: new insights and properties. 42. T. J. W and J. S. B 1993 Journal of Sound and Vibration 163, 101–122. The application of the Wigner distribution to the identification of structure-borne noise components. 43. N. K and J. K. H 1993 Journal of Mechanical Systems and Signal Processing 7(5), 425–449. Application of cepstral techniques for the determination of reflection coefficients for dispersive systems: part I theory and numerical results; part II, comparison between theory and experiment. 44. J. C. M and J. K. H 1994 Mechanical Systems and Signal Processing 8(3), 243–258. A comparison between the modified spectrogram and the pseudo-Wigner–Ville distribution with and without modification. 45. F. A and P. F 1995 IEEE Transactions on Signal Processing 43(5), 1068–1089. Improving the readability of time–frequency and time–scale representations by the reassignment method. 46. Ph. G and P. R. W 1995 IEEE Signal Processing Chapter U.K. Symposium on Applications of Time–Frequency and Time–Scale Methods, University of Warwick, Coventry. Wavelet transforms for the analysis of dispersive systems. 47. Y. H. T 1983 Ph.D. Thesis, University of Southampton. Aspects of evolutionary spectral analysis with applications to problems in acoustics. 48. R. F. H 1983 Ph.D. Thesis, University of Southampton. The nonstationary response of vehicles on rough ground. 49. J. S. L 1989 Ph.D. Thesis, University of Southampton. Time-varying filter modelling and time–frequency characterisation of nonstationary sound fields due to a moving source. 50. J. S. L and J. K. H 1987 Proceedings of the ICASSP, Dallas, 41.8.1–4.1.8.4. Estimation of the directionality pattern of a moving acoustic source. 51. R. F. H and J. K. H 1986 Journal of Sound and Vibration 107, 29–38. Evolutionary (frequency/time) spectral analysis of the response of vehicles moving on rough ground by using covariance equivalent modelling. 52. J. K. H and R. F. H 1985 Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing, Tampa, Florida, 1025–1027. Wigner–Ville and evolutionary spectra for covariance equivalent non-stationary random processes. 53. S. N and G. V. R 1992 Vehicle System Dynamics 21(2), Active control of nonstationary response with nonlinear suspensions.