A Low Complexity DOA Estimation Method of CD Sources in Impulsive Noise

Direction of arrival (DOA) estimation, as the main technology of passive radio monitoring and positioning, has been deeply investigated. However, the DOA for distributed sources is challenging to estimate in environments with impulsive noise. Although many methods have been proposed for DOA estimation, most of them assume that array output signals contain Gaussian noise. Therefore, the performance of these methods is often poor for alpha-stable distributed impulsive noise. Furthermore, subspace-based DOA estimation methods for distributed sources require a two-dimensional (2D) peak search, which increases the consumption of system computing resources. In this paper, a Q-function-based kernel function is proposed, and its properties are derived. On this basis, a novel DOA estimation method is proposed for coherently distributed (CD) sources in impulsive noise. To reduce computational complexity, a Lagrangian quadratic optimization function is derived by approximating the generalized array manifold of the CD source. By solving this optimization function, a 2D peak search can be reduced to several onedimensional (1D) peak searches. The simulation results illustrate that the accuracy and robustness of the proposed method outperform those of existing methods.


I. INTRODUCTION
Considerable progress has been made in direction of arrival (DOA) estimation, an important research field of passive radio monitoring and positioning, in recent decades, driven by the needs of military and civil fields. Consequently, numerous branches of DOA estimation have been derived for different engineering and technical fields [1][2][3][4]. As a milestone of super-resolution DOA estimation, Schmidt introduced the subspace technology into DOA estimation and proposed the multiple signal classification (MUSIC) algorithm [5]. To reduce the computational complexity of MUSIC due to spectral peak search, Roy proposed the estimating signal parameter via rotational invariance techniques (ESPRIT) algorithm [6]. Since the rotation invariance of the signal subspace is utilized, ESPRIT avoids a peak search, thereby improving efficiency. Further, many improved algorithms based on MUSIC and ESPRIT have been proposed in [7][8][9][10][11]. Additionally, these methods focus on the "point source", which assumes that there is only a direct wave with concentrated energy between the source and the array without angular spread [12][13][14]. In most scenarios, the point source model is reasonable and can produce correct DOA estimation results. Nevertheless, as the space electromagnetic environment becomes increasingly complex, multipath propagation generally exists in the process of signal propagation [15]. Many realistic examples, in which the point source assumption is unreasonable, such as undersea echo detection and localization of acoustic sources by a microphone array, have been identified. For subspace-based algorithms, noise subspace is occupied by the signal subspace multipath propagation results. In other words, the noise subspace cannot be effectively obtained. This disadvantage causes the performance of the MUSIC to degrade significantly or even makes it invalid.
To avoid this drawback, researchers have proposed a spatially distributed source model. From the perspective of timevarying signal changes, distributed sources are classified into coherently distributed sources and incoherently distributed sources. If the signal components of a source incident to the receiving array through multiple paths are fully correlated, the source is a coherently distributed (CD) source. In contrast, if these components are uncorrelated, the source is an incoherently distributed (ID) source.
By extending the subspace technology to the distributed source model, [17] proposed the distributed signal parameter estimator (DSPE) to realize DOA estimation. Theoretically, the DSPE can be considered an extension of MUSIC for distributed sources. To estimate DOAs, the DSPE and MUSIC need to perform a two-dimensional (2D) and onedimensional (1D) peak search, respectively. Therefore, the execution efficiency of the DSPE is much lower than that of MUSIC. This disadvantage limits the application of the DSPE in practice. Bengtsson proposed a low-complexity root-MUSIC algorithm [18], which uses two point sources to approximate a distributed source and uses low-complexity algorithms to estimate the DOAs of the two point sources. The mean of the DOAs of the two point sources is considered the center DOA. Zoubir proposed a robust Capon beamforming method for ID sources [19]. This method calculates the autocorrelation matrix of the received signals and searches peaks to estimate the center DOA. [20] proposed a beamforming algorithm, which performs Cholesky decomposition on the autocorrelation matrix of the distributed source signals and applies peak search to estimate the center DOA. By simplifying the steering vector of the CD source to the Schur-Hadamard product between the steering vector of the point source and the real vector and deriving the rotation invariant structure of the subarrays, [21] proposed a generalized ESPRIT algorithm for DOA estimation. Distributed sources have also been addressed in [22][23][24][25][26][27][28].
In practical applications and theoretical analysis, it is usually assumed that a wireless channel contains Gaussian noise. However, studies have revealed that some natural or human-made noise is typically impulsive, for example, lowfrequency atmospheric noise, noise due to sea waves and mountain discontinuities, and underwater acoustic signals [29]. In this case, impulsive noise lacks a finite covariance; therefore, the performance of subspace-based DOA estimation algorithms [30][31][32][33] degrades significantly and cannot produce reasonable results in impulsive noise environments. In contrast to Gaussian noise, the probability density function (PDF) of the alpha-stable distribution has heavy tails, so it is suitable for describing impulsive noise. Many theories have been proposed to improve the performance of DOA estimation in impulsive noise environments. Among them, fractional lower-order statistics (FLOS) [34][35][36] and correntropy [37][38][39][40][41] are the most representative. Based on these theories, fractional lower-order moment-based MU-SIC (FLOM-MUSIC) [42], phased fractional lower-order moment-based MUSIC (PFLOM-MUSIC) [43], correntropybased correlation MUSIC (CRCO-MUSIC) [44], and generalized autocorrentropy-based DSPE (GCO-DSPE) [45] algorithms are derived. FLOM-MUSIC and PFLOM-MUSIC adopt the fractional exponential to ensure the boundedness of the covariance of sample data. CRCO-MUSIC and the GCO-DSPE use the attenuation characteristics of the exponential operation to effectively eliminate the influence of impulsive noise and improve the accuracy and robustness of DOA estimation. However, the selection of the exponential of FLOS and the kernel size of the correntropy requires prior knowledge of arrays and signals, which is not easy to obtain. Another disadvantage is that to ensure the accuracy of DOA estimation, these algorithms require a large amount of sample data. Therefore, these algorithms are subject to considerable restrictions in practice.
In this paper, the estimates of the spread angle and central DOA for distributed sources in impulsive noise are investigated. The main contributions of this work are summarized as follows: • To improve the accuracy and effectiveness of DOA estimation of CD sources under impulsive noise conditions, a new Q-function-based impulsive noise suppression operator is deduced, and its properties are analyzed and proved. Further, a novel DOA estimation method is proposed by combining this operator with the DSPE. • Since the proposed algorithm requires a 2D peak search to obtain the estimates of the spread angle and central DOA, which entails high computational complexity, this paper approximates the generalized array manifold of CD sources and constructs a Lagrangian quadratic optimization function. By solving this quadratic optimization function, the estimates of the spread angle and central DOA are simplified as a 1D peak search. This paper is organized as follows: The alpha-stable distribution, signal model of the CD source, and DSPE algorithm are presented in Section II. A new Q-function-based impulsive noise suppression operator is deduced in Section III based on the characteristics of the Q-function, and its properties are analyzed and proved. Further, to reduce the computational complexity of the proposed algorithm, in Section IV, we approximate the generalized array manifold of CD sources, and the 2D peak search is transformed into several 1D peak searches via quadratic optimization. Section V presents simulation results that verify the performance of different algorithms, and Section VI gives some conclusions.

A. SIGNAL MODEL
The multipath scattering energy of the distributed source presents a certain spatial distribution around the central DOA. Therefore, the distributed source can be regarded as a beam of densely distributed point sources in space. Considering L distributed sources incident into a uniform linear array (ULA) with M elements, and assuming that the scattering point source is continuous, the array output observation can be expressed as where s i (θ, ψ i , t) is the angular signal density of the ith source, ψ i involves the spread angle and central DOA, a(θ) = 1, e −j2πd/λ sin(θ) , · · · , e −j2π(M−1)d/λ sin(θ) T is the steering vector of the ULA, n (t) is the noise vector, and the range of the observation angle is Θ = [−π/2, π/2]. For CD sources, the angular signal density can be represented as where χ (θ, ψ i , t) is a complex-valued deterministic angular signal density function, and ζ i is a random variable. Eq. (1) can be rewritten as in which, b(ψ i ), namely, a generalized steering vector, is given by where i = 1, 2, · · · , L. As expressed by Eq. (4), the generalized steering vector b(ψ i ) is the integral of a(θ) over θ.
If we can further define the generalized array manifold matrix as , under the assumption that the noise and signal are uncorrelated, the correlation matrix of x(t) is given by where R ss (ψ) is a noise-free correlation matrix, and R nn is noise correlation matrix. Γ will be diagonal [46] and its (i, j)th entry is defined as E{ζ i (t)ζ j * (t)}.
If we define the angular cross-correlation kernel (ACCK) as then R ss (ψ) can be re-expressed as Assume that the array output signals from different CD sources are uncorrelated; then, the ACCK can be simplified to where δ i j denotes the Kronecker delta and is the angular autocorrelation kernel (AACK) of the ith source. Thus, Eq. (7) can be rewritten as

B. THE DSPE ALGORITHM
Consider conventional super-resolution subspace-based techniques, where the covariance matrix of the array output signals can be decomposed into noise and signal subspaces via singular value decomposition (SVD). If we denote the noise subspace and signal subspaces as U n and U s , respectively, then, the 2D space spectrum of the distributed source can be expressed as The spread angle and central DOA can be estimated by the following criterion.
where tr(·) denotes the trace of a matrix, and B(ψ) can be expressed as For CD sources, the 2D spatial spectrum can be expressed as The criterion for estimating the angular spread and central DOA can be represented aŝ This algorithm is referred to as the DSPE [17]. However, to estimate the unknown parameter vector ψ, a 2D peak search must be performed, which leads to high computational complexity and additional economic resource consumption.
In practical applications, since the joint probability density function is unknown and only a finite number of data samples are available, R only can be obtained by employed N data samples as followsR

C. ALPHA-STABLE DISTRIBUTION
The alpha-stable distribution, whose PDF has heavy tails, is a continuous probability distribution developed by Paul Pierre Levy [47]. As an effective theoretical tool, the alpha-stable distribution is usually applied to describe the impulsive noise. The thickness of the tails of the alpha-stable distribution can be controlled via the parameter α; thus, the distribution is a flexible modeling tool. Alpha-stable distribution is usually described by its char- VOLUME 4, 2016 acteristic function [48,49] and where 0 < α ≤ 2 is the characteristic exponent. −1 ≤ β ≤ 1 is the symmetric parameter used to determine the sign and degree of asymmetry. −∞ < µ < +∞ is the location parameter. In particular, if β = 0, this distribution is referred to as the symmetric alpha-stable (SαS) distribution, which is symmetric about µ. γ > 0 is the dispersion, which can be applied to measure the dispersion of sample data. The alphastable distribution has two important properties that have a crucial role in the modeling of uncertainty [29].
1) Generalized central limit theorem: Let Y 1 ,Y 2 , · · · ,Y N be independent and identical distribution (i.i.d.) random variables and N → ∞, Y is the limit in the distribution of normalized sums of the form if and only if the distribution of Y is stable. Further, have finite variance and are i.i.d., the limit distribution is Gaussian. 2) Stability property: Let Y 1 and Y 2 be i.i.d. random variables with the same distribution as Y ; for arbitrary constants d 1 and d 2 , there are constants a 1 and a 2 such that Eq. (21) can be interpreted as a 1 Y + a 2 and d 1 Y 1 + d 2 Y 2 have the same distribution. Fig. 1 shows the comparison of alpha-stable distributed impulsive noise and Gaussian noise. Distinct spikes in the impulsive noise are observed, but the Gaussian noise fluctuates within a small range of the mean value.

III. PROPOSED SOLUTION TO ADDRESS IMPULSIVE NOISE
Since the impulsive noise lacks a finite covariance, the performance of the super-resolution subspace-based DOA estimation method will be severely degraded, and reasonable estimation results cannot be obtained. Therefore, this section focuses on suppressing the impulsive noise of array output signals and improving the accuracy of DOA estimation for CD sources.

A. Q-FUNCTION AND ITS PROPERTIES
Considering the standard normal random variable X, the Qfunction [50] Q(x) can be given by Formally, the Q-function is also usually expressed in the form of the following integral Further, the Q-function can be calculated by the error function where erf(·) is the error function. It can be deduced that the domain of the Q-function is (−∞, +∞), and the range is [0, 1]. Therefore, if x contains outliers (impulsive noise), the Q-function can map the amplitude of the sample data to the limited interval [0, 1], which means that the Q-function can suppress impulsive noise. However, from the perspective of signal parameter estimation theory, the impulsive noise suppression effect of the Q-function also has shortcomings: for negative sample data with large outliers, the outputs of the Q-function are nearly 1, and the outputs of the Q-function corresponding to small negative sample data are nearly 0. In this case, the characteristics of the effective signals contained in the sample data cannot be reasonably shown. The Q-function is asymmetric: if the moduli of two sample data are equal but their signs are opposite, and after the Q-function is applied, the values of the data are not equal. This result will affect the performance of parameter estimation.

B. Q-FUNCTION-BASED KERNEL FUNCTION
Considering that the Q-function has the ability to suppress impulsive noise, this paper expands the Q-function as the Q-function-based (QFB) impulsive noise suppression operator. The QFB operator is defined by Inspired by the Hampel identifier [51], in Eq. (25), h(x) is expressed as |x | πmed( |x |) , and med(·) is the median operator. Additionally, we list eight main properties of the QFB operator.
(1) The QFB operator is an even function symmetric with respect to x = 0, and it reaches the maximum at x = 0.
(2) The QFB operator is bounded: The QFB operator is a kernel function. proof: According to property (2), f QFB (x) is bounded, so it is a positive definite function. Furthermore, from property (1), we know that f QFB (x) is symmetric. Therefore, f QFB (x) is a kernel function that is referred to as the QFB kernel.
(4) Consider that X(t) is a random process, then is a reproducing kernel. proof: is also a symmetric function. Moreover, since f QFB (X(t 1 ) − X(t 2 )) satisfies the positive definite condition, for any set of real numbers {d 1 , d 2 , · · · , d n } that are not all zero Because the strictly positive definite function for two random variables is greater than zero, Therefore, C(t 1 , t 2 ) satisfies not only positive definiteness but also boundedness. According to the Moore-Aronszajn theorem [52], for each real symmetric positive definite function of two real-valued random variables, there is a unique Hilbert space with the kernel function as the reproducing kernel [53,54]. Therefore, C(t 1 , t 2 ) is a reproducing kernel.
(5) Let X and Y be two random variables; C(X,Y ) is a second-order statistic of the mapped feature space data.
proof: For the QFB kernel, we can obtain where ϕ(·) denotes a nonlinear mapping, which is induced by the QFB kernel.
Further analysis shows that We can obtain a mean square cost function induced by the QFB kernel in the feature space Let N sample data {(x n , y n )} N n=1 be subject to the joint probability density f XY , and u n = x n − y n . The estimate of J(X,Y ) can be written aŝ It can be deduced that J(X,Y ) is concave. Notably, J(X,Y ) can be applied for adaptive filtering and matrix low-rank approximation in impulsive noise environments.
(6) If we define u = |x| /[πmed(|x|)], then the QFB operator can be represented as a continued fraction By truncating the second item in Eq. (33), an upper bound can be expressed as (7) Since the calculation of the QFB operator involves an integral operation, which has high computational complexity, in practical applications, the following approximate method can be employed (8) Let X be a random variable; the expression X · f QFB (X) is bounded. proof: where u = |x| /[πmed(|x|)]. Therefore, X · f QFB (X) is bounded.

VOLUME 4, 2016
Furthermore, if X · f QFB (X) is regarded as a new random variable obtained by the weighting between the random variable X and a factor related to the statistical characteristics of X, this new random variable is bounded. Therefore, this new random variable can be utilized in subspace-based parameter estimation methods.

C. QFB-DSPE ALGORITHM
Since the alpha-stable distributed random variables lack a finite variance, the performance of subspace-based DOA estimation algorithms will be greatly affected. Therefore, a novel QFB-based DOA estimation algorithm, which is referred to as the QFB-DSPE, is proposed by extending the DSPE from Gaussian noise to impulsive noise.
To realize DOA estimation by subspace technology, the correlation matrix of the array output signals must be constructed. However, if the signals contain alpha-stable distributed impulsive noise, the correlation matrix will be unbounded, which will seriously reduce the accuracy and reliability of DOA estimation. To avoid this disadvantage, an M × M pseudo-correlation matrix R is proposed, and its (i, j)th entry can be written as Furthermore, to better identify and suppress impulsive noise, Eq. (37) can be rewritten as where x i and x j are the ith component of x and jth component of x, respectively. The implementation procedures of the QFB-DSPE are expressed as follows: 1) Compute the M × M matrixR based on Eq. (38).
2) Perform SVD onR to obtain the signal and noise subspaces.
3) Construct the 2D spatial spectrum of the QFB-DSPE in Eq. (61). 4) Search the 2D local peaks and obtain the estimates of the spread angle and central DOA.

IV. A FAST PEAK SEARCH METHOD
From the perspective of algorithm execution efficiency, a 2D peak search requires considerable computing resources and has high complexity. The fourth step of the QFB-DSPE algorithm requires a 2D peak search, which limits the application ULA target σ σ ϑ FIGURE 2: Spatial structure of the distributed source of the QFB-DSPE algorithm in situations where real-time requirements are high. Therefore, this section derives a fast method to obtain the spread angle and central DOA for CD sources.

A. PROPOSED SOLUTION FOR THE FAST PEAK SEARCH
According to the spatial structure of the distributed source, as shown in Fig. 2, the energy of multipath signals exists in the spatial area with the central DOA ϑ as the center and spread angle σ as the radius. Further, the energy of multipath signals follows a certain distribution in this area, and the multipath signals outside this area do not contribute to the DOA estimation of the distributed source. Therefore, the deterministic angular signal density function in Eq. (2) is described as a unimodal conjugate symmetric function, which is always modeled as a Gaussian distribution or uniform distribution. If the deterministic angular signal density function is defined as χ (θ − ϑ), then χ (θ − ϑ) = χ (ϑ − θ). That is, the deterministic angular signal density function is symmetric about ϑ.
Assume that the deterministic angular signal density be modeled as a uniform distribution Then, [b (ψ)] m can be computed by Eq. (4): For a small σ, that is, the angular spread is small, [b (ψ)] m can be obtained in a straightforward manner [17] by The generalized steering vector can be rewritten as The generalized steering vector b(ϑ, σ) has the following form where diag(·) denotes a diagonal matrix.
By substituting Eq. (45) into Eq. (14) and according to [55,56], the 2D spatial spectrum of the CD source can be expressed as Eq. (46) can be further rewritten as Therefore, Eq. (46) can be simplified to For the sake of discussion, we define a new function Eq. (50) is a quadratic programming problem. To eliminate the trivial solution h(σ = 0) M×1 when solving ϑ and σ, the following constraint is considered where c is a constant greater than 0. Therefore, the estimation problem of the spread angle and center DOA for CD sources can be converted to solve a quadratic programming problem with equality constraints as follows In the following section, this quadratic programming problem is solved by the Lagrange method.
We construct the Lagrangian optimization function as follows where is a constant. By calculating the partial derivative of L(ϑ, σ) with respect to h(σ), the following formula can be obtained It can be observed from Eq. (54) that h(σ) and can be calculated as follows Establishing the equations in Eq. (56) and solving them, we can obtain Substituting Eq. (57) into Eq. (50) yields Substituting Eq. (58) into Eq. (49), the 1D spatial spectral function containing only the central DOA can be obtained as follows Therefore, the center DOA can be estimated bŷ The spread angle will be estimated by substituting Eq. (60) into Eq. (49).σ = arg max

B. COMPLEXITY ANALYSIS
To simultaneously determine the DOA of L CD sources, a one-time 1D peak search is required to perform the central DOA estimation, and L 1D peak searches are required to estimate the spread angle. Let the angle extension range and step size of a peak search be ∆ and ω, respectively. For the DSPE, a 2D peak search requires ∆π/ω 2 +(∆+π)/ω+1 steps. In the performance analysis of DOA estimation algorithms under Gaussian noise environments, the signal-to-noise ratio (SNR) is considered the main parameter. However, the alphastable distributed impulsive noise lacks a finite second-order statistic for α < 2. Therefore, to evaluate the ratio of the signal power over noise dispersion γ, the generalized signalto-noise ratio (GSNR) [43] is defined by Since the QFB-DSPE algorithm obtains the spread angle and central DOA by searching local peaks, the changes in these two parameters will show similar characteristics with the changes in the GSNR, number of snapshots, and characteristic exponent.
The root mean square error (RMSE) in Eq. (63) is employed to evaluate the estimation accuracy of the spread angle and central DOA [57].
where N is the number of all Monte Carlo (MC) trials and ϑ l (n) andσ l (n) are the estimates of ϑ l and σ, respectively, for the lth MC trial. If the resolution criterion of Eq. (64) holds the two arrival signals are considered resolvable, where ϑ 1 and ϑ 2 are the angles of the two arrival signals and ϑ a is the mid-angle between them. Additionally, p(·) is the reciprocal of the spatial spectrum. The probability of resolution [58][59][60] is the probability of δ(ϑ 1 , ϑ 2 ) > 0.

A. THE EFFECT OF THE GSNR
We define one data sampling of all elements in the ULA as a snapshot, and in these experiments, the number of snapshots is N = 600. For the SαS impulsive noise, the characteristic exponent α is fixed at 1.3 in Fig. 3 and at 1.8 in Fig. 4, which correspond to the cases of strong impulsive noise and weak impulsive noise, respectively. and DSPE for a wide range of GSNRs. As shown in Fig. 3(b), as the GSNR increases from 0 dB to 20 dB, the RMSEs of these algorithms decrease. The RMSE curves of FLOM and PFLOM are similar, and that of CRCO is below them because the correntropy suppresses impulsive noise better than the FLOS. The DSPE algorithm cannot suppress impulsive noise, so its RMSE is the largest. However, the QFB-DSPE always has the lowest RMSE. Fig. 3(a) shows that as the GSNR increases, the ability of the FLOM, PFLOM, CRCO, and DSPE algorithms to distinguish the incident sources is improved. The probability of resolution of CRCO reaches 1 the fastest, and that of DSPE reaches 1 the slowest. However, when the GSNR is kept in [0, 20] dB, the QFB-DSPE can completely distinguish the incident sources. In Fig. 4, we observe that if the characteristic exponent is set to α = 1.8, which is equivalent to moderate impulsive noise, the DOA estimation accuracy and resolution of the incident sources of these algorithms are improved, and the QFB-DSPE still outperforms the other algorithms.

B. EFFECT OF THE CHARACTERISTIC EXPONENT α
The effect of the characteristic exponent α on the performance of DOA estimation, including the RMSE and probability of resolution, is investigated in this section. The GSNR is 5 dB, and the number of snapshots is N = 600. The characteristic exponent is increased from α = 1.0 to α = 2.0. As shown in Fig. 5(a), in α ∈ [0, 1], the CRCO and QFB-DSPE can completely distinguish the incident sources, but the FLOM and PFLOM can only completely distinguish the incident sources when α > 1.5. For the DSPE, if and only if α > 1.8, can the incident sources be distinguished. In terms of the source discrimination ability, the CRCO and QFB-DSPE have the same performance and outperform the other algorithms. As shown in Fig. 5(b), the RMSEs of the DSPE, PFLOM, FLOM, and CRCO increase significantly as the impulsivity of the noise increases. The increase in impulsive components in noise will reduce the accuracy of DOA estimation of these algorithms. However, the RMSE of QFB-DSPE is minimal and does not change; therefore, the QFB-DSPE has not only high DOA estimation accuracy but also high robustness to the variation in impulse noise intensity.

C. EFFECT OF THE NUMBER OF SNAPSHOTS
In Fig. 6, the impact of the change in the number of snapshots on the QFB-DSPE and the other algorithms is shown. The characteristic exponent of the SαS impulsive noise is α = 1.5, and the GSNR is 5 dB. This phenomenon shows that as the quantity of sample data increases, the impulsive components and effective signal components increase simultaneously. Since the DSPE cannot suppress impulsive noise, the increase in effective signals cannot offset the impact of impulsive noise on performance. Therefore, an increase in the number of snapshots does not improve the resolution of probability and RMSE of the DSPE. For the CRCO, QFB-DSPE, FLOM, and PFLOM algorithms, the correntropy, Qfunction, and FLOS are utilized to suppress impulsive noise; so the increase in the number of effective signals helps to improve the performance of these algorithms. However, the QFB-DSPE still yields better performance than the other algorithms, so in terms of suppressing impulsive noise, the QFB is better than the FLOS and correntropy.

D. EFFICIENCY OF DIFFERENT ALGORITHMS
To evaluate the execution efficiency of different algorithms, the GSNR is set to 10 dB; the characteristic exponent of the alpha-stable distributed impulsive noise is set to α = 1.5; the number of snapshots is set to N = 600; and the number of elements in the ULA is set to M = 10 and M = 20, respectively. Each algorithm performs 500 iterations, and the average consumption time is shown in TABLE 1.
As shown in TABLE 1, for a certain algorithm, with an increase in the number of elements in the ULA, the consumption time for this algorithm increases due to the increase in received data with the increase in elements. When the number of elements is fixed, since the DSPE does not include  an impulsive noise suppression mechanism, it consumes less time than the FLOM, PFLOM, and CRCO. However, the time consumption of the QFB-DSPE is the least, which not only shows that the QFB-DSPE has the lowest complexity but also shows that the 2D spectral peak search in the DOA estimation algorithm of distributed sources has the highest complexity.

VI. CONCLUSIONS
We present a new DOA estimation of the CD source algorithm with low computational complexity in impulsive noise. Based on a comprehensive analysis of the Q-function, the QFB impulsive noise suppression operator is deduced, and its properties are analyzed and proved. By combining the QFB operator with the subspace technology, this paper proposes the QFB-DSPE algorithm to estimate the spread angle and central DOA of the CD source. The QFB-DSPE requires a 2D peak search, which consumes a substantial amount of computing resources and has low efficiency. To improve the execution efficiency of the QFB-DSPE to achieve a spread angle and central DOA estimation, this paper proposes a fast algorithm by transforming the 2D peak search to multiple 1D peak searches. Experimental results illustrate that the QFB-DSPE is better than the other algorithms in terms of both source separation capability and estimation accuracy. Moreover, the QFB-DSPE shows high robustness in environments with different intensities of impulsive noise.