Extension of the Capon's spectral estimator to time–frequency analysis and to the analysis of polynomial-phase signals

Signal Processing 83 (2003) 575 – 592 www.elsevier.com/locate/sigpro

Extension of the Capon’s spectral estimator to time–frequency analysis and to the analysis of polynomial-phase signals ∗ , Mehmet Tankut Ozgen Department of Electrical and Electronics Engineering, Anadolu University, TR-26470, Eskis!ehir, Turkey Received 14 September 2001; received in revised form 4 August 2002

Abstract Incorporation of the linear time-varying 6lter and its Zadeh’s generalized transfer function concepts to the derivation of the Capon’s minimum-variance spectral estimator leads to a new, bilinear, cross-term suppressed and alias-free time-frequency representation (TFR) that has a higher resolution than the spectrogram with the same window width. Time-variant autocorrelation function of the nonstationary signal of interest is employed in this proposed TFR. By adopting an approximation for time-variant autocorrelation functions, we obtain another new, bilinear, parameterized TFR related to the spectrogram, the frequency resolution of which can be adjusted by varying its parameter for a 6xed window width. We compare resolution and cross-term suppression properties of these proposed TFR’s with other basic bilinear TFR’s, via simulations on synthesized signals. Then, by incorporating polynomial-phase kernel functions to the Capon’s estimator, we propose new bilinear signal representations for the analysis of constant-amplitude polynomial-phase signals and apply them to interference excision in direct-sequence spread-spectrum communications. ? 2002 Elsevier Science B.V. All rights reserved. Keywords: Time–frequency representations; Capon’s minimum-variance spectral estimator; Spectrogram; Wigner distribution; Polynomial-phase signals; Polynomial spectrogram; Direct-sequence spread-spectrum communications; Interference excision; Notch 6ltering

1. Introduction Capon’s minimum-variance spectral estimator is a 6lter bank approach based on a data-adaptive bandpass 6lter centered at the analysis frequency [10,20,26,30]. In general, it has a higher-frequency resolution than nonparametric spectral estimators, attributed to its lower variance [20,26,30]. Hence, a time-frequency version of Capon’s method is expected to have a ∗

Present address: 16. Sokak 18B-4, Bahcelievler, Ankara 06490, Turkey. Tel.: +90-222-335-0580x6463; fax: +90-222-323-9501. E-mail address: [email protected] , (M.T. Ozgen).

better time-frequency resolution than the spectrogram (magnitude squared short-time Fourier transform) [15], which is the time-variant extension of a nonparametric spectral estimator, the periodogram [30]. The evolutionary periodogram (EP) [18] and the data-adaptive evolutionary (DAE) spectral estimator [19] are recently proposed bilinear, time-dependent spectral estimators for nonstationary signals, i.e., time-frequency representations (TFR’s), that are derived analogously to Capon’s method in the stationary signal case. They can both be interpreted as the energy of the output of a time-varying bandpass 6lter centered at the analysis frequency [18,19]. In both approaches, the signal is modelled, at a given frequency,

0165-1684/03/$ - see front matter ? 2002 Elsevier Science B.V. All rights reserved. doi:10.1016/S0165-1684(02)00487-5

576

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

Nomenclature (·)T (·)H

denotes transpose denotes conjugate (Hermitian) transpose of matrices and vectors (·)−1 is the inverse of a matrix (·)∗ is the complex conjugate E{·} expectation operation (n; w) discrete time and frequency variables

as a sinusoid with a time-varying amplitude. This amplitude is represented as an expansion of orthonormal functions over the time interval of interest; and, its minimum mean-square error estimate is computed and used in estimating the time-varying spectrum at that frequency, while minimizing the interference from signal components at other frequencies [18,19]. To derive the EP, interfering signal components at other frequencies are assumed to be uncorrelated or white [18]. This assumption is discarded, in the derivation of the DAE estimator, to obtain a better time-frequency resolution than the EP [19]. In this paper, we follow a diKerent and simpler approach than the above work, in order to extend Capon’s estimator to time-frequency analysis. We directly incorporate the linear time-varying 6lter and its Zadeh’s generalized transfer function concepts to the derivation of the Capon’s minimum-variance estimator. This results in a new, bilinear, cross-term suppressed and alias-free TFR that has a higher time-frequency resolution than both the spectrogram with the same window width and the DAE estimator. This new, high-resolution TFR, which we name as Capon’s TFR, has a diKerent form than the DAE estimator; but, they both reduce to Capon’s method when the signal of interest is stationary. Short-time windowing inherent in the proposed Capon’s TFR suppresses cross-terms in the time direction, and, bandpass 6ltering at the analysis frequency inherent in Capon’s method eliminates cross-terms in the frequency direction. Thus, the proposed Capon’s TFR, although bilinear like the Wigner distribution (WD) [11], has cross-term suppression, whereas the WD suKers from the presence of cross-terms for multicomponent signals.

 = (1 ; 2 ; : : : ; R ) Rx;p (n) Rˆ x;p (n) ap (w) and ap (n; )

phase-coeIcient variable vector time-dependent autocorrelation matrix of the short-time signal vector its estimate kernel function vectors

Our proposed Capon’s TFR diKers from the Capon’s estimator for stationary signals only in that it is based on time-dependent autocorrelation sample estimates of the nonstationary signal of interest. The autocorrelation estimator we use in Capon’s TFR is given by Eq. (59) of [19], which is basically a diKerent way of windowing the lag product of the signal. By adopting an approximation for the above autocorrelation estimate based on the outer product of the signal vector, we develop another new, bilinear, parameterized TFR that can be viewed as the spectrogram passed through a data-adaptive, nonlinear steepness function. Hence, for a 6xed window width in time, the frequency resolution of this new TFR can be improved by adjusting its steepness parameter. Then, by replacing complex exponential kernel functions with higher-order polynomial-phase kernel functions in these two proposed TFR’s based on Capon’s method, we obtain two diKerent, bilinear time/phase-coeIcient representations for the analysis of constant-amplitude polynomial-phase signals. Some of the other tools for the polynomial-phase signal analysis are the maximum-likelihood (ML) estimation algorithm for estimating phase coeIcients [7,23,24,12], the high-order ambiguity function (HAF) [23–25], the polynomial Wigner–Ville distribution (PWVD) [8,9,27], and more general polynomial time-frequency distributions [6]. High-order nonlinear transforms indicated above are also applied to multicomponent signals successfully [9,5,28,34]. Moreover, [5] reports improved robustness to noise; SNR threshold of the product HAF is drawn to as low as 0 dB. The feature that distinguishes our proposed time/phase-coeIcient representations from the above

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

work is that their derivation is based on a notch 6ltering operation in the phase-coeIcient domain. Hence, they are aimed at and suitable for simultaneous detection and excision of polynomial-phase interference components in direct-sequence spread-spectrum (DS-SS) communication systems. The paper is organized as follows. Section 2 derives the proposed Capon’s TFR. Section 3 presents the autocorrelation estimator [19] used in this TFR from a windowing perspective, and also brieMy reviews other estimators for time-dependent correlation functions of nonstationary signals. Then, another parameterized TFR is developed by approximating the used estimator. Section 4 compares proposed TFR’s with the alias-free discrete WD [22], the spectrogram with identical window width, the spectrogram with reassignment [3] and the DAE estimator in terms of resolution and cross-term suppression, via simulation examples. Section 5 derives new bilinear signal representations for the analysis of polynomial-phase signals, based on Capon’s estimator. It highlights notch 6ltering inherent in this derivation and applies them to interference suppression in DS-SS communications. Section 6 presents simulation results on polynomial-phase signal detection and notch 6ltering. Section 7 is the conclusion, which also summarizes the comparisons.

2. Derivation of Capon’s TFR The derivation given below is in parallel with that of the Capon’s estimator for the stationary signal case [20,26,30]. Let x(n); 0 6 n 6 N − 1 be the discrete-time, nonstationary signal of interest, where N is the data length. This signal is passed through a linear, causal, time-varying, bandpass 6lter centered at the analysis frequency w: y(n) =

n


T

h(n; m)x(m) = h (n)x(n)

(1)

m=n−p

for 0 6 n 6 N − 1, where h(n; m) is the impulse response of the 6lter and p is the 6lter order and the integer-valued window-width parameter T

h(n) = [h(n; n); h(n; n − 1); : : : ; h(n; n − p)]

(2)

577

is the impulse response vector and x(n) = [x(n); x(n − 1); : : : ; x(n − p)]T

(3)

is the short-time signal vector, at time n, above. From Eq. (1), the power of the 6lter output, at time n, can be written as E{|y(n)|2 } = hT (n)Rx;p (n)h∗ (n);

(4)

where the expectation above is an ensemble average at time n and, Rx;p (n) = E{x(n)xH (n)}

(5)

is the time-dependent autocorrelation matrix of the signal vector. The time-varying frequency response of the above time-varying 6lter is given by the Zadeh’s generalized transfer function of the 6lter evaluated on the unit circle [18,19] H (n; w) =

n


h(n; m)e−jw(n−m) = hT (n)ap (w);

m=n−p

(6) where ap (w) = [1; e−jw ; : : : ; e−jpw ]T

(7)

is the kernel function vector composed of complex exponentials at the analysis frequency. To make the 6lter as selective as possible for a frequency band around the analysis frequency w, at time n, we minimize the power at the 6lter output in Eq. (4) subject to the constraint that the signal component at the frequency w is passed undistorted, i.e., H (n; w)=1: min hT (n)Rx;p (n)h∗ (n) h(n)

subject to hT (n)ap (w) = 1;

(8)

so that components of x(n) at other frequencies will be suppressed. The solution to (8) can be found, e.g., by means of quadratic minimization [30], as ∗ H −1 hopt (n) = R−1 x;p (n)ap (w)=[ap (w)Rx;p (n)ap (w)]:

(9)

Inserting Eq. (9) into Eq. (4) gives Px (n; w) =

1 apH (w)R−1 x;p (n)ap (w)

(10)

578

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

as the 6lter output power, which can be viewed as the power of x(n) in a passband centered at the frequency w and in a time interval around n. By varying w over the whole radian frequency range − 6 w 6 , and n over the whole time range 0 6 n 6 N − 1, a TFR (or a time-varying spectral estimator) for the signal x(n) is obtained as a direct extension of Capon’s method. Eq. (10), together with Eqs. (3), (5) and (7), constitutes Capon’s TFR, as we name it. Short-time windowing imposed in the above derivation suppresses cross-terms in the time direction, and, bandpass 6ltering tuned to the analysis frequency eliminates cross-terms in the frequency direction. Thus, the proposed Capon’s TFR has cross-term suppression. It is also an alias-free TFR, in case of Nyquist rate sampling, due to the discrete-time nature of its derivation. The feature of Eq. (10) that distinguishes it from the Capon’s estimator for stationary signals, i.e., that provides time variation in Eq. (10) is the use of time-dependent autocorrelation of the nonstationary x(n), Rx;p (n). Its estimation is essential for Capon’s TFR and is discussed in the next section. The autocorrelation estimation is also important in that the singularity of the autocorrelation matrix estimator Rˆ x;p (n) should be avoided in, e.g., time intervals where the signal lies at zero.

3. Autocorrelation estimation and a parameterized version of the spectrogram 3.1. Time-dependent autocorrelation estimation Kayhan et al. [19] propose an estimator for the autocorrelation of a nonstationary signal, based on representing time-varying amplitudes of frequency components of the signal as expansions of orthonormal functions and estimating and substituting these amplitudes in the expectation. This estimator is given by Eq. (59) in [19], in its 6nal form. By substituting the Fourier expansion functions, i (n) = ej(2=N )in ; 0 6 n 6 N − 1, 0 6 i 6 M − 1, into Eq. (59) of [19], we obtain the autocorrelation estimator used in Capon’s TFR. M is the number of the expansion functions, above. The estimate of the autocorrelation rx (n; m) = E{x(n)x∗ (m)} is then

given by rˆx (n; m) =

N −1 N
x(k) x∗ (k + m − n) M k=0

×

sin2 [(=N )M (k − n)] sin2 [(=N )(k − n)]

(11)

for 0 6 n; m 6 N − 1. This is the autocorrelation estimator used when estimating the elements of Rx;p (n) in Eq. (10), in this paper. It is actually a diKerent form of windowing the lag product of the signal, and, M controls the window extent. Hence, there are two parameters of Capon’s TFR given by Eqs. (10) and (11); the short-time window-width parameter p, and, the window-width parameter M involved in autocorrelation estimation. Together, they control the resolution and the robustness of resulting time-frequency patterns. The window in the above estimator is w(k − n) =

sin2 [(=N )M (k − n)] sin2 [(=N )(k − n)]

around the time sample n. For M =1, w(k −n)=1 and Eq. (11) reduces to the standard biased autocorrelation estimate for stationary signals. For M = N , w(k − n) = N 2 (k − n) ((·) denoting the Kronecker delta) and Eq. (11) reduces to the outer product of the signal, making the estimate of Rx;p (n) singular. As M increases and approaches to the data length N , the main lobe of the window w(k) narrows and its side-lobes are eKectively suppressed. Since the entire signal is used in the estimator of Eq. (11), Rˆ x;p (n) is prevented from becoming singular in time intervals where the signal is zero. We brieMy review other methods for the estimation of time-variant correlation functions of nonstationary signals, below. Computing a sliding window average of the lag product of the signal, as done above, and an orthogonal series expansion in time of the lag-product time-series with the lag treated as a parameter are two main approaches [14]. As alternative estimators, we have tried sliding Gaussian window with w(k − n) = exp[ − (k − n)2 =W ] and sliding rectangular window with w(k − n) = rect[(k − n)=W ], W being the window-width parameter, in time-variant autocorrelation estimation needed in Eq. (10). The orthogonal series expansion approach such as that in [21] also reduces to a form

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

579

of windowing the lag product. We have substituted Fourier series in this method and used the resulting window function. As a third approach, we have tried the signal correlation matrix estimator used in the exponentially weighted recursive least-squares (RLS) adaptive algorithm for direct-form FIR 6lters, where signal outer product terms are weighted by an exponential forgetting factor [26]. Simulations on monocomponent and multicomponent constant-amplitude quadratic-phase signals indicate that estimators of Eqs. (11), (12) and that in (14) (the latter two will be presented in the next subsection yield, in general, sharper and more robust linear tracks in time-frequency patterns. Hence, in this paper, simulation results on Capon’s method presented in Section 4 and 6 are obtained using these autocorrelation estimators.

It follows from the Schwarz inequality that the above TFR is positive. It can be viewed as the spectrogram passed through a data-adaptive, nonlinear steepness function. Hence, for a 6xed window width p + 1, the frequency resolution of this new TFR can be improved by adjusting its steepness parameter . As  is decreased towards 0, the nonlinearity becomes steeper resulting in a sharper TFR. This is veri6ed by the argument given below. If  is in the approximate range 0:01x2 —0:1x2 ,

3.2. A parameterized version of the spectrogram

|apH (w)x(n)|2 6 ap (w)2 x(n)2 = (p + 1)x(n)2

The autocorrelation estimator Rˆ x;p (n) resulting from Eq. (11) is dominated by the outer product x(n)xH (n), for M = N − 1. This motivates the following approximation to the time-dependent autocorrelation matrix, as a diKerent way of perturbing the outer product:

with equality if and only if x = Kap where K is a complex-valued constant. Hence,

Rˆ x;p (n) = I + x(n)xH (n);

(12)

where I is the (p + 1) × (p + 1) identity matrix and  is a small, positive parameter. By substituting Eq. (12) into Eq. (10), and using the matrix inversion lemma [17] as given by   xxH 1 H −1 I− [I + xx ] =   + xH x for inverting Rˆ x;p (n), we obtain a new, bilinear TFR related to the spectrogram with rectangular window and parameterized by : Px (n; w; ) =

 ; (13) (p + 1) − [|apH (w)x(n)|2 =( + x(n)2 )]

where x(n) is the short-time signal vector given by Eq. (3),  ·  denotes the Euclidean norm of a vector, ap (w) is the kernel function vector given by Eq. (7), |apH (w)x(n)|2 is the spectrogram with rectangular window, and p + 1 is the width of this rectangular window.

Px (n; w; ) ≈

 (p + 1) − [|apH (w)x(n)|2 =x(n)2 ]

which is actually the spectrogram written in a slightly diKerent form. Suppose that we decrease  further. From Schwarz inequality,

lim Px (n; w; ) = 0

→0

= |K|2

if x(n) = Kap (w) if x(n) = Kap (w):

This shows that as  gets smaller and approaches to 0, the resolution of the TFR given by Eq. (13) improves and Px (n; w; ) becomes binary in the limiting case. As indicated above, the threshold of the nonlinearity is determined by the match x(n) = Kap (w), in the limit. When  is further decreased beyond this highresolution range, time-frequency patterns are lost, unless the signals are stationary, as our simulations have indicated. Since there is a mismatch between the short-time signal vector x(n) and the kernel function vector ap (w) unless the signal of interest is short-time stationary or stationary, the binary nature of Px (n; w; ) in the limiting case fails to track the nonstationarity of the signal. The threshold of the nonlinear steepness function will be above the maximum of the spectrogram to which it is applied, in that case. Hence, the resolution cannot be improved inde6nitely; there is a limit to the improvement. For multicomponent quadratic-phase signals, Eq. (13) provides only a slight resolution improvement over the spectrogram, according to our simulations,

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

580

Table 1 3-dB bandwidths, in radians, of the proposed time-frequency representations given by Eqs. (10) and (13) respectively, and those of the spectrogram with rectangular window of the same width p + 1

Signal model

Proposed TFR’s 

48 p(p+1)(p+2)

x(n) = exp(jw0 n)

W =

x(n) = A exp(jw0 n) + w(n)

WCapon =



PQ x (n; w; ) = U

m=−L

WSpect =

48 p(p+1)(p+2)SNR

since for small  it cannot resolve the nonstationarity of such signals. Hence, to improve the ability of tracking the nonstationarity, we propose an extended version of Eq. (13), by incorporating additional shifted outer product terms to the autocorrelation estimator of Eq. (12). This new TFR is obtained as,

apH (w)[I +

Spectrogram

1 : x(n + m)xH (n + m)]−1 ap (w) (14)

Incorporation of additional outer product terms provides a trade-oK between the resolution and the ability of tracking the nonstationarity of the signal (robustness of patterns), and results in improved resolution for constant-amplitude multicomponent quadratic-phase signals. 3.3. Frequency resolution improvement Derivations of 3-dB bandwidths presented in Table 1 in this subsection in parallel with those presented in [20]. To quantify the resolution improvement of Px (n; w; ) given by Eq. (13) over the spectrogram, we present 3-dB bandwidths W and WSpect of Px (n; w; ) and the spectrogram with rectangular window of the same width p + 1, respectively, in the second row of Table 1 for the indicated single-tone signal. As indicated there, the frequency resolution of Px (n; w; ) improves inde6nitely with square root of , as it approaches to zero, for a sinusoidal signal. Px (n; w; ) has better resolution than the spectrogram

WSpect =

 

24 p(p+2) 24 p(p+2)

+

24 p(p+1)(p+2)SNR

as long as  is smaller than half of the window width, which is satis6ed by the small- range of Eq. (13). 3-dB bandwidths WCapon and WSpect of Capon’s TFR given by Eq. (10) and the spectrogram with rectangular window of the same width p + 1, respectively, are presented in the third row of Table 1 for a noisy sinusoid, where the additive noise signal w(n) is zero-mean, white, circular Gaussian noise with variance 2 . The signal-to-noise ratio is SNR = A2 =2 . Comparison of those 6gures reveals that the proposed Capon’s TFR has better resolution than the spectrogram as long as SNR is greater than 1=(p + 1), the inverse of the window width. Hence, for a given SNR level, choosing the short-time window width greater than the inverse of this level maintains a better frequency resolution for the part of Capon’s TFR. However, if SNR1 time resolution would be lost by that choice. 3.4. Computational cost Capon’s TFR derived in Section 2 is computationally intensive. Its computational cost is determined by the computation of autocorrelation sample estimates using Eq. (11) and by the inversion of the non-Toeplitz (p + 1) × (p + 1) matrix Rˆ x;p (n) in Eq. (10) for 0 6 n 6 N − 1. O(N 3 ) computations are required for the entire time-frequency pattern. The parameterized spectrogram given by Eq. (13) requires O(N 2 log N ) computations for the entire time-frequency plane, when FFT algorithms are used. The extended version of the parameterized spectrogram given by Eq. (14) requires O(N 2 log N ) computations due to the evaluation of the quadratic form in the denominator for diKerent frequencies using FFT

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

algorithms, or O(N (p + 1)3 ) computations due to the inversion of the time-variant autocorrelation matrix estimator for each time sample, whichever is dominant for particular N and p values. 4. Simulations and comparisons with other bilinear TFR’s In the time-frequency patterns presented in this section, the horizontal axis is the discrete-time axis; 0 6 n 6 N − 1. The vertical axis is the discrete radian frequency axis; wm =(2=N )m, −N=2 6 m 6 N=2−1. N is the signal length. In the 6rst example, we consider the signal 2

x(n) = ej(=N )n ;

0 6 n 6 N=2 − 1 2

= ej[=(2N )](n−N=2) ;

N=2 6 n 6 N − 1;

where N = 128. Fig. 1(a) shows the alias-free discrete WD for this signal. Fig. 1(b) shows the spectrogram with rectangular window, with the width parameter p = 5. Fig. 1(c) presents the spectrogram with reassignment, with p = 5 again. Fig. 1(d) presents the DAE spectrum using the Fourier expansion functions with M = 14. M is the number of expansion functions, here. Fig. 1(e) displays Capon’s TFR with the same short-time window-width parameter p = 5 and with the window parameter M = N − 1 in autocorrelation estimation. Fig. 1(f) displays the parameterized version of the spectrogram given by Eq. (13) with p = 5 and the steepness parameter  = 10−6 Ex , where Ex is the signal energy. We compare these 6gures. Capon’s TFR and the TFR of Eq. (13) have almost the same time-frequency resolution, better than those of the spectrogram with the same window and the DAE estimator. The parameter of the DAE estimator is chosen to give the best resolution for it, in this comparison. The WD has cross-terms, whereas Capon’s TFR and the TFR of Eq. (13) do not. The spectrogram with reassignment has slightly better resolution than our proposed TFR’s for this signal, but its performance becomes worse for multicomponent signals. The reassignment operation splits time-frequency tracks for such signals. Moreover, it makes the TFR singular in time-intervals where the signal lies at zero. These issues are illustrated in other examples below.

581

In Example 2, we consider the signal   n2 ; 0 6 n 6 N=2 − 1 x(n) = cos N = 0;

N=2 6 n 6 3N=4 − 1   = cos n ; 3N=4 6 n 6 N − 1; 4 where N = 128 again. Fig. 2(a) displays the alias-free WD for this signal. Fig. 2(b) shows the spectrogram with rectangular window, with p = 5 again. The reassignment stage causes division by zero in the zero-interval of the signal. Fig. 2(c) shows the DAE spectrum using the Fourier expansion functions with M = 14 again. Fig. 2(d) presents Capon’s TFR with parameters p = 5 and M = N=2. Fig. 2(e) presents the TFR given by Eq. (14) with parameters p = 5;  = 10−4 Ex ; U = 0 and L = 1. Small- form of Eq. (13) could not track the nonstationarity of the signal above; hence, both the outer product and its delayed version by one sample are used in the autocorrelation estimator to resolve the nonstationarity. By comparing these 6gures, we see that our proposed TFR’s have better resolution than the spectrogram and the DAE estimator. Their resolution is almost as good as that of the WD for this signal, even better in the third segment of the signal. In Example 3, we compare the spectrogram with reassignment and our proposed TFR’s for x(n) = cos[(=128)n2 ], 0 6 n 6 127. Fig. 3(a) is the spectrogram with rectangular window and with reassignment stage, for this signal. The window width parameter is p = 5. The degradation in the time-frequency pattern is clearly seen in this 6gure, for this multicomponent signal. Fig. 3(b) displays Capon’s TFR with p = 5 and M = N=2. Fig. 3(c) represents the parameterized spectrogram proposed by Eq. (14) in this paper, with parameters p = 5,  = 10−3 Ex and U = L = 2. Linear time-frequency tracks are much more robust in them. We have also computed these bilinear TFR’s for the above signals embedded in additive, zero-mean, white, circular Gaussian noise, and, by increasing the noise variance gradually, we have attempted to visually identify the linear time-frequency tracks of these quadratic-phase signals for each time, until no such identi6cation has been possible. No signi6cant diKerence in regard to robustness to noise has been

582

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

Fig. 1. (a) Alias-free discrete Wigner distribution for the signal in Example 1. (b) Spectrogram with rectangular window. (c) Spectrogram with rectangular window and with the reassignment stage. (d) Data-adaptive evolutionary spectrum. (e) Capon’s time-frequency representation. (f) Parameterized spectrogram given by Eq. (13).

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

583

Fig. 2. (a) Alias-free discrete Wigner distribution for the signal in Example 2. (b) Spectrogram with rectangular window. (c) Data-adaptive evolutionary spectrum. (d) Capon’s time-frequency representation. (e) Parameterized spectrogram given by Eq. (14).

584

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

Fig. 3. (a) Spectrogram with rectangular window and with the reassignment stage for the signal in Example 3. (b) Capon’s time-frequency representation. (c) Parameterized spectrogram given by Eq. (14).

found among these compared TFR’s. Hence, we do not present these noisy time–frequency patterns here.

A multicomponent polynomial-phase signal to be analyzed in this section is given by

5. Capon’s method for the polynomial-phase signal analysis

x(n) =

5.1. Modifying Capon’s method

where each component has a constant amplitude Ak and a phase k (n) modelled as a polynomial function of time of order Rk . Let the above signal be passed through a timevarying 6lter with order p, as given by Eq. (1), in order to select a particular polynomial-phase component in Eq. (15) and suppress the others. If a phase term, exp[j (n)], is given as the input to this 6lter, the output of the 6lter will be this phase term multiplied with

In this subsection, we propose two new bilinear signal representations for the analysis of constant-amplitude polynomial-phase signals, based on the Capon’s estimator. They are obtained by replacing complex exponential kernel functions with higher order polynomial-phase kernel functions in the two TFR’s developed in Sections 2 and 3.

K


A k ej

k (n)

;

0 6 n 6 N − 1;

(15)

k=1

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

the factor

analogous to Eq. (13), as n


H (n; ) =

h(n; m)e−j[

Px (n; ; )

(n)− (m)]

m=n−p T

= h (n)ap (n; );

(16)

where h(n) is the impulse response vector of the 6lter given by Eq. (2) and ap (n; ) = [1; e

−j[ (n)− (n−1)]

;:::;e

−j[ (n)− (n−p)] T

] (17)

is the kernel function vector composed of polynomialphase kernel functions.  = (1 ; 2 ; : : : ; R ) is the coeIcient vector for the polynomial phase (n) =

585

R


r nr ;

(18)

r=0

excluding the zeroth-order coeIcient. Eq. (16), together with Eqs. (17) and (18), can be viewed as the extension of the Zadeh’s generalized transfer function of the 6lter to polynomial-phase signals. The frequency variable w in Eq. (6) is replaced by the phase-coeIcient variable vector  above. The kernel function vector ap (n; ) becomes time dependent, as can be seen from Eqs. (17) and (18). In order to select the polynomial-phase component with the coeIcient vector  of the signal in Eq. (15), and, suppress the other components; we minimize the power at the 6lter output subject to the constraint that the signal component with phase coeIcients  is passed undistorted, i.e., H (n; ) = 1. As the solution to this constrained minimization, analogous to the development given by Eqs. (8)–(10), we obtain the following time/phase-coeIcient representation 1 ; (19) Px (n; ) = −1 H ap (n; )Rx;p (n)ap (n; ) where Rx;p (n) is the time-dependent autocorrelation matrix of the short-time signal vector x(n) in Eq. (3). ap (n; ) is the kernel function vector given by Eqs. (17) and (18). By replacing Rx;p (n) above with the autocorrelation matrix estimator Rˆ x;p (n) obtained by Eq. (11), we obtain our 6rst bilinear representation proposed for the polynomial-phase signal analysis in this paper. By substituting the estimator given by Eq. (12) into Eq. (19), we obtain our second bilinear representation,

=

 : (p + 1) − [|apH (n; )x(n)|2 =( + x(n)2 )] (20)

 is, again, the steepness parameter above. By incorporating additional shifted outer product terms to the autocorrelation estimator of Eq. (12), an extended version of the above representation can be obtained, as given by substituting the kernel function vector ap (n; ) into Eq. (14). Short-time windowing imposed in Eqs. (19) and (20) suppresses cross-terms in the time direction, and, 6ltering tuned to a particular polynomial-phase component with phase coeIcients k eliminates cross-terms in the phase-coeIcient direction. Thus, these two representations have cross-term suppression. As will be illustrated by the simulations in Section 6; the representations given by Eqs. (19) and (20) are not energy or power distributions of the signal in the time/phase-coeIcient space (n; ), since an arbitrary signal cannot be decomposed into polynomial-phase components of order greater than or equal to 2, as it usually can in the case of complex exponentials by means of the Fourier transform. Eqs. (19) and (20) are merely tools for identifying the components of a multicomponent polynomialphase signal as in Eq. (15), when the phase (n) in the kernel function vector matches to any of the phases k (n) of the components. Each polynomial-phase component is displayed as a ridge by these representations, along a line segment parallel to the time axis and located at  = k , in the (n; ) space. 5.2. Capon notch ;lters in the phase-coe
586

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

for simultaneous detection and suppression of interference signals in DS-SS communication systems. Amin [1] uses bilinear TFR’s in detecting instantaneous frequency (IF) tracks of interference components and suppresses those components by a time-varying, open-loop, adaptive notch 6lter with zeros placed at identi6ed IF’s. However, such a notch 6ltering is not successful in suppressing multicomponent interferences [2]. Multicomponent interference excision is achieved successfully by projecting the received signal onto the interference-orthogonal subspace [2,33,4]. However, this operation requires accurate IF estimation for each interference component. When their phase orders are greater than two, i.e., for nonlinear FM laws, IF estimation based on bilinear TFR’s may be severely biased [29]. Hence, other IF estimation methods have been proposed [16,13,29]; but, [13,29] again use FIR notch 6lters tuned to estimated IF for interference suppression, which doesn’t perform well for multicomponent jammers. Ikram et al. [16] uses the HAF for the monocomponent jammer case only. Assuming a parametric interference model given by Eq. (15), Eq. (19) or (20) can be used to estimate interference parameters, and the underlying notch 6lter can be used to excise the identi6ed component simultaneously. This constitutes an alternative model-based interference excision technique. According to the used signal model, after being demodulated and sampled at chip rate, the signal for one bit at the receiver of the DS-SS system becomes x(n) = p(n) + i(n) + w(n);

0 6 n 6 N − 1; (21)

where p(n) is the pseudo-noise (PN) chip sequence with length N for the transmitted bit, i(n) is the interference signal given by Eq. (15) in general, w(n) is the complex additive white Gaussian noise (AWGN) with zero mean and variance 2 both for in-phase and quadrature noise components. The general time/phase-coeIcient representation given by Eq. (19) is the signal power at the output of the notch 6lter obtained by substituting the kernel function vector ap (n; ) into Eq. (9). This notch 6lter passes a polynomial-phase component with phase-coeIcient vector  with unity gain, and suppresses other components. It can be used to suppress detected polynomial-phase components in the received signal x(n) of Eq. (21), as follows.

A monocomponent interference signal i(n) in Eq. (21) is detected as n = arg max Px (n; ) 

= arg max 

1

(22)

apH (n; )R−1 x;p (n)ap (n; )

at each time instant n, where Px (n; ) is Capon’s representation given by Eq. (19) or (20) computed for the received signal x(n). This is because p(n) and w(n) are spread over the entire time/phase-coeIcient domain, whereas i(n) is localized to robust line segments. n is the identi6ed phase-coeIcient vector of the polynomial-phase component i(n) at time n. i(n) is simultaneously excised by the following notch 6ltering operation: y(n) = x(n) −

apH (n; n )R−1 x;p (n)x(n)

apH (n; n )R−1 x;p (n)ap (n; n )

;

(23)

where x(n) is the received signal sample at time n, x(n) is the short-time received signal vector as in Eq. (3), and y(n) is the notch 6lter output before despreading in the DS-SS receiver. If the interference signal component i(n) in the received signal x(n) of Eq. (21) is multicomponent as in Eq. (15), detected peaks of Px (n; ) at each time instant n reveal phase-coeIcient vectors of polynomial-phase components for that time instant: n; 1 ; : : : ; n; k ; : : : ; n; K . K is the number of detected components at time n. Derivation of the notch 6lter that excises these multiple components is in parallel with that of Capon’s method. We construct the 6lter that minimizes its output power subject to the constraint that identi6ed polynomial-phase components are passed with unity gain, by using quadratic minimization [30]. Then, selected components are subtracted from the received signal to excise them out. The resulting notch 6lter is −1 −1 y(n) = x(n) − 11×K [AH n Rx;p (n)An ] −1 ×AH n Rx;p (n)x(n);

(24)

where 11×K is an all-1 row vector of size K and An = [ap (n; n; 1 ); : : : ; ap (n; n; k ); : : : ; ap (n; n; K )] is the kernel matrix whose columns are kernel function vectors evaluated for K phase-coeIcient vectors identi6ed from peaks of Px (n; ) at time n.

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

Eqs. (23) and (24) are Capon notch 6lters for monocomponent and multicomponent cases, respectively. If phase-coeIcient variable vector  is replaced by the frequency variable w, they reduce to time-varying, open-loop, adaptive notch 6lters in the frequency domain. By embedding peak detection and notch 6ltering operations given by Eqs. (22)–(24) into the computation of Capon’s representation, we obtain a new model-based interference excision technique for DS-SS systems. It will be illustrated by simulations in Section 6. 5.3. Polynomial spectrogram and the ML estimation algorithm In this subsection, we mention the polynomial spectrogram since it is closely related to the representations proposed in Section 5.1 and given by Eqs. (19) and (20). Eq. (20) can be viewed as a variant of it. These three transforms are compared in the simulations presented in Section 6. The polynomial spectrogram used in the ML estimation algorithm for estimating phase coeIcients is [7,23,24]: 2  N −1  
 −j (k)  Sx (n; ) =  x(k)w(k − n)e (25)    k=0

for 0 6 n 6 N − 1, where the phase (k) is given by Eq. (18), =(1 ; 2 ; : : : ; R ) is the coeIcient vector of the phase, and w(·) is the short-time window centered at time sample n. The ML estimates of the phase coeIcients of a constant-amplitude monocomponent polynomial-phase signal in additive, white Gaussian noise are those coef6cients that maximize Eq. (25) over the whole parameter space, as shown in [7]. For a constant-amplitude multicomponent signal as in Eq. (15), computing Eq. (25) and identifying its maxima as the signal components constitute the approximate ML algorithm derived in [24] for estimating the phase coeIcients. 5.4. Selection of the signal length Eqs. (19), (20) and (25) in this section employ polynomial-phase kernel functions with phase orders greater than or equal to 2. Hence, their ability of detecting and identifying polynomial-phase components

587

of a signal as in Eq. (15) correctly depends on the quasi-orthogonality property of such kernel functions. This property is strengthened when the signal length N is a prime integer, as proved for quadratic-phase chirp signals in [32]. We are unable to prove a similar fact for higher-order polynomial-phase functions; but, our simulations con6rm that such a quasi-orthogonality property is also strengthened for them when the signal length N is a prime. Therefore, N should be selected as a prime number for optimal performance. Xia [32] also shows that the number of quadraticphase chirp √ components detectable by Eq. (25) is limited to O( N =2) for a 6xed prime N , when the signal is noiseless. Our simulations con6rm that there are also upper bounds for the number of higher order polynomial-phase signal components detectable by Eqs. (19), (20) and (25), for a 6xed signal length N . As the number of signal components increases, local maxima that identify the signal components become less dominant and the performance rapidly deteriorates. This will be illustrated in the next section. 6. Simulations for polynomial-phase signals In the 6gures of this section, the discrete-time axis spans the entire signal duration, 0 6 n 6 N − 1. The discrete phase-coeIcient axis spans the range, m = (=N )m, −N 6 m 6 N − 1, in radians. N is the signal length. In Example 4, we consider the following monocomponent signal with no noise 4

x(n) = ej(10=127)n ;

0 6 n 6 63 4

= e−j(25=127)n ;

64 6 n 6 126:

Fig. 4(a) shows the polynomial spectrogram of Eq. (25) for this signal, with the window-width parameter p = 10. Fig. 4(b) displays Capon’s time/phase-coeIcient representation given by Eq. (19), with parameters p = 10 and M = N − 1 = 126. Fig. 4(c) shows the parameterized version of the polynomial spectrogram (or the second version of Capon’s representation) proposed as Eq. (20) by us, with parameters p = 10 and  = 10−6 Ex . Kernel function phases are k (n) = (=N )kn4 ; −N 6 k 6 N − 1 for all three representations.

588

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

Fig. 4. (a) Polynomial spectrogram given by Eq. (25) for the signal in Example 4. (b) Capon’s time/chirp-rate representation given by Eq. (19). (c) Parameterized version of the polynomial spectrogram, proposed by Eq. (20).

These 6gures show that proposed representations given by Eqs. (19) and (20) represent the monocomponent signal as clearer line segments than the polynomial spectrogram. Example 5 is related to the following two-component, noiseless signal; x(n) = sin[(25=N )n4 ], 0 6 n 6 N − 1 with N = 127. Fig. 5(a) presents the polynomial spectrogram for this signal. Figs. 5(b) and 5(c) show proposed time/phase-coeIcient representations given by Eqs. (19) and (20), respectively. The parameter values and the kernel functions are the same with those of Example 4. Performance degradation of all three representations are evident in this two-component signal case, as predicted in Section 5.4. Horizontal lines that represent signal components are not as clear as those in Fig. 4.

In Example 6, we consider a signal modelled by Eq. (21). x(n) can be viewed as a received signal for one bit, in a DS-SS communication system, with 64 chips per bit and sampled at a rate of 2 samples/chip. Signal length is N = 128 samples. p(n) is the PN chip sequence representing the bit with 64 ± 1 chips. i(n) = 10 cos[(30=N )n3 ] is the interference component. w(n) is real-valued, zero-mean AWGN with standard deviation  = 0:25. Two-component interference signal is detected as peaks exhibited by Capon’s representation proposed as Eq. (20), and is simultaneously excised by the Capon notch 6lter given by Eq. (24) and embedded in Eq. (20). Figs. 6(a) and 6(b) show the PN sequence to be recovered and the received signal for the corresponding bit, respectively. Fig. 6(c) gives the output of the Capon notch 6lter with parameters p = 20 and

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

589

Fig. 5. (a) Polynomial spectrogram given by Eq. (25) for the signal in Example 5. (b) Capon’s time/chirp-rate representation given by Eq. (19). (c) Parameterized version of the polynomial spectrogram, proposed by Eq. (20).

 = 10−1 Ex . Autocorrelation matrix estimator used in Eq. (24) is that of Eq. (12). By passing this output through a hardlimiter with zero threshold, we obtain Fig. 6(d) as the recovered chip sequence. Comparison with Fig. 6(a) reveals that the original PN chip sequence is recovered with slight initial degradation. 7. Conclusion For the simulations presented in Section 4, proposed eKectively cross-term-free TFR’s based on Capon’s method have better frequency resolution than the spectrogram with the same window width, as quanti6ed by 3-dB bandwidth 6gures presented in Table 1. Simulations in Section 4 also show their superior resolution to the DAE estimator [19].

In the second part of the paper, we extend Capon’s method to develop two diKerent, bilinear signal transforms for the analysis of constant-amplitude polynomial-phase signals. They are given by Eqs. (19) and (20), and should be viewed as alternatives to the polynomial spectrogram of Eq. (25) used in the ML estimation algorithm for estimating phase coeIcients [7,23,24]. Eq. (20) is actually a variant of the polynomial spectrogram and gives better results than it for noiseless monocomponent signals, and similar results to it otherwise, as demonstrated by the simulations in Section 6. Interference excision technique for DS-SS communication systems proposed in Section 2 is only an alternative to other existing model-based approaches [2,33,4], and does not replace them. Polynomial

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

590 2

10

RECEIVED SIGNAL

PN CHIP SEQUENCE

1.5 1 0.5 0 −0.5 −1 −1.5 −2

5

0

−5 −10

0

50

(a)

100

0

50

(b)

TIME

100 TIME

RECOVERED PN SEQUENCE

NOTCH FILTER OUTPUT

2 2 1 0 −1 −2 −3 0

(c)

50

1 0.5 0 −0.5 −1 −1.5 −2 0

100 TIME

1.5

(d)

50

100 TIME

Fig. 6. (a) A PN chip sequence representing a bit in a DS-SS communication system. (b) Received signal for that bit according to the signal model in Example 6. (c) Output of the Capon notch 6lter of Eq. (24) used in conjunction with the representation of Eq. (20). (d) Recovered chip sequence by hardlimiting the sequence in (c).

spectrogram or generalized Wigner-Hough transform as in [4] are other model-based tools for identifying interference components. Once they are identi6ed by these or other methods mentioned in Section 5.2, projection 6ltering techniques [2,33,4] oKer excellent alternatives for their excision. Our simulations indicate that proposed Capon notch 6lters of Eqs. (23) and (24) give similar results to short-time versions of those projection 6lters in [2,4]. Suleesathira and Chaparro [31] reports another interference excision technique that does not require parametric models for interference IF’s. IF estimation is based on recursive correction of a piecewise linear estimate given by the evolutionary-Hough transform

initially. Then, the synthesized interference signal is subtracted from the received signal. Time/phase-coeIcient space with a large dimension results in high computational cost and diIcult search procedure for maxima, hence, limits practical use of the proposed polynomial-phase signal analysis and interference excision techniques. Therefore, phase orders should be kept small and the window length should be chosen so that the signal’s IF law could be approximated by polynomials of chosen order within the window, as suggested in [7]. While using Eqs. (19), (20) and (25) in polynomialphase signal identi6cation, quasi-orthogonality of polynomial-phase kernels may cause spurious peaks

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

in them that may be misinterpreted as false signal components. Hence, exhibited peaks above appropriately chosen thresholds should be taken as signal components only. Threshold levels are determined by the polynomial order and the signal length.

[14]

References

[16]

[1] M.G. Amin, Interference mitigation in spread spectrum communication systems using time-frequency distributions, IEEE Trans. Signal Process. 45 (1) (January 1997) 90–101. [2] M.G. Amin, G.R. Mandapati, Nonstationary interference excision in spread spectrum communications using projection 6ltering methods, Proceedings of the 32nd Annual Asilomar Conference on Signals, Systems, Computers, Paci6c Grove, CA, November 1998, pp. 827–831. [3] F. Auger, P. Flandrin, Improving the readability of time-frequency and time-scale representations by the reassignment method, IEEE Trans. Signal Process. 43 (5) (May 1995) 1068–1089. [4] S. Barbarossa, A. Scaglione, Adaptive time-varying cancellation of wideband interferences in spread-spectrum communications based on time-frequency distributions, IEEE Trans. Signal Process. 47 (4) (April 1999) 957–965. [5] S. Barbarossa, A. Scaglione, G.B. Giannakis, Product high-order ambiguity function for multicomponent polynomial-phase signal modeling, IEEE Trans. Signal Process. 46 (3) (March 1998) 691–708. [6] M. Benidir, A. Ouldali, Polynomial phase signal analysis based on the polynomial derivatives decompositions, IEEE Trans. Signal Process. 47 (7) (July 1999) 1954–1965. [7] B. Boashash, Estimating and interpreting the instantaneous frequency of a signal—Part 2: algorithms and applications, Proc. IEEE 80 (4) (April 1992) 540–568. [8] B. Boashash, P. O’Shea, Polynomial Wigner–Ville distributions and their relationship to time-varying higher order spectra, IEEE Trans. Signal Process. 42 (1) (January 1994) 216–220. [9] B. Boashash, B. Ristic, Polynomial time-frequency distributions and time-varying higher order spectra: application to the analysis of multicomponent FM signals and to the treatment of multiplicative noise, Signal Process. 67 (1) (January 1998) 1–23. [10] J. Capon, High-resolution frequency-wavenumber spectrum analysis, Proc. IEEE 57 (8) (August 1969) 1408–1418. [11] T.A.C.M. Claasen, W.F.G. Mecklenbrauker, The Wigner distribution—a tool for time-frequency signal analysis. Part 1: continuous-time signals, Part 2: discrete-time signals, Philips J. Res. 35 (4) (1980) 217–250, 276 –300. [12] B. Friedlander, J.M. Francos, Estimation of amplitude and phase parameters of multicomponent signals, IEEE Trans. Signal Process. 43 (4) (April 1995) 917–926. [13] H.W. Fung, A.C. Kot, K.H. Li, A novel approach to FM interference suppression in DS spread-spectrum communication systems, Proceedings of the IEEE

[15]

[17] [18] [19] [20] [21] [22]

[23] [24] [25] [26] [27]

[28]

[29]

591

International Conference on Acoustics, Speech, Signal Processing, June 2000, pp. 2829 –2832. W.A. Gardner, Correlation estimation and time-series modeling for nonstationary processes, Signal Process. 15 (1) (January 1988) 31–41. F. Hlawatsch, G.F. Boudreaux-Bartels, Linear and quadratic time-frequency signal representations, IEEE Signal Process. Mag. 9 (4) (April 1992) 21–67. M.Z. Ikram, A. Belouchrani, K. Abed-Meraim, D. Gesbert, Parametric estimation and suppression of nonstationary interference in DS-spread spectrum communications, Proceedings of the 32nd Annual Asilomar Conference on Signals, Systems, Computers, Paci6c Grove, CA, November 1998, pp. 1401–1405. T. Kailath, A.H. Sayed, B. Hassibi, Linear Estimation, Prentice-Hall, Upper Saddle River, NJ, 2000, Appendix A, 729pp. A.S. Kayhan, A. El-Jaroudi, L.F. Chaparro, Evolutionary periodogram for nonstationary signals, IEEE Trans. Signal Process. 42 (6) (June 1994) 1527–1536. A.S. Kayhan, A. El-Jaroudi, L.F. Chaparro, Data-adaptive evolutionary spectral estimation, IEEE Trans. Signal Process. 43 (1) (January 1995) 204–213. R.T. Lacoss, Data adaptive spectral analysis methods, Geophysics 36 (4) (August 1971) 661–675. V.Z. Marmarelis, Practical estimation of correlation functions of nonstationary Gaussian processes, IEEE Trans. Inform. Theory 29 (1983) 937–938. A.H. Nuttall, Alias-free Wigner distribution function and complex ambiguity function for discrete-time samples, NUSC Technical Report 8533, Naval Underwater Systems Center, April 1989. S. Peleg, B. Friedlander, The discrete polynomial-phase transform, IEEE Trans. Signal Process. 43 (8) (August 1995) 1901–1914. S. Peleg, B. Friedlander, Multicomponent signal analysis using the polynomial-phase transform, IEEE Trans. Aerosp. Electron. Syst. 32 (1) (January 1996) 378–386. S. Peleg, B. Porat, Estimation and classi6cation of polynomial-phase signals, IEEE Trans. Inform. Theory 37 (2) (March 1991) 422–430. J.G. Proakis, C.M. Rader, F. Ling, C.L. Nikias, Advanced Digital Signal Processing, Macmillan Publishing Company, NY, 1992 (Chapter 8), pp. 521–525, Chapter 6, pp. 352–353. B. Ristic, B. Boashash, Instantaneous frequency estimation of quadratic and cubic FM signals using the cross polynomial Wigner–Ville distribution, IEEE Trans. Signal Process. 44 (6) (June 1996) 1549–1553. S. Shamsunder, G.B. Giannakis, B. Friedlander, Estimating random amplitude polynomial phase signals: a cyclostationary approach, IEEE Trans. Signal Process. 43 (2) (February 1995) 493–505. P. Shan, A.A. Beex, FM interference suppression in spread spectrum communications using time-varying autoregressive model based instantaneous frequency estimation, Proceedings of the IEEE International Conference on Acoustics, Speech, Signal Processing, March 1999, pp. 2559 –2562.

592

( M.T. Ozgen / Signal Processing 83 (2003) 575 – 592

[30] P. Stoica, R.L. Moses, Introduction to Spectral Analysis, Prentice-Hall, Upper Saddle River, NJ, 1997, Chapter 5, pp. 196 –206, Chapter 2, pp. 23–27, Appendix A, pp. 283–284. [31] R. Suleesathira, L.F. Chaparro, Interference mitigation in spread spectrum using discrete evolutionary and Hough transforms, Proceedings of the IEEE International Conference on Acoustics, Speech, Signal Processing, June 2000, pp. 2821–2824. [32] X.-G. Xia, Discrete chirp-Fourier transform and its application to chirp rate estimation, IEEE Trans. Signal Process. 48 (11) (November 2000) 3122–3133.

[33] Y. Zhang, M.G. Amin, Array processing for nonstationary interference suppression in DS/SS communications using subspace projection techniques, IEEE Trans. Signal Process. 49 (12) (December 2001) 3005–3014. [34] G.T. Zhou, Y. Wang, Exploring lag diversity in the high-order ambiguity function for polynomial phase signals, IEEE Signal Process. Lett. 4 (8) (August 1997) 240–242.