Channel Covariance Matrix Estimation via Dimension Reduction for Hybrid MIMO MmWave Communication Systems

TL;DR

This paper proposes a low-complexity method for estimating channel covariance matrices in hybrid MIMO mmWave systems, addressing the challenge of lower-dimensional observations due to hybrid structures.

Contribution

It formulates the covariance estimation as a structured low-rank matrix sensing problem and introduces an efficient algorithm to solve it.

Findings

Effective covariance estimation demonstrated with ULA and USPA arrays.

Proposed method outperforms traditional approaches in accuracy and complexity.

Numerical results validate the approach's practicality in mmWave MIMO systems.

Abstract

Hybrid massive MIMO structures with lower hardware complexity and power consumption have been considered as a potential candidate for millimeter wave (mmWave) communications. Channel covariance information can be used for designing transmitter precoders, receiver combiners, channel estimators, etc. However, hybrid structures allow only a lower-dimensional signal to be observed, which adds difficulties for channel covariance matrix estimation. In this paper, we formulate the channel covariance estimation as a structured low-rank matrix sensing problem via Kronecker product expansion and use a low-complexity algorithm to solve this problem. Numerical results with uniform linear arrays (ULA) and uniform squared planar arrays (USPA) are provided to demonstrate the effectiveness of our proposed method.

Tables

Algorithm 1: the GCG-Alt Algorithm for Estimating
$\widehat{𝐂}$ of (28)

Equations

$$\mathbf{H}=\frac{1}{\sqrt{L}}\sum _{\_}{k=1}^K\sum _{\_}{l=1}^L{g}_{\_}kl{\mathbf{a}}_{\_}r({\varphi }_{\_}^{r}kl,{\theta }_{\_}^{r}kl){\mathbf{a}}_{\_}{t}^H({\varphi }_{\_}^{t}kl,{\theta }_{\_}^{t}kl)\in {\mathbb{C}}^{N_{\_}r\times N_{\_}t},$$

$$\mathbf{R}$$

$$=\frac{1}{L}\sum _{\_}{k=1}^K{\gamma }_{\_}^{2}k\sum _{\_}{l=1}^L{\tilde{\mathbf{T}}}_{\_}^{t}kl\otimes {\tilde{\mathbf{T}}}_{\_}^{r}kl\in {\mathbb{C}}^{N_{\_}r{N}_{\_}t\times N_{\_}r{N}_{\_}t},$$

$${\tilde{\mathbf{T}}}_{\_}^{t}kl\triangleq {\mathbf{a}}_{\_}^{\ast }t({\varphi }_{\_}^{t}kl,{\theta }_{\_}^{t}kl){\mathbf{a}}_{\_}^{T}t({\varphi }_{\_}^{t}kl,{\theta }_{\_}^{t}kl)\in {\mathbb{C}}^{N_{\_}t\times N_{\_}t},$$

$${\tilde{\mathbf{T}}}_{\_}^{r}kl\triangleq {\mathbf{a}}_{\_}r({\varphi }_{\_}^{r}kl,{\theta }_{\_}^{r}kl){\mathbf{a}}_{\_}^{H}r({\varphi }_{\_}^{r}kl,{\theta }_{\_}^{r}kl)\in {\mathbb{C}}^{N_{\_}r\times N_{\_}r}.$$

$$\mathbf{a}({\varphi }_{\_}kl)=\frac{1}{\sqrt{N_{\_}a}}[1,{\mathrm{e}}^{j\frac{2\pi }{{\lambda }_{\_}c}d\sin ({\varphi }_{\_}kl)},\cdots ,{\mathrm{e}}^{j(N_{\_}a-1)\frac{2\pi }{{\lambda }_{\_}c}d\sin ({\varphi }_{\_}kl)}{]}^T,$$

$${\tilde{\mathbf{T}}}_{\_}^{t}kl={\mathbf{a}}_{\_}^{\ast }t({\varphi }_{\_}^{t}kl){\mathbf{a}}_{\_}^{T}t({\varphi }_{\_}^{t}kl)$$

$${\tilde{\mathbf{T}}}_{\_}^{r}kl={\mathbf{a}}_{\_}r({\varphi }_{\_}^{t}kl){\mathbf{a}}_{\_}^{H}r({\varphi }_{\_}^{t}kl),$$

$$\mathbf{y}={\mathbf{W}}^H\mathbf{H}\mathbf{f}s+{\mathbf{W}}^H\mathbf{n},$$

$${\mathbf{T}}_{\_}^{t}kl=\frac{{\gamma }_{\_}k}{\sqrt{L}}{\tilde{\mathbf{T}}}_{\_}^{t}kl\in {\mathbb{C}}^{N_{\_}t\times N_{\_}t}$$

$${\mathbf{T}}_{\_}^{r}kl=\frac{{\gamma }_{\_}k}{\sqrt{L}}{\tilde{\mathbf{T}}}_{\_}^{r}kl\in {\mathbb{C}}^{N_{\_}r\times N_{\_}r}$$

$$\mathbf{R}=\sum _{\_}{k=1}^K\sum _{\_}{l=1}^L{\mathbf{T}}_{\_}^{t}kl\otimes {\mathbf{T}}_{\_}^{r}kl\in {\mathbb{C}}^{N_{\_}t{N}_{\_}r\times N_{\_}t{N}_{\_}r},$$

$${\mathbf{R}}_{\_}mn\triangleq [\mathbf{R}{]}_{\_}((m-1)N_{\_}r+1):m{N}_{\_}r,((n-1)N_{\_}r+1):n{N}_{\_}r,$$

$${\mathbf{R}}_{\_}p=\mathcal{P}(\mathbf{R})$$

$$[\mathcal{P}(\mathbf{R}){]}_{\_}m+(n-1)N_{\_}t,:={\text{vec}}^T({\mathbf{R}}_{\_}mn)\in {\mathbb{C}}^{1\times N_{\_}^{2}r}.$$

$${\mathbf{t}}_{\_}^{t}kl=\text{vec}({\mathbf{T}}_{\_}^{t}kl)\in {\mathbb{C}}^{N_{\_}^{2}t\times 1},$$

$${\mathbf{t}}_{\_}^{r}kl=\text{vec}({\mathbf{T}}_{\_}^{r}kl)\in {\mathbb{C}}^{N_{\_}^{2}r\times 1}.$$

$${\mathbf{R}}_{\_}p=\sum _{\_}{k=1}^K\sum _{\_}{l=1}^L{\mathbf{t}}_{\_}^{t}kl({\mathbf{t}}_{\_}^{r}kl{)}^T\in {\mathbb{C}}^{N_{\_}{t}^2\times N_{\_}{r}^2}.$$

$${\mathbf{t}}_{\_}^{t}kl$$

$${\mathrm{e}}^{-j(N_{\_}t-1)\frac{2\pi }{{\lambda }_{\_}c}d\sin ({\varphi }_{\_}^{t}kl)}{\mathbf{a}}_{\_}^{T}t({\varphi }_{\_}^{t}kl){]}^T$$

$$=\frac{{\gamma }_{\_}k}{\sqrt{L}}{\mathbf{a}}_{\_}^{\ast }t({\varphi }_{\_}^{t}kl)\otimes {\mathbf{a}}_{\_}t({\varphi }_{\_}^{t}kl),$$

$${\mathcal{C}}_{\_}t=\{{\mathbf{a}}_{\_}^{\ast }t({\varphi }_{\_}^{t}kl)\otimes {\mathbf{a}}_{\_}t({\varphi }_{\_}^{t}kl),1\le k\le K,1\le l\le L\}.$$

$$p_{\_}e\stackrel{\Delta }{=}\frac{{\sum }_{\_}^{r_{\_}\mathrm{sub}}j=1{\sigma }_{\_}^{2}j}{{\sum }_{\_}^{r_{\_}\mathrm{ch}}i=1{\sigma }_{\_}^{2}i}$$

$${\varphi }_{\_}^{t}kl\sim \mathcal{U}({\varphi }_{\_}^{t}k-{\upsilon }_{\_}^{t}h,{\varphi }_{\_}^{t}k+{\upsilon }_{\_}^{t}h),l=1,2,\cdots ,L,$$

$${\varphi }_{\_}^{r}kl\sim \mathcal{U}({\varphi }_{\_}^{r}k-{\upsilon }_{\_}^{r}h,{\varphi }_{\_}^{r}k+{\upsilon }_{\_}^{r}h),l=1,2,\cdots ,L,$$

$${\mathbf{a}}_{\_}t({\varphi }_{\_}^{t}kl)=\frac{1}{\sqrt{3}}[1,{\mathrm{e}}^{j\frac{2\pi }{{\lambda }_{\_}c}d\sin ({\varphi }_{\_}^{k}kl)},{\mathrm{e}}^{j2\frac{2\pi }{{\lambda }_{\_}c}d\sin ({\varphi }_{\_}^{k}kl)}{]}^T.$$

$${\mathbf{t}}_{\_}^{t}kl=\frac{\sqrt{L}}{{\gamma }_{\_}k}{\mathbf{a}}_{\_}^{\ast }t({\varphi }_{\_}^{t}kl)\otimes {\mathbf{a}}_{\_}t({\varphi }_{\_}^{t}kl)=\frac{\sqrt{L}}{3{\gamma }_{\_}k}\left[\begin{array}{c}1\\ {\mathrm{e}}^{j\frac{2\pi }{{\lambda }_{\_}c}d\sin ({\varphi }_{\_}^{k}kl)}\\ {\mathrm{e}}^{j2\frac{2\pi }{{\lambda }_{\_}c}d\sin ({\varphi }_{\_}^{k}kl)}\\ {\mathrm{e}}^{-j\frac{2\pi }{{\lambda }_{\_}c}d\sin ({\varphi }_{\_}^{k}kl)}\\ 1\\ {\mathrm{e}}^{j\frac{2\pi }{{\lambda }_{\_}c}d\sin ({\varphi }_{\_}^{k}kl)}\\ {\mathrm{e}}^{-j2\frac{2\pi }{{\lambda }_{\_}c}d\sin ({\varphi }_{\_}^{k}kl)}\\ {\mathrm{e}}^{-j\frac{2\pi }{{\lambda }_{\_}c}d\sin ({\varphi }_{\_}^{k}kl)}\\ 1\end{array}\right].$$

$${\mathbf{a}}_{\_}kl=\frac{\sqrt{L}}{3{\gamma }_{\_}k}[{\mathbf{a}}_{\_}^{T}t({\varphi }_{\_}^{t}kl),{\mathrm{e}}^{-j\frac{2\pi }{{\lambda }_{\_}c}d\sin ({\varphi }_{\_}^{k}kl)},{\mathrm{e}}^{-j2\frac{2\pi }{{\lambda }_{\_}c}d\sin ({\varphi }_{\_}^{k}kl)}{]}^T\in {\mathbb{C}}^{5\times 1},$$

$${\mathbf{t}}_{\_}^{t}kl=\left[\begin{array}{c}[{\mathbf{I}}_{\_}3,0_{\_}3\times 2]\\ [0_{\_}1\times 3,1,0]\\ [{\mathbf{I}}_{\_}2,0_{\_}2\times 3]\\ [0_{\_}1\times 4,1]\\ [0_{\_}1\times 3,1,0]\\ [1,0_{\_}1\times 4]\end{array}\right]{\mathbf{a}}_{\_}kl.$$

$${\mathbf{a}}_{\_}kl(x+N_{\_}t-1)=({\mathbf{a}}_{\_}kl(x){)}^{\ast },2\le x\le N_{\_}t,$$

$${\mathbf{b}}_{\_}kl(y+N_{\_}r-1)=({\mathbf{b}}_{\_}kl(y){)}^{\ast },2\le y\le N_{\_}r.$$

$${\mathbf{t}}_{\_}^{t}kl={\mathbf{\Gamma }}_{\_}u{\mathbf{a}}_{\_}kl,\text{and}{\mathbf{t}}_{\_}^{r}kl={\mathbf{\Gamma }}_{\_}v{\mathbf{b}}_{\_}kl,$$

Full text

I Introduction

Millimeter wave (mmWave) communications are promising for future-generation wireless communications for their advantages such as large bandwidths, narrow beams, and secure transmissions [1, 2]. Large-scale multiple-input multiple-output (MIMO) hybrid structures equipped with only a few RF chains have generated great interests for mmWave systems due to their low complexity and near-optimal performance [3, 4]. To exploit the potential of large-scale MIMO hybrid systems, e.g., to achieve high data transmission rates, precoders and combiners must be carefully designed. They can be designed based on the instantaneous channel matrix [5, 6], which may be estimated by using channel estimation techniques [3, 7]. However, the instantaneous channel can vary fast at mmWave frequencies and the precoder/combiner have to be redesigned once the instantaneous channel changes.

Although the instantaneous mmWave channel can change very fast, the long-term channel statistics, e.g., the angular power spectrum, can be stationary for tens to hundreds of coherence blocks [9]. Recently, the channel covariance information has been utilized to design the analog precoders/combiners [8, 9], which remain fixed when the covariance matrix is unchanged. The effective digital system has a reduced dimensionality, which greatly reduces the cost for acquiring the instantaneous channel state information (CSI) and simplifies the optimization of the digital precoders and combines. To realize the designs in [8, 9], the channel covariance matrix should be firstly estimated. With large antenna arrays, the channel covariance matrix has a large dimensionality, which demands a large number of observations to be used when traditional covariance matrix estimators are adopted. Meanwhile, the hybrid structure only allows a reduced number of observations to be acquired at the receiver, which makes the channel covariance estimation task challenging. In order to address this challenge, [10] proposes several compressive sensing (CS) based channel covariance estimators, which explore the relations between the angle of departure (AoD)/angle of arrival (AoA) and the channel covariance matrix. Their methods need a dictionary for searching the AoD/AoA and the resulting performance improves when the resolution of the dictionary becomes higher. However, high-resolution dictionary yields high computational complexity. Moreover, these CS-based estimators require the number of paths in the channel to be known as a prior. In [11], an analytical expression of the channel covariance matrix is derived and computed through the information obtained from one instantaneous channel realization, which can be estimated from low-dimensional observations. However, the analytical expression is given only for vector channels, which may not be suitable for mmWave communications.

In this paper, we investigate the mmWave channel covariance matrix estimation problem for hybrid mmWave communication systems that are equipped with uniform linear arrays (ULA) or uniform square planar arrays (USPA). The main contributions are as follows:

We show that the mmWave MIMO channel covariance matrix follows a Kronecker product expansion model [12]. Following [12, 13, 14], we show that this model can be used for reducing the effective dimension of the large-dimensional channel covariance matrices in mmWave MIMO systems. We further show that permutation can reduce the rank of the mmWave channel covariance matrix, which admits an expression of the summation of vector outer products. We thus formulate the channel covariance matrix estimation problem as a low-rank matrix sensing problem. 2. 2.

Although the aforementioned low-rank matrix sensing problem has a lower size than the original problem, the complexity can still be high when the numbers of the transmitter/receiver antennas are large. In order to reduce the complexity, we further exploit the structures of the ULA or USPA to reduce the dimensionality of the problem and formulate the problem as a structured low-rank matrix sensing problem. We adapt the recently proposed generalized conditional gradient and alternating minimization (GCG-Alt) algorithm [17], which has a low computational complexity, to find the solution. Numerical results with ULA and USPA suggest that our proposed estimator is effective in estimating the mmWave channel covariance matrix.

The rest of this paper is organized as follows. We introduce the spatial channel model and the hybrid system in Section II. In Section III, we formulate the channel covariance estimation problem as a structured low-rank matrix sensing problem and present the solution. We show the simulation results in Section IV and conclude the paper in Section V.

Notations: Bold uppercase ${\mathbf{A}}$ denotes a matrix and bold lowercase ${\mathbf{a}}$ denotes a column vector. ${\mathbf{A}}^{\ast}$, ${\mathbf{A}}^{T}$, and ${\mathbf{A}}^{H}$ denote the conjugate, transpose, and conjugate transpose of matrix ${\mathbf{A}}$, respectively. ${\mathbf{a}}(i)$ denotes the $i$-th element of vector ${\mathbf{a}}$. $[{\mathbf{A}}]_{\_}{a:b,:}$ denotes the submatrix of ${\mathbf{A}}$ made of its $a$-th to $b$-th rows. $[{\mathbf{A}}]_{\_}{a:b,c:d}$ denotes the submatrix of ${\mathbf{A}}$ defined by its $a$-th to $b$-th rows and $c$-th to $d$-th columns. $\|{\mathbf{A}}\|_{\_}F$ and $\|{\mathbf{A}}\|_{\_}{\ast}$ are the Frobenius norm and the nuclear norm of ${\mathbf{A}}$. For ${\mathbf{A}}\in\mathbb{C}^{M\times N}$, $\text{vec}({\mathbf{A}})\in\mathbb{C}^{MN\times 1}$ is a column vector obtained through the vectorization of ${\mathbf{A}}$ and $\text{vec}^{-1}({\mathbf{A}})\in\mathbb{C}^{M\times N}$ is a matrix obtained by the inverse of vectorization. For matrices ${\mathbf{A}}$ and ${\mathbf{B}}$, ${\mathbf{A}}\otimes{\mathbf{B}}$ denotes the Kronecker product of ${\mathbf{A}}$ and ${\mathbf{B}}$. $\mathcal{CN}(a,b^{2})$ represents complex Gaussian distribution with mean $a$ and variance $b^{2}$. $\mathcal{U}(a,b)$ represents uniform distribution with support $[a,b]$.

II Spatial Channel Model

Consider point-to-point mmWave transmissions, where the transmitter has $N_{\_}t$ antennas and the receiver has $N_{\_}r$ antennas. We assume the following spatial channel [15]:

[TABLE]

where $K$ is the number of clusters, and $L$ is the number of rays within each cluster. As reported in [15], the number of clusters is often small, e.g., $K=1,2$, but the number of rays inside each cluster can be large, e.g., $L=30$. ${\mathbf{a}}_{\_}r(\phi^{r}_{\_}{kl},\theta^{r}_{\_}{kl})$ and ${\mathbf{a}}_{\_}t(\phi^{t}_{\_}{kl},\theta^{t}_{\_}{kl})$ are the array response vectors at the receiver and transmitter, respectively, where $\phi^{r}_{\_}{kl},\theta^{r}_{\_}{kl},\phi^{t}_{\_}{kl}$, and $\theta^{t}_{\_}{kl}$ are the azimuth AoA, elevation AoA, azimuth AoD, and elevation AoD on the $l$-th ray of the $k$-th cluster, respectively. These angles can be characterized by cluster center angles and angular spreads: Each cluster covers a range of angles and the angular spread describes the span of each cluster. The angular spread in the mmWave propagation environment is considered to be small [11]. Measurements of the angular spread taken in the urban area of New York City are presented in [15] in terms of the root-mean-square (rms) of all the measurements. At the carrier frequency $f_{\_}c=28$ GHz, example angular spreads of $15.5^{\circ},6^{\circ},10.2^{\circ}$, and $0^{\circ}$ are reported for $\phi^{r}_{\_}{kl},\theta^{r}_{\_}{kl},\phi^{t}_{\_}{kl}$, and $\theta^{t}_{\_}{kl}$, respectively. The small-scale fading coefficient $g_{\_}{kl}$ is assumed complex Gaussian, i.e., $g_{\_}{kl}\sim\mathcal{CN}(0,\gamma^{2}_{\_}k)$, where $\gamma^{2}_{\_}{k}$ is the fraction power of the $k$-th cluster [15, Eq. (7)].

As discussed in [9], though the small-scale fading gains $\{g_{\_}{kl}\}$ change fast, the AoDs/AoAs and $\gamma^{2}_{\_}k$ may remain stationary over tens to hundreds of coherence blocks. Assume that $\{g_{\_}{kl},\forall k,\forall l\}$ are mutually independent. The channel covariance matrix can then be modeled as

[TABLE]

where

[TABLE]

and

[TABLE]

Note that expression (II) is the same as the channel covariance expression in [9] when $L=1$. In the following, we first present our proposed covariance matrix estimation method for systems equipped with the ULA and then discuss its adaptation to systems that adopt the USPA.

For the ULA, the array responses ${\mathbf{a}}_{\_}t(\phi^{t}_{\_}{kl},\theta^{t}_{\_}{kl})$ and ${\mathbf{a}}_{\_}r(\phi^{r}_{\_}{kl},\theta^{r}_{\_}{kl})$ are independent of the elevation angles. They can thus be abbreviated as ${\mathbf{a}}_{\_}t(\phi^{t}_{\_}{kl})$ and ${\mathbf{a}}_{\_}r(\phi^{r}_{\_}{kl})$. For an $N_{\_}a$-element ULA with distance $d$ between adjacent antennas, the array response is

[TABLE]

where $\lambda_{\_}c$ is the carrier wavelength and $N_{\_}a=N_{\_}t$ or $N_{\_}r$ is the number of antennas at the transmitter or receiver. Accordingly, $\widetilde{{\mathbf{T}}}^{t}_{\_}{{kl}}$ of (3) and $\widetilde{{\mathbf{T}}}^{r}_{\_}{{kl}}$ of (4) become

[TABLE]

and

[TABLE]

respectively, which are Toeplitz-Hermitian. Since the Kronecker product of two Toeplitz-Hermitian matrices is block-Toeplitz-Hermitian [16], the channel covariance matrix ${\mathbf{R}}$ defined in (II) is block-Toeplitz-Hermitian.

We next discuss the hybrid system. We assume phase shifter-based hybrid transceivers [17] shown in Fig. 1, where the antennas and analog phase shifters at the transmitter or receiver are fully connected. Assume that there are $K_{\_}t\ll N_{\_}t$ radio frequency (RF) chains at the transmitter and $K_{\_}r\ll N_{\_}r$ RF chains at the receiver. For single-stream transmissions with one symbol $s$ transmitted, the received signal is written as

[TABLE]

where ${\mathbf{W}}$ and ${\mathbf{f}}$ are the receiving processing matrix and transmitting processing vector, respectively, and ${\mathbf{n}}$ is the noise vector. Up to $K_{\_}r$ digital symbols can be observed at the receiver after each transmission. In hybrid transceivers, we have ${\mathbf{W}}={\mathbf{W}}_{\_}{\rm{RF}}{\mathbf{W}}_{\_}{\rm{BB}}$ and ${\mathbf{f}}={\mathbf{F}}_{\_}{\rm{RF}}{\mathbf{f}}_{\_}{\rm{BB}}$, where ${\mathbf{W}}_{\_}{\rm{RF}}$ and ${\mathbf{F}}_{\_}{\rm{RF}}$ are the analog combiner and precoder, respectively, and ${\mathbf{W}}_{\_}{\rm{BB}}$ and ${\mathbf{f}}_{\_}{\rm{BB}}$ are the digital combiner and precoder, respectively. In addition, due to the constraints of the phase shifters in the RF combiner and precoder, the entries in ${\mathbf{W}}_{\_}{\rm{RF}}$ and ${\mathbf{F}}_{\_}{\rm{RF}}$ have constant modulus. Note that using single-stream transmissions during the channel training avoids the interferences caused by transmitting multiple symbols simultaneously and this has been widely considered [3, 4, 7].

When $N_{\_}t$ and $N_{\_}r$ are large, the dimension of the channel covariance matrix ${\mathbf{R}}$ is large. In this case, estimating ${\mathbf{R}}$ can be difficult when only a small number of observations available, which is typical in the hybrid system. From (II), ${\mathbf{R}}$ follows the Kronecker product expansion model [12]. In the following, we explore this property and the block-Toeplitz-Hermitian structure of ${\mathbf{R}}$ to reduce the dimensionality of the problem of estimating ${\mathbf{R}}$, and formulate the channel covariance matrix estimation problem as a structured low-rank matrix sensing problem.

III Structured Low-Rank Covariance Matrix Sensing

III-A Rank Reduction by Permutation

Define

[TABLE]

and

[TABLE]

respectively, where $1\leq l\leq L$ and $1\leq k\leq K$. Then ${\mathbf{R}}$ of (II) can be written compactly as

[TABLE]

where the summation involves $KL$ terms. Note that ${\mathbf{T}}^{t}_{\_}{kl}$ and ${\mathbf{T}}^{r}_{\_}{kl}$ are Toeplitz-Hermitian. Denote the following $N_{\_}r\times N_{\_}r$ submatrix of ${\mathbf{R}}$ as

[TABLE]

where $1\leq m\leq N_{\_}t$ and $1\leq n\leq N_{\_}t$. Define a permutation operator $\mathcal{P}(\cdot)$ that permutes the $N_{\_}tN_{\_}r\times N_{\_}tN_{\_}r$ matrix ${\mathbf{R}}$ into a $N_{\_}t^{2}\times N_{\_}r^{2}$ matrix

[TABLE]

by stacking each submatrix ${\mathbf{R}}_{\_}{mn}$ into a row vector as

[TABLE]

We write

[TABLE]

and

[TABLE]

Then based on the Kronecker product expansion property [16], [18], ${\mathbf{R}}_{\_}p$ can be written as a sum of vector outer products

[TABLE]

Note that if we have ${\mathbf{R}}_{\_}p$, we can obtain ${\mathbf{R}}$ as $\mathcal{P}^{-1}({\mathbf{R}}_{\_}p)$.

From (12), the column space of ${\mathbf{R}}_{\_}p$ is spanned by $\{{\mathbf{t}}^{t}_{\_}{kl}\}$ and the row space of ${\mathbf{R}}_{\_}{p}$ is spanned by $\{{\mathbf{t}}^{r}_{\_}{kl}\}$. Recall that ${\mathbf{t}}^{t}_{\_}{kl}=\text{vec}({\mathbf{T}}^{t}_{\_}{kl})$ and by using the relation between ${\mathbf{T}}^{t}_{\_}{kl}$ and the transmitter array response ${\mathbf{a}}_{\_}t(\phi^{t}_{\_}{kl})$ shown in (5) and (8), ${\mathbf{t}}^{t}_{\_}{kl}$ can be written as

[TABLE]

where $\phi^{t}_{\_}{kl}$ is the azimuth AoD. We can see that ${\mathbf{t}}^{t}_{\_}{kl}$ consists of the array response vector ${\mathbf{a}}_{\_}{t}(\phi^{t}_{\_}{kl})$ and the column space of ${\mathbf{R}}_{\_}p$ is determined by the set

[TABLE]

As introduced earlier, small angular spreads are observed in the mmWave propagation environment, which indicates that the AoDs inside a cluster are closely spaced and their corresponding array response vectors are highly correlated. Therefore, for the $k$-th cluster, though the number of rays $L$ inside can be large, the space spanned by $\{{\mathbf{a}}^{\ast}_{\_}t(\phi^{t}_{\_}{kl})\otimes{\mathbf{a}}_{\_}t(\phi^{t}_{\_}{kl}),1\leq l\leq L\}$ may be well approximated by a low-rank space. In addition, since the number of clusters $K$ is generally small (e.g., $K=1$ or $2$), both $\bf{\mathcal{C}}_{\_}t$ and ${\mathbf{R}}$ can be low-rank. This is similar to the low-rankness of the mmWave channel ${\mathbf{H}}$, which has been validated by the experimental and simulation results in [15].

The low-rank property of ${\mathbf{R}}_{\_}p$ can be shown numerically. Denote by $r_{\_}{\rm{ch}}$ the rank of ${\mathbf{R}}_{\_}p$ or ${\mathbf{R}}$, and let $\sigma_{\_}1>\sigma_{\_}2>\ldots>\sigma_{\_}{r_{\_}{\rm{ch}}}$ be the singular values of ${\mathbf{R}}_{\_}p$ or ${\mathbf{R}}$. We may use

[TABLE]

to measure the energy captured by a rank-$r_{\_}{\rm{sub}}$ approximation of ${\mathbf{R}}_{\_}p$ or ${\mathbf{R}}$, where $r_{\_}{\rm{sub}}$ is the rank of the subspace of ${\mathbf{R}}_{\_}p$ or ${\mathbf{R}}$.

Fig. 2 shows an example of a ULA system with $K=\{1,2,3,4\},L=30$, $N_{\_}t=64$, and $N_{\_}r=16$. The covariance matrix ${\mathbf{R}}$ and its permuted version ${\mathbf{R}}_{\_}p$ have sizes of $1024\times 1024$ and $4096\times 256$, respectively, which shows that ${\mathbf{R}}_{\_}p$ is a taller matrix. The horizontal AoDs

[TABLE]

where the center angles $\phi^{t}_{\_}k$ are distributed uniformly in $[0,2\pi]$ and separated by at least one angular spread $\upsilon^{t}_{\_}h=10.2^{\circ}$. Similarly, the horizontal AoAs

[TABLE]

where $\upsilon^{r}_{\_}h=15.5^{\circ}$. The cluster powers are generated following [15, Tab. I]. It can be seen from Fig. 2 that for capturing a majority of the total energy, e.g., with $p_{\_}e=0.95,0.99$, the required $r_{\_}{\rm{sub}}$ for ${\mathbf{R}}_{\_}p$ is generally much smaller than $\text{min}(N^{2}_{\_}t,N^{2}_{\_}r,KL)$ and is also much smaller than that for ${\mathbf{R}}$. In the following, we use $r_{\_}{p}$ as the $r_{\_}{\rm{sub}}$ of ${\mathbf{R}}_{\_}p$ and $r_{\_}{\rm{R}}$ as the $r_{\_}{\rm{sub}}$ of ${\mathbf{R}}$ for a certain $p_{\_}e$. As such, ${\mathbf{R}}_{\_}p$ may be well approximated as a rank-$r_{\_}p$ matrix. One may use low-rank matrix recovery methods, e.g., matrix completion methods, to estimate the best rank-$r_{\_}p$ approximation of ${\mathbf{R}}_{\_}p$ from a small amount of observations. However, when $N_{\_}t$ and $N_{\_}r$ are large, which is the case in mmWave communications, the number of parameters required by the rank-$r_{\_}p$ approximation of ${\mathbf{R}}_{\_}p$, i.e., $(N^{2}_{\_}t+N_{\_}t^{2})\times r_{\_}p$, is still large. Therefore, estimating the subspaces of ${\mathbf{R}}_{\_}{p}$ can be computationally expensive.

III-B Dimension Reduction By Exploiting the Toeplitz-Hermitian Structure

Recall that ${\mathbf{R}}$ is block-Toeplitz-Hermitian and ${\mathbf{R}}_{\_}p=\mathcal{P}({\mathbf{R}})$ is a permutation of ${\mathbf{R}}$. From (12) and (III-A), we can see that ${\mathbf{R}}_{\_}p$ is also specially structured: ${\mathbf{R}}_{\_}p$ is the summation of the outer products of ${\mathbf{t}}^{t}_{\_}{kl}$ and ${\mathbf{t}}^{r}_{\_}{kl}$, where ${\mathbf{t}}^{t}_{\_}{kl}$ and ${\mathbf{t}}^{r}_{\_}{kl}$ are the vectorizations of Toeplitz-Hermitian matrices ${\mathbf{T}}^{t}_{\_}{kl}$ and ${\mathbf{T}}^{r}_{\_}{kl}$, respectively. Since the Toeplitz-Hermitian matrix ${\mathbf{T}}^{t}_{\_}{kl}\in\mathbb{C}^{N^{2}_{\_}t\times N^{2}_{\_}t}$ is determined by its first column and first row (its first row is the conjugate transpose of its first column), we can represent ${\mathbf{t}}^{t}_{\_}{kl}$ in terms of the entries in the first column and first row of ${\mathbf{T}}^{t}_{\_}{kl}$. We can represent ${\mathbf{t}}^{r}_{\_}{kl}$ in the same way. Therefore, the total numbers of unknowns in ${\mathbf{t}}^{t}_{\_}{kl}$ and ${\mathbf{t}}^{r}_{\_}{kl}$ are $2N_{\_}t-1$ and $2N_{\_}r-1$, respectively. Then we can reduce the problem size of $(N^{2}_{\_}t+N_{\_}r^{2})\times r_{\_}p$ to $2(N_{\_}t+N_{\_}r-1)\times r_{\_}p$. In the following, we show how the problem size can be reduced.

First, let us use an example with $N_{\_}t=3$ to illustrate the structure of ${\mathbf{t}}^{t}_{\_}{{kl}}$. The array response

[TABLE]

Then according to (III-A), we have

[TABLE]

We can see that all the $9$ elements in ${\mathbf{t}}^{t}_{\_}{kl}$ can be represented by the elements in ${\mathbf{a}}_{\_}t(\phi^{t}_{\_}{kl})$ and ${\mathbf{a}}^{\ast}_{\_}t(\phi^{t}_{\_}{kl})$. Now construct a vector

[TABLE]

then ${\mathbf{t}}^{t}_{\_}{kl}$ can be rewritten as

[TABLE]

In fact, ${\mathbf{a}}_{\_}{kl}(4)=({\mathbf{a}}_{\_}{kl}(2))^{\ast}$ and ${\mathbf{a}}_{\_}{kl}(5)=({\mathbf{a}}_{\_}{kl}(3))^{\ast}$. Therefore, ${\mathbf{t}}^{t}_{\_}{kl}$ can be expressed as a product of a weight matrix and a vector ${\mathbf{a}}_{\_}{kl}$. Furthermore, the weight matrix depends only on the structure of the antenna array and is independent of the path angles.

Similarly, for the general cases, we can express ${\mathbf{t}}^{t}_{\_}{kl}$ with a weight matrix ${\mathbf{\Gamma}}_{\_}u\in\mathbb{C}^{N^{2}_{\_}t\times(2N_{\_}t-1)}$ and a vector ${\mathbf{a}}_{\_}{kl}\in\mathbb{C}^{(2N_{\_}t-1)\times 1}$, and express ${\mathbf{t}}^{r}_{\_}{kl}$ with a weight matrix ${\mathbf{\Gamma}}_{\_}v\in\mathbb{C}^{N^{2}_{\_}r\times(2N_{\_}r-1)}$ and a vector ${\mathbf{b}}_{\_}{kl}\in\mathbb{C}^{(2N_{\_}r-1)\times 1}$. We require

[TABLE]

and

[TABLE]

We then have

[TABLE]

where ${\mathbf{\Gamma}}_{\_}u=[{\mathbf{\Gamma}}_{\_}{u1},{\mathbf{\Gamma}}_{\_}{u2}]$ with ${\mathbf{\Gamma}}_{\_}{u1}\in\mathbb{C}^{N^{2}_{\_}{t}\times N_{\_}{t}}$, ${\mathbf{\Gamma}}_{\_}{u2}\in\mathbb{C}^{N^{2}_{\_}{t}\times(N_{\_}{t}-1)}$, and

[TABLE]

with ${\mathbf{e}}_{\_}i\in\mathbb{C}^{(N_{\_}{t}-1)\times 1}$ being a vector whose $i$-th entry is $1$ and other entries are zero. ${\mathbf{\Gamma}}_{\_}v$ is constructed similarly as ${\mathbf{\Gamma}}_{\_}u$, and ${\mathbf{\Gamma}}_{\_}u$ and ${\mathbf{\Gamma}}_{\_}v$ are both full-rank. This is because ${\mathbf{\Gamma}}_{\_}u$ and ${\mathbf{\Gamma}}_{\_}v$ consist of $1$’s and [math]’s, and there is only one $1$ in each row of ${\mathbf{\Gamma}}_{\_}u$ and ${\mathbf{\Gamma}}_{\_}v$. Therefore, (12) can be rewritten as

[TABLE]

where ${\mathbf{C}}=\sum_{\_}{k=1}^{K}\sum_{\_}{l=1}^{L}{\mathbf{a}}_{\_}{kl}{\mathbf{b}}^{T}_{\_}{kl}$.

As shown above, ${\mathbf{R}}_{\_}p$ is approximately low-rank. Since the fixed weight matrices ${\mathbf{\Gamma}}_{\_}u$ and ${\mathbf{\Gamma}}_{\_}v$ are full-rank, ${\mathbf{C}}$ is approximately low-rank. Hence estimating a low-rank approximation of ${\mathbf{R}}_{\_}{p}$ is equivalent to estimating a low-rank approximation of ${\mathbf{C}}$. Note that ${\mathbf{C}}\in\mathbb{C}^{(2N_{\_}t-1)\times(2N_{\_}r-1)}$ is much smaller than ${\mathbf{R}}_{\_}p\in\mathbb{C}^{N_{\_}t^{2}\times N_{\_}r^{2}}$ and this can greatly reduce the complexity of the problem.

III-C Training

We assume that the channel matrix ${\mathbf{H}}$ remains static during a snapshot and suppose we have $T$ snapshots. For different snapshots, we assume that the AoAs/AoDs and the fraction power $\gamma_{\_}k^{2}$ remain unchanged, but the small-scale fading gain $g_{\_}{kl}\sim\mathcal{CN}(0,\gamma^{2}_{\_}k)$ can change [9]. Suppose the transmitter sends out $S$ training beams during each snapshot. For the $s$-th training beam of the $t$-th snapshot, we employ the transmitting vector ${\mathbf{f}}_{\_}{{t,s}}\in\mathbb{C}^{N_{\_}{t}}$ and the receiving matrix ${\mathbf{W}}_{\_}{{t,s}}\in\mathbb{C}^{N_{\_}{r}\times K_{\_}{r}}$. Therefore, in each snapshot, after the transmitter sends out $S$ training beams, the receiver receives $SK_{\_}r$ symbols and the sampling ratio is $SK_{\_}r/N_{\_}rN_{\_}t$. We design ${\mathbf{f}}_{\_}{t,s}$ and ${\mathbf{W}}_{\_}{t,s}$ and their corresponding ${\mathbf{F}}_{\_}{\rm{RF}},{\mathbf{f}}_{\_}{\rm{BB}},{\mathbf{W}}_{\_}{\rm{RF}}$, and ${\mathbf{W}}_{\_}{\rm{BB}}$ realizations for the hybrid structure according to the training scheme in [17, Section III.D]. For the $s$-th training beam of the $t$-th snapshot, the received signal is

[TABLE]

where ${\mathbf{n}}_{\_}{t,s}$ is the noise vector and ${{\mathbf{H}}_{\_}t}$ is the channel matrix at snapshot $t$. Without loss of generality, assume identical training symbols $s=\sqrt{P}$. By setting $\|{\mathbf{f}}_{\_}{t,s}\|^{2}_{\_}F=1$, the total transmitting power is $\|{\mathbf{f}}_{\_}{t,s}s\|_{\_}F^{2}=P$ and the pilot-to-noise ratio $\rm{(PNR)}$ is defined as

[TABLE]

where the noise is assumed to be an additive white Gaussian noise (AWGN) with variance $\sigma^{2}$. In the $t$-th snapshot and after the transmitter sends out all the $S$ training beams, stack the received signals as

[TABLE]

where

[TABLE]

Suppose the trainings are the same for different snapshots, i.e., ${\mathbf{f}}_{\_}{1,s}={\mathbf{f}}_{\_}{2,s}=\ldots={\mathbf{f}}_{\_}{T,s}={\mathbf{f}}_{\_}{s}$, ${\mathbf{W}}_{\_}{1,s}={\mathbf{W}}_{\_}{2,s}=\ldots={\mathbf{W}}_{\_}{T,s}={\mathbf{W}}_{\_}{s}$. We then have

[TABLE]

and

[TABLE]

where ${\mathbf{\Sigma}}$ and ${\mathbf{\Sigma}}_{\_}{n}$ represent the covariance matrices of the received signal ${\mathbf{y}}_{\_}t$ and the noise ${\mathbf{n}}_{\_}t$, respectively. After $T$ snapshots, we can compute the dimension-reduced sample covariance matrix (SCM) of ${\mathbf{y}}_{\_}t$ as

[TABLE]

We permute ${\mathbf{S}}$ into ${\mathbf{S}}_{\_}{p}\in\mathbb{C}^{S^{2}\times K^{2}_{\_}{r}}$ in a similar procedure as ${\mathbf{R}}$ is permuted into ${\mathbf{R}}_{\_}p$.

III-D Low-Rank Matrix Sensing Problem

We can now formulate the channel covariance estimation problem as a low-rank matrix sensing problem [20]:

[TABLE]

where $\widehat{{\mathbf{R}}}_{\_}{p}$ is the estimate of ${\mathbf{R}}_{\_}{p}$, $\mathcal{A}:\mathbb{C}^{N^{2}_{\_}t\times N^{2}_{\_}r}\rightarrow\mathbb{C}^{S^{2}K^{2}_{\_}r\times 1}$ is an appropriate linear map, and $\zeta^{2}$ is a constant to account for the fitting error. Replacing $\widehat{{\mathbf{R}}}_{\_}{p}$ with (III-B), we can reformulate (26) as

[TABLE]

where $\widehat{{\mathbf{C}}}$ is the estimate of ${\mathbf{C}}$.

In general, problem (27) is a nonconvex optimization problem and difficult to solve. In this paper, we solve the relaxed version of problem (27) [19]:

[TABLE]

where

[TABLE]

and $\mu>0$ is a regularization coefficient. After some manipulations, we have $\mathcal{A}({\mathbf{\Gamma}}_{\_}u\widehat{{\mathbf{C}}}{\mathbf{\Gamma}}^{T}_{\_}v)={\mathbf{Q}}\text{vec}(\widehat{{\mathbf{C}}})$, where

[TABLE]

and ${\mathbf{Q}}\in\mathbb{C}^{S^{2}K^{2}_{\_}r\times(2N_{\_}t-1)(2N_{\_}r-1)}$. The direct evaluation of $\|\widehat{{\mathbf{C}}}\|_{\_}{\ast}$, which is the nuclear norm (i.e., the summation of the singular values) of $\widehat{{\mathbf{C}}}$, is computationally expensive. Following [17], $\|\widehat{{\mathbf{C}}}\|_{\_}{\ast}$ can be written as

[TABLE]

Therefore, finding a $\widehat{{\mathbf{C}}}$ to minimize the objective function in (28) becomes finding a pair of $({\mathbf{U}},{\mathbf{V}})$ to minimize

[TABLE]

A similar low-rank recovery problem is recently studied in [17] for instantaneous mmWave channel estimation, where a training scheme is designed such that the channel can be estimated by solving a matrix completion (MC) problem. A generalized conditional gradient and alternating minimization (GCG-Alt) algorithm is developed, which is shown to be able to provide accurate low-rank solutions at a low complexity. In this work, we adapt the GCG-Alt algorithm to solve (III-D) for our covariance matrix estimation problem. In the following, we discuss the key steps of the GCG-Alt algorithm for solving (III-D) and refer the reader to [17] for more detailed treatments.

The GCG-Alt algorithm consists of a relaxed GCG algorithm and an AltMin algorithm. Let $\widehat{{\mathbf{C}}}_{\_}{k-1}$ be the solution to ${\mathbf{C}}$ at the $(k-1)$-th GCG iteration. The relaxed GCG algorithm first produces an output

[TABLE]

where ${\mathbf{Z}}_{\_}k$ is the outer product of the top singular vector pair of $\text{vec}^{-1}(-\nabla f(\widehat{{\mathbf{C}}}_{\_}k))$. The calculations of $\text{vec}^{-1}(-\nabla f(\widehat{{\mathbf{C}}}_{\_}k))$ and the parameter $\theta_{\_}k$ here are different from those in [17]. For problem (III-D), we calculate $\text{vec}^{-1}(-\nabla f(\widehat{{\mathbf{C}}}_{\_}k))$ as

[TABLE]

and the parameter $\theta_{\_}k$ as

[TABLE]

where ${\mathbf{q}}_{\_}{z_{\_}k}={\mathbf{Q}}\text{vec}({\mathbf{Z}}_{\_}k)$ and $\mathcal{R}(\cdot)$ denotes the real part of a number. Since $\widehat{{\mathbf{C}}}_{\_}k={\mathbf{U}}_{\_}k{\mathbf{V}}^{T}_{\_}k$, updating $\widehat{{\mathbf{C}}}_{\_}k$ is equivalent to updating ${\mathbf{U}}_{\_}k=[\sqrt{1-\eta_{\_}k}{\mathbf{U}}_{\_}{k-1},\sqrt{\theta_{\_}k}{\mathbf{u}}_{\_}k]$ and ${\mathbf{V}}_{\_}k=[\sqrt{1-\eta_{\_}k}{\mathbf{V}}_{\_}{k-1},\sqrt{\theta_{\_}k}{\mathbf{v}}_{\_}k]$. Then the obtained ${\mathbf{U}}_{\_}k$ and ${\mathbf{V}}_{\_}k$ are used as the initial input of the AltMin algorithm, i.e., ${\mathbf{U}}^{0}_{\_}k\leftarrow{\mathbf{U}}_{\_}k,{\mathbf{V}}^{0}_{\_}k\leftarrow{\mathbf{V}}_{\_}k$. After $I_{\_}{a}$ iterations of the AltMin algorithm, update ${\mathbf{U}}_{\_}k={\mathbf{U}}^{I_{\_}{a}}_{\_}k$ and ${\mathbf{V}}_{\_}k={\mathbf{V}}^{I_{\_}{a}}_{\_}k$. For completeness, we summarize the GCG-Alt algorithm in Algorithm I. After obtaining $\widehat{{\mathbf{C}}}$, we have $\widehat{{\mathbf{R}}}^{r}_{\_}{p}={\mathbf{\Gamma}}_{\_}u\widehat{{\mathbf{C}}}{\mathbf{\Gamma}}^{T}_{\_}u$ and $\widehat{{\mathbf{R}}}=\mathcal{P}^{-1}(\widehat{{\mathbf{R}}}^{r}_{\_}{p})$.

III-E Computational Complexity

Define a flop as an operation of real-valued numbers. Let $M=SK_{\_}r$ be the number of received symbols during each snapshot. Following the computational complexity analysis in [17], the computational complexity of the GCG-Alt estimator is about $8r_{\_}{\rm{est}}(I_{\_}ar_{\_}{\rm{est}}+I_{\_}a+1)(2N_{\_}t-1)^{2}(2N_{\_}r-1)^{2}+8/3I_{\_}ar_{\_}{\rm{est}}(r_{\_}{\rm{est}}+1)(2r_{\_}{\rm{est}}+1)(2N_{\_}t-1)(2N_{\_}r-1)(N_{\_}r+N_{\_}t-1)+I_{\_}ar^{2}_{\_}{\rm{est}}(r_{\_}{\rm{est}}+1)^{2}((2N_{\_}r-1)^{3}+(2N_{\_}t-1)^{3})+16r_{\_}{\rm{est}}(2N_{\_}r-1)(2N_{\_}t-1)M^{2}$, where $I_{\_}a$ is the number of iterations of the AltMin algorithm and $r_{\_}{\rm{est}}$ is the estimated rank of $\widehat{{\mathbf{C}}}$ by the GCG-Alt estimator. Later in Section IV, we show the computational complexity of the GCG-Alt estimator with specific examples.

III-F Extension to the USPA System

We now follow the same process introduced in Section III. A-D to estimate the channel covariance matrix for USPA systems. To account for the different array structure of the USPA, the weight matrices of (17) are redesigned. For a $\sqrt{N}_{\_}a\times\sqrt{N}_{\_}a$ USPA placed on the $yz$ plane with distance $d$ between adjacent antennas, the array response is

[TABLE]

where

[TABLE]

[TABLE]

is the array response along the $y$ axis and

[TABLE]

is the array response along the $z$ axis. We design the weight matrices by examining the structure of $\widetilde{{\mathbf{T}}}^{t}_{\_}{kl}$ defined in (3) which is written as

[TABLE]

where ${\mathbf{a}}_{\_}{t_{\_}{y}}(\phi^{t}_{\_}{kl},\theta^{t}_{\_}{kl})$ and ${\mathbf{a}}_{\_}{t_{\_}{z}}(\theta^{t}_{\_}{kl})$ are the transmitter array response vectors along the $y$ axis and $z$ axis, respectively. Note that $\widetilde{{\mathbf{T}}}^{t}_{\_}{{kl}}$ is block-Toeplitz-Hermitian. Let

[TABLE]

and

[TABLE]

we can verify that

[TABLE]

and $\widetilde{{\mathbf{T}}}^{y}_{\_}{{kl}}$ and $\widetilde{{\mathbf{T}}}^{z}_{\_}{{kl}}$ are Toeplitz-Hermitian. Then for the USPA, ${\mathbf{T}}^{t}_{\_}{kl}$ of (8) can be written as

[TABLE]

In Section III. B, we have expressed ${\mathbf{t}}^{t}_{\_}{kl}=\text{vec}^{-1}({\mathbf{T}}^{t}_{\_}{kl})$, where ${\mathbf{T}}^{t}_{\_}{kl}$ of (8) is Toeplitz-Hermitian matrix, in terms of a weight matrix and a vector. We have similar expressions for the vectorizations of the Toeplitz-Hermitian matrices $\widetilde{{\mathbf{T}}}^{y}_{\_}{{kl}}$ and $\widetilde{{\mathbf{T}}}^{z}_{\_}{{kl}}$. Let

[TABLE]

where ${\mathbf{\Gamma}}_{\_}y\in\mathbb{C}^{N_{\_}t\times(2\sqrt{N_{\_}t}-1)}$ is the weight matrix and ${\mathbf{a}}^{y}_{\_}{kl}\in\mathbb{C}^{(2\sqrt{N_{\_}t}-1)\times 1}$, and

[TABLE]

where ${\mathbf{\Gamma}}_{\_}{z}\in\mathbb{C}^{N_{\_}t\times(2\sqrt{N_{\_}t}-1)}$ is the weight matrix and ${\mathbf{a}}^{z}_{\_}{kl}\in\mathbb{C}^{(2\sqrt{N_{\_}t}-1)\times 1}$. Let

[TABLE]

and

[TABLE]

By exploring the matrix vectorization process, we have

[TABLE]

where

[TABLE]

is the weight matrix and

[TABLE]

is a vector. Then for the USPA system, ${\mathbf{\Gamma}}_{\_}u$ of (17) becomes (41) and ${\mathbf{\Gamma}}_{\_}v$ of (17) is constructed similarly as (41); the sizes of vectors ${\mathbf{a}}_{\_}{kl}$ and ${\mathbf{b}}_{\_}{kl}$ of (17) have changed: ${\mathbf{a}}_{\_}{kl}\in\mathbb{C}^{(2\sqrt{N_{\_}t}-1)^{2}\times 1}$ and ${\mathbf{b}}_{\_}{kl}\in\mathbb{C}^{(2\sqrt{N_{\_}r}-1)^{2}\times 1}$, and consequently, the size for matrix ${\mathbf{C}}$ of (III-B) has changed: ${\mathbf{C}}\in\mathbb{C}^{(2\sqrt{N_{\_}t}-1)^{2}\times(2\sqrt{N_{\_}r}-1)^{2}}$. After obtaining the weight matrices, we can follow the processes in Section III. C-D to estimate ${\mathbf{C}}$ and then have the channel covariance matrix estimated as $\widehat{{\mathbf{R}}}=\mathcal{P}^{-1}({\mathbf{\Gamma}}_{\_}u\widehat{{\mathbf{C}}}{\mathbf{\Gamma}}^{T}_{\_}v)$.

IV Simulations

We now evaluate the performance of our proposed design for fully connected hybrid transceivers with the ULA and USPA.

IV-A The ULA system

We assume a carrier frequency of $f_{\_}c=28$ GHz. For the ULA system, $N_{\_}t=64,N_{\_}r=16,K_{\_}t=16$, and $K_{\_}r=4$. The number of clusters $K=\{1,2\}$ and there are $L=30$ rays in each cluster. The horizontal AoDs and AoAs are generated as (15) with $\upsilon^{t}_{\_}h=10.2^{\circ}$ and as (16) with $\upsilon^{r}_{\_}h=15.5^{\circ}$, respectively. The cluster powers are generated following [15, Tab. I]. We compare the GCG-Alt estimator with the DCOMP estimator in [10], which has varying receiving processing matrix ${\mathbf{W}}_{\_}{t,s}$ and transmitting processing vector ${\mathbf{f}}_{\_}{t,s}$ during training and has the best performance among other estimators in [10]. The DCOMP estimator needs a dictionary matrix with $G_{\_}t$ grid points that is associated with AoD and a dictionary matrix with $G_{\_}r$ grid points that is associated with AoA. Let $L_{\_}p$ be the number of paths in the channel, the DCOMP estimator assumes that $L_{\_}p$ is known. For the DCOMP estimator, we set $G_{\_}t=2N_{\_}t=128,G_{\_}r=2N_{\_}r=32$, and $L_{\_}p=r_{\_}{\rm{R}}$. Based on Fig. 2, for $p_{\_}e=0.99$, $r_{\_}{\rm{R}}=18$ and $24$ for $K=1$ and $2$, respectively. For the GCG-Alt estimator, we set $\mu=\sigma^{2}$, $\epsilon=0.003$, and $\epsilon_{\_}a=0.1$. The performance metric $\eta$ [10]

[TABLE]

is used to measure how close the subspace of $\widehat{{\mathbf{R}}}$ is to the subspace of ${\mathbf{R}}$, where $\widehat{{\mathbf{M}}}\in\mathbb{C}^{N_{\_}tN_{\_}r\times r_{\_}{\rm{R}}}$ and ${\mathbf{M}}\in\mathbb{C}^{N_{\_}tN_{\_}r\times r_{\_}{\rm{R}}}$ are the singular vector matrices of $\widehat{{\mathbf{R}}}$ and ${\mathbf{R}}$, respectively. We also use the average of the normalized mean square error

[TABLE]

to measure their performance.

We set ${\rm{PNR}}=10$ dB and the number of training beams $S=32$, and compare the GCG-Alt estimator with the DCOMP estimator under different $T$. With $S=32$ per snapshot, the sampling ratio at each snapshot is $SK_{\_}r/N_{\_}rN_{\_}t=12.5\%$. The comparison result shown in Fig. 3 suggests that when the sampling ratio per snapshot is $12.5\%$, our proposed estimator requires fewer snapshots to obtain a $\widehat{{\mathbf{R}}}$ whose subspace is close to that of ${\mathbf{R}}$, as compared to the DCOMP estimator. The $\rm{NMSE}$ result shown in Fig. 4 suggests that our proposed GCG-Alt estimator can obtain a more accurate covariance matrix estimate.

We also compare the computational complexity of the GCG-Alt estimator and the DCOMP estimator. The computational complexity of the DCOMP estimator is about $8TL_{\_}pG_{\_}tG_{\_}r(M^{2}+M)$ flops, where $M=SK_{\_}r$. For the GCG-Alt estimator, based on our observations, the number of iterations of the AltMin algorithm $I_{\_}a\leq 2$, the estimated rank $r_{\_}{\rm{est}}\approx 4$ when $K=1$ and $r_{\_}{\rm{est}}\approx 5$ when $K=2$. Fig. 5 shows the comparison results with different $T$. We can see that the computational complexity of the GCG-Alt estimator is lower than the DCOMP estimator. Also, the computational complexity of the GCG-Alt estimator does not increase as $T$ increases. This is because we use ${\mathbf{S}}_{\_}p\in\mathbb{C}^{S^{2}\times K^{2}_{\_}r}$, which is the permutation of the SCM of ${\mathbf{y}}_{\_}t$ shown in (25), and its size is irrelevant to $T$.

Then we set the number of snapshots $T=40$, and compare the GCG-Alt estimator with the DCOMP estimator under different $S$. The result shown in Fig. 6 suggests that when $T=40$, the GCG-Alt estimator can obtain a more accurate subspace estimation than the DCOMP estimator when the number of training beams $S\geq 24$ per snapshot. Note that $S=24$ corresponds to a sampling ratio of $9.375\%$ per snapshot.

The GCG-Alt estimator explores both the Kronecker structure and the block-Toeplitz-Hermitian structure of ${\mathbf{R}}$ while the DCOMP estimator only considers the Hermitian structure of ${\mathbf{R}}$, so the GCG-Alt estimator can reach an accurate subspace estimation of ${\mathbf{R}}$ with fewer snapshots. We use the same training for different snapshots while the DCOMP estimator uses different trainings per snapshot (i.e., varying ${\mathbf{W}}_{\_}{t,s}$ and ${\mathbf{f}}_{\_}{t,s}$ ). When $S$ is small (e.g., $S\leq 16$), the DCOMP estimator outperforms the GCG-Alt estimator. However, the GCG-Alt estimator performs better when $S$ becomes larger (e.g., $S\geq 24$). Note that for the DCOMP estimator, estimating more paths (i.e., $L_{\_}p$ is large) yields better performance, but its computational complexity also increases.

IV-B The USPA system

We next consider the system with the USPA at the transmitter and receiver. The parameters $f_{\_}c,K,L,\phi^{t}_{\_}{kl},\phi^{r}_{\_}{kl}$ are assumed the same as in the ULA system. The transmitter has an $8\times 8$ USPA (i.e., $N_{\_}t=64$) and $K_{\_}t=16$ RF chains, and the receiver has a $4\times 4$ USPA (i.e., $N_{\_}r=16$) and $K_{\_}r=4$ RF chains. We assume the elevation AoD angular spread $\upsilon^{t}_{\_}v=0^{\circ}$ and the elevation AoA angular spread $\upsilon^{r}_{\_}v=6^{\circ}$ based on the measurement results in [15]. The elevation AoDs and AoAs are distributed as

[TABLE]

[TABLE]

with the elevation center angles $\theta^{t}_{\_}k$ and $\theta^{r}_{\_}k$ being generated in the same manner as the azimuth center angles in the ULA system. For the DCOMP estimator, we set $G_{\_}t=2\sqrt{N_{\_}t}\times 2\sqrt{N_{\_}t}=256$ and $G_{\_}r=2\sqrt{N_{\_}r}\times 2\sqrt{N_{\_}r}=64$. The parameters $L_{\_}p,\mu,\epsilon,\epsilon_{\_}a$ for the GCG-Alt estimator and DCOMP estimator are the same as in the ULA system.

We set ${\rm{PNR}}=10$ dB. The performance comparison with $S=32$ under different $T$ is shown in Fig. 7 and the performance comparison with $T=40$ under different $S$ is shown in Fig. 8. We can see that both of the GCG-Alt estimator and the DCOMP estimator achieve higher $\eta$ for the USPA system. One reason for this is that the USPA system has lower resolution than the ULA system in the azimuth direction even though they have the same number of transmitter and receiver antennas. For the USPA system, the azimuth AoD is resolved by an $\sqrt{N_{\_}t}=8$-element antenna array and the azimuth AoA is resolved by a $\sqrt{N_{\_}r}=4$-element antenna array; while for the ULA system, the azimuth AoD is resolved by a $N_{\_}t=64$-element antenna array and the azimuth AoA is resolved by a $N_{\_}r=16$-element antenna array. Therefore, for the same angular spread, the USPA system resolves fewer paths than the ULA system, which results in a lower rank.

We also show the effects of angular spreads on the performance of the estimators. We set $\upsilon^{t}_{\_}v=0^{\circ},\upsilon^{r}_{\_}v=6^{\circ},K=1,{\rm{PNR}}=10$ dB, $S=16$, and $T=16$. The estimators’ performance under different angular spreads for the azimuth AoD/AoA (i.e., different $\upsilon^{t}_{\_}h$ and $\upsilon^{r}_{\_}h$) shown in Fig. 9 suggests that the estimators achieve lower $\eta$ when $\upsilon^{t}_{\_}h$ and $\upsilon^{r}_{\_}h$ are larger.

V Conclusions

We have formulated the channel covariance estimation problem for hybrid mmWave systems as a structured low-rank matrix sensing problem by exploiting Kronecker product expansion and the structures of the ULA/USPA. The formulated problem has a reduced dimensionality and is solved by using a low-complexity GCG-Alt algorithm. The computational complexity analysis and numerical results suggest that our proposed method is effective in estimating the mmWave channel covariance matrix.