Signal representation using adaptive normalized Gaussian functions

SIGNAL PROCESSING Signal Processing 36 (1994) 1- l I

ELSEVIER

Signal representation using adaptive normalized Gaussian functions Shie

Qian*,

Dapang

Chen

DSP Group, National Instruments, 6504 Bridge Point Parkway, Austin, TX 78730-5039, USA

Received 17 February 1992 ; revised 3 August 1992 and 19 March 1993

Abstract In this paper, a new joint time frequency signal representation, the adaptive Gaussian basis representation (AGR), is presented . Unlike the Gabor expansion and the wavelet decomposition, the bandwidth and time frequency centers of the localized Gaussian elementary functions h p (t) used in the AGR can be adjusted to best match the analyzed signal . Each expansion coefficient B, is defined as the inner product s,(t) and h,(t), where s i lt) is the remainder of the orthogonal projection of s,_, (t) onto h,_ 1 (t). Consequently, the AGR not only accurately captures signal local behavior, but also has a monotonically decreasing reconstruction error I s,(t) 2 By combining the AGR and the Wigner Ville distribution, we further develop an adaptive spectrogram that is non-negative, cross-term free, and of high resolution . Finally, an efficient numerical algorithm to compute the optimal Gaussian elementary functions h,(t) is discussed . I',

.

Zusammentassung Im folgenden Beitrag wird eine neue Zeit-Frequenz-Signaldarstellung vorgestellt, ndmlich die adaptive GauBBasis-Darstellung (AGR). .Anders als bei der Gabor- oder Waveletzerlegung konnen die Bandbreiten and ZeitFrequenz-Zentren der lokalisierten GauB'schen Elementarfunktionen h i lt), die in der AGR verwendet werden, so angepal3t werden, daB sie dem analysierten Signal bestmoglich entsprechen . Jeder Entwicklungskoeffizient B, wird als inneres Produkt von s,(t) and h,(t) definiert, wobei s pit) der Rest der Orthogonalprojektion von s, t (t) auf h, 1 (q ist . Daraus folgt, daB die AGR nicht nor das lokale Signalverhalten in exakter Weise erfaRt . sondern auch einen monoton abnehmenden Rekonstruktionsfehler SP(t) II 2 aufweist . Durch die Kombination von AGR and Wigner-Ville-Verteilung entwickeln wir daruber hinaus ein adaptives Spektrogramm, das nichtnegativ ist . frei von Kreuztermen and hochauflosend . SchlieBlich wird ein effizienter numerischer Algorithmus zur Berechnung der optimalen GauB'schen Elementarfunktionen h,(t) diskutiert. Resume Cet article presente one nouvelle representation conjointe temps frequence d'un signal : la representation adaptative de base gaussienne (AGR) . A la difference de ('expansion de Gabor et de la decomposition en ondellette, la largeur de bande et les centres temps-frequences des functions gaussiennes elementaires h,(t) utilisees dans le AGR pcuvent titre ajustes pour s'ajuster au mieux au signal . Chaque coefficient d'expansion B, est defini comme le product scalaire de s,(t) et de

*Corresponding author . 0165-1684/94,1$7 .00 'c", 1994 Elsevier Science B .V . All rights reserved SSDI0165-1684(93)E0061-0



2

S. Qian, D . Chen / Signal Processing 36 (1994) 1-/i

on s p (t) est le reste de la projection orthogonale de s,_ ,(t) sur h,-,(t) . En consequence, non seulement l'AGR reproduit efficacement les caracteristiques locales du signal, mais en plus, I'erreur de reconstruction s p (r)i decroit de maniere monotone, En combinant I'AGR et la distribution de Wigner-Ville, nous developpons ensuite on spectrogramme adaptatif non-negatif, sans terme croise et de haute resolution . Finallement, tin algorithme numerique efficace permettant de calculer les fonctions gaussiennes elementaires de maniere-optimal est discute h,(t),

Key words:

Gahor expansion ; Gaussian function ; Inner product ; Orthogonal; Wavelet; Wigner-Ville distribution

1 . Introduction Because the majority of signals encountered in the real world are of a time-varying nature, it is important to explore their local behavior simultaneously in time and frequency . The most popular current methods for characterizing time-varying signals are the Gabor expansion [5] and the wavelet decomposition [4,6] . While the bandwidths of the analysis and synthesis functions are fixed in the Gabor expansion,' the time and frequency resolutions of the elementary functions in the wavelet decomposition are regulated in a uniform manner . Both the Gabor expansion and the wavelet decomposition are defined on a regular sampling grid . In many signal processing applications, however, it is desirable to have flexible elementary functions, in order to accommodate different components of the analyzed signals . A decomposed signal can be perfectly reconstructed from the coefficients of its Gabor expansion or wavelet decomposition . However, unless the elementary functions used constitute an orthonormal base, such representations may be redundant . In some applications, such as data compression and noise reduction, accurate representation is not necessary, and a certain degree of distortion is acceptable . The major objective may be to efficiently represent the interesting components (according to some criterion) in the signal, while retaining sufficient information to reconstruct the signal with a limited amount of distortion . In this paper, we present a signal representation scheme that uses adaptive normalized Gaussian functions, which we term the adaptive Gaussian

basis representation (AGR) . Unlike the Gabor expansion and the wavelet decomposition, in the AGR the time and frequency resolution and time-frequency centers of the normalized Gaussian elementary functions hp (t) are adjusted to best match the signal . The coefficients B o are defined as inner products of the signal and Gaussian elementary functions, thus reflecting the signal behavior at the time-frequency instant directly associated with the localized elementary functions . Numerical simulations indicate that the AGR is economical in representation and produces reliable results in the presence of random noise . An independent study by Mallat and Zhang has obtained a very similar result using a matching pursuit algorithm [8] . A dictionary of a large set of elementary functions (atoms) is first constructed . At each stage, a match is obtained between one of the atoms and the residual signal . The process continues until the residual signal becomes negligent . In Section 2, the framework of the AGR is introduced, followed by a demonstration that the reconstruction error of the AGR monotonically decreases and converges to zero . In Section 3, we combine the AGR with the Wigner-Ville distribution to construct an adaptive spectrogram - a new signal energy distribution in the joint time -frequency domain . This new distribution is non-negative, cross-term interference free, and of high resolution . Section 4 is devoted to a fast numerical algorithm which estimates the optimal Gaussian elementary functions . Finally, numerical simulations are presented to demonstrate the effectiveness of the approach .

2. Adaptive Gaussian basis representation 'The methods of determining the Gabor coefficients C,„,,, are generally not unique . In this paper, C,,, is computed by the biorthogonal method [1] .

The goal of the AGR is to map a signal of interest s(t) into a space constituted by a class of localized



S. Qian, D . Chen , Signal Processing 36 (1994) 1 11

elementary functions h,(t) : s(t) = Y B,h,(t) .

(I)

P -0

In order to best characterize the time-varying nature of s(t). we let the elementary function h,(t) be the normalized Gaussian function with an adjustable variance a P and time frequency center (t,, f,) :

reflect the signal local behavior . In the following, we describe a procedure for computing B,, a„ t, and f, . Begin with stage p = 0. The parameters a n , t, and J P are chosen such that h p (t) is most similar to s,(0, i .e . . B,1 2 = max,,,,, .r,

2

_ (7[x

)! P ER"

25

exp

t,,frclB .

-

t2 2a, ' (2)

I TF = 4,c'

where' t 2 1 gp(t)1 2 dt

J and (3)

h,

,(t):

sp(t) =

s,-1

(t)-B,-Ih,(5)

Since h p(t) is subject to normalization ( II h,(t) = 1 . s, t (1)11

12 =

1),

(6) t 1 2, where II - 2 denotes the Euclidean norm . Based on Eqs . (5) and (6), one can perform the following multistage decomposition . At each stage, the coefficients Be are determined according to Eq . (4), see Table 1, Such a decomposition is depicted in Fig . 1 . After Q stages of decomposition, the following relationships hold : 2

J

and where Qt) denotes the Fourier transform of g,(t) . The quantities T and F are traditional time and frequency resolution measures (energy variances) of g,(t), respectively. In Eq . (2), the variance a, determines a tradeoff between the time and frequency resolutions of g,(t) as governed by Eq . (3). In our approach, the parameters a„ t, and f, are chosen such that h,(t) is most similar to s(t) on local basis. Since the time- and frequency-shifted Gaussian functions do not constitute an orthogonal basis, the selection of the coefficients B, is generally not unique . To characterize the time-varying nature of s,(t), the elementary functions h,(t) must, as discussed in detail later, be localized in the joint time-frequency domain, and the coefficients B, must

(4)

where s n (t) = At) . For p >_ I, s,(t) is the remainder after the orthogonal projection of .s,-,(t) onto

Is ,(t)1

F 2 = J f 2 I Gg(t)2 I dJ,

t,/'PER .

x,ER' ,

The Gaussian functions g,(t) are not compact in either time or frequency, but they do possess the best joint time-frequency resolution in the sense of having the minimum time-bandwidth product

T2 =

f S,(t)hp(t)dt

=maxap,,, .f,,

h,(t) = qp(t - t,) exp{j2ttJp t},

g,(t)

IFP1 2

2 - I Bp

Q

s(t)

_ L, B,h,(t) +

sQ+,

(7)

(t)

,-a

and Q Is(1)11

2

=~,IB,I

p=o

2

+I':sQ+1(t)11

(8)

2

In the following, we will show that the norm of the residual 5 Q _ 1 (t) monotonically decreases and converges to zero . Let 0, be the angle between s,(t) and h,(t) . Then cos O,

,,, h,) _ _ s, f'

BP

(9)

IIs,II

Substituting Eq . (9) into Eq . (6) yields I

s,(t)11 2 =

( sin 0,_ 1 ) 2

Obviously, 11 s,(t) 1 terms of is,-201

2

II s,-1(1)11 2 .

can further be represented in

2 , e .g .,

2

Unless otherwise indicated, the limits of integration are assumed to extend from - :z, to +x .

I s,(t)

)I 2 =

( sin 0,

1) 2 (sin

Op

2) 2

I s,-

2(t)1' 2 .

(10)



4

S . Qian, D. Chen Signal Processing 36 (1994, 1-11

Table I Stage

Remainder

0

s o (t) = s(t) s,(t) -

P

son)

Remainder energy 1ISoy) - B,ho(i)

s,tt)=sp-,(t)-Bp-,ho-,(t)

i

t

Fig . l . The coefficient B p is the orthogonal projection of s, on h, . s,„ is the remainder after the projection of s, onto h, .

Therefore, it can be shown that after Q stages, Q SQ+1(t)11 2- Iso(t)Il i

]I

(sin 9 ,)2

v=t

11S(t)11 2 (sinOma .) 2Q ,

(11)

where (sin 0 ma „) 2 >(sin O,) 2 , Vp. Since (sin 0,)' is strictly less than one (otherwise, h,(t) is orthogonal to s,(t) and Be = 0), the sequence U s(t)II '(sin Oma,) 2Q monotonically decreases and converges to zero as Q increases without bound . Thus, Eq . (8) can be written as IIs(t)11 2

= I

IB,1 2 .

P=o

In summary, for signal s(t), the error-free AGR is

B,h,(t),

s(t) _ P=o

B,= f s,(t)It (t) dt, s,(t)

=s ,-,

(t) - Bp-t h,-, (t),

1!1s1(t)11 211

(13)

11s,(t) I ' =

Is(tyI

2

11so(t)11' - IBoI'

1 2 - iB,,-,1 2

Is,-,( )11

where the optimal Gaussian elementary functions h,(t) are determined according to Eq . (4). The main features of the AGR are summarized below . 1 . The AGR faithfully depicts the signal local behavior, because : (a) The elementary functions h,(t), adaptive normalized Gaussian functions, are optimally concentrated in the joint time-frequency domain ; (b) The coefficient B„ the inner product between s,(t) and h,(t), is an indication of signal behavior at the time-frequency instant, [t, - T„ t, + T,] x [f, - F, J, + F,], directly characterized by h,(t) . (Although the {h,(t)} are not even linearly independent, the AGR is similar to the orthonormal representation . The merit of this feature will be further explored when the AGR is used to construct the spectrogram in the next section . 2 . By regulating the variance a„ the constituent function h,(t) can be made to track a variety of signals. In the AGR, the bandwidth of g,(t) is automatically adjusted to best match the signal . The difficulty in selecting a window bandwidth that is encountered in the Gabor expansion as well as the wavelet decomposition is avoided . 3 . The error term 11sQ+111 2 monotonically decreases as Q increases, which makes it easier to balance the error tolerance and the computation complexity . In this scheme, adding the new basis function h Q+t (t) does not affect the previously selected parameters, B, and h,(t), for 0 < p <, Q . 4 . Because of its adaptive nature, a signal in the AGR is generally concentrated in a very small subspace . Consequently, we can use a subspace representation Eq. (7) to approximate the signal s(t) . If random noise tends to spread 'evenly' in



S . Qian . D . (hen Signal Processing 36 1994) 1-11

Fig . 2 .a Frequency Hopper

Fig . 3 .a FM Signal

the entire space, the subspace representation will effectively improve the SNR . Although this subspace noise reduction principle is generally true for many transformations, it is more efficient in AGR due to the highly concentrated signal energy . Figs . 2(a) and 3(a) illustrate a 256-point fourtone frequency-hop signal and a sampled FM signal, respectively . Figs . 2(b) and 3(h) depict the relationship between the number of elementary functions Q, and the reconstruction error

S

I S(t) II` NRQ=10log dB . s Q +, (t) I As expected . the reconstruction error monotonically decreases as the number of elementary functions Q increases . For the four-tone frequency-hop signal, using the first four elementary functions to approximate the original signal yields an SNR as high as 20 dB .

a

Fig . 2 .b SNR (dh) vs . Number of Elementary Functions

Fig . 3 .6 SNR (dh) vs . Number of Elementary Functions

In the next section, we will further explore the advantages of the AGR by a'nplying it in a WVD decomposition, which yields an adaptive spectrogram - a new signal energy distribution in the joint time-frequency domain .

3 . Energy distribution in the joint time-frequency domain The Gabor expansion, the wavelet decomposition, and the AGR are all linear transformations of the signal s(t) . There is another very useful timelocalization method, the WVD E12,13], which is quadratic in the signal s(t) . The WVD for the signal s(t) is defined as W VD 5 (t, f 1 = Is(t+r/2)s*(t-z,12)e -)21 `da .

(14)



6

S . Qian, D. Chen / Signal Processing 36 (1994) 1-11

Although it possesses many useful properties for signal processing [2, 3], the WVD suffers crossterm interference and problems with negativity that greatly impede its useful application . Recently, the decomposition of the Wigner-Ville distribution in terms of the WVD of a complete set of Gaussian functions was investigated [10] . That work was motivated by the fact that the WVD of the Gaussian function is non-negative and well localized in the joint time-frequency domain . The basic procedure is first to decompose the signal s(t) into a linear combination of time- and frequencyshifted Gaussian functions gp (t), such as s(t)_

Y

Bp9 P (t-t4)e 2ifpl .

( 15)

p=0

Then, taking the WVD with respect to both sides of Eq . (15), yields WVD,(t,f)=

IB P I 2 WVD g ,(t - tp,J

Y

- f')

a signal s(t), the Gabor expansion is defined as' K

s(t)

_ Y

."= . C

1]

mAT)ei 2 "" Apt

Cm."9(t -

f s(t)y*(t-,AT)e - 2u " AP, dt. i

(18)

The existence of Eq . (17) is only possible when AT AF <, 1 [4] . In the orthogonal-like Gabor expansion, the analysis function y(t) is chosen to be similar to the synthesis function g(t), i .e., y(t) z flg(t), where /1 is a real constant . Consequently, the Gabor expansion Eq . (17) has an orthogonal-like form, although the functions {g(t-mAT)e are not even linearly independent . Replacing t P by m AT and fp by n AF in Eq . (16), the CDR based on the orthogonal-like Gabor expansion has the following closed form [10] : 11 CDR(t,f)= Y E ICm,"I 2

P=0

x WVDq (t-mAT, f- nAF) Re(1B 0 8g1

+2 Y

J~ P-oq>

x WVD,"9"

=2

`p

I

t - tP 2 tq .l -

'p 2 fq)l J

(16)

The first sum on the right-hand side of Eq . (16) is a linear combination of time-frequency shifted Gaussian functions WVD, which are non-negative and concentrated at (t,, fL ) . The second sum consists of the cross-Gaussian functions WVD, which are oscillating at (](t, + 1,), (f, + J,)) and do not have physical interpretations. If energy contained in the cross-terms is negligible, the auto-terms alone can reasonably be used to represent a signal spectrogram. This is called a cross-term deleted representation (CDR). The CDR is non-negative and cross-term free . In general, the time- and frequency-shifted Gaussian functions do not constitute an orthonormal basis . Therefore, the choice of B P in Eq . (15) is not unique . Indeed, the selection of B,, remains an open question in designing the CDR . One possible scheme is to decompose s(t) via the orthogonal-like Gabor expansion [9] . For

f

xexp

ICm ."I 3 - (t-m

AT) 2

R

-(2nQ) 2 (f-nAF) c ',

(19)

where I W VD g (t, f) 1 2 = 11911' = t . The necessity of seeking y(t) x #g(t) is demonstrated by a simple example that follows . Assume that the analyzed signal s(t) is concentrated around (m e AT,n 0 AF). It is natural to expect that IC m ,,," P Z , the weight of WVD g (tm o AT, f- n o AF), is large . If the analysis function y(t) significantly differs from the synthesis function 9(t), however, ICm,,," q 1 2 = I < S(t), ym" .""(t) > Ir may not be large . The peak may appear elsewhere . Consequently, the resulting CDR may fail to reflect the signal energy distribution .

'The index p is dropped from g e(t) here because the variance a„ = a', of the Gaussian functions is fixed in the Gabor expansion.



7

S . Qian . D . C'hen) Signal Processing 36 '1994) 1 11

(a)

(b)

(c)

(d)

(e)

Fig . 4 . (a) Power spectrum of frequency hopper, (h) W igner-Ville Distribution (5 .98 s), (c) squared STFT with Gaussian window (6 .33 s), (d) CDR (1 .01 sl. (e) ADS (4.24 s). For frequency hopper signal, the adaptive spectrogram has best resolution in the joint time-frequency domain .

As our extensive experiments have indicated, the CDR based on the orthogonal-like Gabor expansion has certain advantages, such as computational economy, it is cross-term interference free, and it has good time frequency resolution . On the other hand, to achieve °:(t) z f y(t), oversampling is required, which causes resolution deterioration . Analogous to the CDR, taking the WVD on

Eq . (13) and removing the cross-term yields the so-called adaptive spectrogram (ADS) : ADS(t,J)=2 j IBphexp

-(t

p°o - (2 n) 2 ap(f C( p E

P8 5 ,

tp E

R,

fp c [0,

lp)Z ar) .

tpt a 1,

, (20)



8

S. Qian, D . Chen / Signal Processing 36 (1994) 1-11 (a)

(b)

(c)

(d)

(e)

Fig. 5 . (a) Random noise corrupted frequency hopper signal, (b) Wigner Ville distribution, (c) squared STFT with Gaussian window, (d) CDR, (e) ADS . In general, the cross-term deleted representation is less sensitive to noise, particularly, the ADS, because in ADS the signal is highly concentrated in a small subspace .

Since the exponential term in Eq . (20) is a WVD for a normalized Gaussian function, it has unit energy . Using Eq . (12), we can readily verify ~ADS (t,1)112= V IBnI2=Is(t)I2 . on

(21)

Thus, the energy contained in ADS(t, f) is identical to the energy contained in the signal s(t) . Hence, the ADS can be considered as a signal energy distribution in the joint time-frequency domain . This new distribution is non-negative, cross-term interference free, and of high resolution .



S . Qian . D . Chen

Signal Proeessin,e 36 (1994) 1- 11

Figs . 4 and 5 illustrate the computation of the spectrogram of a clean and a noise-corrupted frequency hopper signal via different approaches . The numbers in the parentheses indicate the computation time incurred using a Macintosh Ilci . The ADS is not only cross-term interference free, but also has very good resolution . Moreover, the ADS is less sensitive to noise, because in the AGR the signal is generally concentrated in a small subspace .

4. Parameter estimation of optimal Gaussian elementary functions As shown earlier, one can implement the AGR in an iterative manner . At the stage p, first find the solution of Eq . (4) : xp , tp and fp ( the algorithm for solving Eq . (4) will be introduced later). Then compute the coefficient B, and the remainder .s, 1 . Check the remainder ' s n ,-, 2 . Terminate the process if the remainder is acceptable, or otherwise go to the next stage p + 1 until the desired reconstruction error is achieved . The major difficulty in implementing the AGR is the computation of the optimal elementary function hp ( t) defined in Eq . (4), which is a multiextrema optimization problem . The error surface of Eq . (4) is rather complicated . So far no analytical solution has been found. Although some optimization methods that deal with multiextrema functions, such as simulated annealing [7, 11], have been proposed, the computational complexity prevents them from being applied to many applications . In the rest of this section, we introduce an efficient numerical algorithm for the determination of the optimal normalized Gaussian elementary functions hp (, t). Fp in Eq . (4) is in fact a Fourier transformation of spy) with the window g,(t - tn ) . B p (the largest Fp ) can be obtained by first computing F p for all possible a p , tp, fa and then choosing the largest Fp . In practice, the parameters 1 . , t p and f;, are usually digitized in order to utilize the digital computers . For example, in our algorithm, the variance x has k n + 1 levels and the time-frequency centers are limited at tp = m AT and f, = n AF, where A T and AF are time and frequency sampling intervals, ni, n e Z . AT and AF define a tradeoff between the

9

accuracy and the computational burden . We let A 7' = T and AF = F. T and F, which are defined in Eq . (3), can be considered as a window effective m n) can width and frequency bandwidth of g(t). Fz,(m be estimated by

Fp (m, n) =

f

sp(t)gp(t -m AT,,)

xexp[-j21rnAFp t ;dt,

(22)

which is similar to the analysis formula of the Gabor expansion Eq .(18) sampled window Fourier transform, where AT, and AFp correspond to the time and frequency sampling intervals at the pth stage . Since ATp AF,= Tp Fp = 1 /4n < 1 (critical sampling rate), FF (m, n) contains enough information to recover the signal s p(t) . The variance a of g(t) determines the effective window width 7 . F Iin, n) roughly represents the signal's time frequency behavior at t„ = mATo and fp = n A F'p . Analogous to zoom processing, we first use a window, g (i), with a larger variance x (thereby poor time resolution) to scan the data and roughly locate the peak (similar to using a wide angle to catch the object of interest from the large area) . Then we gradually reduce the effective window width (by using a smaller variance a, similar to using a smaller angle) and the search area, until the peak is out of focus . Since at first, the variance a is large, the time increment AT is also large, and therefore a few samples in the time domain can determine the area of interest . As the variance a decreases, the search becomes more detailed in the time domain but the search range also shrinks at the same time . In the simulations, the variance a was determined by (the selection of kp depends on the application) ak =

4` 221k o kl,

k=0,1,2, . . .,k 0 .

R/

This particular choice is very similar to a wavelet representation where the scale is halved in each step . The normalized Gaussian elementary function (window function) at kth level is gk(t)_(n)!k)-O .2 5exp

-

tz 27.

(23)



10

S . Qian, D . Chen / Signal Processing 36 (1994) 1-11

j 1,IAA& '1. , O AT O

I,

m,dT,

m2AT2 Fig. 6 . Matching s,(d with multi-resolution windows g &(t) .

The time increment at kth level is set to be ° k l , ATk=Tk=\0 .5ak=2 k which is halved as k is increased by one . Eq . (4) shows that the estimation of the parameters aa , t, and f, is equivalent to matching the signal s,(t) with all possible elementary functions h(t) in the time-frequency domain . The computation is done by comparing the squared Fourier coefficients, IF,(m, n)J 2 , of sp (t)gk (t - mAT,) at all levels k . Fig . 6 illustrates how the search is done and how the parameters are estimated, where a bold line indicates a Gaussian function gk(t1 k ATk ) . Such that the maximum squared Fourier coefficient is found . For k = 0. the Fourier transforms for s,(t)gn (t - m ATo ) must be evaluated at all times instances m ATo . On the other hand, since ATo > ATk for k > 0, only a few time sampling steps will cover the entire time record . Suppose that a match is found at the time instance m r, ATn , it can be reasonably expected that at the next level k = 1, the best match will occur in the vicinity of m 0 ATo because g n (t) is localized . Consequently, while the searching in time domain goes to more detail (TI = i To ), the search at the level k = 1 is limited to a smaller period around t=m 0 ATp . This process is continued until the

maximum squared Fourier coefficients no longer increase (say at k = K). Then, the variance and time-frequency center corresponding to the maximum squared Fourier coefficient of s,(t)gk(t) are the solutions a,,, t o and f, . After h,(t), the best match, is found, the remainder s o+ ,(t) is given by sp+1(t) = sv(t) - hp(t) .

Note that s, + , (t) differs from s,(t) only for t close to t„ the time center of h,(t), because of the localization of h,(t) . Once s,+ 1 (t) is found, in step p + 1, % p+ 1, t, +1 and fp+1 can be determined in the same way . Continue this process until the error, s Q (t) 1 I 2 , is acceptable. It should be pointed out that the sampled window Fourier transform of s,(t) can be computed by FFT if s,(t) and 9k(t) are digitized as s,(i) and g k (i), respectively_ If 9k(i) is considered as an L-point window (regardless of the effective window width Tk ), then the computation required to determine each optimal elementary function is approximately [9] L log, k =3L(K+2){log 2 L-z(K+1)}, k_o 2

K+1

3L

where K s k o is the number of levels needed to find the optimal elementary function h,(t) .



S . Qian . 1) . Chen ! Signal Processing 36 (1994) 1 -11

The efficient algorithm presented in this section produces an approximate solution of Eq . (4), but the results, as our simulations indicate, were very close to those given by more complicated algorithms [7, 11] . In addition, the algorithm presented here guarantees convergence and takes a finite number of steps to terminate as opposed to simulated annealing algorithms, whose convergence is not always guaranteed and for which the execution time is very sensitive to the signal . Most strikingly, a simulated annealing algorithm would require almost an hour of VAX CPU time to process a 256-point frequency-hop signal, whereas the algorithm presented here uses less than three seconds on a Macintosh Ilci computer!

5 . Conclusions Because the Gaussian function is optimally concentrated in the joint time-frequency domain, it is natural to use data-adaptive Gaussian functions as the elementary functions to represent the timevarying signal-' Time-shifted and frequencymodulated Gaussian functions, however, do not constitute an orthonormal base . It has been an open question as to how to select the coefficients and to what extent the resulting coefficients represent the analyzed signal . In the AGR, the data-adaptive localized Gaussian functions h,(t) are used as the elementary functions, and the coefficients B, are defined as the inner product of s,(r) and h,(t), where s,(t) is the remainder from the orthogonal projection of s,_,(t) onto h,_ t (r). Consequently, the AGR not only captures local signal behavior, but also has a monotonically decreasing reconstruction error. Its„,(t)I12 . Based on the AGR and Wigner Ville distribution, we further developed an adaptive spectrogram-a new signal energy

'As one reviewer and our experience indicated, the elementary function could further accommodate the chirp type signal . if the variance a, is allowed to be complex . Such innovation, however, would tremendously increase the computation needed to determine the optimal elementary functions .

11

distribution in the joint time-frequency domain, which is non-negative, cross-term interference free, and of high resolution . Finally, an efficient algorithm to estimate the optimal Gaussian elementary functions was discussed .

References [1] MT - Bastiaans, "On the sliding-window representation in digital signal processing", IEEE Traits, Acousr . Speech Signal Process., Vol. ASSP-33, No. 4, August 1985 . pp .868 873 . [2] T .A .C .M . Claasen and W .F .G. Mecklcnhrduker, The Wigner distribution A tool for time-frequency signal analysis", Philips J . Res ., Vol . 35 . 1980 . pp 217-250. 276-300, 1067 1072 . [3] L . Cohen, 6 Time-frequency distribution A review', Proc. IEEE, V o l . 77 . No. 7, July 1989 . pp. 941-981 . L4J 1 . Datibechies, 'The wavelet transform, time frequency localization and signal analysis", IEEE Trans . Inform . Theart . September 1990, pp . 961-1005 . [5] D . Gabor . "Theory of communication" . .1 IEE ILondon), Vol . 93. No . ill, November 1946, pp . 429-457 . [6] P. Goupillaud, A . Grossmann and J . Modet,'Cycle-octave and related transformation in seismic signal analysis", Geoexplorarion, Vol . 23, 1984, pp-85-102 . [7] S . Kirkpatrick .. C .D. Gclatt and M .P . Vecchi . "Optimization by simulated annealing", Science . Vol . 220, 13 May, 1983 . pp . 671-680. [8] S . Mallat and Z. Zhang, "Matching pursuit with time frequency dictionaries", IEEE Trans. Signal Process., December 1993 . [9] S . Qian and D. Chen, "Discrete Gabor Transform", IEEE Trans . Signal Process .. Vol . 41, No . 7, July 1993, pp . 2429 2439 . [10] S . Qian and J .M . Morris . "Wigner distribution decomposition and cross-term deleted representation", Signal Proce .ssing, Vol . 27 . No . 2 . May 1992, pp. 125-144. [I I] H . Szu and R . Hartley, "Fast simulated annealing", PhDs . Left ., Vol . 122, Nos . 3, 4, 8 June 1987, pp . 157--162 . [t2'] l . Ville, "Theorie et applications de la notion de signal analylique" . Cables et Transmission, Vol . 2, 1948, pp . 61-74 (in French) . [13] E.P . Wigner, "On the quantum correction for thermodynamic equilibrium", Phis . Ree ., Vol . 40, 1932, pp . 7496