# Two-Dimensional DOA Estimation for Uniform Rectangular Array Using Reduced-Dimension Propagator Method

Authors: Ming Zhou; Xiaofei Zhang; Xiaofeng Qiu; Chenghua Wang
Journal: International Journal of Antennas and Propagation
 (2015)

---

## Abstract

A novel algorithm is proposed for two-dimensional direction of arrival (2D-DOA) estimation with uniform rectangular array using reduced-dimension propagator method (RD-PM). The proposed algorithm requires no eigenvalue decomposition of the covariance matrix of the receive data and simplifies two-dimensional global searching in two-dimensional PM (2D-PM) to one-dimensional local searching. The complexity of the proposed algorithm is much lower than that of 2D-PM. The angle estimation performance of the proposed algorithm is better than that of estimation of signal parameters via rotational invariance techniques (ESPRIT) algorithm and conventional PM algorithms, also very close to 2D-PM. The angle estimation error and Cramér-Rao bound (CRB) are derived in this paper. Furthermore, the proposed algorithm can achieve automatically paired 2D-DOA estimation. The simulation results verify the effectiveness of the algorithm.


---

## Body

## 1. Introduction

Direction-of-arrival (DOA) estimation is a fundamental problem in array signal processing and has been widely used in many fields [1–5]. Till now, many algorithms have been proposed for DOA estimation with uniform linear array. Among them, multiple signal classification (MUSIC) algorithm [6] and estimation of signal parameters via rotational invariance technique (ESPRIT) algorithm [7] are widely used superresolution methods. The two-dimensional DOA (2D-DOA) estimation problem, which plays an important role in array signal processing, has received more attention. This problem is usually considered with rectangular array, two parallel uniform linear arrays, L-shape array, and so forth. Also, many algorithms have been considered for 2D-DOA estimation, which include 2D-MUSIC algorithm [8], 2D Unitary ESPRIT algorithm [9], modified 2D-ESPRIT algorithm [10], matrix pencil methods [11, 12], maximum likelihood method [13], parallel factor (PARAFAC) algorithm [14], and high order cumulant method [15].Propagator method, which is known as a low complexity method without eigenvalue decomposition (EVD) of the covariance matrix of the received data, has been proposed for DOA estimation through peak searching [17–20]. Due to its low complexity, PM algorithms are widely used for 2D-DOA estimation. In [21], a PM-based DOA estimation algorithm is proposed for two parallel uniform linear arrays via the rotational invariance property of propagator matrix, which requires extra pairing match. Reference [22] presents an efficient 2D-DOA estimation algorithm with two parallel uniform linear arrays using PM. In [23], an improved PM algorithm is proposed for 2D-DOA estimation, which has better angle estimation performance than the algorithms in [21, 22]. The above-mentioned PM algorithms can be extended to the rectangular array for 2D-DOA estimation; however, they only employ the rotational invariance property of propagator matrices, and as a result of that, the accuracy of angle estimation performance is not sufficient in some cases, especially with low signal to noise ratio (SNR). Two-dimensional PM (2D-PM) algorithm through peak searching can be extended for 2D-DOA estimation; however, the high computational complexity caused by 2D peak searching makes it inefficient.In this paper, we derive a reduced-dimension PM (RD-PM) algorithm, which reduces the high complexity for 2D-DOA estimation with uniform rectangular array compared with 2D-PM algorithm. The proposed algorithm applies the rotational invariance property of propagator matrix to get the initial angle estimation and then employs one-dimensional local searching to get more accurate angle and finally obtains the other angle via least square (LS) method and estimate azimuth and elevation angles. The proposed algorithm has the following advantages:(
1
) the proposed algorithm has lower computational complexity than 2D-PM algorithm since it requires no EVD of the covariance matrix of the receive data and it only requires one-dimensional local searching; (
2
) the angle estimation performance of the proposed algorithm is better than that of ESPRIT algorithm and the PM algorithms in [21, 23], also very close to that of 2D-PM algorithm; (
3
) it can obtain automatically paired parameters estimation.The remainder of this paper is structured as follows: Section2 shows the data model for uniform rectangular array, while Section 3 proposes the RD-PM algorithm for 2D-DOA estimation. The angle estimation error is derived as well as CRB in Section 4. In Section 5, the simulation results verify the feasibility and effectiveness of the proposed algorithm, and the conclusions are showed in Section 6.Notion. ·
T, ·
H, ·
-
1, and ·
+ denote transpose, conjugate-transpose, inverse, and pseudoinverse operations, respectively; d
i
a
g
(
v
) stands for diagonal matrix whose diagonal element is a vector v; 1
P, I
P, and 0
P denote a P
×
1 vector of ones, a P
×
P identity matrix, and a P
×
1 vector of zeros, respectively. ⊗, ∘, and ⊙ are the Kronecker product, the Khatri-Rao product, and the Hadamard product, respectively; a
n
g
l
e
· is to get the phase angle.


## 2. Data Model

As illustrated in Figure1, consider a uniform rectangular array having M
×
N sensors; M and N are the numbers of sensors in x-axis and y-axis, respectively. The distance between the two adjacent elements is d in both x-axis and y-axis. Assume that there are K uncorrelated sources, and θ
k
,
ϕ
k are the elevation and azimuth angles of the kth source.Figure 1
The structure of uniform rectangular array [16].The received signal of the first subarrayx
1
(
t
) can be expressed as(1)
x
1
t
=
A
x
s
t
+
n
1
t
,where A
x
=
[
a
x
(
ϕ
1
,
θ
1
)
,
a
x
(
ϕ
2
,
θ
2
)
,
…
,
a
x
(
ϕ
k
,
θ
k
)
], a
x
(
ϕ
k
,
θ
k
)
=
[
1
,
e
-
j
2
π
d
s
i
n
θ
k
c
o
s
ϕ
k
/
λ
,
… , e
-
j
2
π
(
M
-
1
)
d
s
i
n
θ
k
c
o
s
ϕ
k
/
λ
]
T, and λ is the wavelength. n
1
(
t
) is the additive white Gaussian noise of the first subarray and s
(
t
)
∈
C
K
×
1 is the source vector.Similarly, the received signal of thenth subarray x
n
(
t
) in the rectangular array is also denoted as (2)
x
n
t
=
A
x
Φ
n
-
1
s
t
+
n
n
t
,where Φ
=
d
i
a
g
(
e
-
j
2
π
d
s
i
n
θ
1
s
i
n
ϕ
1
/
λ
,
…
,
e
-
j
2
π
d
s
i
n
θ
k
s
i
n
ϕ
k
/
λ
) and n
n
(
t
) is the additive white Gaussian noise of the nth subarray. So the integrated received signal of the rectangular array is shown as follows:(3)
x
t
=
x
1
t
x
2
t
⋮
x
N
t
=
A
x
A
x
Φ
⋮
A
x
Φ
N
-
1
s
t
+
n
1
t
n
2
t
⋮
n
N
t
,where x
i
(
t
)
(
i
=
1
,
…
,
N
) stands for the received signal of the ith subarray. Equation (3) is also denoted as(4)
x
t
=
A
y
∘
A
x
s
t
+
n
t
=
A
s
t
+
n
t
,where A
=
A
y
∘
A
x, A
y
=
[
a
y
(
ϕ
1
,
θ
1
)
,
a
y
(
ϕ
2
,
θ
2
)
,
…
,
a
y
(
ϕ
k
,
θ
k
)
], a
y
(
ϕ
k
,
θ
k
)
=
[
1
,
e
-
j
2
π
d
s
i
n
θ
k
s
i
n
ϕ
k
/
λ
,
…
,
e
-
j
2
π
(
M
-
1
)
d
s
i
n
θ
k
s
i
n
ϕ
k
/
λ
]
T, and n
(
t
)
=
[
n
1
(
t
)
T
,
n
2
(
t
)
T
,
…
,
n
N
(
t
)
T
]
T.For the signal model in (4), we collect J snapshots to get the covariance matrix R
^
x as(5)
R
^
x
=
1
J
∑
j
=
1
J
x
t
j
x
H
t
j
.


## 3. Two-Dimensional DOA Estimation Algorithm

### 3.1.2
D-PM Algorithm

In this section, we give a brief introduction of 2D-PM algorithm.Partition the matrixA as (6)
A
=
A
1
A
2
,where A
1
∈
C
K
×
K and A
2
∈
C
(
M
N
-
K
)
×
K. Since A
1 is a K
×
K nonsingular matrix, A
2 is a linear transformation of A
1 as [17](7)
A
2
=
P
H
A
1
,where P
∈
C
K
×
M
N
-
K is the propagator matrix. We define the matrix Q as (8)
Q
=
P
-
I
M
N
-
K
.According to (6) and (7), we have (9)
Q
H
A
=
0
M
N
-
K
×
K
.Then, we construct the function of the 2D-PM spatial spectrum as (10)
F
2D-PM
θ
,
ϕ
=
1
a
y
ϕ
⊗
a
x
θ
H
Q
Q
H
a
y
ϕ
⊗
a
x
θ
.By searching ϕ and θ, the K largest peaks of the costing function in (10) correspond to the angles information. However, the high computational cost caused by the 2D peak searching makes this algorithm inefficient.For estimating the propagator matrixP, we partition the covariance matrix R
^
x as(11)
R
^
x
=
G
^
,
H
^
,where G
^
∈
C
M
N
×
K
,
H
^
∈
C
M
N
×
M
N
-
K. The propagator matrix P
^ meets the minimization problem J
c
s
m
(
P
^
)
=
H
^
-
G
^
P
^
2, and the estimation of P
^ is (12)
P
^
=
G
^
H
G
^
-
1
G
^
H
H
^
.


### 3.2. Reduced-Dimension Propagator Method for 2D-DOA Estimation

#### 3.2.1. Initial Estimation

Define(13)
u
≜
s
i
n
θ
c
o
s
ϕ
,
v
≜
s
i
n
θ
s
i
n
ϕ
.Then, a
x
(
ϕ
,
θ
) and a
y
(
ϕ
,
θ
) can be rewritten as(14)
a
x
v
=
1
,
e
-
j
2
π
d
v
/
λ
,
…
,
e
-
j
2
π
N
-
1
d
v
/
λ
T
a
y
u
=
1
,
e
-
j
2
π
d
u
/
λ
,
…
,
e
-
j
2
π
M
-
1
d
u
/
λ
T
.Construct matrix P
c as (15)
P
c
=
I
K
P
H
.Divide P
c into N blocks as P
c
=
[
P
c
1
T
,
P
c
2
T
,
…
,
P
c
N
T
]
T, and every block contains M
×
K elements. Construct P
1
=
[
P
c
1
T
,
P
c
2
T
,
…
,
P
c
(
N
-
1
)
T
]
T and P
2
=
[
P
c
2
T
,
P
c
3
T
,
…
,
P
c
N
T
]
T; then we have(16)
P
2
=
P
1
T
-
1
Φ
T
,where Φ
=
d
i
a
g
(
e
-
j
2
π
d
u
1
/
λ
,
…
,
e
-
j
2
π
d
u
k
/
λ
). Use EVD of P
1
+
P
2 and obtain the kth eigenvalue as p
k; then the initial estimation of u
^
k
i
n
i is (17)
u
^
k
i
n
i
=
-
a
n
g
l
e
p
k
λ
2
π
d
.


#### 3.2.2. 2D-DOA Estimation

Define(18)
V
u
,
v
=
a
y
u
⊗
a
x
v
H
Q
Q
H
a
y
u
⊗
a
x
v
.Equation (18) can also be denoted as(19)
V
u
,
v
=
a
x
v
H
a
y
u
⊗
I
M
H
Q
Q
H
a
y
u
⊗
I
M
a
x
v
=
a
x
v
H
C
u
a
x
v
,where C
(
u
)
=
[
a
y
(
u
)
⊗
I
M
]
H
Q
Q
H
[
a
y
(
v
)
⊗
I
M
]. Equation (19) is the problem of quadratic optimization. We consider the constraint e
1
H
a
x
(
v
)
=
1, where e
1
=
[
1,0
,
…
,
0
]
T
∈
R
M
×
1. The optimization problem can be reconstructed via linearly constrained minimum variance as [24] (20)
m
i
n
u
,
v
a
x
v
H
C
u
a
x
v
,
s
.
t
.
l
l
l
l
l
l
e
1
H
a
x
v
=
1
.The costing function is (21)
L
u
,
v
=
a
x
v
H
C
u
a
x
v
-
ω
e
1
H
a
x
v
-
1
,where ω is a constant. We have (22)
∂
∂
a
x
v
L
u
,
v
=
2
C
u
a
x
v
+
ω
e
1
=
0
.According to (22), we obtain a
x
(
v
)
=
μ
C
-
1
(
u
)
e
1, where μ is also a constant. As e
1
H
a
x
(
v
)
=
1, μ
=
1
/
{
e
1
H
C
-
1
(
u
)
e
1
}. Then, a
x
(
v
) can be expressed as (23)
a
x
v
=
C
-
1
u
e
1
e
1
H
C
-
1
u
e
1
.Combining (20) and (23), the estimation of u is (24)
u
^
=
a
r
g
m
i
n
u
1
e
1
H
C
u
-
1
e
1
=
a
r
g
m
a
x
u
e
1
H
C
-
1
u
e
1
.Equation (24) also can be present as (25)
u
^
k
=
a
r
g
m
a
x
e
1
H
C
-
1
u
e
1
,
k
=
1,2
,
…
,
K
.We find the K largest peaks of the (
1,1
) element of C
-
1
(
u
) through searching u
∈
[
u
^
k
i
n
i
-
Δ
u
,
u
^
k
i
n
i
+
Δ
u
], where Δ
u is a small value. The K largest peaks (
u
1
,
u
2
,
…
,
u
K
) correspond to the estimations of s
i
n
θ
k
s
i
n
ϕ
k
(
k
=
1,2
,
…
,
K
); then we obtain Kvectors of a
^
x
(
v
1
)
,
a
^
x
(
v
2
)
,
…
,
a
^
x
(
v
K
) via (23). Since a
x
(
v
k
)
=
[
1
,
e
-
j
2
π
d
v
k
/
λ
,
…
,
e
-
j
2
π
(
N
-
1
)
d
v
k
/
λ
]
T, we have(26)
g
k
=
-
a
n
g
l
e
a
x
v
k
=
0
,
2
π
d
v
k
λ
,
…
,
2
π
M
-
1
d
v
k
λ
T
=
v
k
q
,where q
=
[
0,2
π
d
/
λ
,
…
,
2
π
(
M
-
1
)
d
/
λ
]
T. Normalize a
^
x
(
v
K
) and use LS method to estimate v
K; the LS fitting is (27)
m
i
n
c
k
E
c
k
-
g
^
k
F
2
,where g
^
k
=
-
a
n
g
l
e
(
a
^
x
(
v
k
)
), c
k
=
[
c
k
0
,
c
k
1
]
T
∈
R
2
×
1, E
=
[
1
M
,
q
], and c
k
1 is the estimation of v
k. The solution to c
k is (28)
c
k
0
,
v
k
T
=
E
T
E
-
1
E
T
g
^
k
.Finally, we get the estimation of θ
k
,
ϕ
k as (29)
θ
^
k
=
s
i
n
-
1
a
b
s
v
^
k
+
j
u
^
k
ϕ
^
k
=
a
n
g
l
e
v
^
k
+
j
u
^
k
.Till now, we have achieved the proposed algorithm for 2D-DOA estimation. The major steps and complexity of the algorithm are shown as follows:(1)
construct the covariance matrixR
^
x via (5);…, O
{
J
M
2
N
2
};
(2)
partition theR
^
x to get P; then construct Q
,
P
c
,
P
1
,
P
2 and obtain the initial estimation of u
^
k
i
n
i from the EVD of P
1
+
P
2;…, O
{
K
M
2
N
2
+
K
2
M
N
+
2
K
2
M
(
N
-
1
)
+
2
K
3
};
(3)
find theK largest peaks of the (
1,1
) element of Q
-
1
(
u
) through searching u
∈
[
u
^
k
i
n
i
-
Δ
u
,
u
^
k
i
n
i
+
Δ
u
] and obtain the estimation of u
^
k;…, O
{
n
1
K
[
(
M
2
N
+
M
2
)
(
M
N
-
K
)
+
M
2
]
};
(4)
computea
^
x
(
v
K
) via (23), and then v
^
k can be estimated through the LS fitting; …
,
O
{
6
M
+
12
};
(5)
estimateθ
^
k
,
ϕ
^
k via (29).Remark 1.
In this paper, we assume that the target numberK is a prior known parameter. If not, we can estimate K by applying the method in [25, 26].Remark 2.
The proposed algorithm gets the estimation ofu
^
k through searching u
∈
[
u
^
k
i
n
i
-
Δ
u
,
u
^
k
i
n
i
+
Δ
u
] and then obtains the corresponding a
^
x
(
v
K
) and estimates corresponding v
^
k. Since u
^
k and v
^
k are correspondence, the proposed algorithm can obtain automatically paired parameters estimation.Remark 3.
The proposed algorithm requiresu
k
,
v
k
(
k
=
1
,
…
,
K
) being different. If the sources have the same u or v, then the direction matrix A
=
A
x
∘
A
y is a matrix of rank loss according to the property of Khatri-Rao product [27]. That will lead many algorithms invalid, such as ESPRIT, PMs in [21, 23], and our algorithm.


### 3.3. Complexity Analysis

The proposed algorithm has much lower complexity than 2D-PM algorithm. The major complexity of the proposed algorithm isO
{
n
1
K
[
(
M
2
N
+
M
2
)
(
M
N
-
K
)
+
M
2
]
+
J
M
2
N
2
+
2
K
2
M
(
N
-
1
)
+
2
K
3
+
M
2
N
2
K
+
K
2
M
N
}, where n
1 is the number of steps within the local searching range; 2D-PM algorithm requires O
{
J
M
2
N
2
+
M
2
N
2
K
+
2
K
3
+
n
2
[
(
M
N
+
1
)
(
M
N
-
K
)
]
}, where n
2 is the number of steps within the global searching range and n
2
≫
n
1, while ESPRIT algorithm costs O
{
J
M
2
N
2
+
M
3
N
3
+
2
K
3
(
M
-
1
)
N
+
6
K
3
+
2
K
2
(
N
-
1
)
M
}.The complexity comparison with different parameters is shown in Figures2 and 3. From Figures 2 and 3, we can find that the proposed algorithm has much lower complexity than 2D-PM algorithm since the proposed algorithm only requires a one-dimensional local searching.Figure 2
Complexity comparison of different algorithms with different values ofM
/
N.Figure 3
Complexity comparison of different algorithms with different values ofK
/
J.


## 3.1.2
D-PM Algorithm

In this section, we give a brief introduction of 2D-PM algorithm.Partition the matrixA as (6)
A
=
A
1
A
2
,where A
1
∈
C
K
×
K and A
2
∈
C
(
M
N
-
K
)
×
K. Since A
1 is a K
×
K nonsingular matrix, A
2 is a linear transformation of A
1 as [17](7)
A
2
=
P
H
A
1
,where P
∈
C
K
×
M
N
-
K is the propagator matrix. We define the matrix Q as (8)
Q
=
P
-
I
M
N
-
K
.According to (6) and (7), we have (9)
Q
H
A
=
0
M
N
-
K
×
K
.Then, we construct the function of the 2D-PM spatial spectrum as (10)
F
2D-PM
θ
,
ϕ
=
1
a
y
ϕ
⊗
a
x
θ
H
Q
Q
H
a
y
ϕ
⊗
a
x
θ
.By searching ϕ and θ, the K largest peaks of the costing function in (10) correspond to the angles information. However, the high computational cost caused by the 2D peak searching makes this algorithm inefficient.For estimating the propagator matrixP, we partition the covariance matrix R
^
x as(11)
R
^
x
=
G
^
,
H
^
,where G
^
∈
C
M
N
×
K
,
H
^
∈
C
M
N
×
M
N
-
K. The propagator matrix P
^ meets the minimization problem J
c
s
m
(
P
^
)
=
H
^
-
G
^
P
^
2, and the estimation of P
^ is (12)
P
^
=
G
^
H
G
^
-
1
G
^
H
H
^
.


## 3.2. Reduced-Dimension Propagator Method for 2D-DOA Estimation

### 3.2.1. Initial Estimation

Define(13)
u
≜
s
i
n
θ
c
o
s
ϕ
,
v
≜
s
i
n
θ
s
i
n
ϕ
.Then, a
x
(
ϕ
,
θ
) and a
y
(
ϕ
,
θ
) can be rewritten as(14)
a
x
v
=
1
,
e
-
j
2
π
d
v
/
λ
,
…
,
e
-
j
2
π
N
-
1
d
v
/
λ
T
a
y
u
=
1
,
e
-
j
2
π
d
u
/
λ
,
…
,
e
-
j
2
π
M
-
1
d
u
/
λ
T
.Construct matrix P
c as (15)
P
c
=
I
K
P
H
.Divide P
c into N blocks as P
c
=
[
P
c
1
T
,
P
c
2
T
,
…
,
P
c
N
T
]
T, and every block contains M
×
K elements. Construct P
1
=
[
P
c
1
T
,
P
c
2
T
,
…
,
P
c
(
N
-
1
)
T
]
T and P
2
=
[
P
c
2
T
,
P
c
3
T
,
…
,
P
c
N
T
]
T; then we have(16)
P
2
=
P
1
T
-
1
Φ
T
,where Φ
=
d
i
a
g
(
e
-
j
2
π
d
u
1
/
λ
,
…
,
e
-
j
2
π
d
u
k
/
λ
). Use EVD of P
1
+
P
2 and obtain the kth eigenvalue as p
k; then the initial estimation of u
^
k
i
n
i is (17)
u
^
k
i
n
i
=
-
a
n
g
l
e
p
k
λ
2
π
d
.


### 3.2.2. 2D-DOA Estimation

Define(18)
V
u
,
v
=
a
y
u
⊗
a
x
v
H
Q
Q
H
a
y
u
⊗
a
x
v
.Equation (18) can also be denoted as(19)
V
u
,
v
=
a
x
v
H
a
y
u
⊗
I
M
H
Q
Q
H
a
y
u
⊗
I
M
a
x
v
=
a
x
v
H
C
u
a
x
v
,where C
(
u
)
=
[
a
y
(
u
)
⊗
I
M
]
H
Q
Q
H
[
a
y
(
v
)
⊗
I
M
]. Equation (19) is the problem of quadratic optimization. We consider the constraint e
1
H
a
x
(
v
)
=
1, where e
1
=
[
1,0
,
…
,
0
]
T
∈
R
M
×
1. The optimization problem can be reconstructed via linearly constrained minimum variance as [24] (20)
m
i
n
u
,
v
a
x
v
H
C
u
a
x
v
,
s
.
t
.
l
l
l
l
l
l
e
1
H
a
x
v
=
1
.The costing function is (21)
L
u
,
v
=
a
x
v
H
C
u
a
x
v
-
ω
e
1
H
a
x
v
-
1
,where ω is a constant. We have (22)
∂
∂
a
x
v
L
u
,
v
=
2
C
u
a
x
v
+
ω
e
1
=
0
.According to (22), we obtain a
x
(
v
)
=
μ
C
-
1
(
u
)
e
1, where μ is also a constant. As e
1
H
a
x
(
v
)
=
1, μ
=
1
/
{
e
1
H
C
-
1
(
u
)
e
1
}. Then, a
x
(
v
) can be expressed as (23)
a
x
v
=
C
-
1
u
e
1
e
1
H
C
-
1
u
e
1
.Combining (20) and (23), the estimation of u is (24)
u
^
=
a
r
g
m
i
n
u
1
e
1
H
C
u
-
1
e
1
=
a
r
g
m
a
x
u
e
1
H
C
-
1
u
e
1
.Equation (24) also can be present as (25)
u
^
k
=
a
r
g
m
a
x
e
1
H
C
-
1
u
e
1
,
k
=
1,2
,
…
,
K
.We find the K largest peaks of the (
1,1
) element of C
-
1
(
u
) through searching u
∈
[
u
^
k
i
n
i
-
Δ
u
,
u
^
k
i
n
i
+
Δ
u
], where Δ
u is a small value. The K largest peaks (
u
1
,
u
2
,
…
,
u
K
) correspond to the estimations of s
i
n
θ
k
s
i
n
ϕ
k
(
k
=
1,2
,
…
,
K
); then we obtain Kvectors of a
^
x
(
v
1
)
,
a
^
x
(
v
2
)
,
…
,
a
^
x
(
v
K
) via (23). Since a
x
(
v
k
)
=
[
1
,
e
-
j
2
π
d
v
k
/
λ
,
…
,
e
-
j
2
π
(
N
-
1
)
d
v
k
/
λ
]
T, we have(26)
g
k
=
-
a
n
g
l
e
a
x
v
k
=
0
,
2
π
d
v
k
λ
,
…
,
2
π
M
-
1
d
v
k
λ
T
=
v
k
q
,where q
=
[
0,2
π
d
/
λ
,
…
,
2
π
(
M
-
1
)
d
/
λ
]
T. Normalize a
^
x
(
v
K
) and use LS method to estimate v
K; the LS fitting is (27)
m
i
n
c
k
E
c
k
-
g
^
k
F
2
,where g
^
k
=
-
a
n
g
l
e
(
a
^
x
(
v
k
)
), c
k
=
[
c
k
0
,
c
k
1
]
T
∈
R
2
×
1, E
=
[
1
M
,
q
], and c
k
1 is the estimation of v
k. The solution to c
k is (28)
c
k
0
,
v
k
T
=
E
T
E
-
1
E
T
g
^
k
.Finally, we get the estimation of θ
k
,
ϕ
k as (29)
θ
^
k
=
s
i
n
-
1
a
b
s
v
^
k
+
j
u
^
k
ϕ
^
k
=
a
n
g
l
e
v
^
k
+
j
u
^
k
.Till now, we have achieved the proposed algorithm for 2D-DOA estimation. The major steps and complexity of the algorithm are shown as follows:(1)
construct the covariance matrixR
^
x via (5);…, O
{
J
M
2
N
2
};
(2)
partition theR
^
x to get P; then construct Q
,
P
c
,
P
1
,
P
2 and obtain the initial estimation of u
^
k
i
n
i from the EVD of P
1
+
P
2;…, O
{
K
M
2
N
2
+
K
2
M
N
+
2
K
2
M
(
N
-
1
)
+
2
K
3
};
(3)
find theK largest peaks of the (
1,1
) element of Q
-
1
(
u
) through searching u
∈
[
u
^
k
i
n
i
-
Δ
u
,
u
^
k
i
n
i
+
Δ
u
] and obtain the estimation of u
^
k;…, O
{
n
1
K
[
(
M
2
N
+
M
2
)
(
M
N
-
K
)
+
M
2
]
};
(4)
computea
^
x
(
v
K
) via (23), and then v
^
k can be estimated through the LS fitting; …
,
O
{
6
M
+
12
};
(5)
estimateθ
^
k
,
ϕ
^
k via (29).Remark 1.
In this paper, we assume that the target numberK is a prior known parameter. If not, we can estimate K by applying the method in [25, 26].Remark 2.
The proposed algorithm gets the estimation ofu
^
k through searching u
∈
[
u
^
k
i
n
i
-
Δ
u
,
u
^
k
i
n
i
+
Δ
u
] and then obtains the corresponding a
^
x
(
v
K
) and estimates corresponding v
^
k. Since u
^
k and v
^
k are correspondence, the proposed algorithm can obtain automatically paired parameters estimation.Remark 3.
The proposed algorithm requiresu
k
,
v
k
(
k
=
1
,
…
,
K
) being different. If the sources have the same u or v, then the direction matrix A
=
A
x
∘
A
y is a matrix of rank loss according to the property of Khatri-Rao product [27]. That will lead many algorithms invalid, such as ESPRIT, PMs in [21, 23], and our algorithm.


## 3.2.1. Initial Estimation

Define(13)
u
≜
s
i
n
θ
c
o
s
ϕ
,
v
≜
s
i
n
θ
s
i
n
ϕ
.Then, a
x
(
ϕ
,
θ
) and a
y
(
ϕ
,
θ
) can be rewritten as(14)
a
x
v
=
1
,
e
-
j
2
π
d
v
/
λ
,
…
,
e
-
j
2
π
N
-
1
d
v
/
λ
T
a
y
u
=
1
,
e
-
j
2
π
d
u
/
λ
,
…
,
e
-
j
2
π
M
-
1
d
u
/
λ
T
.Construct matrix P
c as (15)
P
c
=
I
K
P
H
.Divide P
c into N blocks as P
c
=
[
P
c
1
T
,
P
c
2
T
,
…
,
P
c
N
T
]
T, and every block contains M
×
K elements. Construct P
1
=
[
P
c
1
T
,
P
c
2
T
,
…
,
P
c
(
N
-
1
)
T
]
T and P
2
=
[
P
c
2
T
,
P
c
3
T
,
…
,
P
c
N
T
]
T; then we have(16)
P
2
=
P
1
T
-
1
Φ
T
,where Φ
=
d
i
a
g
(
e
-
j
2
π
d
u
1
/
λ
,
…
,
e
-
j
2
π
d
u
k
/
λ
). Use EVD of P
1
+
P
2 and obtain the kth eigenvalue as p
k; then the initial estimation of u
^
k
i
n
i is (17)
u
^
k
i
n
i
=
-
a
n
g
l
e
p
k
λ
2
π
d
.


## 3.2.2. 2D-DOA Estimation

Define(18)
V
u
,
v
=
a
y
u
⊗
a
x
v
H
Q
Q
H
a
y
u
⊗
a
x
v
.Equation (18) can also be denoted as(19)
V
u
,
v
=
a
x
v
H
a
y
u
⊗
I
M
H
Q
Q
H
a
y
u
⊗
I
M
a
x
v
=
a
x
v
H
C
u
a
x
v
,where C
(
u
)
=
[
a
y
(
u
)
⊗
I
M
]
H
Q
Q
H
[
a
y
(
v
)
⊗
I
M
]. Equation (19) is the problem of quadratic optimization. We consider the constraint e
1
H
a
x
(
v
)
=
1, where e
1
=
[
1,0
,
…
,
0
]
T
∈
R
M
×
1. The optimization problem can be reconstructed via linearly constrained minimum variance as [24] (20)
m
i
n
u
,
v
a
x
v
H
C
u
a
x
v
,
s
.
t
.
l
l
l
l
l
l
e
1
H
a
x
v
=
1
.The costing function is (21)
L
u
,
v
=
a
x
v
H
C
u
a
x
v
-
ω
e
1
H
a
x
v
-
1
,where ω is a constant. We have (22)
∂
∂
a
x
v
L
u
,
v
=
2
C
u
a
x
v
+
ω
e
1
=
0
.According to (22), we obtain a
x
(
v
)
=
μ
C
-
1
(
u
)
e
1, where μ is also a constant. As e
1
H
a
x
(
v
)
=
1, μ
=
1
/
{
e
1
H
C
-
1
(
u
)
e
1
}. Then, a
x
(
v
) can be expressed as (23)
a
x
v
=
C
-
1
u
e
1
e
1
H
C
-
1
u
e
1
.Combining (20) and (23), the estimation of u is (24)
u
^
=
a
r
g
m
i
n
u
1
e
1
H
C
u
-
1
e
1
=
a
r
g
m
a
x
u
e
1
H
C
-
1
u
e
1
.Equation (24) also can be present as (25)
u
^
k
=
a
r
g
m
a
x
e
1
H
C
-
1
u
e
1
,
k
=
1,2
,
…
,
K
.We find the K largest peaks of the (
1,1
) element of C
-
1
(
u
) through searching u
∈
[
u
^
k
i
n
i
-
Δ
u
,
u
^
k
i
n
i
+
Δ
u
], where Δ
u is a small value. The K largest peaks (
u
1
,
u
2
,
…
,
u
K
) correspond to the estimations of s
i
n
θ
k
s
i
n
ϕ
k
(
k
=
1,2
,
…
,
K
); then we obtain Kvectors of a
^
x
(
v
1
)
,
a
^
x
(
v
2
)
,
…
,
a
^
x
(
v
K
) via (23). Since a
x
(
v
k
)
=
[
1
,
e
-
j
2
π
d
v
k
/
λ
,
…
,
e
-
j
2
π
(
N
-
1
)
d
v
k
/
λ
]
T, we have(26)
g
k
=
-
a
n
g
l
e
a
x
v
k
=
0
,
2
π
d
v
k
λ
,
…
,
2
π
M
-
1
d
v
k
λ
T
=
v
k
q
,where q
=
[
0,2
π
d
/
λ
,
…
,
2
π
(
M
-
1
)
d
/
λ
]
T. Normalize a
^
x
(
v
K
) and use LS method to estimate v
K; the LS fitting is (27)
m
i
n
c
k
E
c
k
-
g
^
k
F
2
,where g
^
k
=
-
a
n
g
l
e
(
a
^
x
(
v
k
)
), c
k
=
[
c
k
0
,
c
k
1
]
T
∈
R
2
×
1, E
=
[
1
M
,
q
], and c
k
1 is the estimation of v
k. The solution to c
k is (28)
c
k
0
,
v
k
T
=
E
T
E
-
1
E
T
g
^
k
.Finally, we get the estimation of θ
k
,
ϕ
k as (29)
θ
^
k
=
s
i
n
-
1
a
b
s
v
^
k
+
j
u
^
k
ϕ
^
k
=
a
n
g
l
e
v
^
k
+
j
u
^
k
.Till now, we have achieved the proposed algorithm for 2D-DOA estimation. The major steps and complexity of the algorithm are shown as follows:(1)
construct the covariance matrixR
^
x via (5);…, O
{
J
M
2
N
2
};
(2)
partition theR
^
x to get P; then construct Q
,
P
c
,
P
1
,
P
2 and obtain the initial estimation of u
^
k
i
n
i from the EVD of P
1
+
P
2;…, O
{
K
M
2
N
2
+
K
2
M
N
+
2
K
2
M
(
N
-
1
)
+
2
K
3
};
(3)
find theK largest peaks of the (
1,1
) element of Q
-
1
(
u
) through searching u
∈
[
u
^
k
i
n
i
-
Δ
u
,
u
^
k
i
n
i
+
Δ
u
] and obtain the estimation of u
^
k;…, O
{
n
1
K
[
(
M
2
N
+
M
2
)
(
M
N
-
K
)
+
M
2
]
};
(4)
computea
^
x
(
v
K
) via (23), and then v
^
k can be estimated through the LS fitting; …
,
O
{
6
M
+
12
};
(5)
estimateθ
^
k
,
ϕ
^
k via (29).Remark 1.
In this paper, we assume that the target numberK is a prior known parameter. If not, we can estimate K by applying the method in [25, 26].Remark 2.
The proposed algorithm gets the estimation ofu
^
k through searching u
∈
[
u
^
k
i
n
i
-
Δ
u
,
u
^
k
i
n
i
+
Δ
u
] and then obtains the corresponding a
^
x
(
v
K
) and estimates corresponding v
^
k. Since u
^
k and v
^
k are correspondence, the proposed algorithm can obtain automatically paired parameters estimation.Remark 3.
The proposed algorithm requiresu
k
,
v
k
(
k
=
1
,
…
,
K
) being different. If the sources have the same u or v, then the direction matrix A
=
A
x
∘
A
y is a matrix of rank loss according to the property of Khatri-Rao product [27]. That will lead many algorithms invalid, such as ESPRIT, PMs in [21, 23], and our algorithm.


## 3.3. Complexity Analysis

The proposed algorithm has much lower complexity than 2D-PM algorithm. The major complexity of the proposed algorithm isO
{
n
1
K
[
(
M
2
N
+
M
2
)
(
M
N
-
K
)
+
M
2
]
+
J
M
2
N
2
+
2
K
2
M
(
N
-
1
)
+
2
K
3
+
M
2
N
2
K
+
K
2
M
N
}, where n
1 is the number of steps within the local searching range; 2D-PM algorithm requires O
{
J
M
2
N
2
+
M
2
N
2
K
+
2
K
3
+
n
2
[
(
M
N
+
1
)
(
M
N
-
K
)
]
}, where n
2 is the number of steps within the global searching range and n
2
≫
n
1, while ESPRIT algorithm costs O
{
J
M
2
N
2
+
M
3
N
3
+
2
K
3
(
M
-
1
)
N
+
6
K
3
+
2
K
2
(
N
-
1
)
M
}.The complexity comparison with different parameters is shown in Figures2 and 3. From Figures 2 and 3, we can find that the proposed algorithm has much lower complexity than 2D-PM algorithm since the proposed algorithm only requires a one-dimensional local searching.Figure 2
Complexity comparison of different algorithms with different values ofM
/
N.Figure 3
Complexity comparison of different algorithms with different values ofK
/
J.


## 4. Performance Analysis

This section aims at analyzing the theoretic estimation error of the proposed algorithm and deriving the Cramér-Rao bound (CRB).


### 4.1. Error Analysis

Assume that the receive data in (1) can be written as (30)
X
~
=
X
+
Δ
X
,where X
~ is the noisy data matrix and Δ
X is the additive perturbation matrix with mean zero and variance σ
2 statistically independent. Since R
^
x
=
1
/
J
∑
j
=
1
J
x
(
t
j
)
x
H
(
t
j
), we have(31)
R
^
x
=
X
+
Δ
X
X
+
Δ
X
H
J
=
X
X
H
J
+
Δ
X
Δ
X
H
J
=
R
x
+
Δ
R
x
,where R
x
=
X
X
H
/
J and Δ
R
x
=
Δ
X
Δ
X
H
/
J. According to (11), we also have (32)
G
^
=
G
+
Δ
G
,
H
^
=
H
+
Δ
H
.The estimation of propagator matrix P can be denoted as (33)
P
^
=
G
+
Δ
G
G
+
Δ
G
H
-
1
G
+
Δ
G
H
+
Δ
H
H
.With a first-order expansion in Δ
G and Δ
H, P
^ is rewritten as (34)
P
^
=
P
+
Δ
P
,where Δ
P
=
(
G
G
H
)
-
1
G
(
Δ
H
H
-
Δ
G
H
P
). We also obtain(35)
Q
^
=
P
^
-
I
M
N
-
K
=
Q
+
Δ
Q
,where Δ
Q
=
T
Δ
X
H
Q with T
=
-
(
G
G
H
)
-
1
G
0
M
N
-
K.SinceC
^
=
Q
^
Q
^
H
=
C
+
Δ
C, using first-order expansion we obtain (36)
Δ
C
=
Q
Δ
Q
H
+
Δ
Q
Q
H
.The proposed algorithm has the same cost function as 2D-PM algorithm, which is a reduced-dimension form of 2D-PM algorithm. For simplification, we derive the theoretical error of the proposed algorithm through the cost function of 2D-PM algorithm.The function from which we search for the minima in order to determine the estimates is(37)
W
u
,
v
=
a
y
u
⊗
a
x
v
H
C
^
H
a
y
u
⊗
a
x
v
.Define r
=
[
u
,
v
]
T, ∇
r
=
[
∂
/
∂
u
,
∂
/
∂
v
]
T, d
u
(
k
)
=
∂
a
k
/
∂
u
k, and d
v
(
k
)
=
∂
a
k
/
∂
v
k. Use r
i
=
[
u
i
,
v
i
] to denote the true angles corresponding to the ith target and r
^
i
=
[
u
^
i
,
v
^
i
] denotes the estimated ones. Then, (38)
∇
r
W
u
,
v
r
^
i
=
0
.After the Taylor expansion around r
i (use W instead of W
(
u
,
v
) for simplification),(39)
0
=
∇
r
W
r
i
+
∇
r
(
∇
r
W
)
T
r
i
ξ
r
^
i
-
r
i
,where r
i
ξ is a vector of some value on the line segment joining r
^
i and r
i.Define(40)
H
i
=
l
i
m
N
→
∞
∇
r
(
∇
r
W
)
T
r
i
.Under large sample, (
r
^
i
-
r
i
) is asymptotically zero-mean Gaussian distributed [17–19], so the associated asymptotic covariance matrix is given by(41)
Φ
i
k
=
l
i
m
N
→
∞
E
r
^
i
-
r
i
r
^
k
-
r
k
T
=
H
i
-
1
l
i
m
N
→
∞
E
∇
r
i
f
∇
r
k
f
T
r
i
,
r
k
H
k
-
1
.According to (36) and (41) and taking into account that a
i
H
Q
=
0, we can determine(42)
Φ
i
k
=
σ
4
2
w
i
w
k
d
v
H
i
P
A
⊥
d
v
i
ρ
i
ρ
i
d
u
H
i
P
A
⊥
d
u
i
×
R
e
a
i
H
T
T
H
a
k
d
u
H
k
P
A
⊥
d
u
i
d
v
H
k
P
A
⊥
d
u
i
d
u
H
k
P
A
⊥
d
v
i
d
v
H
k
P
A
⊥
d
v
i
×
d
v
H
k
P
A
⊥
d
v
k
ρ
k
ρ
k
d
u
H
k
P
A
⊥
d
u
k
,where w
(
i
)
=
d
u
H
(
i
)
P
A
⊥
d
u
(
i
)
d
v
H
(
i
)
P
A
⊥
d
v
(
i
)
-
1
/
4
d
u
H
(
i
)
P
A
⊥
d
v
(
i
)
+
d
v
H
(
i
)
P
A
⊥
d
u
(
i
)
2, P
A
⊥
=
I
-
A
(
A
H
A
)
-
1
A
H, and ρ
i
=
-
d
u
H
i
P
A
⊥
d
v
i
+
d
v
H
i
P
A
⊥
d
u
i
/
2.That is,E
[
(
∂
u
k
)
2
]
=
Φ
k
k
(
1,1
), E
[
(
∂
v
k
)
2
]
=
Φ
k
k
(
2,2
), and E
[
∂
v
k
∂
u
k
]
=
Φ
k
k
(
2,1
).According to (29), the mean square error (MSE) of azimuth and elevation can be derived as(43)
E
∂
ϕ
k
2
=
1
2
s
i
n
2
θ
k
E
∂
u
k
2
+
E
∂
v
k
2
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
-
R
e
E
∂
v
k
2
+
2
j
E
∂
v
k
∂
u
k
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
-
E
∂
u
k
2
e
-
j
2
ϕ
k
=
1
2
s
i
n
2
θ
k
c
o
s
2
ϕ
k
Φ
k
k
1,1
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
+
s
i
n
2
ϕ
k
Φ
k
k
2,2
-
s
i
n
2
ϕ
k
Φ
k
k
2,1
E
∂
θ
k
2
=
1
s
i
n
2
θ
k
c
o
s
2
θ
k
×
u
2
E
∂
u
2
+
v
2
E
∂
v
2
+
2
u
v
E
∂
v
∂
u
=
1
c
o
s
2
θ
k
s
i
n
2
ϕ
k
Φ
k
k
1,1
+
c
o
s
2
ϕ
k
Φ
k
k
2,2
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
+
s
i
n
2
ϕ
k
Φ
k
k
2,1

s
i
n
2
ϕ
k
Φ
k
k
.The comparison of theoretical estimation error and simulation results will be shown in Section 5.


### 4.2. CRB

According to [28], we can derive the CRB for uniform rectangular array (44)
C
R
B
=
σ
2
2
J
R
e
D
H
Π
A
⊥
D
⊙
P
T
-
1
,where D
=
[
∂
a
1
/
∂
θ
1
,
∂
a
2
/
∂
θ
2
,
…
,
∂
a
K
/
∂
θ
K
,
∂
a
1
/
∂
ϕ
1, ∂
a
2
/
∂
ϕ
2
,
…
,
∂
a
K
/
∂
ϕ
K
] and a
k is the kth column of A; P
=
P
s
P
s
P
s
P
s, P
^
s
=
1
/
J
∑
t
=
1
J
b
t
b
H
t, Π
A
⊥
=
I
M
N
-
A
A
H
A
-
1
A
H, and σ
2 is the covariance of the noise. ⊙ is the Hadamard product.


## 4.1. Error Analysis

Assume that the receive data in (1) can be written as (30)
X
~
=
X
+
Δ
X
,where X
~ is the noisy data matrix and Δ
X is the additive perturbation matrix with mean zero and variance σ
2 statistically independent. Since R
^
x
=
1
/
J
∑
j
=
1
J
x
(
t
j
)
x
H
(
t
j
), we have(31)
R
^
x
=
X
+
Δ
X
X
+
Δ
X
H
J
=
X
X
H
J
+
Δ
X
Δ
X
H
J
=
R
x
+
Δ
R
x
,where R
x
=
X
X
H
/
J and Δ
R
x
=
Δ
X
Δ
X
H
/
J. According to (11), we also have (32)
G
^
=
G
+
Δ
G
,
H
^
=
H
+
Δ
H
.The estimation of propagator matrix P can be denoted as (33)
P
^
=
G
+
Δ
G
G
+
Δ
G
H
-
1
G
+
Δ
G
H
+
Δ
H
H
.With a first-order expansion in Δ
G and Δ
H, P
^ is rewritten as (34)
P
^
=
P
+
Δ
P
,where Δ
P
=
(
G
G
H
)
-
1
G
(
Δ
H
H
-
Δ
G
H
P
). We also obtain(35)
Q
^
=
P
^
-
I
M
N
-
K
=
Q
+
Δ
Q
,where Δ
Q
=
T
Δ
X
H
Q with T
=
-
(
G
G
H
)
-
1
G
0
M
N
-
K.SinceC
^
=
Q
^
Q
^
H
=
C
+
Δ
C, using first-order expansion we obtain (36)
Δ
C
=
Q
Δ
Q
H
+
Δ
Q
Q
H
.The proposed algorithm has the same cost function as 2D-PM algorithm, which is a reduced-dimension form of 2D-PM algorithm. For simplification, we derive the theoretical error of the proposed algorithm through the cost function of 2D-PM algorithm.The function from which we search for the minima in order to determine the estimates is(37)
W
u
,
v
=
a
y
u
⊗
a
x
v
H
C
^
H
a
y
u
⊗
a
x
v
.Define r
=
[
u
,
v
]
T, ∇
r
=
[
∂
/
∂
u
,
∂
/
∂
v
]
T, d
u
(
k
)
=
∂
a
k
/
∂
u
k, and d
v
(
k
)
=
∂
a
k
/
∂
v
k. Use r
i
=
[
u
i
,
v
i
] to denote the true angles corresponding to the ith target and r
^
i
=
[
u
^
i
,
v
^
i
] denotes the estimated ones. Then, (38)
∇
r
W
u
,
v
r
^
i
=
0
.After the Taylor expansion around r
i (use W instead of W
(
u
,
v
) for simplification),(39)
0
=
∇
r
W
r
i
+
∇
r
(
∇
r
W
)
T
r
i
ξ
r
^
i
-
r
i
,where r
i
ξ is a vector of some value on the line segment joining r
^
i and r
i.Define(40)
H
i
=
l
i
m
N
→
∞
∇
r
(
∇
r
W
)
T
r
i
.Under large sample, (
r
^
i
-
r
i
) is asymptotically zero-mean Gaussian distributed [17–19], so the associated asymptotic covariance matrix is given by(41)
Φ
i
k
=
l
i
m
N
→
∞
E
r
^
i
-
r
i
r
^
k
-
r
k
T
=
H
i
-
1
l
i
m
N
→
∞
E
∇
r
i
f
∇
r
k
f
T
r
i
,
r
k
H
k
-
1
.According to (36) and (41) and taking into account that a
i
H
Q
=
0, we can determine(42)
Φ
i
k
=
σ
4
2
w
i
w
k
d
v
H
i
P
A
⊥
d
v
i
ρ
i
ρ
i
d
u
H
i
P
A
⊥
d
u
i
×
R
e
a
i
H
T
T
H
a
k
d
u
H
k
P
A
⊥
d
u
i
d
v
H
k
P
A
⊥
d
u
i
d
u
H
k
P
A
⊥
d
v
i
d
v
H
k
P
A
⊥
d
v
i
×
d
v
H
k
P
A
⊥
d
v
k
ρ
k
ρ
k
d
u
H
k
P
A
⊥
d
u
k
,where w
(
i
)
=
d
u
H
(
i
)
P
A
⊥
d
u
(
i
)
d
v
H
(
i
)
P
A
⊥
d
v
(
i
)
-
1
/
4
d
u
H
(
i
)
P
A
⊥
d
v
(
i
)
+
d
v
H
(
i
)
P
A
⊥
d
u
(
i
)
2, P
A
⊥
=
I
-
A
(
A
H
A
)
-
1
A
H, and ρ
i
=
-
d
u
H
i
P
A
⊥
d
v
i
+
d
v
H
i
P
A
⊥
d
u
i
/
2.That is,E
[
(
∂
u
k
)
2
]
=
Φ
k
k
(
1,1
), E
[
(
∂
v
k
)
2
]
=
Φ
k
k
(
2,2
), and E
[
∂
v
k
∂
u
k
]
=
Φ
k
k
(
2,1
).According to (29), the mean square error (MSE) of azimuth and elevation can be derived as(43)
E
∂
ϕ
k
2
=
1
2
s
i
n
2
θ
k
E
∂
u
k
2
+
E
∂
v
k
2
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
-
R
e
E
∂
v
k
2
+
2
j
E
∂
v
k
∂
u
k
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
-
E
∂
u
k
2
e
-
j
2
ϕ
k
=
1
2
s
i
n
2
θ
k
c
o
s
2
ϕ
k
Φ
k
k
1,1
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
+
s
i
n
2
ϕ
k
Φ
k
k
2,2
-
s
i
n
2
ϕ
k
Φ
k
k
2,1
E
∂
θ
k
2
=
1
s
i
n
2
θ
k
c
o
s
2
θ
k
×
u
2
E
∂
u
2
+
v
2
E
∂
v
2
+
2
u
v
E
∂
v
∂
u
=
1
c
o
s
2
θ
k
s
i
n
2
ϕ
k
Φ
k
k
1,1
+
c
o
s
2
ϕ
k
Φ
k
k
2,2
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
l
+
s
i
n
2
ϕ
k
Φ
k
k
2,1

s
i
n
2
ϕ
k
Φ
k
k
.The comparison of theoretical estimation error and simulation results will be shown in Section 5.


## 4.2. CRB

According to [28], we can derive the CRB for uniform rectangular array (44)
C
R
B
=
σ
2
2
J
R
e
D
H
Π
A
⊥
D
⊙
P
T
-
1
,where D
=
[
∂
a
1
/
∂
θ
1
,
∂
a
2
/
∂
θ
2
,
…
,
∂
a
K
/
∂
θ
K
,
∂
a
1
/
∂
ϕ
1, ∂
a
2
/
∂
ϕ
2
,
…
,
∂
a
K
/
∂
ϕ
K
] and a
k is the kth column of A; P
=
P
s
P
s
P
s
P
s, P
^
s
=
1
/
J
∑
t
=
1
J
b
t
b
H
t, Π
A
⊥
=
I
M
N
-
A
A
H
A
-
1
A
H, and σ
2 is the covariance of the noise. ⊙ is the Hadamard product.


## 5. Simulation and Analysis

### 5.1. Simulation Results

We present 1000 Monte Carlo simulations to assess the angle estimation performance of the RD-PM algorithm. Define root mean square error (RMSE):(45)
R
M
S
E
=
1
K
∑
k
=
1
K
1
1000
∑
l
=
1
1000
θ
^
k
,
l
-
θ
k
2
+
ϕ
^
k
,
l
-
ϕ
k
2
,where θ
^
k
,
l
,
ϕ
^
k
,
l are the estimates of θ
k
,
ϕ
k in the lth Monte Carlo trial, respectively. We assume that there are K
=
3 noncoherent sources with (
θ
1
,
ϕ
1
)
=
(
1
0
∘
,
1
5
∘
), (
θ
2
,
ϕ
2
)
=
(
2
0
∘
,
2
5
∘
), and (
θ
3
,
ϕ
3
)
=
(
3
0
∘
,
3
5
∘
); the spacing between two adjacent elements is d
=
λ
/
2. Note that M and N are the numbers of sensors in x
-axis and y-axis, respectively. J is the number of snapshots.Simulation 1.
Figure4 depicts angle estimation performance of the proposed algorithm with M
=
8, N
=
10, J
=
100, and S
N
R
=
20 dB. The results show that elevation and azimuth angles can be clearly observed.Figure 4
Angle estimation performance of the proposed algorithm with SNR = 20 dB.Simulation 2.
Figures5 and 6 present angle estimation performance comparison among the proposed algorithm, 2D-PM algorithm, the PM algorithm in [21], the PM algorithm in [23], ESPRIT algorithm, RD-MUSIC algorithm, CRB, and the theoretic estimation error of RD-PM. We can see that the estimation performance of the proposed algorithm is better than ESPRIT algorithm and the PM algorithms in [21, 23] and is very close to 2D-PM algorithm and RD-MUSIC algorithm. Compared with 2D-PM algorithm, the proposed algorithm has the same costing function and simplifies two-dimensional global searching to one-dimensional local searching. As a result of that, our algorithm has very close angle estimation performance to 2D-PM algorithm and greatly reduces complexity. Furthermore, Figures 5 and 6 also show that the RMSE of the proposed algorithm is almost the same as the theoretic estimation error.Figure 5
Angle estimation performance comparison (M
=
8, N
=
8, J
=
100, and K
=
3).Figure 6
Angle estimation performance comparison (M
=
8, N
=
10, J
=
120, and K
=
2).Simulation 3.
Figure7 presents angle estimation performance of the proposed algorithm with M
=
8, N
=
8, K
=
3, and different values of J. We can find that the angle estimation performance is improved when J increases. When J increases, we get more samples to estimate more accurate propagator matrix P, and the angle estimation performance is enhanced.Figure 7
Angle estimation performance with different values ofJ.Simulation 4.
Figure8 shows angle estimation performance of the proposed algorithm with M
=
8, N
=
8, J
=
100, and different values of K. It is shown that the angle estimation performance of the proposed algorithm is improved when K decreases. More sources enhance the interference, so the angle estimation performance of the proposed algorithm is affected.Figure 8
Angle estimation performance with different values ofK.Simulation 5.
Figures9 and 10 investigate angle estimation performance of the proposed algorithm with K
=
3, J
=
100, and different values of M or N,respectively. It is indicated that the performance of the proposed algorithm is improved when M or N increases. Multiple sensors enhance the aperture of the array as well as diversity gain, so the angle estimation performance is improved.Figure 9
Angle estimation performance withN
=
8 and different values of M.Figure 10
Angle estimation performance withM
=
8 and different values of N.Simulation 6.
Figure11 depicts angle estimation performance of the proposed algorithm when the two sources are closely spaced with M
=
8, N
=
8, J
=
100, and SNR = 20 dB. It is shown in Figure 11 that our algorithm works well for the closely spaced sources.Figure 11
Angle estimation performance of the proposed algorithm with the closely spaced sources.


## 5.1. Simulation Results

We present 1000 Monte Carlo simulations to assess the angle estimation performance of the RD-PM algorithm. Define root mean square error (RMSE):(45)
R
M
S
E
=
1
K
∑
k
=
1
K
1
1000
∑
l
=
1
1000
θ
^
k
,
l
-
θ
k
2
+
ϕ
^
k
,
l
-
ϕ
k
2
,where θ
^
k
,
l
,
ϕ
^
k
,
l are the estimates of θ
k
,
ϕ
k in the lth Monte Carlo trial, respectively. We assume that there are K
=
3 noncoherent sources with (
θ
1
,
ϕ
1
)
=
(
1
0
∘
,
1
5
∘
), (
θ
2
,
ϕ
2
)
=
(
2
0
∘
,
2
5
∘
), and (
θ
3
,
ϕ
3
)
=
(
3
0
∘
,
3
5
∘
); the spacing between two adjacent elements is d
=
λ
/
2. Note that M and N are the numbers of sensors in x
-axis and y-axis, respectively. J is the number of snapshots.Simulation 1.
Figure4 depicts angle estimation performance of the proposed algorithm with M
=
8, N
=
10, J
=
100, and S
N
R
=
20 dB. The results show that elevation and azimuth angles can be clearly observed.Figure 4
Angle estimation performance of the proposed algorithm with SNR = 20 dB.Simulation 2.
Figures5 and 6 present angle estimation performance comparison among the proposed algorithm, 2D-PM algorithm, the PM algorithm in [21], the PM algorithm in [23], ESPRIT algorithm, RD-MUSIC algorithm, CRB, and the theoretic estimation error of RD-PM. We can see that the estimation performance of the proposed algorithm is better than ESPRIT algorithm and the PM algorithms in [21, 23] and is very close to 2D-PM algorithm and RD-MUSIC algorithm. Compared with 2D-PM algorithm, the proposed algorithm has the same costing function and simplifies two-dimensional global searching to one-dimensional local searching. As a result of that, our algorithm has very close angle estimation performance to 2D-PM algorithm and greatly reduces complexity. Furthermore, Figures 5 and 6 also show that the RMSE of the proposed algorithm is almost the same as the theoretic estimation error.Figure 5
Angle estimation performance comparison (M
=
8, N
=
8, J
=
100, and K
=
3).Figure 6
Angle estimation performance comparison (M
=
8, N
=
10, J
=
120, and K
=
2).Simulation 3.
Figure7 presents angle estimation performance of the proposed algorithm with M
=
8, N
=
8, K
=
3, and different values of J. We can find that the angle estimation performance is improved when J increases. When J increases, we get more samples to estimate more accurate propagator matrix P, and the angle estimation performance is enhanced.Figure 7
Angle estimation performance with different values ofJ.Simulation 4.
Figure8 shows angle estimation performance of the proposed algorithm with M
=
8, N
=
8, J
=
100, and different values of K. It is shown that the angle estimation performance of the proposed algorithm is improved when K decreases. More sources enhance the interference, so the angle estimation performance of the proposed algorithm is affected.Figure 8
Angle estimation performance with different values ofK.Simulation 5.
Figures9 and 10 investigate angle estimation performance of the proposed algorithm with K
=
3, J
=
100, and different values of M or N,respectively. It is indicated that the performance of the proposed algorithm is improved when M or N increases. Multiple sensors enhance the aperture of the array as well as diversity gain, so the angle estimation performance is improved.Figure 9
Angle estimation performance withN
=
8 and different values of M.Figure 10
Angle estimation performance withM
=
8 and different values of N.Simulation 6.
Figure11 depicts angle estimation performance of the proposed algorithm when the two sources are closely spaced with M
=
8, N
=
8, J
=
100, and SNR = 20 dB. It is shown in Figure 11 that our algorithm works well for the closely spaced sources.Figure 11
Angle estimation performance of the proposed algorithm with the closely spaced sources.


## 6. Conclusions

In this paper, we have presented a novel algorithm for 2D-DOA estimation for uniform rectangular array using RD-PM algorithm. The proposed algorithm, which only requires one-dimensional local searching and avoids the EVD of the covariance matrix of the receive data, has a much lower complexity than 2D-PM algorithm. Simulation results show that the angle estimation performance of the proposed algorithm is better than that of ESPRIT algorithm and conventional PM algorithms, also very close to that of 2D-PM algorithm. Furthermore, the proposed algorithm can achieve paired 2D-DOA estimation automatically.



---

*Source: 485351-2015-05-24.xml*