Message passing-based link configuration in short range millimeter wave systems

TL;DR

This paper introduces a message passing framework for short-range millimeter wave systems that exploits channel structure to improve beam alignment with fewer measurements, addressing limitations of far-field assumptions.

Contribution

A novel message passing-based method that models short-range mmWave channels by dividing them into subchannels, leveraging structure and antenna geometry for efficient beam alignment.

Findings

Achieves better beam alignment with fewer measurements.

Outperforms standard compressed sensing techniques.

Effectively exploits channel structure across subchannels.

Abstract

Millimeter wave (mmWave) communication in typical wearable and data center settings is short range. As the distance between the transmitter and the receiver in short range scenarios can be comparable to the length of the antenna arrays, the common far field approximation for the channel may not be applicable. As a result, dictionaries that result in a sparse channel representation in the far field setting may not be appropriate for short distances. In this paper, we develop a novel framework to exploit the structure in short range mmWave channels. The proposed method splits the channel into several subchannels for which the far field approximation can be applied. Then, the structure within and across different subchannels is leveraged using message passing. We show how information about the antenna array geometry can be used to design message passing factors that incorporate structure across successive subchannels. Simulation results indicate that our framework can be used to achieve better beam alignment with fewer channel measurements when compared to standard compressed sensing-based techniques that do not exploit structure across subchannels.

Equations

$$y_k[m]={\mathbf{w}}_k^T[m]{\mathbf{H}}_k\mathbf{f}[m]+v_k[m],\forall k\in \{1,2,...N_{\mathrm{RF}}\},$$

$$\mathbf{H}(i,j)=\frac{\lambda }{4\pi d_{ij}}{e}^{-\mathsf{j}2\pi d_{ij}/\lambda }.$$

$${\mathbf{a}}_{{}_N}(\Delta )=[1,e^{-\mathsf{j}\pi \sin \Delta },e^{-\mathsf{j}2\pi \sin \Delta },..,e^{-\mathsf{j}(N-1)\pi \sin \Delta }{]}^T.$$

$${\mathbf{H}}_k={\mathbf{U}}_N{\mathbf{X}}_k{\mathbf{U}}_N.$$

$$z\left[n\right]$$

$$y_k[m]=({\mathbf{J}}_{r_k[m]}\mathbf{z}{)}^T{\mathbf{U}}_N{\mathbf{X}}_k{\mathbf{U}}_N{\mathbf{J}}_{c[m]}\mathbf{z}+v_k[m].$$

$$y_k[m]={\mathbf{e}}_{r_k[m]}^T{\mathbf{U}}_N{\mathbf{S}}_k{\mathbf{U}}_N{\mathbf{e}}_{c[m]}+v_k[m].$$

$${\mathbf{A}}_k(m,:)=({\mathbf{e}}_{c[m]}^T{\mathbf{U}}_N)\otimes ({\mathbf{e}}_{r_k[m]}^T{\mathbf{U}}_N).$$

$$y_k[m]={\mathbf{A}}_k{\mathbf{s}}_k+v_k[m].$$

$$p(y_k[m]|{\mathsf{s}}_k)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{|y_k[m]-{\mathbf{A}}_k(m,:){\mathsf{s}}_k{|}^2}{2{\sigma }^2}\right).$$

$${\rho }_{k,\mathrm{D}}=\frac{\mathrm{SNR}}{[({\mathbf{H}}_{\mathrm{UL},\mathrm{D}}^{\ast }{\mathbf{H}}_{\mathrm{UL},\mathrm{D}}+\mathbf{I}/\mathrm{SNR}{)}^{-1}{]}_{k,k}}-1.$$

$$y_{\mathrm{AP},k}[m]={\alpha }_{\mathrm{AP},k}{\mathbf{w}}_k^T[m]{\mathbf{a}}_N({\omega }_k)+v_k[m].$$

$$p({\omega }_k)=\exp \left(-|{\hat{\alpha }}_{\mathrm{AP},k}{|}^2\Vert {\tilde{\mathbf{y}}}_{\mathrm{AP},k}-{\mathbf{\Psi }}_k{\mathbf{a}}_N({\omega }_k){\Vert }_2^2/{\sigma }^2\right).$$

$$(d_1\cos {\omega }_1{)}^2+(d_1\sin {\omega }_1+{\ell }_1{)}^2=d^2$$

$${\omega }_2={\tan }^{-1}\left(\frac{{\ell }_1-{\ell }_2+d_1\sin {\omega }_1}{d_1\cos {\omega }_1}\right).$$

$$p_{\mathrm{g}}({\omega }_2|{\omega }_1)=\frac{1}{D}{\int }_{d_{\min }}^{d_{\max }}\delta ({\omega }_2-\mathcal{G}({\omega }_1,r))dr,$$

$$P_{\mathrm{in}}^{\mathrm{fwd}}({\omega }_2)={\int }_{-\pi /2}^{\pi /2}{P}_{\mathrm{out}}^{\mathrm{fwd}}({\omega }_1)p_{\mathrm{g}}({\omega }_2|{\omega }_1)d{\omega }_1.$$

$$P_{\mathrm{out}}^{\mathrm{fwd}}({\omega }_2)=p({\omega }_2)P_{\mathrm{in}}^{\mathrm{fwd}}({\omega }_2).$$

$$P_{\mathrm{in}}^{\mathrm{bwd}}({\omega }_{k-1})$$

$$P_{\mathrm{out}}^{\mathrm{bwd}}({\omega }_k)$$

$$p_{{}_{\mathrm{GMP}}}({\omega }_k)=p({\omega }_k)P_{\mathrm{in}}^{\mathrm{fwd}}({\omega }_k)P_{\mathrm{in}}^{\mathrm{bwd}}({\omega }_k).$$

Full text

I Introduction

Millimeter wave (mmWave) systems can support high data rates and enable emerging applications like augmented reality and virtual reality [2]. Integrating mmWave radios in small scale devices like wearables can be easier as large antenna arrays occupy a small footprint at mmWave when compared to lower frequency systems [3]. Challenges like cost and power consumption in mmWave radios can be addressed by using analog beamforming-based radio designs or low-resolution receivers [4]. To support high data rates in mmWave systems, the wireless link between the transmitter and the receiver must be configured properly.

Establishing the link between the transmitting and receiving radios can be challenging at mmWave due to hardware constraints [4]. As typical mmWave beamforming systems have a large number of antennas and fewer radio frequency (RF) chains, conventional beam alignment techniques may require a lot of training overhead when applied to such systems [5]. Algorithms that exploit structure in mmWave channels, like sparsity in an appropriate dictionary or low rank, can be useful in mmWave settings [4]. Structure-aware algorithms, which are based on compressed sensing (CS) or matrix completion, can learn the mmWave channel or perform beam alignment with fewer channel measurements when compared to standard brute-force search approaches [6, 7]. As a result, structure-aware algorithms can enable rapid mmWave link configuration. Most sparsity- or low rank-aware algorithms, however, make the far field assumption to model the multiple-input multiple-output (MIMO) channel [6, 7]. In typical mmWave wearable settings where the distance between the transmitting and the receiving radios can be comparable to the length of the antenna arrays, the far field approximation may cease to hold [8]. Therefore, it is necessary to develop new channel estimation or beam alignment techniques that are tailored to short range mmWave systems.

In this paper, we propose a new framework to leverage structure in short range mmWave communication channels. Our framework assumes a subarray-based hybrid beamforming system [4], and is based on two key observations. First, the far field approximation can be applied for a subchannel although it may fail for a full channel. As a result, it is reasonable to assume that the subchannels have a sparse representation at mmWave, in an appropriate dictionary. Second, the channel entries across appropriately chosen subchannels are correlated if the antenna arrays corresponding to the subchannels are co-located. Our framework splits the MIMO channel into several subchannels to exploit structure within and across the subchannels.

The correlation across the sparse subchannel representations can be used to regularize the channel estimation or beam alignment problems. Prior work on dynamic CS has shown how signal recovery techniques that account for both correlation and sparsity of a group of vectors can outperform standard CS that solves for the vectors independently [9, 10]. Our framework allows applying such dynamic CS techniques for channel estimation and beam alignment in short range settings. Dynamic CS-based channel estimation techniques were developed for far field systems in [6, 11, 12, 13] to exploit structure along time, frequency or spatial dimensions. For example, the far field channel was split into a group of subchannels for compressive subspace estimation [14]. The splitting in [14] was performed along the frequency dimension to deal with beam squint in wideband systems. To the best of our knowledge, however, short range link configuration techniques that are based on dynamic CS have not been developed. Furthermore, the proposed geometry-aware message passing technique, which incorporates information about the antenna geometry, is novel. We summarize the main contributions of our work as follows.

• •

We model a short range mmWave MIMO channel as a collection of several subchannels that satisfy the far field approximation. We use an appropriate sparsifying dictionary for the subchannels, and show that the support and amplitude of the sparse coefficients exhibit smooth dynamics across subchannels.

• •

We show how dynamic compressed sensing-approximate message passing (DCS-AMP) [10], a dynamic CS algorithm, can be used to leverage the sparsity and correlation in subchannels for channel estimation. We also provide a low complexity implementation of DCS-AMP by establishing an equivalence between subchannel estimation and dynamic magnetic resonance imaging (MRI) CS.

• •

We develop an alternative message passing technique that uses information about the antenna geometry and the range of the transceiver distance, for compressive beam alignment. The proposed algorithm solves for the angles-of-arrival in subarrays, and can exploit subchannel structure beyond sparsity and smooth variation in subchannel coefficients.

• •

We evaluate the performance of DCS-AMP and the proposed geometry-aided message passing algorithm, for a short range mmWave scenario in a multi-user setting. Our results indicate that the DCS-AMP-based approach is useful in the low SNR regime, while the proposed algorithm achieves better beam alignment at high SNR.

The link configuration techniques based on DCS-AMP and geometry-aided message passing use a different pilot structure, and exploit subchannel dynamics in different ways. On the one hand, pilot transmissions in the DCS-AMP-based approach are only performed in the uplink, to estimate the MIMO channel. On the other hand, the proposed geometry-aided message passing-based method uses both uplink and downlink pilot transmissions, as in the IEEE 802.11ad standard [5]. Furthermore, the geometry-based method does not intend to estimate the MIMO channel, but solves for the local angles-of-arrivals (AoAs) that are defined for each subarray. As a result, the proposed technique requires low complexity and limited feedback when compared to the DCS-AMP-based method. Our geometry-aided message passing algorithm, however, can only recover the line-of-sight (LoS) component, unlike the DCS-AMP-based approach that can also recover the non-LoS components in the channel. The DCS-AMP-based technique accounts for the smooth variation in subchannels by assuming that the subchannel coefficients follow a Gauss-Markov dynamics. The proposed geometry-aided message passing algorithm explicitly captures the variation in subchannels by using message passing factors that are computed using the array geometry.

Notation : $\mathbf{A}$ is a matrix, $\mathbf{a}$ is a column vector and $a,A$ denote scalars. Using this notation $\mathbf{A}^{T},\overline{\mathbf{A}}$ and $\mathbf{A}^{\ast}$ represent the transpose, conjugate and conjugate transpose of $\mathbf{A}$. The matrix $[\mathbf{A};\mathbf{B}]$ is obtained by vertically stacking $\mathbf{A}$ and $\mathbf{B}$. The spectral norm and the Frobenius norm of $\mathbf{A}$ are denoted by $\|\mathbf{A}\|_{2}$ and $\|\mathbf{A}\|_{F}$. The scalar $a\left[m\right]$ denotes the $m^{\mathrm{th}}$ element of $\mathbf{a}$, and $\mathbf{A}\left(k,\ell\right)$ is the entry of $\mathbf{A}$ in the $k^{\mathrm{th}}$ row and ${\ell}^{\mathrm{th}}$ column. The $k^{\mathrm{th}}$ row of $\mathbf{A}$ is denoted as $\mathbf{A}(k,:)$. The matrix $|\mathbf{A}|$ contains the element-wise magnitude of $\mathbf{A}$. The symbol $\odot$ is used for the Hadamard product. The matrices $\mathbf{U}_{N}$ and $\mathbf{I}_{N}$ in $\mathbb{C}^{N\times N}$ denote the unitary Discrete Fourier Transform (DFT) matrix and the identity matrix. We use $\mathbf{e}_{k}$ to represent the $(k+1)^{\mathrm{th}}$ canonical basis vector. The set $\mathcal{I}_{N}=\{1,2,3,.\,.\,.,N\}$. $\mathcal{N}_{c}(\mu,v_{a})$ denotes the complex Gaussian distribution with mean $\mu$ and variance $v_{a}$. We define $\mathsf{j}=\sqrt{-1}$.

II System and channel model

In this section, we explain the system and channel model in a short range mmWave setting. For simplicity of exposition, we assume a point-to-point link in a narrowband setting. The simulation results in Section V consider a more complex multi-user scenario.

We consider a mmWave setting in Fig. 1 where the access point (AP) has a subarray-based hybrid beamforming architecture with $N_{\mathrm{RF}}$ RF chains. Each of the $N_{\mathrm{RF}}$ subarrays at the AP is a half-wavelength spaced uniform linear array (ULA) with $N$ antennas. The $N$ antennas of a subarray are connected to a single RF chain through phase shifters. Therefore, each subarray consists of a single phased array. We use $\lambda$ to denote the wavelength corresponding to the mmWave carrier frequency. The resolution of all the $NN_{\mathrm{RF}}$ phase shifters at the AP is assumed to be $q$ bits. The use of low resolution phase shifters results in a lower power consumption, and introduces new challenges for beam alignment. These challenges can be addressed by a careful design of codebooks [15]. The alphabet for phase shifts is defined as $\mathbb{Q}=\{e^{\mathsf{j}2\pi\ell/2^{q}}/\sqrt{N}:\ell\in\mathcal{I}_{2^{q}}\}$. The spacing between consecutive subarrays is assumed to be constant; this spacing can be larger than $\lambda/2$. The antennas at the AP are placed along the horizontal axis and occupy a length of $L_{\mathrm{AP}}$ units. Consider a single station (STA) in Fig. 1 equipped with a half-wavelength spaced ULA of $N$ elements. We assume that the antenna arrays at the AP and the STA operate at the same polarization and are coplanar. Let $d$ denote the distance between the midpoints of the antenna arrays at the AP and the STA, and $L_{\mathrm{STA}}$ denote the length of the antenna array at the STA. We define the angle $\gamma$ as the orientation of the STA relative to the normal of the AP array, and $\theta$ as the angle made by the STA array with the horizontal. The ULA at the STA is equipped with phase shifters and a single RF chain. Our framework can also be extended to other antenna architectures like uniform planar arrays, by using appropriate array response vectors in the formulation.

Now, we describe the system model used to obtain uplink channel measurements. We assume a narrowband channel for simplicity. Extending our framework to incorporate frequency selectiveness of channels is an interesting research direction. We consider pilot transmissions in the uplink as the AP may have better computational capabilities than the STA, and can acquire $N_{\mathrm{RF}}$ channel measurements in parallel. Let $\mathbf{H}_{k}\in\mathbb{C}^{N\times N}$ be the baseband equivalent of the short range MIMO channel matrix between the STA and the $k^{\mathrm{th}}$ subarray at the AP. The matrix $\mathbf{H}_{k}$ is defined as the $k^{\mathrm{th}}$ subchannel in the MIMO system. We define the full MIMO channel matrix as $\mathbf{H}\in\mathbb{C}^{NN_{\mathrm{RF}}\times N}$. The channel $\mathbf{H}$ is defined as $\mathbf{H}=[\mathbf{H}_{1};\mathbf{H}_{2};\cdots;\mathbf{H}_{N_{\mathrm{RF}}}]$. Channel measurements in the MIMO system are acquired by applying beam training vectors to the phased arrays at the STA and the AP. Let ${\mathbf{f}}[m]\in\mathbb{Q}^{N}$ be the transmit beam training vector applied to the phased array at the STA in the $m^{\mathrm{th}}$ training slot. During this slot, the AP acquires channel measurements by applying ${\mathbf{w}_{k}}[m]\in\mathbb{Q}^{N}$ as the receive beam training vector to its $k^{\mathrm{th}}$ subarray. The channel measurement acquired by the $k^{\mathrm{th}}$ subarray in the $m^{\mathrm{th}}$ slot is defined as

[TABLE]

where $v_{k}[m]\sim\mathcal{N}_{c}(0,\sigma^{2})$ is additive white Gaussian noise. It can be noticed from (1) that $N_{\mathrm{RF}}$ projections of the channel matrix can be obtained in a given training slot. As $\mathbf{H}$ consists of $N_{\mathrm{RF}}N^{2}$ entries, beam alignment through conventional MIMO channel estimation requires $N^{2}$ training slots. In this paper, we show that a reasonable approximation of the MIMO channel $\mathbf{H}$ can be estimated in far fewer than $N^{2}$ slots by exploiting structure in the short range mmWave communication channel. The channel approximation obtained is then used to configure the antenna arrays at the AP and the STA.

To explain the key ideas underlying our framework, we consider a geometric LoS channel model in a short range setting. The simulation results in Section V, however, consider a more realistic short range scenario that also includes non-LoS components. We use $d_{ij}$ to denote the distance between the $i^{\mathrm{th}}$ antenna element at the AP and the $j^{\mathrm{th}}$ element at the STA. Here, $i\in\mathcal{I}_{NN_{\mathrm{RF}}}$ and $j\in\mathcal{I}_{N}$. The discrete time baseband equivalent of the MIMO channel matrix can be expressed as [8]

[TABLE]

The channel matrix with entries defined by (2) can be approximated as a rank-one matrix when $d$ is significantly larger than the length of the antenna arrays at the AP and the STA, i.e., $d\gg\mathrm{max}(L,(N-1)\lambda/2)$ [16]. In settings where $d$ can be comparable to the length of the antenna array, the full channel matrix, i.e., $\mathbf{H}$, cannot always be approximated as a rank-one matrix.

Now, we show how subchannels exhibit a rank-one characteristic when compared to the full channel for reasonably short distances. To explain our argument, we define an energy metric that measures the rank-one nature of a channel. We define $E_{\mathrm{F}}(d)=\|\mathbf{H}\|^{2}_{2}/\|\mathbf{H}\|^{2}_{F}$ and $E_{\mathrm{S}}(d)=\sum^{N_{\mathrm{RF}}}_{k=1}\|\mathbf{H}_{k}\|^{2}_{2}/(\|\mathbf{H}_{k}\|^{2}_{F}N_{\mathrm{RF}})$ as the energy metrics for the full channel and the subchannels. While $E_{\mathrm{F}}(d)$ represents the fraction of energy contained in the rank-one approximation of the full channel, $E_{\mathrm{S}}(d)$ denotes the average of the fractions computed for each subchannel.

For the short range communication setting in Fig. 1, it can be observed from Fig. 2 that the energy metric increases as a function of distance $d$. The observation is consistent with the fact that LoS MIMO channels can be well approximated as rank-one in the far field setting. An important observation from Fig. 2 is that the energy metric for subchannels, i.e., $E_{\mathrm{S}}(d)$, is higher than that for the full channel. Therefore, subchannels can be better approximated as rank-one matrices when compared to the full channel. The approximation is reasonable as the far field assumption is likely to be valid for arrays of shorter length, i.e., subarrays, when compared to the full array.

The array response vectors in a subchannel can be modeled as Vandermonde vectors when the rank-one approximation is valid for the subchannel [16]. For the half wavelength-spaced ULA, we define the array response vector of dimension $N$ as

[TABLE]

We use $\alpha_{k}$ to denote the path gain of the ray between the STA and the $k^{\mathrm{th}}$ subarray at the AP. The angle of departure and the angle of arrival of the ray, relative to the boresight, are denoted by $\phi_{k}$ and $\omega_{k}$. Note that a unique angle of arrival and departure can be defined for an array only when the far field approximation, equivalently the rank-one approximation, is valid. The angles $\{\omega_{k}\}_{k=1}^{N_{\mathrm{RF}}}$ are called the local AoAs, as they are defined for subarrays or sections of the AP array. Under the rank-one approximation of each subchannel, we can write $\mathbf{H}_{k}\approx\alpha_{k}\mathbf{a}_{{}_{N}}(\omega_{k})\mathbf{a}^{T}_{{}_{N}}(\phi_{k})$ for every $k\in\mathcal{I}_{N_{\mathrm{RF}}}$. Due to the Vandermonde structure of the array response vectors, $\mathbf{H}_{k}$ is compressible in the 2D-DFT dictionary. We use $\mathbf{X}_{k}$ to denote the inverse 2D-DFT of $\mathbf{H}_{k}$ such that

[TABLE]

The matrix $\mathbf{X}_{k}$, defined as the beamspace subchannel, is approximately sparse [4]. It is not exactly sparse as the angles of departure and arrival may not be aligned with those corresponding to the DFT dictionary [6].

Now, we argue that the non-zero coefficients in $\{\mathbf{X}_{k}\}_{k=1}^{N_{\mathrm{RF}}}$ can exhibit smooth dynamics across the subchannel index $k$ in short range settings. The angles of arrival and departure determine the support of the non-zero coefficients in the beamspace subchannel [17]. As the subarrays at the AP are co-located, the angles of departure and arrival, i.e., $\{\phi_{k}\}_{k=1}^{N_{\mathrm{RF}}}$ and $\{\omega_{k}\}_{k=1}^{N_{\mathrm{RF}}}$, vary slowly with $k$. Therefore, the support of the non-zero coefficients in $\mathbf{X}_{k}$ drifts slowly with $k$. For the extreme case when $d\gg\mathrm{max}(L,(N-1)\lambda/2)$, it can be shown that the drift is negligible. An illustration of the dynamics of the coefficients in $\{\mathbf{X}_{k}\}_{k=1}^{N_{\mathrm{RF}}}$ is provided in Fig. 3 for the short range communication setting in Fig. 1. The beamspace subchannels in the short range system can be estimated from sub-Nyquist channel measurements by exploiting the sparse and correlated nature of subchannels.

III Short range channel estimation using dynamic compressed sensing

In this section, we transform the beamspace subchannel estimation problem into an equivalent masked beamspace subchannel estimation problem by using Zadoff-Chu sequences in the antenna domain [18]. Under such a transformation, we show that short range subchannel estimation is exactly equivalent to magnetic resonance imaging of a sequence of correlated angiogram images [10]. Then, we show how DCS-AMP algorithm [10] can be used to estimate the masked beamspace subchannels with $M\ll N^{2}$ pilot transmissions.

III-A Zadoff-Chu-based training and the spectral mask concept

Zadoff-Chu (ZC) sequences have been extensively applied along the time and frequency dimensions in wireless systems due to their constant amplitude and zero autocorrelation (CAZAC) properties [19]. Recent work has demonstrated the advantages of using ZC sequences over the common random phase shift sequences for CS-based beam alignment [18, 20, 21]. In this section, we use ZC sequences along the spatial dimension for compressive channel estimation in the subarray-based hybrid beamforming architecture shown in Fig. 1. We define $\mathbf{z}\in\mathbb{C}^{N}$ as a unit norm ZC sequence of root $t\in\mathcal{I}_{N}$. The root $t$ must be co-prime with $N$. The $n^{\mathrm{th}}$ element of $\mathbf{z}$ is defined as [19]

[TABLE]

The constant amplitude of the entries of the ZC sequence allows realizing $\mathbf{z}$ in phased arrays. The constant amplitude of the DFT of $\mathbf{z}$ allows efficient compressed sensing-based beam alignment of 2D-DFT sparse channels [20]. For an extensive treatment on the optimality of ZC-based channel acquisition for CS-based beam alignment, we refer the interested reader to [20] and [21].

The channel measurements for DCS-AMP-based channel estimation are acquired by applying random circulant shifts of the ZC sequence to the phased arrays at the AP and the STA. For the $m^{\mathrm{th}}$ transmission in a sequence of $M$ uplink pilot transmissions, the STA applies a $c[m]\in\mathcal{I}_{N}$ circulantly shifted version of $\mathbf{z}$ to its phased array. At the same time, the AP uses an $r_{k}[m]\in\mathcal{I}_{N}$ circulantly shifted version of $\mathbf{z}$ as the beam training vector for its $k^{\mathrm{th}}$ subarray. The pair of ZC sequences used at the STA and the $k^{\mathrm{th}}$ subarray at the AP determine the channel measurement corresponding to $\mathbf{H}_{k}$. As the STA and an AP subarray can independently apply $N$ different circulant shifts to their phased array, a maximum of $N^{2}$ spatially distinct subchannel measurements can be obtained with the ZC-based training. In this paper, we focus on the sub-Nyquist setting, i.e., $M<N^{2}$ measurements are acquired per subchannel. For each $m$, the STA chooses $c[m]\in\mathcal{I}_{N}$ at random with replacement. The AP chooses $\{r_{k}[m]\}^{N_{\mathrm{RF}}}_{k=1}$ at random such that none of the coordinates in $\{(r_{k}[m],c[m])\}^{N_{\mathrm{RF}}}_{k=1}$ were sampled in the previous random draws. Each of the $N_{\mathrm{RF}}$ subchannels is projected onto a set of $M$ linearly independent matrices that are determined by the random sampling scheme.

Now, we define masked beamspace subchannels using the spectral mask concept in [18]. We define $\mathbf{J}\in\mathbb{R}^{N\times N}$ as a circulant matrix with its first row as $\left(0,1,0,...,0\right)$. The subsequent rows of $\mathbf{J}$ are generated by circulantly shifting the previous row to the right by one unit. We define the $\ell$-circulant delay matrix as $\mathbf{J}_{\ell}=\mathbf{J}\mathbf{J}\cdot\cdot\cdot\mathbf{J}\,(\mathrm{\ell\,times})$. The matrix $\mathbf{J}_{0}$ is defined as the $N\times N$ identity matrix. Using these definitions, it can be observed that the beam training vector used by the STA for the $m^{\mathrm{th}}$ pilot transmission is $\mathbf{f}[m]=\mathbf{J}_{c[m]}\mathbf{z}$. Similarly, the AP applies $\mathbf{w}_{k}[m]=\mathbf{J}_{r_{k}[m]}\mathbf{z}$ to its $k^{\mathrm{th}}$ subarray, to acquire a channel measurement for the $m^{\mathrm{th}}$ pilot transmission. From the system model in (1), the channel measurement $y_{k}[m]$ can be expressed as

[TABLE]

We define $\boldsymbol{\Lambda}_{\mathbf{z}}=\mathrm{diag}(\sqrt{N}\mathbf{U}_{N}\mathbf{z})$ as a diagonal matrix that contains the scaled DFT of $\mathbf{z}$ along its diagonal. The matrix $\boldsymbol{\Lambda}_{\mathbf{z}}$ is unimodular along its diagonal by the zero autocorrelation property of the ZC sequence [20]. For each $k$, the masked beamspace subchannel associated with $\mathbf{X}_{k}$ is defined as $\mathbf{S}_{k}=\boldsymbol{\Lambda}_{\mathbf{z}}\mathbf{X}_{k}\boldsymbol{\Lambda}_{\mathbf{z}}$. Using the spectral mask concept [18], (6) can be rewritten as

[TABLE]

The spectral mask concept [18] shows that estimating the masked beamspace subchannel $\mathbf{S}_{k}$ is equivalent to estimating the true beamspace subchannel $\mathbf{X}_{k}$. To recover the sparse vector $\mathbf{s}_{k}=\mathrm{vec}(\mathbf{S}_{k})$, the CS matrix corresponding to the linear model in (7) is defined as $\mathbf{A}_{k}\in\mathbb{C}^{M\times N^{2}}$, where

[TABLE]

The channel measurement $y_{k}[m]$ can be expressed using (7) and (8) as

[TABLE]

From (7), the channel measurement in (6) can be interpreted as a noisy version of the $(r_{k}[m],c[m])$ sample of the 2D-DFT of the masked subchannel $\mathbf{S}_{k}$. As a result, the CS problem for masked beamspace subchannel recovery is a partial 2D-DFT CS problem [18].

The map between $\mathbf{X}_{k}=\boldsymbol{\Lambda}^{-1}_{\mathbf{z}}\mathbf{S}_{k}\boldsymbol{\Lambda}^{-1}_{\mathbf{z}}$ and $\mathbf{S}_{k}$ is perfectly conditioned and invertible, due to the unimodular nature of the diagonal of $\boldsymbol{\Lambda}_{\mathbf{z}}$. Specifically, $|\mathbf{S}_{k}(r,c)|=|\mathbf{X}_{k}(r,c)|$ for every $k$, $r$ and $c$. Therefore, the masked beamspace subchannel $\mathbf{S}_{k}$ is sparse with the same locations of sparsity as that of the beamspace subchannel $\mathbf{X}_{k}$. Furthermore, the dynamics across the masked beamspace subchannels is smooth as the coefficients in the beamspace subchannels exhibit smooth variation across the subchannel index. Due to the structure preserving nature of the spectral mask, it is reasonable to solve for the masked beamspace subchannels instead of the original beamspace subchannels. From a computational complexity perspective, solving for $\{\mathbf{S}_{k}\}_{k=1}^{N_{\mathrm{RF}}}$ from the channel measurements is advantageous because CS matrix multiplications corresponding to (7) can be efficiently implemented using the fast Fourier transform (FFT) [18].

The problem of recovering masked beamspace subchannels is similar to that of recovering a sequence of temporally correlated angiogram images in magnetic resonance imaging (MRI) [10]. The sequence of sparse matrices $\{\mathbf{S}_{k}\}_{k=1}^{N_{\mathrm{RF}}}$ is analogous to the sequence of sparse angiogram images that are sampled across time. In both settings, the matrices exhibit smooth dynamics across the sequence index. Furthermore, image acquisition in MR-based angiogram imaging and the ZC-based channel acquisition obtain samples of the 2D-DFT of a sparse matrix. The equivalence between the two problems, and the success of DCS-AMP in recovering a sequence of smoothly varying MR images motivates its application for short range mmWave channel estimation.

III-B Masked beamspace recovery using DCS-AMP

In this section, we explain how the sequence $\{\mathbf{s}_{k}\}_{k=1}^{N_{\mathrm{RF}}}$, i.e., the vector versions of the masked beamspace subchannels, can be recovered from the channel measurements $\{\mathbf{y}_{k}\}_{k=1}^{N_{\mathrm{RF}}}$ using DCS-AMP [10]. Standard CS algorithms like approximate message passing (AMP) [22], when applied independently over the $N_{\mathrm{RF}}$ sub-problems, i.e., each sub-problem solves for a subchannel, do not leverage the smooth variation across $\{\mathbf{s}_{k}\}_{k=1}^{N_{\mathrm{RF}}}$ for channel estimation. The correlation across these vectors can be used to regularize the problem, and dynamic CS techniques like DCS-AMP [10] that account for such correlation can outperform standard AMP-based CS.

In DCS-AMP [10], the vectors $\{\mathbf{y}_{k}\}_{k=1}^{N_{\mathrm{RF}}}$ and $\{\mathbf{s}_{k}\}_{k=1}^{N_{\mathrm{RF}}}$ are considered as the realizations of random vectors $\{\boldsymbol{\mathsf{y}}_{k}\}_{k=1}^{N_{\mathrm{RF}}}$ and $\{\boldsymbol{\mathsf{s}}_{k}\}_{k=1}^{N_{\mathrm{RF}}}$. The random variables $\beta_{k}[n]\in\{0,1\}$ and $\eta_{k}[n]\in\mathbb{C}$ are used to model the support and amplitude of $\mathsf{s}_{k}[n]$, i.e., $\mathsf{s}_{k}[n]=\beta_{k}[n]\eta_{k}[n],\,\forall\,k,n$. The factor graph that illustrates the dependencies between the random variables is shown in Fig. 4. The elements of the random vectors are modeled using circular nodes, also called as variable nodes, in the factor graph. The rectangular nodes in Fig. 4 are called the factors. Factors in Fig. 4 are used to incorporate information about the random variables like sparsity and structural dependencies. For example, the dependency between $\mathsf{s}_{k}[n]$, $\beta_{k}[n]$ and $\eta_{k}[n]$ is captured using the factor $f_{k}[n]=\delta(\mathsf{s}_{k}[n]-\beta_{k}[n]\eta_{k}[n])$, where $\delta(\cdot)$ denotes the Dirac delta function [10]. The Dirac delta factor $f_{k}[n]$ ensures that the random variables satisfy $\mathsf{s}_{k}[n]=\beta_{k}[n]\eta_{k}[n]$. Each plane in the factor graph contains the variable nodes and factors corresponding to a particular subchannel. DCS-AMP performs information flows, i.e., probability distribution flows, within and across different planes to obtain subchannel estimates.

Channel measurements acquired by the subarrays provide information about the random vectors $\{\boldsymbol{\mathsf{s}}_{k}\}_{k=1}^{N_{\mathrm{RF}}}$. It can be observed from (9) that $y_{k}[m]$ is a realization of $\mathbf{A}_{k}(m,:)\boldsymbol{\mathsf{s}}_{k}$, i.e., the projection of the random vector $\boldsymbol{\mathsf{s}}_{k}$. The factor corresponding to $y_{k}[m]$ is denoted by $g_{k}[m]$ in Fig. 4. The factors $\{g_{k}[m]\}_{k,m}$ ensure fidelity of the message passing solution to the observed channel measurements. It can be observed from (7) and (9) that $y_{k}[m]$ depends on all the entries of $\mathbf{s}_{k}$ through a 2D-DFT equation. Therefore, $g_{k}[m]$, the factor corresponding to $y_{k}[m]$, is linked to all the $N^{2}$ elements of the random vector $\boldsymbol{\mathsf{s}}_{k}$ through a 2D-DFT. Specifically, $g_{k}[m]$ contains the likelihood function

[TABLE]

The factor graph built until this point is applicable to any generic linear regression problem.

Now, we describe the key components in DCS-AMP that incorporate sparsity and structural dependencies in the beamspace subchannels. For every $n\in\mathcal{I}_{N^{2}}$, the support correlation across the subchannels is incorporated by assuming that $\{\beta_{k}[n]\}^{N_{\mathrm{RF}}}_{k=1}$ follow a binary Markov process [10]. The $N^{2}$ different random processes are assumed to be independent with the same transition probabilities that characterize a binary Markov process. The correlation in the amplitudes of the beamspace subchannels is taken into account by assuming that $\{\theta_{k}[n]\}^{N_{\mathrm{RF}}}_{k=1}$ follows a Gauss-Markov process [10]. Similar to the binary Markov case, the $N^{2}$ different Gauss-Markov processes are assumed to be independent with the same correlation and variance parameters. The factors corresponding to the $n^{\mathrm{th}}$ binary Markov process and the Gauss-Markov processed are modeled using $\{b_{k}[n]\}_{k=1}^{N_{\mathrm{RF}}}$ and $\{q_{k}[n]\}_{k=1}^{N_{\mathrm{RF}}}$ in the factor graph. DCS-AMP assumes that both the Markov processes are in steady state. The steady state probability that $\beta_{k}[n]=1$ is denoted by $\epsilon$. The scalars $\kappa$, $\zeta$ and $\rho$ are used to denote the correlation coefficient, the mean and the variance of the steady state Gauss-Markov process corresponding to an $\eta_{k}[n]$. As $\mathsf{s}_{k}[n]=\beta_{k}[n]\eta_{k}[n]$, it can be observed that $\mathsf{s}_{k}[n]$ is [math] with a probability of $1-\epsilon$, or is distributed as a Gaussian with mean $\zeta$ and variance $\rho$ with a probability of $\epsilon$. As a result, $\mathsf{s}_{k}[n]$ follows a Bernoulli-Gaussian (BG) prior, which is a reasonable distribution to model sparse signals. For detailed mathematical expressions of the factors in Fig. 4, we refer the interested reader to [10].

Similar to standard message passing, DCS-AMP performs iterative message flows between the variable nodes and the factors to estimate the marginal posteriors of $\{\mathsf{s}_{k}[n]\}_{k,n}$ for the channel measurements $\{y_{k}[m]\}_{k,m}$. We explain the sequence of message flows using $f_{1}[1]$. The factors $b_{1}[1]$ and $q_{1}[1]$ send their messages, i.e., probability distributions, to $f_{1}[1]$ through $\beta_{1}[1]$ and $\eta_{1}[1]$. Then, $f_{1}[1]$ combines the received messages to generate a BG distribution. Standard AMP iterations between $\{g_{1}[m]\}^{M}_{m=1}$ and $\{\mathsf{s}_{1}[n]\}^{N^{2}}_{n=1}$ are performed for subchannel $1$ using the BG distributions received from the factors $\{f_{1}[n]\}^{N^{2}}_{n=1}$. At the end of these iterations, $f_{1}[1]$ receives Gaussian distributions from the factors $\{g_{1}[m]\}^{M}_{m=1}$. These messages are appropriately combined and sent to the next layer, i.e., factor graph for subchannel 2, through $b_{2}[1]$ and $q_{2}[1]$. Then, message passing is performed in layer $2$. The process of message flows within and across the layers is stopped when the $N_{\mathrm{RF}}^{\mathrm{th}}$ layer is updated. Such sequence of flows across the $N_{\mathrm{RF}}$ layers in Fig. 4 is called as a single forward pass in DCS-AMP. A backward pass starts from the $N_{\mathrm{RF}}^{\mathrm{th}}$ layer and performs message flows till the first layer is reached. After multiple forward and backward passes, DCS-AMP is expected to converge. At this point, the marginal posterior of $\mathsf{s}_{k}[n]$ for the channel measurements $\{\mathbf{y}_{k}\}^{N_{\mathrm{RF}}}_{k=1}$ is estimated by multiplying all the messages, i.e., probability densities, received by $\mathsf{s}_{k}[n]$. The mean of the estimated marginal posteriors gives the minimum mean-squared estimate of the sparse masked beamspace subchannels.

In practice, the prior parameters $\epsilon$, $\kappa$, $\xi$, $\zeta$ and $\rho$ are unkown apriori. To overcome this issue, DCS-AMP embeds the message passing algorithm within an Expectation-Maximization (EM) block that learns the parameters online [10]. For better subchannel recovery, the identical prior assumption can be waived to learn different set of parameters for different groups of subchannel locations. Inspired by the MRI example in [10], we develop a strategy to group the subchannel locations into two sets that correspond to potentially active and inactive coefficients. Let $\{\hat{\mathbf{S}}_{k,\mathrm{A}}\}^{N_{\mathrm{RF}}}_{k=1}$ be the subchannel estimates obtained using standard AMP-based CS over each subchannel, i.e., $\hat{\mathbf{S}}_{k,\mathrm{A}}$ is estimated from just $\mathbf{y}_{k}$. The factor graph corresponding to standard AMP is same as that in Fig. 4, but without the inter-subchannel factors $\{b_{k}[n]\}_{k,n}$ and $\{q_{k}[n]\}_{k,n}$. As a result, standard AMP-based CS does not exploit the smooth variation across subchannels and may result in poor subchannel estimates. Our DCS-AMP-based approach uses the subchannel estimates from standard AMP to determine the potentially active and inactive locations. We define an energy matrix corresponding to the subchannels recovered through standard AMP as $\hat{\mathbf{S}}_{\mathrm{E,A}}=\sum_{k=1}^{N_{\mathrm{RF}}}|\hat{\mathbf{S}}_{k,\mathrm{A}}|^{2}$. We define $\mathcal{S}_{1}\subset\mathcal{I}_{N}\times\mathcal{I}_{N}$ as a set that contains the locations of the largest coefficients in $\hat{\mathbf{S}}_{\mathrm{E,A}}$ that contribute to a $\delta_{\mathrm{E}}$ fraction of the total sum in $\hat{\mathbf{S}}_{\mathrm{E,A}}$. The set $\mathcal{S}_{2}$ contains the locations of the remaining coefficients such that $\mathcal{S}_{1}\cup\mathcal{S}_{2}=\mathcal{I}_{N}\times\mathcal{I}_{N}$. Note that the sets $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are invariant with the subchannel index $k$. For $\delta_{\mathrm{E}}$ close to $1$, it can be observed that $\mathcal{S}_{1}$ represents the set of potentially active coefficients. In this paper, we choose $\delta_{\mathrm{E}}=0.9$. The prior parameters in DCS-AMP are learned independently for the masked beamspace locations in $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$, using the EM algorithm [10].111An implementation of DCS-AMP-based subchannel recovery can be found on our page [23]; the code is based on the MRI example in [10].

We now discuss about the complexity of DCS-AMP with ZC-based subchannel acquisition. In a single pass of DCS-AMP, the complexity of message passing iterations corresponding to each subchannel is determined by matrix-vector multiplications involving $\mathbf{A}_{k}$ and $\mathbf{A}^{\ast}_{k}$ [10]. It can be observed from (8) that the matrix-vector products can be efficiently implemented using the fast Fourier transform, which has a complexity of $\mathcal{O}(N\mathrm{log}N)$. Let $T$ be the number of forward and backward passes that are carried out in DCS-AMP. In each pass of DCS-AMP, message passing iterations are performed for $N_{\mathrm{RF}}$ subchannels. As a result, the overall complexity of DCS-AMP is $\mathcal{O}(TN_{\mathrm{RF}}N\mathrm{log}N)$.

III-C Beam alignment and achievable rate

We now describe how to estimate subchannels using the estimates obtained with DCS-AMP. Let $\{\hat{\mathbf{S}}_{k,\mathrm{D}}\}^{N_{\mathrm{RF}}}_{k=1}$ be the masked beamspace estimates obtained using DCS-AMP. For a ZC-based spectral mask, the beamspace subchannel estimates $\{\hat{\mathbf{X}}_{k,\mathrm{D}}\}^{N_{\mathrm{RF}}}_{k=1}$ are found by inverting the spectral mask as $\hat{\mathbf{X}}_{k,\mathrm{D}}=\boldsymbol{\Lambda}^{-1}_{\mathbf{z}}\hat{\mathbf{S}}_{k,\mathrm{D}}\boldsymbol{\Lambda}^{-1}_{\mathbf{z}}$ [18]. The $k^{\mathrm{th}}$ subchannel estimate derived using DCS-AMP is then $\hat{\mathbf{H}}_{k,\mathrm{D}}=\mathbf{U}_{N}\hat{\mathbf{X}}_{k,\mathrm{D}}\mathbf{U}_{N}$.

To evaluate the performance of our DCS-AMP based approach, we consider a setting with $N_{\mathrm{RF}}$ STAs placed at distinct coordinates around the AP [24]. For simplicity, we assume that all the $N_{\mathrm{RF}}$ STAs are at the same distance of $d$ from the AP. We use $\mathbf{H}_{k,u}\in\mathbb{C}^{N\times N}$ to denote the subchannel matrix between the $u^{\mathrm{th}}$ STA and the $k^{\mathrm{th}}$ subarray at the AP. Under the perfect synchronization assumption, the $N_{\mathrm{RF}}$ STAs perform uplink pilot transmissions one at a time for channel estimation. We use $\hat{\mathbf{H}}_{k,u,\mathrm{D}}$ to denote the subchannel estimate corresponding to $\mathbf{H}_{k,u}$, that is obtained with DCS-AMP. At the end of the channel estimation phase, the AP allocates a single RF chain to a STA for data transmission. Without loss of generality, we assume that the $k^{\mathrm{th}}$ RF chain at the AP is allocated to the $k^{\mathrm{th}}$ STA. A graphical illustration of the simulation setting is shown in Fig. 5, where directional uplink data transmission is perfomed between a STA and a single associated RF chain at the AP.

Now, we explain how the directional beams at the AP and the STAs are determined from the estimated channel. For DCS-AMP-based channel estimation, the beamforming vectors at the $k^{\mathrm{th}}$ STA and the $k^{\mathrm{th}}$ subarray of the AP are computed from $\hat{\mathbf{H}}_{k,k,\mathrm{D}}$. We define $\mathbf{f}_{k,\mathrm{D}}$ as the $q$-bit phase quantized version of the top right singular vector of $\hat{\mathbf{H}}_{k,k,\mathrm{D}}$, i.e., the singular vector that corresponds to the maximum singular value of $\hat{\mathbf{H}}_{k,k,\mathrm{D}}$ [20]. Similarly, $\mathbf{w}_{k,\mathrm{D}}$ is defined as the conjugate of the $q$-bit phase quantized version of the top left singular vector of $\hat{\mathbf{H}}_{k,k,\mathrm{D}}$. As the beamformers are computed at the AP, the AP must signal the quantized phases to be used at the STAs using a control channel. Beamforming using a predefined codebook like the DFT can help reduce the signalling overhead. In this paper, however, we consider beamforming based on the quantized SVD-based approach to study the best case performance of DCSAMP. During uplink transmission, the STAs simultaneously transmit their data streams to the AP using $\{\mathbf{f}_{k,\mathrm{D}}\}^{N_{\mathrm{RF}}}_{k=1}$. We use $\mathbf{H}_{\mathrm{UL,D}}\in\mathbb{C}^{N_{\mathrm{RF}}\times N_{\mathrm{RF}}}$ to denote the effective multi-user uplink channel seen after analog combining with $\{\mathbf{w}_{k,\mathrm{D}}\}^{N_{\mathrm{RF}}}_{k=1}$. The $(i,j)^{\mathrm{th}}$ entry of $\mathbf{H}_{\mathrm{UL,D}}$ is given by $\mathbf{H}_{\mathrm{UL,D}}(i,j)=\mathbf{w}^{T}_{i,\mathrm{D}}\mathbf{H}_{i,j}\mathbf{f}_{j,\mathrm{D}}$. In the data transmission phase, the received signal at the $k^{\mathrm{th}}$ AP subarray comprises of the signal from the $k^{\mathrm{th}}$ STA in addition to interference from the remaining $N_{\mathrm{RF}}-1$ STAs. The off-diagonal terms in $\mathbf{H}_{\mathrm{UL,D}}$ indicate the amount of interference after analog beamforming using the DCS-AMP-based channel estimates. The interference can be mitigated by appropriate combining in the digital baseband of the AP.

In this paper, the classical MMSE beamformer is used at the AP to mitigate interference from other STAs. We assume that all the $N_{\mathrm{RF}}$ STAs transmit with the same amount of power. In such case, it is reasonable to assume that the received SNR, ignoring interference, is the same at each subarray. For DCS-AMP-based channel recovery, the signal-to-interference-plus-noise ratio (SINR) of the received signal from the $k^{\mathrm{th}}$ STA, after MMSE beamforming, is given by [24]

[TABLE]

The achievable rate corresponding to the $k^{\mathrm{th}}$ user is defined as $R_{k,\mathrm{D}}=\mathrm{log}_{2}(1+\rho_{k,\mathrm{D}})$.

IV Short range link configuration through geometry-aided message passing

In typical short range settings, information about the range of the transceiver distance, i.e., $d$ in Fig. 1, may be available. For example, it is reasonable to assume that $d\in[30\,\mathrm{cm},130\,\mathrm{cm}]$ in on-body wearable communication settings. In this section, we show that the array geometry of the transceiver together with the distribution on $d$, determine the correlation between the subchannels. The DCS-AMP-based approach, however, does not explicitly incorporate such geometry-based correlation. We construct factors to model such correlation and develop a parametric message passing algorithm that solves for the local AoAs.

To explain the idea underlying our method, we consider a single STA setting shown in Fig. 1. Beam alignment with the proposed approach is performed in two stages. In the first stage, pilot transmission is performed in the uplink by using a fixed transmit beam training vector at the STA. During this stage, the subarrays at the AP use different beam training vectors to acquire channel measurements for local AoA estimation. In the second stage, one of the subarrays at the AP performs directional pilot transmission using the local AoA estimated in the first stage. The STA uses different beam training vectors for AoA estimation in the second stage. Finally, the estimated AoA is used for beam alignment at the STA.

IV-A Geometry-aware local AoA estimation at the AP

The proposed geometry-aided message passing algorithm models the local AoAs, i.e., $\{\omega_{k}\}^{N_{\mathrm{RF}}}_{k=1}$, as realizations of random variables $\{\Omega_{k}\}^{N_{\mathrm{RF}}}_{k=1}$. Our message passing algorithm uses information from the channel measurements through angle-based likelihoods, and also models the correlation among the angles.

IV-A1 Angle-based likelihoods from channel measurements

Channel measurements for local AoA estimation are obtained by using different beam training vectors at the AP. The STA, however, uses a fixed ZC sequence $\mathbf{z}$ as the beam training vector throughout the local AoA estimation process at the AP. The channel response seen by the $k^{\mathrm{th}}$ subarray of the AP is then $\mathbf{h}_{\mathrm{AP},k}=\mathbf{H}_{k}\mathbf{z}$. As $\mathbf{H}_{k}=\alpha_{k}\mathbf{a}_{N}(\omega_{k})\mathbf{a}^{T}_{N}(\phi_{k})$ under the far field assumption for the subchannel, we can write $\mathbf{h}_{\mathrm{AP},k}=\alpha_{\mathrm{AP},k}\mathbf{a}_{N}(\omega_{k})$, where $\alpha_{\mathrm{AP},k}=\alpha_{k}\mathbf{a}^{T}_{N}(\phi_{k})\mathbf{z}$ is an unknown constant. The channel measurement acquired when the $k^{\mathrm{th}}$ subarray applies a beam training vector $\mathbf{w}_{k}[m]$ can be expressed as

[TABLE]

The channel measurement in (12) is a noisy projection of $\mathbf{a}_{N}(\omega_{k})$, up to an unknown scaling. The unknown scaling $\alpha_{\mathrm{AP},k}$ in $y_{\mathrm{AP},k}[m]$ is invariant with $m$, due to the use of a fixed ZC sequence at the STA array.

We describe the set of beam training vectors used at the AP to obtain channel measurements. We define $M_{\mathrm{AP}}$ as the number of channel measurements acquired by an AP subarray for local AoA estimation. The first $M_{\mathrm{AP}}-1$ projections are acquired by applying distinct circulant shifts $\{r_{k}[m]\}^{M_{\mathrm{AP}}-1}_{m=1}$ of the ZC sequence $\mathbf{z}$ to the $k^{\mathrm{th}}$ subarray of the AP, i.e., $\mathbf{w}_{k}[m]=\mathbf{J}_{r_{k}[m]}\mathbf{z}$. Then, the $M_{\mathrm{AP}}^{\mathrm{th}}$ channel measurement is acquired using $\mathbf{w}_{k}[M_{\mathrm{AP}}]=\mathbf{w}_{k}[1]\odot(1,-1,-1,\cdots-1)^{T}$. The $M_{\mathrm{AP}}^{\mathrm{th}}$ beam training vector is chosen differently to estimate the unknown complex gain $\alpha_{\mathrm{AP},k}$ in (12). With $\tilde{w}_{k}$ defined as the first entry of $\mathbf{w}_{k}[1]$, it can be observed from (12) that $y_{\mathrm{AP},k}[1]+y_{\mathrm{AP},k}[M_{\mathrm{AP}}]$ is a noisy version of $2\alpha_{\mathrm{AP},k}\tilde{w}_{k}$. The unknown gain $\alpha_{\mathrm{AP},k}$ is estimated as $\hat{\alpha}_{\mathrm{AP},k}=(y_{\mathrm{AP},k}[1]+y_{\mathrm{AP},k}[M_{\mathrm{AP}}])/2\tilde{w}_{k}$. We define a projection matrix $\boldsymbol{\Psi}_{k}\in\mathbb{C}^{M_{\mathrm{AP}}\times N}$ such that $\boldsymbol{\Psi}_{k}(m,:)=\mathbf{w}^{T}_{k}[m]$. The vector of channel measurements in (12) is then $\mathbf{y}_{\mathrm{AP},k}=\alpha_{\mathrm{AP},k}\boldsymbol{\Psi}_{k}\mathbf{a}_{N}(\omega_{k})+\mathbf{v}_{k}$. The gain compensated channel measurements are defined as $\tilde{\mathbf{y}}_{\mathrm{AP},k}=\mathbf{y}_{\mathrm{AP},k}/\hat{\alpha}_{\mathrm{AP},k}$. We ignore the error in estimating $\alpha_{\mathrm{AP},k}$ to conclude that $\tilde{\mathbf{y}}_{\mathrm{AP},k}$ is a realization of the random variable $\mathcal{N}_{c}(\boldsymbol{\Psi}_{k}\mathbf{a}_{N}(\omega_{k}),\sigma^{2}\mathbf{I}/|\hat{\alpha}_{\mathrm{AP},k}|^{2})$.

The likelihood of the local AoA $\omega_{k}$ is computed from the gain compensated channel measurements as

[TABLE]

A reasonable approach to estimate the local AoAs $\{\omega_{k}\}_{k=1}^{N_{\mathrm{RF}}}$ is to maximize the likelihood functions $\{p(\omega_{k})\}_{k=1}^{N_{\mathrm{RF}}}$ independently. We define the maximum likelihood (ML)-based local AoA estimate as $\hat{\omega}_{k,\mathrm{ML}}=\mathrm{arg\,max}\,p(\omega_{k})$. The ML-based approach, however, does not exploit the joint dependency among the local AoAs. As a result, it is possible that such an approach can result in local AoA estimates that are inconsistent with the geometry of the problem. For example, it can be observed from Fig. 1 that the combination $\omega_{1}=\pi/4$ and $\omega_{2}=-\pi/4$ is infeasible according to the geometry. Therefore, it is important to account for the dependencies among the local AoAs, that arise from the geometry, in the estimation technique.

IV-A2 Geometry-based dependencies for local AoA estimation

In this section, we model the dependencies among the local AoAs by constructing conditional probability densities of $\{\Omega_{k}\}_{k=1}^{N_{\mathrm{RF}}}$. Without loss of generality, we consider the local AoAs $\omega_{1}$ and $\omega_{2}$ to explain our construction. The distance between the midpoints of the STA and the first AP subarray is defined as $d_{1}$. We use $\ell_{k}$ to denote the distance between the midpoint of the $k^{\mathrm{th}}$ AP subarray and the midpoint of the AP array. For a particular $d$, the distance $d_{1}$ in Fig. 1 can be obtained from

[TABLE]

and, subsequently, the local AoA $\omega_{2}$ can be expressed as

[TABLE]

The relation between $\omega_{1}$ and $\omega_{2}$ is determined by $\omega_{2}=\mathcal{G}(\omega_{1},d)$, where $\mathcal{G}$ is a function defined by (14) and (15). It can be observed from (14) and (15) that there is a unique $\omega_{2}$ for a particular $\omega_{1}$ and $d$.

We construct factors for a realistic scenario in which the distance $d$ is unknown. Under the assumption that $d\in[d_{\mathrm{min}},d_{\mathrm{max}}]$, the range of $\omega_{2}$ for a given $\omega_{1}$ can be determined using (14) and (15). As the local AoAs are considered as realizations of random variables, we capture the local AoA dependencies using conditional distributions. We use $p_{\mathrm{g}}(\omega_{2}|\omega_{1})$ to denote the distribution of $\Omega_{2}$ conditioned on $\Omega_{1}=\omega_{1}$. We assume that $d$ is uniformly distributed in $[d_{\mathrm{min}},d_{\mathrm{max}}]$ to write

[TABLE]

where $D=d_{\mathrm{max}}-d_{\mathrm{min}}$, and the Dirac-delta function in (16) indicates $\omega_{2}=\mathcal{G}(\omega_{1},r)$. The procedure to determine $p_{\mathrm{g}}(\omega_{2}|\omega_{1})$, can be used to compute other conditional distributions, i.e., $\{p_{\mathrm{g}}(\omega_{k}|\omega_{n})\}_{k,n}$. It is important to note that the conditional distributions only depend on the geometry of the AP array and prior information about the transceiver distance. Therefore, the conditional distributions $\{p_{\mathrm{g}}(\omega_{k}|\omega_{n})\}_{k,n}$ can be computed offline, i.e., before acquiring channel measurements.

IV-A3 Message passing algorithm at the AP

We now describe the factor graph in Fig. 6 for our geometry-aided message passing algorithm. The factor graph consists of circular nodes, called as variable nodes, that represent the random variables $\{\Omega_{k}\}_{k=1}^{N_{\mathrm{RF}}}$. The rectangular nodes in Fig. 6, called as factors, contain the likelihoods or the conditional local AoA probabilities. Our algorithm performs message flows between the variable nodes and the factors to estimate the local AoAs. The message flows include a single forward pass and a single backward pass, and have a similar structure to the flows in DCS-AMP.

Now, we explain how messages are constructed in the forward pass. In the first step of the forward pass, the local AoA likelihood $p(\omega_{1})$ is sent to node $\Omega_{1}$. The message received by node $\Omega_{1}$ can be thought of a probability distribution of $\Omega_{1}$ that is believed by the factor based on the channel measurements. The node $\Omega_{1}$ then forwards $P^{\mathrm{fwd}}_{\mathrm{out}}(\omega_{1})=p(\omega_{1})$ to the geometry factor that contains $p_{\mathrm{g}}(\omega_{2}|\omega_{1})$, so that the belief about $\Omega_{2}$ can be computed. The message $P^{\mathrm{fwd}}_{\mathrm{in}}(\omega_{2})$, sent from the geometry factor to node $\Omega_{2}$, is defined as

[TABLE]

It can be observed from (17) that $P^{\mathrm{fwd}}_{\mathrm{in}}(\omega_{2})$ is $p_{\mathrm{g}}(\omega_{2}|\omega_{1})$ marginalized over $\omega_{1}$ using the belief from node $\Omega_{1}$, i.e., $P^{\mathrm{fwd}}_{\mathrm{out}}(\omega_{1})$. The message $P^{\mathrm{fwd}}_{\mathrm{in}}(\omega_{2})$ can be interpreted as a side information about $\Omega_{2}$ as it does not depend on the measurements acquired by the second subarray of the AP. This side information is generated using the channel measurements in the first subarray of the AP, and geometry-based dependency between $\omega_{1}$ and $\omega_{2}$. The node $\Omega_{2}$ combines information from the channel measurements, i.e., $p(\omega_{2})$, and the side information $P^{\mathrm{fwd}}_{\mathrm{in}}(\omega_{2})$ by constructing

[TABLE]

The message in (18) is then sent to the factor containing $p_{\mathrm{g}}(\omega_{3}|\omega_{2})$, which computes the belief about $\Omega_{3}$ using $p_{\mathrm{g}}(\omega_{3}|\omega_{2})$ and $P^{\mathrm{fwd}}_{\mathrm{out}}(\omega_{2})$. The message sent by the factor, i.e., $P^{\mathrm{fwd}}_{\mathrm{in}}(\omega_{3})$, is computed using an expression similar to (17). The process of message flows is executed until the last node $\Omega_{N_{\mathrm{RF}}}$ is reached. The messages in (17) and (18) can be generalized, by setting $\omega_{1}$ to $\omega_{k-1}$ and $\omega_{2}$ to $\omega_{k}$, to obtain recursive equations. The forward pass computes the message inflow $\{P^{\mathrm{fwd}}_{\mathrm{in}}(\omega_{k})\}_{k=2}^{N_{\mathrm{RF}}}$ using the generalized recursive equations.

The forward pass performs message flows in a particular direction, i.e., increasing $k$. As a result, it does not exploit information from subarrays $k,k+1,\cdots,N_{\mathrm{RF}}$, to generate side information about $\Omega_{k-1}$. Now, we describe the second component of our message passing algorithm, i.e., the backward pass. Message flows in the backward pass begin at node $\Omega_{N_{\mathrm{RF}}}$, and are performed until node $\Omega_{1}$ is reached. In the first step of the backward pass, the likelihood derived from the channel measurements, i.e., $P^{\mathrm{bwd}}_{\mathrm{out}}(\omega_{N_{\mathrm{RF}}})=p(\omega_{N_{\mathrm{RF}}})$, flows into the geometry factor containing $p_{\mathrm{g}}(\omega_{N_{\mathrm{RF}}-1}|\omega_{N_{\mathrm{RF}}})$. Then, the geometry factor computes a belief about $\Omega_{N_{\mathrm{RF}}-1}$ using $P^{\mathrm{bwd}}_{\mathrm{out}}(\omega_{N_{\mathrm{RF}}})$ and $p_{\mathrm{g}}(\omega_{N_{\mathrm{RF}}-1}|\omega_{N_{\mathrm{RF}}})$. The belief about $\Omega_{N_{\mathrm{RF}}-1}$ is denoted by $P^{\mathrm{bwd}}_{\mathrm{in}}(\omega_{N_{\mathrm{RF}}-1})$. The messages in the backward pass are computed using the recursive equations

[TABLE]

It can be observed from Fig. 6 that the nodes $\Omega_{1}$ and $\Omega_{N_{\mathrm{RF}}}$ do not receive any side information in the forward and backward passes. To validate the recursive equations, we define $P^{\mathrm{fwd}}_{\mathrm{in}}(\omega_{1})=\mathcal{U}(\omega_{1})$ and $P^{\mathrm{bwd}}_{\mathrm{in}}(\omega_{N_{\mathrm{RF}}})=\mathcal{U}(\omega_{N_{\mathrm{RF}}})$, where $\mathcal{U}(\omega)$ denotes a uniform distribution over $\omega$.

The final component of our message passing algorithm combines information from the forward and backward passes in a manner that is similar to DCS-AMP. It can be observed that node $\Omega_{k}$ receives side information in the form of $P^{\mathrm{fwd}}_{\mathrm{in}}(\omega_{k})$ and $P^{\mathrm{bwd}}_{\mathrm{in}}(\omega_{k})$ in the forward and backward passes. Furthermore, it also has access to the likelihood $p(\omega_{k})$. These three sources of information about $\Omega_{k}$ are combined by defining a new belief

[TABLE]

The local AoA estimate, using our geometry-aided message passing algorithm, is defined as $\hat{\omega}_{k}=\mathrm{arg\,max}\,p_{{}_{\mathrm{GMP}}}(\omega_{k})$.

Our implementation considers a discrete support for local AoAs to compute the message flows in the forward and backward passes. We assume an angular resolution of $\delta_{\omega}\pi$, for each of the local AoAs, to compute the integrals in (16), (17), and (19) using a discrete sum. The complexity of our algorithm is determined by the integration step in message passing. For example, it can be observed that the discrete version of (19) requires summing up $\mathcal{O}(1/\delta_{\omega})$ terms for each candidate $\omega_{k-1}$. As there are $\mathcal{O}(1/\delta_{\omega})$ candidates for $\omega_{k-1}$, the complexity of the integration step is $\mathcal{O}(1/\delta^{2}_{\omega})$. As $N_{\mathrm{RF}}-1$ such integrations are performed in each of the forward and backward passes, the overall complexity of our algorithm is $\mathcal{O}(N_{\mathrm{RF}}/\delta^{2}_{\omega})$. We would like to highlight that the proposed technique requires a single forward and backward pass, unlike DCS-AMP that performs multiple passes. The factor graph in both these techniques only model the correlation between adjacent subchannels. For example, the proposed technique does not include factors like $p_{\mathrm{g}}(\omega_{k}|\omega_{\ell})$ for $|k-\ell|>1$. Incorporating such factors can result in short cycles that may not be desirable from a message passing perspective [25].

IV-B Link configuration in a multi-user setting

We explain how our geometry-aided message passing algorithm is applied in a multi-user setting in Fig.  5. The signal structure used for our algorithm is similar to the one in the IEEE 802.11ad standard [5]. First, the STAs transmit pilots so that the AP can estimate the local AoAs using the proposed algorithm. A single STA is assumed to be active at any point of time during the uplink pilot transmission phase. Such an assumption can be relaxed by using orthogonal pilot sequences, in the time or frequency dimensions, at different STAs. Coordination among the STAs for uplink signalling can be achieved with scheduling using the control channel. We use $\omega_{k,u}$ to denote the local AoA at the $k^{\mathrm{th}}$ AP subarray corresponding to the $u^{\mathrm{th}}$ STA. For the geometry-aided message passing algorithm, the local AoA estimate corresponding to $\omega_{k,u}$ is denoted by $\hat{\omega}_{k,u,{\mathrm{GMP}}}$. The AP dedicates its $k^{\mathrm{th}}$ subarray to the $k^{\mathrm{th}}$ STA for downlink signalling and data communication. To achieve directional beamforming at the AP, a $q$-bit phase quantized version of $\overline{\mathbf{a}_{N}}(\hat{\omega}_{k,k,{\mathrm{GMP}}})$ is used as the beamforming vector at the $k^{\mathrm{th}}$ subarray of the AP. The subarrays at the AP perform downlink pilot transmissions using these directional beams so that the AoA at the STAs can be determined.

In the downlink pilot transmission phase, each STA acquires $M_{\mathrm{STA}}$ channel measurements by applying distinct random circulant shifts of $\mathbf{z}$ to its phased array. We use $\phi_{k,k}$ to denote the AoA at the $k^{\mathrm{th}}$ STA that is associated with the $k^{\mathrm{th}}$ subarray at the AP. Similar to the procedure in Sec. IV-A1, the $k^{\mathrm{th}}$ STA computes the likelihood function $p(\phi_{k,k})$ with the received channel measurements. The AoA estimate at the $k^{\mathrm{th}}$ STA is given by $\hat{\phi}_{k,k,\mathrm{ML}}=\mathrm{arg\,max}\,p(\phi_{k,k})$. Using the estimated local AoA, the $k^{\mathrm{th}}$ STA performs beam alignment by applying a $q$-bit phase quantized version of $\overline{\mathbf{a}_{N}}(\hat{\phi}_{k,k,\mathrm{ML}})$ to its phased array. It is important to note that the AP must feedback information about the best beamformers to the STAs in the DCS-AMP-based approach. The proposed method, however, does not require any such feedback as the beamformers are computed at the receiving end. A graphical illustration of multi-user data transmission with the estimated beamformers at the AP and the STAs is shown in Fig.  5. We define the $\mathbf{H}_{\mathrm{UL,GMP}}\in\mathbb{C}^{N_{\mathrm{RF}}\times N_{\mathrm{RF}}}$ as the multi-user uplink channel seen after configuring the phased arrays at the AP and the STA. The $(i,j)^{\mathrm{th}}$ entry of $\mathbf{H}_{\mathrm{UL,GMP}}$ is given by $\mathbf{H}_{\mathrm{UL,GMP}}(k,u)=\mathbf{a}^{\ast}_{N}(\hat{\omega}_{k,k,{\mathrm{GMP}}})\mathbf{H}_{k,u}\overline{\mathbf{a}_{N}}(\hat{\phi}_{u,u,\mathrm{ML}})$. Finally, an MMSE beamformer is used at the AP to cancel the interference that arises in the multi-user setting. The SINR and the achievable rate after MMSE beamforming are computed by using $\mathbf{H}_{\mathrm{UL,GMP}}$ in (11).

V Simulations

In this section, we evaluate the proposed techniques in a realistic communication setting that also has non-LoS components. We consider a hardware architecture in Fig. 1, where the AP is equipped with $N_{\mathrm{RF}}=4$ subarrays. We assume that each subarray at the AP has $N=16$ antennas. The resolution of the phase shifters at the AP and the STA is set as $q=2\,\mathrm{bits}$. The carrier frequency of the mmWave system is set to $60\,\mathrm{GHz}$, that corresponds to a wavelength of $\lambda=5\,\mathrm{mm}$. It can be observed that each subarray is of length $3.75\,\mathrm{cm}$. We assume that the length of the antenna array at the AP is $L_{\mathrm{AP}}=20\,\mathrm{cm}$. It can be shown that the spacing between consecutive subarrays at the AP is about $1\,\mathrm{cm}$. The AP considered in this setting is practical for wearable applications as its length, i.e., $L_{\mathrm{AP}}$, is comparable to the typical length of augmented reality or virtual reality headsets.

We consider $4$ STAs, where each STA is equipped with a phased array of $N=16$ elements. The antenna array at the STA occupies a length of $L_{\mathrm{STA}}=4\,\mathrm{cm}$, and can be realized in wearables like smart watches and fitness trackers. We assume that AP and the STA arrays are coplanar, and operate at the same polarization. For the $u^{\mathrm{th}}$ STA, the angles corresponding to $\gamma$ and $\theta$ in Fig. 1 are labeled as $\gamma_{u}$ and $\theta_{u}$. The location of the $u^{\mathrm{th}}$ STA relative to the AP is determined by the triplet $(d,\gamma_{u},\theta_{u})$. We assume that $|\gamma_{u}|\leq 75^{\mathrm{o}}$ for any $u$. The angles $\{\gamma_{u}\}_{u=1}^{N_{\mathrm{RF}}}$ are randomly chosen without repetition from the discrete set $\{-75^{\mathrm{o}}+15^{\mathrm{o}}i:i\in\mathcal{I}_{9}\}$. The orientation of the STAs relative to the AP array, i.e., $\{\theta_{u}\}_{u=1}^{N_{\mathrm{RF}}}$, are chosen at random from the discrete set $\{10^{\mathrm{o}}i:i\in\mathcal{I}_{18}\}$. The simulation results we report are for $1000$ random configurations of the STAs chosen according to this distribution; this corresponds to $1000$ different channel realizations. A graphical illustration of the short range mmWave setting is shown in Fig. 5 for a particular realization of $\{(\gamma_{u},\theta_{u})\}_{u=1}^{N_{\mathrm{RF}}}$.

Now, we describe the indoor short range communication scenario that is considered for our simulations. The AP and the STAs are placed in a room of floor area $5\,\mathrm{m}\,\times 5\,\mathrm{m}$ and a height of $3\,\mathrm{m}$, as shown in Fig. 7. The AP is placed on the sidewall given by $y=5\,\mathrm{m}$ at a height of $1.5\,\mathrm{m}$, and is oriented along the $x$-axis. All the $N_{\mathrm{RF}}$ STAs are placed at a distance of $d=80\,\mathrm{cm}$ from the AP, on the horizontal plane given by $z=1.5\,\mathrm{m}$. The center of the AP array is at $(2.5\,\mathrm{m},\,5\,\mathrm{m},\,1.5\,\mathrm{m})$. The $NN_{\mathrm{RF}}\times N$ channel matrix $\mathbf{H}$ is computed using the geometric-ray-based approach in [26]. Investigating the performance of our algorithm with more realistic channel models, which account for blockage and spatio-temporal channel dynamics [27], is an interesting research direction. For each combination of antennas at the AP and the STA, the geometric approach in [26] considers three kinds of propagation rays. The first kind is due to the direct path between the antennas, and the channel entry that results from such a ray is given in (2). The second type of propagation is determined by one bounce reflections in the environment. It is assumed that the floor is carpeted, and reflections from the floor are ignored due to high scattering by carpeted surfaces [28]. As a result, there are three one bounce reflections between a pair of antennas, i.e., one due to the ceiling and two due to the sidewalls. The third kind of propagation is determined by two bounce reflections and consists of $4$ rays. These rays arise due to reflections from a pair of sidewalls, and reflections from ceiling and sidewalls. In our simulations, the reflection parameters of the walls and the ceiling were set according to [26]. For a detailed description on the geometric short range channel model, we refer the reader to [26]. We observed that the coefficients of the beamspace subchannels corresponding to this indoor channel model exhibit sparsity and smooth variation across the subchannel index.

Channel measurements in the DCS-AMP-based approach and the geometry-aided message passing approach are acquired using circulant shifts of a ZC sequence $\mathbf{z}$. The root of the ZC sequence, i.e., $t$ in (5), is chosen as $9$ to ensure that $\mathbf{z}$ can be realized in the 2-bit phased array [20]. In the DCS-AMP-based method, the AP acquires $MN_{\mathrm{RF}}$ ZC-based projections of the channel matrix associated with each STA. Pilot transmissions are only performed in the uplink for link configuration using DCS-AMP. For the geometry-aided message passing appproach, $M_{\mathrm{AP}}$ and $M_{\mathrm{STA}}$ pilot transmissions are perfomed in the uplink and the downlink for each STA. In this paper, we use $M_{\mathrm{AP}}=M/2$ and $M_{\mathrm{STA}}=M/2$ so that the total number of channel measurements corresponding to each STA is same as that in the DCS-AMP-based approach. The resolution in the AoA domain, for the geometry-aided technique, is chosen as $1^{\mathrm{o}}$, i.e., $\delta_{\omega}=1/180$. Furthermore, the geometry factors were computed by assuming $d_{\mathrm{min}}=30\,\mathrm{cm}$ and $d_{\mathrm{max}}=130\,\mathrm{cm}$. The number of forward and backward passes in DCS-AMP is set to $T=12$. The reported results are for several random realizations of random circulant shifts at the AP and the STA, for each channel realization. The performance of the two techniques is evaluated in terms of the achievable rate per STA.

We compare the performance of the proposed link configuration algorithms with three different benchmarks. The first one is based on standard AMP-based CS, that solves for each sparse masked beamspace subchannel independently using the EM-BG-AMP algorithm [29]. The same ZC-based CS matrices are used for both DCS-AMP and standard AMP-based techniques. The second benchmark is based on maximum likelihood-based local AoA estimation at the AP. In this approach, the local AoA in a subarray is determined only from the channel measurements acquired in that subarray. The two benchmarks allow us to study the use of inter-subchannel factors in DCS-AMP, and geometry factors in the proposed geometry-aided message passing techniques. Finally, the third benchmark is based on the perfect CSI scenario in which the beamformers are derived from $\mathbf{H}_{k,k}$. The rate achieved with each of the three benchmarks is computed using the procedure in Sec. III-C.

As each subchannel has dimensions of $16\times 16$, it can be observed that beam alignment through conventional MIMO channel estimation requires $M=256$ pilot transmissions from every STA. The proposed techniques can be used to perform beam alignment with fewer channel measurements, i.e., $M<256$. In Fig. 8, we plot the empirical cumulative distribution function (CDF) of the achievable rate obtained with the proposed methods and the benchmarks. The achievable rate considered in our results is a concatenated version of the per-user rates. For $M=16$ pilot transmissions per STA and an SNR of $10\,\mathrm{dB}$, the CDFs in Fig. 8 indicate that the message passing-based techniques can be used to achieve a reasonable rate with sub-Nyquist channel measurements. DCS-AMP achieves a better rate over standard AMP as it exploits the smooth subchannel variation in addition to sparsity. It can be observed from Fig. 8 that the proposed geometry-aided message passing algorithm with $M_{\mathrm{AP}}=M_{\mathrm{STA}}=8$ performs better beam alignment when compared to the maximum likelihood-based approach, as it exploits additional local AoA dependencies using the geometry factors. The plot in Fig. 8 shows that the proposed message passing method outperforms DCS-AMP for short range link configuration. It is because our algorithm exploits the structure across the subchannels explicity using geometry factors, unlike the DCS-AMP-based approach that models subchannel variation using a Gauss-Markov process. For a particular STA and a channel realization, the likelihoods $\{p(\omega_{k})\}_{k=1}^{N_{\mathrm{RF}}}$ and the distributions $\{p_{{}_{\mathrm{GMP}}}(\omega_{k})\}_{k=1}^{N_{\mathrm{RF}}}$ are shown in Fig. 9. It can be noticed from Fig. 9 that geometry-aided message passing reduces the uncertainity in local AoAs when compared to the ML-based approach.

Now, we study the performance of the short range link configuration algorithms as a function of the SNR. The results in Fig. 10 are for $M=16$. An interesting observation from Fig. 10 is that DCS-AMP algorithm performs better than the proposed geometry-aided algorithm for SNRs of $0\,\mathrm{dB}$ and $5\,\mathrm{dB}$. At low SNRs, the poor performance of our geometry-based algorithm is mainly due to two reasons. First, the estimation errors in $\alpha_{\mathrm{AP},k}$, at low SNRs, impact the likelihoods in (13) and our message passing algorithm. Second, the proposed geometry-based technique uses a fixed transmit beam training vector at the STA, i.e., $\mathbf{z}$, for local AoA estimation at the AP. As a result, it is susceptible to nulls in the quasi-omnidirectional beam at the STA, especially at low SNRs. The DCS-AMP-based approach is robust to such nulls due to the use of different shifted ZC sequences $\mathbf{z}$ at the STA, i.e., the location of the nulls vary across different beam configurations at the STA.

We investigate how the achievable rate varies with the number of pilot transmissions using Fig. 11. It can be observed from Fig. 11 that DCS-AMP results in poor performance when compared to standard AMP for $M<14$. The transition effect in Fig. 11 is best understood by studying the impact of introducing the inter-subchannel factors to the factor graph of standard AMP; such factors ensure that the neighbouring subchannels are not significantly different from each other. Similar to any CS algorithm, the probability that standard AMP fails to recover some subchannels is higher for a lower $M$. When the number of subchannels that cannot be recovered by standard AMP is significantly smaller than $N_{\mathrm{RF}}$, the inter-subchannel factors in DCS-AMP “correct” the poor subchannel estimates obtained with standard AMP, by using information from other subchannels. The correction is likely to be successful if the neighbouring subchannels can be successfully recovered with standard AMP. For a small enough $M$ for which standard AMP fails to recover a significant number of subchannels, introducing the inter-subchannel factors in DCS-AMP can corrupt the successfully recovered subchannels. As a result, DCS-AMP performs poorly when compared to standard AMP for a sufficiently small $M$.

From Fig. 11, it can be observed that the proposed geometry-aided message passing algorithm achieves near optimal rates with fewer pilot transmissions, i.e., $M\approx 10$. The proposed algorithm has a lower complexity than DCS-AMP due to two reasons. First, our algorithm operates on a lower dimensional representation of the MIMO channel, i.e., using local AoAs, when compared to DCS-AMP that performs optimization over the full MIMO channel. Second, a single forward and a single backward pass is carried out in our algorithm, in comparison to DCS-AMP that uses multiple passes. As seen in Fig. 11, the geometry-based approach performs better than the ML-based approach as it exploits the dependencies among the local AoAs. It is important to note, however, that the geometry-aided message passing algorithm assumes a LoS scenario and may perform poorly in channels that have strong non-LoS components. In such settings, DCS-AMP algorithm can perform better than the geometry-based approach as it can also handle the non-LoS components.

VI Conclusions and Future work

Short range mmWave channels exhibit different structural properties when compared to the commonly studied far field channels. For instance, a typical short range LoS MIMO channel can have a rank that is higher than one. The differences in structural properties motivate the need to develop new channel structure-aware beam alignment techniques for short range channels. In this paper, we modeled the structure in short range mmWave channels by splitting the channel into several subchannels. We showed how exploiting subchannel correlation in addition to the sparse representation of each subchannel, can result in better beam alignment when compared to standard CS algorithms. We have also developed a low complexity message passing algorithm that ensures faithfulness of the solution to the problem geometry. Our results indicate that the algorithms based on dynamic CS and message passing are promising candidates to exploit structure in short range mmWave systems. Investigating the impact of antenna polarization and incorporating frequency selectiveness of channels are interesting research directions. In our future work, we will design efficient CS matrices for geometry aided-message passing.