# Gridless DOA Estimation for Minimum-Redundancy Linear Array in Nonuniform Noise

Authors: Changyun Qi; Gong Zhang; Jiawen Yuan
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

A gridless direction-of-arrival (DOA) estimation method to improve the estimation accuracy and resolution in nonuniform noise is proposed in this paper. This algorithm adopts the structure of minimum-redundancy linear array (MRA) and can be composed of two stages. In the first stage, by minimizing the rank of the covariance matrix of the true signal, the covariance matrix that filters out nonuniform noise is obtained, and then a gridless residual energy constraint scheme is designed to reconstruct the signal covariance matrix of the Hermitian Toeplitz structure. Finally, the unknown DOAs can be determined from the recovered covariance matrix, and the number of sources can be acquired as a byproduct. The proposed algorithm can be regarded as a gridless version method based on sparsity. Simulation results indicate that the proposed method has higher estimation accuracy and resolution compared with existing algorithms.


---

## Body

## 1. Introduction

Direction-of-arrival (DOA) estimation is one of the most important topics in the field of array signal processing [1]; its main purpose is to estimate the angle information of unknown signals spatially based on the array sensor model. Many existing DOA estimation algorithms have good estimation accuracy and resolution [2–5] in additive white Gaussian noise. However, in practical applications, the noise is usually considered to be nonuniform due to the misalignment of the antenna array or the nonidealities of the receiving channels. This nonuniformity leads to a decrease in the estimation accuracy of many DOA estimation algorithms based on the Gaussian white noise.In recent decades, many measures have been taken to improve the performance of DOA estimation in nonuniform noise. He et al. [6] proposed a DOA estimation algorithm based on the sparse representation of the covariance vector. A linear transformation matrix is designed to remove the nonuniform noise power component, and then it is converted to solve the sparse reconstruction problem based on L1 norm. This algorithm is an improvement of the L1 norm-based array covariance vector sparse representation (L1-SRACV) algorithm [7], so it is also called improved L1-SRACV (IL1-SRACV) algorithm. Yet the diagonal term information of the signal covariance is lost in the process of removing nonuniform noise power. Yang et al. [8] proposed a sparse parameter estimation technique (SPA) without discretization, using the same covariance matching criteria as sparse iterative covariance-based estimation (SPICE) [9] to transform the covariance fitting problem into a positive semidefinite programming (SDP) [10] problem to solve, and then the postprocessing technology is used to estimate the target position, noise variance, and other information from the covariance matrix. This algorithm has many significant advantages, but its performance is poor under low signal-to-noise ratio (SNR). Liao et al. [11] proposed a rank and trace minimization (RTM) algorithm, which transforms the rank minimization problem of the covariance matrix into the noise power maximization problem and then uses the received signal and the nonuniform noise covariance matrix to perform difference to achieve DOA estimation. The algorithm is simple to implement and has low computational complexity, but the structural characteristics of the signal covariance matrix are not taken into account, and the estimation accuracy is low.Besides, most of the algorithms are based on uniform linear arrays, which results in the resolution being affected by the array aperture. Using well-designed sparse arrays, such as minimum-redundant arrays (MRAs) [12], nested arrays [13], and coprime arrays [14], fewer array elements can be used to obtain larger antenna apertures, and therefore, it becomes an effective way to improve resolution. Among them, MRA is defined as the minimum redundancy of the array element position difference under the condition that the number of array element position difference is continuous. Given the number of array elements, the relative position can be easily determined by looking up the minimum-redundant array structure given in reference [15].In order to overcome the defects of the existing DOA estimation algorithms in nonuniform noise and combine the advantages of sparse arrays, this paper proposes a gridless DOA estimation algorithm for MRA in nonuniform noise. First, the matrix minimization problem is solved to obtain the noise-removed signal covariance matrix. Since the actual covariance matrix is obtained from the finite snapshot, even if the noise is suppressed in the early stage, the covariance matrix estimation error caused by the finite snapshot still exists. In order to reduce the estimation error and improve the estimation accuracy, a gridless residual constraint scheme that does not rely on noise parameters is designed, and the noiseless Toeplitz matrix is reconstructed. By operating on this matrix, we can not only obtain target angle information but also get the number of sources as a byproduct. The proposed algorithm can be viewed as a gridless version of method based on sparsity, which overcomes the problem of basis mismatch in compressed sensing algorithms.


## 2. Signal Model and Assumptions

Suppose that a sensor array composed ofM elements receives K sources in unknown directions, taking the first array element as the reference array element and the position information of the array element D=d1,d2,...,dM=d×0,c1,c2,...,cM−1, where d is usually taken as the half wavelength of the signal. The schematic diagram of the array reception of the sparse linear array is shown in Figure 1.Figure 1
Schematic diagram of sparse linear array reception.AssumingK < M, the received data can be formulated as follows:(1)yt=Aθst+nt,t=1,2,…,L,where t is the snapshot index, L stands for the number of snapshots, yt=y1t,y2t,…,yMtT∈ℂM denotes the observation vector, st=s1t,s2t,…,sKtT∈ℂK represents the source signal vector, andnt=n1t,n2t,…,nMtT∈ℂM is the noise vector, in which the superscript T means transpose operation. ℂ represents the set of complex numbers, θ=θ1,θ2,…,θKT∈ℂK,whereθk∈−90∘,90∘ denotes unknown source directions, and Aθ=aθ1,aθ2,…,aθKT∈ℂM×K denotes the array manifold matrix, where the steering vector of the kth signal aθK satisfies the following equation:(2)aθk=1,ejπc1sinθk,…,ejπcM−1sinθkT∈ℂM.It is easy to see the sparse array can be regarded as a uniform linear array with missing elements. In other words, the received signal can be further expressed as follows:(3)yt=ΦA1θst+Φn1t,t=1,2,…,L,where A1θ∈ℂC×K refers to the array manifold matrix of a uniform linear array with C sensors, where C=cM−1+1, Φ∈0,1M×C represents a selection matrix, and in the ith row of which, only the position of the column corresponding to (ci−1+1) is 1, and the others are 0. Since the aperture of the array has been expanded from (M − 1)d to (C − 1)d, a sparse array composed of the same number of elements has a higher resolution.Some standard assumptions are formulated here for the solution of the problem. We assume that the signal and the nonuniform noise are considered to be uncorrelated spatially and temporarily and independent of each other:(4)Est1sHt2=Pδt1,t2,Ent1nHt2=Qδt1,t2,Est1nHt2=Ent1sHt2=0,Q=diagσ,σ=σ12,σ22,...,σM2T,where E{∙} denotes the mathematical expectation and (∙)H represents conjugate transpose operation. δt1,t2 represents the Kronecker delta function, the value of which is one when t1=t2, and otherwise, it is zero. P∈ℂK×KandQ∈ℂM×M represent the signal correlation matrix and noise correlation matrix, respectively, diag(∙) denotes the diagonal matrix formed by the vectors in brackets, and σm2 represents the noise variance superimposed on the mth sensor.Based on the above assumptions, the covariance matrix of the signal received by an array can be modelled as follows:(5)R=EytyHt=AθPAHθ+Q=ΦA1θPA1HθΦH+ΦQ1ΦH.It should be noted that in practical applications, the covariance matrix is usually replaced by the sampling covariance matrixRe:(6)Re=1L∑t=1NytyHt.


## 3. Rank Minimization-Based Gridless DOA Estimation

The gridless DOA estimation algorithm refers to the angle estimation in the continuous domain to avoid discretization of the angle [16], which effectively overcomes the base mismatch problem in the compressed sensing algorithm [17] and provides a fresh way to improve the estimation accuracy and resolution.


### 3.1. Denoising Method Based on Rank Minimization

For noncoherent sources,rankRe−Q=rankΦA1PA1HΦH=K, where rank(∙) denotes the rank operation of the matrix in brackets. For any diagonal matrix Λ composed of real numbers, rankRe−Λ=rankΦA1θPA1HθΦH+ΦQ1ΦH−Λ≥K, and the equality holds if and only if Λ=ΦQ1ΦH. On account of K<M, the problem of estimating the noise covariance matrix Q can be formulated as follows:(7)minimizeQrankRe−ΦQ1ΦH,s.t.Re−ΦQ1ΦH≥0,Q1∈Z+,where ≥0 means matrix Re−ΦQ1ΦH is positive semidefinite and Z+ represents a collection of positive definite matrices. Since rank minimization problem is nonconvex and not easy to solve [18], a common method is to use nuclear norm minimization to approximate rank minimization through convex relaxation; since the matrix Re−ΦQ1ΦH is a Hermitian matrix and has the characteristic of positive semidefinite, trace Re−ΦQ1ΦH is able to be used to take the place of Re−ΦQ1ΦH∗. Consequently, the following expression can be obtained:(8)minimizeQtraceRe−ΦQ1ΦH,s.t.Re−ΦQ1ΦH≥0,Q1∈Z+,where trace(∙) represents the trace operation of the matrix which returns the sum of the elements on the main diagonal of the matrix. According to the properties of the matrix trace and the real diagonal characteristics of nonuniform noise, we can easily obtain(9)traceRe−ΦQ1ΦH=traceRe−traceΦQ1ΦH=traceRe−EMσ,where EM∈0,1C×1 denotes a vector, in which only the elements at the ci+1 position are 1 and the others are 0. Based on equations (8) and (9), the minimization problem in equation (7) can be further transformed as follows:(10)maximizeσEMσ,s.t.Re−ΦdiagσΦH≥0,σm2≥0,m=1,2,...,C,where diagσ signifies the diagonal matrix where the diagonal terms are in turn the elements in the vector σ. The optimal estimate of the noise power Qe can be acquired by using the CVX optimization toolbox [19] to solve the above positive semidefinite problem. Then, through the difference of the covariance matrices:(11)Rw=Re−ΦQeΦH,the signal covariance matrix with the noise component removed can be obtained.


### 3.2. Residual Constraint Scheme Based on Covariance Matching Criteria

Since the actual covariance matrix is obtained from the finite snapshot, even if the noise is suppressed in the early stage, the estimation error of the covariance matrix caused by the finite snapshot still exists. To make this problem solved, this paper puts forward a gridless residual constraint scheme that does not depend on the statistical parameters of noise and combines the structural priors of the covariance to achieve the reconstruction of the covariance matrix.Inspired by the prior information of the low-rank and semidefinite Toeplitz structure of the noise-free signal covariance matrixT(u) and based on the covariance fitting criterion CMC2 [20], the following low-rank matrix reconstruction model can be established:(12)minimizeurankTu,Re−1/2Re−ΦRΦHRe−1/2F2≤β2,s.t.ΦTuΦH≥0,where ·F2 represents the square of Frobenius norm of the matrix, and let Δr=vecRe−ΦRΦH, in which vec(∙) is a vectorized operator. It can be obtained from reference [20] that the residuals follow an asymptotic normal distribution, that is,(13)Δr∼AsN0M2,Q,where Q=RT⊗R/L can be replaced by its approximate estimate Qe=ReT⊗Re/L, (∙)T represents transpose operation, AsNμ,σ2 represents an asymptotic normal distribution whose mean and variance are μ and σ2, respectively, and ⊗ refers to the Kronecker inner product of the matrix, and it can be deduced further as follows:(14)Qe−1/2Δr22∼Asχ2M2,where Asχ2M2 stands for asymptotic chi-square distribution with M2 degrees of freedom and ·22 represents the square of Euclidean norm of the matrix. Combining the properties of matrix traces and vectorized operators, we can obtain(15)Re−1/2Re−ΦRΦHRe−1/2F2=trRe−ΦRΦHRe−1/2Re−ΦRΦHRe−1/2=vecHRe−ΦRΦHvecRe−1Re−ΦRΦHRe−1=vecHRe−ΦRΦHRe−T⊗Re−1vecRe−ΦRΦH=Qe−1/2vecRe−ΦRΦH22=Qe−1/2Δr22.LetRw0=R−ΦQeΦH, then Re−R=Rw−Rw0; it can be further acquired that Qe−1/2vecRw−Rw022 also follows the chi-square distribution of M2. As a result, the value of β becomes simple and convenient via MATLAB function chi2inv(1 − η, M2) [7], where η is a very small number, which means that qualification in equation (12) is satisfied with (1 − η) probability. Thus, a specific energy constraint solution can be obtained. Model (12) can be transformed into the following:(16)minimizeurankTu,s.t.Qe−1/2vecE22≤β2,Tu≥0,where E=Rw−ΦTuΦH. Same as the previous section, the trace function is also used here to replace the rank function, and as a consequence, model (16) can be equivalent to(17)minimizeutraceTu,s.t.Qe−1/2vecE22≤β2,Tu≥0.The above model can be solved by SDPT3 in the CVX optimization toolbox; after the optimal solutionue of u is obtained, T(ue) can be used to substitute the estimated value of the noise-free covariance matrix T(u). Since the rank of T(u) is equal to the number of incident sources in theory, an estimation of the number of incident sources can be obtained based on the calculated rank of T(ue). However, errors inevitably exist in the actual system, which cause T(ue) usually full rank matrix; the M-K smaller singular values acquired by the singular value decomposition are not strictly equal to zero. At this time, it is unreasonable to obtain the number of sources based on the number of nonzero singular values, but we can set a threshold to judge(18)K^=m:λ^m≥κ,m∈M,where λ^m denotes the mth singular value of the matrix T(ue) and κ is the threshold. Reference [21] proves that the value of κ is preferably equal to 0.05 λ^max through extensive simulation experiments. Finally, we can possess the estimated values of the arrival angle and power of the incident sources by performing Vandermonde decomposition [8] on T(ue). Besides, since the estimated value of the number of sources is known, subspace-based algorithms (e.g., MUSIC) can also be used to further estimate the direction of the incident signal.


## 3.1. Denoising Method Based on Rank Minimization

For noncoherent sources,rankRe−Q=rankΦA1PA1HΦH=K, where rank(∙) denotes the rank operation of the matrix in brackets. For any diagonal matrix Λ composed of real numbers, rankRe−Λ=rankΦA1θPA1HθΦH+ΦQ1ΦH−Λ≥K, and the equality holds if and only if Λ=ΦQ1ΦH. On account of K<M, the problem of estimating the noise covariance matrix Q can be formulated as follows:(7)minimizeQrankRe−ΦQ1ΦH,s.t.Re−ΦQ1ΦH≥0,Q1∈Z+,where ≥0 means matrix Re−ΦQ1ΦH is positive semidefinite and Z+ represents a collection of positive definite matrices. Since rank minimization problem is nonconvex and not easy to solve [18], a common method is to use nuclear norm minimization to approximate rank minimization through convex relaxation; since the matrix Re−ΦQ1ΦH is a Hermitian matrix and has the characteristic of positive semidefinite, trace Re−ΦQ1ΦH is able to be used to take the place of Re−ΦQ1ΦH∗. Consequently, the following expression can be obtained:(8)minimizeQtraceRe−ΦQ1ΦH,s.t.Re−ΦQ1ΦH≥0,Q1∈Z+,where trace(∙) represents the trace operation of the matrix which returns the sum of the elements on the main diagonal of the matrix. According to the properties of the matrix trace and the real diagonal characteristics of nonuniform noise, we can easily obtain(9)traceRe−ΦQ1ΦH=traceRe−traceΦQ1ΦH=traceRe−EMσ,where EM∈0,1C×1 denotes a vector, in which only the elements at the ci+1 position are 1 and the others are 0. Based on equations (8) and (9), the minimization problem in equation (7) can be further transformed as follows:(10)maximizeσEMσ,s.t.Re−ΦdiagσΦH≥0,σm2≥0,m=1,2,...,C,where diagσ signifies the diagonal matrix where the diagonal terms are in turn the elements in the vector σ. The optimal estimate of the noise power Qe can be acquired by using the CVX optimization toolbox [19] to solve the above positive semidefinite problem. Then, through the difference of the covariance matrices:(11)Rw=Re−ΦQeΦH,the signal covariance matrix with the noise component removed can be obtained.


## 3.2. Residual Constraint Scheme Based on Covariance Matching Criteria

Since the actual covariance matrix is obtained from the finite snapshot, even if the noise is suppressed in the early stage, the estimation error of the covariance matrix caused by the finite snapshot still exists. To make this problem solved, this paper puts forward a gridless residual constraint scheme that does not depend on the statistical parameters of noise and combines the structural priors of the covariance to achieve the reconstruction of the covariance matrix.Inspired by the prior information of the low-rank and semidefinite Toeplitz structure of the noise-free signal covariance matrixT(u) and based on the covariance fitting criterion CMC2 [20], the following low-rank matrix reconstruction model can be established:(12)minimizeurankTu,Re−1/2Re−ΦRΦHRe−1/2F2≤β2,s.t.ΦTuΦH≥0,where ·F2 represents the square of Frobenius norm of the matrix, and let Δr=vecRe−ΦRΦH, in which vec(∙) is a vectorized operator. It can be obtained from reference [20] that the residuals follow an asymptotic normal distribution, that is,(13)Δr∼AsN0M2,Q,where Q=RT⊗R/L can be replaced by its approximate estimate Qe=ReT⊗Re/L, (∙)T represents transpose operation, AsNμ,σ2 represents an asymptotic normal distribution whose mean and variance are μ and σ2, respectively, and ⊗ refers to the Kronecker inner product of the matrix, and it can be deduced further as follows:(14)Qe−1/2Δr22∼Asχ2M2,where Asχ2M2 stands for asymptotic chi-square distribution with M2 degrees of freedom and ·22 represents the square of Euclidean norm of the matrix. Combining the properties of matrix traces and vectorized operators, we can obtain(15)Re−1/2Re−ΦRΦHRe−1/2F2=trRe−ΦRΦHRe−1/2Re−ΦRΦHRe−1/2=vecHRe−ΦRΦHvecRe−1Re−ΦRΦHRe−1=vecHRe−ΦRΦHRe−T⊗Re−1vecRe−ΦRΦH=Qe−1/2vecRe−ΦRΦH22=Qe−1/2Δr22.LetRw0=R−ΦQeΦH, then Re−R=Rw−Rw0; it can be further acquired that Qe−1/2vecRw−Rw022 also follows the chi-square distribution of M2. As a result, the value of β becomes simple and convenient via MATLAB function chi2inv(1 − η, M2) [7], where η is a very small number, which means that qualification in equation (12) is satisfied with (1 − η) probability. Thus, a specific energy constraint solution can be obtained. Model (12) can be transformed into the following:(16)minimizeurankTu,s.t.Qe−1/2vecE22≤β2,Tu≥0,where E=Rw−ΦTuΦH. Same as the previous section, the trace function is also used here to replace the rank function, and as a consequence, model (16) can be equivalent to(17)minimizeutraceTu,s.t.Qe−1/2vecE22≤β2,Tu≥0.The above model can be solved by SDPT3 in the CVX optimization toolbox; after the optimal solutionue of u is obtained, T(ue) can be used to substitute the estimated value of the noise-free covariance matrix T(u). Since the rank of T(u) is equal to the number of incident sources in theory, an estimation of the number of incident sources can be obtained based on the calculated rank of T(ue). However, errors inevitably exist in the actual system, which cause T(ue) usually full rank matrix; the M-K smaller singular values acquired by the singular value decomposition are not strictly equal to zero. At this time, it is unreasonable to obtain the number of sources based on the number of nonzero singular values, but we can set a threshold to judge(18)K^=m:λ^m≥κ,m∈M,where λ^m denotes the mth singular value of the matrix T(ue) and κ is the threshold. Reference [21] proves that the value of κ is preferably equal to 0.05 λ^max through extensive simulation experiments. Finally, we can possess the estimated values of the arrival angle and power of the incident sources by performing Vandermonde decomposition [8] on T(ue). Besides, since the estimated value of the number of sources is known, subspace-based algorithms (e.g., MUSIC) can also be used to further estimate the direction of the incident signal.


## 4. Algorithm Analysis

### 4.1. Estimated Accuracy and Resolution

The proposed gridless DOA estimation algorithm makes use of the diagonal characteristics of the nonuniform noise covariance matrix and the Toeplitz and Hermitian properties of the noise-free signal covariance matrix to effectively improve the estimation accuracy of the algorithm and obtain the number of sources as a byproduct. The array aperture is expanded by using MRA, and the resolution is further improved compared with the conventional uniform linear array.


### 4.2. Identifiability

It can be seen from Section3.1 that the proposed algorithm is mainly based on the knowledge that the rank of the correlation matrix of incoherent signals is equal to the number of sources in the stage of removing nonuniform noise; therefore, this algorithm is no longer applicable when the sources are coherent. In addition, the number of signal sources should be as small as possible than the number of sensors to ensure the accuracy of the noise power estimation. Empirically, the number of sources that can be identified by the algorithm is M/2.


### 4.3. Computational Complexity

This article mainly measures the computational complexity of the algorithm by comparing the running time of the CPU. When the proposed algorithm is solved with CVX, in the scenario of experiment 6, the number of iterations is basically about 13, while IL1-SRACV needs about 28 iterations under the same error tolerance. Compared with SPA, the proposed algorithm requires a longer running time due to the two-step operation, but in return, the estimation accuracy is higher.All the experimental processes were completed on a desktop computer with Win10, 3.6 GHz, 8-core processor, and the MATLAB software version used was R2019b.


### 4.4. Deterministic CRBs

As we all know, the Cramer–Rao bound (CRB) provides a lowest lower bound for arbitrary unbiased parameter estimation [22, 23]. The CRB derivation used in the subsequent simulation is mainly based on reference [24], and the specific formula can be written as follows:(19)CRB=σ22LReD˜HPA˜⊥D˜⊙P^T−1,where A˜=Q−1/2ΦA1θ, D=da1θ/θθ=θ1,da1θ/θθ=θ2,…,da1θ/θθ=θk, D˜=Q−1/2ΦD, P^=1/N∑i=1NsisiH, PA˜⊥ represents the orthogonal projection of matrix A˜ on P, and ⊙ means the Hadamard product of a matrix.


## 4.1. Estimated Accuracy and Resolution

The proposed gridless DOA estimation algorithm makes use of the diagonal characteristics of the nonuniform noise covariance matrix and the Toeplitz and Hermitian properties of the noise-free signal covariance matrix to effectively improve the estimation accuracy of the algorithm and obtain the number of sources as a byproduct. The array aperture is expanded by using MRA, and the resolution is further improved compared with the conventional uniform linear array.


## 4.2. Identifiability

It can be seen from Section3.1 that the proposed algorithm is mainly based on the knowledge that the rank of the correlation matrix of incoherent signals is equal to the number of sources in the stage of removing nonuniform noise; therefore, this algorithm is no longer applicable when the sources are coherent. In addition, the number of signal sources should be as small as possible than the number of sensors to ensure the accuracy of the noise power estimation. Empirically, the number of sources that can be identified by the algorithm is M/2.


## 4.3. Computational Complexity

This article mainly measures the computational complexity of the algorithm by comparing the running time of the CPU. When the proposed algorithm is solved with CVX, in the scenario of experiment 6, the number of iterations is basically about 13, while IL1-SRACV needs about 28 iterations under the same error tolerance. Compared with SPA, the proposed algorithm requires a longer running time due to the two-step operation, but in return, the estimation accuracy is higher.All the experimental processes were completed on a desktop computer with Win10, 3.6 GHz, 8-core processor, and the MATLAB software version used was R2019b.


## 4.4. Deterministic CRBs

As we all know, the Cramer–Rao bound (CRB) provides a lowest lower bound for arbitrary unbiased parameter estimation [22, 23]. The CRB derivation used in the subsequent simulation is mainly based on reference [24], and the specific formula can be written as follows:(19)CRB=σ22LReD˜HPA˜⊥D˜⊙P^T−1,where A˜=Q−1/2ΦA1θ, D=da1θ/θθ=θ1,da1θ/θθ=θ2,…,da1θ/θθ=θk, D˜=Q−1/2ΦD, P^=1/N∑i=1NsisiH, PA˜⊥ represents the orthogonal projection of matrix A˜ on P, and ⊙ means the Hadamard product of a matrix.


## 5. Numerical Simulations

In this part, we conducted simulation experiments to test the performance of the proposed two-stage method and compared it with subspace-based method RTM, sparse representation-based method IL1-SRACV, and gridless method SPA. To be fair, we extend the algorithms listed above to MRA. Unless other specified, we assume that the incident sources are independent with the same powerp, the array is the MRA with 8 elements, and the relative position of MRA is D=0,1,4,10,16,18,21,23×d. The covariance matrix of nonuniform noise satisfies(20)Q=diag10.0,2.0,8.5,10.5,5.0,10.0,10.5,3.0.The SNR is defined as follows [22]:(21)SNR=10log10pM×∑m=1M1/qmdB.


### 5.1. Experiment 1

We assume four independent sources from the direction ofθ=−5.0°,0°,5.0°,10.0°. The snapshot number L = 300, SNR = 0 dB, and η = 0.001. Figure 2 presents the spatial spectrum of the proposed algorithm, SPA, and RTM. It is easy to see that the proposed method has a relatively sharp peak and is easier to distinguish under the same conditions.Figure 2
Spatial spectrum of the proposed method (M = 8, SNR = 0 dB, L = 300, and θ=−5.0°,0°,5.0°,10.0°).


### 5.2. Experiment 2

Let two independent sources are relatively close, whereθ=−1.0°,1.0°; other experimental conditions are consistent with experiment 1. Figure 3 shows the spatial spectrum of the proposed algorithm in ULA and MRA composed of 8 array elements. As can be seen from Figure 3, the spatial spectrum of the proposed method in MRA can better distinguish two sources that are closer together. In other words, the adoption of MRA can effectively improve the resolution.Figure 3
Spatial spectrum of the proposed method (M = 8, SNR = 0 dB, L = 300, and θ=−1.0°,1.0°).


### 5.3. Experiment 3

Letθ=−3.41°,6.37° and L = 500; the change of the SNR is increased from −10 dB in steps of 2 dB to 10 dB. 300 Monte Carlo experiments were conducted to reduce the contingency of the experiment. The estimation accuracy of the algorithm is mainly measured by root mean square error (RMSE), which can be expressed as follows:(22)RMSE=∑l=1L∑k=1Kθ^kl−θk2L×K.Figure4 shows the plot of estimation error and success probability with SNR, where Figure 4(a) shows the change curve of RMSE and Figure 4(b) shows the change curve of success probability. Here, the success probability refers to the ratio of the number of experiments with the error between the estimated angles and the real angles not more than 0.2∘ and the total number of Monte Carlo experiments. It can be seen from Figure 4(a) that the RMSEs of all methods decrease with the increase of SNR. Compared with other algorithms, the proposed method has the smallest RMSE and gradually fits to the CRB with the increase of SNR. As can be seen from Figure 4(b), the proposed method can guarantee the estimation error within 0.2∘ at SNR=−6dB. From the above results, we can conclude that the proposed method has obvious advantages over other methods in estimation accuracy.Figure 4
The RMSE and the success probability of different methods vary with SNR (M = 8 and L = 500). (a) The change curve of RMSE. (b) The change curve of success probability.
(a)(b)


### 5.4. Experiment 4

We set the SNR to 0 dB, and the number of snapshots increased from 100 to 1000 with an interval of 100, andθ=−3.41°,5.37°. 300 independent Monte Carlo experiments were also conducted. Figure 5 shows the change curve of RMSE with the number of snapshots. It can be summarized from the figure that the estimation error of the proposed method is smaller than that of other algorithms overall. With the increase in the number of snapshots, the RMSE of the proposed algorithm is closer to the CRB low bound.Figure 5
The RMSE of different methods varies with snapshot (M = 8 and SNR = 0 dB).


### 5.5. Experiment 5

LetL = 500; other experimental conditions are the same as experiment 4, and the change of the SNR is increased from −10 dB in steps of 2 dB to 10 dB. Figure 6 shows the success rate of the source estimation of the proposed algorithm compared with that of Akaike’s information criterion (AIC), minimum description length (MDL) [25], and Gershgorin disk estimator (GDE) [26]. The success rate is defined as the ratio of the number of accurate estimations of all sources to the total number of Monte Carlo experiments. It can be summarized from Figure 6 that the proposed method has better estimation accuracy than other algorithms under low SNR and has a higher success probability overall.Figure 6
The success probability of sources varies with SNR (M = 8 and L = 500).


### 5.6. Experiment 6

Letθ=−4.78°,2.56° and L = 400; other experimental conditions are the same as experiment 3. Figure 7 shows the curve of the average running time of different algorithms versus SNR. It can be seen from the figure that RTM has low computational complexity and takes the shortest time. The running time of the proposed algorithm is much lower than that of IL1-SRACV and slightly higher than that of SPA. The average running time of the algorithms increases slightly with the increase of the SNR.Figure 7
The average running time varies with SNR (M = 8 and L = 400).


## 5.1. Experiment 1

We assume four independent sources from the direction ofθ=−5.0°,0°,5.0°,10.0°. The snapshot number L = 300, SNR = 0 dB, and η = 0.001. Figure 2 presents the spatial spectrum of the proposed algorithm, SPA, and RTM. It is easy to see that the proposed method has a relatively sharp peak and is easier to distinguish under the same conditions.Figure 2
Spatial spectrum of the proposed method (M = 8, SNR = 0 dB, L = 300, and θ=−5.0°,0°,5.0°,10.0°).


## 5.2. Experiment 2

Let two independent sources are relatively close, whereθ=−1.0°,1.0°; other experimental conditions are consistent with experiment 1. Figure 3 shows the spatial spectrum of the proposed algorithm in ULA and MRA composed of 8 array elements. As can be seen from Figure 3, the spatial spectrum of the proposed method in MRA can better distinguish two sources that are closer together. In other words, the adoption of MRA can effectively improve the resolution.Figure 3
Spatial spectrum of the proposed method (M = 8, SNR = 0 dB, L = 300, and θ=−1.0°,1.0°).


## 5.3. Experiment 3

Letθ=−3.41°,6.37° and L = 500; the change of the SNR is increased from −10 dB in steps of 2 dB to 10 dB. 300 Monte Carlo experiments were conducted to reduce the contingency of the experiment. The estimation accuracy of the algorithm is mainly measured by root mean square error (RMSE), which can be expressed as follows:(22)RMSE=∑l=1L∑k=1Kθ^kl−θk2L×K.Figure4 shows the plot of estimation error and success probability with SNR, where Figure 4(a) shows the change curve of RMSE and Figure 4(b) shows the change curve of success probability. Here, the success probability refers to the ratio of the number of experiments with the error between the estimated angles and the real angles not more than 0.2∘ and the total number of Monte Carlo experiments. It can be seen from Figure 4(a) that the RMSEs of all methods decrease with the increase of SNR. Compared with other algorithms, the proposed method has the smallest RMSE and gradually fits to the CRB with the increase of SNR. As can be seen from Figure 4(b), the proposed method can guarantee the estimation error within 0.2∘ at SNR=−6dB. From the above results, we can conclude that the proposed method has obvious advantages over other methods in estimation accuracy.Figure 4
The RMSE and the success probability of different methods vary with SNR (M = 8 and L = 500). (a) The change curve of RMSE. (b) The change curve of success probability.
(a)(b)


## 5.4. Experiment 4

We set the SNR to 0 dB, and the number of snapshots increased from 100 to 1000 with an interval of 100, andθ=−3.41°,5.37°. 300 independent Monte Carlo experiments were also conducted. Figure 5 shows the change curve of RMSE with the number of snapshots. It can be summarized from the figure that the estimation error of the proposed method is smaller than that of other algorithms overall. With the increase in the number of snapshots, the RMSE of the proposed algorithm is closer to the CRB low bound.Figure 5
The RMSE of different methods varies with snapshot (M = 8 and SNR = 0 dB).


## 5.5. Experiment 5

LetL = 500; other experimental conditions are the same as experiment 4, and the change of the SNR is increased from −10 dB in steps of 2 dB to 10 dB. Figure 6 shows the success rate of the source estimation of the proposed algorithm compared with that of Akaike’s information criterion (AIC), minimum description length (MDL) [25], and Gershgorin disk estimator (GDE) [26]. The success rate is defined as the ratio of the number of accurate estimations of all sources to the total number of Monte Carlo experiments. It can be summarized from Figure 6 that the proposed method has better estimation accuracy than other algorithms under low SNR and has a higher success probability overall.Figure 6
The success probability of sources varies with SNR (M = 8 and L = 500).


## 5.6. Experiment 6

Letθ=−4.78°,2.56° and L = 400; other experimental conditions are the same as experiment 3. Figure 7 shows the curve of the average running time of different algorithms versus SNR. It can be seen from the figure that RTM has low computational complexity and takes the shortest time. The running time of the proposed algorithm is much lower than that of IL1-SRACV and slightly higher than that of SPA. The average running time of the algorithms increases slightly with the increase of the SNR.Figure 7
The average running time varies with SNR (M = 8 and L = 400).


## 6. Conclusions

This paper proposes a two-stage gridless DOA estimation algorithm for MRA in nonuniform noise. The proposed method can not only effectively eliminate the effect of nonuniform noise but also further improve the estimation accuracy and resolution. Moreover, the number of signals is also able to be acquired as a byproduct. Simulation experiments prove the performance of the proposed method.



---

*Source: 1580391-2020-10-24.xml*