Channel Estimation and Hybrid Precoding for Distributed Phased Arrays Based MIMO Wireless Communications

TL;DR

This paper introduces a novel DPA-MIMO architecture enabling efficient channel estimation and hybrid precoding in mmWave systems, leveraging sparse recovery algorithms and sub-array grouping to improve performance and reduce training overhead.

Contribution

It proposes a coordinated pilot design, two sparse channel recovery algorithms, and a SIC-based hybrid precoding method tailored for DPA-MIMO systems, advancing the state-of-the-art in reconfigurable mmWave MIMO.

Findings

The proposed estimators effectively exploit channel sparsity.

Hybrid precoding achieves good performance with low complexity.

Reduced training overhead compared to traditional methods.

Abstract

Distributed phased arrays based multiple-input multiple-output (DPA-MIMO) is a newly introduced architecture that enables both spatial multiplexing and beamforming while facilitating highly reconfigurable hardware implementation in millimeter-wave (mmWave) frequency bands. With a DPA-MIMO system, we focus on channel state information (CSI) acquisition and hybrid precoding. As benefited from a coordinated and open-loop pilot beam pattern design, all the sub-arrays can perform channel sounding with less training overhead compared with the traditional orthogonal operation of each sub-array. Furthermore, two sparse channel recovery algorithms, known as joint orthogonal matching pursuit (JOMP) and joint sparse Bayesian learning with $\ell_2$ reweighting (JSBL-$\ell_2$), are proposed to exploit the hidden structured sparsity in the beam-domain channel vector. Finally, successive interference cancellation (SIC) based hybrid precoding through sub-array grouping is illustrated for the DPA-MIMO system, which decomposes the joint sub-array RF beamformer design into an interactive per-sub-array-group handle. Simulation results show that the proposed two channel estimators fully take advantage of the partial coupling characteristic of DPA-MIMO channels to perform channel recovery, and the proposed hybrid precoding algorithm is suitable for such array-of-sub-arrays architecture with satisfactory performance and low complexity.

Equations

$${\mathbf{H}}_{m,n}=\sqrt{\frac{N_t^{\mathrm{sub}}{N}_r^{\mathrm{sub}}}{L_{m,n}}}\sum _{i=0}^{L_{m,n}-1}{\alpha }_{m,n}^{\left(i\right)}{\mathbf{a}}_r\left({\vartheta }_{m,n}^{\left(i\right)}\right){\mathbf{a}}_t^H\left({\psi }_{m,n}^{\left(i\right)}\right),$$

$$\mathbf{G}=\left[\begin{array}{ccc}{\mathbf{G}}_{1,1}& \cdots & {\mathbf{G}}_{1,M_t}\\ ⋮& \ddots & ⋮\\ {\mathbf{G}}_{M_r,1}& \cdots & {\mathbf{G}}_{M_r,M_t}\end{array}\right].$$

$${\Omega }_c\triangleq \bigcap _{m=1}^{M_r}\bigcap _{n=1}^{M_t}\mathrm{supp}\left\{{\mathbf{G}}_{m,n}\right\}.$$

$${\Omega }_{m,n}\triangleq \mathrm{supp}\left\{{\mathbf{G}}_{m,n}\right\}\backslash {\Omega }_c.$$

$$y_{q,p}={\mathbf{w}}_q^H\mathbf{H}{\mathbf{f}}_p{x}_p+{\mathbf{w}}_q^H{\tilde{\mathbf{z}}}_p,$$

$${\mathbf{y}}_p={\mathbf{W}}^H\mathbf{H}{\mathbf{f}}_p{x}_p+{\mathbf{z}}_p,$$

$$\mathbf{Y}={\mathbf{W}}^H\mathbf{HF}{\mathbf{X}}_{\mathrm{p}}+\mathbf{Z},$$

$$\mathbf{Y}=\sqrt{P_{\mathrm{p}}}{\mathbf{W}}_{\mathrm{D}}^H{\mathbf{W}}_{\mathrm{R}}^H\mathbf{H}{\mathbf{F}}_{\mathrm{R}}{\mathbf{F}}_{\mathrm{D}}+\mathbf{Z},$$

$$\mathbf{F}=$$

$$=$$

$$\cdot \underset{{\mathbf{F}}_{\text{D}}}{\underbrace{\text{blkdiag}\left\{{\mathbf{F}}_{\text{D},1},\cdots ,{\mathbf{F}}_{\text{D},b_1},\cdots ,{\mathbf{F}}_{\text{D},N_t^{\text{sub}}}\right\}}},$$

$${\mathbf{W}}_{\text{R}}$$

$${\mathbf{W}}_{\text{D}}$$

$$\mathbf{y}$$

$$\stackrel{\left(b\right)}{=}\mathbf{Q}\mathrm{vec}\left(\mathbf{G}\right)+\mathbf{z},$$

$$\mathbf{Q}\triangleq \left({\mathbf{F}}_{\mathrm{D}}^T{\mathbf{F}}_{\mathrm{R}}^T{\mathbf{A}}_t^{\ast }\right)\otimes \left({\mathbf{W}}_{\mathrm{D}}^H{\mathbf{W}}_{\mathrm{R}}^H{\mathbf{A}}_r\right).$$

$$\mathbf{x}\triangleq$$

$$=$$

$$\displaystyle\underbrace{\mathrm{vec}^{T}\left(\mathbf{G}_{1,M_{t}}\right)}_{\left(M_{r}\left(M_{t}-1\right)+1\right)\mathrm{th}\;\mathrm{block}},\cdots,\underbrace{\mathrm{vec}^{T}\left(\mathbf{G}_{M_{r},M_{t}}\right)}_{M_{t}M_{r}\mathrm{th}\;\mathrm{block}}\Big{]}^{T},$$

$$\underset{\mathbf{x}}{\min }$$

$$\displaystyle\mathbf{x}\;\textrm{satisfies the joint sparsity model in Definition~{}{\ref{def_dpa_mimo_chl}}}.$$

$${\mu }^{\mathrm{tot}}\left(\mathbf{\Phi }\right)$$

$$\stackrel{\left(b\right)}{\le }{\mu }^{\mathrm{tot}}\left({\mathbf{F}}_{\mathrm{D}}^T{\mathbf{F}}_{\mathrm{R}}^T{\mathbf{A}}_t^{\ast }\right)\cdot {\mu }^{\mathrm{tot}}\left({\mathbf{W}}_{\mathrm{D}}^H{\mathbf{W}}_{\mathrm{R}}^H{\mathbf{A}}_r\right),$$

$$\begin{array}{rl}{\mu }^{\mathrm{tot}}\left({\mathbf{F}}_{\mathrm{D}}^T{\mathbf{F}}_{\mathrm{R}}^T{\mathbf{A}}_t^{\ast }\right)& \stackrel{\left(a\right)}{=}{\Vert {\left({\mathbf{F}}_{\mathrm{D}}^T{\mathbf{F}}_{\mathrm{R}}^T{\mathbf{A}}_t^{\ast }\right)}^H{\mathbf{F}}_{\mathrm{D}}^T{\mathbf{F}}_{\mathrm{R}}^T{\mathbf{A}}_t^{\ast }-{\mathbf{I}}_{N_t^{\mathrm{tot}}}\Vert }_F^2\\ & \stackrel{\left(b\right)}{=}{\Vert {\mathbf{F}}_{\mathrm{D}}^T{\mathbf{F}}_{\mathrm{R}}^T{\mathbf{A}}_t^{\ast }{\mathbf{A}}_t^T{\mathbf{F}}_{\mathrm{R}}^{\ast }{\mathbf{F}}_{\mathrm{D}}^{\ast }-{\mathbf{I}}_{N_t^{\mathrm{beam}}}\Vert }_F^2+d_0\\ & \stackrel{\left(c\right)}{=}{\Vert {\mathbf{F}}_{\mathrm{D}}^T{\mathbf{F}}_{\mathrm{D}}^{\ast }-{\mathbf{I}}_{N_t^{\mathrm{beam}}}\Vert }_F^2+d_0,\end{array}$$

$$\underset{{\Vert {\left[{\mathbf{F}}_{\text{D},b_1}\right]}_{:,m}\Vert }_2^2=1}{\min }$$

$$\mathbf{F}^{\star}_{\mathrm{D},b_{1}}=\bar{\mathbf{U}}_{t}\Big{[}\mathbf{I}_{\frac{N_{t}^{\mathrm{beam}}}{N_{t}^{\mathrm{sub}}}},\mathbf{0}_{\frac{N_{t}^{\mathrm{beam}}}{N_{t}^{\mathrm{sub}}},M_{t}-\frac{N_{t}^{\mathrm{beam}}}{N_{t}^{\mathrm{sub}}}}\Big{]}^{T}\bar{\mathbf{V}}_{t}^{H},$$

$$\left[{\left({\mathbf{f}}_{\mathrm{R},1}^{\left[m\right]}\right)}^{\star },\cdots ,{\left({\mathbf{f}}_{\mathrm{R},N_t^{\mathrm{sub}}}^{\left[m\right]}\right)}^{\star }\right]=\mathrm{circshift}\left({\mathbf{F}}_{N_t^{\mathrm{sub}}},m-1\right),$$

$$\mathbf{W}^{\star}_{\mathrm{D},b_{2}}=\bar{\mathbf{U}}_{r}\Big{[}\mathbf{I}_{\frac{N_{r}^{\mathrm{beam}}}{N_{r}^{\mathrm{sub}}}},\mathbf{0}_{\frac{N_{r}^{\mathrm{beam}}}{N_{r}^{\mathrm{sub}}},M_{r}-\frac{N_{r}^{\mathrm{beam}}}{N_{r}^{\mathrm{sub}}}}\Big{]}^{T}\bar{\mathbf{V}}_{r}^{H},$$

$$\left[{\left({\mathbf{w}}_{\mathrm{R},1}^{\left[n\right]}\right)}^{\star },\cdots ,{\left({\mathbf{w}}_{\mathrm{R},N_r^{\mathrm{sub}}}^{\left[n\right]}\right)}^{\star }\right]=\mathrm{circshift}\left({\mathbf{F}}_{N_r^{\mathrm{sub}}},n-1\right),$$

$${\mathbf{\gamma }}_{\left(II\right)}$$

$$=\text{arg}\underset{\mathbf{\gamma }⪰0}{\min }{\mathbf{y}}^H{\mathbf{\Sigma }}_y^{-1}\mathbf{y}+\ln \left|{\mathbf{\Sigma }}_y\right|,$$

$$\begin{gathered} {\mathbf{x}}_{\left(II\right)}=\mathbb{E}\left[\mathbf{x}|\mathbf{y};{\mathbf{\gamma }}_{\left(II\right)}\right] \\ ={\mathbf{\Gamma }}_{\left(II\right)}{\mathbf{\Phi }}^H{\left(\lambda \mathbf{I}+\mathbf{\Phi }{\mathbf{\Gamma }}_{\left(II\right)}{\mathbf{\Phi }}^H\right)}^{-1}\mathbf{y}. \end{gathered}$$

$$\underset{\mathbf{C},\mathbf{S}}{\min }{\Vert \mathbf{y}-\mathbf{\Phi }\text{vec}\left(\mathbf{C}+\mathbf{S}\right)\Vert }_2^2+{\beta }_1\sum _{k=1}^K{\Vert {\mathbf{s}}_k\Vert }_1+{\beta }_2{\Vert \mathbf{C}\Vert }_{1,2}$$

$${\Vert {\mathbf{x}}_k\Vert }_1$$

$${\Vert \mathbf{X}\Vert }_{1,2}$$

Full text

I Introduction

Driven by the tremendous growth in demand for wireless data, many new technologies have been proposed for fifth generation cellular communications (5G) to enable orders of magnitude increases in the network capacity [1, 2]. In the physical layer, the exploration of new spectrum in the so-called 5G upper bands, for example, from 6 GHz up to 100 GHz, including the millimeter-wave (mmWave) frequencies, has made multi-gigabit-per-second wireless communications more promising and feasible [3]. However, applying mmWave communications to commercial cellular networks is very challenging mainly due to, first, much higher propagation losses compared with those at lower microwave frequencies; second, strict constraints on hardware designs and implementations which include but are not limited to, antenna performance and dimension, circuits and systems integration challenges [4], power consumption and power supply, form factor (particularly critical for a mobile handset device), etc., according to [5]. Fortunately, massive antenna elements working at mmWave bands can be accommodated into a limited hardware area due to shorter wavelength, facilitating large beamforming gain for combating large pathloss and establishing stable links with reasonable signal-to-noise ratio (SNR) values.

To enable massive multiple-input multiple-output (MIMO) communication with less ratio frequency (RF) chains, a hybrid analog-digital solution [6, 7] was proposed. In this architecture, the signal processing in conventional MIMO is divided into low-dimensional digital beamforming, and high-dimensional analog beamforming that is implemented with low-cost phase shifters. This hybrid transceiver topology is further categorized into the fully-connected and sub-array based structures in terms of how RF chains are mapped to antennas. In a fully-connected structure, each RF chain enables full array gain through individual connection to all antennas [8]; while for the latter structure, each RF chain is only connected to partial antennas, which reduces complexity at the penalty of degrading beamforming gain [9]. In practice, the array-of-sub-array structure has drawn great attention due to its low-complexity implementation, high energy efficiency, and flexible configurations [10, 11, 12].

Recently, a distributed phased arrays based MIMO (DPA-MIMO) architecture which can be easily applied to both base station (BS) and user equipment (UE) designs, has been proposed for practical system and hardware designs [5]. One of its key characteristics distinguishing from the traditional array-of-sub-array structure is to deploy sub-arrays in separate locations similar to distributed antenna systems [13] with relatively small separation. More importantly, the underlying mechanism of DPA-MIMO requires the sub-arrays to separate from each other for several reasons [4, 5, 14]: 1) The practical coupling effects that cause spatial interference (such as the unavoidable side lobes) can be mitigated to guarantee the isolation and independence of each sub-array; 2) The heat dissipation capability and thermal performance can be enhanced by separating sub-arrays with some distance; 3) In particular, at the UE end, separating sub-arrays is critical to overcome the human body (hand) blockage issue. This highly reconfigurable architecture facilitates the multi-beam multi-stream based 5G system and hardware designs under some realistic constraints and resources limitation, which enables appealing advanced features for both academic research and industrial applications [15, 16, 17, 18].

Owing to the geographically separated sub-array structure, the DPA-MIMO based channel modeling requires an in-depth exploration. The basic idea has been inspired by the channel measurements for outdoor BS composed of a very large array [19]. Instead, this extra-large MIMO channel cannot be characterized as a wide sense stationary (WSS) model in the spatial domain. Following the analysis in [19, 20], this spatially non-WSS channel feature results from that, some clusters are only visible to a part of the large array due to some practical factors, such as the cluster sizes, the array aperture, and the spacing between these clusters and the array. Furthermore, this spatial non-stationarity in massive MIMO channels was visualized by introducing the concepts of partially visible clusters and wholly visible clusters which are categorized according to their visibility regions (VRs) [21, 22]. A recent measurement campaign has further verified this spatially non-stationary characteristic over a $40\times 40$ planar RX array at frequencies from 13 to 17 GHz [23]. Based on the above discussion, increasing the spacing between adjacent sub-arrays inevitably triggers an effect that independent scatterers appear over different sub-arrays in a DPA-MIMO system due to their VRs, which has also occurred in indoor THz communications [24]. Meanwhile, as constrained by hardware dimension and power consumption, a reasonable distance between adjacent sub-arrays should be set [5]. As a result, the channels between any transmitter (TX) and receiver (RX) sub-arrays are partial coupling, i.e., both common and local scatterers exist in a DPA-MIMO transmission environment. Due to the inherent sparsity of mmWave channels, we consider modeling the DPA-MIMO channel using the virtual angular domain representation [25]. In this way, the channel coupling relationship can be conveniently characterized as structured sparsity in the beam domain [26].

Channel state information (CSI) acquisition in mmWave systems is challenging due to high dimensional channels, and low point-to-point SNR before beamforming. This makes conventional channel estimators such as the least square (LS) approach infeasible. Exploiting the mmWave channel sparsity to reduce the training overhead is expected to address this challenge [27, 26, 28]. By appropriately choosing the RF and digital precoder matrices, the work in [29] developed a spatial grid based orthogonal matching pursuit (OMP) method to estimate the channel of hybrid MIMO systems. A sparse Bayesian learning (SBL) based channel estimator in [30] was further demonstrated to achieve better performance than the greedy method that is sensitive to the choice of the dictionary matrix in [29]. By leveraging the common sparsity over multiple measurement vectors (MMV), the authors in [31] applied an SBL based approach to mmWave hybrid MIMO systems for accurate support detection. This type of block sparsity, e.g., the common angle-domain channel sparsity over all pilot subcarriers, has also been observed and exploited in [26, 32] to enhance broadband channel estimation when orthogonal frequency-division multiplexing (OFDM) modulation is utilized. However, different channels between any pair of the TX and RX sub-arrays are jointly correlated due to the shared common scatterers in DPA-MIMO systems. Therefore, it is highly desirable to exploit both common and innovation sparsity to decrease the pilot overhead and improve the accuracy of the DPA-MIMO channel estimate.

With the obtained channel, hybrid precoding should be performed to facilitate directional data transmission in DPA-MIMO systems. Some works have been devoted to hybrid precoding in energy-efficient sub-array architectures [11, 12, 33, 9]. In [33], Yu et al. proposed an alternating optimization based method to minimize the Euclidean distance between the fully digital precoder and the hybrid precoders. Therein, during each iteration, the optimal digital precoder is obtained by solving a semidefinite relaxation (SDR) problem with a heavy computational burden while the optimal RF precoder has a closed-form expression. The work in [9] designed hybrid precoders by creatively introducing the mechanism of successive interference cancellation (SIC) in multi-user detection. Thus, a total achievable rate maximization problem, with nonconvex constant amplitude constraints of phase shifters, is decomposed into a series of subrate optimization problems each of which handles one sub-array. However, this algorithm only focuses on the TX design with a simple diagonal baseband precoding matrix under the case that the number of data streams equals that of RF chains. Inspired by [9], we extend the SIC idea to hybrid precoding in DPA-MIMO systems.

In this paper, we consider cooperative multi-sub-array based channel estimation and hybrid precoding for DPA-MIMO systems. The main contributions are summarized below.

• •

We exploit joint channel sparsity among distributed sub-arrays. The inter-sub-array coupling channel model motivates us to take advantage of a multi-sub-array coordinated channel sounding strategy which undoubtedly decreases the training overhead. Based on this strategy, we formulate the DPA-MIMO channel estimation problem as a structured single measurement vector (SMV) recovery problem in compressed sensing (CS) [34]. Instead of using traditional random pilots at the cost of complicated RF hardware, we design deterministic pilot beam patterns by minimizing the total coherence of the equivalent measurement matrix [35], which have successful applications in fully-connected hybrid MIMO [29, 31] and lens arrays [36].

• •

We propose two customized algorithms to find the optimal sparse channel vector with equi-length structured blocks each of which has both the common and innovation supports. Inspired by [37], we divide the proposed joint orthogonal matching pursuit (JOMP) algorithm into two intuitive parts, namely the common support identification and the innovation support identification. Another one is called the joint SBL (JSBL)-$\ell_{2}$ algorithm which adapts the SBL framework [38, 39, 40] to the structured DPA-MIMO channel estimation problem by capitalizing on a dual-space transform.

• •

We propose a low-complexity SIC-based hybrid precoding scheme through sub-array grouping for the array-of-sub-arrays architecture. For the design of the RF beamformers, the idea of SIC is used to decompose the spectral efficiency (SE) maximization problem into several subproblems each of which is only related to one group of sub-arrays, thereby facilitating efficient handling of the sub-arrays group by group.

Organization: The rest of the paper is organized as follows. Section II introduces a jointly sparse DPA-MIMO channel model and a scheme of cooperative multi-sub-array beam training. Section III formulates the DPA-MIMO channel estimation problem and proposes the sub-array based pilot beam pattern design. Section IV presents two channel recovery algorithms based on the structured channel sparsity. Section V specifies the SIC-based hybrid precoding design through sub-array grouping, with simulation results analyzed in Section VI. Finally, Section VII concludes this paper.

Notation: bold uppercase $\mathbf{A}$ (bold lowercase $\mathbf{a}$) denotes a matrix (a vector). We denote $\left[\mathbf{A}\right]_{i,j}$ and $\left[\mathbf{A}\right]_{:,j}$ as its $\left(i,j\right)$th element and $j$th column, respectively. ${\mathrm{vec}}\left(\cdot\right)$ stacks the columns of a matrix into a tall vector, and ${\mathrm{Tr}}\left\{\cdot\right\}$ stands for the matrix trace operation. ${{\mathbf{I}}_{N}}$ and ${{\mathbf{0}}_{M,N}}$ denote the $N\times N$ dimensional identity matrix and the $M\times N$ dimensional all-zero matrix, respectively. ${\left(\cdot\right)^{H}}$, ${\left(\cdot\right)^{T}}$ , ${\left(\cdot\right)^{*}}$, ${\left(\cdot\right)^{-1}}$ and ${\left(\cdot\right)^{\dagger}}$ stand for the conjugate transpose, transpose, conjugate, inverse and pseudo-inverse, respectively. ${\mathrm{diag}}\left\{{\mathbf{a}}\right\}$, ${\mathrm{diag}}\left\{{\mathbf{A}}\right\}$ and ${\mathrm{blkdiag}}\left\{{\mathbf{A}_{1}},\cdots,{\mathbf{A}_{N}}\right\}$ represent a diagonal matrix with $\mathbf{a}$ along its main diagonal, a vector constructed by the main diagonal of the matrix $\mathbf{A}$, a block diagonal matrix , respectively. $\otimes$ denotes the Kronecker product of two matrices. $\ell_{0}$ , $\ell_{1}$ and $\ell_{2}$ norm of vectors are denoted by $\left\|\cdot\right\|_{0}$, $\left\|\cdot\right\|_{1}$ and $\left\|\cdot\right\|_{2}$, respectively. ${\left\|\mathbf{A}\right\|_{\mathrm{F}}}$ denotes the Frobenius norm and the mixed $\ell_{1}/\ell_{2}$ norm is defined as ${\left\|{\mathbf{A}}\right\|_{1,2}}\triangleq\sum\nolimits_{i}{\sqrt{\sum\nolimits_{j}{{{|{{{\left[{\mathbf{A}}\right]}_{i,j}}}|}^{2}}}}}$. ${\mathcal{CN}}\left(\bm{\mu},\mathbf{R}\right)$ denotes the complex Gaussian distribution with mean $\bm{\mu}$ and covariance matrix $\mathbf{R}$. ${\mathbb{E}}\left\{\cdot\right\}$ is the expectation operator. Finally, $\backslash$ denotes the set subtraction operation.

II System Model

In this section, we present the jointly sparse DPA-MIMO channel model and the cooperative multi-sub-array beam training design in the DPA-MIMO system.

II-A Joint Channel Sparsity Model

Consider a $M_{r}\times M_{t}$ DPA-MIMO system shown in Fig. 1, where a TX with $M_{t}$ sub-arrays communicates $N_{s}$ data streams to a RX with $M_{r}$ subrrays. We denote by $N_{t}^{\rm tot}$ ($N_{r}^{\rm tot}$) the total number of antennas at the TX/RX end. Note that $N_{t}^{\rm tot}=M_{t}N_{t}^{\rm sub}$ and $N_{r}^{\rm tot}=M_{r}N_{r}^{\rm sub}$. Furthermore, we assume that each sub-array is a uniform linear array (ULA)111The ULA is widely applied to the existing cellular UE and BS, e.g., the performance comparison of single-user precoding schemes with multi-user precoding schemes has been conducted based on ULA configuration in the high-impact work [8]. Moreover, [41] studied the efficacy of different beamforming approaches for initial UE discovery in mmWave MIMO systems, also based on a typical example of ULA setup case. On the other hand, admittedly, the uniform planar arrays (UPAs) are believed to play a promising and crucial role at both BS [5] and UE [42] end, particularly for mmWave cellular communications, since they can accommodate many more antenna elements into a sub-array in two-dimensional (2-D) configuration to enable 3-D beamforming. The proposed DPA-MIMO application can be also extended to the UPA setup by using the Kronecker product based codebook [43]. We expect to extend and discuss this topic carefully in the future work., and all the sub-arrays are equally spaced and are arranged in the same axis at both TX and RX ends. Fig. 2 shows that $d_{e}$ is antenna element spacing inside each sub-array and $d_{a}$ defines the inter-sub-array spacing. In order to avoid causing too serious grating lobes and significant channel capacity degradation in the DPA-MIMO system, the antenna spacing and inter-sub-array spacing usually satisfy $d_{e}=0.5\lambda_{c}$ and $d_{a}\geq 1.5\lambda_{c}$ respectively [5]. Furthermore, $d_{a}>>d_{e}$ can be easily satisfied at mmWave frequencies. For example, by adopting the DPA-MIMO architecture on a unmanned aerial vehicles (UAV) where $d_{a}$ is larger than 20 times the free-space wavelength, or 40 times the $d_{e}$, Gbps data-rate communication for multi-user scenarios without interference can be enabled and has been preliminarily verified from field tests [16, 17].

Compared with the rich scattering channel model often used for microwave frequencies, mmWave channels are better characterized by a limited number of scattering clusters [41]. Thus, the $L_{m,n}$-path narrowband channel matrix between the $n$th TX sub-array and the $m$th RX sub-array is formulated as

[TABLE]

where $\alpha_{m,n}^{\left(0\right)}$ is the complex gain of the line-of-sight (LoS) component with $\vartheta_{m,n}^{\left(0\right)}$ and $\psi_{m,n}^{\left(0\right)}$ representing its spatial directions composed of an angle of arrival (AoA) and an angle of departure (AoD), respectively. For $i=1,2,\cdots,L_{m,n}-1$, $\alpha_{m,n}^{\left(i\right)}$ is the complex gain of the $i$th non-LoS (NLoS) component with $\vartheta_{m,n}^{\left(i\right)}$ and $\psi_{m,n}^{\left(i\right)}$ denoting its spatial directions composed of an AoA and an AoD, respectively. Note that some AoDs/AoAs are correlated along different sub-arrays, which will be investigated in the proposed jointly sparse DPA-MIMO channel model defined in Definition 1. The path amplitudes are assumed to be Rayleigh distributed, i.e., $\alpha_{m,n}^{\left(0\right)}\sim\mathcal{CN}\left(0,\sigma_{\mathrm{LoS}}^{2}\right)$ and $\alpha_{m,n}^{\left(i\right)}\sim\mathcal{CN}\left(0,\sigma_{\mathrm{NLoS}}^{2}\right)$, where $\sigma_{\mathrm{LoS}}^{2}$ and $\sigma_{\mathrm{NLoS}}^{2}$ are the variances of the LoS and NLoS path gain, respectively [44]. For an ULA with $N$ antennas, the array response vector is ${\mathbf{a}}\left(\psi\right)\triangleq\frac{1}{{\sqrt{N}}}{\left[{1,{e^{-j\frac{{2\pi{d_{e}}}}{{{\lambda_{c}}}}\psi}},\cdots,{e^{-j\frac{{2\pi{d_{e}}}}{{{\lambda_{c}}}}\left({N-1}\right)\psi}}}\right]^{T}}$ and the spatial direction is defined as $\psi=\cos\theta$ where $\theta$ is the physical direction and $\lambda_{c}$ is the carrier wavelength. We use $\mathbf{a}_{t}\left(\cdot\right)$ and $\mathbf{a}_{r}\left(\cdot\right)$ to denote the array response vectors for the TX and RX sub-arrays, respectively, and define the entire channel matrix between all TX and RX sub-arrays using $\mathbf{H}$.

In a DPA-MIMO system, since the AoAs or AoDs for different sub-arrays are partially overlapped as shown in Fig. 1, the spatial resolution is actually determined by the number of sub-array antennas. Assuming in an extreme case that two sub-arrays have different AoAs, the dimension of the channel should equal the number of sub-array antennas [45]. The physical spatial domain and the beam domain are related through a spatial unitary transform matrix, which contains the array steering vectors of uniformly spaced orthogonal spatial directions covering the entire space, e.g., $\mathbf{U}=\left[\mathbf{a}\left(\bar{\psi}_{1}\right),\cdots,\mathbf{a}\left(\bar{\psi}_{N_{G}}\right)\right]$, where $\bar{\psi}_{g}=\frac{2}{{{N_{G}}}}\left({g-1}\right)-1$ for $g=1,\cdots,N_{G}$ with $N_{G}\in\left\{{N_{t}^{{\text{sub}}},N_{r}^{{\text{sub}}}}\right\}$ [25]. Note that finite-resolution discrete dictionary may cause the power leakage problem. One way to handle this problem is to use a redundant dictionary matrix. Interested readers can refer to [36] for detailed operations. We use $\mathbf{U}_{t}$ ($\mathbf{U}_{r}$) to denote the spatial transform matrix for each TX/RX sub-array. Thus, the beam-domain channel matrix between the $n$th TX sub-array and the $m$th RX sub-array can be represented as $\mathbf{H}_{m,n}=\mathbf{U}_{r}\mathbf{G}_{m,n}\mathbf{U}_{t}^{H}$. Subsequently, we can express the relationship between the entire spatial channel and the entire beam-domain channel for the DPA-MIMO system as ${\bf{H}}=\mathbf{A}_{r}{\bf{G}}\mathbf{A}_{t}^{H}$, where $\mathbf{A}_{t}={{{\bf{I}}_{{M_{t}}}}\otimes{\bf{U}}_{t}}$ and $\mathbf{A}_{r}={{{\bf{I}}_{{M_{r}}}}\otimes{{\bf{U}}_{r}}}$ constitute the beam-domain transform matrix for the TX and the RX, respectively, and $\mathbf{G}$ is the entire beam-domain channel matrix that has the following form

[TABLE]

Since spatially correlated mmWave channels are expected to have limited scattering, the beam-domain subchannels $\mathbf{G}_{m,n}$’s are sparse by neglecting the subtle grid quantization error of the AoAs and AoDs [25, 46]. Furthermore, the entire beam-domain channel $\mathbf{G}$ is composed of all the subchannels $\mathbf{G}_{m,n}$’s leading to its sparsity equal to the sum of all the components, which is vividly shown in the left part of Fig. 3. In order to decrease the training overhead in the DPA-MIMO system working on mmWave bands, it is more efficient to estimate the sparse beam-domain channel $\mathbf{G}$ instead of the original physical channel $\mathbf{H}$ [44, 2].

As illustrated in Fig. 1, there exist both the common and local scatterers between any pair of the TX and RX sub-arrays due to large sub-array spacing. In this figure, the VR of the common scatterer covers all the sub-arrays while the VR of any local scatterer only illuminates certain sub-array at the TX/RX end. This result is analogous to that of the newly studied extra-large antenna arrays [20, 22, 47]. In the DPA-MIMO architecture, the structure composed of all sub-arrays at the TX/RX end can be treated as another special version of an extra-large ULA where segments of antenna elements are removed at uniform intervals. Therefore, the spatial non-stationarity of DPA-MIMO means that different scatterers may be observed by different sub-arrays, which is parallel to the spatially non-WSS assumption made on the large array [48]. Due to the inherent channel sparsity at mmWave frequencies, it is more convenient to investigate this kind of spatially non-WSS channels from the perspective of virtual angular domain [37]. Therefore, we conclude the following assumption on the beam-domain channel matrices in the DPA-MIMO system:

Definition 1** (Jointly Sparse DPA-MIMO Channel).**

The channel matrices $\left\{\mathbf{G}_{m,n}\right\}$ have the following properties:

• •

Common sparsity* due to both the LoS path and common scattering: Denote $\mathrm{supp}\left\{\mathbf{A}\right\}$ as the index set of non-zero entries of the matrix $\mathbf{A}$. Then, $\left\{\mathbf{G}_{m,n}\right\}$ are simultaneously sparse. Different $\left\{\mathbf{G}_{m,n}\right\}$ share a common support, i.e., the index set $\Omega_{c}$ which satisfies*

[TABLE]

• •

Innovation sparsity* due to local scattering: There exist unique components for $\mathbf{G}_{m,n}$, i.e.,*

[TABLE]

Remark 1**.**

Common sparsity consists of the paths and sub-rays that lead to the same AoA and AoD for different sub-arrays. The exact components of the common sparsity set depend on the TX and RX positions, sub-array spacings, and the relative locations and orientations of the scattering surfaces. Moreover, the property of common sparsity is based on the prerequisite that all the sub-arrays at the TX (RX) are installed in parallel, e.g., all the sub-arrays composed of ULAs at the TX (RX) are deployed in the same line [49] or all the sub-arrays composed of UPAs at the TX (RX) are parallelly placed on the same plane [11, 5]. Some examples of demonstrations are shown in Fig. 5 and Fig. 13 of [5]. Furthermore, take mmWave UPAs as an example, the spacing between sub-arrays is proposed to be larger than 1.5 times the free-space wavelength (one wavelength is 10.7 mm when the carrier frequency is 28 GHz) [5], in order to avoid electromagnetic mutual coupling [50]. Thus, several free-space wavelengths that are equal to several centimeters are still relatively small compared to the transmission distance of the LoS communication path and the transmission distance between any common scatterer and any sub-array. As a result, the LoS path or any path caused by common scattering produces almost same AoA and AoD upon all the sub-arrays based on the far-field approximation [51], which leads to a common support in the beam domain.

It is observed from Definition 1 that the DPA-MIMO channel sparsity support is parametrized by the set $\left\{\Omega_{c},\left\{\Omega_{m,n}\right\}\right\}$, where $\Omega_{m,n}$ and $\Omega_{c}$ determine the innovation sparsity support and the shared common sparsity support, respectively. In the most common case of $\Omega_{m,n}=\varnothing$ for $\forall m,n$, all sub-arrays at the TX/RX have the same AoDs (AoAs) [49, 11, 9]. In another case of $\Omega_{c}=\varnothing$, independent scatterers are present for each sub-array [24]. Furthermore, Fig. 3 is given for a better understanding of the structured sparsity in the DPA-MIMO channel. In this example, we assume a $2\times 2$ DPA-MIMO system with $N_{t}^{\mathrm{sub}}=N_{r}^{\mathrm{sub}}=4$, and there are the LoS path and one common cluster among all sub-arrays and one local cluster within each transceiver sub-array pair. By applying the vectorized operation to each beam-domain channel matrix $\mathbf{G}_{m,n}$ and putting the obtained channel vectors together, we can formulate a new structured sparse matrix $\left[\mathrm{vec}\left(\mathbf{G}_{1,1}\right),\mathrm{vec}\left(\mathbf{G}_{2,1}\right),\mathrm{vec}\left(\mathbf{G}_{1,2}\right),\mathrm{vec}\left(\mathbf{G}_{2,2}\right)\right]$ that is row-sparse plus element-sparse. In the following, for both the LoS path and the shared paths formed by the common clusters, we call them the common paths indiscriminately.

II-B Cooperative Multi-Sub-Array Beam Training

Spatial scanning based beam training approaches are widely adopted for mmWave channel estimation due to its simplicity and high performance [52, 53, 54]. However, for the non-cooperative spatial scanning based training process, each sub-array at the TX should individually spend $N_{t}^{\mathrm{sub}}N_{r}^{\mathrm{sub}}$ training beams defined in discrete Fourier transform (DFT) based RF codebooks [11, 55]. In order to decrease the probing overhead, we adopt a cooperative multi-sub-array beam training scheme222The hierarchical codebook based scheme [46, 41, 56, 45] only requires $SL^{2}\left\lceil SL/M_{r}\right\rceil\log_{S}\left(N^{\mathrm{tot}}/L\right)$ TX training beams, where $L_{m,n}=L$ for $\forall m,n$, $N_{t}^{\mathrm{tot}}=N_{r}^{\mathrm{tot}}=N^{\mathrm{tot}}$ and $S$ is a design parameter that is usually set to be 2. However, this low overhead scheme is not suitable for the partial coupling DPA-MIMO channel characterized in Definition 1.. For CSI acquisition, the TX uses $N_{t}^{\mathrm{beam}}$ ($N_{t}^{\mathrm{beam}}\leq N_{t}^{\mathrm{tot}}$) pilot beam patterns denoted as \big{\{}\mathbf{f}_{p}\in\mathbb{C}^{N_{t}^{\mathrm{tot}}\times 1}:{\left\|\mathbf{f}_{p}\right\|_{2}^{2}=1},{p=1,\cdots,N_{t}^{\mathrm{beam}}}\big{\}}, and the RX adopts $N_{r}^{\mathrm{beam}}$ ($N_{r}^{\mathrm{beam}}\leq N_{r}^{\mathrm{tot}}$) pilot beam patterns denoted as \big{\{}\mathbf{w}_{q}\in\mathbb{C}^{N_{r}^{\mathrm{tot}}\times 1}:{\left\|\mathbf{w}_{q}\right\|_{2}^{2}=1},{q=1,\cdots,N_{r}^{\mathrm{beam}}}\big{\}}.

During the training period, the TX successively sends its training beam patterns $\left\{\mathbf{f}_{p}\right\}$ which are received by the RX through its beam patterns $\left\{\mathbf{w}_{q}\right\}$. The $\left({q,p}\right)$th received sample for the $p$th TX beam pattern is given by [31]

[TABLE]

where $x_{p}$ is the transmitted pilot symbol and $\tilde{\mathbf{z}}_{p}\in\mathbb{C}^{N_{r}^{\mathrm{tot}}\times 1}$ is a noise vector with $\mathcal{CN}\left(\mathbf{0},\sigma_{\mathrm{z}}^{2}\mathbf{I}_{N_{r}^{\mathrm{tot}}}\right)$. Collecting $y_{{q},p}$ for ${q}\in\left\{1,\cdots,N_{r}^{\mathrm{beam}}\right\}$, we have $\mathbf{y}_{p}\in\mathbb{C}^{N_{r}^{\mathrm{beam}}\times 1}$ given by

[TABLE]

where $\mathbf{W}=\left[\mathbf{w}_{1},\cdots,\mathbf{w}_{{N_{r}^{\mathrm{beam}}}}\right]\in\mathbb{C}^{N_{r}^{\mathrm{tot}}\times N_{r}^{\mathrm{beam}}}$ and ${\mathbf{z}}_{p}={{\mathbf{W}}^{H}}{\tilde{\mathbf{z}}_{p}}\in\mathbb{C}^{N_{r}^{\mathrm{beam}}\times 1}$. To represent the received signals for all TX beam patterns, we collect $\mathbf{y}_{p}$ for $p\in\left\{1,\cdots,N_{t}^{\mathrm{beam}}\right\}$ yielding

[TABLE]

where the concatenated received signal matrix $\mathbf{Y}=\left[\mathbf{y}_{1},\cdots,\mathbf{y}_{N_{t}^{\mathrm{beam}}}\right]\in\mathbb{C}^{N_{r}^{\mathrm{beam}}\times N_{t}^{\mathrm{beam}}}$, the complete TX processing matrix $\mathbf{F}=\left[\mathbf{f}_{1},\cdots,\mathbf{f}_{N_{t}^{\mathrm{beam}}}\right]\in\mathbb{C}^{N_{t}^{\mathrm{tot}}\times N_{t}^{\mathrm{beam}}}$ and the concatenated noise matrix $\mathbf{Z}=\left[{{{\mathbf{z}}_{1}},\cdots,{{\mathbf{z}}_{N_{t}^{{\text{beam}}}}}}\right]\in\mathbb{C}^{N_{r}^{\mathrm{beam}}\times N_{t}^{\mathrm{beam}}}$. In general, we choose the pilot matrix $\mathbf{X}_{\mathrm{p}}=\sqrt{P_{\mathrm{p}}}\mathbf{I}_{N_{t}^{\mathrm{beam}}}$ where $P_{\mathrm{p}}$ is the pilot power per transmission.

In our array-of-sub-arrays architecture, the TX and RX processing matrices are decomposed as $\mathbf{F}=\mathbf{F}_{\rm{R}}\mathbf{F}_{\rm{D}}$ and $\mathbf{W}=\mathbf{W}_{\rm{R}}\mathbf{W}_{\rm{D}}$. Thus, the multiple measurements in (7) is expressed as

[TABLE]

where $\mathbf{F}_{\mathrm{R}}\in\mathbb{C}^{N_{t}^{\mathrm{tot}}\times N_{t}^{\mathrm{tot}}}$ and $\mathbf{W}_{\mathrm{R}}\in\mathbb{C}^{N_{r}^{\mathrm{tot}}\times N_{r}^{\mathrm{tot}}}$ denote RF beamforming matrices at the TX and the RX, respectively, while $\mathbf{F}_{\mathrm{D}}\in\mathbb{C}^{N_{t}^{\mathrm{tot}}\times N_{t}^{\mathrm{beam}}}$ and $\mathbf{W}_{\mathrm{D}}\in\mathbb{C}^{N_{r}^{\mathrm{tot}}\times N_{r}^{\mathrm{beam}}}$ denote the TX and the RX baseband processing matrices, respectively. At the TX, $\mathbf{F}$ is further partitioned into ${{N_{t}^{{\text{tot}}}}\mathord{\left/{\vphantom{{N_{t}^{{\text{tot}}}}{{M_{t}}}}}\right.\kern-1.2pt}{{M_{t}}}}=N_{t}^{{\text{sub}}}$ blocks, that is,

[TABLE]

where $\mathbf{F}_{\mathrm{R},b_{1}}\in\mathbb{C}^{N_{t}^{\mathrm{tot}}\times M_{t}}$ represents the shared analog processing matrix for all the beam patterns at block $b_{1}$ and each column of $\mathbf{F}_{\mathrm{D},b_{1}}\in\mathbb{C}^{M_{t}\times\frac{N_{t}^{\mathrm{beam}}}{N_{t}^{\mathrm{sub}}}}$ denotes one digital pilot beamforming vector for the corresponding beam pattern at block $b_{1}$ [29, 31, 36]. In the same way, the RX admits the following processing by blocks

[TABLE]

where $\mathbf{W}_{\mathrm{R},b_{2}}\in\mathbb{C}^{N_{r}^{\mathrm{tot}}\times M_{r}}$ and $\mathbf{W}_{\mathrm{D},b_{2}}\in\mathbb{C}^{M_{r}\times\frac{N_{r}^{\mathrm{beam}}}{N_{r}^{\mathrm{sub}}}}$. Moreover, each column of the RF processing matrix is zero except for a continuous block of non-zero entries (consisting of the beamforming weights used on the corresponding sub-array), i.e., ${{\mathbf{F}}_{{\text{R}},{b_{1}}}}={\text{blkdiag}}\left\{{{\mathbf{f}}_{{\text{R}},{b_{1}}}^{\left[1\right]},\cdots,{\mathbf{f}}_{{\text{R}},{b_{1}}}^{\left[n\right]},\cdots,{\mathbf{f}}_{{\text{R}},{b_{1}}}^{\left[{{M_{t}}}\right]}}\right\}$ and ${{\mathbf{W}}_{{\text{R}},{b_{2}}}}={\text{blkdiag}}\left\{{{\mathbf{w}}_{{\text{R}},{b_{2}}}^{\left[1\right]},\cdots{\mathbf{w}}_{{\text{R}},{b_{2}}}^{\left[m\right]},\cdots,{\mathbf{w}}_{{\text{R}},{b_{2}}}^{\left[{{M_{r}}}\right]}}\right\}$, where ${\mathbf{f}_{\mathrm{R},b_{1}}^{\left[n\right]}}\in\mathbb{C}^{N_{t}^{\mathrm{sub}}\times 1}$ and ${\mathbf{w}_{\mathrm{R},b_{2}}^{\left[m\right]}}\in\mathbb{C}^{N_{r}^{\mathrm{sub}}\times 1}$ denote the RF beamforming weights for the $n$th TX sub-array and the $m$th RX sub-array, respectively. Furthermore, the amplitude of each element of ${\mathbf{f}_{\mathrm{R},b_{1}}^{\left[n\right]}}$ and ${\mathbf{w}_{\mathrm{R},b_{2}}^{\left[m\right]}}$ equals $\frac{1}{{\sqrt{N_{t}^{{\text{sub}}}}}}$ and $\frac{1}{{\sqrt{N_{r}^{{\text{sub}}}}}}$, respectively.

III Formulation of Channel Estimation Problem and Pilot Beam Pattern Design

In this section, we first exploit the jointly sparse nature of the DPA-MIMO channel, and formulate its channel estimation problem as a structured sparse vector recovery problem. Then, we propose a deterministic beam training scheme.

III-A Formulation of Channel Estimation Problem

To exploit the sparse nature of the DPA-MIMO channel, it is necessary to vectorize the received signal matrix $\mathbf{Y}$ in (8). After denoting $\frac{{{\text{vec}}\left({\mathbf{Y}}\right)}}{{\sqrt{{P_{\text{p}}}}}}$ by $\mathbf{y}\in\mathbb{C}^{N_{t}^{\mathrm{beam}}N_{r}^{\mathrm{beam}}\times 1}$, we have

[TABLE]

where (a) follows from the equivalent noise vector $\mathbf{z}\triangleq\frac{1}{{\sqrt{{P_{\text{p}}}}}}\left[\mathbf{z}_{1}^{T},\cdots,\mathbf{z}_{N_{t}^{\mathrm{beam}}}^{T}\right]^{T}\in\mathbb{C}^{N_{t}^{\mathrm{beam}}N_{r}^{\mathrm{beam}}\times 1}$ and the properties of Kronecker product, $\mathrm{vec}\left(\mathbf{ABC}\right)=\left(\mathbf{C}^{T}\otimes\mathbf{A}\right)\mathrm{vec}\left(\mathbf{B}\right)$ and $\left(\mathbf{A}\otimes\mathbf{B}\right)^{T}=\mathbf{A}^{T}\otimes\mathbf{B}^{T}$ [57], and (b) follows from $\mathrm{vec}\left(\mathbf{H}\right)=\left(\mathbf{A}_{t}^{\ast}\otimes\mathbf{A}_{r}\right)\mathrm{vec}\left(\mathbf{G}\right)$ and $\left(\mathbf{A}\otimes\mathbf{B}\right)\left(\mathbf{C}\otimes\mathbf{D}\right)=\left(\mathbf{AC}\right)\otimes\left(\mathbf{BD}\right)$. The equivalent sensing matrix $\mathbf{Q}\in\mathbb{C}^{N_{t}^{\mathrm{beam}}N_{r}^{\mathrm{beam}}\times N_{t}^{\mathrm{tot}}N_{r}^{\mathrm{tot}}}$ can be defined as

[TABLE]

The formulation of the vectorized received signal in (11) represents a sparse formulation of the channel estimation problem as $\mathrm{vec}\left(\mathbf{G}\right)$ has only $N_{0}=\left|\Omega_{c}\right|+{\sum\nolimits_{m=1}^{M_{r}}\sum\nolimits_{n=1}^{M_{t}}{\left|\Omega_{m,n}\right|}}$ non-zero elements and $N_{0}\ll N_{t}^{\mathrm{tot}}N_{r}^{\mathrm{tot}}$. This implies that the number of required measurements $N_{t}^{\mathrm{beam}}N_{r}^{\mathrm{beam}}$ to detect the non-zero elements can be much less than $N_{t}^{\mathrm{tot}}N_{r}^{\mathrm{tot}}$. We expect to exploit the hidden joint sparsity in the beam-domain channel to reduce the required training and improve the performance of channel estimation. For convenience, we exchange the order of elements in $\mathrm{vec}\left({\mathbf{G}}\right)$ to get a new vector with $M_{t}M_{r}$ equi-length blocks as

[TABLE]

where the block size is $N_{t}^{\mathrm{sub}}N_{r}^{\mathrm{sub}}$. This manipulation is also shown in Fig. 3. In this case, the corresponding equivalent measurement matrix $\bm{\Phi}\triangleq\mathbf{Q}\bm{\Pi}$ is obtained by exchanging the column order of $\mathbf{Q}$ where $\bm{\Pi}$ is a column permutation matrix, such that $\bm{\Phi}{\mathbf{x}}=\mathbf{Q}\mathrm{vec}\left({\mathbf{G}}\right)$. As a result, the problem of DPA-MIMO channel recovery at the RX can be formulated as

[TABLE]

However, problem (14) is very challenging due to the common and innovation sparsity requirement in the constraint which is quite different from the conventional CS-recovery problem with a simple sparsity ($\ell_{0}$-norm) constraint. In addition, the equivalent measurement matrix has to be carefully designed to guarantee the recovery of the non-zero elements of the vector with high probability by using a small number of measurements.

III-B Open-Loop Pilot Beam Pattern Design

Instead of randomized sensing matrices frequently used for CS-based channel estimation, a deterministic measurement matrix designed by minimizing its total coherence can improve the recovery performance [35, 26]. This strategy has been recently applied to pilot beam pattern design for the fully-connected structure [29, 31]. Due to its excellent performance improvement and zero feedback overhead, this strategy is further applied to channel sounding in the DPA-MIMO system.

Since the total coherence of $\bm{\Phi}$ is defined by $\mu^{\mathrm{tot}}\left(\bm{\Phi}\right)\triangleq\sum_{k=1}^{N_{t}^{\mathrm{tot}}}\sum_{l\neq k}^{N_{t}^{\mathrm{tot}}}\left(\left[\bm{\Phi}\right]_{:,k}^{H}\left[\bm{\Phi}\right]_{:,l}\right)^{2}$, we have

[TABLE]

where $\left(a\right)$ follows from the definition of the total coherence, and $\left(b\right)$ can be derived in a similar way as [29, Lemma 7] and this upper-bound plays an important role in the decoupling of pilot beam design at both ends. Therefore, we can decompose the design problem of minimizing $\mu^{\mathrm{tot}}\left(\bm{\Phi}\right)$ into two separate designs, namely the design of $\mathbf{F}_{\mathrm{D}}$ and $\mathbf{F}_{\mathrm{R}}$ via minimizing $\mu^{\mathrm{tot}}\left({\bf{F}}_{{\rm{D}}}^{T}{\bf{F}}_{{\rm{R}}}^{T}\mathbf{A}_{t}^{*}\right)$ and the design of $\mathbf{W}_{\mathrm{D}}$ and $\mathbf{W}_{\mathrm{R}}$ via minimizing $\mu^{\mathrm{tot}}\left({\bf{W}}_{{\rm{D}}}^{H}{\bf{W}}_{{\rm{R}}}^{H}\mathbf{A}_{r}\right)$.

Sequentially, we can simplify $\mu^{\mathrm{tot}}\left({\bf{F}}_{{\rm{D}}}^{T}{\bf{F}}_{{\rm{R}}}^{T}\mathbf{A}_{t}^{*}\right)$ as

[TABLE]

where $d_{0}=N_{t}^{\mathrm{tot}}-N_{t}^{\mathrm{beam}}$, $\left(a\right)$ is to make the equivalent measurement matrix approximate an identity matrix [35], $\left(b\right)$ results from the relationship between the Frobenius norm and the trace, and $\left(c\right)$ comes from $\mathbf{A}_{t}\mathbf{A}_{t}^{H}=\mathbf{I}_{N_{t}^{\mathrm{tot}}}$ and $\mathbf{F}_{\text{R}}^{T}\mathbf{F}_{\text{R}}^{*}=\mathbf{I}_{N_{t}^{\text{tot}}}$. Thus, we can optimize $\mathbf{F}_{\mathrm{D}}$ by minimizing $\left\|{\bf{F}}_{{\rm{D}}}^{T}{\bf{F}}_{{\rm{D}}}^{*}-\mathbf{I}_{N_{t}^{\mathrm{beam}}}\right\|_{F}^{2}$, which can be further transformed into several parallel subproblems for ${b_{1}}=1,2,\cdots,N_{t}^{{\text{sub}}}$

[TABLE]

Via relaxing the individual power constraints and then using the method of Lagrange multipliers [29, Theorem 2], the optimal baseband precoder of the $b_{1}$th block is given by

[TABLE]

where $\bar{\mathbf{U}}_{t}\in\mathbb{C}^{M_{t}\times M_{t}}$ and $\bar{\mathbf{V}}_{t}\in\mathbb{C}^{\frac{N_{t}^{\mathrm{beam}}}{N_{t}^{\mathrm{sub}}}\times\frac{N_{t}^{\mathrm{beam}}}{N_{t}^{\mathrm{sub}}}}$ are arbitrary unitary matrices, e.g., unitary DFT matrices.

In order to make RF pilot beams cover a full range of AoDs, we choose the unitary DFT matrix as the solution of TX sub-array $m$

[TABLE]

where $\mathbf{F}_{N}$ denotes the $N$-dimensional unitary DFT matrix and $\mathrm{circshift}\left(\mathbf{A},m\right)$ represents moving the columns of a matrix $\mathbf{A}$ to the right for $\left(m\right)$ columns in a circular manner. Different from the previous pilot beam pattern design for the fully-connected structure [29, 31], the dimension of the DFT matrix equals the number of sub-array antennas in a DPA-MIMO system due to the joint sparse channels. In this way, the entire TX array simultaneously probes different spatial directions using RF beams.

Similar operation can be applied to the RX, leading to the optimal combiners as

[TABLE]

[TABLE]

where $\bar{\mathbf{U}}_{r}\in\mathbb{C}^{M_{r}\times M_{r}}$ and $\bar{\mathbf{V}}_{r}\in\mathbb{C}^{\frac{N_{r}^{\mathrm{beam}}}{N_{r}^{\mathrm{sub}}}\times\frac{N_{r}^{\mathrm{beam}}}{N_{r}^{\mathrm{sub}}}}$ are arbitrary unitary matrices. Furthermore, the noise vector after the designed RF and baseband processing remains i.i.d. Gaussian with $\mathcal{CN}\left(\mathbf{0},\frac{{\sigma_{\text{z}}^{2}}}{{{P_{\text{p}}}}}\mathbf{I}_{N_{t}^{\mathrm{beam}}N_{r}^{\mathrm{beam}}}\right)$ without prewhitening, which provides much convenience for the further design of CSI recovery algorithms in the next section.

IV DPA-MIMO Channel Estimation Algorithms

To solve problem (14), we present two customized algorithms, i.e., an OMP based greedy algorithm with low complexity and a SBL inspired algorithm with excellent accuracy. For notational simplicity, we define the number of training overhead $N\triangleq N_{t}^{\mathrm{beam}}N_{r}^{\mathrm{beam}}$, the number of resolvable spatial directions $B\triangleq N_{t}^{\mathrm{sub}}N_{r}^{\mathrm{sub}}$ and $K\triangleq M_{t}M_{r}$, respectively.

IV-A Proposed JOMP Algorithm

The details of the proposed JOMP algorithm are described in Algorithm 1 that is inspired by [37], where not only the per-link channel sparsity but also the joint sparsity structure are exploited to enable a distributed compressive CSI estimation scheme for multi-user massive MIMO systems. First of all, $\delta_{1}$ and $\delta_{2}$ are defined as the predetermined thresholds to control the residual error of different loops within reasonable ranges. In principle, they are related to the noise level $\sigma_{\mathrm{z}}^{2}$ and the measurement vector length $N$. Meanwhile, $T_{\mathrm{1,JOMP}}$ and $T_{\mathrm{2,JOMP}}$ are the maximal numbers of iterations to guarantee the convergence. It is obvious that the ideal values of the two iteration numbers should be set as $T_{\mathrm{1,JOMP}}=\left|\Omega_{c}\right|$ and $T_{\mathrm{2,JOMP}}=\sum_{m=1}^{M_{r}}\sum_{n=1}^{M_{t}}\left|\Omega_{m,n}\right|$. Although it is difficult to acquire the true common and innovation sparsity levels ($\mathbb{S}\triangleq\left\{|\Omega_{c}|,\left\{|\Omega_{m,n}|\right\}\right\}$), their statistic bounds change over a very long time scale and can be easily obtained from long-term stochastic learning and estimation as in [37]. Thus, we can choose a lower bound of the common sparsity $|\Omega_{c}|$ as the iteration number for the common support identification, while rendering an upper bound of the sum of all the individual sparsity as the iteration number for the innovation support identification $|\Omega_{m,n}|$., i.e., $T_{\mathrm{1,JOMP}}\leq\left|\Omega_{c}\right|$ and $T_{\mathrm{2,JOMP}}\geq\sum_{m=1}^{M_{r}}\sum_{n=1}^{M_{t}}\left|\Omega_{m,n}\right|$.

The proposed JOMP algorithm is further divided into two parts, where the first part aims at common support identification, and the second part continues the innovation support identification. Observe that the estimation target $\mathbf{x}$ in (13) has non-zero elements at the same positions of each block. Therefore, motivated by the simultaneous sparse approximation algorithm proposed for MMV problems in [58], we wish to find a group of equi-spaced atoms in the equivalent measurement matrix $\bm{\Phi}$ by maximizing the sum of their absolute correlations with the residual $\mathbf{r}_{1}$. This procedure is done in step 1 and the absolute sum has the equivalent expression as \sum_{k=1}^{K}\left\|\bm{\Phi}_{:,b+\left(k-1\right)B}^{H}\mathbf{r}_{1}\right\|_{2}^{2}=\Big{\|}\left(\bm{\Phi}_{\Omega_{b}^{K}}\right)^{H}\mathbf{r}_{1}\Big{\|}_{2}^{2}. After the common support $\Omega_{c}^{e}$ is detected, the standard OMP method in [59] is used to identify the innovation support $\Omega_{i}^{e}$ as realized from step 5 to 8. Depending on the estimated support index, the LS method is finally used to recovery the channel vector.

IV-B Proposed JSBL-$\ell_{2}$ Algorithm

IV-B1 Introduction to SBL

Consider the classical sparse recovery model without the structured sparsity $\mathbf{y}=\bm{\Phi}\mathbf{x}+\mathbf{z}$, where $\mathbf{z}$ is a noise vector with $\mathcal{CN}\left(\mathbf{0},\lambda\mathbf{I}_{N}\right)$ and $\lambda$ is the known equivalent noise variance. Thus we have the Gaussian likelihood model $p\left(\mathbf{y}|\mathbf{x}\right)=\mathcal{CN}\left(\bm{\Phi}\mathbf{x},\lambda\mathbf{I}_{N}\right)$. Assume the parametrized Gaussian prior $p\left(\mathbf{x}\right)=\mathcal{CN}\left(\mathbf{0},\bm{\Gamma}\right)$, where $\bm{\Gamma}=\mathrm{diag}\left\{\bm{\gamma}\right\}$ with a vector of hyper-parameters $\bm{\gamma}$ governing the prior variances of the elements in $\mathbf{x}$. This latent variable based Gaussian assumption is reasonable from the perspective of variational approximations, which is verified in [38, Sec. V]. For a fixed $\bm{\gamma}$, using the Bayesian rules we can obtain the Gaussian posterior density of $\mathbf{x}$ as $p\left(\mathbf{x}|\mathbf{y}\right)=\mathcal{CN}\left(\bm{\mu}_{x},\bm{\Sigma}_{x}\right)$, where $\bm{\mu}_{x}=\bm{\Gamma}\bm{\Phi}^{H}\left(\lambda\mathbf{I}_{N}+\bm{\Phi}\bm{\Gamma}\bm{\Phi}^{H}\right)^{-1}\mathbf{y}$ and $\bm{\Sigma}_{x}=\bm{\Gamma}-\bm{\Gamma}\bm{\Phi}^{H}\left(\lambda\mathbf{I}_{N}+\bm{\Phi}\bm{\Gamma}\bm{\Phi}^{H}\right)^{-1}\bm{\Phi}\bm{\Gamma}$. The next key task is to estimate the latent variables $\bm{\gamma}$. By treating $\mathbf{x}$ as the hidden variables and integrating them out [38], we obtain the maximum a posterior (MAP) estimate on $\bm{\gamma}$ as

[TABLE]

where the covariance matrix of $\mathbf{y}$ denotes $\bm{\Sigma}_{y}=\lambda\mathbf{I}_{N}+\bm{\Phi}\bm{\Gamma}\bm{\Phi}^{H}$. Once $\bm{\gamma}_{\left(II\right)}$ is obtained, a commonly accepted point estimate for $\mathbf{x}$ naturally emerges as

[TABLE]

This procedure is referred to Type II estimation, also called empirical Bayesian. From (23), it can be observed that a sparse $\bm{\gamma}_{\left(II\right)}$ leads to a corresponding sparse estimate $\mathbf{x}_{\left(II\right)}$. Note that the logarithm term $\ln\left|\bm{\Sigma}_{y}\right|$ in (22) is a concave function with respect to $\bm{\gamma}$ according to [38, Lemma 1], thereby favoring a sparse $\bm{\gamma}$, which further results in a sparse $\mathbf{x}$ through (23). The traditional SBL algorithm assumes independent priors for an estimated signal, which however fails to consider the joint sparsity property in the DPA-MIMO channel.

IV-B2 SBL-Inspired Cost Function

By reshaping the vector $\mathbf{x}$, we define a new matrix $\mathbf{X}\triangleq\left[\mathbf{x}_{1},\cdots,\mathbf{x}_{K}\right]\in\mathbb{C}^{B\times K}$ which is both row-sparse and element-sparse, as shown in Fig. 3. In order to promote such a structure, $\mathbf{X}$ can be viewed as the summation of an element-sparse matrix $\mathbf{S}\triangleq\left[\mathbf{s}_{1},\cdots,\mathbf{s}_{K}\right]\in\mathbb{C}^{B\times K}$ and a row-sparse matrix $\mathbf{C}\triangleq\left[\mathbf{c}_{1},\cdots,\mathbf{c}_{K}\right]\in\mathbb{C}^{B\times K}$ [40]. Furthermore, by using convex approximation, problem (14) can be transformed into the following convex optimization problem333In fact, we just provide one of the intuitive transformations to tackle the challenging problem in (14). As SBL-like algorithms adapt to the measurement matrix and promote a sparse solution, we further address the formulated optimization problem from an SBL perspective. Provided that one finds a special structured prior that contains both common and individual sparsity, turbo-type message passing algorithms may further improve the estimation performance [60, 28]. This meaningful topic requires further investigation.:

[TABLE]

where $\beta_{1}\geq 0$ and $\beta_{2}\geq 0$ are weights regarding element-sparsity and row-sparsity respectively.

In order to promote sparsity of the solution, we transform the cost function of (24) in $x$-space to the SBL-like cost function in $\gamma$-space by using a dual-space view [39], where the following variational representations are used [61] :

[TABLE]

where $x_{bk}\triangleq\left[\mathbf{X}\right]_{b,k}$, $\gamma_{b}^{c}$ and $\gamma_{bk}^{s}$ are scalars, $\bm{\gamma}^{c}\triangleq\left[\gamma_{1}^{c},\cdots,\gamma_{B}^{c}\right]^{T}$ is a vector common to all columns of $\mathbf{X}$, and $\bm{\gamma}^{s}\triangleq\left[\gamma_{1,1}^{s},\cdots,\gamma_{B,1}^{s},\cdots,\gamma_{1,K}^{s},\cdots,\gamma_{B,K}^{s}\right]^{T}$ is a vector with each element corresponding to that of $\mathrm{vec}\left(\mathbf{X}\right)$. By using the identity (its derivation is given in Appendix A)

[TABLE]

where $\mathbf{s}\triangleq\left[\mathbf{s}_{1}^{T},\cdots,\mathbf{s}_{K}^{T}\right]^{T}\in\mathbb{C}^{BK\times 1}$, $\mathbf{c}\triangleq\left[\mathbf{c}_{1}^{T},\cdots,\mathbf{c}_{K}^{T}\right]^{T}\in\mathbb{C}^{BK\times 1}$, $\bm{\Gamma}^{c}\triangleq\mathrm{diag}\left\{\bm{\gamma}^{c}\right\}$, $\bm{\Gamma}^{s}\triangleq\mathrm{diag}\left\{\bm{\gamma}^{s}\right\}$, and

[TABLE]

we can further express the convex cost function of (24) in $\gamma$-space as

[TABLE]

Comparing the data-related term $\mathbf{y}^{H}\left(\bm{\Sigma}^{sc}\right)^{-1}\mathbf{y}$ in (28) and that of the SBL cost function in (22), we can observe that the common component $\bm{\gamma}^{c}$ and the innovation component $\bm{\gamma}^{s}$ interact with each other in a manner like $\bm{\Gamma}=\bm{\Gamma}^{s}+\mathbf{I}_{K}\otimes\bm{\Gamma}^{c}$.

Following the innovative decoupling idea in [40], we can replace the convex penalties in the existing models with the SBL counterpoints to obtain some of the corresponding benefits, even without any formal probabilistic model for this derivation444As a matter of fact, [40] was the first attempt to apply simultaneous SBL approximation to the two multitask structured sparse models with real variables, i.e., row-sparse with embedded element-sparse and row-sparse plus element-sparse. In this paper, we extend the SBL framework to perform the recovery of the complex structured SMV defined in (13) for DPA-MIMO channel estimation.. Therefore, we put forth a new cost function in $\gamma$-space from (28) as

[TABLE]

where

[TABLE]

Since the log-determinant function is concave and non-decreasing, the term $\ln\left|\bm{\Sigma}^{c}\right|$ and the term $\ln\left|\bm{\Sigma}^{s}\right|$ promote a sparse common component $\bm{\gamma}^{c}$ and a sparse innovation component $\bm{\gamma}^{s}$, respectively. Moreover, the weight $\beta$ regarded as a tradeoff between row sparsity and element sparsity in the defined matrix $\mathbf{X}$, should be tuned with training data or given with prior information. Note that this weight depends on the large-scale properties of the scattering environment and changes over a long time scale.

IV-B3 $\ell_{2}$ Reweighting Scheme

Following an extension of the duality space analysis for the basic SBL framework [39], we can transform the cost function of (29) from $\gamma$-space to $x$-space. First, via using the identity (26) and standard determinant identities, we can upper-bound (29) by

[TABLE]

where we define two concave functions, namely

[TABLE]

Due to the duality of concave conjugate functions, we have the following upper bounds given by

[TABLE]

By using (31) and (33), we can then perform block coordinate descent (BCD) optimization over the following approximation with irrelevant terms dropped:

[TABLE]

where $\mathcal{Z}\triangleq\{\mathbf{c},\mathbf{s},\bm{\gamma}^{c},\bm{\gamma}^{s},\mathbf{z}^{c},\mathbf{z}^{s}\}$.

By fixing the other variables, we first calculate the optimal values of $\mathbf{s}$ and $\mathbf{c}$. In this way, the optimization problem (34) is equivalent to problem (26). The optimal solutions are expressed as

[TABLE]

which are also proved in Appendix A.

We then optimize $\mathbf{z}^{c}$ and $\mathbf{z}^{s}$. According to the duality relationship in (33), their optimal values are obtained as

[TABLE]

where the Moore-Penrose pseudo-inverse [57] is used for computation reduction and $\bm{\Xi}\triangleq\left[\mathbf{I}_{B},\cdots,\mathbf{I}_{B}\right]\in\mathbb{C}^{B\times BK}$. Sequentially, the optimal hyperparameters are given by

[TABLE]

with the other variables fixed. Finally, by alternately minimizing and repeatedly updating the upper-bound function (34), we obtain the reweighted algorithm described in Algorithm 2. As the objective function in (34) decreases or keeps unchanged in each iteration, the proposed JSBL-$\ell_{2}$ algorithm can promise a local minimum of problem (34).

V SIC-Based Hybrid Precoding through sub-array grouping

As a result of CSI acquisition in Section IV, we now consider the channel is known at both TX and RX ends. The task of this section is to design the hybrid precoding and combining matrices for the DPA-MIMO system. The processed received signal after combining is given by

[TABLE]

where $\mathbf{s}_{\mathrm{d}}\in\mathbb{C}^{N_{s}\times 1}$ is the data vector such that $\mathbb{E}\left\{\mathbf{s}_{\mathrm{d}}\mathbf{s}_{\mathrm{d}}^{H}\right\}=\frac{{{P_{\text{d}}}}}{{{N_{s}}}}{{\mathbf{I}}_{{N_{s}}}}$, $P_{\mathrm{d}}$ is the average transmitting power, $\mathbf{z}_{\mathrm{d}}$ is a Gaussian noise vector with $\mathcal{CN}\left(\mathbf{0},\sigma_{\mathrm{z}}^{2}\mathbf{I}_{N_{r}^{\mathrm{tot}}}\right)$, $\mathbf{V}_{\mathrm{R}}={\mathrm{blkdiag}}\left\{{\mathbf{v}_{\mathrm{R}}^{\left[1\right]}},\cdots,{\mathbf{v}_{\mathrm{R}}^{\left[M_{t}\right]}}\right\}$ and $\mathbf{U}_{\mathrm{R}}={\mathrm{blkdiag}}\left\{{\mathbf{u}_{\mathrm{R}}^{\left[1\right]}},\cdots,{\mathbf{u}_{\mathrm{R}}^{\left[M_{r}\right]}}\right\}$. Furthermore, we assume the amplitude of each element of ${\mathbf{v}_{\mathrm{R}}^{\left[n\right]}}$ and ${\mathbf{u}_{\mathrm{R}}^{\left[m\right]}}$ equals $\frac{1}{{\sqrt{N_{t}^{{\text{sub}}}}}}$ and $\frac{1}{{\sqrt{N_{r}^{{\text{sub}}}}}}$, respectively. With Gaussian signaling employed at the TX, the instantaneous achievable SE is

[TABLE]

where ${{\mathbf{R}}_{\text{z}}}\triangleq\sigma_{\text{z}}^{2}{\mathbf{U}}_{{\text{D}}}^{H}{\mathbf{U}}_{{\text{R}}}^{H}{{\mathbf{U}}_{{\text{R}}}}{{\mathbf{U}}_{{\text{D}}}}$. Thus we determine the hybrid precoders and combiners by maximizing the SE defined in (39) with a transmitting power constraint $\left\|\mathbf{V}_{\mathrm{R}}\mathbf{V}_{\mathrm{D}}\right\|_{F}^{2}=N_{s}$.

V-A Design of $\mathbf{V}_{\mathrm{D}}$ and $\mathbf{V}_{\mathrm{R}}$

In lieu of maximizing the SE, we design $\mathbf{V}_{\mathrm{R}}$ and $\mathbf{V}_{\mathrm{D}}$ to maximize the mutual information achieved by Gaussian signaling [6]. Then, the hybrid precoder design problem can be written as

[TABLE]

where ${\mathcal{V}}_{t}$ denotes the set of non-zero element of ${\mathbf{V}}_{\text{R}}$.

V-A1 Design of $\mathbf{V}_{\mathrm{D}}$

We first group the sub-arrays at the TX according to the number of data streams $N_{s}$ satisfying $N_{s}\leq M_{t}$. We simply assign successive equal number of sub-arrays to each data stream if $M_{t}$ is multiples of $N_{s}$, otherwise, redundant sub-arrays are all assigned to any data stream. Note that this homogeneous grouping strategy is not a special case of the hybridly connected structure based partition strategy in [62]. We use the vector $\mathbf{d}_{t}\triangleq\left[d_{1},\cdots,d_{N_{s}}\right]$ to indicate the number of sub-arrays assigned to each data stream. The optimal grouping strategy for the DPA-MIMO system is scheduled for future research. Furthermore, assisted by the combination of sub-array grouping and variable decoupling, the proposed scheme can handle both the TX and the RX design with a flexible number of data steams in comparison with [9]. Thus, we have the precoders with the new structures as $\bar{\mathbf{V}}_{\mathrm{D}}\in\mathbb{C}^{N_{s}\times N_{s}}$ and $\bar{\mathbf{V}}_{\mathrm{R}}=\mathrm{blkdiag}\left\{\bar{\mathbf{v}}_{\mathrm{R}}^{\left[1\right]},\cdots,\bar{\mathbf{v}}_{\mathrm{R}}^{\left[N_{s}\right]}\right\}\in\mathbb{C}^{N_{t}^{\mathrm{tot}}\times N_{s}}$ where

[TABLE]

In Fig. 4, we give a example of sub-array grouping with $N_{s}=2$ and $M_{t}=4$. In this scenario, two data symbols are first processed by a $2\times 2$ digital precoding matrix. Then, each precoded data is allocated to two sub-arrays. Such a grouping operation facilitates the efficient design of RF precoding using SIC, as will be shown in the following.

With the defined variables $\bar{\mathbf{V}}_{\mathrm{R}}$ and $\bar{\mathbf{V}}_{\mathrm{D}}$ substituted into (40), we follow the common decoupling procedure [12] that given a fixed RF precoder $\bar{\mathbf{V}}_{\mathrm{R}}$ and an equivalent channel matrix $\mathbf{H}_{\mathrm{eq}}\triangleq\mathbf{H}\bar{\mathbf{V}}_{\mathrm{R}}$, the optimal digital precoder has a closed-form water-filling solution as

[TABLE]

where $\mathbf{D}_{t}=\bar{\mathbf{V}}_{\mathrm{R}}^{H}\bar{\mathbf{V}}_{\mathrm{R}}=\mathrm{diag}\left\{\mathbf{d}_{t}\right\}$, $\mathbf{U}_{\mathrm{e}}$ is the set of right singular vectors corresponding to the $N_{s}$ largest singular values of $\mathbf{H}_{\mathrm{eq}}\mathbf{D}_{t}^{-1/2}$ and $\bm{\Lambda}_{\mathrm{e}}$ is a diagonal matrix with the allocated powers to each data stream on its main diagonal.

V-A2 Design of $\mathbf{V}_{\mathrm{R}}$

With an equal power allocation scheme, i.e., $\bm{\Lambda}_{\mathrm{e}}\approx\mathbf{I}_{N_{s}}$, which shows a little loss in performance for moderate and high SNR regimes, we obtain the approximately optimal digital precoder $\bar{\mathbf{V}}_{\mathrm{D}}\approx\mathbf{D}_{t}^{-1/2}\mathbf{U}_{\mathrm{e}}$. Thus, the RF precoder can be obtained by solving the following problem

[TABLE]

where $\tilde{\mathbf{V}}_{\mathrm{R}}\triangleq\bar{\mathbf{V}}_{\mathrm{R}}\mathbf{D}_{t}^{-1/2}\in\mathbb{C}^{N_{t}^{\mathrm{tot}}\times N_{s}}$ and $\tilde{\mathcal{V}}_{t}$ denotes the set of non-zero element of $\tilde{\mathbf{V}}_{\text{R}}$.

According to the block structure of $\tilde{\mathbf{V}}_{\mathrm{R}}$, it is observed that problem (43) with nonconvex constraints can be decomposed into a series of simple subproblems, each of which only considers one specific group of sub-arrays [9, 62]. In particular, we can divide the matrix $\tilde{\mathbf{V}}_{\mathrm{R}}$ as $\tilde{\mathbf{V}}_{\mathrm{R}}=\left[\tilde{\mathbf{V}}_{N_{s}-1}\tilde{\mathbf{v}}_{N_{s}}\right]$, where $\tilde{\mathbf{v}}_{N_{s}}$ is the $N_{s}$th column and $\tilde{\mathbf{V}}_{N_{s}-1}$ is a matrix containing the first $N_{s}-1$ columns of $\tilde{\mathbf{V}}_{\mathrm{R}}$, respectively. Thus, the cost function of (43) can be written as

[TABLE]

where $\mathbf{R}_{n-1}=\mathbf{I}_{N_{r}^{\mathrm{tot}}}+P_{\mathrm{d}}\left(N_{s}\sigma_{\mathrm{z}}^{2}\right)^{-1}\mathbf{H}\tilde{\mathbf{V}}_{n-1}\tilde{\mathbf{V}}_{n-1}^{H}\mathbf{H}^{H}$ and $\mathbf{R}_{0}=\mathbf{I}_{N_{r}^{\mathrm{tot}}}$. This equality follows from $\left|\mathbf{I}+\mathbf{A}\mathbf{B}\right|=\left|\mathbf{I}+\mathbf{B}\mathbf{A}\right|$ and $\left|\mathbf{A}\mathbf{B}\right|=\left|\mathbf{A}\right|\left|\mathbf{B}\right|$ [57].

In (44), the total achievable rate can be a summation of the subrates of all data streams. Motivated by the idea of SIC for multi-user detection, we can first optimize the achievable subrate of the first data stream and then update $\mathbf{R}_{1}$. After some iterations, the optimal RF precoder design subproblem for the $n$th group of sub-arrays can be equivalently represented by

[TABLE]

where $\mathbf{T}_{n-1}=\mathbf{H}^{H}\mathbf{R}_{n-1}^{-1}\mathbf{H}$. Note that $\mathbf{T}_{n}$ can be iteratively obtained without matrix inverse, with its derivation shown in Appendix B. Due to the special structure of $\tilde{\mathbf{v}}_{n}$, (45) can be further simplified to

[TABLE]

where $\hat{\mathbf{T}}_{n-1}$ is a $d_{n}N_{t}^{\mathrm{sub}}\times d_{n}N_{t}^{\mathrm{sub}}$ Hermitian matrix formed as a submatrix of matrix $\mathbf{T}_{n-1}$ by taking the $\left(N_{t}^{\mathrm{sub}}\sum_{i=1}^{n-1}d_{i}+1\right)$th row and column to the $\left(N_{t}^{\mathrm{sub}}\sum_{i=1}^{n}d_{i}\right)$th row and column of $\mathbf{T}_{n-1}$. Since each element of $\bar{\mathbf{v}}_{\mathrm{R}}^{\left[n\right]}$ can be separated in (46), its optimal solution can be obtained by iteratively updating the following equation until convergence

[TABLE]

where the symbol $\measuredangle\left(.\right)$ extracts the corresponding phases of the element. Note that this iterative procedure is guaranteed to converge to a local optimum since the objective function of (46) increases in each iteration.

V-B Design of $\mathbf{U}_{\mathrm{D}}$ and $\mathbf{U}_{\mathrm{R}}$

V-B1 Design of $\mathbf{U}_{\mathrm{R}}$

We group the sub-arrays at the RX in the same way as the TX. Furthermore, we decouple the design of $\bar{\mathbf{U}}_{\mathrm{R}}$ and $\bar{\mathbf{U}}_{\mathrm{D}}$ by first optimizing the RF combiner with assumed ideal digital combiner and then finding the optimal digital combiner for the obtained RF combiner. Since $\bar{\mathbf{U}}_{\mathrm{R}}^{H}\bar{\mathbf{U}}_{\mathrm{R}}=\mathbf{D}_{r}=\mathrm{diag}\left\{\mathbf{d}_{r}\right\}$ where $\mathbf{d}_{r}\triangleq\left[\bar{d}_{1},\cdots,\bar{d}_{N_{s}}\right]$, we can define $\tilde{\mathbf{U}}_{\mathrm{R}}\triangleq\bar{\mathbf{U}}_{\mathrm{R}}{\mathbf{D}}_{r}^{-1/2}$, which results in the similar problem as (43):

[TABLE]

where $\mathcal{V}_{r}$ denotes the set of non-zero element of $\tilde{\mathbf{U}}_{\text{R}}$.

V-B2 Design of $\mathbf{U}_{\mathrm{D}}$

Assuming all other beamformers are fixed, the optimal digital combiner based on the minimum mean-square error (MMSE) criterion is formulated as [6]

[TABLE]

where $\mathbf{J}=\frac{P_{\mathrm{d}}}{N_{s}}\bar{\mathbf{U}}_{\mathrm{R}}^{H}\mathbf{H}\bar{\mathbf{V}}_{t}\bar{\mathbf{V}}_{t}^{H}\mathbf{H}^{H}\bar{\mathbf{U}}_{\mathrm{R}}+\sigma_{\mathrm{z}}^{2}\mathbf{D}_{r}\in\mathbb{C}^{N_{s}\times N_{s}}$ and $\bar{\mathbf{V}}_{t}=\bar{\mathbf{V}}_{\mathrm{R}}^{\star}\bar{\mathbf{V}}_{\mathrm{D}}^{\star}$. Finally, the proposed hybrid precoding approach is summarized in Algorithm 3.

VI Simulation Results

The performance of the proposed algorithms is evaluated through simulation with the following parameters. The TX and RX sub-arrays are ULAs of half-wavelength antenna spacing. For the sub-array spacing, we simply set $d_{a}=9\lambda_{c}$ for algorithm verification. Note that the practical value of $d_{a}$ has relationship with joint sparsity and it is worth pointing out that this relationship needs to be validated through extensive measurement and study of mmWave channels [5, 4]. Since the amplitude of the LoS components is typically 5 to $10\;\text{dB}$ stronger than that of the NLoS components at mmWave frequencies [63], the channel coefficients are generated through (1) with the variances of the channel paths as $\sigma_{\mathrm{LoS}}^{2}=1$ and $\sigma_{\mathrm{NLoS}}^{2}=10^{-0.5}$ [44]. We denote the channel common sparsity as $L_{c}=\left|\Omega_{c}\right|\geq 1$ and assume the equal channel individual sparsity among different sub-arrays as $L=\left|\Omega_{m,n}\right|+\left|\Omega_{c}\right|=5$ for simulation convenience. For the generation of the proposed DPA-MIMO channel, we first generate $L_{c}$ common paths one of which is the LoS path for all the sub-arrays. For any common path, the phase variation across the sub-arrays at the TX/RX induced by the sub-array spacing should be considered [11], i.e., for common path $i$ with $\vartheta_{m,n}^{\left(i\right)}={\vartheta^{\left(i\right)}}$ and $\psi_{m,n}^{\left(i\right)}={\psi^{\left(i\right)}}$, its complex amplitude is calculated by $\alpha_{m,n}^{\left(i\right)}=\alpha_{1,1}^{\left(i\right)}{e^{-j\left({\varphi_{r,m}^{\left(i\right)}+\varphi_{t,n}^{\left(i\right)}}\right)}}$ where $\varphi_{r,m}^{\left(i\right)}=\frac{{2\pi}}{{{\lambda_{c}}}}\left({m-1}\right)\left({{d_{a}}+\left({N_{r}^{{\text{sub}}}-1}\right){d_{e}}}\right)\cos\left({{\vartheta^{\left(i\right)}}}\right)$ and $\varphi_{t,n}^{\left(i\right)}=\frac{{2\pi}}{{{\lambda_{c}}}}\left({n-1}\right)\left({{d_{a}}+\left({N_{t}^{{\text{sub}}}-1}\right){d_{e}}}\right)\cos\left({{\psi^{\left(i\right)}}}\right)$. Then, we generate independent $\left(L-L_{c}\right)$ paths for each pair of the TX and RX sub-arrays.

In Algorithm 1, the threshold parameters are set to be $\delta_{1}=N\sigma_{\mathrm{z}}^{2}$ and $\delta_{2}=0.1N\sigma_{\mathrm{z}}^{2}$; the maximal iteration numbers are chosen as $T_{\mathrm{1,JOMP}}=L_{c}$ and $T_{\mathrm{2,JOMP}}=\left(L+2-T_{\mathrm{1,JOMP}}\right)K$. In Algorithm 2, we set the weight $\beta=3.3$, the iteration number $T_{\mathrm{JSBL}}=80$ and the error tolerance $\epsilon=10^{-4}$. In the following, two types of SNRs are considered: one is the pilot-to-noise ratio (PNR) defined as $10\log_{10}\left(P_{\mathrm{p}}/\sigma_{\mathrm{z}}^{2}\right)$, and the other is the data-to-noise ratio (DNR) defined as $10\log_{10}\left(P_{\mathrm{d}}/\sigma_{\mathrm{z}}^{2}\right)$. The performance metric for channel estimation is the normalized MSE (NMSE) defined as {\text{NMSE}}\triangleq 10{\log_{10}}\Big{(}{\mathbb{E}\Big{\{}{\frac{1}{{N_{t}^{{\text{tot}}}N_{r}^{{\text{tot}}}}}\Big{\|}{{\mathbf{H}}-\hat{\mathbf{H}}}\Big{\|}_{F}^{2}}\Big{\}}}\Big{)}. The hybrid precoding schemes based on the channel estimates $\hat{\mathbf{H}}$ are evaluated through the SE defined in (39).

This section consists of the following two parts:

• •

In the first part, we compare the NMSE of the proposed JOMP and JSBL-$\ell_{2}$ estimators with the conventional OMP and SBL estimators by employing the designed training beam patterns. Additionally, the non-cooperative DFT codebook based channel estimator, named as DFTB, is served as a benchmark. The comparison of pilot beam overhead and computational complexity for the channel estimators is summarized in Table I. On the one hand, the SBL and the proposed JSBL-$\ell_{2}$ algorithms require higher complexity than the DFTB, the OMP and the proposed JOMP algorithms. On the other hand, the CS-based algorithms embrace less pilot beam overhead than the DFTB algorithm when $N_{t}^{{\text{beam}}}N_{r}^{{\text{beam}}}<N_{t}^{{\text{tot}}}N_{r}^{{\text{tot}}}$.

• •

In the second part, we investigate the performance of the proposed hybrid precoding scheme termed as GSIC, the SDR based alternating optimization (SDR-AO) scheme [33], the hybrid precoding scheme [12] and the optimal fully-digital (FD) precoding scheme. We further compare the SE realized through the proposed hybrid beamformers based on the channel estimates obtained in the previous part.

VI-A Performance Evaluation of Channel Estimation

In Fig. 5, we compare the NMSE versus PNR for a DPA-MIMO setup with $M_{t}=M_{r}=4$, $N_{t}^{\mathrm{sub}}=N_{r}^{\mathrm{sub}}=12$, $N_{t}^{\mathrm{beam}}=24$, $N_{r}^{\mathrm{beam}}=36$ (partial-training case), $L_{c}=3$ and $L=5$. It is notable that the JSBL-$\ell_{2}$ estimator substantially outperforms the other estimators at low PNR. We further observe that the greedy JOMP estimator almost achieves the same channel estimation performance as the SBL estimator. At high PNR, the DFTB estimator with the full-training overhead achieves better performance over the other CS-based methods.

In Fig. 6, we show the NMSE versus the common sparsity $L_{c}$ varying from 1 to $L$ for a DPA-MIMO setup with $M_{t}=M_{r}=4$, $N_{t}^{\mathrm{sub}}=N_{r}^{\mathrm{sub}}=12$, $N_{t}^{\mathrm{beam}}=24$, $N_{r}^{\mathrm{beam}}=36$ (partial-training case) and ${\text{PNR}}=\;5\;{\text{dB}}$. For both the proposed JOMP and JSBL-$\ell_{2}$ estimators, better channel estimation performance is obtained with an increasing number of the common support $L_{c}$, while the OMP, SBL and DFTB estimators keep the constant NMSE. This is because the two customized estimators take advantage of the jointly sparse characteristic of the DPA-MIMO channel to enhance the quality of estimated channels. Moreover, for $L_{c}\geq 3$, the JSBL-$\ell_{2}$ estimator surpasses the DFTB estimator in the accuracy of channel estimates.

In Fig. 7, we further investigate that with how much training the proposed two CS-based methods can be competitive to the full-training based DFTB estimator for different PNRs. The number of training beams at both ends is assumed to the same as $N^{\mathrm{beam}}$ for a DPA-MIMO setup with $M_{t}=M_{r}=4$, $N_{t}^{\mathrm{sub}}=N_{r}^{\mathrm{sub}}=12$, $L_{c}=3$ and $L=5$. In Fig. LABEL:sub@fig_nmse_vs_beam_5dB, with an increasing $N^{\mathrm{beam}}$, the NMSE of the proposed two estimators decreases monotonically at ${\text{PNR}}=\;5\;{\text{dB}}$. More specifically, the JOMP and JSBL-$\ell_{2}$ estimators can realize the nearly same reconstruction accuracy of the DFTB estimator with its $\frac{3}{4}\times\frac{3}{4}=56.25\%$ and $\frac{1}{2}\times\frac{1}{2}=25\%$ training overhead, respectively. As shown in Fig. LABEL:sub@fig_nmse_vs_beam_20dB at ${\text{PNR}}=\;20\;{\text{dB}}$, the JOMP estimator finally approaches the same NMSE value as the DFTB estimator with $N^{\mathrm{beam}}=48$, while the JSBL-$\ell_{2}$ estimator still maintains the performance advantage with a high PNR.

VI-B Performance Evaluation of Hybrid Precoding

First, we compare the SE versus DNR for a various number of sub-arrays at the TX/RX with perfect CSI and a fixed number of total antennas. We set the parameters as $N_{t}^{\mathrm{tot}}=N_{r}^{\mathrm{tot}}=90$, $L_{c}=3$ and $L=5$. It is observed from Fig. LABEL:sub@fig_hybrid_precoding_perfect_csi_M_3_Ns_3 and Fig. LABEL:sub@fig_hybrid_precoding_perfect_csi_M_6_Ns_3 that there are distinct gaps between the FD scheme and the other three hybrid precoding schemes, which results from the beamforming gain loss of the sub-array based structure compared with the FD structure [33]. We further find that the proposed GSIC scheme can successfully surpass the SDR-AO and SVD-based schemes. The efficiency of the proposed GSIC scheme is verified through average central processing unit processing time. By using the same configuration as the Matlab experiment demonstrated in Fig. LABEL:sub@fig_hybrid_precoding_perfect_csi_M_6_Ns_3, the proposed GSIC scheme almost run 115 times faster than the SDR-AO scheme [33] while only about 1.5 times slower than the hybrid precoding scheme in [12]. Therefore, in the following, we only employ the proposed GSIC scheme to design the hybrid beamformers based on the estimated channels obtained in Subsection VI-A.

In Fig. 9, we investigate the SE versus DNR for a DPA-MIMO setup with $M_{t}=M_{r}=4$, $N_{t}^{\text{sub}}=N_{r}^{\text{sub}}=12$, $N_{s}=4$, $L=5$, $L_{c}=3$, $\text{PNR}=5\;\text{dB}$, $N_{t}^{\text{beam}}=24$ and $N_{r}^{\text{beam}}=36$. Obviously, the improved channel estimation accuracy of the proposed algorithms are reflected in their enhanced SEs. Furthermore, the SE curve of the proposed JSBL-$\ell_{2}$ estimator nearly coincides with that of the DFTB estimator. In addition, the proposed less-complexity JOMP estimator can provide the same SE as the SBL estimator.

In Fig. 10, we present the SE versus the common sparsity $L_{c}$ varing from 1 to $L$ for a DPA-MIMO setup with $M_{t}=M_{r}=4$, $N_{t}^{\text{sub}}=N_{r}^{\text{sub}}=12$, $N_{s}=4$, $L=5$, $\text{PNR}=5\;\text{dB}$, $N_{t}^{\text{beam}}=24$ and $N_{r}^{\text{beam}}=36$. It is interesting that the SE of any precoder decreases with an increasing number of common supports. This mainly comes from that the spatial degrees of freedom (DoFs) are reduced due to less independent scatterers of the DPA-MIMO channels. Additionally, the SE curves of the proposed JOMP and JSBL-$\ell_{2}$ estimators approach that of the DFTB estimator as the number of the common supports increases.

In Fig. 11, we depict the SE versus the number of training beams for a DPA-MIMO setup with $M_{t}=M_{r}=4$, $N_{t}^{\text{sub}}=N_{r}^{\text{sub}}=12$, $N_{s}=4$, $L=5$, $L_{c}=3$, $\text{DNR}=5\;\text{dB}$ and $\text{PNR}=5\;\text{dB}$. The SE of the proposed JOMP and JSBL-$\ell_{2}$ estimators evidently boost towards the perfect CSI case with an increasing number of training beams. Furthermore, the SE of the proposed two estimators approximates that of the DFTB estimator (full-training case) with only its $\frac{3}{4}\times\frac{3}{4}=56.25\%$ training overhead.

VII Conclusion

In this paper, we have focused on modeling and analysis of the narrowband DPA-MIMO based transceiver system, and designed efficient channel estimation and hybrid precoding schemes for such distributed array-of-sub-arrays architecture working on mmWave bands. Based on the reasonable analysis in Section II, the DPA-MIMO channel has high probability to manifest a hidden structured sparsity in the beam-domain channel vector due to the partially shared scatterers among the distributed sub-arrays at mmWave frequencies at the TX/RX. In light of this characteristic, we have formulated a structured SMV problem that estimates the AoDs, AoAs and the corresponding gain of significant paths. In order to guarantee the good recovery performance and decrease the training feedback overhead, the open-loop training beam patterns were designed through minimizing the total coherence of the equivalent measurement matrix. The simulation and comparison results have demonstrated that the proposed channel estimators can better exploit the structured channel properties defined in Definition 1 than the existing CS-based estimators such as the OMP and SBL estimators, and the proposed hybrid precoding method enjoys the low-complexity while achieving good performance.

More realistic channel modeling for DPA-MIMO is expected for future research. Specifically, the channel parameters including common and local scattering components can be investigated by using ray-tracing tools, and the scatter evolution on both sub-array and time axes should be also taken into account [64]. In addition, interesting and practical topics in DPA-MIMO applications cover many diversified situations, such as multi-user channel acquisition and hybrid precoding schemes [8], and optimal sub-array placement for preventing blockage [42].

Appendix A Derivation of (26)

Let $\mathbf{s}=\mathbf{x}-\mathbf{c}$, $\mathbf{A}=\left(\bm{\Gamma}^{s}\right)^{-1}$ and $\mathbf{B}=\left(\mathbf{I}_{K}\otimes\bm{\Gamma}^{c}\right)^{-1}$, thus we can transform the objective function in (26) to the following

[TABLE]

For the fixed $\mathbf{x}$, we have an unconstrained quadratic function only with respect to $\mathbf{c}$ and get its optimal solution as $\mathbf{c}^{\star}=\left(\mathbf{A}+\mathbf{B}\right)^{-1}\mathbf{A}\mathbf{x}$. After submitting the optimal $\mathbf{c}^{\star}$ into (50), we have

[TABLE]

where (a) follows from the Woodbury identity [57], and (b) follows from the identity [39]

[TABLE]

with $\bm{\Gamma}=\mathbf{A}^{-1}+\mathbf{B}^{-1}$.

Obviously, the optimal value of $\mathbf{x}$ is expressed as

[TABLE]

Given the following matrix identity [57]

[TABLE]

we further obtain

[TABLE]

Appendix B Iterative computation of $\mathbf{T}_{n}$ in (46)

In order to compute $\mathbf{T}_{n}$ iteratively, we further express $\mathbf{R}_{n}$ as $\mathbf{R}_{n}=\mathbf{R}_{n-1}+\frac{P_{\mathrm{d}}}{N_{s}\sigma_{\mathrm{z}}^{2}}\mathbf{H}\tilde{\mathbf{v}}_{n}\tilde{\mathbf{v}}_{n}^{H}\mathbf{H}^{H}.$ With the Sherman-Morrison formula [57]

[TABLE]

we have

[TABLE]

where $\rho_{n}=\frac{P_{\mathrm{d}}}{N_{s}\sigma_{\mathrm{z}}^{2}}\tilde{\mathbf{v}}_{n}^{H}\mathbf{T}_{n-1}\tilde{\mathbf{v}}_{n}$. Thus, the matrix $\mathbf{T}_{n}$ can be iteratively obtained by

[TABLE]

Acknowledgment

The authors would like to thank the Associate Editor and the anonymous reviewers for their valuable comments and helpful suggestions.