# Reweighted Covariance Fitting Based on Nonconvex Schatten-p Minimization for Gridless Direction of Arrival Estimation

Authors: Lama Zien Alabideen; Oumayma Al-Dakkak; Khaldoun Khorzom
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this paper, we reformulate the gridless direction of arrival (DoA) estimation problem in a novel reweighted covariance fitting (CF) method. The proposed method promotes joint sparsity among different snapshots by means of nonconvex Schatten-p quasi-norm penalty. Furthermore, for more tractable and scalable optimization problem, we apply the unified surrogate for Schatten-p quasi-norm with two-factor matrix norms. Then, a locally convergent iterative reweighted minimization method is derived and solved efficiently via a semidefinite program using the optimization toolbox. Finally, numerical simulations are carried out in the background of unknown nonuniform noise and under the consideration of coprime array (CPA) structure. The results illustrate the superiority of the proposed method in terms of resolution, robustness against nonuniform noise, and correlations of sources, in addition to its applicability in a limited number of snapshots.


---

## Body

## 1. Introduction

DoA estimation is a major task in array signal processing. It has wide applications in radar, sonar, wireless communications, underwater acoustics, and seismology [1–3]. The main goal of DoA estimation is to retrieve the direction information of several sources from the outputs of the antenna array. The study of DoA estimation has a long history. One of an outstanding class of methods is designated as subspace-based methods such as root multiple signal classification (RootMUSIC) [4] and the estimation of parameters by rotational invariant techniques (ESPRIT) [5], and their variants. In general, these methods have high estimation accuracy with low computational complexity, but they suffer from considerable limitations. For example, they require a priori knowledge of the number of sources that may be difficult to acquire in practice. Moreover, they need a sufficient number of data snapshots to estimate the data covariance matrix precisely. In addition, they are sensitive to source correlations that may lead to a rank deficiency in the sample data covariance matrix.With the development of sparse representation and compressed sensing (CS) [6, 7], the last decade has received remarkable advances resulting in a number of sparse methods for DoA estimation. Sparse methods can be applied in several challenging scenarios, including cases with no knowledge of the number of sources, a limited number of snapshots, and highly correlated sources. The early grid-based sparse methods have been focused on discrete linear systems, meaning that the DoAs have to lie on a set of predefined grid points so that the collected snapshots can be sparsely represented under a finite discrete dictionary. However, in reality, the DoAs are continuously valued and nonlinear with the measured data; hence, there is no guarantee that the true DoAs will be on the grid. This results in a grid or basis mismatch problem [8, 9]. To address this issue, several off-grid methods have been proposed [10–12]. However, these methods bring complexity issue while the grid selection remains a major problem since gridding still exists. Recently, gridless sparse methods have been developed [13, 14]. These methods eliminate the requirement of gridding by working directly in a continuous infinite dictionary and recover the DoAs by solving a semidefinite programming (SDP) problem.


## 2. Related Works

The gridless sparse methods in DoA estimation are categorized into two methodologies. The first one is based on the concept of atomic norm minimization (ANM) [15]. While, the second methodology makes use of sparse representation of second-order statistics of the array output, which is designated as gridless CF.ANM has been successfully applied in estimating a spectrally-sparse signal without any statistical assumptions in [12, 16]. It has been shown that the frequencies can be exactly estimated if they are sufficiently separated by at least 4/n, where n denotes the size of full data. Later, this condition is reduced to 2.52/n in [17]. Initially, the atomic norm soft-thresholding (AST) method was introduced in [18] for full data case in the presence of stochastic noise. Then, the study was extended to the case of compressive data for a single measurement vector (SMV), where the problem was referred to as continuous CS (CCS) in [12]. Next, CCS was expanded to the case of multiple measurement vectors (MMVs) in [19]. In order to alleviate the limitation of estimation resolution due to the effect of separation condition, the reweighted atomic-norm minimization (RAM) was proposed in [20] with applicability to DoA estimation.The gridless CF technique in DoA estimation is not new. The covariance matching estimation technique (COMET) in [21] and the asymptotic maximum likelihood (AML) technique in [22] are both consistent gridless DoA estimators as the signal-to-noise ratio (SNR) increases. In addition, the continuous sparse recovery (CSR) method [23] and the nuclear norm minimization method followed by MUSIC (NNM-MUSIC) [13], do not require the number of sources to estimate the DoAs precisely, while the knowledge of error bounds is essential. The gridless sparse iterative covariance-based estimation (GL-SPICE or GLS) method [24] (also referred to as the sparse and parametric approach (SPA) [25]) has many merits. Mainly, GLS is a hyperparameter-free method that can estimate the noise variance(s) consistently, also it is applicable to an arbitrary number of snapshots and robust against moderate sources correlation. However, in the case of finite snapshots, GLS suffers from the separation condition due to the theoretical relation with ANM. Recently, the covariance matrix reconstruction approach (CMRA) was developed in [26] as a super-resolution method. Nevertheless, CMRA faces difficulties in highly correlated sources scenarios. Furthermore, the study in [27] has proposed an improved version of CMRA (ICMRA) by applying a family of nonconvex penalties on the singular values of the covariance matrix. ICMRA has shown an enhancement in sparsity and resolution compared to CMRA. On the contrary, ICMRA suffers from unsatisfactory performance in the presence of low SNR and low number of snapshots. A clearer comparison of the abovementioned methods is provided in Table 1.Table 1
Gridless DoA estimation methods and their comparisons.
Gridless sparse methodology
Methods
Specifications
Atomic norm
Covariance fitting
AST
[12]
RAM
MRAM
SPA/GLS
CMRA
ICMRA
AML
COMET
CSR
NNM-MUSIC
Need for error bound
√
√
√
√
√
√
Need for model order
√
√
√
Angle separation condition
√
√
√
√
√
√
√
√
Sparsity promoting
√
√
√
√
√
√
√
√
Noise estimation
√
√
√
√
Feasibility to low number of snapshots
√
√
√
√
√
Immunity to sources correlation
√
√
√
√
√
High detection probability in low SNR
√
√
√
√Concerning noise, many sparse methods assume that the antenna array noises are spatially uniform white with equal variances. This assumption helps in choosing a suitable regularization parameter in a sparse representation optimization formula. In several practical applications, this may be violated since the noise levels among all antennas exhibit some perturbations due to the imperfection of the array calibration, the nonidealities of the communication channels, or the mutual coupling between antennas. In many sparse methods, the antennas are sparsely placed so their noises are spatially uncorrelated and modeled as a diagonal covariance matrix [28, 29], but the antennas noise variances are no longer identical (i.e., nonuniform noise). Such a noise model is relevant in situations with CF-based methods [14], while ANM-based methods are favored when the noise is uniform with bounded power or when the sources are highly correlated. In fact, the nonuniform noise is a challenging problem in sparse estimation methods, where the sample covariance matrix is no longer low-rank. To eliminate this problem the recent research studies get rid of noise in two major methods, the first one is by using linear transform primarily and then solving noise-free optimization problem as in [30, 31]. However, this method needs to design a linear transform able to separate the noise from the signal, which is pointless in the general cases, specifically when the transformation requires prior knowledge of the noise parameters. The second method treats the noise variance as a variable to be optimized in the optimization problem [14, 25]. In this paper, the latter method is adopted.Despite the convexity of the preliminary gridless optimization models, it has been shown that the estimation performance of such convex relaxation will degrade in the presence of measurement noise or when some assumptions are violated (e.g., spatial sources separation limit and signals correlation). Hence, the solution may significantly deviate from the original (i.e., rank minimization). In order to achieve a better approximation of the original solution, many nonconvex surrogate functions have been proposed, includingLp norm, Logarithm, Laplace, and Schatten-p norm which are well-studied in the literature [32, 33]. In [27], the DoA problem has been formulated in a nonconvex minimization problem using a family of nonconvex penalties on the singular values of the covariance matrix resulting in an enhancement in both sparsity and resolution. Recently, we have proposed a high-resolution gridless method for noiseless joint sparse frequency estimation named modified reweighted atomic norm minimization (MRAM) [34]. MRAM was based on nonconvex relaxation of ANM by Schatten-p norm. This approximation has offered improved performance compared to trace norm relaxation. Hence, due to the one-to-one relation between frequency and direction information, we look forward to achieve a better gridless DoA estimation using Schatten-p quasi-norm penalty.


## 3. Main Contributions

Since, gridless CF is preferable in the presence of a nonuniform noise background, in this paper we develop a new gridless sparse method for DoA estimation based on recently developed gridless CF methods. By exploiting the positive semidefinite (PSD) and Toeplitz structure of the parameterized data covariance matrix, a gridless CF criterion is formulated in an efficient joint sparsity optimization problem. More specifically, we propose a novel objective function based on introducing the nonconvex Schatten-p quasi-norm as a sparsity metric. Then, for a more tractable and scalable optimization problem, the unified surrogate for Schatten-p quasi-norm with two-factor matrix norms is applied. Consequently, an iterative reweighted minimization method is derived with a SDP formulation. Regarding the major role of the array geometry and in order to acquire a high degree of freedom with less mutual coupling and a reasonable number of physical antennas, the coprime array (CPA) structure [23, 35] is adopted. Then, the new method has been tested, where the numerical results have shown an improved performance as compared to the state-of-art methods.Notations 1.
Throughout this paper, the bold letters are reserved for matrices and vectors. The transpose is denoted by⋅T, and the complex conjugate or Hermitian is denoted by ⋅H. tr⋅ and rank⋅ represent matrix trace and rank, respectively. E⌈⌉ stands for the expectation. ⋅2, ⋅∗, and ⋅F refer to l2, nuclear, and Forbenius norms, respectively. A≥0 means that A is a positive semidefinite matrix.


## 4. Problem Formulation

Assume there areK uncorrelated narrow-band far-field source signals sk, impinge from directions θ=θ1,θ2,…,θKT on a sparse linear CPA array. In particular, for a coprime number pair N,M with N>M. CPA consists of two uniform linear subarrays, the first consists of M elements having interelement spacing of Nd, whereas the other consists of N elements having interelement spacing of Md, where d is chosen to be half of the wavelength λ/2. According to the properties of coprime integers, the elements of the subarrays do not overlap with each other when they are aligned except the first element. Therefore, CPA consists of M+N−1 physical antennas in total. The CPA output data matrix can be considered as the output of virtual N×M elements uniform linear array (ULA) by retaining the antennas indexed by Ω=r1,r2,…,rM+N−1⊂1,…,NM. The observation model of CPA at t snapshot can be modeled as follows:(1)yΩt=∑k=1KaΩθkskt+nΩt=AΩθst+nΩt,t=1,…,L,where L is the number of snapshots that are jointly sparse, in the sense that they share the same support where the DoAs are encoded, and yΩt∈ℂM+N−1×1, st∈ℂK×1, and nΩt∈ℂM+N−1×1 are the array output, the source signal vector, and the noise powers vector, respectively. AΩθ=aΩθ1,…,aΩθK∈ℂM+N−1×K denotes the steering matrix, where aΩθk is the steering vector of the kth source which is determined by the geometry of the antenna array and is given by(2)aΩθk=ei2πr1d/λsinθk,ei2πr2d/λsinθk,…,ei2πrM+N−1d/λsinθkT,where ri is the physical positions of the ith CPA element. Furthermore, it is assumed that s1K are spatially and temporarily uncorrelated with independently and uniformly distributed phases. It follows that EssH=diagP, where P=P1,P2,…,PK∈ℝ+1×K is the vector of the powers of sources.Moreover, suppose the noisenΩ1M+N−1 is independent of s1K and satisfies EnΩnΩH=diagσΩ, where σΩ∈ℝ+M+N−1×1 denotes the noise variances vector whose elements are assumed to be unequal in this study. Briefly, the array output for L snapshots can be modeled as(3)YΩ=yΩ1,…,yΩL=AΩθS+NΩ∈ℂM+N−1×L,where S=s1,…,sL∈ℂK×L is the sources signals matrix and NΩ=nΩ1,…,nΩL∈ℂM+N−1×L is the measurement noise matrix.Concerning the main goal of DoA estimation, in this paper, we seek to estimate the unknown parametersθ,P,σΩ given the measurement matrix YΩ.


## 5. Gridless CF with Nonconvex Schatten-p Minimization

In this section, we propose a new gridless method to estimate the parametersθ,P,σΩ by introducing the nonconvex Schatten-p norm regularization into the same CF criterion as in Sparse Iterative Covariance-based Estimation (SPICE) method [36] and its gridless version method (GL-SPICE or GLS) [24, 25]. Mainly, the proposed method inherits the merits of GLS, with enhanced performance as it is verified by the numerical simulations in Section 6.


### 5.1. Gridless CF

The key idea of gridless CF is exploiting the structure of the data covariance matrix and performing the estimation under ungrided parameterized domain with specific statistical considerations (i.e., uncorrelated sources and reasonable number of snapshots). Based on the assumptions in Section4, the uncorrelated data yΩ have the data covariance matrix RΩ∈ℂM+N−1×M+N−1 given by(4)RΩ=EyΩtyΩHt=AΩθdiagPAΩHθ+diagσΩ.Note thatRΩ−1 exists in the presence of noise for σΩ>0. RΩ is a highly nonlinear function of the DoAs parameters. To linearize the problem, a reparameterization of RΩ using a PSD Toeplitz matrix Tu is performed. Hence, according to Vandermonde decomposition lemma (see, e.g., [3]), RΩ can be linearly reparametrized as(5)RΩu,σΩ=Tu,Ω+diagσΩ,Tu≥0,σΩ>0,where Tu,Ω=ΓΩTuΓΩT=AΩθdiagPAΩHθ for Tu∈ℂNM×NM is a PSD Toeplitz matrix formed by the vector u∈ℂ1×NM. ΓΩ∈0,1M+N−1×NM denotes a selection matrix, such that the ith row of ΓΩ contains all zeros but a single 1 at the ri position. The importance of the parameterization in equation (5) is that Tu is linear in the vector u, which comprises the interesting parameters.Recall, the sample covariance matrix is set asR˜Ω=1/LYΩYΩH. In this paper, we are mainly interested in the case, where L≥M+N−1 (i.e., R˜Ω−1 exists with probability one) and K≤NM−1 (i.e., Tu is rank-deficient, then θ,P can be uniquely determined [25]). As well, we consider the same CF criterion as in the SPICE method, which is given by(6)fu,σΩ=RΩ−1/2R˜Ω−RΩR˜Ω−1/2F2=trRΩ−1R˜Ω+trR˜Ω−1RΩ−2NM.By using the linear parameterization in equation (5), the CF can be accomplished by minimizing f leading to the following equivalences:(7)minu,σ≥0trR˜Ω1/2Tu,Ω+diagσΩ−1R˜Ω1/2+trR˜Ω−1Tu,Ω+trR˜Ω−1diagσΩsubject toTu≥0,(8)minX,u,σ≥0trX+trR˜Ω−1Tu,Ω+trR˜Ω−1diagσΩsubject toTu≥0,X≥R˜Ω1/2Tu,Ω+diagσΩ−1R˜Ω1/2,(9)minX,u,σ≥0trX+trR˜Ω−1Tu,Ω+trR˜Ω−1diagσΩsubject toXR˜Ω−1/2R˜Ω−1/2Tu,Ω+diagσΩ≥0.Hence, iff minimization is casted as SDP, then it is convex. Problem (9) was derived under the development of the GLS method. In particular, the GLS method transforms the direction estimation problem into the estimation of a PSD Toeplitz matrix in which the parameters of the directions are encoded. Therefore, once the SDP in problem (9) is solved, the parameters of interest can be obtained from the estimated covariance matrix by either using the conventional subspace methods or applying the Vandermonde decomposition lemma.In spite of the strong theoretical guarantees of the GLS method, it suffers from considerable limitations according to Table1. In fact, we believe that its performance can be further improved by employing the nonconvex penalty, which to the best of our knowledge, has not been studied in the literature. Note that the nonconvex Schatten-p quasi-norm regularization is able to achieve a more accurate and sparser solution than the trace norm. In this paper, we reformulate the minimization problem of f in a new objective function by introducing Schatten-p quasi-norm. Then, the resultant optimization problem is modeled as a SDP formula.


### 5.2. Schatten-p Norm

Let⋅P be the Schatten-p norm of the matrix M∈ℂn×m, which is defined as the lp norm of its singular values σkM, i.e.,(10)Mp=∑k=1minn,mσkpM1/p=trMTMp/21/p,0<p<∞.Recently, Schatten-p norm has received significant attention from researchers in various domains [27, 33], due to its outstanding properties that make it a good choice in many optimization problems (e.g., rank function relaxation). Specifically, Mp equals to the nuclear norm when p=1. As well, if 0<p<1, then the quasi-norm (for 0<p<1, Mpis a quasi-norm because it does not satisfy the triangle inequality)Mp tends to rankM. Moreover, Zhang et al. [37] verified theoretically that Mp with small p requires significantly fewer measurements, with superior performance over the conventional nuclear norm in the majority of the problems.


### 5.3. The Proposed Method

Motivated by the properties of Schatten-p norm and GLS method, we develop a new method by replacing thetrR˜Ω−1Tu,Ω term in problem (8) by its Schatten-p quasi-norm, leading to the following optimization problem:(11)minX,u,σ≥0trX+λZupp+trR˜Ω−1diagσΩ,0<p<1subject toTu≥0,X≥R˜Ω1/2Tu,Ω+diagσΩ1R˜Ω1/2,where λ>0 is a regularization parameter and Zu=R˜Ω−1Tu,Ω is a PSD matrix clearly:(12)Zupp=∑k=1M+N−1σkpZu,0<p<1.According to equation (10), Zupp is given aswhere σkZu≥01M+N−1 are descendingly sorted singular values of Zu.Since0<p<1, then for any p,p1,p2>0 satisfying(13)1p=1p1+1p2,the bi-Schatten-p norm surrogate can be applied successfully [38]. In other words, the Schatten-p quasi-norm of Zu matrix is equivalent to minimizing the product of the Schatten-p1 norm (or quasi-norm) and Schatten-p2 norm (or quasi-norm) of its two factor matrices. So, as an alternative to the nonconvex quasi-norm in (12), the following efficient factorization is used:(14)1pZupp=minR˜Ω−1,Tu,ΩZ=R˜Ω−1Tu,Ω1p1R˜Ω−1p1p1+1p2Tu,Ωp2p2.Notice thatR˜Ω−1 is fixed, so whatever the value of p1 is on minimizing Problem (14), one only needs to minimize 1/p2Tu,Ωp2p2. Substituting the spectral regularization in Problem (11) with Problem (14) and setting λ=1/p, we infer the following problem:(15)minX,u,σ≥0trX+1p2Tu,Ωp2p2+trR˜Ω−1diagσΩsubject toTu≥0,X≥R˜Ω1/2Tu,Ω+diagσΩ−1R˜Ω1/2.With respect to equation (13), if p2≥1, Problem (15) is strongly convex; thus, it can be easily solved by off-the-shelf algorithms. When p2<1, then Problem (15) becomes nonconvex and must be relaxed. However, we empirically found that the later problem has superior properties in terms of estimation performance and solution sparsity. Therefore, in the following, we consider only the nonconvex model and work on the relaxation issue of Problem (15).Let’s setp2=p for simplicity. Since the proposed objective function in Problem (15) is nonconvex, then it must be linearized. A popular locally convergent approach to the minimization of such a problem is the majorization maximization (MM) method [39]. Accordingly, at the lth iteration, the Schatten-p quasi-norm term will be replaced by its tangent plane, which is given as(16)σkTu,Ωl+εp+pσkTu,Ωl+ε1−pσkTu,Ω−σkTu,Ωl,where ε >0 is a smoothing parameter. Notice that Tu,Ωl is constant. Consequently, by applying equation (16) into Problem (15), the resultant optimization formula will be rewritten as(17)minX,u,σ≥0trX+1pWlTu,Ω⋆+trR˜Ω−1diagσΩsubject toTu≥0,X≥R˜Ω1/2Tu,Ω+diagσΩ−1R˜Ω1/2.In Problem (17), the Schatten-p quasi-norm term is reduced to the weighted nuclear norm (WNN), where WlTu,Ω⋆=trWlTu,Ω=∑k=1M+N−1wklσkTu,Ω and Wl=D⋅diagw1l,w2l,…,wN+M−1l⋅DH is the reweighting matrix in the lth iteration, whose elements are wkl=p/σkTu,Ωl+ε1−p. Therefore, we name our method as reweighted p-GLS (Rwp-GLS). The matrix Tu,Ωl is updated based on the latest solution of the problem, and D is obtained by using the eigen decomposition of Tu,Ωl so that DDH=I.SinceTu,Ωpp is a concave function when 0<p<1, then at each iteration its value decreases significantly, faster than the decrease of its tangent plane. It follows that, by iteratively solving Problem (17), the objective function monotonically decreases and converges to a local minimum. Following the same steps in deriving Problem (9), at each iteration of Rwp-GLS, we need to solve the following SDP:(18)minX,u,σ≥0trX+1pWlTu,Ω⋆+trR˜Ω−1diagσΩsubject toconstraint of9.By solving the SDP in Problem (18) by off-the-shelf SDP solvers, e.g., SDPT3 [40], R^Ω=Tu^+diagσ^Ω is obtained given the solution u^,σ^Ω. Once R^Ω is obtained, the following task is to estimate the parameters θ^,P^. Note that Rwp-GLS does not assume any knowledge of K, which must be determined in order to acquire the final estimation step. Recall that rankTu^=K; hence, K determination can be formulated as rankTu^ estimation. In Rwp-GLS, based on the eigenvalue decomposition of Tu^, K is set as the number of eigenvalues that are greater than a predefined threshold. Considering the assumptions in Section 5.1, we acquire the parameters θ^,P^ from Tu^ uniquely by applying the Vandermonde decomposition lemma as Tu^=AΩθ^diagP^AΩHθ^ (readers are referred to [25], for more information).


### 5.4. Convergence Analysis

As shown, Rwp-GLS method converges to a local minimum, here we analyze its convergence from other side. It is well known that WNN minimization problem is nonconvex in general, and no efficient methods in theory can guarantee the global minimum. According to [41], WNN converges weakly to the solution if the weights are nondescendingly ordered, through a weighted soft-thresholding operator. By observing that σkTu,Ω1M+N−1 are nonnegative monotonically decreasing, then the ascending order of the resultant weights wk1NM will be kept throughout the reweighting process. Hence, we can conclude that Rwp-GLS method converges weakly to the solution.


### 5.5. Computational Complexity

The overall computational complexity of Rwp-GLS has two dominated parts, the first one is related to the computation of the reweighted CF, which is an iterative process, and the second is the cost ofθ^,P^ estimation. To make the following discussion more sense, let’s set M¯=M+N−1, and N¯=MN. According to the iterative part, the overall computational complexity equals the cost per iteration times the number of iterations. In each iteration, Rwp-GLS requires computing R˜Ω1/2 and Wl, where the numbers of required floating-point operations (flops) are OM¯2L+M¯3 and OM¯3, respectively. Then, solving the SDP must be performed, where the interior-point based SDPT3 solver requires OM¯4N¯2.5+N¯4.5 flops [40], in the case of N¯+M¯2 variable size, and N¯+M¯ as the dimension of the PSD matrix. Finally, applying the Vandermonde decomposition to estimate θ^,P^ consumes up to ON¯3.


## 5.1. Gridless CF

The key idea of gridless CF is exploiting the structure of the data covariance matrix and performing the estimation under ungrided parameterized domain with specific statistical considerations (i.e., uncorrelated sources and reasonable number of snapshots). Based on the assumptions in Section4, the uncorrelated data yΩ have the data covariance matrix RΩ∈ℂM+N−1×M+N−1 given by(4)RΩ=EyΩtyΩHt=AΩθdiagPAΩHθ+diagσΩ.Note thatRΩ−1 exists in the presence of noise for σΩ>0. RΩ is a highly nonlinear function of the DoAs parameters. To linearize the problem, a reparameterization of RΩ using a PSD Toeplitz matrix Tu is performed. Hence, according to Vandermonde decomposition lemma (see, e.g., [3]), RΩ can be linearly reparametrized as(5)RΩu,σΩ=Tu,Ω+diagσΩ,Tu≥0,σΩ>0,where Tu,Ω=ΓΩTuΓΩT=AΩθdiagPAΩHθ for Tu∈ℂNM×NM is a PSD Toeplitz matrix formed by the vector u∈ℂ1×NM. ΓΩ∈0,1M+N−1×NM denotes a selection matrix, such that the ith row of ΓΩ contains all zeros but a single 1 at the ri position. The importance of the parameterization in equation (5) is that Tu is linear in the vector u, which comprises the interesting parameters.Recall, the sample covariance matrix is set asR˜Ω=1/LYΩYΩH. In this paper, we are mainly interested in the case, where L≥M+N−1 (i.e., R˜Ω−1 exists with probability one) and K≤NM−1 (i.e., Tu is rank-deficient, then θ,P can be uniquely determined [25]). As well, we consider the same CF criterion as in the SPICE method, which is given by(6)fu,σΩ=RΩ−1/2R˜Ω−RΩR˜Ω−1/2F2=trRΩ−1R˜Ω+trR˜Ω−1RΩ−2NM.By using the linear parameterization in equation (5), the CF can be accomplished by minimizing f leading to the following equivalences:(7)minu,σ≥0trR˜Ω1/2Tu,Ω+diagσΩ−1R˜Ω1/2+trR˜Ω−1Tu,Ω+trR˜Ω−1diagσΩsubject toTu≥0,(8)minX,u,σ≥0trX+trR˜Ω−1Tu,Ω+trR˜Ω−1diagσΩsubject toTu≥0,X≥R˜Ω1/2Tu,Ω+diagσΩ−1R˜Ω1/2,(9)minX,u,σ≥0trX+trR˜Ω−1Tu,Ω+trR˜Ω−1diagσΩsubject toXR˜Ω−1/2R˜Ω−1/2Tu,Ω+diagσΩ≥0.Hence, iff minimization is casted as SDP, then it is convex. Problem (9) was derived under the development of the GLS method. In particular, the GLS method transforms the direction estimation problem into the estimation of a PSD Toeplitz matrix in which the parameters of the directions are encoded. Therefore, once the SDP in problem (9) is solved, the parameters of interest can be obtained from the estimated covariance matrix by either using the conventional subspace methods or applying the Vandermonde decomposition lemma.In spite of the strong theoretical guarantees of the GLS method, it suffers from considerable limitations according to Table1. In fact, we believe that its performance can be further improved by employing the nonconvex penalty, which to the best of our knowledge, has not been studied in the literature. Note that the nonconvex Schatten-p quasi-norm regularization is able to achieve a more accurate and sparser solution than the trace norm. In this paper, we reformulate the minimization problem of f in a new objective function by introducing Schatten-p quasi-norm. Then, the resultant optimization problem is modeled as a SDP formula.


## 5.2. Schatten-p Norm

Let⋅P be the Schatten-p norm of the matrix M∈ℂn×m, which is defined as the lp norm of its singular values σkM, i.e.,(10)Mp=∑k=1minn,mσkpM1/p=trMTMp/21/p,0<p<∞.Recently, Schatten-p norm has received significant attention from researchers in various domains [27, 33], due to its outstanding properties that make it a good choice in many optimization problems (e.g., rank function relaxation). Specifically, Mp equals to the nuclear norm when p=1. As well, if 0<p<1, then the quasi-norm (for 0<p<1, Mpis a quasi-norm because it does not satisfy the triangle inequality)Mp tends to rankM. Moreover, Zhang et al. [37] verified theoretically that Mp with small p requires significantly fewer measurements, with superior performance over the conventional nuclear norm in the majority of the problems.


## 5.3. The Proposed Method

Motivated by the properties of Schatten-p norm and GLS method, we develop a new method by replacing thetrR˜Ω−1Tu,Ω term in problem (8) by its Schatten-p quasi-norm, leading to the following optimization problem:(11)minX,u,σ≥0trX+λZupp+trR˜Ω−1diagσΩ,0<p<1subject toTu≥0,X≥R˜Ω1/2Tu,Ω+diagσΩ1R˜Ω1/2,where λ>0 is a regularization parameter and Zu=R˜Ω−1Tu,Ω is a PSD matrix clearly:(12)Zupp=∑k=1M+N−1σkpZu,0<p<1.According to equation (10), Zupp is given aswhere σkZu≥01M+N−1 are descendingly sorted singular values of Zu.Since0<p<1, then for any p,p1,p2>0 satisfying(13)1p=1p1+1p2,the bi-Schatten-p norm surrogate can be applied successfully [38]. In other words, the Schatten-p quasi-norm of Zu matrix is equivalent to minimizing the product of the Schatten-p1 norm (or quasi-norm) and Schatten-p2 norm (or quasi-norm) of its two factor matrices. So, as an alternative to the nonconvex quasi-norm in (12), the following efficient factorization is used:(14)1pZupp=minR˜Ω−1,Tu,ΩZ=R˜Ω−1Tu,Ω1p1R˜Ω−1p1p1+1p2Tu,Ωp2p2.Notice thatR˜Ω−1 is fixed, so whatever the value of p1 is on minimizing Problem (14), one only needs to minimize 1/p2Tu,Ωp2p2. Substituting the spectral regularization in Problem (11) with Problem (14) and setting λ=1/p, we infer the following problem:(15)minX,u,σ≥0trX+1p2Tu,Ωp2p2+trR˜Ω−1diagσΩsubject toTu≥0,X≥R˜Ω1/2Tu,Ω+diagσΩ−1R˜Ω1/2.With respect to equation (13), if p2≥1, Problem (15) is strongly convex; thus, it can be easily solved by off-the-shelf algorithms. When p2<1, then Problem (15) becomes nonconvex and must be relaxed. However, we empirically found that the later problem has superior properties in terms of estimation performance and solution sparsity. Therefore, in the following, we consider only the nonconvex model and work on the relaxation issue of Problem (15).Let’s setp2=p for simplicity. Since the proposed objective function in Problem (15) is nonconvex, then it must be linearized. A popular locally convergent approach to the minimization of such a problem is the majorization maximization (MM) method [39]. Accordingly, at the lth iteration, the Schatten-p quasi-norm term will be replaced by its tangent plane, which is given as(16)σkTu,Ωl+εp+pσkTu,Ωl+ε1−pσkTu,Ω−σkTu,Ωl,where ε >0 is a smoothing parameter. Notice that Tu,Ωl is constant. Consequently, by applying equation (16) into Problem (15), the resultant optimization formula will be rewritten as(17)minX,u,σ≥0trX+1pWlTu,Ω⋆+trR˜Ω−1diagσΩsubject toTu≥0,X≥R˜Ω1/2Tu,Ω+diagσΩ−1R˜Ω1/2.In Problem (17), the Schatten-p quasi-norm term is reduced to the weighted nuclear norm (WNN), where WlTu,Ω⋆=trWlTu,Ω=∑k=1M+N−1wklσkTu,Ω and Wl=D⋅diagw1l,w2l,…,wN+M−1l⋅DH is the reweighting matrix in the lth iteration, whose elements are wkl=p/σkTu,Ωl+ε1−p. Therefore, we name our method as reweighted p-GLS (Rwp-GLS). The matrix Tu,Ωl is updated based on the latest solution of the problem, and D is obtained by using the eigen decomposition of Tu,Ωl so that DDH=I.SinceTu,Ωpp is a concave function when 0<p<1, then at each iteration its value decreases significantly, faster than the decrease of its tangent plane. It follows that, by iteratively solving Problem (17), the objective function monotonically decreases and converges to a local minimum. Following the same steps in deriving Problem (9), at each iteration of Rwp-GLS, we need to solve the following SDP:(18)minX,u,σ≥0trX+1pWlTu,Ω⋆+trR˜Ω−1diagσΩsubject toconstraint of9.By solving the SDP in Problem (18) by off-the-shelf SDP solvers, e.g., SDPT3 [40], R^Ω=Tu^+diagσ^Ω is obtained given the solution u^,σ^Ω. Once R^Ω is obtained, the following task is to estimate the parameters θ^,P^. Note that Rwp-GLS does not assume any knowledge of K, which must be determined in order to acquire the final estimation step. Recall that rankTu^=K; hence, K determination can be formulated as rankTu^ estimation. In Rwp-GLS, based on the eigenvalue decomposition of Tu^, K is set as the number of eigenvalues that are greater than a predefined threshold. Considering the assumptions in Section 5.1, we acquire the parameters θ^,P^ from Tu^ uniquely by applying the Vandermonde decomposition lemma as Tu^=AΩθ^diagP^AΩHθ^ (readers are referred to [25], for more information).


## 5.4. Convergence Analysis

As shown, Rwp-GLS method converges to a local minimum, here we analyze its convergence from other side. It is well known that WNN minimization problem is nonconvex in general, and no efficient methods in theory can guarantee the global minimum. According to [41], WNN converges weakly to the solution if the weights are nondescendingly ordered, through a weighted soft-thresholding operator. By observing that σkTu,Ω1M+N−1 are nonnegative monotonically decreasing, then the ascending order of the resultant weights wk1NM will be kept throughout the reweighting process. Hence, we can conclude that Rwp-GLS method converges weakly to the solution.


## 5.5. Computational Complexity

The overall computational complexity of Rwp-GLS has two dominated parts, the first one is related to the computation of the reweighted CF, which is an iterative process, and the second is the cost ofθ^,P^ estimation. To make the following discussion more sense, let’s set M¯=M+N−1, and N¯=MN. According to the iterative part, the overall computational complexity equals the cost per iteration times the number of iterations. In each iteration, Rwp-GLS requires computing R˜Ω1/2 and Wl, where the numbers of required floating-point operations (flops) are OM¯2L+M¯3 and OM¯3, respectively. Then, solving the SDP must be performed, where the interior-point based SDPT3 solver requires OM¯4N¯2.5+N¯4.5 flops [40], in the case of N¯+M¯2 variable size, and N¯+M¯ as the dimension of the PSD matrix. Finally, applying the Vandermonde decomposition to estimate θ^,P^ consumes up to ON¯3.


## 6. Simulations and Results

In this section, the performance of Rwp-GLS method is tested and compared with state-of-art methods via numerical simulations. The implementation of Rwp-GLS is done using MATLAB v.14. a by SDPT3 solver of CVX toolbox [42], on a PC with 2 GHz CPU.


### 6.1. Simulation Setups

All the simulations are carried out for two equal-power uncorrelated signals (i.e.K=2), impinge onto 7-elements CPA with M=3, and N=5 with half-wavelength element spacing. Consider both signals and noise are i.i.d. Gaussian. The root mean squared error (RMSE) of the DoA estimation is computed as 1/TK∑i=1T∑k=1Kθk−θ^k,i2, where θkk=1K is the set of true directions, θ^kk=1K indicates the set of estimated directions and T is the number of Monte Carlo runs. The SNR is defined for equal-power sources P1=P2=…=PK=P as [28],(19)SNR=10log10PM+N−1∑i=1M+N−1σΩi−1,where σΩii=1M+N−1 are spatially uncorrelated white Gaussian random. Furthermore, to measure the variation of the noise on the antenna array, we use the worst noise power ratio (WNPR) as WNPR=σmax/σmin [28].At first, we discuss the choice ofε and p parameters and study the influence of applying different values of both on the performance of the proposed method. Since Problem (15) puts a penalty on 1/p∑k=1σkTu+εp, we plot the function hpσ=1/pσ+εp together with l0, and l1 norms in Figure 1. As shown in Figure 1(a), hpσ converges faster toward l0 whenever p⟶0 and yields to a sparser optimizer. While, by fixing p as small as 0.1 and using different values of ε, then hpσ will get close to l0 norm as ε⟶∞, and approach l1 norm when ε⟶0 as depicted in Figure 1(b).Figure 1
Property ofhpσ when (a)ε=1 and different values of p, (b)p=0.1 and different values of ε.
(a)
(b)However, after linearization, we expect that the resultant reweighting function in Problem (18) will play the same role as hpσ when ε⟶0 and p⟶0. Hence, throughout the iterative process, ε is initiated by 10−2 and reduced with a factor η in each iteration. Regarding the iterative process, the iteration is stopped if the relative change of u^ at two consecutive iterations is less than 10−4, or the maximum number of iterations which sets to 10 is reached. Finally, the results in each experiment are evaluated after 100 Monte Carlo runs.


### 6.2. DoA Estimation Performance Metrics

Consider two uncorrelated sources with spatial separation limitΔθ is defined as [23],(20)Δθ=minθk,θlθk≠θlsinθk−sinθl.Taking for exampleθ=17.46°,27.9° (i.e. Δθ=2.52/NM), SNR=0dB dB with σΩ=1,4,7.3,10,8,9,5 (i.e. WNPR=10) and L=200 snapshots.Initially, the success probability (SP) performance as a function of the number of iterations is examined. In Figure2(a), we plot SP for two different values of p∈0.1,0.5, while each is repeated for three values of factor η∈2,5,10. As shown, SP of Rwp-GLS is the best all over the number of iterations when p=0.1 regardless of η value. In contrast to the case of p=0.5, where SP slightly varies with η, and further it gets worse when the number of iteration becomes larger. However, using larger η is somewhat better. This indicates again the outperformance of Rwp-GLS with remarkable convergence characteristics for small values of p where the impact of η becomes negligible.Figure 2
Success probability performance versus number of iterations and various: (a)p and η, (b) reweighting function for p=0.1 and η=10.
(a)
(b)Figure2(b) points to the influence of the reweighting function choice that is applied in the first iteration. In the case of p=0.1, it is observed that SP is better when a nontrivial weighting based on the sample covariance matrix of the measurements is used rather than the weightless first iteration.Further, we test the performance of Rwp-GLS method for different values of0<p<1 with respect to different settings of WNPR and SNR. Figure 3(a) shows that the RMSE of Rwp-GLS perform the best, over the range 1≤WNPR≤200 for 0.1≤p≤0.2, where SNR = 0 dB. In Figure 3(b), we fix WNPR=10 and repeat the test for −20≤SNRdB≤20. The results show that the best performance is realized for 0.1≤p≤0.4.Figure 3
RMSE of Rwp-GLS with respect to p and (a) WNPR for SNR = 0 dB and (b) SNR for WNPR = 10.
(a)
(b)Regarding the results of the previous experiments, in the following simulations, we report the results considering the following settings:p=0.1, ε=10−2, and η=10, besides to adopt the weighted first iteration as the one tested in Figure 2(b).


### 6.3. Comparisons with Prior Arts

In this section, the DoA estimation performance of Rwp-GLS is evaluated and compared with the conventional RootMUSIC method, which is applied directly to the aforementioned CPA. In addition, we consider two remarkable gridless CF methods: SPA (we compare with SPA because it uses the same gridless CF procedure as GLS and agrees with the comparative methods in the model order estimation technique, while GLS adopts different technique) and ICMRA with Schatten-p penalty (ICMRAlp) to perform this comparison. These methods are derived under similar assumptions on sources and noise and applied with their best settings. We fix WNPR=10 and L=200 all over the following simulations unless otherwise is mentioned. Moreover, to deliver a benchmark for the evaluation, the stochastic Cramer–Rao lower bound (CRLB) for nonuniform noise given in [30] is included.Suppose two signals impinge onto the CPA from directionsθ=5.74°,15.55°. By varying the SNR from −15 dB to 25 dB, Figure 4 shows that the performance of Rwp-GLS is the best, and it improves persistently with increasing SNR and follows the CRLB when SNR ≥ −5 dB. Thus, Rwp-GLS can separate the two incident signals precisely when SNR is as low as −5 dB, in contrast to RootMUSIC, which suffers fixed error at high SNR due to the limited number of array elements.Figure 4
RMSEs of RootMUSIC, SPA, ICMRAlp, Rwp-GLS, and CRLB with SNR.Figure5 presents the CPU times of the experiment associated with Figure 4. Rwp-GLS is an iterative method with significant convergence characteristics due to well initialization and its computational complexity in one iteration has a similar complexity as SPA; then, Rwp-GLS costs a similar CPU time to SPA and differs in weighting computing. Note that the time usages of both Rwp-GLS and SPA become slightly greater when increasing the SNR because of the σΩ estimation. It is worth mentioning, ICMRAlp is an iterative method with convergence characteristics as Rwp-GLS, but without the ability to estimate the noise variances, which explains the slight decrease in CPU time-consuming.Figure 5
CPU times comparisons of RootMUSIC, SPA, ICMRAlp, and Rwp-GLS.In the following simulation, we evaluate the estimation performance of the Rwp-GLS versus the number of snapshots varying from 10 to 250. Consider the settings ofθ=5.74°,15.55°, SNR = 0 dB, and σΩ=1,8,2,10,3,6,5. Figure 6 shows that the estimation performance of Rwp-GLS outperforms the other methods and it follows well the CRLB as L≥50. In contrast to ICMRAlp, whose performance degrades significantly when L≤150.Figure 6
RMSEs of RootMUSIC, SPA, ICMRAlp, Rwp-GLS, and CRLB with number of snapshots L.Next, we examine the performance of DoA estimation versus angle separation. Assume two signals impinge onto the CPA from directions1°,1°+Δθ. The RMSEs of the comparative methods are shown in Figure 7, with Δθ varying from 0.5° to 5°, SNR=15 dB, and σΩ=1,8,2,10,3,6,5. It can be observed that Rwp-GLS surpasses the other methods in localizing spatially adjacent signals. In particular, Rwp-GLS is able to coincide with the CRLB when the sources are only 1° apart.Figure 7
RMSEs of RootMUSIC, SPA, ICMRAlp, Rwp-GLS, and CRLB with sources separation.In Figure8, the RMSEs of DoA estimation using the aforementioned methods are plotted against the variation of correlation factor 0.1≤ρ≤0.9. Furthermore, we carry out this simulation for two closely spaced sources, for example, let θ=5.74°,11.54° (i.e., Δθ=1.5/NM), where, SNR=0dB and σΩ=3.5,5,1,8,10,4.5,2.2 In Figure 8, we observe that, as the correlation factor increases the performances of the compared methods degrade, in particular, ICMRAlp is damaged significantly. However, the Rwp-GLS method has the best performance and exhibits immunity to the incorporated correlations among sources.Figure 8
RMSEs of RootMUSIC, SPA, ICMRAlp, and Rwp-GLS with different correlation factors.To further validate the robustness of the Rwp-GLS against the nonuniform noise, we present the RMSEs versus the WNPR in Figure9, where WNPR is varied from 20 to 200, under the settings of SNR = 0 dB and θ=5.74°,11.54° for two uncorrelated signals. It can be seen that, as WNPR increments, Rwp-GLS has the lowest RMSE over the studied WNPR range. This indicates again that the proposed method is robust against the nonuniform noise.Figure 9
RMSEs of RootMUSIC, SPA, ICMRAlp, and Rwp-GLS with WNPR.


### 6.4. Power Estimation

Finally, we attempt to estimate the DoAs and powers of 7 equal-power signals impinged onto the CPA from directionsθ=2°,8°,19°,29°,41°,55°,77° simultaneously. SNR is set as 10 dB and σΩ=1,4,3,8,5,7,10(i.e., WNRP = 10). We carry out 500 independent trials and show the results in Figure 10. Black circles denote the positions of the true powers and corresponding DoAs, while the colored dots denote the estimated ones. As shown, the Rwp-GLS method can correctly locate both DoAs and powers for all sources.Figure 10
DoAs and power estimation of Rwp-GLS forθ=2°,8°,19°,29°,41°,55°,77°.


## 6.1. Simulation Setups

All the simulations are carried out for two equal-power uncorrelated signals (i.e.K=2), impinge onto 7-elements CPA with M=3, and N=5 with half-wavelength element spacing. Consider both signals and noise are i.i.d. Gaussian. The root mean squared error (RMSE) of the DoA estimation is computed as 1/TK∑i=1T∑k=1Kθk−θ^k,i2, where θkk=1K is the set of true directions, θ^kk=1K indicates the set of estimated directions and T is the number of Monte Carlo runs. The SNR is defined for equal-power sources P1=P2=…=PK=P as [28],(19)SNR=10log10PM+N−1∑i=1M+N−1σΩi−1,where σΩii=1M+N−1 are spatially uncorrelated white Gaussian random. Furthermore, to measure the variation of the noise on the antenna array, we use the worst noise power ratio (WNPR) as WNPR=σmax/σmin [28].At first, we discuss the choice ofε and p parameters and study the influence of applying different values of both on the performance of the proposed method. Since Problem (15) puts a penalty on 1/p∑k=1σkTu+εp, we plot the function hpσ=1/pσ+εp together with l0, and l1 norms in Figure 1. As shown in Figure 1(a), hpσ converges faster toward l0 whenever p⟶0 and yields to a sparser optimizer. While, by fixing p as small as 0.1 and using different values of ε, then hpσ will get close to l0 norm as ε⟶∞, and approach l1 norm when ε⟶0 as depicted in Figure 1(b).Figure 1
Property ofhpσ when (a)ε=1 and different values of p, (b)p=0.1 and different values of ε.
(a)
(b)However, after linearization, we expect that the resultant reweighting function in Problem (18) will play the same role as hpσ when ε⟶0 and p⟶0. Hence, throughout the iterative process, ε is initiated by 10−2 and reduced with a factor η in each iteration. Regarding the iterative process, the iteration is stopped if the relative change of u^ at two consecutive iterations is less than 10−4, or the maximum number of iterations which sets to 10 is reached. Finally, the results in each experiment are evaluated after 100 Monte Carlo runs.


## 6.2. DoA Estimation Performance Metrics

Consider two uncorrelated sources with spatial separation limitΔθ is defined as [23],(20)Δθ=minθk,θlθk≠θlsinθk−sinθl.Taking for exampleθ=17.46°,27.9° (i.e. Δθ=2.52/NM), SNR=0dB dB with σΩ=1,4,7.3,10,8,9,5 (i.e. WNPR=10) and L=200 snapshots.Initially, the success probability (SP) performance as a function of the number of iterations is examined. In Figure2(a), we plot SP for two different values of p∈0.1,0.5, while each is repeated for three values of factor η∈2,5,10. As shown, SP of Rwp-GLS is the best all over the number of iterations when p=0.1 regardless of η value. In contrast to the case of p=0.5, where SP slightly varies with η, and further it gets worse when the number of iteration becomes larger. However, using larger η is somewhat better. This indicates again the outperformance of Rwp-GLS with remarkable convergence characteristics for small values of p where the impact of η becomes negligible.Figure 2
Success probability performance versus number of iterations and various: (a)p and η, (b) reweighting function for p=0.1 and η=10.
(a)
(b)Figure2(b) points to the influence of the reweighting function choice that is applied in the first iteration. In the case of p=0.1, it is observed that SP is better when a nontrivial weighting based on the sample covariance matrix of the measurements is used rather than the weightless first iteration.Further, we test the performance of Rwp-GLS method for different values of0<p<1 with respect to different settings of WNPR and SNR. Figure 3(a) shows that the RMSE of Rwp-GLS perform the best, over the range 1≤WNPR≤200 for 0.1≤p≤0.2, where SNR = 0 dB. In Figure 3(b), we fix WNPR=10 and repeat the test for −20≤SNRdB≤20. The results show that the best performance is realized for 0.1≤p≤0.4.Figure 3
RMSE of Rwp-GLS with respect to p and (a) WNPR for SNR = 0 dB and (b) SNR for WNPR = 10.
(a)
(b)Regarding the results of the previous experiments, in the following simulations, we report the results considering the following settings:p=0.1, ε=10−2, and η=10, besides to adopt the weighted first iteration as the one tested in Figure 2(b).


## 6.3. Comparisons with Prior Arts

In this section, the DoA estimation performance of Rwp-GLS is evaluated and compared with the conventional RootMUSIC method, which is applied directly to the aforementioned CPA. In addition, we consider two remarkable gridless CF methods: SPA (we compare with SPA because it uses the same gridless CF procedure as GLS and agrees with the comparative methods in the model order estimation technique, while GLS adopts different technique) and ICMRA with Schatten-p penalty (ICMRAlp) to perform this comparison. These methods are derived under similar assumptions on sources and noise and applied with their best settings. We fix WNPR=10 and L=200 all over the following simulations unless otherwise is mentioned. Moreover, to deliver a benchmark for the evaluation, the stochastic Cramer–Rao lower bound (CRLB) for nonuniform noise given in [30] is included.Suppose two signals impinge onto the CPA from directionsθ=5.74°,15.55°. By varying the SNR from −15 dB to 25 dB, Figure 4 shows that the performance of Rwp-GLS is the best, and it improves persistently with increasing SNR and follows the CRLB when SNR ≥ −5 dB. Thus, Rwp-GLS can separate the two incident signals precisely when SNR is as low as −5 dB, in contrast to RootMUSIC, which suffers fixed error at high SNR due to the limited number of array elements.Figure 4
RMSEs of RootMUSIC, SPA, ICMRAlp, Rwp-GLS, and CRLB with SNR.Figure5 presents the CPU times of the experiment associated with Figure 4. Rwp-GLS is an iterative method with significant convergence characteristics due to well initialization and its computational complexity in one iteration has a similar complexity as SPA; then, Rwp-GLS costs a similar CPU time to SPA and differs in weighting computing. Note that the time usages of both Rwp-GLS and SPA become slightly greater when increasing the SNR because of the σΩ estimation. It is worth mentioning, ICMRAlp is an iterative method with convergence characteristics as Rwp-GLS, but without the ability to estimate the noise variances, which explains the slight decrease in CPU time-consuming.Figure 5
CPU times comparisons of RootMUSIC, SPA, ICMRAlp, and Rwp-GLS.In the following simulation, we evaluate the estimation performance of the Rwp-GLS versus the number of snapshots varying from 10 to 250. Consider the settings ofθ=5.74°,15.55°, SNR = 0 dB, and σΩ=1,8,2,10,3,6,5. Figure 6 shows that the estimation performance of Rwp-GLS outperforms the other methods and it follows well the CRLB as L≥50. In contrast to ICMRAlp, whose performance degrades significantly when L≤150.Figure 6
RMSEs of RootMUSIC, SPA, ICMRAlp, Rwp-GLS, and CRLB with number of snapshots L.Next, we examine the performance of DoA estimation versus angle separation. Assume two signals impinge onto the CPA from directions1°,1°+Δθ. The RMSEs of the comparative methods are shown in Figure 7, with Δθ varying from 0.5° to 5°, SNR=15 dB, and σΩ=1,8,2,10,3,6,5. It can be observed that Rwp-GLS surpasses the other methods in localizing spatially adjacent signals. In particular, Rwp-GLS is able to coincide with the CRLB when the sources are only 1° apart.Figure 7
RMSEs of RootMUSIC, SPA, ICMRAlp, Rwp-GLS, and CRLB with sources separation.In Figure8, the RMSEs of DoA estimation using the aforementioned methods are plotted against the variation of correlation factor 0.1≤ρ≤0.9. Furthermore, we carry out this simulation for two closely spaced sources, for example, let θ=5.74°,11.54° (i.e., Δθ=1.5/NM), where, SNR=0dB and σΩ=3.5,5,1,8,10,4.5,2.2 In Figure 8, we observe that, as the correlation factor increases the performances of the compared methods degrade, in particular, ICMRAlp is damaged significantly. However, the Rwp-GLS method has the best performance and exhibits immunity to the incorporated correlations among sources.Figure 8
RMSEs of RootMUSIC, SPA, ICMRAlp, and Rwp-GLS with different correlation factors.To further validate the robustness of the Rwp-GLS against the nonuniform noise, we present the RMSEs versus the WNPR in Figure9, where WNPR is varied from 20 to 200, under the settings of SNR = 0 dB and θ=5.74°,11.54° for two uncorrelated signals. It can be seen that, as WNPR increments, Rwp-GLS has the lowest RMSE over the studied WNPR range. This indicates again that the proposed method is robust against the nonuniform noise.Figure 9
RMSEs of RootMUSIC, SPA, ICMRAlp, and Rwp-GLS with WNPR.


## 6.4. Power Estimation

Finally, we attempt to estimate the DoAs and powers of 7 equal-power signals impinged onto the CPA from directionsθ=2°,8°,19°,29°,41°,55°,77° simultaneously. SNR is set as 10 dB and σΩ=1,4,3,8,5,7,10(i.e., WNRP = 10). We carry out 500 independent trials and show the results in Figure 10. Black circles denote the positions of the true powers and corresponding DoAs, while the colored dots denote the estimated ones. As shown, the Rwp-GLS method can correctly locate both DoAs and powers for all sources.Figure 10
DoAs and power estimation of Rwp-GLS forθ=2°,8°,19°,29°,41°,55°,77°.


## 7. Conclusions

In this paper, we have proposed a gridless sparse method for DoA estimation named Rwp-GLS. The new method reformulates the DoA estimation, by introducing the unified surrogate of the nonconvex Schatten-p quasi-norm to penalize the parametrized covariance matrix in the gridless CF criterion. Then, an iterative reweighted minimization method was derived with a SDP formulation. The effectiveness of Rwp-GLS is verified through numerical simulations under the consideration of CPA structure and unknown nonuniform noise background. The results show the superiority of Rwp-GLS in terms of resolution, applicability in a limited number of snapshots, and immunity against nonuniform noise and correlation of sources. In the future, we look forward to reasonable computational efficiency implementation and arbitrary array geometries incorporation.



---

*Source: 3012952-2020-04-27.xml*