Evolutionary chirp representation of non-stationary signals via Gabor transform

Signal Processing 81 (2001) 2429–2436

www.elsevier.com/locate/sigpro

Short communication

Evolutionary chirp representation of non-stationary signals via Gabor transform Ayd)n Akana ;1 , Luis F. Chaparrob ; ∗ a Department

b Department

of Electronics Engineering, University of Istanbul Avcilar 34850, Istanbul Turkey of Electrical Engineering, University of Pittsburgh, 348 Benedum Hall, Pittsburgh, PA 15261, USA Received 2 August 1999; received in revised form 28 March 2001

Abstract In this paper, we propose a chirp time–frequency representation for non-stationary signals, and associate with it— via a multi-window Gabor expansion—the corresponding evolutionary spectra. Representations based on rectangular time–frequency plane tilings give poor time and frequency localization in the spectrum, especially when the signal is not modeled well by 7xed bandwidth analysis. We propose a representation that uses scaled and translated windows modulated by chirps as bases. Considering a chirp-based Wold–Cramer model, the signal evolutionary spectrum with improved time and frequency resolutions is obtained from the kernel of the representation. The chirp representation optimally chooses scales and linear chirp slopes by maximizing a local energy concentration measure. Parsimonious signal representation and well-localized evolutionary spectrum are obtained simultaneously. As an application of our representation, we consider the excision of broad-band jammers in spread spectrum communications. Examples illustrating the improvement in the time and frequency resolution of the signal spectrum using our procedure are given. ? 2001 Elsevier Science B.V. All rights reserved. Keywords: Evolutionary spectrum; Time–frequency analysis; Discrete Gabor expansion; Time–frequency signal representation; Spread spectrum

1. Introduction Time–frequency representations characterize signals in terms of joint time and frequency content, and are useful in dealing with signals with ∗

Corresponding author. Tel.: +1-412-624-9665; fax: +1-412-624-8003. E-mail addresses: [email protected] (A. Akan), [email protected] (L.F. Chaparro). 1 This work was partially supported by The Research Fund C of The University of Istanbul, Project numbers: O-436=060598 and B-428=13042000.

time-dependent spectra, such as speech, audio and biomedical signals. When developing a representation for a non-stationary signal, it is important not only to select an appropriate basis for it, but also to associate with the representation a time-dependent spectrum. Such is the case of the short-time Fourier transform (STFT) [5] and the Spectrogram associated with it, or of the Wold–Cramer representation [14] and Priestley’s evolutionary spectrum (ES) [13,14]. Recently, the Gabor representation [15,17], of which the STFT is a special case, has been connected with the evolutionary spectrum [1].

0165-1684/01/$ - see front matter ? 2001 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 0 1 ) 0 0 1 3 1 - 1

2430

A. Akan, L.F. Chaparro / Signal Processing 81 (2001) 2429–2436

The Wold–Cramer representation [14] of a random non-stationary signal v(n) consists of an in7nite sum of sinusoids with random and time-varying magnitudes and phases
 A(n; !) ej!n d Z(!): (1) v(n) = −

The evolutionary spectrum of v(n) [13,14] is then given by S(n; !) = (1=2)|A(n; !)|2 , for which different estimators have been proposed [9,10]. An analogous representation for a deterministic signal x(n) with a time-varying spectrum can be de7ned as x(n) =

K−1 

A(n; !k ) e j!k n ;

(2)

k=0

where !k = 2k=K, i.e., a sum of time-varying sinusoids at discrete frequencies. Most of the existing time–frequency representations use sinusoidal bases. Since these representations use the same time and frequency resolutions for the whole time–frequency plane, resulting in a rectangular tiling, they are not adequate for signals with wide-band components. These representations are parsimonious only when the signal is composed of narrow-band signals [2]. Improvement in the spectral resolutions can be achieved by combining rectangular tiling representations [1]. As an extension of the traditional Gabor expansion, the multi-window representation of a 7nite support signal x(n) is given by x(n) =

K−1 I −1 M −1 1   ai; m; k h˜i (n − mL) e j!k n I k=0 i=0 m=0

0 6 n 6 N − 1;

(3)

where I is the number of scaled windows, M and K are the number of time and frequency samples considered, L and L are the time-step and frequency-step, and !k = 2k=K. We also have that ML = KL = N , and L and L are chosen so that LL 6 N , or L 6 K, for numerical stability. The case when L = K is called critically sampled, while when L ¡ K is the oversampled case [1,15,17]. The window h˜i (n) is a periodic extension of the synthesis window generated by contracting a mother

Gaussian window g(n): hi (n) = 2i=2 g(2i n);

0 6 n 6 N − 1;

(4)

for the diLerent scales i ∈ I . The Gabor coeMcients are given by ai; m; k =

N −1 

x(n)˜∗i (n − mL) e−j!k n ;

(5)

n=0

where the analysis window ˜i (n) is bio-orthogonal to h˜i (n) [17]. In the critically sampled case, ˜i (n) is unique but not necessarily well-localized, while in the oversampled case it is not unique, but can be chosen well-localized. The evolutionary spectrum of x(n) is found [1] as 1 SE (n; !k ) = |A(n; !k )|2 ; K where the evolutionary kernel A(n; !k ) is A(n; !k ) =

N −1 

x(‘)w(n; ‘) e−j!k ‘ ;

‘=0

and the time-varying window w(n; ‘) depends on the analysis and synthesis windows and is much better localized than either of them. Choosing a non-sinusoidal basis, the time-varying frequency content of the signal can be captured by a representation with fewer coeMcients than using sinusoidal bases. This parsimony comes at the expense of selecting appropriate chirp bases, a diMcult problem. In previous solutions to this problem, it has been argued that this parsimony can be achieved by using linear-chirp-modulated bases [2,4,8,11]. The best matched bases in some of these methods are searched in a redundant dictionary of basis functions, requiring computationally expensive iterations. Mann and Haykin [12] have proposed a continuous chirplet transform (CCT) with parameters time-shift, frequency-shift, scale, chirp-rate and dispersion rate to represent a continuous-time signal. In this paper we propose to use scaled, time and frequency shifted and chirp modulated Gaussian windows to obtain a representation for non-stationary discrete signals similar to the short-time Fourier transform with time-varying windows and chirp bases. A proposed optimal procedure permits us to search locally for appropriate scales and chirp-rates that best 7t the signal.

A. Akan, L.F. Chaparro / Signal Processing 81 (2001) 2429–2436

The rest of the paper is organized as follows. In Section 2, we propose a Wold–Cramer representation for deterministic signals in terms of chirps. In Section 3 we consider an optimal approach to obtain the chirp-based representation and propose an algorithm to choose the scales and chirp slopes. In Section 4 we propose the application of our chirp representation to jammer excision in spread spectrum communication systems and show examples that illustrate the evolutionary spectral analysis. Conclusions and discussion of the results follow. 2. Wold--Cramer chirp representation

p

 K−1  p

Ap (n; !k ) ej(!k n+p (n)) ;

tion is achieved by the non-rectangular tiling provided by the chirp basis. Calculation of the Gabor coeMcients given by aqs; l =

N −1 

x(n) ˜∗ (n − sL) e−j(n!l +q (n)) ;

(8)

n=0

˜ } be requires again that analysis windows {(:) ˜ }. bi-orthogonal to the synthesis windows {h(:) Comparing the Gabor representation (7) with the chirp Wold–Cramer representation (6), the evolutionary kernel is given by  p ˜ − mL) am; k h(n (9) Ap (n; !k ) = m

An extension of the sinusoidal Wold–Cramer representation for x(n), can then be given by  xp (n) ejp (n) x(n) = =

2431

(6)

k=0

where again !k = 2k=K and {p (n)} are phase functions corresponding to instantaneous frequencies (IFs) {!p (n)}, or slopes, of the chirp components xp (n). Thus, x(n) is a sum of chirp components xp (n) with slopes {!p (n)} and sinusoidally represented. When {p (n) = 0} this representation coincides with the sinusoidal. Although, whenever x(n) is composed of wide-band signals, the representation in (6) would be more parsimonious than the sinusoidal representation, the problem is in de7ning the chirp basis. Since attempting to estimate the IFs [3,5] of a multi-component signal x(n) is a very diMcult task, we will try a direct representation in terms of chirps. Consider then the extension of the sinusoidal Gabor representation given by (3). Assuming for simplicity only one scale, the discrete, 7nite-length signal x(n) can be represented by the following chirp Gabor representation: P−1 K−1 −1   M ˜ − mL) ej(n!k +p (n)) ; apm; k h(n x(n) = p=0 k=0 m=0

(7) where P is a large number of chirp slopes over the time–frequency plane. EMciency in the representa-

and so the evolutionary spectrum of x(n) is  2  1   SE (n; !k ) =  Ap (n; !k − !p (n)) :  K p

(10)

Although parsimony can be achieved with this representation, its implementation is diMcult because it is not clear how to de7ne the instantaneous phases of the chirps. One could, however, attempt to obtain such a representation by searching for regions in the TF plane where the dechirped signal x(n) e−jq (n) and the sinusoidal component xq (n) are very close. That is, to 7nd regions Rq in the time–frequency plane where the error xq (n) − x(n) e−jq (n) is minimal. By choosing for x(n) the representation of xq (n) ejq (n) in each of these regions, an optimal representation is obtained. This is the approach we consider in the following section. 3. Optimal chirp representation Our objective is to search for regions {Rp } in the time–frequency plane where evolutionary kernels, or corresponding Gabor coeMcients, can be found to make x(n) coincide with xp (n) ejp (n) for IFs !p (n) = !k + d p (n)= d n, p = 0; : : : ; P − 1 and k = 0; : : : ; K − 1. To implement this, we 7rst select a set of chirps with incrementally changing slopes that provide a good covering of the time– frequency plane. A sinusoidal Gabor representation of the dechirped signal x(n) ejp (n) is obtained for diLerent scales i = 0; : : : ; I − 1.

2432

A. Akan, L.F. Chaparro / Signal Processing 81 (2001) 2429–2436

We now optimally search for a set of Gabor coeMcients, and the corresponding evolutionary kernel, that best represents the signal. The optimization consist in choosing coeMcients that maximize an energy concentration measure in small regions of the time–frequency plane. After considering each of the I scales and the P chirp slopes, we obtain an array of Gabor coeMcients {ai;m;pk ; m = 0; : : : ; M − 1; k = 0; : : : ; K − 1}. The basic assumption in this search is that, at each time–frequency location, there is one set of these coeMcients that provides a better representation of the signal than any other set. Repeating this search over all regions in the time–frequency plane, we obtain the best representation possible. An optimal search method for multi-component signals requires a local measure that considers time– frequency localization in small regions. Jones and Parks [8] used a concentration measure to obtain the Gaussian window parameters which maximize the concentration of the spectrogram. This energy concentration measure is the ratio of the L4 -norm and the squared L2 -norm of the magnitude of the STFT. It is similar to the “kurtosis” and other concentration measures used in statistics. The following algorithm measures the energy concentration of Gabor coeMcients, {ai;m;pk ; m = 0; : : : ; M − 1; k = 0; : : : ; K − 1}, instead of the spectrum. The optimal coeMcients are used to synthesize the 7nal signal representation and to calculate its evolutionary spectrum. The energy concentration measure permits us to determine the most appropriate slope and scale parameters in small regions of the time–frequency plane. The optimization algorithm has the following steps: (i) Compute the Gabor coeMcient sets {ai;m;pk } for a set of chirp slopes {!p (n); p = 0; 1; : : : ; P − 1} and a set of window scales {2i ; i = 0; 1; : : : ; I − 1}. Divide the M × K coeMcient plane into blocks of size S × Q. (ii) In each of the S × Q blocks, compute the following local energy concentration measure: Ei; p (m; k) = S=2−1 Q=2−1

#i;2 p (s − mL; q − k)|ai;s;pq |4 2 ;  S=2−1 Q=2−1 2 (s − mL; q − k)|ai; p |2 # s; q i; p s=−S=2 q=−Q=2 s=−S=2

q=−Q=2

(11)

for each of the Gabor coeMcient sets {ai;m;pk }, p = 0; 1; : : : ; P − 1 and i = 0; 1; : : : ; I − 1. (iii) Obtain the Gabor coeMcient subset that gives the maximum energy concentration in each of the S × Q blocks. (iv) Use the optimal Gabor coeMcient to compute the evolutionary kernel Ap (n; !k ) using (9). Using this kernel the signal and the evolutionary spectrum can be calculated as indicated in equations (6) and (10). In (11) the time–frequency localization window #(·) is the Wigner distribution [7,8] of the analysis window (·) used to calculate the Gabor coef7cient ai;m;pk . As indicated in [7,8], #(·) is largest at its center and decays monotonically in all radial directions. Thus, only neighboring coeMcients aLect the concentration measure. We will show in the simulations that by using linear chirp basis functions for the time–frequency representation and combining the coeMcients that maximize the above concentration measure, we are able to obtain well localized evolutionary spectra. 4. Experimental results Example 1. As an application of the above representation, we consider the excision of wide-band chirp jammers in Spread Spectrum communications via time–frequency masking. Obtaining a well-localized spectrum for the received signal signi7cantly enhances the excision of the jamming signal. Direct sequence spread spectrum communication systems are widely used in transmission environments where high power jamming interference or multipath problems exist [16]. In such systems, the transmitted signal occupies a bandwidth much wider than necessary to send the message. Spreading of the message spectrum is accomplished by modulating the message with a spreading signal before transmission. In the receiver, the received signal is demodulated by a synchronized replica of the spreading signal. This permits us to recover the message and reduces the energy level of jammers and channel noise. This robustness of the system is

A. Akan, L.F. Chaparro / Signal Processing 81 (2001) 2429–2436

2433

Fig. 1. Evolutionary spectrum of (a) message, (b) corrupted message, (c) recovered (solid line), and corrupted (dashed line) signals, (d) error signal.

not strong for non-stationary, wide-band jammers. Jammer excision in SS communication systems is an important problem especially in military communications. Excising the jammers prior to despreading in the receiver increases the signal-to-noise ratio [16]. In our method, the evolutionary kernel of the corrupted received signal r(n) is computed with the proposed method. If the jamming signal is a combination of chirps with diLerent slopes, the evolutionary spectrum of the received signal will

be well-localized so that masking can be used to get rid of it. The masked signal can be obtained from Ay (n; !k ) = Ar (n; !k ) M (n; !k ); where Ay (n; !k ) corresponds to the kernel of the masked signal. The mask M (n; !k ) is obtained by using as threshold the average & of the peaks of Ar (n; !k ) in the absence of the jammer (this value can be computed in the regions where the jammer

2434

A. Akan, L.F. Chaparro / Signal Processing 81 (2001) 2429–2436

Fig. 2. (a) Evolutionary spectrum, (b) multi-window only evolutionary spectrum, (c) frequency marginal (dashed line), and magnitude-squared spectrum (solid line), (d) time marginal (dashed line) and magnitude-squared signal (solid line).

is not present)  0 Ar (n; !k ) ¿ &; M (n; !k ) = 1 otherwise:

(12)

The output can then be synthesized using Ay (n; !k ). Notice that the above procedure does not depend on the shape or number of components of the jammer, since our representation adapts locally to the signal. As an example, consider a 16-bit message sequence, and as a spreading function a pseudo-noise sequence with 16 chips. White noise (SNR of 2 dB), and a linear chirp jammer, with a

jammer-to-signal ratio of 18 dB, is added to the signal. Figs. 1(a) and (b) show the evolutionary spectra (in a dB scale) of the spread message and of the corrupted receive signal using parameters N = 256; L = 2; K = 256; I = 3, and P = 5. With no excision at the receiver, the bit error is 4 out of 16 bits. After masking the jammer as indicated above, the output signal is obtained and given (solid line) together with the distorted signal (dashed line) in Fig. 1(c). The error between the original and the recovered signals is given in Fig. 1(d). It is shown from the 7gures that linear-chirp jammer

A. Akan, L.F. Chaparro / Signal Processing 81 (2001) 2429–2436

2435

Fig. 3. (a) Evolutionary spectrum of speech, (b) combination of spectrograms in Jones and Baraniuk’s adaptive method, (c) frequency marginal (dashed line), and magnitude-squared signal spectrum (solid line), (d) time marginal (dashed line), and magnitude-squared signal (solid line).

is removed by the proposed method, improving the bit error performance (bit error = 0). Example 2. Consider a signal of length N = 128, composed of two linear chirps, one increasing and the other decreasing in frequency, a sinusoid and an impulse function. Most of the existing time– frequency analysis methods will not be able to adapt to the narrow- and wide-band frequency changes of this signal. Doing the analysis using I = 5 scales and P = 7 slopes, L = 1; K = 128 sampling inter-

vals, and the local energy concentration measure to combine the Gabor coeMcients, Fig. 2(a) shows the evolutionary spectrum obtained. As a comparison, Fig. 2(b) shows the evolutionary spectrum obtained with multi-window Gabor analysis with sinusoidal basis using the same Gabor parameters. Thus, by using a linear chirp signal model we are able to adapt to the time-varying frequency content of this signal. The chirp Gabor spectrum also approximates quite well the marginals. In Fig. 2(c) we compare the magnitude-squared spectrum (solid

2436

A. Akan, L.F. Chaparro / Signal Processing 81 (2001) 2429–2436

line) with the frequency marginal (dashed line) with excellent results. Likewise, the magnitude-squared of the signal (solid line) and the time marginal of the spectrum (dashed line) shown in Fig. 2(d) compare quite well. Example 3. We apply our chirp method to 100 ms. of speech signal sampled at 8 kHz, using I = 5; P = 7; L = 8; K = 2000. Fig. 3(a) displays the evolutionary spectrum obtained with our method. We compare this result with the ones obtained using an adaptive combination of spectrograms by Jones and Baraniuk [6] and displayed in Fig. 3(b). As shown, our method provides a better localized spectrum. The magnitude-squared signal spectrum (solid line), and the frequency marginal of the evolutionary spectrum (dashed line) are shown in Fig. 3(c). The magnitude-squared of the signal (solid line) and the time marginal of the evolutionary spectrum (dashed line) are shown in Fig. 3(d). As shown, the frequency marginal of the spectrum obtained with our method approximates the correct marginal due to its high localization. 5. Conclusions In this paper, we have developed a high resolution evolutionary spectral method which uses a linear chirp Gabor expansion. The fact that we are using a Sexible, non-rectangular time–frequency plane tiling allows us to obtain a compact Gabor representation and, at the same time, an evolutionary spectral estimate with high resolution. A given signal is analyzed using a set of analysis windows and chirp rates, and using a local energy concentration measure an appropriate set of Gabor coeMcients are selected. The resulting representation is well-localized in the time–frequency plane. The procedure can be applied very eMciently in the excision of wide-band jammers in spread spectrum communications.

References [1] A. Akan, L.F. Chaparro, Multi-window Gabor expansion for evolutionary spectral analysis, Signal Processing 63 (December 1997) 249–262. [2] R.G. Baraniuk, D.L. Jones, Shear madness: new orthonormal bases and frames using chirp functions, IEEE Trans. Signal Process. 41 (12) (December 1993) 3543– 3549. [3] B. Boashash, Estimating and interpreting the instantaneous frequency of a signal—Part 1: fundamentals, Proc. IEEE 80 (4) (April 1992) 520–539. [4] A. Bultan, A four-parameter atomic decomposition of chirplets, IEEE Trans. Signal Process. 47 (1999) 731–745. [5] L. Cohen, Time-Frequency Analysis, Prentice-Hall, Englewood CliLs, NJ, 1995. [6] D.L. Jones, R.G. Baraniuk, A simple scheme for adapting time–frequency representations, IEEE Trans. Signal Process. 42 (12) (December 1994) 3530–3535. [7] D.L. Jones, T.W. Parks, Time–frequency window leakage in the short-time Fourier transform, Circuits Signals Systems 6 (3) (1987) 263–286. [8] D.L. Jones, T.W. Parks, A high resolution data-adaptive time–frequency representation, IEEE Trans. Signal Process. 38 (12) (December 1990) 2127–2135. [9] A.S. Kayhan, A. El-Jaroudi, L.F. Chaparro, Evolutionary periodogram for non-stationary signals, IEEE Trans. Signal Process. 42 (6) (June 1994) 1527–1536. [10] A.S. Kayhan, A. El-Jaroudi, L.F. Chaparro, Dataadaptive evolutionary spectral estimation, IEEE Trans. Signal Process. 43 (1) (January 1995) 204–213. [11] S. Mallat, Z. Zhang, Matching pursuit with time– frequency dictionaries, IEEE Trans. Signal Process. 41 (12) (December 1993) 3397–3415. [12] S. Mann, S. Haykin, The chirplet transform: physical considerations, IEEE Trans. Signal Process. 43 (1995) 2745–2761. [13] G. Melard, A.H. Schutter, Contributions to evolutionary spectral theory, J. Time Ser. Anal. 10 (January 1989) 41–63. [14] M.B. Priestley, Non-linear and Non-Stationary Time Series Analysis, Academic Press, London, 1988. [15] S. Qian, D. Chen, Discrete Gabor transform, IEEE Trans. Signal Process. 41 (7) (July 1993) 2429–2439. [16] M.K. Simon et al., Spread Spectrum Communications, Computer Science, New York, 1985. [17] J. Wexler, S. Raz, Discrete Gabor expansions, Signal Processing 21 (3) (November 1990) 207–220.