Robust spatial time-frequency distribution matrix estimation with application to direction-of-arrival estimation

Signal Processing 91 (2011) 2630–2638

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Robust spatial time-frequency distribution matrix estimation with application to direction-of-arrival estimation W. Sharif , Y. Chakhchoukh, A.M. Zoubir Signal Processing Group, Technische Universit¨ at Darmstadt, Merckstr. 25, Darmstadt 64283, Germany

a r t i c l e in f o

abstract

Article history: Received 12 February 2011 Received in revised form 20 May 2011 Accepted 23 May 2011 Available online 12 June 2011

We consider the problem of estimating spatial time-frequency distribution (STFD) matrices in the presence of impulsive noise. STFD matrices are widely used in sensor array processing for direction-of-arrival estimation and blind source separation of nonstationary sources. Conventional methods fail when the noise is non-Gaussian or impulsive. We propose robust techniques for STFD estimation which are based on pre-processing, robust position based estimation and robust covariance based estimation. The proposed methods are compared in terms of direction-of-arrival estimation performance. & 2011 Elsevier B.V. All rights reserved.

Keywords: Time-frequency analysis Direction-of-arrival estimation Robust estimation Impulsive noise

1. Introduction Direction-of-arrival (DOA) estimation is important for source localization and interference suppression in many real-world applications such as: radar, sonar and mobile communications [9]. In this context, we encounter sources which are generally non-stationary, e.g., frequency modulated (FM) signals. To exploit this non-stationarity, the spatial time-frequency distribution (STFD) introduced in [1], constitutes a common and powerful tool. Indeed, the STFD enhances the spatial resolution and estimation accuracy [1–3,5] by capturing well the different source signatures present in the time-frequency (TF) plane. Using techniques such as beamforming or MUSIC [10] on STFD matrices provide the final DOA estimates. Existing techniques for STFD based DOA estimation consider only the sensors’ thermal noise, which is assumed to be the Gaussian. However, in practical radio environments, noise is often non-Gaussian and impulsive due to  Corresponding author.

E-mail addresses: [email protected] (W. Sharif), [email protected] (Y. Chakhchoukh), [email protected] (A.M. Zoubir). 0165-1684/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2011.05.022

many natural and man-made phenomena [14,15]. Hence, the estimation is performed in a contaminated environment where the observations contain outliers and deviate from the Gaussian assumption. In [12,13,22], it has been shown that standard timefrequency distributions (TFDs) are severely degraded in the presence of impulsive noise. Since the STFD matrices are obtained by computing TFDs across the sensor array, they suffer from a lack of robustness as well. Consequently, the performance of DOA estimation based on standard STFD matrices degrades significantly. Our goal consists of proposing DOA estimation methods based on robust STFD estimates. A few techniques for robust TFD estimation have been suggested in the literature, e.g., in [12,13]. In [12], Djeddi et al. proposed a fractional-lower-order-moment based TFD. In [13], Barkat et al. proposed a median-based TFD estimation. However, their applicability for STFD based DOA estimation has not been investigated so far. This aspect is also treated here. Furthermore, these basic techniques do not offer a good theoretical trade-off between efficiency and robustness [17]. We propose new methods for STFD estimation which follow performant statistical robust procedures [17,18].

W. Sharif et al. / Signal Processing 91 (2011) 2630–2638

The proposed robust STFD estimation techniques are classified as follows:

 pre-processing based STFD estimation,  robust position based STFD estimation, and  non-iterative robust estimation techniques, e.g., robust covariance based estimation. These robust STFDs are then used to design robust DOA estimation methods. Finally, a comparison analysis of the different approaches in view of DOA estimation quality is performed. This paper is organized as follows. In Section 2, we introduce the signal model. A brief definition of the problem of STFD based DOA estimation is given in Section 3. In Section 4, robust STFD estimation based on preprocessing is discussed. In Sections 5 and 6, we describe robust STFD estimation based on robust position and noniterative robust techniques, respectively. Finally, simulation examples are presented in Section 7 and conclusions are drawn in Section 8. 2. Signal model In narrowband array processing, the baseband signal model for K signals, impinging on an m-element sensor array is given by xðtÞ ¼ AsðtÞ þ nðtÞ,

t 2 f1,2, . . . ,Ng

ð1Þ

Here, A is the m  K mixing matrix, x(t) is the m  1 array output vector and sðtÞ is the K  1 source signal vector at time t. The mixing matrix takes the form A ¼ ½aðy1 Þ, . . . , aðyK ÞT , where aðyk Þ, k 2 f1,Kg is the kth array steering vector and y1 , . . . , yK are the source directions of arrival (DOAs). In (1), nðtÞ ¼ ½n1 ðtÞ,n2 ðtÞ, . . . ,nm ðtÞT is the additive noise vector at time t. The noise is assumed to be impulsive. Moreover, we assume that the noise is spatially and temporally independently and identically distributed (i.i.d). The source signals we considered are constant amplitude FM signals of the form sk ðtÞ ¼ ej2pFk ðtÞ ,

k ¼ 1, . . . ,K

and their corresponding instantaneous frequency (IF) is given by fk ðtÞ ¼ dFk ðtÞ=dt. We assume the IFs of the sources to be known. This does not limit our work since the IF can be independently estimated and numerous methods exist for its estimation [4,20,5], also in impulsive noise environments, e.g., in [13,6]. The source waveforms are assumed to lie within the same baseband frequency range, and have been demodulated using a common carrier frequency. The objective is to estimate the DOAs yk ,k 2 f1,Kg, given N snapshots of xðtÞ and exploiting the knowledge of the instantaneous frequencies of the sources, IFk ,k 2 f1,Kg. 3. STFD matrices and direction finding In this section, we provide a brief introduction to STFD based DOA estimation, where we use the spatial TFD matrix instead of the spatial covariance matrix as in conventional direction finding. The spatial TFD matrix is defined in terms

2631

of auto- and cross-TFD of the sensor signals given by ½Dxx ðt,f Þpq ¼ Dxp xq ðt,f ; jÞ p,q 2 f1,mg,

ð2Þ

ðt,f Þ 2 R

where Dxp xq ðt,f ; jÞ is assumed to be a bilinear TFD of Cohen’s class, for which the kernel function is j. In (2), xp and xq for p,q 2 f1,mg denote the pth and qth sensor signal, respectively. We use a discrete-time form of Cohen’s class of TFDs as can be found in [8] Dxp xq ðt,f , jÞ ¼

1 X

jðm,lÞxp ðt þm þ lÞxnq ðt þmlÞej4pfl

l,m ¼ 1

ð3Þ where jðm,lÞ is the kernel function in the time-lag domain and xn denotes the complex conjugate of x. In the sequel, we use the Wigner–Ville distribution (WVD) for which jðm,lÞ ¼ dðmÞ. By averaging the STFD matrix over a subset of signal IF signatures, the source direction can be estimated using, e.g., MUSIC [10] or beamforming techniques. Let ðtn ,fn Þ 2 S k denote the TF points belonging to the kth source. Then the averaged STFD matrix for the kth source is computed as Dk ¼

X 1 Dxx ðtn ,fn Þ #S k ðt ,f Þ2S n n

ð4Þ

k

where #S k is the number of TF points in S k . We estimate the direction of the source signal represented by the IF segment S k by applying the MUSIC algorithm [10] to the averaged STFD matrix Dk . The averaged STFD matrix provides an effective improvement in SNR by amplification of the source eigenvalues with respect to the noise eigenvalues [3]. It can be seen from (3) that an outlier in the observations xp(t) or xq(t) can have an unbounded influence on the estimate of the ðp,qÞth element of the STFD matrix. Hence, the averaged STFD matrix Dk is not a robust estimate and will lead to wrong estimates of DOAs in the presence of even a small fraction of outlying observations xðtÞ. Therefore, we seek a robust estimate of the STFD matrix. We propose two different approaches for robust estimation of STFD matrices. The first approach is based on suppressing the outlying observations before applying standard STFD estimation, i.e., pre-processing. The second approach is to estimate the STFD matrix robustly. For the problem of direction finding, these two schemes are illustrated in Fig. 1. In the next sections, we will describe the pre-processing and the robust STFD estimation steps and will demonstrate their application to direction finding. 4. Pre-processing based robust estimation There are pre-processing techniques that have been widely used to suppress impulsive interference in many different applications and in particular in conventional direction finding techniques, e.g., in [16]. The objective of pre-processing techniques is to suppress impulsive noise from the observations. The pre-processing based methods can be classified into two different types, i.e., transformation based and diagnostic based. Transformation based

2632

W. Sharif et al. / Signal Processing 91 (2011) 2630–2638

Fig. 1. Proposed schemes for STFD based DOA estimation in impulsive noise.

pre-processing techniques use impulsive resistant transformations on the observations. While diagnostic based pre-processing techniques require detection of the presence of an impulsive contaminated sample and then trim that observation. In this section, we propose some preprocessing techniques which can be used for robust estimation of STFD matrices.

(1) Compute a robust estimate of scale for Rðxp Þ and I ðxp Þ, e.g. MADN. Let s p and s~ p denote the robust scale estimates of Rðxp Þ and I ðxp Þ, respectively (2) Record the indices of Rðxp Þ and I ðxp Þ for which jRðxp Þj4 3s p and jI ðxp Þj 4 3s~ p . Let v p and v~ p denote the indices (3) Clip the observations in xp to 0 that are indexed by vp ¼ v p [ v~ p (4) Repeat the above steps for all p 2 f1,mg (5) Apply classical STFD matrix estimation on the new observations.

4.1. Spatial sign function This transformation based pre-processing technique has been used for direction finding in impulsive noise. It was first proposed by Visuri et al. [16]. The transformed observations after applying the spatial sign function are given by 8 < xðtÞ if xðtÞa0 SfxðtÞg ¼ JxðtÞJ ð5Þ : 0 otherwise pffiffiffiffiffiffiffiffiffiffiffiffi where JxJ ¼ ðxH xÞ. Then the STFD estimation procedure is applied on SfxðtÞg instead of xðtÞ. 4.2. 3-s rejection pre-processing

4.3. Normalization technique This is another transformation based technique where normalization is applied to the observed data to suppress the influence of outliers. Each sensor observation is divided by its magnitude. This technique is different from that in Section 4.1, where all the observations at time t, xðtÞ 2 Cm are normalized by their norm. Herein, each sensor observation is individually normalized by its magnitude. The observations after applying the normalization are given by zp ðtÞ ¼

This diagnostic based pre-processing technique detects and trims the contaminated samples. The detection of a contaminated sample is done by comparing the magnitude of the observation with 3s, where s is the scale estimate. If the sample magnitude is greater than 3s then the sample is considered as impulsive and is discarded. We use the normalized median absolute deviation (MADN) [18] as a robust estimate of scale. MADN is defined as the normalized median of absolute deviations from the median. The normalization is done to obtain a consistent estimate of standard deviation under the Gaussian data. Given the observations y ¼ fy1 ,y2 , . . . ,yn g of a random variable, the MADN is given by MedjyMedðyÞj MADNfyg ¼ 0:6745

Table 1 3-s rejection pre-processing.

ð6Þ

Let Rðxp Þ and I ðxp Þ denote the vectors of real and imaginary parts of the pth sensor observations xp 2 CN , respectively, then the 3-s rejection pre-processing is performed as in Table 1.

xp ðtÞ jxp ðtÞj

ð7Þ

where xp(t) denotes the pth sensor signal at time t. Once having the normalized observations, standard techniques for STFD estimation can be applied to get the direction estimates. 5. Robust position based STFD estimation In this section, we describe robust STFD matrix estimation based on robust estimation of position. Let us consider only the ðp,qÞth element of the STFD matrix at a given (t,f) TF point. In the case of a pseudo-WVD (PWVD), we have ^ x x ðt,f Þ ¼ D p q

L=2 X

xp ðt þlÞxnq ðtlÞej4pfl

ð8Þ

l ¼ L=2

where xp ðt þlÞ denotes the pth sensor signal at time t þ l. Let us denote the summand xp ðt þ lÞxnq ðtlÞej4pfl by Gpq ðt,f ,lÞ. Furthermore, remark that we get fGpq ðt,f ,lÞ, l 2 fL=2,L=2gg 2 CL þ 1 , where L þ1 is the window length.

W. Sharif et al. / Signal Processing 91 (2011) 2630–2638

5.1. Median-based STFD estimation Robust position based TFD estimation described in [13] is based on the reformulation of TFD estimation as an estimation of position of a complex vector fGpq ðt,f ,lÞ l 2 ½L=2,L=2g of length L þ1. It can also be defined as an optimization problem of the following type. ^ x x ðt,f Þ ¼ argminQpq ðt,f ,l,Dpq Þ D p q Dpq

ð9Þ

where Qpq ¼

L=2 X

rðepq ðt,f ,l,Dpq ÞÞ

ð10Þ

l ¼ L=2

and epq ðt,f ,lÞ ¼ Gpq ðt,f ,lÞDpq is the error function depending on ðt,f ,lÞ and where Dpq 2 C denotes the quantity to be estimated. The standard TFD is obtained by taking rðepq ðt,f ,l,Dpq ÞÞ ¼ jepq ðt,f ,l,Dpq Þj2 . The resulting TFD is the mean of the vector fGpq ðt,f ,lÞ, l 2 ½L=2, . . . ,L=2g. The estimators based on jepq ðt,f ,lÞj2 are highly sensitive to outliers and yield biased estimates even in the presence of a small fraction of outliers. To get a robust estimate, other cost functions are required. Barkat et al. [13] and Sahmoudi et al. [19], proposed to use the L1-norm based cost function, i.e., rðeðÞÞ ¼ jeðÞj. The estimator corresponding to minimizing an L1-norm is the median which is highly robust against outliers. Optimizing the L1-norm based cost function requires to solve @Qpq =@Dnpq ¼ 0. It leads to non-linear equations and Dpq is obtained iteratively [7]. As Dpq 2 C, the iterative procedure requires complex operations. In [7], it is proposed to estimate the real and imaginary parts of Dpq separately. This separate estimation requires independence of the real and imaginary parts of Gpq ðt,f ,lÞ which does not necessarily hold. In order to overcome this problem and also to provide more efficient estimators for STFD, we generalize the problem by proposing multidimensional position M-estimation. 5.2. M-estimation based STFD In this section, we take into account the existing correlation between the real and imaginary parts of Gpq ðt,f ,lÞ and treat the complex numbers in the sample fGpq ðt,f ,lÞ, l 2 fL=2, . . . ,L=2gg as two-dimensional vectors Gpq ðt,f ,lÞ ¼ ½R ðGpq ðt,f ,lÞÞ,I ðGpq ðt,f ,lÞÞT , where ½T is the transpose. Let ^ x x ðt,f ÞÞ,I ðD ^ x x ðt,f ÞÞT denotes the vector d^ x x ðt,f Þ ¼ ½RðD

2633

rðÞ function. The most common ones are Huber’s function and Tukey’s biweight function [17]. Here, we choose Huber’s r-function defined by ( d2 =2 ifjdj rc rHuber ðdÞ ¼ ð13Þ cjdjc2 =2 ifjdj 4c where c 2 R þ is an appropriately chosen constant and acts as a threshold that ensures larger weights are applied to the normal observations than the outliers. It also provides a trade-off between robustness and efficiency [17]. It can be either chosen empirically or computed, e.g., based on the distribution of Mahalanobis’ distances. The majority of Mahalanobis’ distance follows a w2p distribution, where p denotes the degree of freedom and is equal to the dimension of the observation [23]. Herein, we have p¼ 2 and the qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi threshold c can be set to w2p,0:975 , where w2p,0:975 denotes the 97:5% quantile of w2p distribution. For the M-estimation based STFD, the goal is to estimate both dpq ðt,f Þ and R.

5.3. S-estimation based STFD Consider again the classical cost function given by L=2 X

Qpq ¼

dðt,f ,l,dpq Þ2

ð14Þ

l ¼ L=2

where STFD estimation can be seen as finding the dpq that minimizes the classical estimate of variance of samples fdðt,f ,l,dpq Þ, l ¼ L=2, . . . ,L=2g. This means that the STFD minimizes the classical estimator of scale. Outliers highly influence this scale estimate. A way to robustify the STFD is to minimize a robust scale instead of the classical scale. The chosen robust scale is an M-estimate of scale [17]. Therefore, the objective of S-estimation based STFD is to find dpq and R which minimize the robust M-scale estimator s^ M of distances fdðt,f ,l,dpq ðt,f ÞÞ, l 2 fL=2, L=2gg. This is given by min s^ M s:t:

  L=2 X dðt,f ,l,dpq Þ 1 r ¼ b0 Lþ 1 l ¼ L=2 s^ M

ð15Þ

ð12Þ

where dðt,f ,l,dpq ðt,f ÞÞ is defined in Eq. (12) and b0 is an appropriately chosen constant which ensures a certain efficiency under the Gaussian distribution and a certain breakdown point (BP). The BP is the maximum percentage of outliers for which the estimator is still robust. Further details for choosing b0 can be found in [17]. Again rðÞ can be chosen according to Huber’s r function or Tukey’s biweight function. Here, we use Tukey’s biweight function given by ( 1½1ðd=kÞ2 3 ifjdj rk rTB ðdÞ ¼ ð16Þ 1 ifjdj 4k

with R 2 R being the covariance matrix of Gpq ðt,f ,lÞ, l 2 fL=2, . . . ,L=2g. There are different ways for choosing the

where k is an appropriately chosen constant that depends on a certain efficiency and a breakdown point.

p q

p q

p q

^ x x ðt,f Þ. Then, containing the real and imaginary parts of D p q STFD based on M-estimation can be written as follows: X d^ xp xq ðt,f Þ ¼ argmin rðdðt,f ,l,dpq ÞÞ ð11Þ dpq

l

where dpq ¼ ½RðDpq Þ,I ðDpq ÞT denotes the quantity to be estimated and dðt,f ,l,dpq Þ is the two-dimensional Mahalanobis distance given by dðt,f ,l,dpq Þ ¼ ðGpq ðt,f ,lÞdpq ÞT R1 ðGpq ðt,f ,lÞdpq Þ 22

2634

W. Sharif et al. / Signal Processing 91 (2011) 2630–2638

be written as follows:

6. Non-iterative robust STFD estimation The robust position based techniques are iterative techniques and require a significant computational effort for the estimation of STFDs. In this section, we describe an existing non-iterative robust TFD estimation method and propose some new non-iterative techniques for robust estimation of STFDs. The proposed techniques are STFD estimation based on robust covariance and STFD estimation based on the weighted mean. 6.1. Fractional-lower-order-moment based STFD Djeddi et al. [12] proposed a non-iterative robust estimation of TFD based on a fractional-lower-order moment (FLOM). The FLOM based ðp,qÞth element of the STFD matrix at TF point (t,f) is given by X aS ðt þlÞx/aS ðtlÞej4pfl ^ x x ðt,f Þ ¼ D x/ ð17Þ p q p q l

xp/aS

ja þ 1 =xn

¼ jxp where p and 0 o a o 1. It is evident from the FLOM operation that it influences only the magnitude of x, while the phase remains unchanged. It is required to tune the value of a to get robust estimates. Also it should be noted that a ¼ 1 gives the standard TFD.

^ x x ðt,f Þ ¼ Covðap ,cq ÞCovðbp ,dq Þ D p q þjðCovðap ,dq Þ þCovðbp ,cq ÞÞ

ð20Þ

where Cov denotes the covariance operator. Once having represented the STFD matrix in terms of simple covariance operations, these standard covariances can be replaced by their robust counterparts to get a robust estimate of the STFD matrix. In this paper, we have used two simple methods for the robust estimation of covariance. The first method is based on a robust estimate of scale. For example, the robustified first term Covðap ,cq Þ in Eq. (20) can be written as  2 s^ r ðap Þs^ r ðcq Þ ap cq þ RCovðap ,cq Þ ¼ s^ r 4 s^ r ðap Þ s^ r ðcq Þ  2 ! ap cq ð21Þ  s^ r s^ r ðap Þ s^ r ðcq Þ where s^ r ðÞ denotes the robust estimator of scale. We use, the MADN as the robust estimate of scale. The second robust covariance estimation is based on a function cðÞ that suppresses the outliers. For such kind of robust covariance, the first term in Eq. (20), for example, can be written as RCovðap ,cq Þ ¼ s^ r ðap Þs^ r ðcq Þ

6.2. Robust covariance based STFD estimation

X ap ðlÞm^ ðap Þ cq ðlÞm^ ðcq Þ r r c c s^ r ðap Þ s^ r ðcq Þ l

ð22Þ Robust estimation of covariance under contaminated observations has been already studied in the statistical literature [18]. The idea here is to define each element of ^ x x ðt,f Þ as the estimation of a pairwise the STFD matrix D p q covariance between sequences of two random variables. Consider again the ðp,qÞth element of STFD matrix, given by L=2 X

^ x x ðt,f Þ ¼ D p q

xp ðt þlÞxnq ðtlÞej4pfl

where c is a bounded monotone or redescending cðÞ function, s^ r and m^ r are robust estimates of scale and location, respectively. Here, for the evaluation of the STFD matrix we use cðÞ ¼ sgnðÞ, the median for robust location and MADN for robust estimation of scale. The algorithm to compute the STFD matrix based on robust covariance estimation is summarized in Table 2.

ð18Þ

l ¼ L=2

6.3. Weighted mean based STFD

Define As described in Section 3, we require an averaged STFD matrix computed over all the TF points belonging to a single source for DOA estimation. The ðp,qÞth element of the averaged STFD matrix can be also written as

xp ¼ fxp ðt þ lÞ, l 2 fL=2, . . . ,L=2gg yq ¼ fxnq ðtlÞej4pfl , l 2 fL=2, . . . ,L=2gg xp ,yq 2 CðL þ 1Þ1 . Let Rðxp Þ, Rðyq Þ and I ðxp Þ, I ðyq Þ denote the real and imaginary parts of xp and yq , respectively. Then (18) can also be written as L=2 X

^ x x ðt,f Þ ¼ D p q

Rðxp ðlÞÞRðyq ðlÞÞ



I ðxp ðlÞÞI ðyq ðlÞÞ þ j

l ¼ L=2

j

L=2 X

L=2 X

Rðxp ðlÞÞI ðyq ðlÞÞ

l ¼ L=2

I ðxp ðlÞÞRðyq ðlÞÞ

n n

Table 2 Robust covariance based STFD estimation.

l ¼ L=2 L=2 X

L=2 X 1 X d^ xp xq ¼ Gpq ðtn ,fn ,lÞ #S ðt ,f Þ2S l ¼ L=2

ð19Þ

l ¼ L=2

where xp ðlÞ denotes the lth element of the vector xp . We can also interpret Rðxp Þ,I ðxp Þ,Rðyq Þ and I ðyq Þ as realizations of random variables ap ,bp ,cq and dq, respectively. Then each of the summations in (19) is the covariance of the corresponding two random variables. Hence (19) can also

For a given (t,f) point and for a certain (p,q): (1) Compute the robust estimates of scale and location for each of the vectors ap, bp, cq and dq. MADN and median can be used as robust estimates of scale and location, respectively (2) Compute the robust estimate of the covariance function either by using Eq. (21) or (22) for each term in Eq. (20) (3) Compute the (p,q)th element of STFD matrix using Eq. (20) (4) Repeat Steps (i)–(iii) for all ðt,f Þ 2 S and all ðp,qÞ 2 f1,mg. S denotes the set of all the TF points belonging to a single source

ð23Þ

W. Sharif et al. / Signal Processing 91 (2011) 2630–2638

where S is the set which contains all the (t,f) points ^ x x Þ,I ðD ^ x x ÞT belonging to a single source, d^ xp xq ¼ ½RðD p q p q and Gpq ðtn ,fn ,lÞ be defined as in Section 5.2. We remarked from the simulations that the elements Gpq ðtn ,fn , . . .Þ for different values of ðtn ,fn Þ 2 S are almost identically distributed when constant amplitude FM signals are considered. Therefore, instead of computing the inner summation in (23) for different ðtn , fn Þ separately and then averaging, we propose to compute the double summation in (23) simultaneously. Let ðt1 ,f1 Þ, . . . ,ðtZ ,fZ Þ 2 S denote the TF points belonging to a single source and Z denote the total number of TF points. Further, let Gpq ¼ ½Gpq ðt1 ,f1 ,L=2Þ, . . . ,Gpq ðt1 ,f1 ,L=2Þ, Gpq ðt2 ,f2 ,L=2Þ, . . . ,Gpq ðt2 ,f2 ,L=2Þ, . . . , Gpq ðtZ ,fZ ,L=2Þ, . . . ,Gpq ðtZ ,fZ ,L=2ÞT

ð24Þ

2635

this approach. The procedure to compute the weighted mean based STFD is summarized in Table 3. In the next section, we provide some simulation examples and discuss the results obtained from the different robust estimation approaches. 7. Simulation results We use a uniform linear array (ULA) of eight sensors and two FM sources impinging from angles y ¼ ½51 41 from broadside. One source is a linear FM signal and the second one is a hyperbolic FM signal. The number N is equal to 128 snapshots. Impulsive noise is generated using an e-contaminated Gaussian mixture, given by np ðtÞ  ð1eÞN C ð0, s2 Þ þ eN C ð0, ks2 Þ

ð29Þ

ZðL þ 1Þ2

denoting the matrix of all the sumwith Gpq 2 R mands in (23). Then, (23) can also be written in the following form: ZðL þ 1Þ 1 X d^ xp xq ¼ Gpq ðzÞ Z z¼1

ð25Þ

In the presence of outliers, summation in (25) is not robust. In order to robustify the above average, we use the weighted mean instead of the mean. The weights are computed based on the outlying measures of the elements of Gpq . The outlying measures can be determined by using the Mahalanobis distance measure given by dðzÞ ¼ ðGpq ðzÞmedðGpq ÞÞT R1 G ðGpq ðzÞmedðGpq ÞÞ

ð26Þ

22

is the robust covariance estimate of Gpq . where RG 2 R The weights are calculated based on Tukey’s biweight function given by ( rðdÞ=d2 if da0 wðdÞ ¼ ð27Þ if d ¼ 0 6=k2 where rðÞ for Tukey’s biweight function is given in Eq. (16), k 40 is an appropriately chosen constant. Let wz , z 2 f1,Z  ðLþ 1Þg denote the weight computed for each vector in (24), then the weighted mean estimate d^ xp xq can be written as follows: PZðL þ 1Þ wz Gpq ðzÞ ¼1 d^ xp xq ¼ zP ð28Þ ZðL þ 1Þ wz z¼1 This is a computationally simple and fast as compared to the iterative robust position based techniques. It is a noniterative technique and is motivated by the requirement of having an averaged STFD matrix for DOA estimation. The results in the simulation section show the effectiveness of Table 3 Weighted mean STFD. (1) For a given (p,q), compute the matrix Gpq as in (24) (2) Compute the median and robust covariance SG of Gpq and the Mahalanobis distance (26) for each element in Gpq (3) Compute the weights for each element in Gpq based on Tukey’s biweights as defined in (27) (4) Compute the weighted mean d^ x x as in Eq. (28) p q

(5) Repeat the above steps for all p,q 2 f1,mg

where np(t) denotes the noise at the pth sensor and N C ðm, s2 Þ denotes complex circular Gaussian noise with mean m and variance s2 . In (29), e denotes the probability of occurrence of an outlier and k is the multiplier to indicate the high variance of the outliers with k b 1. A window length Lþ 1 of 37 is used for the computation of the STFD matrix. The robust position based methods require fixing some tuning parameters. For the S-estimator, the required parameters b0 ,k depend on the desired efficiency under the Gaussian noise and the BP of the estimator. We adjust the efficiency to 95% and the BP to 30%. The parameter c for M-estimation is chosen to be equal to 2. The robust covariance based techniques do not require any tuning parameter. The weighted mean based approach requires adjusting the parameter k which is equal to 5 in this simulation. For FLOM STFDs, we set a ¼ 0:3. In the first simulations setup, we use k ¼ 75 and e ¼ 0:2 for (29). The snapshots are generated with a varying signalto-noise ratio (SNR). STFD matrices are computed using the standard and the proposed techniques, namely pre-processing, robust position and non-iterative robust techniques. The resulting averaged STFD matrices are then used to estimate the DOAs using the MUSIC algorithm [10]. The root mean square error (RMSE) of the estimated DOA for the first source is plotted against the SNR for the preprocessing, the robust position based and the non-iterative robust STFD estimators. Table 4 summarizes the robust STFD estimators which are assessed in the simulations. Fig. 2 depicts the RMSE of pre-processing techniques against SNR. The spatial sign function yields biased estimates when the SNR increases. The normalization and the 3-s methods yield similar results with an increasing SNR. The pre-processing techniques provide a lower RMSE than the one achieved by the standard STFD estimation technique. The Crame´r–Rao bound (CRB) is also included to Table 4 Robust methods for STFD matrices estimation. Pre-proc.

Position based

Non-iterative

Spatial sign 3-s Normalization

Median M-estimation S-estimation

FLOM Rob. cov. Weighted mean

2636

W. Sharif et al. / Signal Processing 91 (2011) 2630–2638

102

102 Standard Robust Covariance I Robust Covaraince II Weighted Mean FLOM CRB

Standard Spatial Sign Normalization 3 Sigma CRB

101 RMSE in deg

RMSE in deg

101

100

10−1

10−1 −20

−15

−10 SNR (dB)

−5

102

−20

0

−15 −10 SNR (dB)

−5

0

Fig. 4. Non-iterative robust STFD estimation techniques under impulsive noise.

Fig. 2. Pre-processing based techniques under impulsive noise.

102

Standard Median M estimation S estimation CRB

Standard Spatial Sign Normalization 3 Sigma CRB

101

101 RMSE in deg

RMSE in deg

100

100

10−1

100

10−1

−20

−15 −10 SNR (dB)

−5

0

Fig. 3. Robust position based estimation techniques under impulsive noise.

give an insight. The CRB derived in [11] is calculated for the Gaussian mixture model of (29). The results show that the computationally simple pre-processing techniques are effective in suppressing the impulsive interference. Fig. 3 shows the comparison of different robust position based estimation techniques. The proposed M-estimation and S-estimation based techniques provide lower RMSE if compared with the already existing iterative median-based technique. The robust position based techniques are computationally expensive. Fig. 4 shows the RMSE of non-iterative robust techniques under impulsive noise. The proposed techniques are the robust covariance and weighted mean based approaches. The robust covariance based STFD and weighted mean

−20

−15

−10 −5 SNR (dB)

0

5

Fig. 5. Robust pre-processing STFD based DOA estimation techniques under the Gaussian noise.

approach provide better RMSE characteristics if compared with the FLOM based approach. The simple weighted mean approach provides the lowest RMSE of all the non-iterative robust techniques. These techniques are computationally simpler than the robust position based estimation techniques. In the second simulation, we compare the proposed estimators when the noise is perfectly Gaussian, i.e., there are no outliers. Figs. 5–7 show the results of RMSE vs. SNR for the pre-processing, the robust position based and the non-iterative robust estimation techniques, respectively. In the case of the pre-processing technique, the preprocessed data give a higher RMSE when compared with the standard estimation technique. The spatial sign covariance function provides a higher bias as the SNR increases.

W. Sharif et al. / Signal Processing 91 (2011) 2630–2638

102

102 Standard Iterative (median) M estimation S estimation CRB

101

Standard Median M−estimation S−estimation Weighted Mean Rob Cov FLOM 3−Sigma

101 RMSE (deg)

RMSE in deg

2637

100

100

10−1 10−1

−20

−15

−10 −5 SNR (dB)

0

5

Fig. 6. Robust position based STFD estimation techniques under the Gaussian noise.

102 Standard Robust Covariance I Robust Covaraince II Weighted Mean FLOM CRB

RMSE in deg

101

100

10−1 −20

−15

−10 −5 SNR (dB)

0

5

Fig. 7. Non-iterative robust STFD estimation techniques under the Gaussian noise.

The 3-s rejection rule yields DOA estimates comparable to the estimates obtained by standard techniques. In the case of robust position based estimates, all the three techniques provide better estimates when compared with the estimates based on standard techniques. This is due to the fact that the noise of the STFD becomes a mixture of the Gaussian and the Laplacian even though the noise in the observations is the Gaussian [22]. This means that robust methods are suitable even if there is no contamination in the signal. This leads to a performance gain for M-estimator and S-estimator compared to the standard and median-based estimators. The median-based technique lacks efficiency in the Gaussian case.

0

0.1

0.2 0.3 epsilon

0.4

0.5

Fig. 8. Sensitivity with respect to percentage of outliers at an SNR of  10 dB and with varying percentage of outliers and k ¼ ½5 10 20 40 50 60 70.

For the non-iterative robust estimation, the resulting RMSE is lower than the one achieved by the standard techniques for low to moderate SNR values. The proposed weighted mean approach provides the lowest RMSE among the non-iterative techniques. Again CRB under the Gaussian noise derived in [21] is also included. In the third simulation, the sensitivity of different estimators is compared with respect to the percentage of outliers. Fig. 8 shows the resulting RMSE vs. e. It can be seen that the standard method departs from the true estimates even for a very small fraction of outliers. The FLOM based method starts worsening afterwards. For the robust position based estimation, the S-estimation based STFD and median based STFDs provide a reasonable RMSE for a contamination of up to 40%. The M-estimation based STFD yields correct estimates of DOA for only 30% of contamination. For the non-iterative robust estimation, the robust covariance based estimation techniques could handle 30% of outliers. The weighted-mean based STFD is able to provide good estimates of DOA even up to 40% of outliers. The 3-s rejection pre-processing departed from the true estimates with 30% of outliers. Hence, most of the proposed robust estimators of STFD are able to handle up to 30% of contamination. In the fourth simulation, the sensitivity of the estimators with respect to the impulsiveness factor k is considered. From Fig. 9, it is evident that all the proposed robust estimators provide a good impulsive rejection over a large range of impulsiveness. From the above results, it is clear that the different proposed robust approaches provide comparable results at lower SNR values. Hence, pre-processing is a suitable and constitutes an effective technique for impulsive noise suppression in low SNR regimes. While in the case of low to moderate SNR values, the proposed weighted mean algorithm is a better alternative than pre-processing techniques.

2638

W. Sharif et al. / Signal Processing 91 (2011) 2630–2638

Universit¨at Darmstadt. We are also thankful to Prof. Stannat from TU Darmstadt for his input and valuable comments.

102

RMSE (deg)

101

Standard Median M−estimation S−estimation Weighted Mean FLOM Rob Cov 3*Sigma Rejection

References

100

10−1 101

102 kappa

Fig. 9. Sensitivity with respect to impulsiveness with an SNR of  5 dB and e ¼ 0:2.

8. Conclusion We have addressed the problem of STFD estimation for frequency modulated (FM) sources in the presence of impulsive noise. We focussed on two different approaches, namely pre-processing and robust methods. We have proposed different pre-processing algorithms which are computationally simple and effective in suppressing the impulsive interference. In the robust position based estimators, we have extended the existing position based estimators by applying robust statistical techniques, such as, M-estimation and S-estimation. The proposed position based techniques provide an opportunity to improve the performance of DOA estimation under the Gaussian noise as well. The position based estimators are computationally expensive in general. Therefore, we have also proposed simpler and non-iterative robust estimation techniques. The proposed non-iterative STFD estimation techniques outperform the existing non-iterative robust techniques. We have considered the sensitivity of the proposed estimators against the percentage of contamination and against varying impulsiveness. The proposed estimators yield reasonable results up to a contamination of 30% and all the estimators show good impulsive rejection capability.

Acknowledgments This work is supported by the Excellence Initiative of the German Federal and State Governments and the Graduate School of Computational Engineering at Technische

[1] M. Amin, Y. Zhang, Direction finding based on spatial timefrequency distribution matrices, Digital Signal Processing 10 (4) (2000) 325–359. [2] Y. Zhang, M. Amin, Spatial averaging of time-frequency distributions for signal recovery in uniform linear arrays, IEEE Transactions on Signal Processing 48 (10) (2000). [3] Y. Zhang, W. Mu, M. Amin, Subspace analysis of spatial timefrequency distribution matrices, IEEE Transactions on Signal Processing 49 (4) (2001) 747–759. [4] L. Cirillo, A.M. Zoubir, M. Amin, Parameter estimation for locally linear FM signals using a time-frequency Hough transform, IEEE Transactions on Signal Processing 56 (9) (2008) 4162–4175. [5] P. Heidenreich, L.A. Cirillo, A.M. Zoubir, Morphological image processing for FM source detection and localization, Signal Processing 89 (2009) 1070–1080. [6] W. Sharif, P. Heidenreich, A.M. Zoubir, Robust direction-of-arrival estimation for FM sources in the presence of impulsive noise, in: Proceedings of 35th IEEE Conference on Acoustic, Speech and Signal Processing (ICASSP), Dallas, USA, March 2010, pp. 3662–3665. [7] B. Boashash, Time Frequency Signal Analysis and Processing—A Comprehensive Reference, Elsevier, 2003. [8] L. Cohen, Time-Frequency Analysis, Prentice-Hall, 1995. [9] H. Krim, M. Viberg, Two decades of array signal processing research, IEEE Signal Processing Magazine 13 (4) (1996) 67–94. [10] R. Schmidt, Multiple emitter location and signal parameter estimation, IEEE Transactions on Antennas Propagation, AP 34 (1986) 276–280. [11] B.M. Sadler, R.J. Kozick, Maximum-likelihood array processing in non-Gaussian noise with Gaussian mixtures, IEEE Transactions on Signal Processing 48 (12) (2000) 3520–3535. [12] M. Djeddi, M. Benidir, Robust polynomial Wigner–Ville distribution for the analysis of polynomial phase signals in a-stable noise, in: IEEE International Conference on Acoustics, Speech and Signal Processing, Montreal, Canada, May 2004, pp. 613–616. [13] B. Barkat, L. Stankovic, Analysis of polynomial FM signals corrupted by heavy-tailed noise, Signal Processing 84 (2004) 69–75. [14] D. Middleton, Non-Gaussian noise models in signal processing for telecommunications: new methods an results for class A and class B noise models, IEEE Transactions on Information Theory 45 (4) (1999) 1129–1149. [15] X. Dong, H. Weng, D.G. Beetner, T.H. Hubing, D.C. Wunsch, M. Noll, ¨ H. Goksu, B. Moss, Detection and identification of vehicles based on their unintended electromagnetic emissions, IEEE Transactions on Electromagnetic Compatibility 48 (4) (2006) 752–759. [16] S. Visuri, H. Oja, V. Koivunen, Subspace-based direction-of-arrival estimation using nonparametric statistics, IEEE Transactions on Signal Processing 49 (9) (2001) 2060–2073. [17] R.A. Maronna, R.D. Martin, V.J. Yohai, Robust Statistics, Theory and Methods, John Wiley & Sons Ltd, 2006 ISBN:0-470-01092-4. [18] P.J. Huber, E.M. Ronchetti, Robust Statistics. Wiley Series in Probability and Statistics, second ed., 2009. [19] M. Sahmoudi, K. Abed-Meraim, Robust quadratic time-frequency distributions for the analysis of multicomponent FM signals in heavytailed noise, in: Proceedings of 13th IEEE Workshop on Statistical Signal Processing, Bordeaux, France, July 2005, pp. 871–876. [20] L. Rankine, M. Mesbah, B. Boashash, IF estimation for multicomponent signals using image processing techniques in the timefrequency domain, Signal Processing 87 (6) (2007) 1234–1250. [21] P. Stoica, A. Nehorai, MUSIC, maximum likelihood, and Crame´r–Rao bound, IEEE Transactions on Acoustics, Speech and Signal Processing 37 (5) (1989) 720–740. ¨ [22] I. Djurovic, L. Stankovic, F. Bohme, Robust L-estimation based forms of signal transforms and time-frequency representations, IEEE Transactions on Signal Processing 51 (7) (2003) 1753–1761. [23] D. Pena, F.J. Prieto, Multivariate outlier detection and robust covariance matrix estimation, Technometrics 43 (3) (2001) 286–310.