Super-Resolution Channel Estimation for Arbitrary Arrays in Hybrid Millimeter-Wave Massive MIMO Systems

TL;DR

This paper introduces a super-resolution channel estimation method for hybrid millimeter-wave massive MIMO systems with arbitrary array geometries, using a generalized array manifold separation and efficient algorithms to improve accuracy and speed.

Contribution

It proposes a novel angle-based super-resolution channel estimation technique that handles arbitrary array geometries and hybrid hardware constraints in mmWave MIMO systems.

Findings

Effective super-resolution angle estimation in arbitrary array geometries.

Significant reduction in sensing time for channel acquisition.

Improved accuracy demonstrated through simulations.

Abstract

This paper develops efficient channel estimation techniques for millimeter-wave (mmWave) massive multiple-input multiple-output (MIMO) systems under practical hardware limitations, including an arbitrary array geometry and a hybrid hardware structure. Taking on an angle-based approach, this work adopts a generalized array manifold separation approach via Jacobi-Anger approximation, which transforms a non-ideal, non-uniform array manifold into a virtual array domain with a desired uniform geometric structure to facilitate super-resolution angle estimation and channel acquisition. Accordingly, structure-based optimization techniques are developed to effectively estimate both the channel covariance and the instantaneous channel state information (CSI) within a short sensing time. In particular, the difference in time-variation of channel path angles and path gains is capitalized to design a two-step CSI estimation scheme that can quickly sense fading channels. Theoretical results are provided on the fundamental limits of the proposed technique in terms of sample efficiency. For computational efficiency, a fast iterative algorithm is developed via the alternating direction method of multipliers. Other related issues such as spurious-peak cancellation in nonuniform linear arrays and extensions to higher-dimensional cases are also discussed. Simulations testify the effectiveness of the proposed approaches in hybrid mmWave massive MIMO systems with arbitrary arrays.

Equations

$$\begin{array}{r}{\mathbf{h}}_t=\sum _{k=1}^K{g}_{k,t}\mathbf{a}({\theta }_k),\forall t,\end{array}$$

$$\begin{array}{r}[\mathbf{a}({\theta }_k){]}_n=e^{j2\pi \frac{r_n}{\lambda }\cos ({\theta }_k-{\varphi }_n)},\end{array}$$

$${\mathbf{h}}_t=\mathbf{A}{\mathbf{g}}_t,$$

$$\begin{array}{r}{\mathbf{x}}_t={\mathbf{h}}_t{z}_t+{\mathbf{w}}_t=\mathbf{A}{\mathbf{g}}_t{z}_t+{\mathbf{w}}_t,\end{array}$$

$$\begin{array}{r}{\mathbf{R}}_x=\mathbb{E}\{{\mathbf{x}}_t{\mathbf{x}}_t^H\}={\mathbf{R}}_h+{\sigma }^2\mathbf{I}=\mathbf{A}{\mathbf{R}}_g{\mathbf{A}}^H+{\sigma }^2\mathbf{I},\end{array}$$

$$\begin{array}{r}{\mathbf{y}}_t=z_t^{\ast }{\mathbf{W}}_t{\mathbf{x}}_t={\mathbf{W}}_t\mathbf{A}{\mathbf{g}}_t+{\mathbf{W}}_t{\mathbf{n}}_t,\end{array}$$

$$\begin{array}{r}[\mathbf{a}({\theta }_k){]}_n=e^{j2\pi (n-1)\frac{d}{\lambda }\cos ({\theta }_k)},\end{array}$$

$$\begin{array}{r}\mathbf{Y}=\mathbf{W}\mathbf{A}\mathbf{G}+\mathbf{W}\mathbf{N}=\mathbf{W}\mathbf{H}+\mathbf{W}\mathbf{N},\end{array}$$

$$\begin{array}{r}{\mathbf{R}}_y=\mathbb{E}\{{\mathbf{y}}_t{\mathbf{y}}_t^H\}=\mathbf{W}({\mathbf{R}}_h+{\sigma }^2\mathbf{I}){\mathbf{W}}^H.\end{array}$$

$$\begin{array}{r}\mathcal{A}=\left\{\mathbf{a}(f){\mathbf{b}}^H|f\in \left(-\frac{1}{2},\frac{1}{2}\right],\mathbf{b}\in {\mathbb{C}}^T,{\Vert \mathbf{b}\Vert }_2=1\right\},\end{array}$$

$$\begin{array}{r}{\Vert \mathbf{H}\Vert }_{\mathcal{A}}=\inf \left\{l>0:\mathbf{H}\in l\text{conv}(\mathcal{A})\right\},\end{array}$$

$$\begin{array}{r}\hat{\mathbf{H}}=\arg \underset{\mathbf{H}}{\min }{\Vert \mathbf{H}\Vert }_{\mathcal{A}}+\frac{\tau }{2}{\Vert \mathbf{Y}-\mathbf{W}\mathbf{H}\Vert }_F^2,\end{array}$$

$$\begin{array}{rl}{\hat{\mathbf{R}}}_h=& \arg \underset{\text{T}({\mathbf{u}}_h)}{\min }\text{tr}\left(\text{T}({\mathbf{u}}_h)\right)+\frac{\tau }{2}{\Vert {\hat{\mathbf{R}}}_y-\mathbf{W}\text{T}({\mathbf{u}}_h){\mathbf{W}}^H\Vert }_F^2\\ & \text{s.t.}\text{T}({\mathbf{u}}_h)⪰0.\end{array}$$

$$\begin{array}{r}[\mathbf{a}({\theta }_k){]}_n=\sum _{i=-\infty }^{+\infty }{j}^i{\text{J}}_i\left(2\pi \frac{r_n}{\lambda }\right)e^{-j{\varphi }_{n}i}{e}^{j{\theta }_{k}i},\end{array}$$

$$\begin{array}{r}I>\frac{2\pi }{\lambda }{r}_{max},\end{array}$$

$$\begin{array}{rl}[\mathbf{a}({\theta }_k){]}_n& \approx \sum _{i=-I}^I{j}^i{\text{J}}_i(2\pi \frac{r_n}{\lambda })e^{-j{\varphi }_{n}i}{e}^{j{\theta }_{k}i}\\ & ={\mathbf{c}}_n^T\mathbf{v}({\theta }_k),\end{array}$$

$$[{\mathbf{c}}_n{]}_i=j^i{\text{J}}_i(2\pi \frac{r_n}{\lambda })e^{-j{\varphi }_{n}i},i=-I,\dots ,I;$$

$$[\mathbf{v}({\theta }_k){]}_i=e^{j{\theta }_{k}i},i=-I,\dots ,I.$$

$$\begin{array}{r}\mathbf{a}({\theta }_k)=[{\mathbf{c}}_1,\dots ,{\mathbf{c}}_N{]}^T\mathbf{v}({\theta }_k)=\mathbf{C}\mathbf{v}({\theta }_k).\end{array}$$

$$\mathbf{V}\left(\mathbf{\theta }\right)=\left[\mathbf{v}({\theta }_1),\dots ,\mathbf{v}({\theta }_K)\right].$$

$$\mathbf{A}=\mathbf{CV}.$$

$$\begin{array}{rl}{\mathbf{h}}_t& =\mathbf{C}\mathbf{V}{\mathbf{g}}_t=\mathbf{C}{\mathbf{d}}_t,\end{array}$$

$$\begin{array}{rl}{\mathbf{R}}_h& =\mathbf{C}\mathbf{V}{\mathbf{R}}_g{\mathbf{V}}^H{\mathbf{C}}^H=\mathbf{C}{\mathbf{R}}_v{\mathbf{C}}^H,\end{array}$$

$${\mathbf{R}}_v=\mathbf{V}{\mathbf{R}}_g{\mathbf{V}}^H.$$

$$\begin{array}{rl}{\mathbf{y}}_t& ={\mathbf{W}}_t\mathbf{C}\mathbf{V}{\mathbf{g}}_t+{\mathbf{W}}_t{\mathbf{n}}_t\\ & ={\mathbf{\Phi }}_t{\mathbf{d}}_t+{\mathbf{W}}_t{\mathbf{n}}_t,\end{array}$$

$$\begin{array}{rl}{\mathbf{R}}_y& =\mathbf{W}\mathbf{C}{\mathbf{R}}_v{\mathbf{C}}^H{\mathbf{W}}^H+{\sigma }^2\mathbf{W}{\mathbf{W}}^H\\ & =\mathbf{\Phi }{\mathbf{R}}_v{\mathbf{\Phi }}^H+{\sigma }^2\mathbf{W}{\mathbf{W}}^H,\end{array}$$

$$\begin{array}{rl}{\hat{\mathbf{R}}}_v=& \arg \underset{\text{T}({\mathbf{u}}_v)}{\min }\text{tr}\left(\text{T}({\mathbf{u}}_v)\right)+\frac{\tau }{2}{\Vert {\hat{\mathbf{R}}}_y-\mathbf{\Phi }\text{T}({\mathbf{u}}_v){\mathbf{\Phi }}^H\Vert }_F^2\\ & \text{s.t.}\text{T}({\mathbf{u}}_v)⪰0.\end{array}$$

$$\begin{array}{r}{\mathcal{A}}^{\prime }=\left\{\mathbf{v}(f){\mathbf{q}}^H|f\in (-\frac{1}{2},\frac{1}{2}],\mathbf{q}\in {\mathbb{C}}^T,{\Vert \mathbf{q}\Vert }_2=1\right\},\end{array}$$

$$\begin{array}{r}\hat{\mathbf{D}}=\arg \underset{\mathbf{D}}{\min }\Vert \mathbf{D}{\Vert }_{{\mathcal{A}}^{\prime }}+\frac{\tau }{2}{\Vert \mathbf{Y}-\mathbf{\Phi }\mathbf{D}\Vert }_F^2.\end{array}$$

$$\begin{array}{r}\hat{\mathbf{D}}=\arg \underset{\mathbf{D}}{\min }\Vert \mathbf{D}{\Vert }_{{\mathcal{A}}^{\prime }}+\frac{\tau }{2}\sum _{t=1}^T{\Vert {\mathbf{y}}_t-{\mathbf{\Phi }}_t{\mathbf{d}}_t\Vert }_F^2.\end{array}$$

$$\displaystyle\begin{split}(\hat{\bm{D}},&\text{T}(\hat{\bm{u}}_{d}),\hat{\bm{Q}})=\\ \arg&\min_{\bm{D},\text{T}(\bm{u}_{d}),\bm{Q}}~{}~{}~{}\frac{1}{2\sqrt{2I+1}}\left(\text{tr}\left(\text{T}\left(\bm{u}_{d}\right)\right)+\text{tr}\left(\bm{Q}\right)\right)\\ &+\frac{\tau}{2}\sum_{t=1}^{T}\left\lVert\bm{y}_{t}-\bm{\Phi}_{t}\bm{d}_{t}\right\rVert_{F}^{2},~{}~{}~{}\text{s.t.}~{}\left[\begin{array}[]{cc}\text{T}(\bm{u}_{d})&\bm{D}\\ \bm{D}^{H}&\bm{Q}\end{array}\right]\succeq 0.\\ \end{split}$$

Full text

I Introduction

In millimeter-wave (mmWave) massive multiple-input multiple-output (MIMO) communications, large antenna gains coupled with the availability of large bandwidths bring many desired benefits such as high throughput, large capacity, and robustness against fading and interference [1, 2], which all hinge on accurate channel knowledge. However, the increase in antennas results in an enlarged channel dimension that gives rise to challenges to traditional channel estimation techniques, in terms of the high signal acquisition cost and the large training overhead [3, 4]. The hinderance in channel estimation is further aggravated by practical hardware limitations. A hybrid analog-digital architecture is widely suggested for massive MIMO transceivers, which reduces the number of radio frequency (RF) chains by balancing between the analog RF part and the digital baseband part [5, 6]. However, under such a hybrid structure, the channel estimator at baseband can only observe a compressed representation of the channel through a few RF chains.

To overcome these challenges, compressive sensing (CS) has been advocated for channel estimation in mmWave massive MIMO systems [7, 8, 9, 10, 11, 12, 13]. These CS-based approaches exploit the channel sparsity that stems from the limited scattering characteristics of mmWave propagation [14, 15, 16, 17]. Through virtual channel modeling [18], the large-dimensional mmWave massive MIMO channels can be represented by only a small number of parameters, including the angles of departure/arrival (AoD/AoA) and the path gains of the sparse scattering paths. Therefore, CS techniques enable channel estimation from a small set of compressively collected training samples. In [7, 8], a sparse multipath channel is formulated as a sparse vector on the angle-delay-Doppler space, and then CS techniques are applied to recover the vectorized sparse channel. In [9], an adaptive CS-based algorithm is proposed to estimate the sparse channel with a hybrid analog-digital hardware architecture. In [10], a hybrid architecture based on phase shifters is proposed to recover the sparse channel via greedy search algorithms. To further reduce the power consumption of phase shifters, a switch-based hybrid architecture is developed in [11], for sparse channel estimation. In [12], the CS-based channel estimation scheme is extended to broadband mmWave MIMO systems. In [13], to reduce the problem complexity, CS-based channel estimation is divided into angle estimation and path gain estimation subproblems, which are solved sequentially.

All the aforementioned techniques aim to estimate the instantaneous channel state information (CSI). Another line of work focuses on estimating the channel statistics, such as the channel covariance [19, 20]. The channel covariance is an important second-order statistic, which remains constant over many channel coherence intervals and therefore can be used for statistics-based design of the precoders, beamformers and linear receivers [21, 22]. To estimate the second-order statistics of the vectorized sparse mmWave MIMO channel, a diagonal-search orthogonal matching pursuit algorithm is developed in [19], which not only utilizes the joint sparsity represented by the available multiple measurement vectors (MMV) but also takes advantage of the Hermitian structure of the channel covariance matrix. In [20], a CS-based channel covariance estimator is proposed by using dynamic sensing schemes and designing dynamic greedy pursuit algorithms for the hybrid architecture.

Existing CS-based channel estimators critically hinge on an on-grid assumption that the values of the AoD/AoA of each propagation path exactly reside on some predefined grid in the angular domain. However, in practice, the AoD/AoA of paths are continuously-valued off grid. As a result, CS-based methods suffer from degraded performance due to the power leakage effect around the recovered discrete grid points, a.k.a., the infamous basis mismatch problem [23]. An angle rotation technique is proposed to alleviate this problem, which is developed upon the spatial basis expansion model [24, 25]. It improves the estimation accuracy in the angular domain, but still experiences finite resolution due to some predefined spatial rotation parameters. A continuous basis pursuit technique is proposed for perturbed CS in [26], which is however limited through the series expansion. On the other hand, classical subspace methods, such as MUSIC and ESPRIT, can achieve super-resolution in angle estimation [27, 28]. But, they require a large number of snapshots for collecting sample statistics, which leads to a long sensing time and consumes large training resources. To circumvent the on-grid assumption required by traditional CS and achieve super-resolution at short sensing time, a gridless CS technique is developed via atomic norm minimization (ANM) in the form of semidefinite programming (SDP) [29, 30, 31]. As a structure-based optimization technique, gridless CS is applied for super-resolution channel estimation in mmWave massive MIMO systems [32, 33, 34], which utilizes not only the sparsity of the channels but also the Vandermonde structure of the antenna arrays.

By capitalizing on the critical Vandermonde structure, gridless CS implicitly assumes the use of an ideal uniform array geometry, that is, the antennas have to be uniformly placed with exactly the same separation distance. However, in practical applications, arbitrary arrays arise in several cases, instead of the perfect uniform arrays. For example, the antenna separation distance is measured in the millimeter range over the mmWave frequency bands. Thus, an ideal uniform array geometry is hard to guarantee due to calibration errors introduced in the manufacturing process and/or antenna installation. Another case of arbitrary arrays appears due to sub-array selection. For example, for the purpose of energy saving in the switch-based hybrid architecture, only a small number of antennas is switched to link the RF chains [35]. Arbitrary arrays no longer present the well-featured Vandermonde structure explicitly in the array manifold, which then excludes the use of a large number of geometric-based channel estimation techniques. To overcome this problem, array manifold separation techniques have been developed in array signal processing [36, 37, 38]. However, the structural feature presented via manifold separation is not efficiently utilized in [36, 37], where conventional subspace methods are used based on large samples. In [38], since the manifold separation through the Fourier series approximation is employed to enable ANM, it suffers from an expensive computational complexity in order to reduce the approximation error in the largely expanded Fourier domain.

In this paper, addressing all the aforementioned challenges cohesively, we seek to design high-performance, low-cost channel estimation solutions for arbitrary arrays in hybrid mmWave massive MIMO systems. Specifically, we propose a super-resolution channel estimation framework that not only utilizes the special channel features of the sparse mmWave massive MIMO propagation, but also fully considers the nonideal array geometry and practical hardware limitations. This framework offers several channel estimation solutions and enables to obtain both the channel statistics and the instantaneous CSI, depending on whether the transceiver design is built on channel covariance [21, 22] or the channel itself [39]. We propose two super-resolution solutions for channel covariance estimation (CCE) through efficient structure-based optimization techniques, with samples collected from multiple snapshots. One is the CCE via the low-rank structured covariance reconstruction (LRSCR), which provides super-resolution accuracy at a low computational cost. The other is the CCE via Dynamic-ANM, which further allows for a dynamic configuration where the hybrid hardware parameters change over time for a better performance. For block fading channels, given the estimate from CCE, the instantaneous CSI can then be estimated in a timely fashion. Since the angles change slowly and can take a long time to acquire accurately from the channel statistics while the path gains vary frequently but are easy to acquire given the estimated angles, the instantaneous CSI estimation is divided by solving two subproblems sequentially, i.e., angle estimation and path gain estimation. In developing these novel super-resolution channel estimation approaches, this work contains the following main contributions.

• •

We leverage a generalized array manifold separation approach to extract the useful geometric structure for a practical system with an arbitrary or imperfect array geometry. In particular, we transform the sparse mmWave massive MIMO channel representation from the physical arbitrary antenna domain to a virtual uniform antenna domain via the Jacobi-Anger approximation [40]. Our method enables gridless CS to exploit the useful Vandermonde structure presented in the virtual uniform array manifold.

• •

This work not only develops super-resolution channel estimation solutions, but also investigates the fundamental limits of gridless CS based channel estimation under the constraints of arbitrary arrays and hybrid structures. Our theoretical results shed light on the minimum number of RF chains required by super-resolution channel estimation, as well as the lower and upper bounds on the mode order selected for the Jacobi-Anger approximation. This leads to a tradeoff between the hardware cost of sparse channel estimation and the approximation accuracy to combat the imperfect array geometry.

• •

To reduce the high computational complexity of the SDP-based channel estimation, we design a fast iterative algorithm through the alternating direction method of multipliers (ADMM) [41]. It provides an efficient first-order algorithm implementation with much lower computational complexity than that of the SDP solvers using the interior-point method.

• •

We tackle several practical implementation issues. Specifically, we overcome the side effect of the array manifold separation operation, by removing the spurious peaks generated by the Jacobi-Anger approximation for nonuniform linear arrays. We also extend our work to more complicated 2-dimensional (2D) scenarios, where both the BS and the MS are equipped with multiple antennas.

Simulation results are provided to testify the advantages of the proposed solutions, which make them attractive for hybrid mmWave massive MIMO systems with arbitrary arrays.

The rest of this paper is organized as follows. Section II presents the system model and problem formulation for sparse channel estimation in hybrid arbitrary arrays. Section III proposes a super-resolution channel estimation framework based on the array manifold separation, in which different channel estimation solutions are developed for obtaining the channel covariance and the instantaneous CSI. Specific issues related to the proposed techniques are discussed in Section IV. Simulation results are presented in Section V, followed by conclusions in Section VI.

Notations: $a$ is a scalar, $\bm{a}$ denotes a vector, $\bm{A}$ is a matrix, and $\mathcal{A}$ represents a set. $(\cdot)^{T}$, $(\cdot)^{*}$, and $(\cdot)^{H}$ are the transpose, conjugate, and conjugate transpose of a vector or matrix, respectively. $\text{conv}(\mathcal{A})$ means the convex hull of a set $\mathcal{A}$. $\text{Real}(\cdot)$ and $\text{Imag}(\cdot)$ compute the real part and the imaginary part of a vector or matrix, respectively. $|a|$ denotes the absolute value of $a$. $\left\lVert\bm{a}\right\rVert_{2}$ is the $\ell_{2}$ norm of $\bm{a}$. $\text{diag}(\bm{a})$ and $\text{diag}(\bm{A}_{1},\bm{A}_{2})$ denote a diagonal matrix with the diagonal elements constructed from $\bm{a}$ and a block diagonal matrix with the submatrices $\bm{A}_{1}$ and $\bm{A}_{2}$, respectively. $\bm{I}$ is an identity matrix and $\bm{I}_{a}$ is an anti-diagonal identity matrix. $\text{T}(\bm{u})$ is a Hermitian Toeplitz matrix with first column being $\bm{u}$. $\|\bm{A}\|_{F}$, $\bm{A}^{\dagger}$, and $\text{tr}(\bm{A})$ are the Frobenius norm, the pseudoinverse, and the trace of $\bm{A}$, respectively. The operation $\text{vec}(\cdot)$ stacks all the columns of a matrix into a vector. $\otimes$ is the Kronecker product of matrices or vectors. $\mathbb{E}\{\cdot\}$ denotes expectation.

II Models and Preliminaries

In this section, we first present the signal model and state the goal of both CCE and CSI estimation. Then, we briefly overview the related prior work on relevant super-resolution techniques that are only applicable for ideal uniform arrays, e.g., the uniform linear array (ULA), and under a fixed hybrid hardware structure.

II-A Channel and Signal Models

Consider a narrowband111In a wideband case with frequency selectivity, the continuous-valued delays of the individual paths of the sparse time-dispersive channels can be estimated via gridless CS to achieve super-resolution accuracy in the time domain [42]. mmWave massive MIMO time division duplex (TDD)222This work can be applied to frequency division duplex (FDD) systems as well, given the angle reciprocity between uplink and downlink [25]. system for channel estimation conducted at the base station (BS). As shown in Fig. 1, the BS has a hybrid structure equipped with $N$ arbitrarily deployed antennas and $M\,(M<N)$ RF chains. For simplicity, we mainly focus on the basic single-antenna case at the mobile station (MS), while we extend to the case of multiple-antenna MS as well in Section IV-D. Noticeably, as shown in Fig. 1, to impose the useful Vandermonde structure in an arbitrary array geometry, a preprocessing block via Jacobi-Anger approximation is added to the channel estimator in the hybrid mmWave massive MIMO system, which will be described in Section III-A.

At the mmWave frequency, the wireless channel experiences limited scattering propagation, which results in a sparse multipath structure [14, 15], as shown in Fig. 1. In this sense, the channel can be described by a geometric model with $K\,(K<M<N)$ scatterers333In this work, we focus on the case of point scatterers. The angle spread issue due to the reflecting areas of shaped scatterers is out of scope of this paper. The impact of angle spreads on the proposed methods will be studied in future work., in which each path is parameterized by the path angle and the path gain. For simplicity, suppose that each scatterer contributes to one propagation path, which can be straightforwardly extended to cluster scattering where each cluster includes multiple scattering paths [39]. Further, in mmWave channels, the angles of the scattering paths remain constant for a relatively long time, while the channel coefficients change very rapidly [15]. Accordingly, the uplink channel $\bm{h}_{t}$ can be expressed as the sum of $K$ paths in the form

[TABLE]

where $g_{k,t}$ denotes the channel gain for the $k$-th scattering path at the $t$-th snapshot, and $\bm{a}(\theta_{k})\in\mathbb{C}^{N}$ is the array manifold vector corresponding to the $k$-th channel path.

In this work, we focus on arbitrary arrays, in which the $n$-th antenna element is placed at a known location $(r_{n},\phi_{n}),n=1,\dots,N$ in polar coordinates. Define $\theta_{k}$ as the angle between the polar axis and the $k$-th path, and take the polar origin as the reference point. Then, the $n$-th component of the array manifold vector for the $k$-th path can be written as

[TABLE]

where $\lambda$ denotes the wavelength. In a compact matrix-vector form, the channel $\bm{h}_{t}$ in (1) can be rewritten as

[TABLE]

where $\bm{g}_{t}=[g_{1,t},\dots,g_{K,t}]^{T}$ and $\bm{A}=[\bm{a}(\theta_{1}),\dots,\bm{a}(\theta_{K})]$.

In uplink channel estimation, the MS sends out training symbols $z_{t}$ which are also known to the BS. For simplicity, let $|z_{t}|=1$ for all snapshots. Then, the received signal at the BS’s antennas can be represented as

[TABLE]

where $\bm{w}_{t}$ denotes additive Gaussian noise distributed as $\mathcal{CN}(\bm{0},\sigma^{2}\bm{I})$. According to (4), the covariance matrices for $\bm{x}_{t},\bm{h}_{t}$ and $\bm{g}_{t}$ have the following linear relationship:

[TABLE]

where $\bm{R}_{h}=\mathbb{E}\{\bm{h}_{t}\bm{h}_{t}^{H}\}$, and $\bm{R}_{g}=\mathbb{E}\{\bm{g}_{t}\bm{g}_{t}^{H}\}$.

The BS adopts a hybrid hardware structure in the form of $\bm{W}_{t}=\bm{W}^{\text{BB}}_{t}\bm{W}^{\text{RF}}_{t}$, where $\bm{W}^{\text{BB}}_{t}\in\mathcal{C}^{M\times M}$ denotes a baseband digital combiner, and $\bm{W}^{\text{RF}}_{t}\in\mathcal{C}^{M\times N}$ is an analog combiner. In this paper, we focus on the case where $\bm{W}^{\text{RF}}_{t}$ is made of a network of random phase shifters, while this work can be applied to other structures such as the switch-based network as well [11]. To further enhance randomness, $\bm{W}^{\text{BB}}_{t}$ can be set as a random Gaussian matrix. After being multiplied with the hybrid combining matrix $\bm{W}_{t}\in\mathcal{C}^{M\times N}$ and the known training symbol $z_{t}^{*}$, the received signal at the lower-dimensional baseband is given by

[TABLE]

where $\bm{n}_{t}=z_{t}^{*}\bm{w}_{t}$. When $\bm{W}_{t}$ in (6) is different snapshot by snapshot, it is a dynamic channel sensing system.

In this paper, we assume a block fading channel, where path gains, and hence the CSI, stay constant within a block but vary from block to block. In contrast, path angles vary much slower, and stay unchanged across blocks, until angle re-calibration is needed. Further, angles can be retrieved from the channel covariance, which is the key idea behind the statistical inference methods for angle estimation. This motivates us to design a two-stage channel estimation framework to obtain both the channel covariance and the instantaneous CSI sequentially, as shown in Fig. 2. In the first stage, we apply CCE to obtain the channel covariance over multiple snapshots. Then, in the second stage, considering the difference in time-variation of path angles and path gains, we design a two-step scheme to do instantaneous CSI estimation, which is further divided into two subproblems: angle estimation and path gain estimation.

Remark 1: When the goal of channel estimation is CCE only, the estimator can terminate upon completing the first stage. The CCE by itself is relevant in two cases: one is to simplify either the channel estimation task or the transceiver design, and the other is when path gains experience fast fading that renders the CSI estimates useless for data transmission. In both cases, CCE-based transceiver design can be adopted [21, 22].

Remark 2: Our task of CCE-based angle estimation is also useful during the system calibration stage for fixed wireless applications, in which case the estimated angles can be used to facilitate several system-level tasks such as user grouping and beam-sectoring.

II-B Prior Art: Efficient Super-Resolution Techniques in Ideal Models

When the block length is not long enough for traditional super-resolution methods to work for massive MIMO when $N$ is very large, we have to focus on those super-resolution techniques that exploit the structural feature of the array manifold to improve the sample efficiency. Therefore, in this subsection, we overview existing efficient super-resolution techniques for channel estimation [32, 33], which critically rely on an assumption that the antenna array has to be a uniform array, e.g., a ULA. Then, the array manifold naturally presents an explicit Vandermonde structure. That is, the $n$-th component of the array manifold vector at the $k$-th path is a special case of (2) with $\phi_{n}=0,\,r_{n}=d(n-1),\forall n$, which has the form

[TABLE]

where $d$ denotes the same separation distance between any two adjacent antennas placed along the polar axis. Without loss of generality, suppose $d=\lambda/2$.

Further, assuming $\bm{W}_{t}=\bm{W},\,\forall t$ as the fixed hybrid structure over time, the received signals $\bm{y}_{t}$ for $t=1,\dots,T$ in (6) can be collected to form a matrix $\bm{Y}$ with $\bm{y}_{t}$ being its columns as

[TABLE]

where the matrices $\bm{H}$, $\bm{G}$ and $\bm{N}$ are similarly defined as $\bm{Y}$ with $\bm{h}_{t}$, $\bm{g}_{t}$ and $\bm{n}_{t}$ being their columns, respectively. From (8), the covariance of $\bm{y}_{t}$ is given as

[TABLE]

Then, an atom set in the MMV case is defined as [43, 44]

[TABLE]

where $\bm{a}(f)\in\mathbb{C}^{N}$ with the $n$-th component being $e^{j2\pi(n-1)f}$. According to the atomic norm theorem [43, 44], the atomic norm of $\bm{H}$ over the atom set $\mathcal{A}$ is defined as

[TABLE]

which seeks the most concise representation of $\bm{H}$ by involving the fewest atoms over $\mathcal{A}$.

From the received signals $\bm{Y}$ in (8), the channel estimation for the instantaneous CSI is conducted by solving the regularized ANM formulation as

[TABLE]

where $\tau$ denotes the regularization coefficient controlling the tradeoff between the ANM and the residual error tolerance to the observations. According to [45], $\tau$ can be set as $\tau={1}/{(\sigma+\frac{\sigma}{\log N}{\sqrt{N\log N+N\log(4\pi\log N)}}\,)}$.

Besides the instantaneous CSI itself, the second-order channel statistics in terms of the channel covariance $\bm{R}_{h}$ also play an important role in precoding design for mmWave massive MIMO channels [21, 22]. Suppose the channel gains of the sparse paths are uncorrelated with each other. Then, $\bm{R}_{h}$ not only presents the low rankness due to the channel sparsity, but also is a well-structured positive semidefinite (PSD) Hermitian Toeplitz matrix thanks to the Vandermonde structure of uniform arrays.

To utilize these useful features in the channel statistics, a structure-based optimization approach named low-rank structured matrix reconstruction (LRSMR) [43, 46], can be applied to do CCE, by imposing the Hermitian Toeplitz structure on $\bm{R}_{h}$ as $\bm{R}_{h}=\text{T}(\bm{u}_{h})$. Thus, from the sample covariance matrix $\hat{\bm{R}}_{y}=\frac{1}{T}\bm{Y}\bm{Y}^{H}$, the CCE can be conducted via LRSCR as,

[TABLE]

Unfortunately, the uniform array assumption cannot be guaranteed in practice, considering the antenna misalignment and subarray selection issues arisen in hybrid mmWave massive MIMO systems. Moreover, the assumption of the fixed combining matrix over all snapshots is ineffective to find out all potential scattering paths in mmWave channel estimation. Regarding these practical situations, two questions arise: 1) can we design new super-resolution channel estimation techniques via LRSCR and ANM for arbitrary arrays? 2) how can we retrieve the desired channel information via the dynamic hybrid structure? In the next section, to fully address these problems, we develop super-resolution and fast channel estimation approaches for obtaining both the channel covariance and the instantaneous CSI.

III Super-Resolution Channel Estimation for Arbitrary Arrays

In this section, we first leverage the Jacobi-Anger approximation to extract the useful Vandermonde structure from a nonideal array geometry. Then, two new CCE methods are developed through structure-based optimization techniques for arbitrary arrays. For efficient estimation of the instantaneous CSI, a two-step solution is proposed, which estimates the path angles and path gains sequentially.

III-A Jacobi-Anger Approximation

The Jacobi-Anger expansion provides a general infinite terms expansion of exponentials of trigonometric functions in the basis of their harmonics [40]. Specifically, (2) in the Jacobi-Anger expansion form is expressed as

[TABLE]

where $\text{J}_{i}(\cdot)$ denotes the $i$-th order Bessel function of the first kind.

Although (14) indicates a summation of infinite series, the value of $|\text{J}_{i}(r)|$ decays very rapidly as the value of $|i|$ increases for any $r>0$, which is a nature of the Bessel function. Thus, the infinite series expansion can be well approximated by keeping only the terms having large absolute values located around the central range of the series, that is, $|i|\leq I$ in (14). To reach a desired precision, the maximum mode order $I$ for the approximation is chosen as [47]

[TABLE]

where $r_{max}$ is the maximum $r_{n}$.

Then, given $I$, (14) can be approximately expressed as

[TABLE]

where $\bm{c}_{n}$ and $\bm{v}(\theta_{k})$ are given by

[TABLE]

According to (16), the $N\times 1$ array manifold vector of an arbitrary array can be approximated as

[TABLE]

Noticeably, it is the Jacobi-Anger expansion that enables to separate the unknown channel-related parameter ($\theta_{k}$) in (17b) from the known array-related configurations ($r_{n}$ and $\phi_{n}$) in (17a), which are otherwise mingled in the original physical antenna domain of (2). Further, due to the exponential form in (17b), the Vandermonde structure is well-presented in the virtual uniform antenna domain, in terms of the virtual array manifold

[TABLE]

Merging (18) and (19), any arbitrary array can be expressed as the multiplication of a Bessel matrix and a Vandermonde matrix, in the form

[TABLE]

Therefore, to appreciate the Vandermonde structure, we reformulate the channel by substituting (18) into (3) as

[TABLE]

where $\bm{d}_{t}=\bm{V}\bm{g}_{t}$ represents the virtual channel in the virtual uniform antenna domain. From (5) and (20), the channel covariance can be rewritten as

[TABLE]

where $\bm{R}_{v}$ denotes the virtual channel covariance matrix given by

[TABLE]

Accordingly, after taking the hybrid structure and arbitrary array into consideration, the received signal in (6) can be expressed as

[TABLE]

where $\bm{\Phi}_{t}=\bm{W}_{t}\bm{C}$ denotes the equivalent sensing matrix.

As the Bessel matrix $\bm{C}$ depends solely on the known array geometry, $\bm{R}_{h}$ in (22) can be estimated as long as $\bm{R}_{v}$ is retrieved. Next, we need to figure out how to estimate $\bm{R}_{v}$ from collected $\{\bm{y}_{t}\}_{t}$.

III-B Channel Covariance Estimation

In this subsection, to estimate the virtual channel covariance $\bm{R}_{v}$, and hence also the actual channel covariance $\bm{R}_{h}$, we develop two super-resolution CCE techniques for the MMV case with multiple snapshots. One is the CCE via LRSCR, and the other is the CCE via ANM. The ANM-based CCE method is applicable to the dynamic hybrid hardware structure.

III-B1 CCE via LRSCR

Suppose the path gains of the fading channels are uncorrelated with each other. Then, $\bm{R}_{v}$ in (23) not only presents the low rankness because of sparse scattering environments, but also is a well-structured PSD Hermitian Toeplitz matrix due to the Vandermonde structure of the virtual uniform array. In this sense, the LRSCR technique can be used to implement the low-rank feature of $\bm{R}_{v}$. Moreover, suppose the hybrid combining matrix $\bm{W}_{t}$ is fixed over time, i.e., $\bm{W}_{t}=\bm{W}$, and $\bm{\Phi}_{t}=\bm{WC}=\bm{\Phi},\forall t$. Substituting (22) into (9), $\bm{R}_{y}$ can then be rewritten as

[TABLE]

which clearly is a function of $\bm{R}_{v}$.

Accordingly, by imposing the Hermitian Toeplitz structure on $\bm{R}_{v}$ in terms of $\bm{R}_{v}=\text{T}(\bm{u}_{v})$, we can describe the LRSCR-based formulation for the estimation of $\bm{R}_{v}$ from the sample covariance $\hat{\bm{R}}_{y}$ as

[TABLE]

Then, given $\hat{\bm{R}}_{v}$ estimated from (26), we finally obtain $\hat{\bm{R}}_{h}$ via (22), with known $\bm{C}$.

Further, to solve (26) with lower computational complexity compared with using off-the-shelf SDP solvers [48], we will develop a fast algorithm via ADMM later in Section IV-B.

It is worth noting that CCE via LRSCR can be done in blind mode from data symbols, since it can work as long as $\hat{\bm{R}}_{y}$ is available to (26).

III-B2 CCE via Dynamic-ANM

Let $\bm{D}=[\bm{d}_{1},\dots,\bm{d}_{T}]$ collect the virtual channels defined in (21) from different time slots. Then, an atom set can be defined in the virtual antenna domain as

[TABLE]

where $\bm{v}(f)\in\mathbb{C}^{2I+1}$ with its components being $e^{j2\pi if},\,i=-I,\ldots,I$. Obviously, $\bm{D}$ is a linear combination of the atoms from the set $\mathcal{A}^{\prime}$. In this sense, the ANM technique can be used to exploit the low rank and Vandermonde features of $\bm{D}$.

Then, with the received signals $\{\bm{y}_{t}\}_{t}$ in the matrix-form $\bm{Y}$, we produce the following ANM formulation:

[TABLE]

Next, we consider a time-varying $\bm{\Phi}_{t}$ in the dynamic hybrid structure case. Since the ANM in our objective function can be maintained, we now only have to formulate the residual error tolerance snapshot by snapshot. Hence, the dynamic version of (28) can be rewritten as

[TABLE]

Further, replacing the atomic norm in (29) by its SDP formulation [43][44], we have

[TABLE]

For solving the problem (30), an ADMM-based fast algorithm can also be designed and iteratively implemented in a similar way as for (26).

In addition, according to the atomic norm definition and its SDP formulation, the reconstructed $\text{T}(\hat{\bm{u}}_{d})$ from (30) can be expressed by a Vandermonde decomposition as

[TABLE]

Finally, using the estimated $\hat{\bm{V}}$ and $\hat{\bm{R}}_{g}$, we can obtain $\hat{\bm{R}}_{v}$ via (23) and then $\hat{\bm{R}}_{h}$ via (22) accordingly.

III-C Two-Step Instantaneous CSI Estimation

As has been mentioned in Section II-A, the spatially sparse mmWave channel $\bm{h}$ is fully determined by parameters: $\bm{\theta}$ and $\bm{g}$. Since the path angles depend only on the relative positions of the BS, the MS, and the scatterers, $\bm{\theta}$ varies slowly, but can take a long time to acquire accurately from $\hat{\bm{R}}_{h}$ according to statistical inference techniques for angle estimation [27, 28, 49]. In contrast, the path gains $\bm{g}$ are easy to acquire given $\hat{\bm{\theta}}$, but vary frequently. Therefore, given the estimated $\hat{\bm{R}}_{h}$, in the second stage, we design a two-step CSI estimation scheme for block transmission. Angle estimates directly result from the CCE output, which stay unchanged for multiple blocks until they change. Given $\hat{\bm{\theta}}$, path gains are estimated at each block, followed by CSI-based data transmission.

III-C1 Angle Estimation

Thanks to the specific matrix structures presented by the second-order statistics of mmWave channels, such as the low rankness due to the sparse scattering propagation and the PSD Toeplitz structure imposed by the SDP formulation through either (26) or (30), the Vandermonde decomposition can be applied to estimate the angles.

• •

If the LRSCR technique is used for CCE in Section III-B1, the recovered $\hat{\bm{R}}_{v}$ from (26) is a low-rank PSD Toeplitz matrix. According to the Vandermonde decomposition lemma [50], $\hat{\bm{R}}_{v}$ can be uniquely expressed as

[TABLE]

This Vandermonde decomposition can be computed efficiently via root finding or by solving a generalized eigenvalue problem [49]. Since the virtual uniform array geometry is solely decided by the angles as in (19), $\hat{\bm{\theta}}$ can be directly extracted from $\hat{\bm{V}}$ according to (17b).

• •

For the CCE based on the Dynamic-ANM technique as developed in Section III-B2, the Vandermonde decomposition can be carried out as in (31). Although the definition of the atom set in (27) leads to a diagonal matrix in the form of $(T\hat{\bm{R}}_{g})^{\frac{1}{2}}$ in (31) that is different from $\hat{\bm{R}}_{g}$ in (32), the common Vandermonde structure of $\hat{\bm{V}}$ still leads to the same estimation results for $\hat{\bm{\theta}}$ via (17b).

III-C2 Path Gain Estimation

Given the obtained angular information $\hat{\bm{\theta}}$ from Section III-C1, next we need to estimate the path gains $\bm{g}$ in a timely fashion as shown in Fig. 2. To this end, we first form the array matrix $\bm{A}$ via (2). Then, following the principle of a matched filter, we tune the precoder as $\bm{W}=\bm{A}^{H}$ for beamforming. Noteworthily, with the obtained $\hat{\bm{\theta}}$ from angle estimation, the $K\times N$ matched filter based $\bm{W}$ applied now for path gain estimation is different from the $M\times N$ random phase shifter based $\bm{W}$ used earlier for CCE when $\hat{\bm{\theta}}$ is unknown. As a result, the estimation of $\bm{g}$ can be expressed as a least squares formulation

[TABLE]

IV Discussion of Related Issues

In this section, we provide detailed discussions on some specific issues related to the proposed solutions. We first provide the theoretical results in terms of fundamental limits for the proposed super-resolution channel estimation in hybrid mmWave massive MIMO with arbitrary arrays. Then, we design a first-order algorithm via ADMM to rapidly implement the super-resolution estimation in lieu of invoking the high-computational SDP. In addition, we study the spurious-peak issue of the Jacobi-Anger approximation as a side effect specific to nonuniform linear arrays, and provide an effective way to solve this problem. Finally, we extend the work to the multiple-antenna MS case, which is developed based on an efficient 2D gridless CS approach.

IV-A Analysis of Fundamental Limits

To get an approximation with certain precision, (15) offers a lower bound on the choice of the maximum mode order $I$ of the Jacobi-Anger expansion. From the view of the Jacobi-Anger approximation, the larger $I$ is, the higher the accuracy the approximation can achieve. On the other hand, from the view of sparse channel estimation, we need to retrieve the sparse virtual channel $\{\bm{d}_{t}\}_{t}$ from the compressed measurements $\{\bm{y}_{t}\}_{t}$ with a high probability. For a given hybrid mmWave massive MIMO system, the dimension of $\bm{y}_{t}$ is fixed and known as $M$. To ensure the proposed gridless CS based methods are feasible and effective, first, it is necessary to study the minimum number of RF chains $M_{min}$ required by our approaches given the approximation and hybrid structure. According to Theorem III.4 in [51], to guarantee successful reconstruction of channel covariance with high probability, $M_{min}$ can be expressed as a function of $I$ and $K$ as

[TABLE]

where $C$ is a numerical constant. Note that although the sensing matrix $\bm{\Phi}_{t}$ in [51] is assumed to be an i.i.d. random Gaussian matrix, it is reasonable to relax (34) to accommodate the analysis in this paper when $\bm{\Phi}_{t}$ is modeled as a random phase shifter based $\bm{W}_{t}$ multiplied by a Bessel matrix $\bm{C}$. Then, replacing $M_{min}$ in (34) by $M$ with $M>M_{min}$, we can obtain an upper bound on $I$. Combining with (15), we have

[TABLE]

It is worth noting that (35) actually reflects the tradeoff between the approximation accuracy to combat the imperfection of the array geometry and the hardware cost required for sparse channel estimation, which thus sheds light on the choices of $I$ and $M$ in practice.

IV-B Fast Implementation via ADMM

To avoid the high computational complexity of the SDP-based solutions, we develop a fast iterative algorithm via ADMM. Next, we mainly discuss the solution for the LRSCR formulation (26) and omit that for the ANM case (30)444The design of an ADMM-based fast algorithm for (30) can be developed in a similar way as for (26). The differences are the two additional variables $\bm{D}$ and $\bm{Q}$ in (30), which can be easily updated by gradient descent in each iteration. Moreover, the algorithm implementation for (30) is actually simpler than that for (26). The reason is that the Toeplitz structured matrix $\text{T}(\bm{u})$ is not included in the least squares term in the objective function in (30), which then simplifies the calculation of the partial derivative of the augmented Lagrangian $\mathcal{L}$ with respect to $\bm{u}^{*}$. As a result, compared with the LRSCR case, the update of $\bm{u}$ in each iteration becomes easier in the ANM case.. To apply ADMM [41], we reformulate (26) as

[TABLE]

whose augmented Lagrangian can be expressed as

[TABLE]

where $\bm{U}$ and $\bm{\Lambda}$ are Hermitian matrices. Then the implementation of ADMM involves the following iterative updates:

[TABLE]

where the superscript $l$ denotes the $l$-th iteration update. In order to implement (38), we take the partial derivative of (37) with respect to $\bm{u}_{v}^{*}$ at the $(l\!+\!1)$-th iteration and force it equal to zero. After taking a series of derivations on $\frac{\partial}{\partial\bm{u}_{v}^{*}}\mathcal{L}(\bm{u}_{v},\bm{U}^{l},\bm{\Lambda}^{l})\left|{}_{\bm{u}_{v}=\bm{u}_{v}^{l+1}}\right.=0$, we obtain

[TABLE]

where $N_{I}=2I+1$ is the column (row) size of $\text{T}(\bm{u}_{v}^{l+1})$, $\bm{e}_{1}$ is the $N_{I}$-length vector with only the first element being one, and $\bm{b}=\mathcal{G}(\bm{B})\in\mathbb{C}^{N_{I}}$ is a mapping from a matrix to a vector where the $n_{I}$-th element of $\bm{b}$ is the sum of all the elements $\bm{B}_{i,j}$ in $\bm{B}$ satisfying $i-j+1=n_{I}$. Moreover, denote $\bm{M}$ as the matrix which satisfies $\bm{b}=\mathcal{G}(\bm{B})=\bm{M}\text{vec}(\bm{B})$ and $\bm{\beta}^{l}=\tau\mathcal{G}(\bm{\Phi}^{H}\hat{\bm{R}}_{y}\bm{\Phi})+\rho\mathcal{G}(\bm{U}^{l}+\rho^{-1}\bm{\Lambda}^{l})$, respectively.

Accordingly, we rewrite (LABEL:eq:partialdiff) as

[TABLE]

where $\bm{u}_{v}^{l+1}=\bm{u}_{\text{R}}^{l+1}+j[0,({\bm{u}_{\text{I}}^{l+1}})^{T}]^{T}$. Since $\text{T}(\bm{u}_{v}^{l+1})$ is only determined by the real and imaginary parts of $\bm{u}_{v}^{l+1}$ as $\bm{u}_{\text{R}}^{l+1}$ and $\bm{u}_{\text{I}}^{l+1}$, respectively, there exists a fixed matrix $\bm{\Gamma}$ satisfying $\text{vec}(\text{T}(\bm{u}_{v}^{l+1}))=\bm{\Gamma}[(\bm{u}_{\text{R}}^{l+1})^{T},(\bm{u}_{\text{I}}^{l+1})^{T}]^{T}$. Moreover, $\bm{\Pi}=\left(\tau\bm{M}\left[(\bm{\Phi}^{H}\bm{\Phi})^{T}\otimes(\bm{\Phi}^{H}\bm{\Phi})\right]+\rho\bm{M}\right)\bm{\Gamma}$. Since $[(\bm{u}_{\text{R}}^{l+1})^{T},(\bm{u}_{\text{I}}^{l+1})^{T}]^{T}\in\mathbb{R}^{2N_{I}-1}$, we can rewrite the $N_{I}$ complex equations of (LABEL:eq:rewrite-partialdiff) into $2N_{I}$ real equations as

[TABLE]

Hence, the update rule for $\bm{u}_{v}$ is given by

[TABLE]

Let $\Xi^{l}=\text{T}(\bm{u}_{v}^{l+1})\!-\!\rho^{-1}\bm{\Lambda}^{l}=\bm{E}^{l}\bm{\Sigma}^{l}{\bm{E}^{l}}^{H}$ be its eigenvalue decomposition, then based on (37) and (39), we have the update of $\bm{U}$ at the $(l\!+\!1)$-th iteration as

[TABLE]

where $\bm{\Sigma}_{+}^{l}$ is obtained by letting all negative eigenvalues of $\bm{\Sigma}^{l}$ be zero.

The iterative algorithm will stop until both primal and dual residuals satisfy the pre-set tolerance level [41].

IV-C Specific Instance of Nonuniform Linear Arrays

Nonuniform linear arrays yield one case of arbitrary arrays. However, unlike general random distributed antennas, the antenna elements of nonuniform linear arrays are distributed along a line, which results in the Bessel matrix used by the Jacobi-Anger approximation being axial symmetric. This axial symmetric characteristic leads to a special issue in the implementation of the proposed Dynamic-ANM and LRSCR methods. Next, we discuss this specific instance in detail.

Suppose the locations of antenna elements of a nonuniform linear array are formed as $(r_{n},0)$, $n=1,\dots,N$ in polar coordinates. Then, (16) is rewritten as

[TABLE]

Based on the property of the Bessel function of the first kind that says

[TABLE]

we have

[TABLE]

From (48), the Bessel matrix $\bm{C}$ explicitly holds

[TABLE]

Accordingly, (18) can be expressed as

[TABLE]

Then, given (50), (21) can be rewritten as

[TABLE]

which indicates

[TABLE]

This means that (52) turns out be an alternative possible solution of (30). Accordingly, different from (31), the estimated $\text{T}(\bm{u}_{d})$ leads to another Vandermonde decomposition as

[TABLE]

Further, since $\bm{R}_{g}$ is a diagonal matrix with positive elements, substituting (50) into (22), we have

[TABLE]

Noticeably, instead of (32), $\text{Real}(\bm{R}_{v})$ is decomposed as

[TABLE]

As a result, beside the estimation of the true AoAs as $\{\hat{\theta}_{k}\}$, spurious results are also generated as $\{-\hat{\theta}_{k}\}$. According to the above analysis that there exist multiple solutions, we need to use the prior knowledge that the feasible domain of AoAs is $[0,\pi)$. Then we can simply remove the spurious results appearing in $(-\pi,0)$.

IV-D Extension to the Multiple-Antenna MS Case

In this subsection, we extend the work to the case where the MS also has a hybrid architecture with multiple RF chains and antennas in arbitrary arrays. Now, the 1D uplink channel model for the case of the single-antenna at the MS in (1) is extended to a 2D uplink channel model

[TABLE]

where $\theta_{\text{MS},k}$ and $\theta_{\text{BS},k}$ denote the continuously-valued AoD and AoA of the $k$-th path at the MS as transmitter and at the BS as receiver, respectively. Accordingly, the vectors $\bm{a}_{\text{MS}}\left(\theta_{\text{MS},k}\right)$ and $\bm{a}_{\text{BS}}\left(\theta_{\text{BS},k}\right)$ represent the array manifold vectors corresponding to the $k$-th path for the $N_{\text{MS}}$-antenna and $N_{\text{BS}}$-antenna arrays, respectively, which both have components in the form of (2). This 2D channel model is general enough to subsume the multi-user case where each column of $\bm{H}$ corresponds to one MS (user) with a single antenna.

Given the hybrid structures and the arbitrary arrays employed at both the MS and BS sides, the received signal can be expressed as

[TABLE]

where $\bm{F}$ denotes the hybrid precoding matrix used at the transmitter side.

Applying the array manifold separation approach described in (20), (57) can be rewritten as

[TABLE]

where $\bm{\Psi}=\bm{V}_{\text{BS}}\,\mbox{diag}(\bm{g})\,\bm{V}_{\text{MS}}^{H}$ denotes the virtual 2D channel that presents the 2D Vandermonde structure in the virtual uniform antenna domain.

To estimate the 2D channel through gridless CS, a straightforward way is to vectorize the 2D formulations and then to cast the 2D Vandermonde structure into a vectorized SDP formulation via a two-level Toeplitz structured matrix, a.k.a., vectorization based ANM (V-ANM) [32, 52, 53]. However, the V-ANM leads to a high computational complexity on the order of $\mathcal{O}(N_{\text{BS}}^{3.5}N_{\text{MS}}^{3.5})$ [54], because of the huge problem scale resulting from the vectorization operation.

To solve this problem of V-ANM, we develop an efficient 2D channel estimation at much lower computational cost, by using a decoupled-ANM (D-ANM) technique [55, 54]. Different from the V-ANM, we introduce a matrix-form atom set $\mathcal{A}_{\text{M}}$ as

[TABLE]

which naturally results in a matrix-form atomic norm as

[TABLE]

Then, given $\bm{Y}$ from (58), the virtual 2D channel $\bm{\Psi}$ can be reconstructed via the following decoupled SDP formulation

[TABLE]

Noticeably, since the PSD constraint in (61) is of size $(N_{\text{BS}}+N_{\text{MS}})\times(N_{\text{BS}}+N_{\text{MS}})$, the D-ANM allows a reduced computational complexity on the order of $\mathcal{O}\left((N_{\text{BS}}+N_{\text{MS}})^{3.5}\right)$ [54], which is much smaller than that of the V-ANM with large arrays.

V Numerical Results

This section presents numerical results to evaluate the channel estimation performance achieved by the proposed methods for arbitrary arrays and a hybrid precoding structure. In each Monte Carlo simulation, the random path angles are generated uniformly from $[0^{\circ},180^{\circ})$. The existing channel covariance estimation methods via covariance orthogonal matching pursuit (COMP) and Dynamic-COMP (DCOMP) and the existing instantaneous CSI estimation methods via simultaneous orthogonal matching pursuit (SOMP) and Dynamic-SOMP (DSOMP) are also simulated as benchmarks for performance comparison [20], where a predefined grid of size $360$ is employed for the grid-based CS technique that leads to an angle-resolution of $0.5^{\circ}$.

V-A Channel Estimation Performance

First, we testify the performance of different channel estimation approaches with an arbitrary planar array, in terms of the normalized mean squared error (NMSE) for the CCE as $\mathbb{E}\{{\|\bm{R}_{h}-\hat{\bm{R}}_{h}\|_{F}^{2}}\}/\mathbb{E}\{{\|\bm{R}_{h}\|_{F}^{2}}\}$ and the NMSE for the CSI estimation as $\mathbb{E}\{{\|\bm{h}-\hat{\bm{h}}\|_{2}^{2}}\}/\mathbb{E}\{{\|\bm{h}\|_{2}^{2}}\}$, respectively. In simulations, our LRSCR-based methods as described in Section III-B1 and the COMP-based and SOMP-based methods in [20] are tested on the fixed combining matrix $\bm{W}$, while our Dynamic-ANM-based methods as developed in Section III-B2 and the DCOMP-based and the DSOMP-based methods in [20] are applied with the dynamic combining matrix $\bm{W}_{t}$. In our proposed two-step CSI estimation scheme, we employ a Vandermonde decomposition in the form of either (31) or (32) to retrieve the path angles in the first step and to estimate the path gains via least squares in (33) in the second step.

Fig. 3 and Fig. 4 present the NMSE of CCE and the NMSE of CSI versus the number of snapshots, respectively. The comparison of the curves indicates that our proposed methods based on LRSCR and Dynamic-ANM outperform the existing methods based on grid-based CS. While the dynamic configurations usually provide a higher sensing accuracy than the fixed counterparts, our LRSCR can even work better than DSOMP as the number of snapshots increases. This is because our methods utilize not only the sparsity of mmWave channels but also the structural feature of the array geometry. Further, our Dynamic-ANM method always achieves the best performance especially given a small number of snapshots, because ANM can efficiently utilize such structures directly from the collected samples. On the other hand, since LRSCR is a statistics-based design, it requires a sufficient number of snapshots for computing an accurate sample covariance. When the number of snapshots becomes small e.g. less than $20$, the finite-sample effect ruins the Toeplitz structure presented in the ideal covariance matrix, which thus degrades the performance of LRSCR. Fig. 5 and Fig. 6 present the NMSE performance of these approaches for different numbers of RF chains, which show the same trends as in Fig. 3 and Fig. 4.

V-B Nonuniform Linear Array

Next, we study the nonuniform linear array as a special case of arbitrary arrays. Fig. 7 shows that the NMSE of the instantaneous CSI estimation via our Dynamic-ANM method can be much smaller than that of the existing ANM without array manifold separation [32], which demonstrates the necessity of the Jacobi-Anger approximation for an imperfect array geometry. Besides, it also indicates that our method proposed for arbitrary arrays can approach the performance of the ideal case of a same-size ULA as the benchmark for the best performance that the proposed techniques can achieve.

Further, we study the side effect of the Jacobi-Anger approximation in terms of the spurious peaks generated in nonuniform linear arrays. Fig. 8 shows the spatial spectra result of LRSCR, where the peaks on the right indicate the true angles. Meanwhile, the spurious peaks appear at the symmetric angles, which can be simply removed given the prior knowledge of the angle range as $[0^{\circ},180^{\circ})$.

V-C Computational Complexity

In addition, we test the computational cost of the proposed fast algorithm via ADMM, compared with the SDP-based solver [56]. By counting the runtime versus the number of antennas, Fig. 9 clearly indicates that as the number of antennas increases the slope of the runtime curve of the ADMM-based solution is much smaller than that of the SDP counterpart. Thus, the proposed fast implementation has low computational complexity and is well suited for large arrays.

V-D 2D Path Angle Estimation

Last but not least, we extend our proposed work to the multiple-antenna MS case where both the BS and MS are equipped with arbitrary multiple antennas, which results in a 2D path angle estimation scenario. As shown in Fig. 10, our proposed D-ANM based on the Jacobi-Anger approximation can precisely retrieve both the AoAs and the AoDs of the sparse scattering paths, which indicates the high performance of our proposed super-resolution 2D channel estimation techniques.

VI Conclusions

Recognizing the imperfect geometry of arbitrary arrays and the hardware constraint of hybrid mmWave massive MIMO systems, this paper has proposed a new super-resolution channel estimation framework that achieves the benefits of the array manifold separation techniques and the structure-based optimization approaches. Through the Jacobi-Anger approximation, the Vandermonde structure is effected in the virtual antenna domain of arbitrary arrays, which enables super-resolution channel estimation based on gridless CS techniques to obtain a high performance at low training costs. In particular, we develop two channel covariance estimation approaches via LRSCR and Dynamic-ANM. Further, considering that angles change relatively slower than path gains, we design a two-step CSI estimation scheme which separates long-term angle estimation from frequent path gain estimation. The theoretical results are provided to investigate the fundamental limits of the proposed super-solution technique in terms of the minimum number of RF chains required for channel estimation and the bounds on the mode order for the Jacobi-Anger approximation. To reduce the computational complexity of the structure-based optimization via SDP, a first-order iterative algorithm is developed through ADMM for fast implementation. To combat the side effect due to the Jacobi-Anger approximation occurring in nonuniform linear arrays, we provide a mechanism to efficiently remove the spurious peaks. Finally, we extend our work to the 2D-angle scenarios, where both the BS and MS are equipped with multiple antennas.