Gridless Variational Bayesian Channel Estimation for Antenna Array Systems with Low Resolution ADCs

TL;DR

This paper introduces a grid-less variational Bayesian method for channel estimation in mm-wave antenna systems with low-resolution ADCs, effectively handling sparse angular channels without grid mismatch issues.

Contribution

It proposes the GL-QVBCE algorithm, a novel grid-less Bayesian approach that outperforms traditional methods and approaches the CRB in low-resolution ADC scenarios.

Findings

GL-QVBCE outperforms least squares in simulations.

The algorithm approaches the Cramér Rao bound asymptotically.

Effective for sparse mm-wave channel estimation with low-res ADCs.

Abstract

Employing low-resolution analog-to-digital converters (ADCs) coupled with large antenna arrays at the receivers has drawn considerable interests in the millimeter wave (mm-wave) system. Since mm-wave channels are sparse in angular dimensions, exploiting the structure could reduce the number of measurements while achieve acceptable performance at the same time. Motivated by the variational Bayesian line spectral estimation (VALSE) algorithm which treats the angles as random parameters, in contrast with previous works which confine the estimate to the set of grid angle points and induce grid mismatch, this paper proposes the grid-less quantized variational Bayesian channel estimation (GL-QVBCE) algorithm for antenna array systems with low resolution ADCs. Compared to the traditional least squares (LS) approach, numerical results show that GL-QVBCE performs significantly better and asymptotically approaches the Cram\`{e}r Rao bound (CRB).

Equations

$${\bar{\mathbf{y}}}_t=\left(\sum _{l=1}^L{\beta }_l\mathbf{a}({\varphi }_l)\right)x_t+{\mathbf{w}}_t,$$

$$\mathbf{a}({\varphi }_l)=[{\mathrm{e}}^{\mathrm{j}{m}_1\pi \sin ({\varphi }_l)},{\mathrm{e}}^{\mathrm{j}{m}_2\pi \sin ({\varphi }_l)},\cdots ,{\mathrm{e}}^{\mathrm{j}{m}_M\pi \sin ({\varphi }_l)}{]}^{\mathrm{T}},$$

$$\mathbf{h}=\sum _{l=1}^L{\beta }_l\mathbf{a}({\varphi }_l)=\mathbf{A}(\mathbf{\varphi })\mathbf{\beta }.$$

$${\mathbf{y}}_t=\mathcal{Q}(\Re \{{\bar{\mathbf{y}}}_t\})+\mathrm{j}\mathcal{Q}(\Im \{{\bar{\mathbf{y}}}_t\}),t=1,\cdots ,T.$$

$$\mathbf{Y}=\mathcal{Q}(\Re \{\mathbf{h}{\mathbf{x}}^{\mathrm{T}}+\mathbf{W}\})+\mathrm{j}\mathcal{Q}(\Im \{\mathbf{h}{\mathbf{x}}^{\mathrm{T}}+\mathbf{W}\}).$$

$$\mathbf{y}\triangleq \mathrm{vec}(\mathbf{Y})=\mathcal{Q}(\Re \{\mathbf{\Phi }\mathbf{h}+\mathbf{w}\})+\mathrm{j}\mathcal{Q}(\Im \{\mathbf{\Phi }\mathbf{h}+\mathbf{w}\}),$$

$$\mathbf{\Phi }=\mathbf{x}\otimes {\mathbf{I}}_M\in {\mathbb{C}}^{TM\times M}.$$

$$\bar{\mathbf{Y}}=\sum _{k=1}^K{\mathbf{h}}_k{\mathbf{x}}_k^{\mathrm{T}}+\mathbf{W}.$$

$$\bar{\mathbf{Y}}{\mathbf{x}}_k^{\ast }$$

$$={\mathbf{h}}_k{\mathbf{x}}_k^{\mathrm{T}}{\mathbf{x}}_k^{\ast }+\mathbf{W}{\mathbf{x}}_k^{\ast },$$

$$\mathbf{h}=\sum _{l=1}^N{\beta }_l\mathbf{a}({\theta }_l)\triangleq \mathbf{A}(\mathbf{\theta })\mathbf{\beta },$$

$$p(s_l)={\rho }^{s_l}(1-\rho {)}^{(1-s_l)},s_l\in \{0,1\}.$$

$$p({\beta }_l|s_l;\tau )=(1-s_l)\delta ({\beta }_l)+s_l\mathcal{C}\mathcal{N}({\beta }_l;0,\tau ).$$

$$p({\theta }_l)=\mathcal{V}\mathcal{M}({\theta }_l;{\mu }_{0,l},{\kappa }_{0,l})=\frac{1}{2\pi I_0({\kappa }_{0,l})}{\mathrm{e}}^{{\kappa }_{0,l}\cos (\theta -{\mu }_{0,l})},$$

$$\mathbf{z}=\mathbf{\Phi }\mathbf{h},$$

$$\mathbf{\Omega }=({\theta }_1,\dots ,{\theta }_N,(\mathbf{w},\mathbf{s})),$$

$$\mathbf{\eta }=\{\rho ,\tau \}$$

$$p(\mathbf{y},\mathbf{z},\mathbf{h},\mathbf{\Omega };\mathbf{\eta })$$

$$\times \prod _{i=1}^{N}p({\theta }_i)p({\beta }_i|s_i)p(s_i).$$

$${\hat{\mathbf{\eta }}}_{\mathrm{ML}}=\underset{\mathbf{\eta }}{\mathrm{argmax}}p(\mathbf{y};\mathbf{\eta })=\int p(\mathbf{y},\mathbf{z},\mathbf{h},\mathbf{\Omega };\mathbf{\eta })\mathrm{d}\mathbf{z}\mathrm{d}\mathbf{h}\mathrm{d}\mathbf{\Omega }.$$

$$(\hat{\mathbf{z}},\hat{\mathbf{h}},\hat{\mathbf{\Omega }})=\mathrm{E}[(\mathbf{z},\mathbf{h},\mathbf{\Omega })|\tilde{\mathbf{y}};{\mathbf{\eta }}_{\mathrm{ML}}],$$

$$p(\mathbf{z},\mathbf{h},\mathbf{\Omega }|\mathbf{y};{\hat{\mathbf{\eta }}}_{\mathrm{ML}})=\frac{p(\mathbf{z},\mathbf{h},\mathbf{\Omega },\mathbf{y};{\hat{\mathbf{\eta }}}_{\mathrm{ML}})}{p(\mathbf{y};{\hat{\mathbf{\eta }}}_{\mathrm{ML}})}$$

$$m_{\mathbf{z}\to {\delta }_{\mathrm{B}}}^t(\mathbf{z})\propto \frac{\mathrm{Proj}[m_{{\delta }_{\mathrm{B}}\to \mathbf{z}}^t(\mathbf{z})p(\mathbf{y}|\mathbf{z})]}{m_{{\delta }_{\mathrm{B}}\to \mathbf{z}}^t(\mathbf{z})}\triangleq \frac{\mathrm{Proj}[q_{\mathrm{C}}^t(\mathbf{z})]}{m_{{\delta }_{\mathrm{B}}\to \mathbf{z}}^t(\mathbf{z})},$$

$${\mathbf{z}}_{\mathrm{C}}^{\mathrm{post}}(t)=\mathrm{E}[\mathbf{z}|q_{\mathrm{C}}^t(\mathbf{z})],$$

$${\mathbf{v}}_{\mathrm{C},\mathrm{z}}^{\mathrm{post}}(t)=\mathrm{Var}[\mathbf{z}|q_{\mathrm{C}}^t(\mathbf{z})],$$

$$\mathrm{Proj}[q_{\mathrm{C}}^t(\mathbf{z})]=\mathcal{C}\mathcal{N}(\mathbf{z};{\mathbf{z}}_{\mathrm{C}}^{\mathrm{post}}(t),\mathrm{diag}({\mathbf{v}}_{\mathrm{C},\mathrm{z}}^{\mathrm{post}}(t))).$$

$$m_{\mathbf{z}\to {\delta }_{\mathrm{B}}}^t(\mathbf{z})\propto \frac{\mathcal{C}\mathcal{N}(\mathbf{z};{\mathbf{z}}_{\mathrm{C}}^{\mathrm{post}}(t),\mathrm{diag}({\mathbf{v}}_{\mathrm{C},\mathrm{z}}^{\mathrm{post}}(\mathrm{t})))}{\mathcal{C}\mathcal{N}(\mathbf{z};{\mathbf{z}}_{\mathrm{B}}^{\mathrm{ext}}(t),\mathrm{diag}({\mathbf{v}}_{\mathrm{B},\mathrm{z}}^{\mathrm{ext}}(\mathrm{t})))}\propto \mathcal{C}\mathcal{N}(\mathbf{z};{\mathbf{z}}_{\mathrm{C}}^{\mathrm{ext}}(t),\mathrm{diag}({\mathbf{v}}_{\mathrm{C},\mathrm{z}}^{\mathrm{ext}}(t))),$$

$${\mathbf{v}}_{\mathrm{C},\mathrm{z}}^{\mathrm{ext}}(t)={\left(\frac{1}{{\mathbf{v}}_{\mathrm{C},\mathrm{z}}^{\mathrm{post}}(t)}-\frac{1}{{\mathbf{v}}_{\mathrm{B},\mathrm{z}}^{\mathrm{ext}}(t)}\right)}^{-1},$$

$${\mathbf{z}}_{\mathrm{C}}^{\mathrm{ext}}(t)={\mathbf{v}}_{\mathrm{C},\mathrm{z}}^{\mathrm{ext}}(t)\odot \left(\frac{{\mathbf{z}}_{\mathrm{C}}^{\mathrm{post}}(t)}{{\mathbf{v}}_{\mathrm{C},\mathrm{z}}^{\mathrm{post}}(t)}-\frac{{\mathbf{z}}_{\mathrm{B}}^{\mathrm{ext}}(t)}{{\mathbf{v}}_{\mathrm{B},\mathrm{z}}^{\mathrm{ext}}(t)}\right),$$

$$\tilde{\mathbf{z}}=\mathbf{\Phi }\mathbf{h}+\tilde{\mathbf{w}},$$

$$m_{{\delta }_{\mathrm{B}}\to \mathbf{h}}^t(\mathbf{h})$$

$$\triangleq \frac{\mathrm{Proj}[q_{\mathrm{B}}^t(\mathbf{h})]}{m_{{\delta }_{\mathrm{A}}\to \mathbf{h}}^t(\mathbf{h})}.$$

$${\mathbf{\Sigma }}_{\mathrm{B},\mathrm{h}}^{\mathrm{post}}={\left({\mathbf{\Phi }}^{\mathrm{H}}\mathrm{diag}({\tilde{\mathbf{\sigma }}}_w^{-2}(t))\mathbf{\Phi }+\mathrm{diag}(1/{\mathbf{v}}_{\mathrm{A},\mathrm{h}}^{\mathrm{ext}}(t))\right)}^{-1},$$

Full text

keywords: Variational Bayesian inference, expectation propagation, quantization, grid-less, MMSE, LMMSE, CRB

I Introduction

As the spectrum available is limited in sub-6 GHz bands, millimeter wave (mm-wave) communications have drawn a great deal of attention as a potential technology for future cellular system [1]. Since mm-wave communication suffers from high attenuation, large antennas are often employed at the receiver side to focus their power. As the bandwidth scales up, the power consumption and implementation complexity of conventional analog-to-digital converters (ADCs) (8-14bits) increase significantly [2]. Consequently, low resolution ADCs (1-3bits) are proposed to be employed at the receiver side [3].

Since low-resolution ADCs incur severe distortion on the original received signal, traditional channel estimation algorithms may suffer significant performance degradation. Consequently, many studies have proposed novel signal processing algorithms to estimate the channel from heavily quantized samples [4, 5]. To reduce the training cost, a joint channel-and-data estimation method is proposed [6], although with high computation complexity.

The mm-wave channel is often characterized by a single dominant path and several other weaker paths [7, 8, 9]. In [10], a grid-based compressive sensing approach which confines the direction of arrivals (DOAs) to the set of grid angle points is proposed. Since the DOAs are continuous parameters, such an approach incurs basis mismatch and becomes more significant as the bit-depth decreases. To address the basis mismatch issue, a grid-less atomic norm based channel estimation approach is proposed under mixed one-bit antenna systems [11]. As shown in [11], exploiting this structure improves the channel estimation performance significantly over conventional channel estimators such as maximum likelihood (ML) and linear minimum mean squared error (LMMSE) estimators. Since the atomic norm based approach involves solving a semidefinite programming (SDP), the computation complexity is high for large antenna systems.

The angular models of the propagation channels for mm-wave communication share the similar structure as the line spectral estimation (LSE) or DOA problem. Consequently, advanced array processing methods may help address the channel estimation problem in wireless communications. In [12], a variational Bayesian LSE (VALSE) is proposed to perform the LSE. The VALSE treats the frequencies as random parameters, and it automatically estimates the model order (number of spectral), noise variance and the nuisance parameters of the prior distribution. In addition, VALSE also outputs the posterior probability density fucntion (PDF) of the frequencies and provides the uncertain degree of frequency estimates, in contrast with the prior work providing only point estimates. Later, multisnapshot VALSE (MVALSE) is developed to deal with the multiple measurement vector (MMV) [13]. In [14], VALSE is extended to solve the LSE from quantized samples. In this work, we derive the VALSE solution to the channel estimation problem.

From the algorithm point of view, many standard Bayesian algorithms [15, 16, 17, 18, 19, 20, 21] such as approximate message passing (AMP) [15], sparse Bayesian learning (SBL) [16], orthogonal approximate message passing (OAMP) [19], vector approximate message passing (VAMP) [21] have been proposed to deal with the standard linear model (SLM). For the generalized linear model (GLM), Bayesian algorithms [27, 29, 25, 28, 30, 31, 32, 26] such as generalized AMP (GAMP) [27], generalized expectation consistent signal recovery (GEC-SR) algorithm [28], and generalized VAMP [29] are developed. Later, a unified Bayesian inference framework is proposed, which demonstrates that the GLM can be solved via exchanging extrinsic information between the SLM and the minimum mean square error (MMSE) module [30]. Then, generalized AMP (Gr-AMP), generalized VAMP (Gr-VAMP), and generalized SBL (Gr-SBL) are developed. In addition, by discreting the DOAs into the grids, Gr-SBL is applied to solve the DOA from one-bit samples [31]. For the inference in deep networks, multi-layer VAMP (ML-VAMP) derived from expectation propagation [34] is proposed, and the mean-squared error performance of ML-VAMP is exactly predicted in a certain large system limit [33]. As shown later, the channel estimation in this paper can be viewed as a problem of estimating the line spectral which undergoes a linear transform followed by a componentwise nonlinear transform. As a result, we could utilize EP to design algorithms by implementing the respective modules and scheduling the messages between these modules 111For more details about EP and its relation to pmessage passing algorithms, please refer to [22, 23, 24, 21, 30, 26].

This paper proposes a grid-less quantized variational Bayesian channel estimation (GL-QVBCE) approach. The GL-QVBCE is designed from the module point of view. To evaluate the performance of the GL-QVBCE, we also use the additive quantization noise model (AQNM) and design GL-VBCE algorithm. Besides, the Cramèr-Rao bound (CRB) is derived under low resolution ADCs to be acted as the benchmarks of the proposed algorithm. Since the DOAs may be fixed across the pilots, while the complex amplitudes may be varied, sequential GL-QVBCE (Seq-GL-QVBCE) algorithm is proposed to track the channel. Finally, substantial numerical experiments are conducted to illustrate the effectiveness of the proposed approach, and investigate the factors on channel estimation performance, including signal to noise ratio (SNR), bit-depth, number of effective antennas, number of pilots, etc..

I-A Notation

For a matrix $\mathbf{A}$, let $[{\mathbf{A}}]_{i,j}$ or $A_{ij}$ denote the $(i,j)$th element of ${\mathbf{A}}$, and ${\rm diag}(\mathbf{A})$ returns a vector with elements being the diagonal elements of $\mathbf{A}$. Let $\Re\{\cdot\}$ and $\Im\{\cdot\}$ denote the real and imaginary part operator, respectively. Let $\otimes$ denote the Kronecker product operator. For a vector ${\mathbf{a}}$, let ${\rm diag}({\mathbf{a}})$ return a diagonal matrix with diagonal elements being $\mathbf{a}$. The symbol $t$ refers to the time or iteration index. Let ${\mathcal{S}}\subset\{1,\cdots,N\}$ be a subset of indices. For a square matrix ${\mathbf{A}}\in{\mathbb{C}}^{N\times N}$, let $\mathbf{A}_{{\mathcal{S}},{\mathcal{S}}}$ denote the submatrix by choosing both the rows and columns of $\mathbf{A}$ indexed by $\mathcal{S}$. Let ${(\cdot)}^{*}_{\mathcal{S}}$, ${(\cdot)}^{\rm T}_{\mathcal{S}}$ and ${(\cdot)}^{\rm H}_{\mathcal{S}}$ be the conjugate, transpose and Hermitian transpose operator of ${(\cdot)}_{\mathcal{S}}$, respectively. For a random vector $\mathbf{x}$ with probability density function (PDF) $p({\mathbf{x}})$, let ${\rm Proj}[p({\mathbf{x}})]$ denote the projection of $p({\mathbf{x}})$ onto Gaussian PDF with diagonal covariance matrix, where the means and variances are matched with that of $p({\mathbf{x}})$.

II System Model

Consider a pilot signal $x_{t}\in{\mathbb{C}}$ that impinges on a linear array with $M$ antennas. The received signal vector $\bar{\mathbf{y}}_{t}$ at the $t$th time-instant is expressed as [11]

[TABLE]

where $L$ is the number of rays, $\beta_{l}\in{\mathbb{C}}$ is the $l$th ray coefficient, ${\mathbf{a}}(\phi_{l})\in{\mathbb{C}}^{M}$ is the array manifold vector with $\phi_{l}\in[-\pi,\pi)$ as the DOAs of the $l$th ray, ${\mathbf{w}}_{t}\in{\mathbb{C}}^{M}$ is the additive white Gaussian noise with mean zero and covariance matrix $\sigma^{2}{\mathbf{I}}_{M}$. The array manifold vector ${\mathbf{a}}(\phi_{l})$ for the linear array is

[TABLE]

where ${\mathcal{M}}=\{m_{1},\cdots,m_{M}\}\subseteq\{0,1,\cdots,N-1\}$. Note that this model can be applied to general linear arrays.

Denote ${\mathbf{A}}({\bm{\phi}})=[{\mathbf{a}}(\phi_{1}),{\mathbf{a}}(\phi_{2}),\cdots,{\mathbf{a}}(\phi_{L})]\in{\mathbb{C}}^{M\times L}$ and ${\bm{\beta}}=[\beta_{1},\beta_{2},\cdots,\beta_{L}]^{\rm T}$. The channel vector ${\mathbf{h}}\in{\mathbb{C}}^{M}$ can be expressed as

[TABLE]

Suppose that ${\mathbf{y}}_{t}$ are quantized into a finite number of bits 222As will be shown later, the algorithm can be easily incorporated and extended to work with the mixed ADC system., i.e.,

[TABLE]

The goal of this paper is to estimate the channel $\mathbf{h}$ from the quantized observations $\{{\mathbf{y}}_{t}\}_{t=1}^{T}$. To estimate $\mathbf{h}$ accurately, the structure of $\mathbf{h}$ is exploited, i.e., the number of paths $L$, the DOAs $\{\phi_{l}\}_{l=1}^{L}$, and their complex gains $\{g_{l}\}_{l=1}^{L}$ are estimated.

Combining (1) and (4), and by defining ${\mathbf{Y}}=[{\mathbf{y}}_{1},\cdots,{\mathbf{y}}_{T}]\in{\mathbb{C}}^{M\times T}$ and the pilot vector ${\mathbf{x}}=[x_{1},\cdots,x_{T}]^{\rm T}\in{\mathbb{C}}^{T}$, one has

[TABLE]

Through vectorization and utilize ${\rm vec}({\mathbf{A}}{\mathbf{B}}{\mathbf{C}})=\left({\mathbf{C}}^{\rm T}\otimes{\mathbf{A}}\right){\rm vec}({\mathbf{B}})$, one obtains

[TABLE]

where

[TABLE]

In the following, efficient gridless quantized variational Bayesian channel estimation (GL-QVBCE) algorithm is designed to estimate the channel vector $\mathbf{h}$. Besides, gridless variational Bayesian channel estimation (GL-VBCE) algorithm can also be designed for the unquantized model.

It is worth noting that the GL-VBCE can be applied to the multiuser scenario. Consider an $K$ user scenario where $K=T$, let ${\mathbf{h}}_{k}\in{\mathbb{C}}^{M}$ denote the channel between the $k$th user and base station (BS) with the same angular structure as (1), let ${\mathbf{x}}_{k}\in{\mathbb{C}}^{T}$ denote the pilot of the $k$th user. Since $K=T$, one can construct orthogonal pilots, i.e., ${\mathbf{x}}_{k_{1}}^{\rm H}{\mathbf{x}}_{k_{2}}=0$, $\forall k_{1}\neq k_{2}$. By defining ${\bar{Y}}=[\bar{\mathbf{y}}_{1},\cdots,\bar{\mathbf{y}}_{T}]\in{\mathbb{C}}^{M\times T}$, we have

[TABLE]

For the $k$th user, we extract the $k$th user information by

[TABLE]

which has the form similar to (1). For the quantized setting, if the additive quantization noise model (AQNM) is adopted (see (62)), then we obtain a single user channel estimation problem similar to equation (9). While if we directly consider the nonlinear quantization effects, the proposed approach can still be applied. As shown in subsection III-B, the nonlinear model is iteratively approximated as a linear model with noise being heteroscedastic (different components having different variance), we still extract the $k$th user information similarly as (II) does and perform the channel estimation. The difference is that the extrinsic message between the modules needs to be redesigned. To avoid the obfuscation of the key features of our algorithms by intricate notations, we focus on the single user model (6) instead.

III GL-QVCBE pAlgorithm

This section develops the algorithm. First, the probabilistic formulation is introduced. Then, the gridless quantized variational Bayesian channel estimation (GL-QVBCE) algorithm consisting of several different modules are designed.

III-A Probabilistic Formulation

To estimate the channel vector $\mathbf{h}$ via exploiting the structures, we borrow the probabilistic formulation from [12]. For completeness, we introduce the probabilistic formulation in the following text.

First, we reparameterize the model by defining $\theta_{i}=\pi\sin\phi_{i}$ and $\theta_{i}\in[-\pi,\pi)$. We term $\theta_{i}$ as the $i$th frequency. Similarly, the number of paths $L$ is termed the number of frequencies. Since the number of frequencies $L$ is usually unknown, the number of frequencies is assumed to be $N$ [12], i.e.,

[TABLE]

where ${\mathbf{A}}({\bm{\theta}})=[{\mathbf{a}}({\theta}_{1}),\cdots,{\mathbf{a}}({\theta}_{N})]$ and $N$ satisfies $N>L$. Since the number of frequencies is $L$, the binary hidden variables ${\mathbf{s}}=[s_{1},...,s_{N}]^{\rm T}$ are introduced, where $s_{l}=1$ means that the $l$th frequency is active, otherwise deactive ($\beta_{l}=0$). The probability mass function of $s_{l}$ is

[TABLE]

Given that $s_{l}=1$, we assume that the channel path coefficient p${\beta}_{l}\sim{\mathcal{CN}}({\beta}_{l};0,\tau)$. Thus $(s_{l},{\beta}_{l})$ follows a Bernoulli-Gaussian distribution, that is

[TABLE]

From (11) and (12), it can be seen that the parameter $\rho$ denotes the probability of the $l$th component being active and $\tau$ is a variance parameter. The variable ${\bm{\theta}}=[\theta_{1},...,\theta_{N}]^{\rm T}$ has the prior PDF $p({\bm{\theta}})=\begin{matrix}\prod_{l=1}^{N}p(\theta_{l})\end{matrix}$. Generally, $p(\theta_{l})$ is encoded through the von Mises distribution [36, p. 36]

[TABLE]

where $\mu_{0,l}$ and $\kappa_{0,l}$ are the mean direction and concentration parameters of the prior of the $l$th frequency $\theta_{l}$, $I_{p}(\cdot)$ is the modified Bessel function of the first kind and the order $p$ [36, p. 348]. Without any knowledge of the frequency $\theta_{l}$, the uninformative prior distribution $p(\theta_{l})={1}/({2\pi})$ is used [12]. For more details please refer to [12, 13].

By defining

[TABLE]

the PDF $p({\mathbf{y}}|{\mathbf{z}})$ of $\mathbf{y}$ conditioned on $\mathbf{z}$ can be easily calculated through (4). Let

[TABLE]

be the set of all random variables and the model parameters, respectively. According to the Bayes rule, the joint PDF $p({\mathbf{y}},{\mathbf{z}},{\mathbf{h}},{\bm{\Omega}};{\bm{\eta}})$ is

[TABLE]

Given the above joint PDF (III-A), the type II maximum likelihood (ML) estimation of the model parameters $\hat{\bm{\eta}}_{\rm ML}$ is

[TABLE]

Then the minimum mean square error (MMSE) estimate of the parameters $({\mathbf{z}},{\mathbf{h}},{\bm{\Omega}})$ is

[TABLE]

where the expectation is taken with respect to

[TABLE]

Directly solving the ML estimate of $\bm{\eta}$ (18) or the MMSE estimate of $({\mathbf{z}},{\mathbf{h}},{\bm{\Omega}})$ (19) are both intractable. As a result, an iterative algorithm is designed in Section III.

III-B Module Design

The additional hidden variables ${\mathbf{z}}$ defined in (14) and two $\delta(\cdot)$ factor nodes in the factor graph are the key to design GL-QVBCE algorithms. pAccording to the factor graph Fig. 1 (a), the problem can be solved by exchanging information between each two modules, i.e., module A, module B and module C, where module A performs the variational line spectral estimation (VALSE) algorithm, module B performs the linear minimum mean square error estimation (LMMSE), and module C performs the componentwise minimum mean square error (MMSE) estimation 333For the mixed one-bit system, the difference is that only the one-bit measurements are needed to be input to the MMSE module C, and the extrinsic information from module C and the unquantized measurements are input to the LMMSE module..

First, we initialize $m_{\delta_{\rm A}\rightarrow{\mathbf{h}}}^{0}({\mathbf{h}})$ and $m_{\delta_{\rm B}\rightarrow{\mathbf{z}}}^{0}({\mathbf{z}})$ referring to Fig. 1, where the subscript is for the index of iteration. For the $t$th iteration with $1\leq t\leq T_{\rm o}$, the algorithm is detailed as follows.

III-B1 MMSE module

For the $t$th iteration, assume $m_{\delta_{\rm B}\rightarrow{\mathbf{z}}}^{t}({\mathbf{z}})={\mathcal{CN}}({\mathbf{z}};{\mathbf{z}}_{\rm B}^{\rm ext}(t);{\rm diag}({\mathbf{v}}_{\rm B,z}^{\rm ext}(t)))$. According to EP [34], the message $m_{{\mathbf{z}}\rightarrow\delta_{\rm B}}^{t}({\mathbf{z}})$ is calculated to be

[TABLE]

where $\propto$ denotes identity up to a normalizing constant. First, the MMSE estimate of $\mathbf{z}$ can be obtained, i.e.,

[TABLE]

where ${\rm E}[\cdot|q_{\rm C}^{t}({\mathbf{z}})]$ and ${\rm Var}[\cdot|q_{\rm C}^{t}({\mathbf{z}})]$ are the mean and variance operations taken componentwise with respect to the distribution $\propto q_{\rm C}^{t}({\mathbf{z}})$. As a result, ${\rm Proj}[q_{\rm C}^{t}({\mathbf{z}})]$ is 444Here the diagonal EP is used. It is numerically shown that the scalar EP which averages the noise variance yields significant performance degradation.

[TABLE]

Substituting (24) in (21), the message $m_{{\mathbf{z}}\rightarrow\delta_{\rm B}}^{t}({\mathbf{z}})$ from the variable node $\mathbf{z}$ to the factor node $\delta_{\rm B}({\mathbf{z}}-{\bm{\Phi}}{\mathbf{h}})$ is calculated as

[TABLE]

where ${\mathbf{z}}_{\rm C}^{\rm ext}(t)$ and ${\mathbf{v}}_{\rm C,z}^{\rm ext}(t)$ are [30]

[TABLE]

where $\odot$ denotes componentwise multiplication.

III-B2 LMMSE module

According to the factor node $\delta_{\rm B}({\mathbf{z}}-{\bm{\Phi}}{\mathbf{h}})$ and the extrinsic message $m_{{\mathbf{z}}\rightarrow\delta_{\rm B}}^{t}({\mathbf{z}})$ transmitted from module C to module B, we obtain a pseudo standard linear model 555For the multiuser model with equation (27), we extract the $k$th user information similarly as (II) does and channel estimation can be performed.

[TABLE]

where $\tilde{\mathbf{z}}={\mathbf{z}}_{\rm C}^{\rm ext}(t)$, $\tilde{\mathbf{w}}\sim{\mathcal{CN}}(\tilde{\mathbf{w}};{\mathbf{0}},{\rm diag}(\tilde{\bm{\sigma}}_{w}^{2}(t)))$, $\tilde{\bm{\sigma}}_{w}^{2}(t)={\mathbf{v}}_{\rm C,z}^{\rm ext}(t)$. Combing the message $m_{\delta_{\rm A}\rightarrow{\mathbf{h}}}^{t}({\mathbf{h}})={\mathcal{CN}}={\mathcal{CN}}({\mathbf{h}};{\mathbf{h}}_{\rm A}^{\rm ext}(t),{\rm diag}({\mathbf{v}}_{\rm A,h}^{\rm ext}(t)))$, we update the message $m_{\delta_{\rm B}\rightarrow{\mathbf{h}}}^{t}({\mathbf{h}})$ as

[TABLE]

First, it can be shown that $q_{\rm B}^{t}({\mathbf{h}})$ is Gaussian distributed with covariance matrix and mean being

[TABLE]

As a result, we project $q_{\rm B}^{t}({\mathbf{h}})$ as

[TABLE]

where

[TABLE]

As a result, $m_{\delta_{\rm B}\rightarrow{\mathbf{h}}}^{t}({\mathbf{h}})$ can be calculated as

[TABLE]

where

[TABLE]

III-B3 VALSE module

According to $m_{\delta_{\rm B}\rightarrow{\mathbf{h}}}^{t}({\mathbf{h}})$, we obtain another pseudo measurement model

[TABLE]

where $\tilde{\mathbf{h}}={\mathbf{h}}_{\rm B}^{\rm ext}(t)$, $\tilde{\mathbf{n}}\sim{\mathcal{CN}}(\tilde{\mathbf{n}};{\mathbf{0}},{\rm diag}({\mathbf{v}}_{\rm B,h}^{\rm ext}(t)))$. Now we run the VALSE algorithm [14] and calculate the posterior mean and variances of ${\mathbf{h}}$ as ${\mathbf{h}}_{\rm A}^{\rm post}(t)$ and ${\mathbf{v}}_{\rm A,h}^{\rm post}(t)$, respectively.

III-B4 From VALSE to LMMSE

We now update the message $m_{\delta_{\rm A}\rightarrow{\mathbf{h}}}^{t+1}({\mathbf{h}})={\mathcal{CN}}({\mathbf{h}};{\mathbf{h}}_{\rm A}^{\rm ext}(t+1),{\rm diag}({\mathbf{v}}_{\rm A,h}^{\rm ext}(t+1)))$, given by

[TABLE]

III-B5 From LMMSE to MMSE

We pupdate the message $m_{\delta_{\rm B}\rightarrow{\mathbf{z}}}^{t+1}({\mathbf{z}})={\mathcal{CN}}({\mathbf{z}};{\mathbf{z}}_{\rm B}^{\rm ext}(t+1);{\rm diag}({\mathbf{v}}_{\rm B}^{\rm ext}(t+1)))$. It can be seen that once this message is updated, the iteration is closed. The message is updated as follows:

[TABLE]

It can be calculated that

[TABLE]

where

[TABLE]

Let ${\mathbf{v}}_{\rm B}^{\rm post}(t+1)={\rm diag}({\bm{\Sigma}}_{\rm B}^{\rm post}(t+1))$. Then

[TABLE]

Utilizing (38), one has

[TABLE]

which closes the loop of the proposed Gridless quantized variational Bayesian channel estimation (GL-QVBCE) algorithm.

The computation complexity of the proposed algorithm is analyzed for a single outer iteration. For the VALSE module, its computation complexity is $O(MN+NL^{3})$ [14]. While for the LMMSE module, the computations are dominated by (29) and (40) involving a $M\times M$ matrix inversion, whose computation complexity is $O(M^{3})$. As for the componentwise MMSE module, its computation complexity is small. The overall computation complexity is $O(M^{3}+MN+NL^{3})$. As a comparison, the atomic norm based algorithm performing cvx [38] has a complexity of $O(M^{4}N^{2.5}+N^{4.5})$ [37]. Since the number of iterations is usually small for the GL-QVBCE, the computation complexity of atomic norm based algorithm is higher than that of the proposed approach.

III-C Extension to Sequential Channel Estimation

In some settings, the DOAs are fixed across the pilots, while the gains and phases are varied. As a result, performing the channel tracking is very important. Similar to [13], sequential GL-QVBCE (Seq-GL-QVBCE) is developed. Since GL-QVBCE outputs the posterior PDF of the frequencies, performing the sequential estimation is very natural.

Suppose that we have performed GL-QVBCE given the first pilot $x_{t}|_{t=1}$ and obtain the posterior PDF $p({\bm{\theta}}_{t}|{t=1})$ of $\bm{\theta}$. For the second pilot $x_{t}|_{t=2}$, we run the GL-QVBCE algorithm. For the VALSE module, we encode $p({\bm{\theta}}_{t}|{t=1})$ as the prior distribution of $\bm{\theta}$. Following the previous steps, Seq-GL-QVBCE can be developed.

The prior distribution obtained from the previous measurements may be too strong and inaccurate, which will deteriorate the estimation performance. As a result, we may use the damping operation to decrease the concentration parameter ${\bm{\kappa}}_{0,t}\in{\mathbb{R}}^{L}$ of the prior distribution $p({\bm{\theta}}_{t})$ 666For large $\kappa$, the variance of the frequency is approximately $1/\kappa$.. For the $(t-1)$th posterior PDF $p({\bm{\theta}}_{t-1})$ with concentration parameter ${\bm{\kappa}}_{t-1}$, we set

[TABLE]

where $0<\lambda\leq 1$.

IV Cramér Rao bound

Before designing the recovery algorithm, the performance bounds of unbiased estimators are derived, i.e., the Cramér Rao bound (CRB). Although the Bayesian algorithm is designed, the CRB acts as the performance benchmark of the algorithm. To derive the CRB, $L$ is assumed to be known, the frequencies ${\bm{\theta}}\in{\mathbb{R}}^{L}$ (or DOAs ${\bm{\phi}}\in{\mathbb{R}}^{L}$) and weights ${\bm{\beta}}\in{\mathbb{C}}^{L}$ are treated as deterministic unknown parameters. As for the quantizer $Q(\cdot)$, the quantization intervals are $\{(t_{b},t_{b+1})\}_{b=0}^{|{\mathcal{D}}|-1}$, where $t_{0}=-\infty$, $t_{{\mathcal{D}}}=\infty$, $\bigcup_{b=0}^{{\mathcal{D}}-1}[t_{b},t_{b+1})={\mathbb{R}}$. Given a real number $a\in[t_{b},t_{b+1})$, the representation is

[TABLE]

Note that for a quantizer with bit-depth $B$, the cardinality of the output of the quantizer is $|{\mathcal{D}}|=2^{B}$. Let $\bm{\kappa}$ denote the set of parameters, i.e., ${\bm{\kappa}}=[{\bm{\theta}}^{\rm T},{\mathbf{g}}^{\rm T},{\bm{\varphi}}^{\rm T}]^{\rm T}\in{\mathbb{R}}^{3L}$, where ${\bm{\beta}}={\mathbf{g}}{\rm e}^{{\rm j}{\bm{\varphi}}}$. The probability mass function (PMF) of the measurements $p({\mathbf{Y}}|{\bm{\kappa}})$ is

[TABLE]

Moreover, the PMFs of $\Re\{Y_{it}\}$ and $\Im\{Y_{it}\}$ are

[TABLE]

where $I_{(\cdot)}$ is the indicator function,

[TABLE]

The CRB is equal to the inverse of the Fisher information matrix (FIM) ${\mathbf{I}}({\bm{\kappa}})\in{\mathbb{R}}^{3L\times 3L}$

[TABLE]

To calculate the FIM, the following Theorem [39] is utilized.

Theorem 1

[39]** The FIM ${\mathbf{I}}({\bm{\kappa}})$ for estimating the unknown parameter ${\bm{\kappa}}$ is

[TABLE]

For a general quantizer, one has

[TABLE]

and

[TABLE]

For the unquantized system, the FIM is

[TABLE]

* *

According to Theorem 1, we need to calculate $\frac{\partial\Re\{Z_{it}\}}{\partial{\bm{\kappa}}}$ and $\frac{\partial\Im\{Z_{it}\}}{\partial{\bm{\kappa}}}$. In our setting, we have

[TABLE]

We need to calculate $\frac{\partial\Re\{Z_{it}\}}{\partial{\bm{\kappa}}}$ and $\frac{\partial\Im\{Z_{it}\}}{\partial{\bm{\kappa}}}$. By defining

[TABLE]

We have, for $l=1,\cdots,L$,

[TABLE]

Substituting (60), (53) and (54) in (52) and (55), respectively, one obtains the FIM for quantized and unquantized settings. The CRBs are the inverse of the corresponding FIMs, i.e., ${\rm CRB}({\bm{\kappa}})={\mathbf{I}}^{-1}({\bm{\kappa}})$ and ${\rm CRB}_{\rm unq}({\bm{\kappa}})={\mathbf{I}}_{\rm unq}^{-1}({\bm{\kappa}})$, respectively. The CRB of the frequencies, gains and phases are $[{\rm CRB}({\bm{\kappa}})]_{1:K,1:K}$, $[{\rm CRB}({\bm{\kappa}})]_{K+1:2K,K+1:2K}$, $[{\rm CRB}({\bm{\kappa}})]_{2K+1:3K,2K+1:3K}$, respectively, which will be used as the performance metrics.

V Numerical Simulation

Performance of the proposed GL-QVBCE is evaluated numerically. The simulation setup is similar to [11]. For the ULA of the channel vector ${\mathbf{h}}={\mathbf{A}}({\bm{\phi}}){\bm{\beta}}$, we generate the multipath channel including a line-of-sight (LoS) path and $L-1$ non-LoS paths. Let $P$ denote the total received power and $P=\sum\limits_{l=1}^{L}{\rm E}[|\beta_{l}|^{2}]$. For each run, the ray gains $\{g_{l}\}_{l=1}^{L}$ are generated from the complex normal distribution, where ${\bm{\beta}}={\mathbf{g}}{\rm e}^{{\rm j}{\bm{\varphi}}}$. In addition, the average powers of the LoS ray gains are set as ${\rm E}[g_{1}^{2}]=0.5P$ and $g_{1}$ is generated from ${\mathcal{N}}(\sqrt{0.45P},0.05P)$, whereas the remaining non-LoS paths gains are ${\rm E}[g_{l}^{2}]=0.5P/(L-1)$, $l=2,\cdots,L$ and are generated from ${\mathcal{N}}(\sqrt{0.45P/(L-1)},0.05P/(L-1))$. The phases $\bm{\varphi}$ are generated uniformly from $(-\pi,\pi)$. The DOAs $\{\phi_{l}\}_{l=1}^{L}$ are chosen uniformly from $[-\pi/2,\pi/2)$. All the results are averaged over $100$ Monte Carlo (MC) trials unless stated otherwise. The signal to noise ratio (SNR) is defined as

[TABLE]

As for the quantizer, zero threshold is chosen for 1-bit quantizer, while a uniform quantizer is chosen for multi-bit quantization. For ${\mathbf{Z}}={\mathbf{h}}{\mathbf{x}}^{\rm T}$, straightforward calculation shows that the variance $\sigma_{z}^{2}$ of the elements of $\mathbf{Z}$ is $\sigma_{z}^{2}=P$. The real and imaginary parts of $\mathbf{Z}$ are quantized separately, and the dynamic range of the quantizer is restricted to be $[-3\sigma_{z}/\sqrt{2},3\sigma_{z}/\sqrt{2}]$. For a uniform quantizer with bit-depth $B$, the quantizer step size $\Delta$ is $\Delta=3\sigma_{z}/2^{B-0.5}$. Then the variance of the additive quantization noise is $\sigma_{q}^{2}=2\Delta^{2}/12=3\sigma_{z}^{2}/2^{2B}$ (including the real and imaginary parts) [40]. It is assumed that the noise variance $\sigma^{2}$ is available.

Three additional algorithms are implemented to make performance comparison. The first is the proposed GL-VBCE which works under unquantized setting. The second is GL-VBCE-AQNM which uses the AQNM model

[TABLE]

where ${\mathbf{N}}_{q}$ denotes the additive quantization noise whose elements are i.i.d. and satisfies $[{\mathbf{N}}_{q}]_{i,t}={\mathcal{CN}}([{\mathbf{N}}_{q}]_{i,t};0,\sigma_{q}^{2})$. Consequently, the AQNM model is the same as that of the unquantized model, and we apply the GL-VBCE with the equivalent noise ${\mathbf{W}}_{\rm eq}$, whose elements are i.i.d. and satisfies $[{\mathbf{W}}_{\rm eq}]_{i,t}={\mathcal{CN}}([{\mathbf{W}}_{\rm eq}]_{i,t};0,\sigma_{w}^{2}+\sigma_{q}^{2})$. The last is the conventional LS approach works under unquantized setting. It can be easily shown that the NMSE of LS is

[TABLE]

i.e., only SNR or the number of pilots can improve the performance of LS. This phenomenon can also be validated by the ensuing numerical simulations.

V-A Channel Estimation Error versus Iteration

At first, the NMSEs of the algorithms versus the iteration are presented and shown in Fig. 2. It can be seen that all the algorithms converge in several iterations. From (63), the NMSE of LS is $-10\log 2\approx-3$ dB due to ${\rm SNR}=0$ dB and $T=2$. All the algorithms utilizing channel structure works better than the traditional LS approach. Under 1-bit quantization, the performance of GL-QVBCE is better than that of GL-VBCE-AQNM. While for 2-bit quantization, the performance of GL-QVBCE is comparable to that of GL-VBCE-AQNM, and is close to GL-VBCE. Since GL-QVBCE involves the nonlinear MMSE, GL-VBCE-AQNM is preferable under 2-bit quantization in this setting.

V-B Comparisons of Channel Estimators versus SNR

This subsection investigates the channel estimation performance versus SNR. Parameters are set as follows: $N=M=64$, $L=2$, $T=2$. The performance of channel estimation error is presented in Fig. 3. The NMSEs of GL-VBCE and LS decreases linearly with $\rm SNR$. For GL-QVBCE under two bit quantization, the NMSE decreases as SNR increases. When SNR increases to $10$ dB, the NMSE is saturated. For the GL-VBCE-AQNM and GL-QVBCE under one bit quantization, the NMSEs first decrease as SNR increases. When the SNR exceeds a certain value, the NMSEs begin to increase. The phenomenon from GL-VBCE-AQNM demonstrates that AQNM is more preferable when SNR is low.

V-C Comparisons with CRB

The performances of the various algorithms against the CRB are investigated. Here we fix the two DOAs as $\phi=[-30^{\circ},60^{\circ}]^{\rm T}$, set the corresponding complex amplitudes as ${\bm{\beta}}=[0.8{\rm e}^{-{\rm j}0.3\pi},0.6{\rm e}^{{\rm j}0.2\pi}]^{\rm T}$. The pilot length is $T=2$. $N=M=96$. Here $P(\hat{L}=L)$ denotes the probability of both the model order $K$ is correctly estimated and ${\rm NMSE}(\hat{\mathbf{h}})\leq-5{\rm dB}$. As for the MSEs of ray-gain estimation, ray-phase estimation and DOA estimation, we only average the case when the two conditions are satisfied. For clarity, the performances of GL-VBCE-AQNM are not presented here. Fig. 4 presents the NMSEs versus SNR under this setting. The phenomena are basically the same as shown in subsection V-B.

The MSEs of ray-gain estimation, ray-phase estimation, DOA estimation against the CRB and the model order probability are presented in Fig. 5. It can be seen that all the algorithms except under one-bit quantization improve as SNR increases. On the one hand, noise is beneficial for recovering the magnitude information as magnitude information is lost under noiseless one-bit quantization. On the other hand, noise deteriorates the estimation performance as noise is unknown and adds uncertainty on the estimation. Thus it is expected that when SNR is low, the NMSE of the gain decreases with SNR, and when SNR is above a threshold, the MSE of the gain increases with SNR, as shown in Fig. 5 under one-bit quantization. For the phase and frequency estimation, the MSEs approach their CRBs under two bit quantization and unquantized setting. For the model order estimation probability, it increases as SNR increases for two bit quantization and unquantized setting. For one-bit quantization, when SNR exceeds $5$ dB, it decreases. The results demonstrate that for one-bit quantization, GL-QVBCE is more preferable under low SNR scenario.

V-D Effects of Bit-depth

This experiment is conducted to illustrate the effect of bit-depth on the channel estimation and results is shown in Fig. 6. It can be seen that increasing the bit-depth $B$ improves the channel estimation performance. As bit-depth increases, the performance gains decreases. From the figure, under 2-bit and 3-bit quantization, the performance gap between GL-QVBCE and GL-VBCE is about $3$ dB and $1$ dB. When $B\geq 4$, both GL-QVBCE and GL-VBCE-AQNM achieve almost identical performance with GL-VBCE under quantized setting.

V-E Effects of Number of Epffective Antennas

Increasing the number of effective antennas will improve the performance of all the algorithms except LS. Fig. 7 shows the NMSE of the channel versus the number of effective antennas $M$. It can be seen that under ${\rm SNR}=0$ dB, the NMSE of the LS approach is still [math] dB as $M$ increases, which is consistent with the results (63). While for the GL-QVBCE and GL-VBCE approach, their performances improves with increasing $M$. Under 1-bit quantization, GL-QVBCE performs better than that of the AQNM model. While for 2-bit quantization, their performances are close.

V-F Effects of Multipath Numbers

The proposed approach utilizes the channel structure and thus improve the performance over traditional channel estimation algorithms. Fig. 8 presents the MSE versus the number of multipath. It can be seen that performance of LS approach is stable irrespective of multipath numbers $L$, as (63) reveals. For all the algorithms exploiting the channel structure, their performance deteriorates as the number of multipath increases.

V-G Pilot Length

An obvious fact is that the performance of the channel estimation algorithms should be improved by increasing the pilot length and Fig. 9 validates this. In addition, it also shows similar phenomenon presented in subsection V-A.

V-H Sequential channel estimation

Consider a scenario where the DOAs are fixed across the pilots, while their complex coefficients are completely independent. In this setting, Seq-GL-QVBCE and Seq-GL-VBCE algorithm are implemented to perform sequential channel estimation. As for the GL-QVBCE and GL-VBCE, we calculate the NMSE of channel for each pilot, separatively. We set $\lambda=0.1$ (44) to damp the concentration parameter of the prior distribution. Results are presented in Fig. 10. It can be seen that both Seq-GL-QVBCE and Seq-GL-VBCE perform better than those of GL-QVBCE and GL-VBCE, respectively.

VI Conclusion

This paper proposes the GL-QVBCE and GL-VBCE algorithm exploiting channel structures to estimate the channel under quantized and unquantized setting. GL-QVBCE treats the DOAs as random parameters and performs the variational Bayesian estimation. When the DOAs are fixed across the pilots, Seq-GL-QVBCE and Seq-GL-VBCE are developed to perform sequential estimation. Numerical results are conducted to illustrate the effectiveness of the proposed algorithm.