# Adaptive Beamforming Based on Compressed Sensing with

Authors: Yubing Han; Jian Wang
Journal: International Journal of Antennas and Propagation
 (2015)

---

## Abstract

An adaptive beamforming based on compressed sensing with smoothedl
0 norm for large-scale sparse receiving array is proposed in this paper. Because of the spatial sparsity of the arriving signal, compressed sensing is applied to sample received signals with a sparse array and reduced channels. The signal of full array is reconstructed by using a compressed sensing reconstruction method based on smoothed l
0 norm. Then an iterative linearly constrained minimum variance beamforming algorithm is adopted to form antenna beam, whose main lobe is steered to the desired direction and nulls to the directions of interferences. Simulation results and Monte Carlo analysis for linear and planar arrays show that the beam performances of our proposed adaptive beamforming are similar to those of full array antenna.


---

## Body

## 1. Introduction

Adaptive beamforming is widely used in array signal processing for enhancing a desired signal while suppressing interference and noise at the output of an array of sensors [1, 2]. In bistatic radar systems, digital beamforming receiving array antennas are urgently needed so that the receiving antenna beams can cover the transmitting antenna beam flexibly. On receiving radar stations, to obtain high antenna gain and high angular measurement accuracy, an antenna array with a large number of elements should be used. So the high cost is a major drawback for a large-scale array. Meanwhile, the computational burden and high-rate data transmission are two bottlenecks in the implementation of an adaptive beamforming algorithm. To reduce the number of radio frequency (RF) front-ends without reducing the antenna aperture, we can use sparse arrays. Compared to the fully populated array, sparse arrays are attractive because they can reduce the number of elements and receiver channels, weight, power consumption, and cost. The sparse equally spaced arrays produce grating lobes inevitably and reduce the scanning range of the beam. In order to avoid grating lobes, sparse arrays are usually designed to be unequally spaced. There are many techniques that have been proposed over the last fifty years to design and synthesize such unequally spaced sparse arrays. The most common is stochastic methods, such as genetic algorithm [3], particle swarm [4], ant colony [5], and simulated annealing [6]. Matrix pencil methods have lately been efficiently applied to reconstruct focused and shaped beam patterns while reducing the number of the array elements [7]. Differential evolution [8] and sparse-periodic hybrid array [9] are also good methods for reducing the side-lobe levels. In [10], a procedure to synthesize sparse arrays with antenna selected via sequential convex optimization is presented. In [11], a method of pattern synthesis with sparse arrays based on Bayesian compressive sampling is proposed. However, these methods only optimize static patterns and are difficult to guarantee the performance of the beam patterns when we need to suppress interference adaptively.Recently, Candes and Donoho reported a novel sampling theory called compressed sensing (CS), also known as compressive sampling [12, 13]. CS theory asserts that one can recover certain signals from far fewer samples or measurements than traditional methods use. To make this possible, CS relies on one principle: sparsity, which pertains to the signals of interest. The theory states that as long as the signal is sparse or compressible, you can reconstruct the original signal from a small number of samples by solving an optimization problem with high probability. The sparsity property of signals has been utilized in a variety of applications including astronomy [14], optical imaging [15], radar [16, 17], and spectrum sensing [18]. In the literatures [19–23], some CS-based DOA estimations have been investigated and compared with each other. In [24], a method of DOA estimation through Bayesian compressive sensing is proposed. To obtain the sparsest solution, we may search for a solution with l
0 norm, that is, a minimum number of nonzero components [13]. It is well known that searching the minimum l
0 norm is an intractable combinatorial problem as the dimension increases, and it is too sensitive to noise [25]. To overcome this drawback, we can relax the problem by substituting the l
1 norm for the l
0 norm. One of the most successful approaches is basis pursuit (BP) [26], which finds the solution by linear programming (LP) methods. However, for the large systems of equations, this algorithm is very slow. Another family of algorithms is iterative reweighted least squares, with FOCUSS as an important member [27]. These methods are faster than BP, but their estimation quality is worse, especially if the number of nonzero elements of the sparsest solution is not very small. Another approach is the greedy algorithms such as matching pursuit (MP) [28] and orthogonal matching pursuit (OMP) [29, 30]; these methods are very fast, but do not provide good estimation of the sources. Different from the above approaches, the literatures [31, 32] proposed an algorithm to minimize the l
0 norm directly. By choosing a continuous function to approximate l
0 norm, the algorithm searches the optimal solution based on gradient ascent method. It is shown that this algorithm is about two to three orders of magnitude faster than BP (based on the interior-point LP solvers), while resulting in the same or better accuracy [32].An adaptive digital beamforming based on CS with smoothedl
0 norm for large-scale sparse receiving array is proposed in this paper. Due to the spatial sparseness of the arriving signal, CS theory is adopted to sample received signals with a sparse array, and the signal of full array is reconstructed by using the CS reconstruction method based on smoothed l
0 norm. Then an iterative linearly constrained minimum variance (LCMV) is used to form antenna beam, whose main lobe is steered to the desired direction and nulls are steered to the directions of interferences. For the proposed adaptive beamforming, it greatly reduces the number of elements and has similar performance to full array, which means that the beams have low side-lobes, deep nulls in the directions of interferences, and no grating lobes. Different from most CS-based DOA estimation in array processing applications, our proposed method applies the CS algorithm to reconstruct the signal of full array and then uses adaptive beamforming to form antenna beam. It cares of the accuracy of the reconstructed signals more than the accuracy of DOA estimations.The paper is organized as follows. In Section2, the signal model and spatial sparsity are discussed. The proposed adaptive beamforming based on CS with smoothed l
0 norm is developed in Section 3. In Section 4, when targets are not on the grid, a coarse-to-fine method is proposed to improve the signal reconstruction results. The numerical simulations under different situations are provided to illustrate the effectiveness of the algorithm in Section 5, and conclusions are given in Section 6.


## 2. Signal Model

Consider an array antenna withN elements and assume that K narrowband signals are incident to the antenna array. One of these signals is the desired signal, and the others are interferences. The received signal can be expressed by(1)
x
t
=
∑
k
=
1
K
‍
s
k
t
a
k
+
v
t
,where x
(
t
)
=
[
x
1
(
t
)
,
…
,
x
N
(
t
)
]
T, v
t is the N
×
1 white Gaussian noise, and s
k
(
t
) and a
k are the complex amplitude and steering vector of the kth incident signal. For an ideal uniform linear array with the spacing of d, we have(2)
a
k
=
1
,
e
j
2
π
d
sin
(
θ
k
)
/
λ
,
…
,
e
j
2
π
(
N
-
1
)
d
sin
(
θ
k
)
/
λ
T
,where λ is the carrier wavelength and θ
k represents the DOA of signal k. For a planar array, the steering vector can be expressed by(3)
a
k
=
e
j
2
π
(
ξ
1
u
k
+
η
1
v
k
)
/
λ
,
…
,
e
j
2
π
(
ξ
N
u
k
+
η
N
v
k
)
/
λ
T
,where u
k
=
sin
(
θ
k
)
cos
⁡
(
φ
k
), v
k
=
sin
(
θ
k
)
sin
(
φ
k
), and (
θ
k
,
φ
k
) are the elevation and azimuth angles and (
ξ
k
,
η
k
) are the coordinates of the kth element. For the 3D array with element position of (
ξ
k
,
η
k
,
ρ
k
), the steering vector is given by(4)
a
k
=
e
j
2
π
(
ξ
1
u
k
+
η
1
v
k
+
ρ
1
b
k
)
/
λ
,
…
,
e
j
2
π
(
ξ
N
u
k
+
η
N
v
k
+
ρ
N
b
k
)
/
λ
T
,where u
k
=
sin
(
θ
k
)
cos
⁡
(
φ
k
), v
k
=
sin
(
θ
k
)
sin
(
φ
k
), and b
k
=
cos
⁡
(
θ
k
).Without loss of generality, we only consider the linear array in this paper unless noted otherwise. For the DOAs in linear array, we assume that there is a partition grid of the entire angle space[
-
9
0
°
,
9
0
°
], such that(5)
-
9
0
°
=
θ
1
<
θ
2
<
⋯
<
θ
Γ
=
9
0
°
,
sin
θ
i
+
1
-
sin
θ
i
=
2
Γ
-
1
,where Γ is the number of grid points. We define a transformation matrix H as(6)
H
=
a
θ
1
,
a
θ
2
,
…
,
a
θ
Γ
.So the received signal can be rewritten as(7)
x
t
=
H
s
t
+
v
t
,where s
(
t
)
=
[
0,0
,
…
,
s
1
(
t
)
,
0
,
…
,
0
,
s
K
(
t
)
,
0
,
…
,
0
]
T is a Γ
×
1 complex vector with K nonzero elements. The indices of nonzero elements in s
(
t
) represent the DOAs of incident signals. In general, the number of incident signal is small; that means s
t
0
=
K
≪
N, so it is reasonable to assume that s
(
t
) is spatially sparse.Next we design aM
×
N dimension measurement matrix Φ with M
≪
N, and the measurement y
(
t
) is obtained by(8)
y
t
=
Φ
x
t
,where M represents the number of receiver channels, and it is smaller than the number of array elements. The measurement matrix Φ means a method of compressed sampling space signals. Substituting (7) into (8) results in(9)
y
t
=
Φ
H
s
(
t
)
+
Φ
v
(
t
)
=
P
s
t
+
Φ
v
t
,where P
=
Φ
H is called the observation matrix. Theoretical studies in [12, 13, 33] have shown that we can accurately reconstruct the “sparse” signal s
(
t
) from the “compressed” vector y
(
t
) when the observation matrix P satisfies the condition of restricted isometry property (RIP). For the selection of Φ, Candes et al. point out that if you want to reconstruct the signal completely, the observation matrix P must guarantee that two different K-sparse signals will not be mapped to the same collection [12]. Here we define a correlation coefficients matrix G with elements of(10)
g
i
j
=
P
i
H
P
j
P
i
2
×
P
j
2
,
i
,
j
=
1
,
…
,
M
,where P
i is the ith row vector of P. For reconstructing the signal completely, we should make the correlation coefficients g
i
j as small as possible. With given H, the optimal correlation coefficients and measurement matrix Φ can be optimized by using genetic algorithm [3], particle swarm optimization [4], or other methods [5, 6]. With the optimization of element positions, it can help reducing the recovery errors. In this paper, we do not discuss the selection of Φ and just randomly select M rows from an identity matrix of size N
×
N for convenience. That is to say, we randomly select M elements from the full array to sample the space signals.Owing to the spatial sparsity of incident signals, we can obtain the estimation ofs
(
t
) from measurement y
(
t
) by solving the CS problem(11)
s
^
(
t
)
=
argmin
s
(
t
)
s
(
t
)
0
s
.
t
.
y
(
t
)
-
P
s
(
t
)
2
2
≤
ɛand reconstruct the received signal using x
^
(
t
)
=
H
s
^
(
t
), where ·
0 is l
0 norm of a vector which stands for the number of nonzeros in a vector; ɛ stands for the noise level. From [25], we can know that solving (11) is both numerically unstable and NP-complete. To solve this problem, we do not relax (11) by substituting the l
1 norm for the l
0 norm such as BP [26], FOCUSS [27], and OMP [29]. The method we present in this paper is based on direct minimization of the l
0 norm. The basic idea is to choose a smoothed continuous function to approximate l
0 norm and search the optimal solution based on gradient ascent and projection method.


## 3. Proposed Adaptive Beamforming

### 3.1. Smoothedl
0 Norm [31, 32]

It is known thatl
0 norm of a vector is a discontinuous function; searching the minimum l
0 norm solution is a combinatorial problem, which is nonconvex and highly nonsmoothed. Furthermore, l
0 norm of a vector is extremely sensitive to noise, and its value will be completely changed even though the noise is very weak. In this section, we introduce a family of smooth complex approximators of l
0 norm, whose minimization can be implemented using gradient based methods. Consider a complex zero-mean Gaussian family of functions with the parameter σ
(12)
f
σ
c
=
exp
⁡
-
c
2
2
σ
2
=
exp
⁡
-
c
r
2
+
c
i
2
2
σ
2
,where |
c
|, c
r, and c
i represent the module and real and imaginary parts of c
∈
C, we have(13)
lim
σ
→
0
f
σ
(
c
)
=
1
,
c
=
0
0
,
c
≠
0or approximately(14)
f
σ
(
c
)
≈
1
,
c
≪
σ
0
,
c
≫
σ
.Then, for a complex vector s
=
(
s
1
,
s
2
,
…
,
s
Γ
)
T, we can define(15)
F
σ
s
=
∑
i
=
1
Γ
‍
f
σ
s
i
.Since the number of entries in s
(
t
) is Γ and the function F
σ
[
s
(
t
)
] is an indicator of the number of zero-entries in s
(
t
), the l
0 norm of s
(
t
) can be approximated by(16)
s
t
0
≈
Γ
-
F
σ
s
tfor small values of σ, and the approximation tends to equality when σ
→
0.Substituting this approximation into (11) yields the CS problem with smoothed l
0 norm(17)
s
^
(
t
)
=
argmin
s
t
Γ
-
F
σ
s
t
=
argmax
s
t
F
σ
s
t
s
.
t
.
y
t
-
P
s
t
2
2
≤
ɛ
.


### 3.2. Compressed Sensing Reconstruction Algorithm

For small values ofσ in (17), F
σ
[
s
(
t
)
] contains a lot of local maxima. Consequently, it is very difficult to directly maximize this function for very small values of σ. However, as the value of σ grows, the function becomes smoother and smoother, and, for sufficiently large values of σ, there is no local maxima. Similar to [31, 32], our approach to solve the problem (17) is then to decrease the value of σ gradually. The underlying thought is to select a σ which ensures that the initial optimization problem is convex and gradually increase the accuracy of the approximation. By careful selection of the sequence of σ, nonconvexity and thereby local minima can be effectively avoided. For each σ, we use some iterations of the steepest ascent algorithm (followed by projection onto the feasible set) to maximize F
σ
[
s
(
t
)
], and the initial value of this steepest ascent algorithm is the maximizer of F
σ
[
s
(
t
)
] obtained for the previous (larger) value of σ. The CS reconstruction algorithm can be described in Algorithm 1. In Algorithm 1, μ is a small positive constant, P
+
=
P
H
(
P
P
H
)
-
1 denotes the Moore-Penrose pseudoinverse of the matrix P, and x
∘
y denotes the Hadamard product (entrywise multiplication) of the vectors x and y.Algorithm 1 (compressed sensing reconstruction algorithm).
Step 1. Initialization.
Step 1.1. Let s
0
(
t
) be equal to the minimum l
2 norm solution of P
s
(
t
)
=
y
(
t
), obtained by s
0
(
t
)
=
P
+
y
(
t
).
Step 1.2. Choose a suitable decreasing sequence [
σ
1
,
…
,
σ
J
].
Step 2. For each σ
j, j
=
1
,
…
,
J, maximize (approximately) the function F
σ
j
[
s
(
t
)
] on the feasible set y
t
-
P
s
t
2
2
≤
ɛ.
Step 2.1. Initialize s
temp
=
s
j
-
1
(
t
).
Step 2.2. Calculate the gradient ∇
F
σ
j
[
s
(
t
)
]
=
(
1
/
σ
j
2
)
s
temp
∘
exp
⁡
{
-
(
s
temp
*
∘
s
temp
)
/
2
σ
j
2
}.
Step 2.3. Carry out the steepest ascent algorithm s
temp
=
s
temp
-
μ
σ
j
2
∇
F
σ
j
[
s
(
t
)
].
Step 2.4. Project onto the feasible set s
temp
=
s
temp
-
P
+
[
P
s
temp
-
y
(
t
)
].
Step 2.5. Assign s
j
(
t
)
=
s
temp.
Step 2.6. If y
t
-
P
s
temp
2
2
>
ɛ and is less than the maximum iteration time, then return to Step 2.2.
Step 3. Obtain the final solution s
^
(
t
)
=
s
J
(
t
).As stated in [31, 32], some remarks are presented here. (1) The internal loop for a fixed σ is repeated a fixed and small number of times. In other words, for increasing the speed, we do not wait for the steepest ascent algorithm to converge. We just need to enter the region near the (global) maximizer of F
σ
[
s
(
t
)
] for escaping from its local maximizers. (2) In Step 1.1 of Algorithm 1, we use the minimum l
2 norm solution (which corresponds to σ
→
∞) as initial estimation of the sparse solution. For the value of σ
1, it may be chosen about two to four times of the maximum absolute value of the elements in s
0
(
t
). Then we use σ
j
=
c
σ
j
-
1, j
≥
2, to determine the next values of σ, where c is usually chosen between 0.5 and 1. For the smallest σ (i.e., σ
J), it can be set about one to two times of (a rough estimation) the standard deviation of the noise. (3) For the selection of μ, the line search can be applied to find the locally optimal μ to minimize F
σ
j
[
s
temp
-
μ
σ
j
2
∇
F
σ
j
(
s
temp
)
]. Performing the line search can be time-consuming, so we just use a fixed small μ in this paper (e.g., μ
=
0.1). It is worth mentioning that we use iteration step-size μ
σ
j
2 in the steepest ascent algorithm, which is proportional to the σ
j
2. The reason is that for smaller values of σ, the function F
σ
[
s
(
t
)
] is more “fluctuating” and hence smaller step-sizes should be used for its maximization.At the end of this section, we take a brief complexity analysis for Algorithm1. The cost of gradient calculation (Step 2.2) and steepest ascent algorithm (Step 2.3) is about O
(
Γ
C
exp
⁡
+
3
Γ
), where C
exp
⁡ is the exponential complexity. For the projection (Step 2.4), it has complexity of O
(
2
M
Γ
+
M
3
) which arises from the matrix inverse and matrix-vector multiplications. So the complexity of Algorithm 1 is about O
{
[
max
⁡
(
C
exp
⁡
,
2
M
)
Γ
+
M
3
]
Q
J
}, where J and Q are the numbers of outer and inner iterations.


### 3.3. Iterative LCMV Beamforming Algorithm

After obtainings
^
(
t
) and reconstructing the echo signals using x
^
(
t
)
=
H
s
^
(
t
), we should apply adaptive beamforming algorithms to enhance the desired signal and suppress interferences and noise with the recovered data. In this paper, we use the iterative LCMV beamforming algorithm to form antenna beam. Of course, other beamforming methods can also be applied, such as subspace projection algorithm, and generalized side-lobe canceller [1].The LCMV algorithm minimizes the total variance of the beamformer output and gets the weight by solving the linear constraint equation [34](18)
w
^
t
=
argmin
w
w
H
R
t
w
s
.
t
.
w
H
a
(
θ
s
)
=
1
,where R
(
t
)
=
E
[
x
^
(
t
)
x
^
H
(
t
)
] is the covariance matrix of reconstructed signal x
^
(
t
) and a
(
θ
s
) is the steering vector in the desired direction. Using the method of Lagrange multiplier, the optimized weight vector can be achieved by(19)
w
^
t
=
R
-
1
t
a
θ
s
a
H
θ
s
R
-
1
t
a
θ
s
-
1
.For online implementation, we can also develop an iterative formula for beamforming weight calculation [35]:(20)
w
^
t
=
A
w
^
t
-
1
-
α
x
^
t
x
^
H
t
w
^
t
-
1
+
F
,where α is the iteration step size, F
=
a
(
θ
s
)
[
a
H
(
θ
s
)
a
(
θ
s
)
]
-
1, and A
=
I
-
F
a
H
(
θ
s
).For the stability of the algorithm,α should satisfy(21)
0
<
α
<
2
λ
max
⁡
,where λ
max
⁡ is the maximum eigenvalue of the covariance matrix R
(
t
)
=
E
[
x
^
(
t
)
x
^
H
(
t
)
]. Since R
(
t
) is positive semidefinite, we have(22)
λ
max
⁡
≤
tr
⁡
R
t
,where tr
⁡
{
·
} is the trace of the matrix argument. Then we have(23)
α
<
2
tr
⁡
R
tas an upper limit for α. When exceeding this limit, α is likely to cause the iterative LCMV algorithm unstability.After obtaining the beamforming weight, we can get the output of adaptive beamformer using(24)
z
t
=
w
^
H
t
x
^
t
.


## 3.1. Smoothedl
0 Norm [31, 32]

It is known thatl
0 norm of a vector is a discontinuous function; searching the minimum l
0 norm solution is a combinatorial problem, which is nonconvex and highly nonsmoothed. Furthermore, l
0 norm of a vector is extremely sensitive to noise, and its value will be completely changed even though the noise is very weak. In this section, we introduce a family of smooth complex approximators of l
0 norm, whose minimization can be implemented using gradient based methods. Consider a complex zero-mean Gaussian family of functions with the parameter σ
(12)
f
σ
c
=
exp
⁡
-
c
2
2
σ
2
=
exp
⁡
-
c
r
2
+
c
i
2
2
σ
2
,where |
c
|, c
r, and c
i represent the module and real and imaginary parts of c
∈
C, we have(13)
lim
σ
→
0
f
σ
(
c
)
=
1
,
c
=
0
0
,
c
≠
0or approximately(14)
f
σ
(
c
)
≈
1
,
c
≪
σ
0
,
c
≫
σ
.Then, for a complex vector s
=
(
s
1
,
s
2
,
…
,
s
Γ
)
T, we can define(15)
F
σ
s
=
∑
i
=
1
Γ
‍
f
σ
s
i
.Since the number of entries in s
(
t
) is Γ and the function F
σ
[
s
(
t
)
] is an indicator of the number of zero-entries in s
(
t
), the l
0 norm of s
(
t
) can be approximated by(16)
s
t
0
≈
Γ
-
F
σ
s
tfor small values of σ, and the approximation tends to equality when σ
→
0.Substituting this approximation into (11) yields the CS problem with smoothed l
0 norm(17)
s
^
(
t
)
=
argmin
s
t
Γ
-
F
σ
s
t
=
argmax
s
t
F
σ
s
t
s
.
t
.
y
t
-
P
s
t
2
2
≤
ɛ
.


## 3.2. Compressed Sensing Reconstruction Algorithm

For small values ofσ in (17), F
σ
[
s
(
t
)
] contains a lot of local maxima. Consequently, it is very difficult to directly maximize this function for very small values of σ. However, as the value of σ grows, the function becomes smoother and smoother, and, for sufficiently large values of σ, there is no local maxima. Similar to [31, 32], our approach to solve the problem (17) is then to decrease the value of σ gradually. The underlying thought is to select a σ which ensures that the initial optimization problem is convex and gradually increase the accuracy of the approximation. By careful selection of the sequence of σ, nonconvexity and thereby local minima can be effectively avoided. For each σ, we use some iterations of the steepest ascent algorithm (followed by projection onto the feasible set) to maximize F
σ
[
s
(
t
)
], and the initial value of this steepest ascent algorithm is the maximizer of F
σ
[
s
(
t
)
] obtained for the previous (larger) value of σ. The CS reconstruction algorithm can be described in Algorithm 1. In Algorithm 1, μ is a small positive constant, P
+
=
P
H
(
P
P
H
)
-
1 denotes the Moore-Penrose pseudoinverse of the matrix P, and x
∘
y denotes the Hadamard product (entrywise multiplication) of the vectors x and y.Algorithm 1 (compressed sensing reconstruction algorithm).
Step 1. Initialization.
Step 1.1. Let s
0
(
t
) be equal to the minimum l
2 norm solution of P
s
(
t
)
=
y
(
t
), obtained by s
0
(
t
)
=
P
+
y
(
t
).
Step 1.2. Choose a suitable decreasing sequence [
σ
1
,
…
,
σ
J
].
Step 2. For each σ
j, j
=
1
,
…
,
J, maximize (approximately) the function F
σ
j
[
s
(
t
)
] on the feasible set y
t
-
P
s
t
2
2
≤
ɛ.
Step 2.1. Initialize s
temp
=
s
j
-
1
(
t
).
Step 2.2. Calculate the gradient ∇
F
σ
j
[
s
(
t
)
]
=
(
1
/
σ
j
2
)
s
temp
∘
exp
⁡
{
-
(
s
temp
*
∘
s
temp
)
/
2
σ
j
2
}.
Step 2.3. Carry out the steepest ascent algorithm s
temp
=
s
temp
-
μ
σ
j
2
∇
F
σ
j
[
s
(
t
)
].
Step 2.4. Project onto the feasible set s
temp
=
s
temp
-
P
+
[
P
s
temp
-
y
(
t
)
].
Step 2.5. Assign s
j
(
t
)
=
s
temp.
Step 2.6. If y
t
-
P
s
temp
2
2
>
ɛ and is less than the maximum iteration time, then return to Step 2.2.
Step 3. Obtain the final solution s
^
(
t
)
=
s
J
(
t
).As stated in [31, 32], some remarks are presented here. (1) The internal loop for a fixed σ is repeated a fixed and small number of times. In other words, for increasing the speed, we do not wait for the steepest ascent algorithm to converge. We just need to enter the region near the (global) maximizer of F
σ
[
s
(
t
)
] for escaping from its local maximizers. (2) In Step 1.1 of Algorithm 1, we use the minimum l
2 norm solution (which corresponds to σ
→
∞) as initial estimation of the sparse solution. For the value of σ
1, it may be chosen about two to four times of the maximum absolute value of the elements in s
0
(
t
). Then we use σ
j
=
c
σ
j
-
1, j
≥
2, to determine the next values of σ, where c is usually chosen between 0.5 and 1. For the smallest σ (i.e., σ
J), it can be set about one to two times of (a rough estimation) the standard deviation of the noise. (3) For the selection of μ, the line search can be applied to find the locally optimal μ to minimize F
σ
j
[
s
temp
-
μ
σ
j
2
∇
F
σ
j
(
s
temp
)
]. Performing the line search can be time-consuming, so we just use a fixed small μ in this paper (e.g., μ
=
0.1). It is worth mentioning that we use iteration step-size μ
σ
j
2 in the steepest ascent algorithm, which is proportional to the σ
j
2. The reason is that for smaller values of σ, the function F
σ
[
s
(
t
)
] is more “fluctuating” and hence smaller step-sizes should be used for its maximization.At the end of this section, we take a brief complexity analysis for Algorithm1. The cost of gradient calculation (Step 2.2) and steepest ascent algorithm (Step 2.3) is about O
(
Γ
C
exp
⁡
+
3
Γ
), where C
exp
⁡ is the exponential complexity. For the projection (Step 2.4), it has complexity of O
(
2
M
Γ
+
M
3
) which arises from the matrix inverse and matrix-vector multiplications. So the complexity of Algorithm 1 is about O
{
[
max
⁡
(
C
exp
⁡
,
2
M
)
Γ
+
M
3
]
Q
J
}, where J and Q are the numbers of outer and inner iterations.


## 3.3. Iterative LCMV Beamforming Algorithm

After obtainings
^
(
t
) and reconstructing the echo signals using x
^
(
t
)
=
H
s
^
(
t
), we should apply adaptive beamforming algorithms to enhance the desired signal and suppress interferences and noise with the recovered data. In this paper, we use the iterative LCMV beamforming algorithm to form antenna beam. Of course, other beamforming methods can also be applied, such as subspace projection algorithm, and generalized side-lobe canceller [1].The LCMV algorithm minimizes the total variance of the beamformer output and gets the weight by solving the linear constraint equation [34](18)
w
^
t
=
argmin
w
w
H
R
t
w
s
.
t
.
w
H
a
(
θ
s
)
=
1
,where R
(
t
)
=
E
[
x
^
(
t
)
x
^
H
(
t
)
] is the covariance matrix of reconstructed signal x
^
(
t
) and a
(
θ
s
) is the steering vector in the desired direction. Using the method of Lagrange multiplier, the optimized weight vector can be achieved by(19)
w
^
t
=
R
-
1
t
a
θ
s
a
H
θ
s
R
-
1
t
a
θ
s
-
1
.For online implementation, we can also develop an iterative formula for beamforming weight calculation [35]:(20)
w
^
t
=
A
w
^
t
-
1
-
α
x
^
t
x
^
H
t
w
^
t
-
1
+
F
,where α is the iteration step size, F
=
a
(
θ
s
)
[
a
H
(
θ
s
)
a
(
θ
s
)
]
-
1, and A
=
I
-
F
a
H
(
θ
s
).For the stability of the algorithm,α should satisfy(21)
0
<
α
<
2
λ
max
⁡
,where λ
max
⁡ is the maximum eigenvalue of the covariance matrix R
(
t
)
=
E
[
x
^
(
t
)
x
^
H
(
t
)
]. Since R
(
t
) is positive semidefinite, we have(22)
λ
max
⁡
≤
tr
⁡
R
t
,where tr
⁡
{
·
} is the trace of the matrix argument. Then we have(23)
α
<
2
tr
⁡
R
tas an upper limit for α. When exceeding this limit, α is likely to cause the iterative LCMV algorithm unstability.After obtaining the beamforming weight, we can get the output of adaptive beamformer using(24)
z
t
=
w
^
H
t
x
^
t
.


## 4. DOAs Not on the Partition Grid

Thus far, in our framework, the estimates of the signal locations are confined to be on the grid of partition. When the signals are not on the grids, the error of the reconstructed signal grows heavily. Increasing the parameterΓ to reduce the grid spacing can solve this problem. But finer uniform grid will increase the computational complexity significantly. Instead of having a uniform fine grid, we make the grid fine only around the regions where signals are present [36]. This requires an approximate knowledge of the locations of the sources, which can be obtained by using a coarse grid first. The steps to form nonuniform refined grid are as follows.(1) Create a coarse partition grid of potential signal locationsθ
~
(
0
)
=
[
-
9
0
°
=
θ
~
1
(
0
)
,
θ
~
2
(
0
)
,
…
,
θ
~
Γ
(
0
)
=
9
0
°
] such that |

sin
(
θ
~
i
+
1
(
0
)
)
-
sin
(
θ
~
i
(
0
)
)
|
=
2
/
(
Γ
-
1
). The grid should not be too rough in order to get the approximated locations estimation of the sources. Here, for example, we initialize Γ
=
200.(2) Form the transformation matrixH related to θ
~
(
0
) and use the CS reconstruction algorithm to reconstruct the sparse signal s
~
(
0
)
(
t
). Find the indices of K largest amplitude of nonzero elements in s
~
(
0
)
(
t
) and get the coarse locations estimation θ
~
j
(
0
), j
=
1
,
…
,
K.(3) As illustrated in Figure1, for each θ
~
j
(
0
), j
=
1
,
…
,
K, divide the interval [
θ
~
j
-
1
(
0
)
,
θ
~
j
+
1
(
0
)
] into 2
γ parts in terms of angular sine representation, where γ is the refined times. Then we can get the refined partition grid around the locations θ
~
j
(
0
), j
=
1
,
…
,
K, and denote the new nonuniform partition grid as θ
~
(
1
).Figure 1
Illustration of grid refinement.


## 5. Numerical Simulation

In this section, we will perform some simulations for linear and planar arrays to demonstrate the effectiveness of our proposed adaptive beamforming.


### 5.1. The Evaluation of DOA and Signal Recovery Error

A uniform linear array withN
=
100 omnidirectional sensors spaced half a wavelength apart is considered. The noise is white complex Gaussian with power σ
n
2
=
1 in the simulations. We select M
=
30 elements from the N array elements randomly and set Γ
=
200 for the angle partition. It is assumed that there is only one incident signal from the directions of sin
(
θ
)
=
0.1 with signal-to-noise ratio (SNR) 10 dB. For comparisons, we compare the OMP, iterative adaptive approach for amplitude and phase estimation (IAA-APES) [37], and propose CS-based DOA methods (called CS-OMP, CS-IAA-APES, and CS-SL0). Figure 2 shows the DOAs of three methods with one incident signal, and Figure 3 gives the signal recovery errors of three different methods. Here we define the signal recovery error as(25)
Error
=
x
^
-
x
2
x
2
,where x
^ is the reconstructed signal and x is the original signal without noise. From Figures 2 and 3, we can see that these three methods have similar DOA performance. However the proposed CS-SL0 outperforms CS-OMP and CS-IAA-APES in terms of signal recovery error. The reason is that the proposed CS-SL0 obtains more accurate complex amplitude vector s
^
(
t
) than the other two CS-based DOA estimations.Figure 2
The DOAs of three different methods with one incident signal (SNR = 10 dB).Figure 3
The signal recovery errors of three different methods with one incident signal.


### 5.2. The Evaluation of Beam Pattern and Output SINR

We also consider a uniform linear array withN
=
100 omnidirectional sensors spaced half a wavelength apart. Assume that there are four space signals from four different directions. One is the desired signal with SNR 10 dB; the other three are interferences with interference-to-noise ratio (INR) 40 dB. The directions of the desired signal and three interferences are set to be 0.1, −0.21, 0.41, and −0.45 in terms of angular sine representation (i.e., sin
(
θ
)), respectively. The other parameters are the same as those in Section 5.1. Three kinds of the beam patterns are shown in Figure 4 by applying iterative LCMV beamforming to the CS-SL0 reconstructed signal, received signal with 30 elements and 100 elements (i.e., full array), respectively. From Figure 4, we can see that the antenna aperture is not reduced when the number of the elements is decreased from 100 to 30. The side-lobe of the beam formed by the signal of 30 elements is −9 dB. The beam formed by the proposed CS-SL0 method has similar performance to the beam formed by the full array (100 elements) signal.Figure 4
The beam patterns of iterative LCMV with different signals.In order to verify the effectiveness of the algorithm in the cases of different SNRs and different DOAs, we randomly select the directions of the desired signal and interferences. The directions of the interferences are assumed to be out of the main lobe. We then do 100 Monte Carlo trials and use a Taylor weight to suppress the level of side-lobes to be lower than −20 dB for practical application. Figure5 shows the average output SINR versus SNR when the input SNR is varied from −10 dB to 30 dB and signal-to-interference ratio (SIR) is set to be −30 dB. From Figure 5, it can be seen that, when the SNR is lower than 0 dB, the desired signal may be lost by using the CS reconstruction (Algorithm 1). It shows that the recovered signal of full array is degraded and the target cannot be detected. When the SNR is higher than 0 dB, the desired signal can be recovered reliably.Figure 5
Output SINR versus SNR of iterative LCMV with CS-SL0 reconstructed signal.To evaluate the performance of side-lobe level and interference suppression, Monte Carlo simulations about side-lobe level and null depth are provided in Figure6 where the input SNR is also from −10 dB to 30 dB and SIR = −30 dB. Obviously, the levels of the side-lobe are lower than −20 dB and almost the same with the change of SNR. The null depths are all lower than −45 dB and decreased with the increase of input SNR. As a result, good performance of the beam pattern can be obtained by using the iterative LCMV algorithm with CS reconstructed signal.Figure 6
The levels of side-lobe and null depths in interference directions.


### 5.3. Simulation of DOAs Not on the Grid

In this subsection, we consider the situation of the DOAs not on the grid. Assume that there are two spatial signals from two different directions. One is the desired signal; the other is the interference. The DOAs of the desired signal and the interference are randomly generated (not on the grid). Set SNR = 0~30 dB, SIR = −30 dB, andΓ
=
200 for the initial angle partition. Then we do Monte Carlo analysis 100 times.For comparisons, we also provide the results of adaptive beamforming using CS-OMP reconstruction method [17]. Figure 7 shows the output SINRs of the coarse and fine grids for CS-SL0 and CS-OMP methods. From Figure 7, we can see that, when the desired signal or the interferences are not on the grid, the performance of the proposed algorithm will degrade seriously. Using the refined method stated in Section 4, we can get two fine nonuniform grids with γ
=
2 and γ
=
5, respectively. For the fine nonuniform grids, the errors of the reconstructed signals have been reduced significantly, which results in higher output SINRs compared to the coarse grid. Also we can see that the performance of CS-SL0 is better than CS-OMP method. For example, when γ
=
5 and SNR = 20 dB, the output SINRs of CS-OMP and CS-SL0 methods are 29.1 dB and 32.2 dB, and there is 3.1 dB improvement. Figure 8 shows the beam pattern with γ
=
5 from one trial, where the directions of the desired signal and the interference are set to be 0.105 and 0.455 (sin
(
θ
)). From Figure 8, we can see that the main lobe is steered to the desired signal direction and the null is deep in the direction of the interference. Therefore, when DOAs are not on the grid, the proposed adaptive beamforming still works with the refined method in Section 4.Figure 7
The output SINRs of the coarse and fine grids for CS-SL0 and CS-OMP.Figure 8
The beam pattern of DOAs not on the grid.


### 5.4. Simulation of 2D Planar Array

Finally, to further illustrate the performance of the proposed adaptive beamformer for 2D planar array, we also test our algorithm for a planar circle array spaced half a wavelength apart. Assume that the radius of the circle array is6
λ. So the number of the arrays is N
=
1051. The entire UV-plane is partitioned into 39
×
39 parts. We randomly select 100 elements from the N array elements. Figure 9 shows the distribution of the 100 elements. The noise power is σ
n
2
=
1, and a desired source signal is from the direction of (
u
,
v
)
=
(
0.1,0
) with SNR = 10 dB. Two spatially fixed interferers arrive from (
u
,
v
)
=
(
0.1,0.5
) and (
u
,
v
)
=
(
0.1
,
-
0.4
) with INR = 40 dB.Figure 9
The distribution of the sparse array element.The experimental results have been shown in Figures10 and 11. Figure 10 is the 3D beam pattern, and Figure 11 is the 2D U-cut beam pattern with u
=
0.1. From Figures 10 and 11, it is also shown that the antenna aperture has not been reduced when the number of the elements is decreased from 1051 to 100. The side-lobe of the beam formed by the recovery data is −40 dB, and the null depths in the direction of the interferences are deep. Therefore, the proposed adaptive beamforming still works for planar circle array.Figure 10
The beam pattern of planar circle array in surface plot.Figure 11
The beam pattern of planar circle array in U-cut plot withu
=
0.1.


## 5.1. The Evaluation of DOA and Signal Recovery Error

A uniform linear array withN
=
100 omnidirectional sensors spaced half a wavelength apart is considered. The noise is white complex Gaussian with power σ
n
2
=
1 in the simulations. We select M
=
30 elements from the N array elements randomly and set Γ
=
200 for the angle partition. It is assumed that there is only one incident signal from the directions of sin
(
θ
)
=
0.1 with signal-to-noise ratio (SNR) 10 dB. For comparisons, we compare the OMP, iterative adaptive approach for amplitude and phase estimation (IAA-APES) [37], and propose CS-based DOA methods (called CS-OMP, CS-IAA-APES, and CS-SL0). Figure 2 shows the DOAs of three methods with one incident signal, and Figure 3 gives the signal recovery errors of three different methods. Here we define the signal recovery error as(25)
Error
=
x
^
-
x
2
x
2
,where x
^ is the reconstructed signal and x is the original signal without noise. From Figures 2 and 3, we can see that these three methods have similar DOA performance. However the proposed CS-SL0 outperforms CS-OMP and CS-IAA-APES in terms of signal recovery error. The reason is that the proposed CS-SL0 obtains more accurate complex amplitude vector s
^
(
t
) than the other two CS-based DOA estimations.Figure 2
The DOAs of three different methods with one incident signal (SNR = 10 dB).Figure 3
The signal recovery errors of three different methods with one incident signal.


## 5.2. The Evaluation of Beam Pattern and Output SINR

We also consider a uniform linear array withN
=
100 omnidirectional sensors spaced half a wavelength apart. Assume that there are four space signals from four different directions. One is the desired signal with SNR 10 dB; the other three are interferences with interference-to-noise ratio (INR) 40 dB. The directions of the desired signal and three interferences are set to be 0.1, −0.21, 0.41, and −0.45 in terms of angular sine representation (i.e., sin
(
θ
)), respectively. The other parameters are the same as those in Section 5.1. Three kinds of the beam patterns are shown in Figure 4 by applying iterative LCMV beamforming to the CS-SL0 reconstructed signal, received signal with 30 elements and 100 elements (i.e., full array), respectively. From Figure 4, we can see that the antenna aperture is not reduced when the number of the elements is decreased from 100 to 30. The side-lobe of the beam formed by the signal of 30 elements is −9 dB. The beam formed by the proposed CS-SL0 method has similar performance to the beam formed by the full array (100 elements) signal.Figure 4
The beam patterns of iterative LCMV with different signals.In order to verify the effectiveness of the algorithm in the cases of different SNRs and different DOAs, we randomly select the directions of the desired signal and interferences. The directions of the interferences are assumed to be out of the main lobe. We then do 100 Monte Carlo trials and use a Taylor weight to suppress the level of side-lobes to be lower than −20 dB for practical application. Figure5 shows the average output SINR versus SNR when the input SNR is varied from −10 dB to 30 dB and signal-to-interference ratio (SIR) is set to be −30 dB. From Figure 5, it can be seen that, when the SNR is lower than 0 dB, the desired signal may be lost by using the CS reconstruction (Algorithm 1). It shows that the recovered signal of full array is degraded and the target cannot be detected. When the SNR is higher than 0 dB, the desired signal can be recovered reliably.Figure 5
Output SINR versus SNR of iterative LCMV with CS-SL0 reconstructed signal.To evaluate the performance of side-lobe level and interference suppression, Monte Carlo simulations about side-lobe level and null depth are provided in Figure6 where the input SNR is also from −10 dB to 30 dB and SIR = −30 dB. Obviously, the levels of the side-lobe are lower than −20 dB and almost the same with the change of SNR. The null depths are all lower than −45 dB and decreased with the increase of input SNR. As a result, good performance of the beam pattern can be obtained by using the iterative LCMV algorithm with CS reconstructed signal.Figure 6
The levels of side-lobe and null depths in interference directions.


## 5.3. Simulation of DOAs Not on the Grid

In this subsection, we consider the situation of the DOAs not on the grid. Assume that there are two spatial signals from two different directions. One is the desired signal; the other is the interference. The DOAs of the desired signal and the interference are randomly generated (not on the grid). Set SNR = 0~30 dB, SIR = −30 dB, andΓ
=
200 for the initial angle partition. Then we do Monte Carlo analysis 100 times.For comparisons, we also provide the results of adaptive beamforming using CS-OMP reconstruction method [17]. Figure 7 shows the output SINRs of the coarse and fine grids for CS-SL0 and CS-OMP methods. From Figure 7, we can see that, when the desired signal or the interferences are not on the grid, the performance of the proposed algorithm will degrade seriously. Using the refined method stated in Section 4, we can get two fine nonuniform grids with γ
=
2 and γ
=
5, respectively. For the fine nonuniform grids, the errors of the reconstructed signals have been reduced significantly, which results in higher output SINRs compared to the coarse grid. Also we can see that the performance of CS-SL0 is better than CS-OMP method. For example, when γ
=
5 and SNR = 20 dB, the output SINRs of CS-OMP and CS-SL0 methods are 29.1 dB and 32.2 dB, and there is 3.1 dB improvement. Figure 8 shows the beam pattern with γ
=
5 from one trial, where the directions of the desired signal and the interference are set to be 0.105 and 0.455 (sin
(
θ
)). From Figure 8, we can see that the main lobe is steered to the desired signal direction and the null is deep in the direction of the interference. Therefore, when DOAs are not on the grid, the proposed adaptive beamforming still works with the refined method in Section 4.Figure 7
The output SINRs of the coarse and fine grids for CS-SL0 and CS-OMP.Figure 8
The beam pattern of DOAs not on the grid.


## 5.4. Simulation of 2D Planar Array

Finally, to further illustrate the performance of the proposed adaptive beamformer for 2D planar array, we also test our algorithm for a planar circle array spaced half a wavelength apart. Assume that the radius of the circle array is6
λ. So the number of the arrays is N
=
1051. The entire UV-plane is partitioned into 39
×
39 parts. We randomly select 100 elements from the N array elements. Figure 9 shows the distribution of the 100 elements. The noise power is σ
n
2
=
1, and a desired source signal is from the direction of (
u
,
v
)
=
(
0.1,0
) with SNR = 10 dB. Two spatially fixed interferers arrive from (
u
,
v
)
=
(
0.1,0.5
) and (
u
,
v
)
=
(
0.1
,
-
0.4
) with INR = 40 dB.Figure 9
The distribution of the sparse array element.The experimental results have been shown in Figures10 and 11. Figure 10 is the 3D beam pattern, and Figure 11 is the 2D U-cut beam pattern with u
=
0.1. From Figures 10 and 11, it is also shown that the antenna aperture has not been reduced when the number of the elements is decreased from 1051 to 100. The side-lobe of the beam formed by the recovery data is −40 dB, and the null depths in the direction of the interferences are deep. Therefore, the proposed adaptive beamforming still works for planar circle array.Figure 10
The beam pattern of planar circle array in surface plot.Figure 11
The beam pattern of planar circle array in U-cut plot withu
=
0.1.


## 6. Conclusions

A new adaptive digital beamforming for large-scale sparse receiving array based on CS with smoothedl
0 norm is presented in this paper. It can greatly reduce the elements of the array and the number of RF front-ends without sacrificing the antenna aperture and deteriorating the performance of the beams. Applying the sparse prior knowledge of the target in space, the whole channel signal of full array can be effectively reconstructed by using the CS reconstruction algorithm with smoothed l
0 norm. Then an adaptive beamformer based on iterative LCMV algorithm is developed with the recovered signal of full array. Simulation results and Monte Carlo analysis for linear and planar arrays demonstrate the effectiveness of the proposed adaptive beamforming, which has similar performance to the beamforming using signal of full array.



---

*Source: 959856-2015-03-18.xml*