Low Complexity Beamspace Super Resolution for DOA Estimation of Linear Array

Abstract

Beamspace processing has become much attractive in recent radar and wireless communication applications, since the advantages of complexity reduction and of performance improvements in array signal processing. In this paper, we concentrate on the beamspace DOA estimation of linear array via atomic norm minimization (ANM). The existed generalized linear spectrum estimation based ANM approaches suffer from the high computational complexity for large scale array, since their complexity depends upon the number of sensors. To deal with this problem, we develop a low dimensional semidefinite programming (SDP) implementation of beamspace atomic norm minimization (BS-ANM) approach for DFT beamspace based on the super resolution theory on the semi-algebraic set. Then, a computational efficient iteration algorithm is proposed based on alternating direction method of multipliers (ADMM) approach. We develop the covariance based DOA estimation methods via BS-ANM and apply the BS-ANM based DOA estimation method to the channel estimation problem for massive MIMO systems. Simulation results demonstrate that the proposed methods exhibit the superior performance compared to the state-of-the-art counterparts.

1. Introduction

Direction-of-Arrive (DOA) estimation is an important topic in many applications, such as radar, sonar, and wireless communication. To lower the hardware cost, reduce the computational burden and improve the performance, the received signals on front-end sensors can be projected into lower dimensional space by digital/analog structure, which is referred to as beamspace processing [1,2,3]. In the past decades, various classic DOA estimation methods have been applied to beamspace [4,5,6,7].

Inspired from the concept of compressed sensing, some discrete sparse representation approaches are extended to the beamspace to improve the DOA estimation performance in low SNR and lack of snapshots’ scenarios [8,9,10]. However, the discrete sparse recovery methods may suffer from the off-grid sources problem, a.k.a. power leakage effect [11].

To deal with the problem, the super resolution methods have recently been proposed for gridless compressed sensing, which is referred to as atomic norm minimization (ANM) [12]. By utilizing the Toeplitz structure of the covariance matrix, ANM approaches are applied to DOA estimation of linear/rectangular arrays [13,14,15,16].

For beamspace processing cases, the existed ANM methods are usually based on the generalized line spectral estimation (GL) framework [17], which regards the beamspace DOA estimation problem as the line spectral estimation problem with linear mapping constraints, which can be solved by the conventional ANM approaches [11,18,19,20]. However, these methods focus on recovering the signal on the sensors of the receiver, which may yield high computational burden (since high dimensional SDP formulation) in the case of large number of sensors, e.g., millimeter-wave massive MIMO system for 5G, even though the dimension of the beamspace may be quite small.

Like the extensive application for gridless super resolution methods with beamspace processing in real-time scenarios, we try to propose the low complexity beamspace super resolution method that is not yet reported in the literature.

We study the beamspace super resolution approaches for discrete Fourier transform (DFT) beamspace via beamspace atomic norm minimization(BS-ANM) in this paper. Firstly, we define the beamspace atomic norm and propose the low dimensional SDP implementations based on the super resolution theory on the semi-algebraic set. Then, we develop the fast algorithm based on the alternating direction method of multipliers (ADMM) method and show the computational complexity analysis. Finally, we apply BS-ANM approaches to the covariance based DOA estimation problem and the channel estimation in millimeter-wave massive MIMO system with lens antenna array. Simulation results indicate that the proposed approaches exhibit the favorable performance and computational cost, as compared to the state-of-the-art methods.

This paper is organized as follows: the signal model and conception of atomic norm are introduced in Section 2. The low dimensional SDP implementation of BS-ANM is proposed in Section 3. In Section 4, we present a fast algorithm based on ADMM for BS-ANM and provide the complexity analysis in Section 5. In Section 6, the BS-ANM based DOA estimation and channel estimation methods are developed and simulations are performed in Section 7. Section 8 concludes this paper.

Notations: ${(•)}^{\mathrm{T}}$ denotes the transpose, ${(•)}^{†}$ denotes pseudo inverse of a matrix, and ${(•)}^{\mathrm{H}}$ denotes conjugate transpose of a matrix or vector. $diag\left(\mathbf{x}\right)$ denotes a diagonal matrix.${∥•∥}_2$ and ${∥•∥}_F$ denote the Euclidean $l_2$ norm and the Frobinus norm. $\mathbb{C},\mathbb{R}$ and $\mathbb{N}$ represent the complex, real and nature number set, respectively. $\mathbf{M}\ge 0$ denotes positive semidifinite matrix. $E\left\{•\right\}$ is the expectation. $deg(•)$ denotes the degree of the polynomial. $vec(•)$ and $\mathrm{Tr}(•)$ are the vectorization and the trace operator. $\inf \left\{•\right\}$ and $\sup \left\{•\right\}$ denote the infimum and supremum of a set, respectively.

2. Signal Model and Atomic Norm

2.1. Signal Model

Considering the uniform linear array (ULA) of N sensors with the $M(M<N)$ dimensional fixed hardware structure beamformer, there are K far field narrowband scatters at angular directions ${\theta }_1,\dots ,{\theta }_K$ and the measurement vector on the sensors at time t can be written as

$$\mathbf{x}\left(t\right)=\mathbf{A}\mathbf{s}\left(t\right)+\mathbf{n}\left(t\right),$$

(1)

where $\mathbf{s}\left(t\right)={[s_1,\cdots ,s_K]}^{\mathrm{T}}\in {\mathbb{C}}^{K\times 1}$ and $\mathbf{n}\left(t\right)\in {\mathbb{C}}^{N\times 1}$ are the source and additional noise component with ${\mathbf{R}}_s=E\left\{\mathbf{s}\left(t\right){\mathbf{s}}^{\mathrm{H}}\left(t\right)\right\}$ and $E\left\{\mathbf{n}\left(t\right){\mathbf{n}}^{\mathrm{H}}\left(t\right)\right\}={\sigma }^2\mathbf{I}$(${\sigma }^2$ is the power of noise). The $N\times K$ array manifold $\mathbf{A}$ is given by

$$\mathbf{A}=[\mathbf{a}\left({\theta }_1\right),\cdots ,\mathbf{a}\left({\theta }_K\right)],$$

(2)

where the column vector of $\mathbf{A}$ is

$$\mathbf{a}{(\theta )=[1,e^{-j\frac{2\pi d}{\lambda }sin\theta },\cdots ,e^{-j\frac{2\pi Nd}{\lambda }sin\theta }]}^{\mathrm{T}}$$

(3)

with ${\theta }_k\in \mathsf{\Xi }$, $\mathsf{\Xi }=(-\pi /2,\pi /2]$, the interspace of elements d, and the wavelength $\lambda$.

Denote the DFT beamspace transform matrix as $\mathbf{W}\in {\mathbb{C}}^{M\times N}$, whose $m\mathrm{th}$ row is [4]

$${\mathbf{w}}_m=\frac{1}{N}{e}^{jm\frac{(N-1)\pi }{2N}}\left[1,e^{-jm\frac{\pi }{N}},\cdots ,e^{-jm\frac{(N-1)\pi }{N}}\right],$$

(4)

the beamspace received signals are modeled as

$$\begin{array}{cc}\hfill \mathbf{y}\left(t\right)& =\mathbf{B}\mathbf{s}\left(t\right)+{\mathbf{n}}_{\mathbf{b}}\left(t\right),\hfill \\ & =\sum _{k=1}^K{s}_k\mathbf{b}\left({\theta }_k\right)+{\mathbf{n}}_{\mathbf{b}}\left(t\right),\hfill \end{array}$$

(5)

where $\mathbf{B}=\mathbf{WA}=[\mathbf{b}\left({\theta }_1\right),\cdots ,\mathbf{b}\left({\theta }_K\right)]$ and ${\mathbf{n}}_{\mathbf{b}}\left(t\right)=\mathbf{Wn}\left(t\right)$ are the beamspace array manifold and noise, respectively. From Equation (5), the covariance matrix in beamspace is

$$\begin{array}{cc}\hfill \mathbf{R}& =E\left\{\mathbf{Y}\left(t\right){\mathbf{Y}}^{\mathrm{H}}\left(t\right)\right\}\hfill \\ & =\mathbf{B}{\mathbf{R}}_s{\mathbf{B}}^{\mathrm{H}}+\mathbf{V}\hfill \end{array}$$

(6)

with $\mathbf{V}={\sigma }^2\mathbf{W}{\mathbf{W}}^{\mathrm{H}}$. Inspired from [21], the covariance matrix in (5) can be represented in the Multiple Measurement Vectors (MMV) model as

$$\begin{array}{cc}\hfill \mathbf{R}& =\mathbf{B}\mathbf{H}+\mathbf{V}\hfill \\ & =\sum _{k=1}^K{c}_k\mathbf{b}\left({\theta }_k\right){\mathbf{h}}_k+\mathbf{V}\hfill \end{array}$$

(7)

where $\mathbf{H}={\mathbf{R}}_s{\mathbf{B}}^{\mathrm{H}}$, $c_k{\mathbf{h}}_k$ is the $k\mathrm{th}$ row vector of $\mathbf{H}$ and ${∥{\mathbf{h}}_k∥}_2=1$.

2.2. Atomic Norm in Beamspace

As shown in [12], the element-space signals of the ULA in Equation (1) are consisted of distinct atoms over the atomic set

$$\mathcal{A}=\left\{{\left[\mathbf{a}\left({\theta }_k\right)\right]}_n=e^{-j\frac{2\pi ndsin\theta }{\lambda }}|{\theta }_k\in \mathsf{\Xi }\right\},$$

(8)

which can be represented by the atomic decomposition as

$$\mathbf{x}=\sum _{k=1}^K{s}_k\mathbf{a}\left({\theta }_k\right),\mathbf{a}\left(\theta \right)\in \mathcal{A}.$$

(9)

The atomic decomposition with the smallest total variation over $\mathcal{A}$ yields the ${\ell }_1$ atomic norm as

$$\begin{array}{c}\hfill {∥\mathbf{x}∥}_{\mathcal{A}}=inf\left\{\sum _{k=1}^K|s_k||\mathbf{x}=\sum _{k=1}^K{s}_k\mathbf{a}\left({\theta }_k\right),\mathbf{a}\left(\theta \right)\in \mathcal{A}\right\}.\end{array}$$

(10)

Followed by generalized line spectral estimation framework in [17], the beamspace signal in Equation (5) can be retrieved by the following atomic norm minimization (ANM) problem:

$$\begin{array}{c}\underset{\mathbf{x}}{min}{∥\mathbf{x}∥}_{\mathcal{A}}\hfill \\ \mathrm{s}.\mathrm{t}.\mathbf{y}=\mathbf{Wx},\hfill \end{array}$$

(11)

if the locations of sources are separated sufficiently. Followed by the Vandermonde decomposition theorem, the SDP implementation of Equation (11) can written as [12]

$$\begin{array}{c}{∥\mathbf{x}∥}_{\mathcal{T}}=\underset{\mathbf{Q},w,\mathbf{z}}{min}\frac{1}{2N}\mathrm{Tr}\left(S\left(\mathbf{Q}\right)\right)+\frac{1}{2}w\hfill \\ \mathrm{s}.\mathrm{t}.\left[\begin{array}{cc}S\left(\mathbf{Q}\right)& \mathbf{x}\\ {\mathbf{x}}^{\mathrm{H}}& w\end{array}\right]\ge 0,\hfill \\ \mathbf{y}=\mathbf{Wx}\hfill \end{array}$$

(12)

where ${∥\mathbf{x}∥}_{\mathcal{T}}={∥\mathbf{x}∥}_{\mathcal{A}}$ and $S\left(\mathbf{Q}\right)$ is a Toeplitz matrix defined by $\mathbf{Q}$.

It is noted that the method in Equation (11) can be applied to arbitrary beamspace design; however, it may result in high computational burden when the number of sensors is large, due to the fact that the computational complexity of SDP based implementation to Equation (12) depends on the dimension of $\mathbf{x}$ [22].

To develop the low complexity super resolution approach for the signal in beamspace, in this paper, we focus on the beamspace atomic norm

$$\begin{array}{c}\hfill {∥\mathbf{y}∥}_{\mathcal{B}}=inf\left\{\sum _{k=1}^K|{s^{\prime }}_k||\mathbf{y}=\sum _{k=1}^K{s^{\prime }}_k{\mathbf{b}}^{\prime }\left({\theta }_k\right),{\mathbf{b}}^{\prime }\left(\theta \right)\in \mathcal{B}\right\}\end{array}$$

(13)

directly, with

$$\mathcal{B}=\left\{{\mathbf{b}}^{\prime }\left(\theta \right)|{\mathbf{b}}^{\prime }\left(\theta \right)=\mathbf{Wa}\left(\theta \right)/{∥\mathbf{Wa}\left(\theta \right)∥}_F,{\theta }_k\in \mathsf{\Xi }\right\}.$$

(14)

We will introduce an approximate low dimensional SDP implementation of the beamspace atomic norm minimization (BS-ANM) problem in Equation (13) for DFT beamspace design in the next section.

3. Low Dimensional Sdp Implementation for Bs-Anm

The generalized line spectral estimation based methods rely on the Vandermonde decomposition theorem. However, the beamspace array manifolds are not Vandermonde, which implies the major difficulty of realizing that the BS-ANM approach is to formulate the solvable convex optimization problem without the Vandermonde decomposition theorem. In this section, we propose an approximate low dimensional SDP implementation to the BS-ANM problem based on the super resolution theory on the semi-algebraic set [23].

We first introduce the dual problem of Equation (13) with Lagrangian analysis,

$$\begin{array}{c}\underset{\mathsf{\Omega }}{max}Re\left\{\mathrm{Tr}\left[{\mathbf{y}}^{\mathrm{H}}\mathsf{\Omega }\right]\right\}\hfill \\ \mathrm{s}.\mathrm{t}.{∥\mathsf{\Omega }∥}_{\mathcal{B}}^{*}\le 1,\hfill \end{array}$$

(15)

which can retrieve the atomic decomposition of the ground truth if there exists the polynomial $\mathbf{q}\left(\theta \right)={\mathbf{b}}^{\prime }\mathrm{H}\left(\theta \right)\mathsf{\Omega }$ satisfying the following dual certificate:

$$\left\{\begin{array}{c}\mathbf{q}\left({\theta }_k\right)=1,\forall k=1,2,\dots ,K\hfill \\ {∥\mathbf{q}\left(\theta \right)∥}_2<1,\forall \theta \ne {\theta }_k\hfill \end{array}\right.$$

(16)

where ${∥•∥}_{\mathcal{B}}^{*}$ is the dual norm corresponding to ${∥•∥}_{\mathcal{B}}$ that

$$\begin{array}{cc}\hfill {∥\mathsf{\Omega }∥}_{\mathcal{B}}^{*}& =\underset{{∥\mathbf{G}∥}_{\mathcal{B}}\le 1}{sup}Re\left\{\mathrm{Tr}\left[{\mathbf{G}}^{\mathrm{H}}\mathsf{\Omega }\right]\right\}\hfill \\ & =\underset{\theta \in \mathsf{\Xi },{∥\rho ∥}_2=1}{\sup }Re\left\{\mathrm{Tr}\left[{\mathsf{\Omega }}^{\mathrm{H}}{\mathbf{b}}^{\prime }\left(\theta \right)\rho \right]\right\}\hfill \\ & =\underset{\theta \in \mathsf{\Xi },{∥\rho ∥}_2=1}{\sup }\left|\mathrm{Tr}\left[\rho {\mathsf{\Omega }}^{\mathrm{H}}{\mathbf{b}}^{\prime }\left(\theta \right)\right]\right|\hfill \\ & =\underset{\theta \in \mathsf{\Xi }}{\sup }{∥{\mathsf{\Omega }}^{\mathrm{H}}{\mathbf{b}}^{\prime }\left(\theta \right)∥}_2.\hfill \end{array}$$

(17)

Letting $z=e^{-j\frac{2\pi ndsin\theta }{\lambda }}$, we have

$$\begin{array}{cc}\hfill & {∥\mathsf{\Omega }∥}_{\mathcal{B}}^{*}\le 1\hfill \\ \hfill \iff & F\left(z\right)=1-{{\mathbf{b}}^{\prime }}^{\mathrm{H}}\left(\theta \right){\mathsf{\Omega }}^{\mathrm{H}}\mathsf{\Omega }{\mathbf{b}}^{\prime }\left(\theta \right)\ge 0,\hfill \end{array}$$

(18)

i.e., $F\left(z\right)$ is a nonnegative trigonometric polynomial. In addition, $F\left(z\right)$ is also a polynomial of $\mathbf{b}$, which can be denoted as $f\left(\mathbf{b}\right)$, where $\mathbf{b}$ denotes $\mathbf{b}\left(\theta \right)$ in Equation (5).

To insight the property of the polynomial $f\left(\mathbf{b}\right)$, we introduce the following theorems:

Theorem 1.

Given a bounded semi-algebraic set

$$\mathcal{D}\left(g\right)=\left\{\left.\mathbf{b}\in {\mathbb{R}}^M\right|g_v\left(\mathbf{b}\right)\ge 0,v=1,2,\cdots ,V\right\},$$

(19)

if there exists the nonnegative polynomial $f\left(\mathbf{b}\right)\ge 0$ for any $\mathbf{b}\in \mathcal{D}\left(g\right)$, we have

$$\begin{array}{c}\hfill f\left(\mathbf{b}\right)=\sum _{\alpha \in {\left\{0,1\right\}}^V}{g}_1^{{\alpha }_1}\cdots g_V^{{\alpha }_V}{s}_{\alpha },\end{array}$$

(20)

with $g_v\left(\mathbf{b}\right)\in \mathbb{R}\left[\mathbf{b}\right]$ and $s_{\alpha }\in \sum \mathbb{R}{\left[\mathbf{b}\right]}^2$, where $\mathbb{R}\left[\mathbf{b}\right]$ and $\sum \mathbb{R}{\left[\mathbf{b}\right]}^2$ denote the set of real polynomials and real sum of squares polynomials on $\mathbf{b}$, respectively.

Proof of Theorem 1.

Theorem 1 can be found in [24] as Theorem 4.1 (See also Corollary 3 in [25]), thus we omit the proof here.  □

Theorem 2.

The beamspace array manifold $\mathcal{B}$ in Equation (14) with the DFT beamspace in Equation (4) is a semi-algebraic set as $\mathcal{D}\left(g\right)$ in Equation (19).

Proof of Theorem 2.

Followed by the invariant relationship between $b_m$ and $b_{m+1}$ in [4] and the property of trigonometric functions, where $b_m$ denotes the $m\mathrm{th}$ element of $\mathbf{b}\left(\theta \right)$, we have

$$\left\{\begin{array}{c}\begin{array}{c}{v}_m(b_m,b_{m+1})sin\left(\frac{u}{2}\right)=p_m(b_m,b_{m+1})cos\left(\frac{u}{2}\right),\hfill \\ {sin}^2\left(\frac{u}{2}\right)+{cos}^2\left(\frac{u}{2}\right)=1,\hfill \\ sin\left(u\right)=2sin\left(\frac{u}{2}\right)cos\left(\frac{u}{2}\right),\hfill \\ cos\left(u\right)=2{cos}^2\left(\frac{u}{2}\right)-1\hfill \end{array}\end{array}\right.$$

(21)

with $u=\frac{2\pi nd}{\lambda }sin\theta$ and

$$\begin{array}{cc}\hfill v_m(b_m,b_{m+1})& =cos\left(\frac{m\pi }{N}\right)b_m+cos\left(\frac{(m+1)\pi }{N}\right)b_{m+1},\hfill \\ \hfill p_m(b_m,b_{m+1})& =sin\left(\frac{m\pi }{N}\right)b_m+sin\left(\frac{(m+1)\pi }{N}\right)b_{m+1}.\hfill \end{array}$$

(22)

Solving the polynomial system in Equation (21) implies that $cos\left(u\right)$ and $sin\left(u\right)$ can be determined by the polynomial of $\left\{b_m,b_{m+1}\right\}$ uniquely, i.e., $cos\left(u\right)=L_m(b_m,b_{m+1})$ and $sin\left(u\right)=J_m(b_m,b_{m+1})$. Accordingly, we can represent $b_m$ as a polynomial of $\mathbf{b}$ such that $b_m=H_m\left(\mathbf{b}\right)$.

Combining the constraints mentioned above together

$$\left\{\begin{array}{c}{L}_1(b_1,b_2)=L_2(b_2,b_3)=\cdots =L_{M-1}(b_{M-1},b_M),\hfill \\ J_1(b_1,b_2)=J_2(b_2,b_3)=\cdots =J_M(b_{M-1},b_M),\hfill \\ L_m^2(b_m,b_{m+1})+J_m^2(b_m,b_{m+1})=1,\hfill \\ H_m\left(\mathbf{b}\right)=b_m,m=1,\cdots ,M,\hfill \end{array}\right.$$

(23)

it implies the solution of the polynomials system Equation (23) is the element of the beamspace atomic set $\mathcal{B}$, i.e., the beamspace array manifold $\mathcal{B}$ is a semi-algebraic set as $\mathcal{D}\left(g\right)$, where $g_v\left(\mathbf{b}\right)$ is constructed from Equation (23) and $g_v\left(\mathbf{b}\right)\equiv 0$ for $\mathbf{b}\in \mathcal{B}$.  □

By utilizing Theorem 1 and Theorem 2, we present the following theorem.

Theorem 3.

Given $f\left(\mathbf{b}\right)$ to be a polynomial of $\mathbf{b}$, if $f\left(\mathbf{b}\right)\ge 0$ on the domain $\mathcal{B}$ with the DFT beamspace in Equation (4), $f\left(\mathbf{b}\right)$ is a sum-of-squares polynomial on $\mathbf{b}$.

Proof of Theorem 3.

Given a beamspace steering vector $\mathbf{b}\in \mathcal{B}$, we have $\left|b_m\right|\le 1$ followed by [4], which implies

$$\chi \left(\mathbf{b}\right)=M-\sum _{m=0}^M{b}_m^2\ge 0,\mathbf{b}\in \mathcal{B},$$

(24)

and thus $\mathcal{D}\left(g\right)$ is bounded, where $\mathcal{D}\left(g\right)$ is the semi-algebraic set as $\mathcal{D}\left(g\right)$ in Equation (19) corresponding to $\mathcal{B}$, followed by Theorem 2.

Using Theorem 1 and the fact that $g_v\left(\mathbf{b}\right)\equiv 0$ for $\mathbf{b}\in \mathcal{B}$, we can say that any $f\left(\mathbf{b}\right)\ge 0$ for $\mathbf{b}\in \mathcal{D}\left(g\right)$ satisfies

$$\begin{array}{cc}\hfill f\left(\mathbf{b}\right)& =\sum _{\alpha \in {\left\{0,1\right\}}^V}{g}_1^{{\alpha }_1}\cdots g_V^{{\alpha }_V}{s}_{\alpha },\hfill \\ & =s_{\mathbf{0}},s_{\alpha }\in \sum \mathbb{R}{\left[\mathbf{b}\right]}^2,\hfill \end{array}$$

(25)

where $s_{\mathbf{0}}$ denotes the $s_{\alpha }$ with ${\left\{{\alpha }_v\right\}}_{v=1,\cdots ,V}=0$, i.e., $f\left(\mathbf{b}\right)$ is sum-of-squares polynomial on $\mathbf{b}$.  □

We note that the result in Equation (25) is with the sum-of-squares relaxation of the polynomial in [26]. If $\mathrm{deg}\left(f\left(\mathbf{b}\right)\right)=\mathrm{deg}\left(s_0\right)$, we can have the SDP implementation to Equation (15) followed by [24] as

$$\begin{array}{c}\underset{\mathbf{P},\mathsf{\Omega }}{min}-Re\left\{\mathrm{Tr}\left[{\mathbf{y}}^{\mathrm{H}}\mathsf{\Omega }\right]\right\}\hfill \\ \mathrm{s}.\mathrm{t}.\mathrm{Tr}\left({\mathsf{\Theta }}_n{\mathbf{W}}^{\mathrm{H}}\mathbf{P}\mathbf{W}\right)=\mathrm{Tr}\left({\mathsf{\Theta }}_n{\mathbf{W}}^{\mathrm{H}}\mathbf{W}\right)\hfill \\ n\in H\hfill \\ \left[\begin{array}{cc}\mathbf{P}& -{\mathsf{\Omega }}^{\mathrm{H}}\\ -\mathsf{\Omega }& \mathbf{I}\end{array}\right]\ge 0\hfill \end{array}$$

(26)

where ${\mathsf{\Theta }}_n\in {\mathbb{R}}^{N\times N}$ is the zeros matrix except ones on the nth diagonal, H is a half space of $[-(N-1),N-1]$. The derivation of Equation (26) can be found in Appendix A.

Therefore, Equation (13) can be implemented by the following SDP problem with the standard Lagrangian analysis if strong duality holds

$$\begin{array}{c}\underset{\mathbf{Q},\mathbf{M}}{min}\frac{1}{2}\mathrm{Tr}\left(\mathbf{W}S\left(\mathbf{Q}\right){\mathbf{W}}^{\mathrm{H}}\right)+\frac{1}{2}\mathrm{Tr}\left(\mathbf{M}\right)\hfill \\ \mathrm{s}.\mathrm{t}.\left[\begin{array}{cc}\mathbf{W}S\left(\mathbf{Q}\right){\mathbf{W}}^{\mathrm{H}}& \mathbf{y}\\ {\mathbf{y}}^{\mathrm{H}}& \mathbf{M}\end{array}\right]\ge 0,\hfill \end{array}$$

(27)

where $S\left(\mathbf{Q}\right)$ is a Hermitian Toeplitz matrix with the first row vector $\mathbf{Q}$. The derivation of Equation (27) can be found in Appendix B.

Regarding the Theorem 3 and its resulting SDP implementations of BS-ANM problem in Equations (26)–(27), we remark as follows:

Remark 1.

It is well known that the positive trigonometric polynomial is sum-of-squares [24], i.e., $f\left(\mathbf{b}\right)\in \sum \mathbb{R}{\left[z\right]}^2$, which formulates the mathematical foundation of conventional ANM approaches, such as [12,17]. In Theorem 3, we further reveal that the positive polynomial $f\left(\mathbf{b}\right)\in \sum \mathbb{R}{\left[\mathbf{b}\right]}^2$ for the DFT beamspace manifold. Based on this theorem, we develop the lower dimensional SDP implementation in Equations (26) and (27), which is our main contribution.

Remark 2.

Compared with the conventional one in Equation (12), the proposed SDP implementation in Equation (27) is with lower dimensional SDP constraints in the beamspace processing scenarios ($N>M$), which results in the low computational complexity as shown in Section 5.

Remark 3.

The SDP formulation in Equation (26) is based on the sum-of-squares relaxation [26], and thus is an approximate algorithm of Equation (15). It can be regarded as the optimal polynomial designing to regression the dual polynomial in Equation (16) under the polynomial order constraint as shown in [27].

Remark 4.

Our work is related with the ANM approaches on the semi-algebraic set [22,23]. In particular, we formulate the beamspace array manifold as the boundary of a semi-algebraic set, which is obviously a semi-algebraic set either, and this results in the low dimensional SDP implementation.

4. BS-ANM via ADMM

Due to the SDP based methods with the high computational burden, we propose a fast implementation of BS-ANM approach via ADMM.

Let $\mathsf{\Lambda }\ge 0$ and $\mathbf{V}\ge 0$ be the Lagrangian multipliers of the constraint of Equation (27), the Lagrangian function of Equation (27) is

$$\begin{array}{c}L(\mathbf{Q},\mathbf{M},\mathsf{\Lambda })=\mathrm{Tr}\left(\mathbf{M}\right)+〈\mathbf{V}-\underset{\mathbf{T}}{\underbrace{\left[\begin{array}{cc}\mathbf{W}S\left(\mathbf{Q}\right){\mathbf{W}}^{\mathrm{H}}& \mathbf{y}\\ {\mathbf{y}}^{\mathrm{H}}& \mathbf{M}\end{array}\right]}},\mathsf{\Lambda }〉\hfill \\ +\frac{\rho }{2}{∥\mathbf{V}-\left[\begin{array}{cc}\mathbf{W}S\left(\mathbf{Q}\right){\mathbf{W}}^{\mathrm{H}}& \mathbf{y}\\ {\mathbf{y}}^{\mathrm{H}}& \mathbf{M}\end{array}\right]∥}_F^2+\mathrm{Tr}\left(\mathbf{W}S\left(\mathbf{Q}\right){\mathbf{W}}^{\mathrm{H}}\right)\hfill \end{array}$$

(28)

where $\mathsf{\Lambda }=\left[\begin{array}{cc}{\mathsf{\Lambda }}_0& \mathsf{\Omega }\\ {\mathsf{\Omega }}^{\mathrm{H}}& {\mathsf{\Lambda }}_1\end{array}\right]$, $\mathbf{V}=\left[\begin{array}{cc}{\mathbf{V}}_0& \mathsf{\Gamma }\\ {\mathsf{\Gamma }}^{\mathrm{H}}& {\mathbf{V}}_1\end{array}\right]$ and the blocks in $\mathsf{\Lambda }$ and $\mathbf{V}$ match the partition of the blocks in $\mathbf{T}$. Then, the partial derivatives of $L(\mathbf{Q},\mathbf{M},\mathsf{\Lambda })$ with respect to $\mathbf{Q}$ and $\mathbf{M}$ are calculated as follows:

$$\begin{array}{c}\frac{\partial L}{\partial \mathbf{M}}=\mathbf{I}-{\mathsf{\Lambda }}_1-\rho ({\mathbf{V}}_1-\mathbf{M}),\hfill \\ \frac{\partial L}{\partial \mathbf{Q}}={\mathbf{e}}_1-\mathfrak{E}({\mathsf{\Lambda }}_0,\mathbf{W})-\frac{\rho }{2}\mathfrak{E}({\mathbf{V}}_0,\mathbf{W})+\frac{\rho }{2}\mathfrak{F}\left(\mathbf{W}\right)\left[\begin{array}{c}{\mathbf{Q}}_{\mathrm{r}}\\ {\mathbf{Q}}_{\mathrm{i}}\end{array}\right],\hfill \end{array}$$

(29)

where ${\mathbf{Q}}_{\mathrm{r}}$ and ${\mathbf{Q}}_{\mathrm{i}}$ denote the real and imaginary part of $\mathbf{Q}$, ${\mathbf{e}}_1$ is the vector with zeros except the first element, $\mathbf{I}$ is the unity matrix. $\mathfrak{E}({\mathbf{V}}_0,\mathbf{W})\in {\mathbb{C}}^{N\times 1}$ is a vector function of ${\mathbf{V}}_0$ and $\mathbf{W}$, whose $n\mathrm{th}$ element is ${\left[\mathfrak{E}({\mathbf{V}}_0,\mathbf{W})\right]}_n=\mathrm{Tr}\left({\mathsf{\Theta }}_n{\mathbf{W}}^{\mathrm{H}}{\mathbf{V}}_0\mathbf{W}\right)$ and $\mathfrak{E}({\mathsf{\Lambda }}_0,\mathbf{W})$ is defined similarly. $\mathfrak{F}\left(\mathbf{W}\right)=\left[{\mathfrak{F}}_r\left(\mathbf{W}\right)-j{\mathfrak{F}}_i\left(\mathbf{W}\right)\right]$ is a matrix function of $\mathbf{W}$, where the element of ${\mathfrak{F}}_r\left(\mathbf{W}\right)\in {\mathbb{C}}^{N\times N}$ on $m\mathrm{th}$ row and $n\mathrm{th}$ column can be described as

$$\begin{array}{c}\hfill {\left[{\mathfrak{F}}_r\left(\mathbf{W}\right)\right]}_{m,n}=\mathrm{Tr}\left(\mathbf{W}{\mathsf{\Theta }}_n{\mathbf{W}}^{\mathrm{H}}{\mathsf{\Pi }}_m\right)\\ \hfill {\left[{\mathfrak{F}}_i\left(\mathbf{W}\right)\right]}_{m,n}=\mathrm{Tr}\left(\mathbf{W}{\mathsf{\Theta }}_n{\mathbf{W}}^{\mathrm{H}}{\mathsf{{\rm Y}}}_m\right)\end{array}$$

(30)

where

$$\begin{array}{cc}\hfill {\mathsf{\Pi }}_m& =\left\{\begin{array}{c}\mathbf{W}{\mathsf{\Theta }}_m{\mathbf{W}}^{\mathrm{H}}+\mathbf{W}{\mathsf{\Theta }}_{-m}{\mathbf{W}}^{\mathrm{H}},m\ne 0\hfill \\ \mathbf{W}{\mathbf{W}}^{\mathrm{H}},m=0,\hfill \end{array}\right.\hfill \\ \hfill {\mathsf{{\rm Y}}}_m& =\left\{\begin{array}{c}-\mathbf{W}{\mathsf{\Theta }}_m{\mathbf{W}}^{\mathrm{H}}+\mathbf{W}{\mathsf{\Theta }}_{-m}{\mathbf{W}}^{\mathrm{H}},m\ne 0\hfill \\ \mathbf{W}{\mathbf{W}}^{\mathrm{H}},m=0.\hfill \end{array}\right.\hfill \end{array}$$

(31)

Therefore, we can give the updates steps for the alternating optimization of the ADMM. Given the results of $l\mathrm{th}$ iteration ${\mathbf{Q}}^l,{\mathsf{\Lambda }}_0^l,{\mathsf{\Lambda }}_1^l,{\mathbf{V}}_1^l,{\mathbf{V}}_0^l$, the updates rules can be described as

$$\begin{array}{cc}\hfill {\mathbf{M}}^{l+1}& =\frac{1}{\rho }{\mathsf{\Lambda }}_1^l+{\mathbf{V}}_1^l-\frac{1}{\rho }\mathbf{I}\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill \left[\begin{array}{c}{\mathbf{Q}}_{\mathrm{r}}^{l+1}\\ {\mathbf{Q}}_{\mathrm{i}}^{l+1}\end{array}\right]& =\frac{2}{\rho }{\mathfrak{F}}^{†}\left(\mathbf{W}\right)\left[\mathfrak{E}({\mathsf{\Lambda }}_0^l,\mathbf{W})+\frac{\rho }{2}\mathfrak{E}({\mathbf{V}}_0^l,\mathbf{W})-{\mathbf{e}}_1\right]\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill {\mathsf{\Lambda }}^{l+1}& ={\mathsf{\Lambda }}^l+\rho ({\mathbf{V}}^{l+1}-{\mathbf{T}}^{l+1})\hfill \end{array}$$

(34)

with

$$\begin{array}{cc}\hfill {\mathbf{T}}^{l+1}& =\left[\begin{array}{cc}\mathbf{W}S\left({\mathbf{Q}}^{l+1}\right){\mathbf{W}}^{\mathrm{H}}& \mathbf{y}\\ {\mathbf{y}}^{\mathrm{H}}& {\mathbf{M}}^{l+1}\end{array}\right]\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill {\mathbf{V}}^{l+1}& ={\left[{\mathbf{T}}^{l+1}-\rho {\mathsf{\Lambda }}^l\right]}_{+}.\hfill \end{array}$$

(36)

In Equation (33), ${\left[•\right]}_{+}$ denotes the orthogonal projection on the positive semidefinite matrices, which is defined as

$$\begin{array}{cc}\hfill {\mathbf{V}}^{l+1}& ={\mathbf{U}}^l{\mathsf{\Sigma }}_{+}{{\mathbf{U}}^l}^{\mathrm{H}}\hfill \end{array}$$

(37)

where ${\mathbf{T}}^{l+1}-\rho {\mathsf{\Lambda }}^l={\mathbf{U}}^l\mathsf{\Sigma }{{\mathbf{U}}^l}^{\mathrm{H}}$ is the eigenvalue decomposition and ${\mathsf{\Sigma }}_{+}$ is $\mathsf{\Sigma }$ with all negative eigenvalue being zeros.

Summarizing above, we completed the iterations steps of ADMM algorithm for the BS-ANM problem in

Table 1. The ADMM based BS-ANM algorithm.

Initialize: ${\mathit{Q}}^0=0,{\mathit{M}}^0=0,{\mathit{\Lambda }}^0=0,{\mathit{V}}^0=0,\mathit{l}=0$
Input: y; Output: ${\mathit{V}}_{\mathit{o}}$
Iteration:
While (the stop condition is not satisfied) do:
1:   Update ${\mathbf{M}}^{l+1}$ and ${\mathbf{Q}}^{l+1}$ using Equations (32) and (33), respectively.
2:   Update ${\mathbf{T}}^{l+1}$ using Equation (35).
3:   Update ${\mathbf{V}}^{l+1}$ using Equation (36).
4:   Update ${\mathsf{\Lambda }}^{l+1}$ using Equation (34).
5:   $l=l+1$;
end
Output:${\mathbf{V}}_o={\mathbf{V}}_l$

. The iteration will stop until converging or achieving the iteration number limits.

5. Complexity Analysis

In this section, we provide the complexity analysis to the implementations of the ANM based methods in beamspace, such as the algorithms in Equations (12) and (27), Table 1 and [28].

5.1. Complexity of SDP Implementations

It is shown in [22] that a primal-dual path following method needs $C\left(ϵ\right)=-O\left(1\right)ln\left(ϵ\right)S$ flops to solve the SDP problem with accuracy $ϵ$, where

$$S={\left\{1+\sum _{i=1}^I{\beta }_i\right\}}^{\frac{1}{2}}\left\{{\gamma }^3+{\gamma }^2\sum _{i=1}^I{\beta }_i^2+\gamma \sum _{i=1}^I{\beta }_i^3\right\}$$

(38)

and the SDP problem is with $\gamma$ independent real variables and I real linear matrix inequalities (LMI) with the $i\mathrm{th}$ LMI size of $\beta \times \beta$.

Firstly, we analyze the GL based SDP implementation in Equation (12). Here, we have $\gamma =4N-2M,I=1$ and ${\beta }_1=2(N+1)$, and thus $S=O\left(N^{4.5}\right)$ for Equation (12), followed by [22].

For the SDP implementation of BS-ANM approach in Equation (27), there are $N-1$ complex variables, two real variables, and a $(M+1)\times (M+1)$ SDP constraint, thus we have $\gamma =2N,I=1$ and ${\beta }_1=2(M+1)$, which gives

$$S=\sqrt{2M+3}\left\{8N^3+16N^2(M+1)+16N{(M+1)}^3\right\}.$$

(39)

Accordingly, for $N>M$, the computational complexity of Equation (27) is $O\left(N{M}^{3.5}\right)$, compared with the computational complexity of $S=O\left(N^{4.5}\right)$ for the GL based method.

5.2. Complexity of ADMM Based Methods

The ADMM based method for generalized line spectral estimation problem has been proposed in [28]. As shown in [28], the step with highest computational burden of the ADMM methods is the projection onto the semidefinite cone as Equation (34), and thus the ADMM based method for GL problem is with the complexity of $O\left(N^3\right)$, since the eigenvalue decomposition requires $O\left(N^3\right)$ operations.

Similarly, the algorithm proposed in Section 4 is with a semidefinite cone projection operator of $M\times M$ matrix, and the complexity of the proposed ADMM method is $O\left(M^3\right)$ accordingly. It is noted that the complexity of the MUSIC [29] is $O\left(M^3\right)$, although the proposed ADMM method has a higher computational cost than the MUSCI method since the iterations.

Therefore, the proposed methods are significantly more computationally efficient than GL based ones, when $N>M$. The complexity gain of the proposed methods is $O{\left(\frac{M}{N}\right)}^{3.5}$ for SDP implementation and $O{\left(\frac{M}{N}\right)}^3$ for ADMM implementation, respectively.

6. DOA Estimation via BS-ANM

Based on the BS-ANM approaches proposed in Section 3 and Section 4, we develop the covariance matrix based DOA estimation algorithm in free space and the channel estimation method for massive MIMO system with lens antennas.

6.1. DOA Estimation with Covariance Matrix

Rewriting the MMV model in Equation (7) that $\mathbf{R}={\displaystyle \sum _{k=1}^K}{c}_k\mathbf{b}\left({\theta }_k\right){\mathbf{h}}_k$, we define the beamspace atomic set as

$${\mathcal{B}}^{\prime }:=\left\{{\mathbf{b}}^{\prime }\left(\theta \right)\mathbf{h}\left|{\mathbf{b}}^{\prime }\left(\theta \right)=\mathbf{Wa}\left(\theta \right)/{∥\mathbf{Wa}\left(\theta \right)∥}_F,\theta \in \mathsf{\Xi },{∥\mathbf{h}∥}_2=1\right.\right\}.$$

(40)

Then, the beamspace atomic ${\ell }_1$ norm of the covariance matrix $\mathbf{R}$ can be defined as

$$\begin{array}{c}\hfill {∥\mathbf{R}∥}_{\mathcal{B}}^{\prime }=inf\left\{\left.\sum _{k=1}^K{c}_k\right|\mathbf{R}=\sum _{k=1}^K{c}_k{\mathbf{b}}^{\prime }\left({\theta }_k\right){\mathbf{h}}_k.\right\}.\end{array}$$

(41)

Followed by Equation (27), the SDP implementation of Equation (41) can be written as

$$\begin{array}{cc}\hfill & \underset{\mathbf{Q},\mathbf{M}}{min}\frac{1}{2}\mathrm{Tr}\left(\mathbf{W}S\left(\mathbf{Q}\right){\mathbf{W}}^{\mathrm{H}}\right)+\frac{1}{2}\mathrm{Tr}\left(\mathbf{M}\right)\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.\left[\begin{array}{cc}\mathbf{W}S\left(\mathbf{Q}\right){\mathbf{W}}^{\mathrm{H}}& \mathbf{R}\\ {\mathbf{R}}^{\mathrm{H}}& \mathbf{M}\end{array}\right]\ge 0.\hfill \end{array}$$

(42)

In practice DOA estimation applications, we usually estimate the noisy covariance matrix as

$$\begin{array}{cc}\hfill \widehat{\mathbf{R}}& =\frac{1}{J}\sum _{t=1}^J\mathbf{Y}\left(t\right){\mathbf{Y}}^{\mathrm{H}}\left(t\right),\hfill \\ \hfill & =\mathbf{B}{\mathbf{R}}_s{\mathbf{B}}^{\mathrm{H}}+\mathit{\eta },\hfill \end{array}$$

(43)

where $\mathit{\eta }$ is the outlier term. Consequently, the BS-ANM approach can be applied to the noisy case via regularization:

$$\begin{array}{cc}\hfill & \underset{\mathbf{Q},\mathbf{M},\mathbf{R}}{min}\frac{1}{2}\mathrm{Tr}\left(\mathbf{W}S\left(\mathbf{Q}\right){\mathbf{W}}^{\mathrm{H}}\right)+\frac{1}{2}\mathrm{Tr}\left(\mathbf{M}\right)\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.\left[\begin{array}{cc}\mathbf{W}S\left(\mathbf{Q}\right){\mathbf{W}}^{\mathrm{H}}& \mathbf{R}\\ {\mathbf{R}}^{\mathrm{H}}& \mathbf{M}\end{array}\right]\ge 0,\hfill \\ \hfill & {∥\widehat{\mathbf{R}}-\mathbf{R}∥}_2^2\le {\epsilon }^2\hfill \end{array}$$

(44)

where $\epsilon$ can be set as $\epsilon =M{\sigma }^2$, followed by [30]. Here, we estimate $\sigma$ with the square root of the smallest eigenvalue of $\widehat{\mathbf{R}}$. Then, the DOAs can be obtained by employing Root-MUSIC [31] to $\mathbf{W}S\left({\mathbf{Q}}^{☆}\right){\mathbf{W}}^{\mathrm{H}}$, where ${\mathbf{Q}}^{☆}$ is the solution of Equation (44).

Summarizing above, we list the steps of the DOA estimation algorithm via the BS-ANM approach in

Table 2. DOA estimation algorithm via BS-ANM.

Inpute: $\mathit{Y}\left(\mathit{t}\right),\mathit{K}$
Output: ${\left\{{\widehat{\mathit{\theta }}}_{\mathit{k}}\right\}}_{\mathit{k}=1,2,\cdots ,\mathit{K}}$
1:   Compute the covariance matrix $\widehat{\mathbf{R}}$ as in Equation (43).
2:   Compute the smallest of eigenvalue of $\widehat{\mathbf{R}}$ to obtain $\epsilon$.
3:   Solve Equation (44) to obtain $\mathbf{W}S\left({\mathbf{Q}}^{☆}\right){\mathbf{W}}^{\mathrm{H}}$.
4:   Apply Root-MUSIC to $\mathbf{W}S\left({\mathbf{Q}}^{☆}\right){\mathbf{W}}^{\mathrm{H}}$ and estimate ${\left\{{\widehat{\theta }}_k\right\}}_{k=1,2,\cdots ,K}$.

.

To apply the ADMM approach to DOA estimation, we propose a subspace based method, inspired from [32]. Consider the covariance matrix in Equation (43)

$$\widehat{\mathbf{R}}={\mathbf{U}}_s{\mathsf{\Sigma }}_s{\mathbf{U}}_s^{\mathrm{H}}+{\sigma }^2{\mathbf{U}}_n{\mathbf{U}}_n^{\mathrm{H}},$$

(45)

where ${\mathbf{U}}_s$ and ${\mathbf{U}}_n$ are the signal and noise subspace, respectively, which is obtained by eigenvalue decomposition. Let $\tilde{\mathbf{y}}={\mathbf{U}}_s{\mathbf{U}}_s^{†}$ be substituted into Equations (32), (33), and (34) to replace $\mathbf{y}$; the proposed ADMM method can be applied to DOA estimation in noisy environments. We list the algorithm in

Table 3. The DOA estimation algorithm via ADMM.

Input: $\mathit{Y}\left(\mathit{t}\right),\mathit{K}$
Output: ${\left\{{\widehat{\mathit{\theta }}}_{\mathit{k}}\right\}}_{\mathit{k}=1,2,\cdots ,\mathit{K}}$
1:   Compute the covariance matrix $\widehat{\mathbf{R}}$ as in Equation (43).
2:   Compute the eigenvalue decomposition of $\widehat{\mathbf{R}}$ and obtain ${\mathbf{U}}_s$ as the eigenvectors corresponding to the large eigenvalues.
3:   Let $\mathbf{y}={\mathbf{U}}_s{\mathbf{U}}_s^{†}$ and compute ${\mathbf{V}}_o$ by the algorithm in Table 1.
4:   Apply Root-MUSIC to ${\mathbf{V}}_o$ to estimate ${\left\{{\widehat{\theta }}_k\right\}}_{k=1,2,\cdots ,K}$.

.

6.2. Channel Estimation with Lens Antenna Array

In the millimeter-wave massive MIMO system, the lens antenna arrays with the DFT beamformer are used for the hardware cost and power reduction [9,10]. In this subsection, we apply the BS-ANM approach to the channel estimation at the base station (BS) in such scenario.

Considering a narrowband millimeter-wave massive MIMO system with N antennas ULA at base station and a single antenna mobile phone, the uplink channel can be expressed as [11]

$$\mathbf{H}\left(t\right)=\sum _{k=1}^K{g}_k\left(t\right)\mathbf{a}\left({\theta }_k\right),$$

(46)

where $\mathbf{a}\left(\theta \right)$ is the steering vector of ULA at BS in Equation (3) and $g_k\left(t\right)$ denotes the channel gain for the $k\mathrm{th}$ propagation path at time block t.

As shown in [9], given a lens antenna array with the M channel DFT beamformer, the received signal after normalizing the known pilot sequence can be given by

$$\begin{array}{cc}\hfill \mathbf{y}\left(t\right)& =\mathbf{B}\mathbf{g}\left(t\right)+{\mathbf{n}}_{\mathbf{b}}\left(t\right)\hfill \end{array}$$

(47)

where $\mathbf{g}\left(t\right)={[g_1\left(t\right),\cdots ,g_K\left(t\right)]}^{\mathrm{T}}$, $\mathbf{B}$ is the beamspace array manifold in Equation (5), then the channel estimation can be regarded as the sparse representation problem with the MMV model.

Thus, the estimated angle of the scatters $\left\{{\widehat{\theta }}_k\right\}$ can be obtained by the algorithm presented in Section 6.1 and the path gain $\left\{{\mathbf{g}}_k\left(t\right)\right\}$ can be given by

$$\begin{array}{cc}\hfill & \widehat{\mathbf{g}}\left(t\right)=\underset{g}{min}{∥\mathbf{y}\left(t\right)-\mathbf{Bg}\left(t\right)∥}_2\hfill \end{array}$$

(48)

Finally, we have the estimated channel state information (CSI) as

$$\widehat{\mathbf{H}}\left(t\right)=\sum _{k=1}^K{\widehat{g}}_k\left(t\right)\mathbf{a}\left({\widehat{\theta }}_k\right).$$

(49)

The steps of the BS-ANM based channel estimation algorithm for the lens antenna array are presented in

Table 4. Channel estimation algorithm via BS-ANM.

Input: $\mathit{y}\left(\mathit{t}\right)$
Output: $\widehat{\mathit{H}}\left(\mathit{t}\right)$
1:   Compute the covariance matrix $\widehat{\mathbf{R}}$ with $\mathbf{y}\left(t\right)$ as in Equation (43).
2:   Solve Equation (44) to obtain $\mathbf{W}S\left({\mathbf{Q}}^{☆}\right){\mathbf{W}}^{\mathrm{H}}$.
3:   Compute the eigenvalue of $\mathbf{W}S\left({\mathbf{Q}}^{☆}\right){\mathbf{W}}^{\mathrm{H}}$ as ${\gamma }_1\ge {\gamma }_2\ge \cdots \ge {\gamma }_M$, and estimate the number of scatters by $K=\min [\eta ,\sup \left\{k:{\gamma }_k\ge 0.01{\gamma }_1\right\}]$.
4:   Apply Root-MUSIC to ${\mathbf{V}}_o$ to estimate ${\left\{{\widehat{\theta }}_k\right\}}_{k=1,2,\cdots ,K}$.
5:   Estimate the path gain $\widehat{\mathbf{g}}\left(t\right)$ using Equation (48).
6:   Estimate $\widehat{\mathbf{H}}\left(t\right)$ by Equation (49).

.

7. Simulations

In this section, we evaluate the performance of the proposed methods through simulations. Here, the proposed SDP and ADMM based methods are referred to as BS-ANM and BS-ADMM, respectively.

Firstly, we compare the resolution capability of the BS-ANM and the generalized line spectral estimation based ANM method (GL-ANM). Consequently, we demonstrate the computational complexity and the DOA estimation performance comparison of the proposed the BS-ANM and BS-ADMM methods with the conventional ANM and subspace based methods, such as Beamspace MUSIC in [29], SPA in [13], and GL-ANM in [17]. As a performance metric of DOA estimation accuracy, the Cramer–Rao bound (CRB) derived in [33] is also presented in simulations. In Section 7.1, Section 7.2 and Section 7.3, simulations perform the DOA estimation in free space, where the signal model is followed from Equation (5). In Section 7.4, we show the performance of the proposed methods for the channel estimation problem in the massive MIMO system as shown in Figure 1.

In simulations of DOA estimation, the root mean square error (RMSE) is utilized to evaluate the performance of algorithms as

$$\mathrm{RMSE}=\sqrt{\frac{1}{TK}\sum _{i=1}^K\sum _{j=1}^T{({\stackrel{⌢}{\theta }}_{ij}-{\theta }_i)}^2}$$

(50)

where ${\stackrel{⌢}{\theta }}_{ij}$ and ${\theta }_i$ denote the estimated and the true DOA of the ith target at jth trial, the number of trials is T and K is the number of sources. Monte Carlo simulations are are performed with $T=200$ to compute RMSE using Intel i7-8700K CPU PC (Santa Clara, CA, USA). All the simulations are implemented by Matlab2017b (Natick, MA, USA). Our code is publicly available at github.com/panda-1982/BS-ANM.

7.1. Comparison of Resolution and Complexity

Consider a ULA with $N=22$ elements and $M=10$ DFT beamformer, $d=\lambda /2$, there are two sources at $\theta ={-20}^{\circ }$ and $\theta ={-20+\delta }^{\circ }$. The RMSE of the BS-ANM method and GL-ANM method in a noiseless case with 100 snapshots are shown in Figure 2 versus angle separation. It can be seen that the BS-ANM method performs similar resolution capacity with the GL-ANM method.

The computing time of GL-ANM, BS-ANM, BS-ADMM, SPA, and MUSIC methods are presented in Figure 3. In Figure 3a, the dimension of beamspace is fixed as $M=10$, and we compare the CPU runtime of evaluated methods versus the number of sensors N. In Figure 3b, the number of sensors is set as $N=64$ and the CPU runtime of evaluated methods is compared versus the dimension of beamspace M. In simulations, the GL-ANM, BS-ANM, and SPA methods are implemented by the CVX software package [34], and the BS-ADMM is implemented in Matlab with the "codegen" feature. As shown from the results, the runtime of GL-ANM is significantly higher than BS-ANM and BS-ADMM, especially when the number of sensors N is large.

7.2. Performance for Uncorrelated Sources

In this subsection, the proposed BS-ANM based methods are evaluated with some baselines in the presence of uncorrelated sources.

Supposing that the array configuration is the same as the setting in Section 7.1, there are two independent sources at $\theta ={-10}^{\circ }$ and $\theta =5^{\circ }$ impinging onto the array and 200 snapshots are collected in each trial. Figure 4 shows the RMSE of DOA estimation of the evaluated algorithms when SNR varies from −15 to 10 dB. It can be seen that the proposed BS-ANM and BS-ADMM methods have almost the identical RMSE to GL-ANM method and perform better than the other competitors. We also note that the performance of MUSIC method degrades significantly compared with the sparsity based methods when SNR is lower than −8 dB.

When the SNR is fixed at SNR = −10 dB, Figure 5 plots the RMSE of DOA estimation of the evaluated algorithms versus a different number of snapshots. The result shows that the GL-ANM method performs the best among the algorithms, and the BS-ANM method has a similar performance to the GL-ANM method. When the number of snapshot is not too low, i.e., larger than 100, the BS-ADMM method also shows the satisfying DOA estimation accuracy. In addition, it is shown that the MUSIC method suffers from the lack of snapshots more seriously than the sparsity based methods.

In Figure 6, we demonstrate the DOA estimation performance versus the number of sources. Considering the ULA with $N=22$ elements and $M=8$ DFT beamformer, there are K uncorrelated sources at $[{-30}^{\circ },{-24}^{\circ },\cdots ,{(6K-36)}^{\circ }]$ and $J=100$ snapshots are used for DOA estimation. The SNR is fixed at $0\mathrm{dB}$. We can see from the simulation results that the proposed BS-ANM and BS-ADMM methods are robust to the multiple sources’ scenarios.

7.3. Performance for Correlated Sources

Next, we study the effect of correlated sources on the evaluated algorithms. Supposing that the array configuration and the source location are set as in Figure 4, the signal sequences consisted of 200 snapshots in each trial with the correlation coefficient of $\zeta =0.5$. We present the RMSE of the evaluated algorithms versus SNR in Figure 7. As we can see from Figure 7, GL-ANM, BS-ANM, and BS-ADMM methods have almost the same performance, which is very close to CRB when SNR is higher than −10 dB. In addition, the SPA and MUSIC methods are suffered from the correlated sources more seriously than the ANM based ones, especially in low SNR.

To evaluate the effect of the number of snapshots in a correlated sources’ scenario, we set $\mathrm{SNR}=0\mathrm{dB}$ and the correlation coefficient of $\zeta =0.7$, and the RMSE of the evaluated algorithms versus the number of snapshots is demonstrated in Figure 8. The simulation results show that the proposed BS-ANM and BS-ADMM methods provide the significant robustness against the correlation of the sources compared with SPA and MUSIC methods and can achieve satisfying DOA estimation performance close to CRB even with a low number of snapshots.

7.4. Performance for Channel Estimation

We verify the channel estimation method proposed in Section 6.2 with lens antenna array for massive MIMO system as shown in Figure 1.

Suppose the $N=64$ ULA with $M=16$ dimensions DFT beamspace as the receiver antennas array and the single antenna transmitter consist of the MIMO system, there are $K=6$ random distributed scatters in the main-lobe region of the receiver array whose locations are satisfying uniform distribution, and the path gains $\left\{g_k\left(t\right)\right\}$ are assumed as independent identical Gaussian distribution. The channel estimation is carried out with $J=10$ snapshots and $\eta =9$ in each trail. The 200 Monte Carlo trials are performed in the simulation. The normalized mean squared error (NMSE) of the CSI estimation

$$\mathrm{NMSE}=\frac{1}{L}\sum _{l=1}^L\frac{{∥{\widehat{\mathbf{H}}}_l\left(t\right)-\mathbf{H}\left(t\right)∥}_F^2}{{∥\mathbf{H}\left(t\right)∥}_F^2}$$

(51)

is utilized to evaluate the performance of the channel estimation methods, where ${\widehat{\mathbf{H}}}_l\left(t\right)$ and $\mathbf{H}\left(t\right)$ are the estimated and ground truth CSI, respectively.

The NMSE of the orthogonal matching pursuit (OMP) based channel estimation method in [8], the GL based channel estimation method in [20] and the proposed method in Section 6.2 versus SNR is shown in Figure 9. The simulation results indicate that the proposed method shows the significant performance improvement over the OMP based method and has almost the same performance as the GL based method.

8. Conclusions

In this paper, we study the low complexity implementation of the beamspace atomic norm minimization for DOA estimation. We propose the beamspace atomic norm and its low dimensional SDP implementation for DFT beamspace. By utilizing this approach, the covariance based BS-ANM DOA estimation methods and BS-ANM based channel estimation method for a massive MIMO system with a lens antenna array are developed. The complexity analysis and simulations indicate that the proposed methods have almost the same performance and significantly lower computational complexity than the generalized line spectral estimation based ANM methods. In addition, the proposed methods demonstrate the performance improvement compared to some state-of-the-art counterparts.