DOA Estimation in Impulsive Noise Based on FISTA Algorithm

: This paper investigates the challenging problem of direction-of-arrival (DOA) estimation in impulsive noise and presents a fast iterative shrinkage-thresholding algorithm (FISTA)-based approach to tackle the difﬁculty. More speciﬁcally, the underlying noise is modelled as the superposition of outliers in the white Gaussian noise. Leveraging on the spot-sparse characteristic of the outlier matrix, the FISTA is conducted on each snapshot of the array output. With the estimated outlier matrix and the coarse on-grid DOA estimates, an alternating optimization method is developed to retrieve the ﬁnal off-grid DOA estimates. Simulation results show that the proposed method outperforms existing methods in terms of resolution capability and estimation accuracy especially in severe noise environments.


Introduction
In the community of array signal processing, direction-of-arrival (DOA) estimation of multiple emitting sources has received considerable attention for its wide range of applications in radar, sonar, and wireless communications [1,2]. The most well-known subspace-based high-resolution DOA estimation algorithms include MUSIC [3] and ES-PRIT [4]. They retrieve DOA estimates by exploiting the orthogonality between the noise subspace and array manifold which are obtained by eigenvalue decomposition (EVD) of the sample covariance matrix of the array output [3][4][5]. It is worth noting that these methods were developed under the assumption of Gaussian noise. However, in some scenarios, sudden bursts may exist in the received data and it is no longer appropriate to model the ambient noise as Gaussian [6]. For instance, the human-made noise created by power lines, heavy current switches, and other sources cannot be assumed to be Gaussian [7]. Instead, this form of noise is described as impulsive noise, whose probability density function (PDF) has heavier tails than the Gaussian distribution, and hence, it has a higher probability of occurrence of abnormally large values [8].
One way to deal with DOA estimation in impulsive noise is modelling the noise as a complex symmetric alpha-stable (SαS) process and making use of the fractional lower-order moments (FLOMs). With this model and assuming that signals and additive noise are jointly SαS, a robust covariation-based MUSIC (ROC-MUSIC) algorithm was derived in [9]. In [10], a FLOM-MUSIC algorithm was developed by applying FLOMs to the scenarios where the signals are circular and the additive noise is SαS with 1 < α < 2. To estimate the parameters of deterministic signals corrupted by impulsive noise, it was proposed in [11] to preprocess the data by passing them through a Gaussian-tailed zero-memory non-linearity, which leads to low-variance estimates and better performance than covariation-based approaches. In addition, robust statistics [12], such as M-estimation, have also been applied to DOA estimation and tracking in impulsive noise [13][14][15]. Unlike the aforementioned methods using either knowledge on the PDF or number of mixtures, a non-parametric statisticsbased DOA estimator was presented in [16] that estimates the signal and noise subspaces from the spatial sign or rank covariance matrix. Note that although the above-mentioned algorithms adopt different strategies, in general they obtain the signal or noise subspace by performing EVD of the covariance or covariance-like matrices. On the contrary, it is shown in [8] that the estimate of signal subspace is obtained by directly solving the minimization of the l p -norm of the residual fitting error matrix. The resulting l p -MUSIC algorithm achieves performance superiority over the traditional subspace-based algorithms.
In general, the above-mentioned subspace-based DOA estimation methods have limited performance in challenging scenarios such as high source correlation, low sample size, and unknown model order. Alternatively, sparse signal recovery (SSR) methodology is preferred in these scenarios and the problem of DOA estimation by employing SSR has been extensively studied in the past decade [17][18][19][20]. Particularly, for SSR-based DOA estimation in impulsive noise, it is proposed in [21] to model the noise as a mixed Gaussian distribution, i.e., a Gaussian process with an extremely high variance representing outliers superimposed on the ambient Gaussian noise, and perform DOA estimation from the perspective of sparse Bayesian learning (SBL). In [22], Lee et al. examined the scenario when faulty elements appear at the sensor array and develop a weighted SSR-based DOA estimation algorithm by introducing weights describing the degree of outliers of each element. One challenge that must be faced for SSR-based DOA estimation methods is the well-known off-grid problem, which is caused by the finite basis obtained by a predetermined grid sampling of a continuous spatial domain, yet the true DOA may lie off the grid [23][24][25].
To address this off-grid problem, in [23], an alternating minimization method is proposed to obtain sparse signal and mismatch simultaneously. In [24], a greedy method based on matching pursuit is developed. In [25], after conducting first-order Taylor expansion on the steering matrix, an interior point method is utilized to gain the corrected DOAs. Besides off-grid issues, computational efficiency is another critical problem that needs to be considered in SSR-based approaches. The heavy computational load often lies in the process of solving the convex optimization problem involved in the algorithm, e.g., using the CVX MATLAB toolbox or employing SBL techniques. To deal with this problem, the fast iterative shrinkage-thresholding algorithm (FISTA) [26] is employed in [27] to solve the l 1 -norm minimization problem involved, and a FISTA-based DOA estimation method is proposed. Benefiting from the advantages of FISTA in terms of computational efficiency, the computational load of this method is significantly reduced in comparison to the traditional SSR-based DOA estimation approaches.
It is worth noting that in impulsive noise environments, the outliers are sparsely distributed in the data. This implies that we may exploit the signal sparsity of the data in the spatial domain and the outlier sparsity in the data matrix. With this concept, a new method to address the problem of DOA estimation in impulsive noise by the methodology of SSR is proposed in this paper. More specifically, we first model the impulsive noise as Gaussian noise superimposed by scattered points with high amplitudes, which represent outliers. Due to the random occurrence of outliers, FISTA cannot be directly conducted on the received array output data. However, from another point of view, these scattered outliers can be considered as having spot-sparse characteristics in the spatial domain. After expansion of the coefficient matrix into a combination of the array manifold matrix and a unit matrix, each snapshot of the array output can be constructed as a sparse signal recovery problem. Performing FISTA on the array output snapshot by snapshot, a rough DOA (i.e., the on-grid part) and an estimate of the complete outlier matrix can be obtained. Finally, with the estimated rough DOA and outlier matrix, an iterative optimization algorithm is developed in this paper to obtain the final off-grid DOA estimates. Compared with the existing DOA estimation methods, the proposed method achieves better DOA estimation performance, especially in highly impulsive noise environments. Moreover, our method does not require a priori knowledge of the number of sources which should be properly determined in advance to implement l p -MUSIC. The computational burden is also significantly reduced in comparison to the method in [21]. Another advantage of our method compared to the method in [27] is that the restriction of uncorrelated sources can be relaxed.
The following notations are utilized throughout this paper. The superscripts (·) T , (·) H , and (·) * denote the transpose, the conjugate transpose, and the conjugate, respectively. ⊗ denotes the Kronecker product of two matrices, and represents the point-wise multiplication of two vectors with the same dimension. · 1 , · 2 , and · F stand for the l 1 -norm, l 2 -norm, and Frobenius norm, respectively. (·) † represents the pseudo-inverse operator. tr(·) is the trace operator. I M is the M × M identity matrix.

Signal Model
Consider Q narrow-band far-field source signals impinging on an M-element uniform linear array (ULA), of which the distance between adjacent sensors is d. The Q signals arrive at the array from distinct directions {θ 1 , ..., θ Q }. The array output of the t-th snapshot can be expressed as vector of data received by the array sensors, s(t) = [s 1 (t), · · · , s Q (t)] T is the Q × 1 vector of the incoming signals, and n(t) = [n 1 (t), · · · , n M (t)] T is the M × 1 vector of the noise. The T snapshots can be expressed in a compact form as Usually, the noise n(t) is assumed to be Gaussian distributed for simplicity. However, in the presence of impulsive noise, sudden bursts can be involved in the array outputs, and these can be considered as outliers [21,22]. In such cases, the following data model can be employed where N is still the Gaussian noise component as in data model (2), and W = [w(1), · · · , w(T)] represents the outlier matrix. Because the outliers are sudden bursts randomly received by the elements of the array, the outlier matrix W possesses spot-sparse characteristics. That is, only few of the elements in W have high magnitudes whereas others are zero.

SSR Model and FISTA Algorithm
If the on-grid case is considered, the sparse representation of the data model (2) can be easily constructed as where A(θ) = [a(θ 1 ), · · · , a(θ N )] ∈ C M×N is an overcomplete basis which is constructed by discretizing the spatial frequency range [−π/2, π/2] at the N grid pointsθ = [θ 1 , · · · ,θ N ]. Typically, N = MK, where K is the oversampling factor satisfying K > M > Q. Generally, a uniform sampling with ϕ n = (n − 1)∆ n , n = 1, 2, · · · N is considered, where ∆ n = π/N is the grid spacing. Under the on-grid assumption, the true spatial frequencies lie exactly on the sampling grid, i.e., where s T q denotes the q-th row in the source symbol matrix S in (2). The row-sparse structure ofS indicates that if the row index n corresponds to the location of the DOAs on the spatial grid, thes T n is non-zero or otherwise zero.
Under the framework of SSR, the sparse signalS can be recovered from the well-known l 1 -norm optimization problem where S 2,1 denotes a mixed norm, i.e., the l 1 -norm of the column vector formed by performing l 2 -norm calculation on each of its rows, and β is a given threshold. The above optimization problem can be solved by standard convex solvers which often suffer from computational inefficiency. In [26], a fast iterative shrinkage-thresholding algorithm (FISTA) is developed to recover the sparse signal s from the following compressed sensing data model: where y ∈ R M is the measurement vector, D ∈ R M×N is the given dictionary matrix, s ∈ R N is the signal of interest, and w ∈ R M is the noise vector. The implementation of FISTA is shown in Algorithm 1, where prox h (·)is Moreau's proximal operator and it is given by and L ∇ f is a Lipschitz continuous gradient. More details of FISTA can be found in [26].
In order to apply FISTA to data model (4), it is assumed in [27] that all the sources are uncorrelated and the noise is white Gaussian noise with a noise power of σ 2 w . Then, the covariance matrix of X is given as where σ 2 n represents the signal power of the source at grid point n. By vectorizing R, we have y =Ās + σ 2 w 1 n whereĀ = [a(θ 1 ) H ⊗ a(θ 1 ), ..., a(θ N ) H ⊗ a(θ N )], and s equals [σ 2 1 , · · · , σ 2 N ] T which is rowsparse. 1 n = [e T 1 , · · · , e T M ] T with e i being a unit vector where the i-th element is one and the remaining are all zeros. Via these means, it is easy to conduct FISTA on (10) to obtain the sparse signal s.
As a matter of fact, the actual spatial frequencies usually do not fall exactly on the discretized gridθ. In this case, the reconstruction accuracy for the signal by data model (4) is deteriorated due to the modelling error for the manifold matrix A(θ). Performing a first-order Taylor expansion on the steering matrix, a more accurate data model can be described as where A(θ) = [a(θ 1 ), · · · , a(θ N )] ∈ C M×N denotes the on-grid part and the symbol ] is the derivative of A(θ),∆ = diag(δ), andδ = θ −θ is the quantization error due to grid mismatch. Let P = ∆S, (11) can be rewritten as where Φ = [A(θ), B(θ)] and C = [S T , P T ] T . By performing FISTA on the vector form of the covariance of Y, the sparse signal is retrieved and the corrected off-grid DOAs are obtained correspondingly.

DOA Estimation in Impulsive Noise
In the case of impulsive noise, the data model (3) is employed. First, we consider the simple on-grid scenario. Using the same overcomplete basis A(θ) ∈ C M×N in (4), we can formulate the sparse representation of (3) as Various robust compressive sensing approaches have been proposed to deal with the outliers in the data. For instance, the Lorentzian-norm is used in [28], and the 1 -norm loss which is able to achieve better performance is reported in [29]. In addition, the p -norm loss which generates the 1 -norm loss is considered in [30]. Specifically, for each snapshot of X, i.e., x(t), the following problem can be formulated to recover the sparse vectors(t): miñ s(t)

x(t) − A(θ)s(t)
p p + µ s(t) 1 (14) where µ is a regularization parameter. Particularly, if p = 1, the problem reduces to the 1 -norm minimization problem in [29]. Note that if p < 1, the problem is non-convex and difficult to solve, and it is also usually not easy to choose the value of p.
Unlike the above-mentioned robust approaches, in this work we exploit the spotsparse characteristics of the outlier matrix. By combining A(θ)S and W, we can rewrite (13) as As described in data model (3), W does not have row-sparse characteristics, which is the feature ofS. This means that, even if we obtain the vector form of the covariance matrix of X in (14), we cannot apply FISTA directly to it to obtain an estimate of the sparse signal. However, if we investigate each snapshot of the array output X, we have Obviously, for the received t-th snapshot x(t), the joint vector ofs(t) and w(t) is row-sparse. For the signal vectors(t), none-zero elements appear at the row indices n which correspond to the location of the DOAs on the spatial grid, though for the outlier vector w(t), non-zero elements appear randomly. Therefore, let us define E = [A(θ) I M ] and v(t) = [s T (t) w T (t)] T ; then, the following problem can be formulated to determine v(t): min v(t) where λ is a regularization parameter. It is known that the above convex problem can be solved via interior point methods. Nevertheless, to be more computationally efficient, we can employ FISTA on x(t) to retrieve the estimate of v(t), and hence, the estimates ofs(t) and w(t). Following this idea, all the signal vectors inS and all the outlier vectors in W could be obtained. Then, by assembling them, the complete signal matrixS and the outlier matrix W could be derived. Finally, by retrieving the non-zero row indices ofS, we can easily obtain the on-grid DOA estimates. However, there are certain concerns that need to be discussed throughout the whole process. First, different from the method in [27] which utilizes the covariance matrix of X, we apply FISTA on one snapshot of X once a time. The complete signal matrix estimateŜ and the outlier matrix estimateŴ are derived by the combination of the estimated vectors of {s(t)} T t=1 and {w(t)} T t=1 . Obviously, it is difficult to ensure estimation accuracy. Second, for the t-th snapshot x(t), there is a large difference in magnitude between the non-zero elements in the signal vectors(t) and the non-zero elements in the outlier vector w(t) when the noise is highly impulsive. In this case, it is hard to obtain an accurate estimate of them simultaneously. In fact, the obtained estimate of w(t) is far more accurate than that of s(t) due to its relatively high magnitude. Finally, only the on-grid case is considered, and the off-grid case needs to be further addressed.
Our scheme can be described as follows. First, we conduct FISTA on all the snapshots of X and obtain all the estimates of {w(t)} T t=1 . By assembling {ŵ(t)} T t=1 , the estimate of the complete outlier matrixŴ could be obtained. In addition, we examine the l 2 -norm value for all the rows inŜ, which is obtained by assembling all the estimates of {s(t)} T t=1 . Then, the row indices with the Q largest l 2 -norm values are marked as the rough on-grid DOA estimates {θ q } Q q=1 . In the next step, we develop an alternate optimization method to obtain the final off-grid DOA estimates based onŴ and {θ q } Q q=1 . Using the first-order Taylor expansion on the steering matrix in (3), we have where A(θ) = [a(θ 1 ), · · · , a(θ Q )] is the on-grid part of the steering matrix, ] is the derivative of A(θ), ∆ = diag(δ), and δ = θ −θ is the quantization error due to grid mismatch. With the estimatedŴ after the FISTA procedure, we can formulate the following optimization problem to obtain the quantization error ∆, Estimating S byŜ = A † (θ)(X −Ŵ), whereθ is the estimate of θ, the quantization error ∆ can be derived as follows Letting the derivative of δ be equal to zero, we havê The main steps of the proposed FISTA-based DOA estimation method for impulsive noise are summarized in Algorithm 2.

Algorithm 2 Proposed FISTA-based method for DOA estimation in impulsive noise
(1) Input: X and A(θ).
Our proposed method does not involve presetting the coefficient p in the l p -norm, as in the l p -MUSIC algorithm, which plays an important role in the estimation accuracy. That is, our method is more robust when the intensity or probability of occurrence of outliers varies. Although the proposed method cannot directly conduct FISTA on the covariance matrix of the array output as in [27] due to impulsive noise, the estimation accuracy of the method can be compensated by the optimization procedure in (17). Moreover, in our method, the uncorrelation restriction for the signals is relaxed. Finally, benefiting from the FISTA technique, the computational load of the proposed method is significantly reduced compared with the method in [21].

Simulation Results
In this section, we present simulations to assess the DOA estimation performance of the proposed method for impulsive noise. In the simulations, we assume two QPSK (quadrature phase shift keying) signals of the same power impinging on a ULA of sensors. The DOAs of the two incoming signals are θ 1 = −10.2 • and θ 2 = 10.3 • , respectively. The number of sensors is M = 12. In each simulation, 100 Monte Carlo trials are run for each algorithm. The probability of resolution and the RMSE (root mean square error) are evaluated for each algorithm for comparison. Among them, a successful resolution is considered when the DOA estimation error for each source is less than 5 • . The RMSE of the DOA estimates is calculated by where K represents the number of the Monte Carlo runs, θ q represents the real DOA of the q-th signal, andθ q,k represents the estimate to θ q in the k-th Monte Carlo run. The l p -MUSIC, the FISTA-based method in [27], and the SBL-based method in [21] are performed for comparison. For l p -MUSIC, as recommended in [8], p = 1.1 is employed.
In the simulations, the variance of the Gaussian noise is assumed to be σ 2 n , and are Bernoulli distributed with P(b i (t) = 1) = p w , that is, p w represents the probability of occurrence of the outliers, and n w (t) is Gaussian distributed with zero-mean and variance σ 2 w . Accordingly, let σ 2 s be the variance of signal, and the signal-to-noise ratio (SNR) is defined as SNR = σ 2 s /σ 2 n , whereas the signal-to-outlier ratio (SOR) is defined as SOR = σ 2 s /σ 2 w . In the first experiment, the SNR is set at 20 dB, and p w is set at 0.2. The number of the data samples is T = 30. Figure 1 shows the performance of the algorithms as a function of the SOR level. Obviously, as the SOR decreases, for each algorithm, the resolution probability deteriorates and the RMSE increases as expected. However, the RMSE of the proposed method demonstrates an improvement over the other three algorithms especially when the SOR is at extremely low levels. Another observation is that the resolution probability of the proposed method remains 1 for all the SOR values, whereas the other three algorithms present a distinct decrease, especially for low SOR values, indicating that our method obtains a more robust DOA estimation performance. Furthermore, this also explains why the RMSEs of the other three algorithms degrade so dramatically at extremely low SOR levels, as these methods yield DOA estimates with large deviations from the true DOAs in some of the trials. In the following two experiments, the SNR and SOR of each algorithm are fixed at 20 dB and −20 dB, whereas p w is set at 0.15 for the former experiment and the number of snapshots for the latter experiment is T = 30. Figures 2 and 3 show the performance of the algorithms as a function of the number of the data samples and p w , respectively. As Figures 2 and 3 indicate, the proposed method presents a more robust performance over the other three methods when few snapshots are available or when p w increases, which indicates a higher occurrence probability of outliers. In the last experiment, we assume that there are two coherent signals impinging on the array. The SNR and SOR are set at 20 dB and −20 dB, respectively. p w is set at 0.2. The number of snapshots is T = 30. Figure 4 shows 100 realizations of the DOA estimates of these four algorithms. As shown in Figure 4, the DOA estimation results of the l p -MUSIC [8], the FISTA-based method [27] and the SBL-based method [21] present a fairly large deviation with the real DOAs in the realizations, and the proposed method gains a much more accurate DOA estimation. This also proves that our proposed method is not only robust for uncorrelated signals but also for coherent signals.

Conclusions
In this paper, we present a new FISTA-based algorithm for DOA estimation in impulsive noise. By modelling the impulsive noise as the superposition of outliers in Gaussian noise, it is convenient to investigate the spot-sparse characteristic of the outlier matrix. For the array output, we implement FISTA via snapshot to obtain the estimates of the outlier matrix and the on-grid DOAs. Based on the estimated outlier matrix and the coarse on-grid DOA estimates, an alternate optimization method is developed to gain the final off-grid DOA estimates. Simulation results demonstrate that the proposed method presents more robust DOA estimation performance in extremely impulsive noise environments. In addition, the proposed method relaxed the uncorrelation assumption for the signals which is needed in the existing FISTA-based DOA estimation method, and its robust DOA estimation performance for coherent signals is demonstrated by simulations.