See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/277165916

Algorithms for Time-Frequency Signal Analysis

Article · January 1992

CITATIONS READS

109 744

2 authors, including:

Boualem Boashash
Qatar University
526 PUBLICATIONS 18,358 CITATIONS

SEE PROFILE

All content following this page was uploaded by Boualem Boashash on 20 May 2016.

The user has requested enhancement of the downloaded file.

Published in Australia, New Zealand, the Pacific, Asia and Africa by
Longman Cheshire Pty Limited
Loilgman House
Kings Gardens
95 Coventry Street
Melbourne 3205 Australia
ISBN 0 582 71286 6
TIME • FREQUENCY SIGNAL ANALYSIS
METHODS AND APPLICATIONS
Offices in Sydney, Brisbane, Adelaide and Perth. Associated companies,
branches and representatives throughout the world.

Copublished in the Western Hemisphere, United Kingdom

and Europe by Halsted Press: an Imprint of
John Wiley & Sons, Inc.
New York Toronto Chichester

Copyright @ Longman Cheshire 1992

First published 1992 Edited by Boualem Boashash
All rights res~ed. Except under the conditions described in the Copyright Act
1968 of Australia and subsequent amendments, no part of this pUblication may
be reproduced, stored in a retrieval system or transmitted in any fonn or by any
means, electronic, mechanical, photocopying, recording or otherwise, without
the prior pennission of the copyright owner.

Designed by Nadia Graziotto

Printed in Hong Kong

National Library of Australia

Catalogulng-in-Publlcatlon data

Time-frequency signal analysis.

Bibliography.
Includ6 index.
ISBN 0-582-71286-6
l.Signal processing. 1. Boashash, Boualem.
621. 38223

Library or Congress Cataloguing-in-Publication data

Time-frequency signal analysis-methods and applications/edited by Boualem

Boashash.
p. en.
Includes bibliographical references and index.
ISBN 0-470-21821-5 ~.~
1. Signal processing. I. Boashash, Boualem. •••
•••
YK5102.T585 1992
621.382'2-dc20 91-32658
Longman Cheshire
CIP
6M WI LEY ~:tisTED
162 Time-Frequency Signal Analysis

[53] Wahba, G., "Practical Approximate Solutions to Linear Operator Chapter 7

Equations when the Data are Noisy", SIAM J. Numer. Anal.,
Volume 14, Number 4, September 1977, pp. 651-667.
[54] Wahba, G., "Smoothing and Ill-Posed Problems", in M. A.
Golberg, editor, Solution Methods for Integral Equations, Plenum
Press, 1978, pp. 183-194. Algorithms for Time-Frequency
[55] West, G., "Eine schnelle Mellin-Transformation", Computing, Signal Analysis
Volume 33, pp. 237-245, 1984. Boualem Boashash and Andrew Reilly
[56] Whitehouse, H. J. and Speiser, J. M., "Algorithms and Architec-
tures for Array Signal Processing", Presented at the Indo- U. S.
Workshop on Systems and Signal Processing, Bangalore, India,
January 8-12, 1988. The conference was sponsored by the U.S.
Office of Naval Research, the Indian Department of Science and
Technology, and the IEEE Bangalore chapter. A book of the ex-
tended abstracts of the conference papers is available, and a con-
ference proceedings is planned.
[57] Whitehouse, H. J. Boashash, B. and Speiser. J. M., "High
Resolution Processing Techniques for Temporal and Spatial Keywords: fast algorithms, efficient implementation, Fortran sub-
Signals" , Presented at the Workshop on High Resolution Methods routine, optimizations, analytic signal, Cohen's class, Wigner-Ville
for Underwater Acoustics, organized by GRETSI, Juan les Pins, distribution
France, June 16, 1989. To appear as a chapter in a volume edited
by M. Bouvet and G. Bienvenu, in the Springer-Verlag Lecture
Notes in Comp'uter Science.
[58] Wiener, Norbert, The Four-ier Integral and Cer-tain of its Appli- 1 Introduction
cations, Dover Publications Inc., New York, p.3.
This chapter presents algorithms which implement time-frequency sig-
[59] Wiener, Norbert, Generalized Har-monic Analysis and Tauber-ian nal analysis techniques on computer systems. Fortran code fragments
Theorems, The M.LT. Press, Massachusetts Institute of Technol- are included.
ogy, Cambridge, Massachusetts, 1964. A generalized framework for time....Jrequency distribution calcu-
[60] Wilcox, C. H. "The Synthesis Problem for Radar Ambiguity lation is provided, and a number of speed-up optimizations are de-
Functions", Mathematics Research Center, MRC Technical Sum- scribed. Algorithms are presented for the calculation of all of the
mary Report No.157, The University of Wisconsin, April 1960. time-frequency representations listed in table 7.1.
[61] Zwick, Philip E. and Imre Kiss, "A New Implementation of the Most popular time-frequency representations can be expressed in
Mellin Transform and its Application to Radar Classification terms of the general bilinear time-frequency distribution representa-
of Ships" IEEE Transaction on Pattem Analysis and Machine tion proposed byL.Cohen [6] (see also chapter 1 in this book). This
Intelligence, Volume PAMI-5, Number 2, pages 191-199, 1983. is shown in equation 1 below: *

Pz(t, f) =/!! e j27rv (u-t) g( v, r )z( u + ~ )z*( u - ~ )e- j27r!r dvdudr (1)

Here z is an analytic signal, (see [11] i.e., a complex signal which

contains no negative frequencies. The function g( v, r) determines
the characteristics of the time-frequency distribution. For example
if .q( v, r) = 1 then the distribution formed is the Wigner-Ville
 Algorithms for Time-Frequency Signal Analysis 165
164 Time-F'requency Signal Analysis

distribution. If the integration with respect to v is performed, then where (~) denotes the discrete time convolution, and fm-+k denotes the
the equation becomes: discrete Fourier transform from time (n) to frequency (k).

Pz(t, f) = JJG(t, r)z(u + ~)z*(u - ~)e-j2'71Ir dudr (2)

2 Analytic signal calculation
This distribution (pz) is related by Fourier transforms to the tin,te- In most practical cases, the signals to be analysed consist only of real
lag representation Rz (t, r) , the Doppler-frequency representatl~n values. In these cases it is necessary, as a first step towards signal
Tz(V,f), and the Doppler--delay representation, Az(v,r) as shown m analysis, to form the corresponding analytic signal, as explained in
figure 7.1. [3].
The direct method of producing the analytic signal is to use its
7Pjl,fl~ definition: a signal with no negative frequency components. That is,
form the Fourier transform of the signal, set the negative frequency
values to zero, and perform the inverse Fourier transform. Unfortu-
R'(I~ ~,fJ nately, finite data length effects (windowing) cause this method to
produce undesirable ripples in the signal. Setting the negative fre-
quency components to zero has the same effect in the time domain as
AAv, '1') any naive frequency domain filtering operation: the resulting response
has sizable ripples for any non-harmonic frequencies.
Another method is to use a Hilbert transform filter to produce
Figure 7.1. Relationships between continuous time representations.
the required complex component of the signal, which when added to
The equivalent discrete time relationships, related by discrete the real part produces the desired analytic signal. This is shown in
Fourier transforms are Rz(n, m), pAn, k), Tz(l, k), AAl, m), as shown equation 5.
in figure 7.2. z(n) = x(n) + jH[x(n)] (5)

)/P'("kJ~ where HO is the Hilbert transform [11], and has the ideal impulse
response shown in equation 6.

R.(n, m) rAI, k)
2 sin 2 (7rn/2)

~A
( )
h n = { 7rn
for n =I 0 (6)
. 0, for n = 0
AAI,m)

Figure 7.2. Relationships between discrete time representations. The Hilbert transform can be implemented as a finite-impulse
response (FIR) digital filter, which may be more efficient than the
The expression to be used for discrete time implementation of these Fourier technique., depending on the length of the signal to be
transforms is the discrete time equivalent of equation 2, and is shown transformed, and the length of the filter impulse response used.
in equation 3. Although an ideal Hilbert transform filter has an infinite impulse
M M
response length, in practise an FIR filter length of 79 samples has
pz(n, k) = L L G(p, m)z(n+p+m)z*(n+p-m)e- j41l'mk/N (3) been shown to provide an adequate approximation[4].
This filter could be a rectangular-windowed version of the infinite
m=-Mp=-M
length filter, or it could be an optimal implementation calculated
This can also be expressed as: using one of the filter design algorithms, such as the Remez-exchange
algorithm. This produces a filter which is optimal in a minimax or
(4) Chebyschev sense [I1J.
166 Time-Frequency Signal Analysis Algorithms for TimecFrequency Signal Analysis 167

3.2 Convolution in the n (time) direction

3 General approach to computation of TFDs
This is simply the matrix multiplication of the weighting function
The discrete-time definition of Cohen's class of time-frequency distri- values G(n,m) with the bilinear kernel values Kz(n,m). This must
butions given in equation 3 forms the basis of the general approach to be performed for each time instant for which output is required.
implementation of time-frequency distributions. This approach can In an actual implementation, either the kernel matrix or the
be expanded into three steps: selection function matrix may be calculated at the point of use (in
1. Form the bilinear product Kz(n, m) = zen + m)z*(n - m). This the convolution). This can save memory space and sometimes time.
is known as the 'bilinear kernel' or just 'kernel' where meaning If the kernel values are only going to be used once, then a considerable
is obvious. memory saving can be made by not storing them at all (remember,
they are complex, so would require a very large array to store).
2. Convolve the kernel with the desired determining function' Conversely, if the selection function is trivial then this need not be
G(n, m) in the n (time) dimension. stored either.
3. Calculate the discrete Fourier transform of this result, to Since the selection function G(n, m) is 8(n) for the Wigner-Ville
produce the time slice of the desired distribution. distribution, the convolution step does not really occur at all.
Most common bilinear transforms exhibit symmetry in both the
These steps will now be elaborated on, and issues arising from m and n directions, and this fact can be used to reduce both
implementation on computing systems mentioned. The next section, computation and storage requirements. In the code fragment shown
section 4, presents an example of some of the optimizations that can in section 3.4, this symmetry is used by multiplying the same G(n, m)
be applied to these techniques for a specific example of Cohen'sdass, value with the required two bilinear kernel values at the same
the Wigner-Ville distribution. time. The delay axis symmetry of G and the corresponding kernel
Hermitian symmetry allow the uppet half of the result matrix R to be
calculated from the lower half through conjugation, rather than direct
3.1 Calculation of the bilinear kernel calculation. Care must be taken that symmetry exists around the R( 1)
element, rather than including it. That is, R(FFTLEN)=conjg(R(2»,
It can be easily shown that the bilinear kernel has Hermitian not conjg(R(1», where FFTLEN is the data length of the Fourier
symmetry: transform to be used.

zen + m)z*(n - m) for m 2: 0

(7)
Kz(n, m) = { K;(n, -m) for m < 0 3.3 Discrete Fourier transform
This means that values of the kernel need only be calculated for The discrete Fourier transform of each time-slice of the distribution
positive time lags. is calculated using a fast-Fourier transform routine (FFT). Examples
It is normally not necessary to store the calculated kernel values of this can be found in [11 J. The most important point to note is that
in an array, as each will only be used once. (Of course there are the use of one of the fast algorithms puts restrictions on the length
exceptions to this rule, such as when a sliding analysis window is of each data' segment (usually to a power of two). This may require
used, with a large overlap between successive windows.) Since the zero-padding of the filtered kernel before performing the transform.
kernel array would be rectangular, twice as long (in the n direction) Another point to note is that these routines operate only on positive
as wide (the m direction), storage requirements for the kernel would time and frequency values, rather than the positive and negative values
increase with the square of the problem size. In a finite storage system expected by simple translation from the continuous domain. That is,
this will reduce the useful analysis window length. rather than representing frequencies from -~fs to ~is, with the zero
The alternative is to calculate the kernel values as they are frequency value appearing in the middle of the array, frequencies from
required, and this is the way the code in section 3.4 works. This o to is are used, with the zero frequency value appearing as the first
code includes the concurrent convolution with the kernel specification element in the array. This has absolutely no effect on the algorithms or
function, G(n, m), which will be described in section 3.2 below, and calculations, due to the cyclic nature of the discrete Fourier transform,
so is best studied after that section has been read. but does affect the way values are referred to. The 'nega.tive' frequency
values are stored at the 'top' of the arra.y, and are usually referenced
168 Time-Frequency Signal Analysis Algorithms for Time-Frequency Signal Analysis 169

with indices like: FFTLEN+1-i. real G(KERNMAX,KERNMAX)

The FFT routine referred to in the code fragments in this chapter integer FFTLEN,mf,hlf,KERNMAX,m,p
performs calculations in place, and requires both length and order
parameters:
do 1 m=O,hlf
subroutine FFT(A,mf,FFTLEN) R(m+1)=G(1,m+1)*z(m)*conjg(z(-m»
c A is an array of complex do 2 p=l,hlf
2 R(m+1)=R(m+1) +G(p+1,m+l)* (
c mf is the calculation order, FFTLEN=2**mf
1 z(p+m)*conjg(z (p-m) + z( -p+m) *conjg(z( --p-m»
c FFTLEN is the length of the array, a power of 2
2 )
if (m. gt. 0) R(FFTLEN-m+1) =conjg (R(m+l) )
The code to find the next power of two greater than or equal to 1 continue
the length of a signal is a simple loop:
c make sure any points not containing calculated values aTe set
c calculate FFTLEN, the minimum (2**mj) .ge. lwin c to zero, as the FFT length may be gTeater than the window:
mf=O
FFTLEN=l do 3 m=hlf+2,FFTLEN-hlf
6 if (FFTLEN.ge.lwin) goto 7 3 R(m)=(O.O,O.O)
mf=mf+1
FFTLEN=FFTLEN+FFTLEN call FFT(R,mf,FFTLEN)
goto 6 return
7 continue end

3.4 Arbitrary transform code fragment 3.5 Notes on implementation techniques

This routine performs the convolution of the bilinear kernel with the The code presented to perform the convolution with the kernel in
distribution specification matrix G, for a signal window of length section 3.4 used negative indices into the signal array z, which ~ay
2 x hlf + 1, centered on element 0 (see below). It then takes the Fourier seem strange, and may not work with a compiler that generates code
transform of the result, producing one time-slice of the distribution. to check array bounds, but is more efficient than the alternatives.
G is known to be symmetrical in both time and lag dimensions, The same effect ca~ be achieved (a little less efficiently) by passing
and so values are only used from the positive quadrant. the whole of the sIgnal (z) array, and also an index i to indicate
It is also known that the resulting spectrum is real, so the result the center of the analysis window, which must then be added into all
of the convolution must have Hermitian symmetry. Thus half of the references to z. That is, the first line of code inside the outer loop in
RO values are computed directly, the others from symmetry. (R is the section 3.4 would change from:
array used to store the kernel products for different lags.) R(m+1)=G(1,m+1)*z(m)*conjg(z(-m»
Finally, it is known that the z array origin is not situated at the to:
beginning of the signal array, but is at least 2 x hlf from either end.
Thus negative array indices may be used, and a substantial complexity R(m+l)=G(1,m+l)*z(i+m)*conjg(z(i-m»
reduction may be realized, as no special cases are introduced at the
beginning and end of the signal (see section 3.5). . Both. implementation techniques have a problem when the analysis
wmdow IS not wholly contained within the extent of the signal, as in
subroutine transform(z,G,FFTLEN,mf,hlf,KERNMAX,R) that case some of the values used in calculations must be zero. There
complex z(1), RCi) are a number of ways to solve this problem:
170 Time-Frequency Signal Analysis Algorithms for Time-Frequency Signal Analysis 171

(I The obvious solution is to check each array index before 9 write(2,*) real(choiCi»
accessing the value, and return zero if the access is outside the 8 continue
array bounds. This is inefficient, since it requires a lot of extra
work in the inner loops of the convolution.
(I Another technique is to move the range checking outside the
loops, so that the indices never exceed the array bounds. This Table 7.1 Some TFDs and their determining functions G ( n, m)
can be quite efficient, especially where the analysis window is
very long. In that case the time required to do the checking, and
the reduced symmetry available may be offset by the reduction
in the number of computations required at the ends of the signal. Time-Frequency Representation G(n,m)

(I The technique used here is to pad the signal in the signal array
with sufficient zero values that the array will not be indexed
ben) mE [-(M-l) (M-l)]
out of bounds. Thus the correct windowing effect is achieved, Windowed Discrete WVD ., 2 ' 2

without any bounds checking within the convolution loops. Here 0 otherwise
the reduced overhead in the 'working' loops makes up for the
extra 'dead' iterations involving computation with zero values. Smoothed WVD using a rectangular window 1 n E [-(P-l) (P-l)]
p 2 ' 2
of odd length P 0 otherwise
This approach is shown in the code fragment below, extracted
frorp. a Choi-Williams distribution program. It reads the data into the Rihaczek-Margenau ~[b(n + m) + b(n - m)]
signal array, converts the signal to an analytic signal (with a call to
SIGANA), and then produces nplts spectra using an analysis window STFT using a Rectangular Window of odd 1 1m + nl ~ (P;-l)
P
of length lwin centered points with a time separation of res samples. length P, 0 otherwise

c make hif max (2*hlf+l .le. lwin) 1

hlf=(lwin+l)/2-1 Born-J ordan-Cohen Iml+l Iml ~Inl
0 otherwise

c read the signal

Choi-Williams (parameter IT) ~ e-un2j4m2
c this is where the zero-padding buffers are built at either 2m

c end ...
do 10 i=l,LWINMAX-l
z(i)=(O.O,O.O)
10 z(n+LWINMAX+i)=(O.O,O.O)
3.6 Code fragments to generate G(n,m)
do 5 i=LWINMAX,LWINMAX+n-l
read (1, *) x Thi~ sectiOI~.shows example code fragments to generate the distri-
5 z(i)=cuaplx(x,O.O) butIOn specIfication matrix G( n, m) for the functions specified in ta-
ble 7.l.
form the analytic signal: (and get nl as power of 2 .ge. n) . Note that for many of these distributions, this computation tech-
c
call SIGANA(n,z(LWINMAX), .false. ,nl) mque ?f con,"olving the bilinear kernel with the G( n, m) array is obvi-
ously mefficlent, because most of the values of G( n, m) are zero, and
th~ others are a constant (e.g., Wigner-Ville distribution smoothed
do 8 ii=O,nplts-l ~lg~er-yille distri?utio~, Rihaczek-Margenau, and STFT). For these
t=ii*res+hlf+LWINMAX dlstnb~tIOns, algorlthm Improvements yield much higher performance.
call transform(z(t),G,FFTLEN,mf,hlf,KERNMAX,choi) The Wlgner-Ville distribution is the only one for which such an im-
do 9 i=i, FFTLEN provement will be demonstrated.
172 Time- Frequency Signal Analysis Algorithms for Time-Frequency Signal Analysis 173

G(n,m) for Choi-Williams distribution wt=1.0/(real(i)+1.0)

do 2 j=O,i
c calculate the distribution specification matrix "0" 2 G(i+1,j+1)= wt
do 2 i=O,hlf do 3 j=i+1,hlf
2 G(i+1,1)=0 3 G(i+1,j+1)= 0.0
G(1,1)=1.0 4 continue

do 4 j=1,hlf G( n,m) for STFT

wt=O.O
do 3 i=O,hlf c calculate the distribution specification matrix "0"
G(i+1,j+1)= exp(-(sigma*i*i)/(4*j*j» wt=1.0/(real(2*hlf-l»
3 wt= wt+2*G(i+1,j+1) do 4 i=O,hlf
wt=wt-G(l, j+1) do 2 j=O ,hlf-i
c normalize array so that we know 2.:n G(n,j) = 1.0 2 G(i+l,j+l)= wt
do 4 i=O,hlf do 3 j=hlf-i+1,hlf
4 G(i+1,j+1)=G(i+1,j+1)/wt 3 G(i+1,j+1)= 0.0
4 continue

In this routine, the variable 'wt' is used to ensure that the sum G( n,m) for Rihaczek-Margenau distribution
along every constant 'j' (parallel to the time axis) is exactly one, a
condition necessary to preserve the marginals of the distribution, as c calculate the distribution specification matrix "0"
described in [2], Property 2-7. These are repeated here: do 2 i=O,hlf
do 2 j=O ,hlf
2 G(i+l,j+1)= 0.0
n
do 3 i=l,hlf
L:pz(n, k) = Iz(n)12 3 G(i+1,i+1)= 0.5
k G(l, 1)=1. 0
This is obtained if

G(n,O) = 8(n) and L: G(n, m) = 1 G(n,m) for smoothed WVD

n
Note: p is half the length of the smoothing region, calculated in a
This condition would not necessarily be met otherwise, due to finite similar manner to the calculation of hlf from lwin.
precision in the calculations. It should be noted that the Choi-
Williams distribution selection function: c calculate the distribution specification matrix "0"
wt=1.0/real(2*p-l)
G( n,m ) - JrJ/7f
- - - e _(y-n2/4m
2 do 4 i=O,hlf
2m
J I do2 j=O ,p-l
I
does not meet these criteria, and does not hold for m = 0, at which 2 G(i+l,j+l)= wt
point the rule G(n, 0) = 8(n) must be imposed. do 3 j=p,hlf
3 G(i+l,j+1)= 0.0
G(n,m) for Born-Jordan-Cohen distribution
4 continue
c calc'ulate the distribution specification matrix "0"
do 4 i=O,hlf
174 Time-Frequency Signal Analysis Algorithms for Time-Frequency Signal Analysis 175

3.7 Effect of windowing 4. Where the distribution will be calculated for every time value
there is another optimization, presented by Eilouti and Khadr~
Time-windowed distributions will often be calculated in preference to [~l. Here a recursive technique is used to calculate the analytic
full-length distributions. Two reasons for this. practice are: sIgnal for successive windows. This will not be described here .
• To reduce unnecessary computation, when it is known that the 5. The formulation of the discrete Wigner-Ville distribution (and
signals of interest have a limited time-extent within the data all of the other distributions in Cohen's class) presented here
collected; or produces a result with a frequency scaling of f = M where
·
th t echmques, . 2M'
• When the data stream is continuous (real-time analysis). For oer such as the short-tIme Fourier transform result
this case, it is impractical to collect all of the data before in a scaling of f = ';;. The consequence of this is that the
analysis, as results are required on early sections before later resulting transforms are twice as long as one might expect. To
sections have been collected. overcome this Sun, Li, Sekhar and Sclabassi [12J have shown
how to pr~duce a transform with the usual length, and a saving
Instead of performing all of the operations on the entire signal, only of a~proxlInately half the computations. Since this technique
sections of the signal are read into the data arrays at once, displacing reqUIres the replacement of the fast Fourier transform with their
previous values, and the calculations are repeated for each 'window'. own 'fast Fourier transform in part' (FFTP) it will not be
If the windows are overlapped, then some data must be retained, and described further here. '
shifted forward within the data arrays by an amount equal to the
time-displacement of the window. Only the code for the outer loop and the transform routine will
be presented. It is assumed that the signal array 'z' and input and
out~ut files have been attended to by the surrounding code in a similar
fashIOn to that presented in section 3.5.
4 Wigner-Ville distribution implementation nplts=n/res
c
Although covered in conjunction with the other bilinear time-frequency
c since we're using the 'do two transforms at once'
distributions in the previous section, the implementation of the Wigner-
Ville distribution will be presented again here, as an example of how c optimization, nplts must be even, make it so (possibly
to optimize the general calculation for special cases. c at the expense of the last plot:
There are a number of optimizations applicable to the calculation c
of the Wigner--Ville distribution (some of these also apply to the nplts=(nplts/2)*2
calculation of other TFDs):
1. The array G( n, m) is trivial, and need not be stored or calcu- c Here is the outer loop of the actual distribution calculation:
lated. Instead, its effect can be incorporated into the algorithm. c Call the distribution routine for each two windows of the
2. The bilinear kernel has Hermitian symmetry, and so only values c signal. The signal array z() is passed with a time ojJ.set, t,
for positive lags need be calculated (This applies for all kernels.) c so that the transform2() routine can operate as though t was always
c zero.
3. Two result spectra can be calculated by each Fourier transform.
This is possible because the spectra are known to be real, the
do 8 ii=O,nplts/2-1
Fourier transform is known to be a linear operator, and the
fast Fourier transform routine being used operates exclusively t=ii*res*2+hlf+LWINMAX
on complex data. call transform2 (z(t) ,z(t+res) ,FFTLEN,mf,hlf,wvd)
do 9 i=l,FFTLEN
The best way to use this optimization is to multiply the second
set of lags by j, so that the resulting spectrum-slice will appear 9 write (2 ,*) real(wvdCi»
in the imaginary part of the Fourier transform result. That is: do 11 i=l,FFTLEN
11 write(2,*) aimag(wvd(i)
8 continue
 Algorithms for Time-Frequency Signal Analysis 177
176 Time-Frequency Signal Analysis
R(m+i)=vi+cnnplx(-ainnag(v2),real(v2))
R(FFTLEN-m+i)=c°nUg(vi)+cnnplx(ainnag(v2),real(v2))
close (2) 1 continue
999 return c make sure any points not containing calculated values are set
end c to zero, as the FFT length may be greater than the window:
c -------------------
do 3 m=hlf+2,FFTLEN-hlf
subroutine transform2(zl,z2,FFTLEN,mf,hlf,R)
3 R(m)=(O.O,O.O)
connplex zi(i), z2(1), R(i)
connplex vi, v2 call FFT(R,mf ,FFTLEN)
integer FFTLEN,mf,hlf,m
return
end
c This routine forms the bilinear lag array for two kernel
c time-slices at once, for lag values up to hlf. (i. e., using
c a window of length 2*hlf+l centered on elements 0 and res.)
c
c The second set of lag values (those around res) are added
c into the lag array in an 'odd 'fashion, as opposed to the
'even' fashion that the first set were inserted.
5 Other· methods
c
c While the general approach presented will allow calculation of most
c Finally, the Fourier transform of the lag array is taken, time-frequency distributions, it does not cope with all possibilities
c forming the two time-slices at once, in the real and imaginary and. fo.r s~me popular cases it is not the most efficient approach:
parts of the array. OptimIzatIOns of the method for the Wigner-Ville distribution have
c
already been noted. If work is to center on one of these other
c
tech~iques, then more computationally efficient approaches may be
c This is achieved by multiplying the second array of lags by
c~nsidered. Some techniques not covered by the general technique
c j before performing the Fourier transform. Since this is a
wIll ~e commented on 15riefly in this section, and references will be
c linear operator, this also rotates the resulting spectrum into provIded for source~ of further information.
c the imaginary plane.
c
c Finally, it is known that the z array origin is not situated 5.1 The wavelet transform
c at the beginning of the signal array, but is at least 2*hlf
This is a linear transform, as opposed to the bilinear distributions
c from either end. Thus negative array indices may be used, and discuss~d in t~is chapter, and has the interesting property of time
c a substantial complexity reduction may be realized, as no resolutIOn whIch is a function of frequency. See [7) for more
c special cases are introduced at the beginning and end of the information.
c signal (see section 3.5)

vi=zl(O)*conUg(zl(O)) 5.2 Cross-Wigner-Ville distribution

v2=z2(O)*conjg(z2(O)) ~he most significant differences between the implementation tech-
R(1)=vl+cnnplx(-ainnag(v2),real(v2)) mques for the Cross-Wigner-Ville distribution and those already de-
do 1 m=i,hlf scribed are:
vl=zi(m)*conjg(zi(-m))
v2=z2(m)*conjg(z2(-m)) • the necessity to read and store two data streams·,
 Algorithms for Time-Frequency Signal Analysis 179
178 Time-Frequency Signal Analysis

• the necessity to provide storage for complex numbers m the

5.5 The Q-distribution
output matrix; A time-frequency technique that will not be covered in depth here is
the Q-distribution of Altes [1], so called because it shares the constant-
• the loss of some of the symmetry optimizations that could be Q or proportional bandwidth property of the Wavelet transform
performed with the bilinear transforms. although it is based on the bilinear Wigner-Ville distribution rathe;
than a linear transform. The paper [1] describes how this ~aybe
Once these are taken into account, all of the presented techniques done efficiently, after which the calculation would proceed as for the
and program structure still apply. The cross Wigner-Ville distribution Wigner-Ville.
is discussed in the chapter by Boles in this volume [5].

5.3 Short-time Fourier transform 6 Description of the TFSA package

The short-time Fourier transform (STFT) is a special case because it
can be efficiently calculated by taking the spectrum of windows of the TFSA is th~'c~frent name of a package of signal analysis tools built
data. The program structure will be somewhat similar, from the point up over a period of time by B. Boashash and his students. Its primary
of view of data input, output and windowing. Instead of calculation function is to produce time-frequency representations and plots from
of the bilinear kernel, and convolution with the relevant determining time series, although it has a number of auxiliary functions. Currently
function, all that is necessary is to form the Fourier transform of each the spectrum analysis tools include the short-time Fourier transform
window of data, and output the squared magnitude of each frequency the Wigner-Ville distribution, an auto-regressive (parametric) model
point. Some care is necessary, if windowing is required. Unlike the based spectrum estimator, Wigner-Ville distribution modified to
bilinear transforms noted, this technique will produce both positive use an auto-regressive model spectrum estimator, Choi-Williams
and negative frequency values (which will be zero for analytic signals), distribution, Born-Jordan-Cohen distribution and ZAM distribution
and so only half of the calculated frequency vector need be output. [14].
Since the Fourier transform approach assumes signal stationarity It also has a suite of test signal generation routines with which
over the period of the window, care must be taken to optimize the a wide variety of test signals, including Gaussian white noise can be
window length for signals that are not stationary. Where linear produced and manipulated.
frequency modulated signals are involved, the optimum window length . The package consists of a large number of stand alone programs
lmked by common file formats and an interactive menu/form based
is (~; )-~. Wi.ndows shorter than this cause unnecessary spreading of front end, which simplifies and speeds the use of these programs
the signal due to the width of the window in the frequency domain. dramatically. The code fragments presented in this chapter were
Longer windows cause a smearing of the signal within the window due extracted from this package, which is constantly being updated. It
to the non-stationarity of the frequency law contained within it. is available for a small cost from the authors, and may be ordered
with the order form which appears at the back of the book.
5.4 Parametric methods
All of the parametric spectrum estimation techniques can be used for
time-frequency spectrum analysis in the same manner that the Fourier 7 Summary
transform is used: by the assumption of short time stationarity. For a
very detailed discussion of these methods, consult the books by Kay In this chapter we have presented a general approach to the calculation
and Marple [9,10]. of the time-frequency energy distributions based on Cohen's class of
The higher resolution offered by these techniques can also be distributions. The principle steps of this approach are:
applied to the bilinear distributions, by replacing the final Fourier
transform with a parametric spectrum calculation. This generally 1. Produce the analytic signal from the real data sequence to be
requires a much greater amount of computation. An algorithm is analysed (section 1).
described by Whitehouse, Boashash and Speiser in [13].
2. Calculate the bilinear kernel values (section
 180 Time- Frequency Signal Analysis Algorithms for Time-Frequency Signal Analysis 181

3. Convolve the bilinear kernel with the distribution specification A rray Processing (Prentice-Hall Signal Processing Series), S.
function G(n, m), (sections 3.2, 3.6). Haykin, editor, Englewood Cliffs, NJ, Prentice-Hall, 1990.
4. Take the discrete Fourier transform with respect to lag (m) for [8J H. H. Eilouti and L. M. Khadra, "Optimized Implementation of
each time instant (section 3.3). Real-Time Discrete Wigner Distribution," Electronics Letters,
25, no. 11, pp. 706,707, 25 May 1989.
This process was demonstrated for the case of the Choi-Williams [9] S. Kay, Modern Spectral Analysis, Englewood Cliffs, NJ, Prentice-
distribution in the code fragments in sections 3.4 and 3.5. Hall, 1989.
This algorithm can be simplified and optimized quite significantly
[10J S. L. Marple, Jr, Digital Spectral Analysis with Applications,
for a number of common members of Cohen's class, due to extra
Englewood Cliffs, NJ, Prentice-Hall, 1987.
properties of the kernel. The Wigner-Ville distribution was presented
as an example of how the kernel selection function G(n, m) can be [11] A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal
incorporated into the algorithm, rather than being stored or calculated Processing, Englewood Cliffs, NJ, Prentice-Hall, 1989.
(see section 4). [12] M. Sun, C. C., Li, L. N. Sekhar and R. J. Sclabassi, "Effi-
References were presentedfdr a number of other popular time- cient Computation of the Disc,r~tet'p'seudo- Wigner Distribution,"
frequency signal analysis techniques which do not fit into this frame- IEEE Trans. Acoust. Speech Signal. Process., 37, no. 11, pp.
work, in section 5. 1735-1742, November 1989.
The final section (6) briefly described a signal analysis package [13] H. J. Whitehouse, B. Boashash and J. M. Speiser, "High Res-
written by the authors and their colleagues which utilizes most of the olution Processing Techniques for Temporal and Spatial Sig-
techniques presented. nals," in High Resolution Techniques in Underwater Acoustics
(L~cture Notes in Control and Information Series), New York-
HeIdelberg-Berlin, Springer-Verlag, 1989.
[14J Y. Zhao, L. E. Atlas and R. J. Marks "The Use of Cone-Shaped
References Kernels for Generalized Time-Frequency Representations of
Nonstationary Signals," IEEE Trans. Acoust. Speech Signal.
[1] R. A. Altes, "Wideband, Proportional Bandwidth Wigner-Ville Process., 38, no. 7, July 1990.
Analysis," IEEE Trans. Acoust. Speech Signal. Process., 38,
June 1990.
[2J B. Boashash, "Time-Frequency Signal Analysis", in Advances
in Spectral Analysis and Array Processing (Prentice-Hall Signal
Processing Series), S. Haykin, editor, Englewood Cliffs, NJ,
Prentice-Hall, 1990.
[3J B. Boashash, "Note on the Use of the Wigner Distribution of
Time-Frequency Signal Analysis", IEEE Trans. Acoust. Speech
Signal. Process., 36, pp. 1518-1521, September 1988.
[4] B. Boashash and P. Black, "An Efficient Real Time Implemen-
tation of the Wigner-Ville Distribution", IEEE Trans. Acoust.
Speech Signal. Process. ASSP 35, no. 11, pp. 1611-1618,
November 1987, CRISSP Reference CSP 87/6.
[5J P. Boles, "Application of the Cross-Wigner-Ville Distribution
to Seismic Surveying" in TFSA Methods and Applications, B.
Boashash, editor, Longman Cheshire, 1991.
[6] L. Cohen, "Time-Frequency Distributions-A Review," IEEE
Proc., 77, no. 7, pp. 941-981, July 1989.
[7J 1. Daubechies, "The Wavelet Transform: a Method for Tirne-
Frequency Localization" in Advances in Spectral Analysis and
View publication stats