# Augmented Lagrange Based on Modified Covariance Matching Criterion Method for DOA Estimation in Compressed Sensing

Authors: Weijian Si; Xinggen Qu; Lutao Liu
Journal: The Scientific World Journal
 (2014)

---

## Abstract

A novel direction of arrival (DOA) estimation method in compressed sensing (CS) is presented, in which DOA estimation is considered as the joint sparse recovery from multiple measurement vectors (MMV). The proposed method is obtained by minimizing the modified-based covariance matching criterion, which is acquired by adding penalties according to the regularization method. This minimization problem is shown to be a semidefinite program (SDP) and transformed into a constrained quadratic programming problem for reducing computational complexity which can be solved by the augmented Lagrange method. The proposed method can significantly improve the performance especially in the scenarios with low signal to noise ratio (SNR), small number of snapshots, and closely spaced correlated sources. In addition, the Cramér-Rao bound (CRB) of the proposed method is developed and the performance guarantee is given according to a version of the restricted isometry property (RIP). The effectiveness and satisfactory performance of the proposed method are illustrated by simulation results.


---

## Body

## 1. Introduction

Direction of arrival (DOA) estimation of multiple narrow-band sources has been an interesting research topic in array signal processing. Its applications span many fields including radar, communication systems, and medical imaging [1, 2]. Many effective algorithms are proposed for DOA estimation, which are mainly classified into three categories. The beamforming algorithms such as MVDR [3] and MEM [4] obtain a nonparametric spatial spectrum by optimizing the filter weights. The subspace algorithms such as MUSIC [5] and ESPRIT [6] and their derivatives [7, 8] exploit the orthogonality of signal subspace and noise subspace for DOA estimation. The subspace fitting algorithms including Maximum Likelihood (ML) [9] and Weighted Subspace Fitting (WSF) [10] solve a multidimensional nonlinear optimization problem to obtain a precise estimation, but a good initialization is required to ensure global convergence and computational complexity is very high. All these algorithms focus on two important issues: resolution (i.e., the ability to correctly resolve two closely spaced sources) and precision (i.e., the deviation from true DOA) which are considered to be the theoretical bases of evaluating certain algorithms. Many high resolution methods suffer from serious performance degradation in the scenarios with low SNR, small number of snapshots, and closely spaced corrected sources. More recently, many research applications involving compressed sensing (CS) [11], especially DOA estimation [12, 13], have become more and more popular in the signal processing. Moreover, the restricted isometry property (RIP) based on the perfect theoretical framework given with modern probability theory by Donoho [11] and Candés et al. [14] provides the performance guarantee in CS. Hence, in this paper, DOA estimation is posed as a joint sparse recovery problem where we recover jointly sparse sources from multiple measurement vectors (MMV) under CS framework.CS is an emerging area which can break through the limit of Nyquist sampling theorem. On one hand, CS can simultaneously capture and store compressed or sparse source at a rate much lower than that of Nyquist sampling. On the other hand, it can recover the original source using nonadaptive linear projection measurements onto a suitable measurement matrix which satisfies the RIP. In CS, the joint sparse recovery aims to find a common support shared by unknown sparse vectors from MMV, which is obtained by the measurement matrix. Support denotes the indices of the nonzero elements in the unknown sparse vectors. A sparse solution can be obtained as long as the support is determined.CS theory has been widely applied to DOA estimation according to the source sparsity, which results from the fact that there are much fewer source directions than all potential directions in the spatial domain. The DOA estimation methods in CS are attractive since they have much better estimation performance than conventional estimation methods. In [15], Gürbüz et al. firstly formulate the DOA estimation problem under CS framework in the time domain. Wang et al. [16] propose a new architecture to estimate DOA by exploiting compressive sampling in the spatial domain. Stoica et al. [17] make full use of covariance matching criterion and present a semiparametric/sparse estimation method and its derivative called LIKES [18] for the separable model. Malioutov et al. [19] present the l
1-SVD algorithm for DOA estimation which combines the SVD of the data matrix with a sparse recovery method based on l
1-norm minimization. A new class of subspace-base algorithms such as compressive MUSIC (CS-MUSIC) [20] and subspace-augmented MUSIC (SA-MUSIC) [21] is proposed in recent years.The RIP and various modified versions of it have been used as a foundation of performance guarantee [21–24] for the joint sparse recovery. The performance guarantee of MUSIC based on joint sparse recovery is given for identifying the unique support in a favorable case [25]. However, in the unfavorable case where the number of measurement vectors is smaller than the sparsity or the covariance matrix tends to lose rank due to correlated sources, performance guarantee fails. Lee et al. [21] propose a new performance guarantee in terms of a version of the RIP under such unfavorable conditions. The performance guarantees of other methods such as greedy algorithms and convex relaxation have been developed to find the sparse solution in [23, 24]. However, the guarantees of these methods cannot be simply extended to the MMV case to obtain a better bound for the sparse solution.In this paper, we propose a novel augmented Lagrange based on modified covariance matching criterion method called AULMC for DOA estimation in CS. This method utilizes the minimization of the modified covariance matching criterion which is acquired by using regularization method to add penalties in order to obtain the stable sparse parameter estimation, especially in the low sparsity case. This minimization problem is shown to be a semidefinite program (SDP) and transformed into the constrained quadratic programming problem for the sake of reducing computational complexity. The augmented Lagrange function is formed to solve this problem by the use of augmented Lagrange method. AULMC has a number of advantages over the other methods. For example, it provides more precise estimation and higher resolution in the scenarios with low SNR, small number of snapshots, and closely spaced correlated sources, and it does not need priori knowledge about the number of sources or to choose the regularization parameter of thel
1-optimization framework which is very difficult to select in the DOA estimation. In addition, we give a detailed derivation process of the closed-form expression for the Cramér-Rao bound (CRB) of the new method and discuss an explicit condition that guarantees performance of the new method. This performance guarantee is given in terms of weaker version of the RIP which is referred to as weak-1 RIP. Simulation results illustrate the performance of the proposed method.It is worth noting that covariance matching criterion has been used for DOA estimation [26]. In [26], a sparse iterative covariance-base estimation method, abbreviated as SPICE, is proposed. Our approach is different from SPICE because it utilizes modified covariance matching criterion which can guarantee the stability of solution even if the source sparsity is rather low. In the future work, we will focus on the application of our approach to data-driven design [27–30]. Now we briefly summarize the contributions of this work as follows.(i)
A modified covariance matching criterion is proposed by adding penalties according to the regularization method.
(ii)
The original SDP problem is transformed into the constrained quadratic programming problem. The motivation to transform the problem is that it can reduce computational complexity.
(iii)
Augmented Lagrange based on modified covariance matching criterion method is devised to solve the resulting programming problem.
(iv)
The CRB and performance guarantee of the new method are given in detail.The rest of this paper is organized as follows. In Section2, we formulate the DOA estimation problem. Section 3 presents a modified covariance matching criterion. A novel augmented Lagrange based on modified covariance matching criterion method for DOA estimation is described in detail in Section 4; the performance of which is analyzed in Section 5. The performance of the proposed method is evaluated by simulations in Section 6. Conclusions are provided in Section 7.


## 2. Problem Formulation

ConsiderL narrow-band sources s
^
1
,
s
^
2
,
…
,
s
^
L impinging on a sensor array that consisted of P omnidirectional sensors from far field. At time instant t, the array received source can be given by(1)
x
(
t
)
=
∑
k
=
1
L
‍
a
(
θ
k
)
s
^
k
(
t
)
+
n
(
t
)
,
where n
(
t
)
∈
𝒞
P
×
1 denotes a noise term, θ
k
∈
Ω is the unknown direction, of the kth source where Ω denotes the entire spatial range and a
(
θ
k
) is the P
×
1 steering vector. Although the DOA estimation of the single snapshot, which is a typical single measurement vector (SMV) model, has its value, the number of snapshots is larger than one in the most practical applications. Correspondingly, the multiple snapshots model is a MMV model.Let{
θ
~
k
}
k
=
1
K denote a fine grid which covers Ω where there exist K
(
K
≫
L
) potential directions of the sources s
1
,
s
2
,
…
,
s
K so that the true directions {
θ
k
}
k
=
1
L are aligned or are close to the grids. This means that if θ
~
k
1
,
θ
~
k
2
,
…
,
θ
~
k
L are equal to θ
1
,
θ
2
,
…
,
θ
L, respectively, we have
(2)
s
k
=
{
s
^
l
,
k
=
k
l
(
l
=
1,2
,
…
,
L
)
,
0
,
elsewhere
.Hence, the multiple snapshots model can be written as the following sparse form:(3)
x
(
t
)
=
∑
k
=
1
K
‍
a
(
θ
~
k
)
s
k
(
t
)
+
n
(
t
)
=
A
s
(
t
)
+
n
(
t
)
,
t
=
t
1
,
t
2
…
t
M
,
where A
=
[
a
(
θ
~
1
)
a
(
θ
~
2
)
⋯
a
(
θ
~
K
)
] is the P
×
K manifold matrix corresponding to all the potential directions which is also referred to as an overcomplete dictionary in CS. s
(
t
)
=
[
s
1
(
t
)
s
2
(
t
)
⋯
s
K
(
t
)
]
T is a L-sparse vector since it has at most L nonzero elements in K elements, and L is defined as sparsity where the operator [
·
]
T denotes transport. {
s
(
t
i
)
}
i
=
1
M are jointly L-sparse if they share a common support. Hence, the matrix S
=
[
s
(
t
1
)
s
(
t
2
)
⋯
s
(
t
M
)
]
∈
𝒞
K
×
M has no more than L nonzero rows in order to be called row L-sparse. The MMV problem is that of identifying the row support of the unknown matrix S from the matrix Y
∈
𝒞
N
×
M that consisted of MMV which is given by
(4)
Y
=
[
y
(
t
1
)
y
(
t
2
)
⋯
y
(
t
M
)
]
=
𝚽
A
S
+
𝚽
N
with a common measurement matrix 𝚽 of the size N
×
P with N
<
P where N is the number of nonadaptive linear projection measurements, such as random Gaussian matrix or random partial Fourier matrix, and noise matrix N.


## 3. Modified Covariance Matching Criterion

In this section, the modified covariance matching criterion is developed in the CS scenario. A conventional covariance matrix of the compressed measurement source with the sizeN
×
N is given by
(5)
R
y
=
E
[
y
(
t
)
y
*
(
t
)
]
=
𝚽
A
R
s
A
*
𝚽
*
+
𝚽
R
n
𝚽
*
,
where R
s
=
E
[
s
(
t
)
s
*
(
t
)
] is a K
×
K covariance matrix of the sparse source whose off-diagonal elements denote the source correlation and diagonal elements denote the source powers. Since the powers are one-to-one corresponding to all the potential directions and our focus is on the source angle parameter estimation, R
s can be reduced to a diagonal matrix R
s
(
θ
~
)
=
diag
⁡
(
θ
~
1
θ
~
2
⋯
θ
~
K
). According to (5), the measurement matrix can change the covariance matrix of the noise to render the noise colored even if the noise is white. Therefore, this adverse factor must be considered before recovering jointly sparse sources. In addition, the operators (
·
)
* and E
[
·
] denote conjugate transpose and expectation, respectively.Since𝚽
R
n
𝚽
* is a positive definite Hermitian matrix, a prewhitened method is given by the Cholesky decomposition. Let B be the Cholesky factor that satisfies
(6)
(
𝚽
R
n
𝚽
*
)
-
1
=
B
*
B
,
where B
∈
𝒞
N
×
N is an upper triangular matrix of positive line. Hence, a prewhitened process is implemented by multiplying y
(
t
) by B in order to obtain a pure source whose covariance matrix is written as
(7)
R
~
y
=
B
𝚽
A
R
s
A
*
𝚽
*
B
*
+
I
N
=
C
R
s
C
*
+
I
N
,
where C
=
[
c
1
,
c
2
,
…
,
c
K
] and I
N is an identity matrix of the size N
×
N. It is important to note that the unknown covariance matrix of the noise is transformed into the known identity matrix after the prewhitened process. Therefore, the prewhitened method improves the robustness to the noise. Then, the covariance matrix of the compressed measurement source denoising is realized by
(8)
R
=
R
~
y
-
I
N
=
C
R
s
C
*
.The parameterθ
~ can be estimated by a class of the covariance matching estimation techniques (COMET) based on covariance matching criterion [31]. This parameter estimation method is well known to be a large-sample approximation of ML method and provide a more attractive solution than ML estimator.The principle of COMET is that of using the right data to minimize its data model by the weighted least-squares (WLS) method. However, the lower the source sparsity is, the more likely it is to be ill-posed for the covariance matrix estimation error meaning that the optimal solution obtained directly by minimizing the conventional covariance matching criterion is instable. Hence, we employ regularization method to add penalties in order to sufficiently exploit prior knowledge to reduce the solution space for determining the stable optimal solution. The modified covariance matching criterion is proposed as the following form:(9)
[
vec
⁡
(
R
^
)
-
vec
⁡
(
R
(
θ
~
)
)
]
*
Γ
*
P
-
1
Γ
[
vec
⁡
(
R
^
)
-
vec
⁡
(
R
(
θ
~
)
)
]
+
α
Tr
⁡
[
R
-
1
(
θ
~
)
-
R
^
-
1
]
+
β
Tr
⁡
[
R
-
1
(
θ
~
)
R
^
-
1
4
]
,
where parameter α
≥
0 controls the solution smoothness (guarantee precision), parameter β
≥
0 controls the solution scale (guarantee sparsity), the inverses of the matrix Γ, the sample covariance matrix R
^, and R
(
θ
~
) are assumed to exist, and the matrix P
-
1 is the inverse of the covariance matrix of the residuals, ɛ
~
=
Γ
vec
⁡
(
R
^
-
R
(
θ
~
)
). Since R
^ is equal to R
(
θ
~
) as M
→
∞, it follows from [32] that vec
⁡
(
R
^
-
R
(
θ
~
)
) satisfies the asymptotic normal distribution with mean zero and variance N
-
1
R
T
(
θ
~
)
⊗
R
^. In addition, the operators Tr
⁡
[
·
], vec
⁡
(
·
), and ⊗ denote matrix trace, column stacking operation, and Kronecker product, respectively. Then, the matrix P can be given by
(10)
P
=
E
[
ɛ
~
(
θ
~
)
ɛ
~
*
(
θ
~
)
]
=
Γ
E
[
vec
⁡
(
R
^
-
R
(
θ
~
)
)
ve
c
*
(
R
^
-
R
(
θ
~
)
)
]
Γ
*
=
Γ
N
-
1
R
T
(
θ
~
)
⊗
R
^
Γ
*
=
Γ
N
-
1
D
Γ
*
.We consider the normalized data model in (9) and choose the regularization parameters α
=
β
=
N based on perceptual criterion [33]. By substituting (10) into (9), we have
(11)
ve
c
*
(
R
^
-
R
(
θ
~
)
)
D
-
1
vec
⁡
(
R
^
-
R
(
θ
~
)
)
+
Tr
⁡
[
R
-
1
(
θ
~
)
-
R
^
-
1
]
+
Tr
⁡
[
R
-
1
(
θ
~
)
R
^
-
1
]
4
.By the properties ofvec
⁡, ⊗, and Tr
⁡, the data model can be further simplified to
(12)
ve
c
*
(
R
^
-
R
(
θ
~
)
)
D
-
1
vec
⁡
(
R
^
-
R
(
θ
~
)
)
+
Tr
⁡
[
R
-
1
(
θ
~
)
(
R
^
-
R
(
θ
~
)
)
R
^
-
1
]
+
Tr
⁡
[
R
-
1
(
θ
~
)
R
^
-
1
]
4
=
∥
R
^
-
1
/
2
(
R
^
-
R
(
θ
~
)
+
I
N
2
)
R
-
1
/
2
(
θ
~
)
∥
2
,
where ∥
·
∥ denotes the Frobenius norm for matrices or the Euclidean norm for vectors. The data model in (12) is considered to be the modified covariance matching criterion. It can be seen from (12) that a positive definite Hermitian matrix is added to the covariance matrix estimation error by applying penalties according to regularization method in order to guarantee the stability of solution.


## 4. DOA Estimation

In this section, we will utilize the minimization of the modified covariance matching criterion to estimate parameterθ
~. Let(13)
θ
~
s
=
arg
min
⁡
θ
~
∥
R
^
-
1
/
2
(
R
^
-
R
(
θ
~
)
+
I
2
)
R
-
1
/
2
(
θ
~
)
∥
2
be the optimal solution of  θ
~ in the structure model of (8). Our goal is to utilize the modified covariance matching criterion for an estimate that is asymptotically equal to θ
~
s. By the properties of the trace and the Hermitian matrix, a derivation process is shown as follows, where we omit the dependence on  θ
~ for notational convenience:
(14)
f
=
Tr
⁡
[
R
^
-
1
(
R
^
-
R
+
I
2
)
R
-
1
(
R
^
-
R
+
I
2
)
]
=
Tr
⁡
[
R
^
-
1
R
-
R
^
-
1
+
(
R
^
1
/
2
+
1
2
R
^
-
1
/
2
)
=
Tr
⁡
[
l
·
R
-
1
/
2
(
R
^
1
/
2
+
1
2
R
^
-
1
/
2
)
R
-
1
/
2
]
.Due to(15)
Tr
⁡
[
R
^
-
1
R
]
=
∑
k
=
1
K
‍
(
c
k
*
R
^
-
1
c
k
)
θ
~
k
we can deduce that the minimization of f is equal to the minimization of h:
(16)
h
=
Tr
⁡
[
(
R
^
1
/
2
+
1
2
R
^
-
1
/
2
)
R
-
1
/
2
(
R
^
1
/
2
+
1
2
R
^
-
1
/
2
)
R
-
1
/
2
]
+
∑
k
=
1
K
‍
(
c
k
*
R
^
-
1
c
k
)
θ
~
k
.Then, we will demonstrate that the minimization ofh in (16) with respect to {
θ
~
k
}
k
=
1
K is an SDP. To prove this fact, One assumes
(17)
R
^
1
/
2
+
1
2
R
^
-
1
/
2
=
[
r
1
r
2
⋯
r
N
]
.Theh in (16) can be rewritten as
(18)
h
=
∑
k
=
1
N
‍
r
k
*
R
-
1
r
k
+
∑
k
=
1
K
‍
(
c
k
*
R
^
-
1
c
k
)
θ
~
k
.By the Schur complement, let{
x
k
}
k
=
1
N be auxiliary variables satisfying
(19)
[
x
k
r
k
*
r
k
R
]
≥
0
.The equivalent form of this minimization problem is expressed as(20)
min
⁡
x
k
,
θ
~
k
∑
k
=
1
N
x
k
+
∑
k
=
1
K
(
c
k
*
R
^
-
1
c
k
)
θ
~
k
s
.
t
.
θ
~
k
≥
0
,
k
=
1,2
,
…
,
K
,
[
x
k
r
k
*
r
k
R
]
≥
0
,
k
=
1,2
,
…
,
N
.It is easy to see that (20) is an SDP [34]. Many software packages can solve an SDP, but solving (16) as an SDP is not a good choice because SDP solvers have generally rather high computational complexity for the values of N, M, and K in the DOA estimation. To solve it effectively, we transform it into the constrained quadratic programming problem for reducing computational complexity, as described next.Since a consistent estimation ofθ
~ is given by (15), we can reformulate the minimization of h in (16) as the following constrained minimization by the Schur inequality of the trace:
(21)
min
⁡
θ
i
≥
0
Tr
⁡
[
(
R
^
1
/
2
+
1
2
R
^
-
1
/
2
)
R
-
1
/
2
(
R
^
1
/
2
+
1
2
R
^
-
1
/
2
)
R
-
1
/
2
]
s
.
t
.
W
z
=
1
K
,
where W is a K
×
K diagonal matrix of [
C
*
R
^
-
1
C
]
-
1
/
2 and z
=
[
z
1
z
2
⋯
z
K
]
T is a K-dimensional vector with z
i
=
θ
i
-
1
/
2
,
i
=
1,2
,
…
,
K. By substituting (8) into (21), the objective function in (21) can be rewritten as
(22)
Tr
⁡
[
(
R
^
1
/
2
+
1
2
R
^
-
1
/
2
)
(
C
*
)
-
1
/
2
R
s
-
1
/
2
C
-
1
/
2
i
i
i
i
i
i
i
·
(
R
^
1
/
2
+
1
2
R
^
-
1
/
2
)
(
C
*
)
-
1
/
2
R
s
-
1
/
2
C
-
1
/
2
]
=
Tr
⁡
[
C
-
1
/
2
(
R
^
1
/
2
+
1
2
R
^
-
1
/
2
)
(
C
*
)
-
1
/
2
R
s
-
1
/
2
i
i
i
i
i
i
i
i
i
i
i
i
i
·
C
-
1
/
2
(
R
^
1
/
2
+
1
2
R
^
-
1
/
2
)
(
C
*
)
-
1
/
2
R
s
-
1
/
2
]
.Based on the following equation:(23)
Tr
⁡
[
M
T
(
d
)
N
T
(
d
)
]
=
d
*
(
M
T
⊙
N
)
d
,
where ⊙ denotes the Schur-Hadamard product, M and N are both K
×
K matrices, T
(
d
)
=
diag
⁡
(
d
1
d
2
⋯
d
K
), and d
=
[
d
1
d
2
⋯
d
K
]
T, the objective function (22) can be written as(24)
Tr
⁡
[
Q
R
s
-
1
/
2
Q
R
s
-
1
/
2
]
=
z
*
V
z
,
where Q
=
C
-
1
/
2
(
R
^
1
/
2
+
(
1
/
2
)
R
^
-
1
/
2
)
(
C
*
)
-
1
/
2 and V
=
Q
T
⊙
Q. Hence, the minimization problem in (21) is transformed into the following form:
(25)
min
⁡
z
f
(
z
)
=
z
*
V
z
s
.
t
.
Π
(
z
)
=
W
z
-
1
K
=
0
K
which is a typical constrained quadratic programming problem.It is well known that an important class of methods for solving the constrained quadratic programming problem is to form the auxiliary function. To solve (25), we form the following augmented Lagrange function with respect to (25):
(26)
p
σ
(
z
,
μ
¯
)
=
f
(
z
)
+
μ
¯
*
Π
(
z
)
+
1
2
σ
∥
Π
(
z
)
∥
2
,
where μ
¯ is the asymptotical solution of the Lagrange multiplier of (25) and σ is a penalty factor.By setting the gradient and the Hessian matrix of (26) with respect to z to zero, we have
(27)
∇
z
p
σ
(
z
,
μ
¯
)
=
∇
z
l
(
z
,
μ
¯
)
+
σ
q
(
z
)
Π
(
z
)
,
∇
z
z
2
p
σ
(
z
,
μ
¯
)
=
∇
z
z
2
l
(
z
,
μ
¯
)
+
σ
∑
i
=
1
K
‍
Π
i
(
z
)
∇
2
Π
(
z
)
+
σ
q
(
z
)
q
*
(
z
)
,
where l
(
z
,
μ
¯
)
=
f
(
z
)
+
μ
¯
*
Π
(
z
) and q
(
z
)
=
∇
Π
*
(
z
)
=
[
∇
Π
1
(
z
)
∇
Π
2
(
z
)
⋯
∇
Π
K
(
z
)
]. One assumes that
(28)
b
(
z
,
μ
¯
)
=
∇
z
z
2
l
(
z
,
μ
¯
)
+
σ
∑
i
=
1
K
‍
Π
i
(
z
)
∇
2
Π
(
z
)
.Applying the Newton method, we obtain(29)
∇
z
z
2
p
σ
k
(
z
k
,
μ
¯
k
)
(
z
k
+
1
-
z
k
)
=
-
∇
z
p
σ
k
(
z
k
,
μ
¯
k
)
.For notational convenience, we assume thatq
k
=
q
(
z
k
), ∇
f
k
=
∇
f
(
z
k
), Π
k
=
Π
(
z
k
), b
k
=
b
(
z
k
,
μ
¯
k
), d
k
=
z
k
+
1
-
z
k and we get the following equation by (29):
(30)
b
k
d
k
+
σ
k
q
k
q
k
*
d
k
+
∇
z
l
(
z
k
,
μ
¯
k
)
+
σ
k
q
k
Π
k
=
0
.Assuming that the inverse ofb
k exists, by left-multiplying (30) by (
1
/
σ
k
)
q
k
*
b
k
-
1, we have
(31)
(
I
K
σ
k
+
q
k
*
b
k
-
1
q
k
)
q
k
*
d
k
=
-
q
k
*
b
k
-
1
q
k
Π
k
-
1
σ
k
q
k
*
b
k
-
1
∇
z
l
(
z
k
,
μ
¯
k
)
.It follows from (31) that(32)
q
k
*
b
k
=
-
Π
k
+
α
k
σ
k
,
where α
k satisfies
(33)
(
I
K
σ
k
+
q
k
*
b
k
-
1
q
k
)
α
k
=
-
q
k
*
b
k
-
1
∇
z
l
(
z
k
,
μ
¯
k
)
+
Π
k
.By substituting (32) into (30), we have
(34)
d
k
=
-
b
k
-
1
(
q
k
α
k
+
∇
z
l
(
z
k
,
μ
¯
k
)
)
.
Both the multiplier factor and the penalty factor are determined with difficulty for utilizing the augmented Lagrange function to solve the constrained quadratic programming problem. Hence, an updated sequence for the multiplier factor is given in terms of Proposition 1.Proposition 1.
One assumes thatμ
(
z
) is the optimal solution of the problem min
⁡
μ
∈
𝒞
K
∥
q
(
z
)
μ
+
∇
f
(
z
)
∥
2. Then, the following equation holds:
(35)
μ
(
z
k
)
=
μ
¯
k
+
α
k
+
β
k
d
k
,
where β
k
=
(
q
k
*
q
k
)
-
1
q
k
*
b
k.Proof.
See AppendixA.It can be deduced from Proposition1 that μ
¯
k
+
α
k is referred to as the next iteration of μ
¯
k. We apply a heuristic update sequence for the penalty factor to achieve a stable solution. If the kth iterative solution z
k is closer to the feasible region than the previous solution z
k
-
1, the penalty factor is decreased. Inversely, we increase the penalty factor when z
k is not closer to the feasible region.The specific steps of solving by the augmented Lagrange method are given as follows.Initialization: setk
=
1, z
1
∈
𝒞
K, σ
1
>
0, τ
∈
(
0,1
), and μ
¯
1
=
-
(
q
1
*
q
1
)
-
1
q
1
*
∇
f
1.(1)
Calculateα
k and d
k in terms of (33) and (34). If ∥
d
k
∥
≃
0, z
k is the KKT point of the problem (25); then stop iteration.
(2)
We use Armijo search method to find the maximum oft
k satisfying
(36)
p
σ
k
(
z
k
+
t
k
d
k
,
μ
¯
k
)
≤
p
σ
k
(
z
k
,
μ
¯
k
)
+
τ
t
k
∇
p
σ
k
*
(
z
k
,
μ
¯
k
)
d
k
.
(3)
Setz
k
+
1
=
z
k
+
t
k
d
k and update the multiplier factor and penalty factor, respectively:
(37)
μ
¯
k
+
1
=
μ
¯
k
+
α
k
,
σ
k
+
1
=
{
2
σ
k
if
∥
Π
(
z
k
+
1
)
∥
>
∥
Π
(
z
k
)
∥
,
1
2
σ
k
if
∥
Π
(
z
k
+
1
)
∥
≤
∥
Π
(
z
k
)
∥
.
(4)
Setk
=
k
+
1 and return to step (1).


## 5. Performance Analysis

### 5.1. Cramér-Rao Bound

In this subsection, the closed-form expression for the CRB of the proposed method with complex white Gaussian noise after the prewhitened process is illustrated. The bound of the noise variance estimation can be computed separately asCRB
n
=
1
/
N
M (see [35]). The remaining parameters consist of an unknown vector φ
=
[
θ
~
T
s
T
]
T. It is not an easy task to get the CRB of the unknown parameters. However, fortunately, [36, 37] have provided a critical inspiration for the derivation in this paper.First, the likelihood function is given by(38)
L
(
y
s
(
t
j
)
,
θ
~
)
=
1
(
2
π
)
M
N
(
1
/
2
)
M
N
·
exp
⁡
{
-
∑
j
=
1
M
‍
(
y
(
t
j
)
-
B
𝚽
A
s
(
t
j
)
)
*
·
exp
⁡
{
l
×
(
y
(
t
j
)
-
B
𝚽
A
s
(
t
j
)
)
∑
j
=
1
M
‍
}
.Thus, the log-likelihood function is(39)
ln
⁡
L
(
y
s
(
t
j
)
,
θ
~
)
=
const
-
∑
j
=
1
M
‍
(
y
(
t
j
)
-
B
𝚽
A
s
(
t
j
)
)
*
const
-
∑
j
=
1
M
l
‍
×
(
y
(
t
j
)
-
B
𝚽
A
s
(
t
j
)
)
.Then, the partial derivatives of (39) with respect to θ
~, s
r
(
t
j
)
=
Re
{
s
(
t
j
)
} and s
i
(
t
j
)
=
Im
⁡
{
s
(
t
j
)
} are given by
(40)
∂
ln
⁡
L
∂
θ
~
=
2
∑
j
=
1
M
‍
Re
{
S
j
*
U
*
𝚽
*
B
*
n
(
t
j
)
}
,
∂
ln
⁡
L
∂
s
r
(
t
j
)
=
2
Re
{
A
*
𝚽
*
B
*
n
(
t
j
)
}
,
∂
ln
⁡
L
∂
s
i
(
t
j
)
=
2
Im
⁡
{
A
*
𝚽
*
B
*
n
(
t
j
)
}
,
where S
j
=
diag
⁡
(
s
(
t
j
)
) and U
=
[
d
a
(
θ
~
1
)
/
d
θ
~
1
d
a
(
θ
~
2
)
/
d
θ
~
2
⋯
d
a
(
θ
~
K
)
/
d
θ
~
K
]. Following [35, 37], we can obtain the Fisher information matrix (FIM) as follows:
(41)
FIM
=
Λ
-
[
Δ
1
r
Δ
1
i
Δ
2
r
⋯
Δ
M
i
Δ
M
r
]
·
[
G
r
-
G
i
0
⋯
0
G
i
G
r
⋱
⋮
0
⋱
0
⋮
⋱
G
r
-
G
i
0
⋯
0
G
i
G
r
]
·
[
Δ
1
r
Δ
1
i
Δ
2
r
⋮
Δ
M
r
Δ
M
i
]
,
where
(42)
E
{
∂
ln
⁡
L
∂
θ
~
(
∂
ln
⁡
L
∂
θ
~
)
T
}
=
Λ
=
2
∑
j
=
1
M
‍
Re
{
S
j
*
U
*
𝚽
*
B
*
B
𝚽
U
S
j
}
,
E
{
∂
ln
⁡
L
∂
s
r
(
t
k
)
(
∂
ln
⁡
L
∂
s
i
(
t
p
)
)
T
}
=
-
G
i
-
1
=
-
Im
⁡
{
G
-
1
}
=
-
Im
⁡
{
2
A
*
𝚽
*
B
*
B
𝚽
A
}
δ
k
p
,
E
{
∂
ln
⁡
L
∂
s
r
(
t
k
)
(
∂
ln
⁡
L
∂
s
r
(
t
p
)
)
T
}
=
E
{
∂
ln
⁡
L
∂
s
i
(
t
k
)
(
∂
ln
⁡
L
∂
s
i
(
t
p
)
)
T
}
=
G
r
-
1
=
Re
{
G
-
1
}
=
Re
{
2
A
*
𝚽
*
B
*
B
𝚽
A
}
δ
k
p
,
E
{
∂
ln
⁡
L
∂
s
r
(
t
j
)
(
∂
ln
⁡
L
∂
θ
~
)
T
}
=
Re
{
Δ
j
}
=
Δ
j
r
=
Re
{
2
A
*
𝚽
*
B
*
B
𝚽
U
S
j
}
,
j
=
1,2
,
…
,
M
,
E
{
∂
ln
⁡
L
∂
s
i
(
t
j
)
(
∂
ln
⁡
L
∂
θ
~
)
T
}
=
Im
⁡
{
Δ
j
}
=
Δ
j
i
=
Im
⁡
{
2
A
*
𝚽
*
B
*
B
𝚽
U
S
j
}
,
j
=
1,2
,
…
,
M
.It is well known that(43)
[
Δ
r
T
Δ
i
T
]
[
G
r
-
G
i
G
i
G
r
]
[
Δ
r
Δ
i
]
=
Re
{
Δ
*
G
Δ
}
.It can be deduced from (41) and (43) that
(44)
CRB
(
θ
~
)
=
FI
M
-
1
=
1
2
{
∑
j
=
1
M
‍
Re
{
S
j
*
U
*
H
*
(
I
-
H
A
(
H
A
)
+
)
H
U
S
j
}
}
-
1
,
where H
=
B
𝚽 is a N
×
P matrix and (
·
)
+ denotes pseudoinverse. Note that the CRB in CS is affected not only by the conventional factors, for example, SNR, array structure and signal relative location, but also by the measurement matrix.


### 5.2. Performance Guarantee

In CS, the RIP has been deeply studied for the joint sparse recovery by minimizing thel
1 norm. We say that the matrix C
∈
𝒞
N
×
K obeys the RIP of the order L if there exists a constant δ
∈
(
0,1
) satisfying
(45)
(
1
-
δ
)
∥
s
∥
2
2
≪
∥
C
s
∥
2
2
≪
(
1
+
δ
)
∥
s
∥
2
2
.Therefore, all submatrices ofC with L columns are uniformly well conditioned. The restricted isometry constant (RIC) of the order L, described as δ
L
(
C
), is the smallest δ that satisfies (45) and δ
L
(
C
) satisfies
(46)
δ
L
(
C
)
=
max
⁡
|
J
′
|
=
L
∥
C
J
′
*
C
J
′
-
I
L
∥
,
where J
′ is a L-dimensional subsupport of J
=
[
1,2
,
…
,
K
] and C
J
′ denotes the submatrix of C with columns indexed by J
′. Note that the condition satisfying the RIP is so demanding that its applications are limited. Therefore, we should make use of a new version of the RIP, which is called the weak-1 RIP [38], to control the size of the recovery error. The weak-1 RIP is given in the following form:
(47)
(
1
-
α
)
∥
v
∥
2
2
≪
∥
C
v
∥
2
2
≪
(
1
+
α
)
∥
v
∥
2
2
for all v supported on U, where the cardinality of the set U is L
+
1. If the matrix C satisfies the weak-1 RIP, it can be deduced that 0
≤
α
≤
ϕ
L
+
1
(
C
U
), where ϕ
k
(
C
U
) denotes the kth largest eigenvalue of C
U. The corresponding weak-1 RIC is given by
(48)
α
L
+
1
w
(
C
;
U
)
=
min
⁡
ϕ
L
+
1
(
C
U
)
.In this paper, when the estimation quality is imperfect, especially in the unfavorable case, the support is no longer identified by the algebraic property ofR. Hence, a new performance guarantee is given in terms of the weak-1 RIP in the following proposition.Proposition 2.
One assumes thatF
=
R
^
1
/
2
+
(
1
/
2
)
R
^
-
1
/
2 and R
-
1
=
𝚿
R
s
-
1
(
θ
~
)
𝚿
* where 𝚿
=
[
ψ
1
ψ
2
…
ψ
K
]. θ
^ is an estimation that is asymptotically equal to θ
~
s such that |
tr
⁡
(
R
-
1
(
θ
^
)
)
-
tr
⁡
(
R
-
1
(
θ
~
s
)
)
|
≤
η for η
∈
(
0
,
1
/
2
∥
F
∥
2
2
). Let J
0 and J
1
=
{
j
i
}
i
=
1
L be the L-dimensional supports that consisted of the indexes of L largest elements in θ
~
s and θ
^, respectively. If the matrix R satisfies
(49)
α
L
+
1
w
(
R
;
J
0
)
>
κ
for
(50)
κ
≥
L
∥
F
∥
2
2
∑
j
∈
J
0
ψ
j
*
ψ
j
1
-
2
η
∥
F
∥
2
2
the support can be identified.Proof.
See AppendixB.It follows from Proposition2 that the performance guarantee of the proposed method requires a mild condition in the unfavorable case. However, in the favorable case, the performance guarantee only requires a much milder condition, α
L
+
1
w
(
R
;
J
0
)
>
0, which is an algebraic condition.


## 5.1. Cramér-Rao Bound

In this subsection, the closed-form expression for the CRB of the proposed method with complex white Gaussian noise after the prewhitened process is illustrated. The bound of the noise variance estimation can be computed separately asCRB
n
=
1
/
N
M (see [35]). The remaining parameters consist of an unknown vector φ
=
[
θ
~
T
s
T
]
T. It is not an easy task to get the CRB of the unknown parameters. However, fortunately, [36, 37] have provided a critical inspiration for the derivation in this paper.First, the likelihood function is given by(38)
L
(
y
s
(
t
j
)
,
θ
~
)
=
1
(
2
π
)
M
N
(
1
/
2
)
M
N
·
exp
⁡
{
-
∑
j
=
1
M
‍
(
y
(
t
j
)
-
B
𝚽
A
s
(
t
j
)
)
*
·
exp
⁡
{
l
×
(
y
(
t
j
)
-
B
𝚽
A
s
(
t
j
)
)
∑
j
=
1
M
‍
}
.Thus, the log-likelihood function is(39)
ln
⁡
L
(
y
s
(
t
j
)
,
θ
~
)
=
const
-
∑
j
=
1
M
‍
(
y
(
t
j
)
-
B
𝚽
A
s
(
t
j
)
)
*
const
-
∑
j
=
1
M
l
‍
×
(
y
(
t
j
)
-
B
𝚽
A
s
(
t
j
)
)
.Then, the partial derivatives of (39) with respect to θ
~, s
r
(
t
j
)
=
Re
{
s
(
t
j
)
} and s
i
(
t
j
)
=
Im
⁡
{
s
(
t
j
)
} are given by
(40)
∂
ln
⁡
L
∂
θ
~
=
2
∑
j
=
1
M
‍
Re
{
S
j
*
U
*
𝚽
*
B
*
n
(
t
j
)
}
,
∂
ln
⁡
L
∂
s
r
(
t
j
)
=
2
Re
{
A
*
𝚽
*
B
*
n
(
t
j
)
}
,
∂
ln
⁡
L
∂
s
i
(
t
j
)
=
2
Im
⁡
{
A
*
𝚽
*
B
*
n
(
t
j
)
}
,
where S
j
=
diag
⁡
(
s
(
t
j
)
) and U
=
[
d
a
(
θ
~
1
)
/
d
θ
~
1
d
a
(
θ
~
2
)
/
d
θ
~
2
⋯
d
a
(
θ
~
K
)
/
d
θ
~
K
]. Following [35, 37], we can obtain the Fisher information matrix (FIM) as follows:
(41)
FIM
=
Λ
-
[
Δ
1
r
Δ
1
i
Δ
2
r
⋯
Δ
M
i
Δ
M
r
]
·
[
G
r
-
G
i
0
⋯
0
G
i
G
r
⋱
⋮
0
⋱
0
⋮
⋱
G
r
-
G
i
0
⋯
0
G
i
G
r
]
·
[
Δ
1
r
Δ
1
i
Δ
2
r
⋮
Δ
M
r
Δ
M
i
]
,
where
(42)
E
{
∂
ln
⁡
L
∂
θ
~
(
∂
ln
⁡
L
∂
θ
~
)
T
}
=
Λ
=
2
∑
j
=
1
M
‍
Re
{
S
j
*
U
*
𝚽
*
B
*
B
𝚽
U
S
j
}
,
E
{
∂
ln
⁡
L
∂
s
r
(
t
k
)
(
∂
ln
⁡
L
∂
s
i
(
t
p
)
)
T
}
=
-
G
i
-
1
=
-
Im
⁡
{
G
-
1
}
=
-
Im
⁡
{
2
A
*
𝚽
*
B
*
B
𝚽
A
}
δ
k
p
,
E
{
∂
ln
⁡
L
∂
s
r
(
t
k
)
(
∂
ln
⁡
L
∂
s
r
(
t
p
)
)
T
}
=
E
{
∂
ln
⁡
L
∂
s
i
(
t
k
)
(
∂
ln
⁡
L
∂
s
i
(
t
p
)
)
T
}
=
G
r
-
1
=
Re
{
G
-
1
}
=
Re
{
2
A
*
𝚽
*
B
*
B
𝚽
A
}
δ
k
p
,
E
{
∂
ln
⁡
L
∂
s
r
(
t
j
)
(
∂
ln
⁡
L
∂
θ
~
)
T
}
=
Re
{
Δ
j
}
=
Δ
j
r
=
Re
{
2
A
*
𝚽
*
B
*
B
𝚽
U
S
j
}
,
j
=
1,2
,
…
,
M
,
E
{
∂
ln
⁡
L
∂
s
i
(
t
j
)
(
∂
ln
⁡
L
∂
θ
~
)
T
}
=
Im
⁡
{
Δ
j
}
=
Δ
j
i
=
Im
⁡
{
2
A
*
𝚽
*
B
*
B
𝚽
U
S
j
}
,
j
=
1,2
,
…
,
M
.It is well known that(43)
[
Δ
r
T
Δ
i
T
]
[
G
r
-
G
i
G
i
G
r
]
[
Δ
r
Δ
i
]
=
Re
{
Δ
*
G
Δ
}
.It can be deduced from (41) and (43) that
(44)
CRB
(
θ
~
)
=
FI
M
-
1
=
1
2
{
∑
j
=
1
M
‍
Re
{
S
j
*
U
*
H
*
(
I
-
H
A
(
H
A
)
+
)
H
U
S
j
}
}
-
1
,
where H
=
B
𝚽 is a N
×
P matrix and (
·
)
+ denotes pseudoinverse. Note that the CRB in CS is affected not only by the conventional factors, for example, SNR, array structure and signal relative location, but also by the measurement matrix.


## 5.2. Performance Guarantee

In CS, the RIP has been deeply studied for the joint sparse recovery by minimizing thel
1 norm. We say that the matrix C
∈
𝒞
N
×
K obeys the RIP of the order L if there exists a constant δ
∈
(
0,1
) satisfying
(45)
(
1
-
δ
)
∥
s
∥
2
2
≪
∥
C
s
∥
2
2
≪
(
1
+
δ
)
∥
s
∥
2
2
.Therefore, all submatrices ofC with L columns are uniformly well conditioned. The restricted isometry constant (RIC) of the order L, described as δ
L
(
C
), is the smallest δ that satisfies (45) and δ
L
(
C
) satisfies
(46)
δ
L
(
C
)
=
max
⁡
|
J
′
|
=
L
∥
C
J
′
*
C
J
′
-
I
L
∥
,
where J
′ is a L-dimensional subsupport of J
=
[
1,2
,
…
,
K
] and C
J
′ denotes the submatrix of C with columns indexed by J
′. Note that the condition satisfying the RIP is so demanding that its applications are limited. Therefore, we should make use of a new version of the RIP, which is called the weak-1 RIP [38], to control the size of the recovery error. The weak-1 RIP is given in the following form:
(47)
(
1
-
α
)
∥
v
∥
2
2
≪
∥
C
v
∥
2
2
≪
(
1
+
α
)
∥
v
∥
2
2
for all v supported on U, where the cardinality of the set U is L
+
1. If the matrix C satisfies the weak-1 RIP, it can be deduced that 0
≤
α
≤
ϕ
L
+
1
(
C
U
), where ϕ
k
(
C
U
) denotes the kth largest eigenvalue of C
U. The corresponding weak-1 RIC is given by
(48)
α
L
+
1
w
(
C
;
U
)
=
min
⁡
ϕ
L
+
1
(
C
U
)
.In this paper, when the estimation quality is imperfect, especially in the unfavorable case, the support is no longer identified by the algebraic property ofR. Hence, a new performance guarantee is given in terms of the weak-1 RIP in the following proposition.Proposition 2.
One assumes thatF
=
R
^
1
/
2
+
(
1
/
2
)
R
^
-
1
/
2 and R
-
1
=
𝚿
R
s
-
1
(
θ
~
)
𝚿
* where 𝚿
=
[
ψ
1
ψ
2
…
ψ
K
]. θ
^ is an estimation that is asymptotically equal to θ
~
s such that |
tr
⁡
(
R
-
1
(
θ
^
)
)
-
tr
⁡
(
R
-
1
(
θ
~
s
)
)
|
≤
η for η
∈
(
0
,
1
/
2
∥
F
∥
2
2
). Let J
0 and J
1
=
{
j
i
}
i
=
1
L be the L-dimensional supports that consisted of the indexes of L largest elements in θ
~
s and θ
^, respectively. If the matrix R satisfies
(49)
α
L
+
1
w
(
R
;
J
0
)
>
κ
for
(50)
κ
≥
L
∥
F
∥
2
2
∑
j
∈
J
0
ψ
j
*
ψ
j
1
-
2
η
∥
F
∥
2
2
the support can be identified.Proof.
See AppendixB.It follows from Proposition2 that the performance guarantee of the proposed method requires a mild condition in the unfavorable case. However, in the favorable case, the performance guarantee only requires a much milder condition, α
L
+
1
w
(
R
;
J
0
)
>
0, which is an algebraic condition.


## 6. Simulation Results

In this section, the performance of the proposed method is illustrated by several simulation results and compared with that of CS-MUISC, SPICE, and CRB under the condition that the number of sources is unknown. We consider the spatial signal impinging on the uniform linear array (ULA) with interspacingλ
/
2 where λ denotes the wavelength of source. In the ULA case, the steering vector corresponding to the DOA equal to θ
k is given by
(51)
a
(
θ
k
)
=
[
1
e
-
j
π
sin
(
θ
k
)
⋯
e
-
j
(
P
-
1
)
π
sin
(
θ
k
)
]
T
,
where the number of the array elements is set to be P
=
8.In the simulation, the average root mean square error (RMSE) of the DOA estimation with 50 Monte Carlo runs is defined as the significant performance index:(52)
RMSE
=
[
∑
m
=
1
50
‍
∑
l
=
1
L
‍
(
θ
¯
l
m
-
θ
l
)
2
50
L
]
1
/
2
,
where θ
¯
l
m is the estimate of θ
l in the mth run.The resolution of the grid is closely related to the precision of the DOA estimation. A coarse grid can lead to poor precision, but a too fine grid increases computational complexity. Therefore, an adaptive grid refinement method is used to balance the tradeoff between precision and computational complexity. In the simulation, we make a coarse grid with1
° step in the range of -
9
0
° to 9
0
° and perform a local fine grid in the vicinity of locations obtained by using the coarse grid.In the first simulation, we display the superimposed spatial spectra of three algorithms in 10 Monte Carlo runs in the scenario with low SNR, small number of snapshots, and five sources impinging from[
-
40
.
3
°
-
2
5
°
5.2
°
15.2
°
38.1
°
], respectively where two most closely spaced sources are correlated and the remaining sources are uncorrelated. The spatial spectra are shown in Figure 1 with 3 dB SNR and 50 snapshots. The following facts can be acquired from Figure 1 as follows: the spatial spectrum obtained by CS-MUSIC suffers from severe interference at the true directions, especially at two most closely spaced correlated sources whose bias is clearly seen from the insert. Two most closely spaced correlated sources can be resolved by SPICE (note that peaks in the spatial spectrum are much larger than 2, but they are cut off at 2 to use the same scale as the other figures in Figure 1), but SPICE can yield false peaks and slight bias in the vicinity of the correlated sources and uncorrelated sources, respectively. The proposed method AULMC yields a nearly ideal spatial spectrum and provides a precise estimation for all the sources. In summary, AULMC outperforms CS-MUSIC and SPICE in terms of the spatial spectrum in the scenario with low SNR, small number of snapshots, and closely spaced correlated sources.Superimposed spatial spectra of CS-MUSIC, SPICE, and AULMC in 10 Monte Carlo runs, where the red vertical dashed lines denote the true DOAs. (a) CS-MUSIC, (b) SPICE and (c) AULMC.
(a)
(b)
(c)We analyze the RMSE of three algorithms under different conditions in the second simulation. The source model is the same as the first simulation. Figure2 shows the RMSE as a function of SNR of all the algorithms and CRB in 50 Monte Carlo runs for the fixed number of snapshots 50, whereas the RMSE versus number of snapshots is shown in Figure 3 for the fixed SNR 5 dB in 50 Monte Carlo runs. Based on Figures 2 and 3, we can draw the conclusions that the RMSE of AULMC is smaller than those of other two algorithms and AULMC has the more significant performance advantages than the other two algorithms, especially in the scenarios with low SNR or small number of snapshots. One possible explanation is that AULMC can give the stable estimation in every Monte Carlo run. It can be also seen that the RMSE is close to the CRB with the increase of SNR and the number of snapshots.Figure 2
RMSE of the DOA estimation versus SNR for 50 snapshots.Figure 3
RMSE of the DOA estimation versus number of snapshots for 5 dB SNR.In Figure4, we display the relation between the RMSE and angel separation of correlated sources which can illustrate the resolving ability. Let two correlated sources at angles 2
0
° and 2
0
°
+
Δ
θ, where the step of Δ
θ is 1
°, be impinged on the ULA. The SNR is 13 dB and the number of snapshots is 100. It can be seen from Figure 4 that when angle separation is 2
°, AULMC fails; however, AULMC can still provide a precise estimation as long as the angle separation is no less than 3
° and has higher resolution than the other two algorithms.Figure 4
RMSE of the DOA estimation versus angle separation, where the SNR is 13 dB and the number of snapshots is 100.Finally, the computation time of different algorithms versus number of snapshots is shown in Table1 for comparing the efficiency of these algorithms and SDP solver. Two correlated sources impinge on the ULA at 2
0
° and 2
5
°. The SNR is fixed at 13 dB. The computation time is obtained by the MATLAB 7.8 (R2009a) on a 2.8 GHz 4 GB PC. For AULMC, the computation time is mainly spent on the iterations of augmented Lagrange.Table 1
Comparison of computation time.
Number of snapshots
Time (sec)
AULMC
CS-MUSIC
SPICE
SDP solver
50
0.0651
0.0569
0.0601
0.7526
75
0.0776
0.0604
0.0714
0.8871
100
0.0860
0.0640
0.0819
1.0268
125
0.0974
0.0691
0.0923
1.3018It can be seen from Table1 that the computation time of SDP solver is the longest, and although the computation time of AULMC is longer than that of other two algorithms, it is comparable. Moreover, it is worth noting that the performance of AULMC is much better than that of CS-MUSIC or SPICE.


## 7. Conclusion

A novel augmented Lagrange based on modified covariance matching criterion method for DOA estimation is proposed in CS. It is proved that the problem of minimizing the modified covariance matching criterion is an SDP, which can be transformed into the constrained quadratic programming problem solved by the augmented Lagrange method. A detailed derivation for the CRB and a theoretical performance guarantee for identifying the support are provided. Simulation results show that AULMC outperforms CS-MUSIC and SPICE in terms of the spatial spectrum and has more precise estimation as well as higher resolution, especially in the scenarios with low SNR, small number of snapshots, and closely spaced correlated sources.



---

*Source: 241469-2014-02-11.xml*