# Multiple Sparse Measurement Gradient Reconstruction Algorithm for DOA Estimation in Compressed Sensing

Authors: Weijian Si; Xinggen Qu; Yilin Jiang; Tao Chen
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

A novel direction of arrival (DOA) estimation method in compressed sensing (CS) is proposed, in which the DOA estimation problem is cast as the joint sparse reconstruction from multiple measurement vectors (MMV). The proposed method is derived through transforming quadratically constrained linear programming (QCLP) into unconstrained convex optimization which overcomes the drawback thatl
1-norm is nondifferentiable when sparse sources are reconstructed by minimizing l
1-norm. The convergence rate and estimation performance of the proposed method can be significantly improved, since the steepest descent step and Barzilai-Borwein step are alternately used as the search step in the unconstrained convex optimization. The proposed method can obtain satisfactory performance especially in these scenarios with low signal to noise ratio (SNR), small number of snapshots, or coherent sources. Simulation results show the superior performance of the proposed method as compared with existing methods.


---

## Body

## 1. Introduction

Direction of arrival (DOA) estimation of multiple narrowband sources is an important research topic in array signal processing. It has been extensively studied in acoustic source localization, radar, and medical imaging [1–3]. Many effective DOA estimation algorithms have been proposed and developed, which mainly include beamforming algorithms such as MVDR [4] and subspace-based algorithms such as MUSIC [5]. To obtain preferable estimation performance, the Nyquist sampling theorem must be used in these conventional methods of data acquisition. However, high-speed sampling rate can impose so enormous pressure on capturing and storing data that requirements on both hardware and software are increased. Moreover, these methods suffer from serious performance degradation in these scenarios with low signal to noise ratio (SNR), small number of snapshots, or coherent sources.Recently, many applications involving compressed sensing (CS) [6–8], especially DOA estimation, have been attracting tremendous research interest in the signal processing. CS is an emerging area, and it can not only capture and store compressed or sparse sources simultaneously at a rate much lower than the Nyquist sampling rate but also reconstruct original sources using nonadaptive linear projection measurements onto a suitable measurement matrix, which satisfies the restricted isometry property (RIP) [9–11]. The sparse reconstruction aims to find the support which is shared by the unknown sparse vectors from multiple measurement vectors (MMV). The support denotes the indices of the nonzero elements in the unknown sparse vectors.CS has been widely applied to DOA estimation since sources are sparse in the spatial domain which results from the fact that there are much fewer true sources directions than all potential directions. Stoica et al. [12] proposed a sparse iterative covariance-based estimation method (SPICE) for array processing which is semiparametric estimation method and can avoid parameter selection. Hyder and Mahata [13] proposed an alternative strategy called joint l
0 approximation (JLZA-DOA) algorithm based on spatial sparsity, which can resolve closely spaced and high correlated sources even if the number of sources is unknown. Figueiredo et al. [14] proposed gradient projection algorithm to solve the bound-constrained quadratic programming formulation. Although it can be simply implemented, the regularization parameter is selected with difficulty and it is only suitable for the single measurement vector (SMV) model, which limits its practical engineering application.In this paper, we propose a novel multiple sparse measurement gradient reconstruction method called MSMGR for DOA estimation in CS. The method is derived through transforming quadratically constrained linear programming (QCLP) into unconstrained convex optimization to overcome the drawback thatl
1-norm is nondifferentiable when minimizing the l
1-norm for the sparse reconstruction. The search steepest descent step [15] and Barzilai-Borwein step [16] are alternately used as the search step to improve the convergence rate and estimation performance significantly. Furthermore, the singular value decomposition (SVD) is incorporated into the proposed method to reduce the computational complexity and the sensitivity to the noise. The proposed method is suitable for both SMV and MMV, and it has higher estimation accuracy and resolution than existing methods especially in these scenarios with low SNR, small number of snapshots, or coherent sources. Simulation results show the superior performance of the proposed method as compared with existing methods.


## 2. Problem Formulation

ConsiderK narrowband far-field sources, s
~
k
t, k
=
1,2
,
…
,
K, impinging on the sensor array consisting of N omnidirectional sensors from different directions, θ
~
i, i
=
1,2
,
…
,
K. The array observation model at tth snapshot can be formulated as(1)
x
t
=
∑
k
=
1
K
a
θ
~
k
s
~
k
t
+
n
t
,where n
(
t
)
∈
C
N
×
1 is a complex Gaussian white noise vector with zero mean and covariance matrix σ
2
I and a
(
θ
k
) is the N
×
1 steering vector of the source from the direction θ
k. Although the DOA estimation based on the single snapshot, which is a typical 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 typical MMV model.In order to cast the DOA estimation as a sparse reconstruction, letθ
k
k
=
1
P denote a fine enough grid which covers the entire spatial domain where there are P
(
P
≫
K
) potential directions of the sources so that the true directions θ
~
k
k
=
1
K are aligned or are close to the grids. It means that there exists θ
k
1
,
θ
k
2
,
…
,
θ
k
K being equal to θ
~
1
,
θ
~
2
,
…
,
θ
~
K, respectively. Thus, we have(2)
s
k
=
s
~
l
k
=
k
l
(
l
=
1,2
,
…
,
K
)
0
elsewhere
.The multiple snapshots model can be written as the following sparse form:(3)
x
t
=
∑
k
=
1
P
a
θ
k
s
k
t
+
n
t
=
A
s
t
+
n
t
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
h
i
h
h
t
=
1,2
,
…
,
L
,where L is the number of snapshots and Α
=
a
θ
1
a
θ
2
⋯
a
(
θ
P
) is the N
×
P array manifold matrix corresponding to all the potential directions which is also defined as an overcomplete dictionary in CS. s
t
=
s
1
(
t
)
s
2
(
t
)
⋯
s
P
(
t
)
T is the P
×
1 sparse vector with K nonzero elements at positions corresponding to the true directions and zero elements at the remaining P
-
K positions, where ·
T denotes the transpose operation. Hence, the matrix S
=
s
(
1
)
s
(
2
)
⋯
s
(
L
)
∈
C
P
×
L has K nonzero rows, that is, row K-sparse, since s
(
t
)
t
=
1
L share the common support. Obviously, the DOA estimation problem of multiple snapshots is that of identifying the row support of the unknown matrix S from the matrix Y
∈
C
M
×
L which is given by(4)
Y
=
y
(
1
)
y
(
2
)
⋯
y
(
L
)
=
Φ
Α
S
+
Φ
N
=
Θ
S
+
Φ
Nwith sensing matrix Θ, noise matrix N, and a common measurement matrix Φ of the size M
×
N with M
<
N where M is the number of nonadaptive linear projection measurements.It is well known that sparse sources can be reconstructed by solving thel
0-norm minimization problem. However, the optimization problem is nonconvex and the optimization method is both numerically unstable and computationally unaccepted [17]. Then, this problem is transformed into the l
1-norm minimization problem [18] so that we can accurately reconstruct the matrix S by solving the following QCLP problem:(5)
min
⁡
S
s
(
2
)
1
s
.
t
.
Y
-
Θ
S
2
2
≤
σ
2
,where s
(
2
) is the P
×
1 unknown sparse vector and ·
2 represents the Frobenius norm of matrices or the Euclidean norm of vectors. The ith entry of s
(
2
) is equal to Euclidean norm of the ith row of S; that is, s
i
(
2
)
=
S
(
i
,
:
)
2.


## 3. DOA Estimation

The SVD is employed on the matrixY to reduce the computational complexity and the sensitivity to the noise. Hence, we have(6)
Y
=
U
Σ
V
H
=
U
S
U
N
Σ
V
H
,where V is the orthogonal matrix and ·
H denotes the conjugate transpose operation. U
S and U
N denote the signal subspace and noise subspace, respectively. The M eigenvalues of the matrix U are arranged from the largest to the smallest; that is, σ
1
≥
σ
2
≥
⋯
σ
K
≫
σ
K
+
1
≥
⋯
≥
σ
M, where K large eigenvalues are dominant. U
S and U
N, respectively, consist of left singular eigenvectors corresponding to K big eigenvalues and M
-
K small eigenvalues. Denote Y
S
=
U
S
=
U
Σ
Β
=
Y
V
Β
∈
C
M
×
K, S
S
=
S
Σ
Β
, and N
S
=
N
Σ
Β where Β
=
I
K
O
K
×
L
-
K
H, I
K is an identity matrix of the size K
×
K, and O
K
×
(
L
-
K
) is a zero matrix of the size K
×
(
L
-
K
) so that we have (7)
Y
S
=
Θ
S
S
+
Φ
N
S
.As one may note, the dimension of the matrixY is reduced from M
×
L to M
×
K which can significantly reduce the computational complexity and the sensitivity to the noise especially in the scenario with large number of snapshots. The essence of the dimension reduction is to keep the signal subspace and discard the noise subspace. By using the SVD decomposition, (5) can be rewritten as the following form:(8)
min
⁡
S
s
S
(
2
)
1
s
.
t
.
Y
S
-
Θ
S
S
2
2
≤
σ
2
,where s
S
(
2
) is also a sparse vector and shares the same support with s
(
2
). To overcome the drawback that l
1-norm is nondifferentiable by minimizing the l
1-norm for solving the sparse sources, (8) is transformed into the unconstrained convex optimization by using Lagrange multiplication [19, 20]:(9)
s
^
S
(
2
)
=
arg
min
⁡
1
2
Y
S
-
Θ
S
S
2
2
+
τ
s
S
(
2
)
1
,where τ is nonnegative and it is called regularized factor that can serve as a tradeoff between the ability of suppressing noise and source sparsity. Note that the search path, which is obtained by projecting the negative gradient of the objective function in (9) onto the feasible set, cannot perform a backtracking line search well. Therefore, we adopt the m
1-norm of matrices to change the search direction which results in(10)
S
^
S
=
arg
min
⁡
1
2
Y
S
-
Θ
S
S
2
2
+
τ
S
S
m
1
,where ·
m
1 denotes the m
1-norm of matrices. A detailed derivation process of using MSMGR for DOA estimation is shown as follows. Assume that r
Γ
k
-
1 denotes the residual of the kth iteration, Γ
k denotes the support of the kth iteration, Θ
Γ
k denotes the submatrix of Θ with columns indexed by Γ
k, and S
^
Γ
k denotes the reconstructed source after the kth iteration. Therefore, the objective function of the kth iteration can be written as(11)
f
S
^
Γ
k
=
1
2
Y
S
-
Θ
Γ
k
S
^
Γ
k
2
2
+
τ
S
^
Γ
k
m
1
.The purpose of the current iteration is to find the sparse reconstructed sourceS
^
Γ
k which can minimize the objective function f
(
S
^
Γ
k
); that is, the residual is minimum after the current iteration. The expansion of (11) can be expressed as(12)
f
S
^
Γ
k
=
1
2
tr
⁡
Y
S
H
Y
S
-
Y
S
H
Θ
Γ
k
S
^
Γ
k
-
S
^
Γ
k

H
Θ
Γ
k
H
Y
S
h
h
h
h
+
S
^
Γ
k

H
Θ
Γ
k
H
Θ
Γ
k
S
^
Γ
k
+
τ
∑
i
=
1
k
∑
l
=
1
K
s
(
i
,
l
)
,where s
(
i
,
l
) denotes the element of the ith row and lth column of the matrix S
^
Γ
k. With further derivation, the minimum of the objective function (12) is equal to the maximum of (13) which can be given as the following form:(13)
f
S
^
Γ
k
=
1
2
tr
⁡
Y
S
H
Θ
Γ
k
S
^
Γ
k
+
S
^
Γ
k

H
Θ
Γ
k
H
Y
S
h
h
h
h
-
S
^
Γ
k

H
Θ
Γ
k
H
Θ
Γ
k
S
^
Γ
k
-
τ
∑
i
=
1
k
∑
l
=
1
K
s
i
,
l
.Based on the properties of matrix trace, (13) can be further simplified to(14)
tr
⁡
Y
S
H
Θ
Γ
k
Z
^
Γ
k
-
1
2
S
^
Γ
k

H
Θ
Γ
k
H
Θ
Γ
k
S
^
Γ
k
-
τ
∑
i
=
1
k
∑
l
=
1
K
s
i
,
l
.Then, the negative gradient is obtained by the partial derivative of (14) with respect to S
^
Γ
k which is given by(15)
d
k
=
Θ
Γ
k
H
Y
S
-
Θ
Γ
k
H
Θ
Γ
k
S
^
Γ
k
-
τ
d
∑
i
=
1
k
∑
l
=
1
K
s
(
i
,
l
)
d
S
^
Γ
k
=
Θ
Γ
k
H
Y
S
-
Θ
Γ
k
H
Θ
Γ
k
S
^
Γ
k
-
τ
D
,where D is referred to as polarity matrix that can judge the polarity of nonzero elements: (16)
d
(
i
,
l
)
=
1
+
j
Re
{
s
(
i
,
l
)
}
>
0
Im
{
s
(
i
,
l
)
}
>
0
1
-
j
Re
{
s
(
i
,
l
)
}
>
0
Im
{
s
(
i
,
l
)
}
<
0
-
1
+
j
Re
{
s
(
i
,
l
)
}
<
0
Im
{
s
(
i
,
l
)
}
>
0
-
1
-
j
Re
{
s
(
i
,
l
)
}
<
0
Im
{
s
(
i
,
l
)
}
<
0
,where d
(
i
,
l
) denotes the element of the ith row and lth column of the matrix D. Since the conventional search step is too small based on the orthogonality, the steepest descent step and Barzilai-Borwein step are alternately exploited as the search step μ
k in order to improve the convergence rate and estimation performance. Then, we have (17)
μ
k
=
μ
SD
k
is
even
μ
BB
k
is
odd
,where μ
SD and μ
BB are the steepest descent step and Barzilai-Borwein step, respectively. The specific steps of MSMGR are given as follows.Initialization. Set the number of iterations k
=
1, r
Γ
0
=
Y
S, Γ
0
=
⌀, Θ
Γ
0
=
⌀, and S
^
Γ
0
=
⌀.Step 1.
Calculate the inner product of the residual of thekth iteration and sensing matrix. Then, update the index p:(18)
p
=
arg
max
⁡
p
r
k
-
1
,
Θ
:
,
p
2
p
=
1,2
,
3
,
…
,
P
.Step 2.
Update the supportΓ
k
=
Γ
k
-
1
∪
p and the corresponding submatrix Θ
Γ
k
=
Θ
Γ
k
-
1
∪
Θ
(
:
,
p
).Step 3.
Calculate the negative gradientd
k and the search step μ
k in terms of (15) and (17), respectively.Step 4.
Update the constructed sourceS
^
Γ
k
=
S
^
Γ
k
-
1
+
d
k
·
μ
k in terms of the negative gradient and the search step.Step 5.
Update the polarity matrixD
=
[
sign
d
k
·
μ
k
;
O
1
×
K
], where O
1
×
K is the zero matrix of the size 1
×
K.Step 6.
Update the residualr
k
=
r
k
-
1
-
Θ
Γ
k
·
d
k
·
μ
k. If the residual satisfies the stopping criterion, stop the iteration; otherwise, set k
=
k
+
1 and return to Step 1.The core of the new method is that of updating the polarity matrix by zero-padding process in Step5 since the dimensions of the support and corresponding submatrix are both expanded in Step 2. Moreover, zero-padding process can guarantee the precision of the DOA estimation. The spectrum of the proposed method is obtained by estimating constructed source power from all potential directions. Like other spectral-based methods, the true directions are estimated by the locations of the highest peaks of the spectrum.


## 4. Simulation Results

In this section, the superior performance of the proposed method is shown as compared with existing JLZA-DOA and MUSIC methods by several numerical simulations. Consider the spatial sources 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 source from the direction θ
k is given by(19)
a
θ
k
=
1
e
-
j
π
sin
⁡
θ
k
⋯
e
-
j
N
-
1
π
sin
⁡
θ
k
T
,where the number of array elements is set to N
=
16. In the simulation, the regularized factor can be chosen as suggested in [21](20)
τ
=
0.1
Θ
H
Y
S
m
∞
.Following [22], it is easy to see that the unique minimum of (10) is the zero matrix for τ
>
Θ
H
Y
S
m
∞. In the simulation, the average root mean square error (RMSE) of the DOA estimation is defined as the significant performance index:(21)
RMSE
=
∑
j
=
1
J
∑
k
=
1
K
θ
^
k
j
-
θ
k
2
J
K
1
/
2
,where J is the number of independent Monte Carlo runs and θ
^
k
j is the estimate of θ
k in the jth run. The resolution of the grid is closely related to the precision of the DOA estimation. A coarse grid leads to poor precision, but a too fine grid increases computational complexity. Therefore, an adaptive grid refinement method is used for the tradeoff between precision and computational complexity. In the simulation, we set a coarse grid with 1
∘ step in the range of -
9
0
∘ to 9
0
∘ and make a local fine grid in the vicinity of the estimated angle.In the first simulation, we show the spatial spectra of three methods in the scenario with low SNR, small number of snapshots, and two uncorrelated sources impinging fromθ
1
=
0
∘ and θ
2
=
2
∘, respectively. The spatial spectra are shown in Figure 1 with SNR = 3 dB and 50 snapshots. The following conclusion can be acquired from Figure 1 that since the spatial spectra obtained by MUSIC and JLZA-DOA have only one peak, MUSIC and JLZA-DOA cannot identify the closely spaced sources accurately. However, the proposed method MSMGR has a nearly ideal spectrum and a precise estimation for the closely spaced sources. Therefore, MSMGR outperforms JLZA-DOA and MUSIC in terms of the spatial spectrum.Figure 1
Spatial spectra of MUSIC, JLZA- DOA, and MSMGR for uncorrelated sources.The RMSEs of three methods are analyzed under different conditions in the second simulation. Consider three sources impinging fromθ
1
=
0
∘
θ
2
=
49
.
5
∘ and θ
3
=
57
.
3
∘, respectively, where the latter two closely spaced sources are correlated and the first source is uncorrelated to them. The forward spatial smoothing method is exploited on the MUSIC called FSS-MUSIC to resolve the correlated sources. Figure 2 shows the RMSE as a function of SNR of three methods in 100 Monte Carlo runs for the fixed number of 50 snapshots whereas the RMSE versus the number of snapshots is shown in Figure 3 for the fixed SNR −6 dB in 100 Monte Carlo runs. The conclusions can be drawn from Figures 2 and 3 that the RMSE of MSMGR is smaller than those of other two methods and MSMGR has better estimation performance than the other two methods, especially in the scenarios with low SNR or small number of snapshots. The reason is that MSMGR can overcome the nondifferentiable drawback and exploit the alternate search step to improve the convergence rate and estimation performance. Moreover, the RMSE of MSMGR is close to those of the other two methods 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 −6 dB SNR.In Figure4, the relation between the RMSE and angle separation of correlated sources is shown, which can illustrate the resolving capability. Let two correlated sources at directions 2
0
∘ and 2
0
∘
+
Δ
θ, where the step of Δ
θ is 1
∘, be impinged on the ULA. It can be seen from Figure 4 that MSMGR suffers from serious performance degradation if the angle separation is 4
∘. However, MSMGR can still provide the most precise estimation as long as the angle separation is no less than 5
∘. Simulation results show that MSMGR has higher resolution than the other two methods.Figure 4
RMSE of the DOA estimation versus angle separation.Finally, we compare the computation time of different methods versus number of snapshots in Table1. Two correlated sources located at 2
0
∘ and 2
7
∘ impinge on the ULA. MSMGR-nosvd denotes that the SVD is not adopted in the process of employing MSMGR for DOA estimation.Table 1
Computation time (sec) of methods.
Estimation methods
Number of snapshots
50
75
100
125
MSMGR
0.685
0.724
0.741
0.779
JLZA-DOA
0.627
0.623
0.645
0.648
FSS-MUSIC
0.423
0.431
0.446
0.451
MSMGR-nosvd
1.521
1.566
1.618
1.654It is easy to see from Table1 that the computation time of MSMGR-nosvd is the longest and the computation time of MSMGR is longer than that of JLZA-DOA and FSS-MUSIC, but it is important to note that the performance of MSMGR is much better than that of these two methods. Moreover, it is proved that the SVD can significantly reduce the computation time.


## 5. Conclusion

In this paper, a novel MSMGR method for DOA estimation is proposed in CS. The proposed method is obtained by transforming QCLP into unconstrained convex optimization to overcome the drawback thatl
1-norm is nondifferentiable when sparse sources are reconstructed by minimizing the l
1-norm. An alternate search step is used to improve the convergence rate and estimation performance. The SVD is used to reduce the computational complexity and the sensitivity to the noise. Simulation results show that MSMGR outperforms JLZA-DOA and MUSIC in terms of the spatial spectrum and has more precise estimation as well as higher resolution; in particular when SNR is low, the number of snapshots is small and sources are coherent.



---

*Source: 152570-2015-08-26.xml*