# Channel Estimation for Quantized Systems based on Conditionally Gaussian
Latent Models

Benedikt Fesl, , Nurettin Turan, ,
Benedikt Böck, , and Wolfgang Utschick, This work was partly funded by Huawei
Sweden Technologies AB, Lund. The authors are with Lehrstuhl für Methoden der
Signalverarbeitung, Technische Universität München, 80333 München, Germany
(e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>utschick@tum.de).

###### Abstract

This work introduces a novel class of channel estimators tailored for coarse
quantization systems. The proposed estimators are founded on conditionally
Gaussian latent generative models, specifically Gaussian mixture models
(GMMs), mixture of factor analyzers (MFAs), and variational autoencoders
(VAEs). These models effectively learn the unknown channel distribution
inherent in radio propagation scenarios, providing valuable prior information.
Conditioning on the latent variable of these generative models yields a
locally Gaussian channel distribution, thus enabling the application of the
well-known Bussgang decomposition. By exploiting the resulting conditional
Bussgang decomposition, we derive parameterized linear minimum mean square
error (MMSE) estimators for the considered generative latent variable models.
In this context, we explore leveraging model-based structural features to
reduce memory and complexity overhead associated with the proposed estimators.
Furthermore, we devise necessary training adaptations, enabling direct
learning of the generative models from quantized pilot observations without
requiring ground-truth channel samples during the training phase. Through
extensive simulations, we demonstrate the superiority of our introduced
estimators over existing state-of-the-art methods for coarsely quantized
systems, as evidenced by significant improvements in mean square error (MSE)
and achievable rate metrics.

###### Index Terms:

Channel estimation, generative latent model, coarse quantization, Bussgang
theorem, covariance recovery.

##  I Introduction

©This work has been submitted to the IEEE for possible publication. Copyright
may be transferred without notice, after which this version may no longer be
accessible.

Massive multiple-input multiple-output (MIMO) and millimeter wave (mmWave)
systems enable the ever-increasing requirements of bandwidth and throughput in
wireless communications. However, deploying a large number of high-precision
analog-to-digital converters (ADCs) for each antenna’s radio frequency (RF)
chain with bandwidths sufficient for mmWave systems is unaffordable in terms
of cost and power consumption [1, 2]. One of the most direct and promising
ways in order to solve the power consumption bottleneck and achieve high
energy efficiency is to use low-resolution ADCs at the base station (BS). In
recent years, considerable research efforts have been devoted to analyzing the
performance of low-resolution quantization systems [3, 4, 5, 6]. Remarkably,
although the low-resolution quantization causes nonlinear distortions at the
receiver, the capacity is not severely reduced, especially at low signal-to-
noise ratios (SNRs) [7].

In order to realize the mentioned favorable characteristics in practical low-
resolution systems of the next generation of cellular systems (6G), accurate
channel estimation is a crucial task. However, the severe nonlinearity of the
ADCs degrades the performance of conventional channel estimation algorithms
[1, 2]; thus, it is necessary to design novel channel estimators that provide
good performance together with reasonable complexity and robustness in
quantized systems.

Various channel estimation algorithms for coarsely quantized systems have been
proposed in recent years. In [8], least squares (LS) estimation is considered,
which is computationally simple but results in rather poor estimation quality.
An iterative channel estimation technique based on the expectation-
maximization (EM) algorithm is proposed in [9], which exhibits limitations due
to high complexity and convergence to local optima. Iterative maximum
likelihood (ML) methods are investigated in [10, 11]; however, they typically
require a large number of pilot signals, resulting in an unaffordable
signaling overhead [2]. MMSE channel estimation approaches in the case of a
Gaussian channel together with uncorrelated channel entries are studied in
[12, 13]; unfortunately, the MMSE estimator has no closed-form in the general
case and is intractable to compute [14]. In [15, 16, 17], joint channel
estimation and decoding is investigated, where payload data is used to assist
in channel estimation. Due to the iterative nature of the approaches, these
methods are considered to have too high complexity for commercial massive MIMO
systems [2].

The works in [18, 19, 20] take into account the sparsity of wireless channels
and utilize compressive sensing (CS) approaches such as iterative hard
thresholding and generalized approximate message passing (GAMP); the main
disadvantages thereby are the sensitivity concerning the (estimated) sparsity
level and the high complexity of the iterative procedure. A well-known linear
MMSE channel estimator is based on Bussgang’s theorem [21]. In [22, 23], the
estimator is derived for the one-bit as well as multi-bit quantization case.
It has the advantage that only a few pilot observations are necessary to
achieve a good channel estimation quality. The critical prerequisite of the
Bussgang estimator is a Gaussian distributed channel with known second-order
statistics, which seriously limits the estimator’s application.

Recently, deep learning-based approaches were investigated for channel
estimation in quantized systems [24, 25, 26, 27]. Although providing good
performance for low numbers of pilots, the approaches generally lack
generalization ability with respect to different numbers of quantization bits,
pilot signals, antennas, and SNRs. Additionally, accumulating a large
representative training dataset consisting of perfect channel state
information (CSI) samples is necessary, which may be costly to acquire in
practical systems, especially with coarse quantization. In recent works,
conditionally Gaussian latent generative models were utilized in order to
learn the underlying unknown channel distribution of a radio propagation
environment and leverage this prior information to design wireless
communication functionalities [28], especially channel estimators for high-
resolution systems [29, 30, 31, 32, 33, 34]. As mentioned before, the highly
nonlinear distortion of low-resolution ADCs makes it impractical to directly
utilize the mentioned channel estimation techniques in coarsely quantized
systems.

Contributions: This article presents novel contributions for channel
estimation in coarsely quantized systems utilizing conditionally Gaussian
latent generative models. Specifically, we investigate the applicability of
three different models: GMM, MFA, and VAE for learning the unknown underlying
channel distribution within a BS cell. Through this, we achieve a
conditionally Gaussian representation of the learned channel distribution. By
inferring the latent variable of the presented generative approaches based on
a quantized pilot observation, we can efficiently determine the parameters of
the conditional Gaussian distribution. By applying Bussgang’s theorem, we
derive a parameterized conditional Bussgang estimator, which serves as the
linear MMSE estimator for the given latent encoding. Importantly, our work
showcases the diverse strengths of the presented models concerning estimation
performance, computational complexity, and memory overhead.

To facilitate practical feasibility, we propose training adaptations that
allow us to learn the channel distribution solely from quantized pilot
observations as training data, which eliminates the need for perfect CSI in
the training stage. This feature enables us to use pilot observations
collected during the regular BS operation for training. In the case of the
GMM, we introduce a novel covariance recovery method that serves as an
unbiased and consistent estimator of the unquantized input covariance matrix,
using only quantized samples. Additionally, for the VAE, we adapt the evidence
lower bound (ELBO) loss function through model-based insights, specifically
accounting for the nonlinear quantization process.

Extensive simulations based on various realistic channel models demonstrate
the outstanding performance of the proposed conditionally Gaussian latent
generative modeling-aided estimators, especially in setups with low numbers of
pilot observations, surpassing state-of-the-art methods.

Notation: The standard Gaussian cumulative distribution function (CDF) is
denoted by
${\operatorname{\Phi}}(x)=\int_{-\infty}^{x}\tfrac{1}{\sqrt{2\pi}}\exp(-\tfrac{x^{2}}{2}){\operatorname{d}}x$
and the error function is given as
$\operatorname{erf}(x)=\tfrac{2}{\sqrt{\pi}}\int_{0}^{x}\exp(-t^{2}){\operatorname{d}}t$.
We denote the indicator function as $\chi(x\in\mathcal{A})$, which returns one
if $x\in\mathcal{A}$ and zero otherwise. The $n$th entry of a vector is
denoted by $[{\bm{x}}]_{n}$. The diagonals and off-diagonals of a matrix are
denoted by $\operatorname{diag}({\bm{A}})$ and
$\operatorname{nondiag}({\bm{A}})={\bm{A}}-\operatorname{diag}({\bm{A}})$.

##  II Preliminaries

####  II-1 System Model

We consider the uplink transmission of $P$ pilot signals from a single-antenna
mobile terminal (MT) to an $N$-antenna BS which operates ADCs with $B$
quantization bits. The quantized receive signal is therefore written as
${\bm{R}}=Q_{B}({\bm{Y}})=Q_{B}({\bm{h}}{\bm{a}}^{\operatorname{T}}+{\bm{N}})$,
where ${\bm{R}}=[{\bm{r}}_{1},\dots,{\bm{r}}_{P}]\in\mathbb{C}^{N\times P}$
contains the $P$ quantized receive signals as columns,
${\bm{Y}}\in\mathbb{C}^{N\times P}$ describes the unquantized receive signal,
${\bm{h}}\in\mathbb{C}^{N}$ denotes the wireless channel,
${\bm{a}}\in\mathbb{C}^{P}$ is the pilot vector which fulfills the power
constraint $\|{\bm{a}}\|_{2}^{2}=P$,
${\bm{N}}=[{\bm{n}}_{1},\dots,{\bm{n}}_{P}]\in\mathbb{C}^{N\times P}$ is
additive white Gaussian noise (AWGN) with
${\bm{n}}_{i}\sim\mathcal{N}_{\mathbb{C}}({\bm{0}},\sigma^{2}\operatorname{\mathbf{I}})$,
and $Q_{B}$ denotes the $B$-bit quantization function, which is discussed
below. By column-wise vectorization, the system model can be written as

$\displaystyle{\bm{r}}=Q_{B}({\bm{y}})=Q_{B}({\bm{A}}{\bm{h}}+{\bm{n}})\in\mathbb{C}^{NP}$
(1)

with ${\bm{r}}=\operatorname{vec}({\bm{R}})$,
${\bm{y}}=\operatorname{vec}({\bm{Y}})$,
${\bm{n}}=\operatorname{vec}({\bm{N}})$, and
${\bm{A}}={\bm{a}}\otimes\operatorname{\mathbf{I}}$. Note that an extension to
a multi-user setup can, in principle, be straightforwardly achieved by
stacking the channels of all users, cf., e.g., [22]; however, the analysis of
multi-user systems is out of the scope of this work. By normalizing the
channels as ${\operatorname{E}}[\|{\bm{h}}\|_{2}^{2}]=N$, the SNR of the
quantizer input is defined as $\text{SNR}=1/\sigma^{2}$.

Typically, several pilot observations are required to achieve reasonable
channel estimation performance in coarsely quantized systems. For the case of
one-bit quantization where the amplitude information is eliminated after the
quantization, it is shown in [14] that a pilot sequence with equidistant phase
shifts in the range $[0,\frac{\pi}{2})$ is MSE-optimal with respect to the
conditional mean estimator (CME) for jointly Gaussian inputs. For the general
case of $B$-bit quantization, the optimization of the pilot sequence is
intractable due to the exponentially increasing number of quantization levels.
Therefore, in this work, we consider pilot sequences that have an equidistant
spacing in both the amplitude and the angle, i.e.,

$\displaystyle[\tilde{{\bm{a}}}]_{i}=\beta_{i}\exp\left({\operatorname{j}}\frac{\pi}{2P}(i-1)\right),~{}i\in\\{1,\dots,P\\},$
(2)

where $\beta_{i}=\frac{1}{2}+\frac{i-1}{2(P-1)}$ is the amplitude spacing. In
order to fulfill the power constraint $\|{\bm{a}}\|_{2}^{2}=P$, the pilot
vector is normalized as
${\bm{a}}=\frac{\sqrt{P}}{\|\tilde{{\bm{a}}}\|_{2}}\tilde{{\bm{a}}}$.

####  II-2 Quantizer Design

In this work, we consider scalar quantizers, i.e., the quantization is
performed elementwise on the input, where the real and imaginary parts are
quantized independently. The quantizer can be described by means of the
$2^{B}$ quantization labels $\ell_{i}$, $i\in\left\\{1,\dots,2^{B}\right\\}$,
and the quantization thresholds $\tau_{i}$,
$i\in\left\\{0,\dots,2^{B}\right\\}$, where $\tau_{0}=-\infty$ and
$\tau_{2^{B}}=\infty$ by definition. The quantization function of the
real/imaginary part of the signal can be denoted as

$\displaystyle Q_{B}(x)=\sum_{i=1}^{2^{B}}\ell_{i}\hskip
1.0pt\chi(\tau_{i-1}\leq x<\tau_{i}).$ (3)

For the case of one-bit quantization $B=1$, the quantization function of the
complex-valued signal ${\bm{y}}$ can be expressed as

$\displaystyle
Q_{1}({\bm{y}})=\frac{1}{\sqrt{2}}\left(\operatorname{sign}(\Re({\bm{y}}))+{\operatorname{j}}\operatorname{sign}(\Im({\bm{y}}))\right).$
(4)

Practicable ADCs usually have uniformly spaced quantization thresholds with a
constant step size $\Delta$, which depends on the input distribution, and the
quantization labels are placed in the middle of two quantization thresholds.
For the case of zero-mean Gaussian input with variance one, the optimal values
for $\Delta$ are computed numerically in [35].

Under the assumption that the elementwise quantizer input is zero-mean
Gaussian distributed with variance $1+\sigma^{2}$ following the considered SNR
definition, we choose the SNR-dependent step size as suggested in [36] as

$\displaystyle\Delta=\sqrt{\tfrac{1}{2}(1+\sigma^{2})}\Delta_{*}$ (5)

where $\Delta_{*}$ is the step size for the standard Gaussian input [35].
Although the quantizer input is generally not Gaussian distributed, this
choice gives a reasonable performance with regard to practical feasibility.
The necessary scaling in (5) is resolved by automatic gain control in
practice. We note that there exist more sophisticated quantizer designs, e.g.,
non-uniform scalar quantization [35, 37]; however, the uniform quantizer is
considered the most practicable choice in wireless communications [5].
Furthermore, the channel estimation techniques proposed in this work are not
limited to uniform quantization but can be straightforwardly extended to non-
uniform quantization.

####  II-3 Channel Models

We work with the 3rd Generation Partnership Project (3GPP) spatial channel
model [38, 39] where channels are modeled conditionally Gaussian:
${\bm{h}}|{\bm{\delta}}\sim\mathcal{N}_{\mathbb{C}}({\bm{0}},{\bm{C}}_{{\bm{\delta}}})$.
The random vector ${\bm{\delta}}$ collects the angles of arrival/departure and
path gains of the main propagation clusters between a MT and the BS. The main
angles are drawn independently and uniformly from the interval $[0,2\pi]$; the
path gains are also drawn uniformly and are subsequently normalized such that
they sum up to one. The BS employs a uniform linear array (ULA) such that the
spatial channel covariance matrix is given by

${\bm{C}}_{{\bm{h}}|{\bm{\delta}}}=\int_{-\pi}^{\pi}\omega(\phi;{\bm{\delta}}){\bm{t}}(\phi){\bm{t}}(\phi)^{\operatorname{H}}{\operatorname{d}}\phi.$
(6)

Here,
${\bm{t}}(\phi)=[1,{\operatorname{e}}^{{\operatorname{j}}\pi\sin(\phi)},\dots,{\operatorname{e}}^{{\operatorname{j}}\pi(N-1)\sin(\phi)}]^{\operatorname{T}}$
is the array steering vector for an angle of arrival $\phi$, and $\omega$ is a
power density consisting of a sum of weighted Laplace densities whose standard
deviations describe the angle spread of the propagation clusters [38]. For
every channel sample, we generate random angles and path gains, combined in
${\bm{\delta}}$, and then draw the sample as
${\bm{h}}\sim\mathcal{N}_{\mathbb{C}}({\bm{0}},{\bm{C}}_{{\bm{h}}|{\bm{\delta}}})$,
which results in an overall non-Gaussian channel distribution [29]. Note that
the conditional Gaussianity of the channel model is not connected to the
conditional Gaussianity of the proposed latent models since the inference of
(6) from a single snapshot is intractable.

To ensure a broader evaluation of different channel models, version 2.4 of the
QuaDRiGa channel simulator [40, 41] is used to generate channel samples. We
simulate an urban macrocell scenario at a center frequency of 6 GHz. The BS’
height is 25 meters, and it covers a $120^{\circ}$ sector. The distances
between the MTs and the BS are in the range of 35–500 meters. We either
consider a pure line-of-sight (LOS) scenario or a mixed LOS/non-line-of-sight
(NLOS) scenario, where in 80% of the cases, the MTs are located indoors at
different floor levels, whereas the MTs’ height is 1.5 meters in the case of
outdoor locations. The BS is equipped with a ULA with $N$ “3GPP-3D” antennas,
and the MTs employ an omnidirectional antenna. The generated channels are
post-processed to remove the effective path gain [41, Sec. 2.7].

####  II-4 Training Datasets

In this work, we consider the channel distribution to be unknown and
arbitrarily complex by means of the sophisticated channel models, cf. Section
II-3. However, we assume the availability of a representative dataset which is
comprised of samples stemming from the respective channel model. In practice,
this means that data samples from the respective BS cell are available. In
this work, we discuss two different setups. First, we assume the availability
of a training dataset consisting of $T$ ground-truth channel samples
$\mathcal{H}=\\{{\bm{h}}_{t}\\}_{t=1}^{T}$. This can be achieved in practice
via measurement campaigns or digital twins, e.g., ray tracing. Afterward, we
consider a training dataset consisting solely of noisy and quantized pilot
observations $\mathcal{R}=\\{{\bm{r}}_{t}\\}_{t=1}^{T}$. Therefore, pilot
observations from the regular BS operation can be utilized, and no ground-
truth channel samples are needed.

##  III Bussgang Estimator

In this section, we briefly revise the linear MMSE estimator based on the
Bussgang decomposition, which is a direct consequence of Bussgang’s theorem
[21]. The Bussgang theorem, and hence also the decomposition, is valid only if
the input to the quantizer is zero-mean Gaussian distributed. However, even if
this is generally not true in the considered setting, the linear MMSE
estimator based thereon is a reasonable baseline for channel estimation in
quantized systems.

In particular, under the assumption of a jointly Gaussian quantizer input, the
Bussgang decomposition implies that the system in (1) can be written as a
linear combination of the desired signal part and an uncorrelated distortion
${\bm{q}}$ as

${\bm{r}}=Q_{B}({\bm{y}})={\bm{B}}{\bm{y}}+{\bm{\eta}}={\bm{B}}{\bm{A}}{\bm{h}}+{\bm{q}},$
(7)

where ${\bm{B}}$ is the Bussgang gain that can be obtained from the linear
MMSE estimation of ${\bm{r}}$ from ${\bm{y}}$ as
${\bm{B}}={\bm{C}}_{{\bm{r}}{\bm{y}}}{\bm{C}}_{{\bm{y}}}^{-1}$, cf. [42, Sec.
9.2], and where the distortion term ${\bm{q}}={\bm{B}}{\bm{n}}+{\bm{\eta}}$
contains both the AWGN ${\bm{n}}$ and the quantization noise ${\bm{\eta}}$.
The Bussgang gain matrix for a uniform quantizer is derived in [43] and is
computed as

$\displaystyle{\bm{B}}=\frac{\Delta}{\sqrt{\pi}}{\bm{D}}_{{\bm{y}}}^{-\frac{1}{2}}\sum_{i=1}^{2^{B}-1}\exp\left(-\Delta^{2}\left(i-2^{B-1}\right)^{2}{\bm{D}}_{{\bm{y}}}^{-1}\right)$
(8)

where ${\bm{D}}_{{\bm{y}}}=\operatorname{diag}({\bm{C}}_{{\bm{y}}})$. In the
case of one-bit quantization $B=1$, by choosing $\Delta=\sqrt{2}$, we get the
well-known solution

$\displaystyle{\bm{B}}=\sqrt{\frac{2}{\pi}}\operatorname{diag}({\bm{C}}_{{\bm{y}}})^{-\frac{1}{2}}.$
(9)

As the statistically equivalent model (7) is linear, one can formulate the
linear MMSE estimator

$\hat{{\bm{h}}}_{\text{Buss}}={\bm{C}}_{{\bm{h}}{\bm{r}}}{\bm{C}}_{{\bm{r}}}^{-1}{\bm{r}}.$
(10)

The cross-correlation matrix between the channel and the received signal is
calculated as
${\bm{C}}_{{\bm{h}}{\bm{r}}}={\operatorname{E}}[{\bm{h}}({\bm{B}}{\bm{A}}{\bm{h}}+{\bm{q}})^{\operatorname{H}}]={\bm{C}}_{{\bm{h}}}{\bm{A}}^{\operatorname{H}}{\bm{B}}^{\operatorname{H}}$
which follows from the fact that the noise term ${\bm{q}}$ is uncorrelated
with the channel ${\bm{h}}$, see [22, Appendix A]. For the one-bit
quantization case, the auto-correlation matrix is equal to the covariance
matrix ${\bm{C}}_{{\bm{r}}}$ due to the elimination of the amplitude
information and can be calculated in closed-form via the so-called arcsine law
[44] as

$\displaystyle{\bm{C}}_{{\bm{r}}}=\frac{2}{\pi}\left(\arcsin\left({\bm{D}}_{{\bm{y}}}^{-\frac{1}{2}}\Re({\bm{C}}_{{\bm{y}}}){\bm{D}}_{{\bm{y}}}^{-\frac{1}{2}}\right)\right.$
(11)
$\displaystyle\left.+{\operatorname{j}}\arcsin\left({\bm{D}}_{{\bm{y}}}^{-\frac{1}{2}}\Im({\bm{C}}_{{\bm{y}}}){\bm{D}}_{{\bm{y}}}^{-\frac{1}{2}}\right)\right).$

Unfortunately, for the multi-bit quantization case, no closed-form expression
for ${\bm{C}}_{{\bm{r}}}$ exists. Besides that, the computation of the
variances after the quantization are no longer non-trivial. As shown in [5,
eq. (2.14)] (adapted for the complex-valued case), for Gaussian input and the
uniform quantizer, the variances are computed as

$[{\bm{C}}_{{\bm{r}}}]_{i,i}=\sum_{i=1}^{2^{B}}2\ell_{i}^{2}\left({\operatorname{\Phi}}\left(\sqrt{2}\tau_{i}c_{i}\right)-{\operatorname{\Phi}}\left(\sqrt{2}\tau_{i-1}c_{i}\right)\right),$
(12)

where $c_{i}=[{\bm{C}}_{{\bm{y}}}]_{i,i}^{-\frac{1}{2}}$. Although there exist
various practicably feasible approximations for the evaluation of the involved
Gaussian CDF, cf. [45, 46], the evaluation of (12) may still be problematic in
time-critical systems; thus, a reasonable approximation is used in this work.
By assuming that the signal’s variance does not change for different antennas,
the Bussgang gain becomes a scaled identity of the form
${\bm{B}}=\rho\operatorname{\mathbf{I}}$, cf. (8). By further neglecting
cross-correlations of the quantization distortion, the quantized covariance
matrix ${\bm{C}}_{{\bm{r}}}$ is well approximated, especially in the low SNR
regime, by, cf. [3],

$\displaystyle{\bm{C}}_{{\bm{r}}}\approx\rho^{2}{\bm{C}}_{{\bm{y}}}+(1-\rho^{2})\operatorname{diag}({\bm{C}}_{{\bm{y}}}).$
(13)

Importantly, the resulting covariance matrix ${\bm{C}}_{{\bm{r}}}$ remains
positive semi-definite (PSD) if $0\leq\rho^{2}\leq 1$. We note that the
expression for ${\bm{C}}_{{\bm{r}}}$ with respect to ${\bm{C}}_{{\bm{y}}}$
generally depends on the quantizer choice and the input distribution, and
useful approximations can be found differently, cf., e.g., [4]. However, the
design of the channel estimation algorithms in this work is not founded upon
the choice in (13), and different approximations can be utilized.

##  IV Parameterized Bussgang Channel Estimators

In this section, we introduce several novel channel estimators which can be
applied in coarse quantization systems for variable numbers of quantization
bits. Therefore, we utilize generative models to approximate the underlying
unknown channel distribution by means of a given dataset. The common
characteristic of the generative models that are used in this context is their
modeling of the data distribution as conditionally Gaussian through a latent
variable model. These models include the GMM, the MFA, and the VAE. Recently,
all three models have been successfully deployed for channel estimation in
high-resolution scenarios, i.e., under the assumption of perfect ADCs [29, 30,
31, 32, 33, 34]. In this work, we exploit the conditional Gaussianity of the
mentioned generative models to parameterize a Bussgang estimator. This is
motivated by the insight that, conditioned on the latent variable, the
underlying channel is Gaussian, which establishes the prerequisite for the
application of Bussgang’s theorem. In this section, it is assumed that a
training dataset $\mathcal{H}$ of ground-truth channel samples is available,
cf. Section II-4.

###  IV-A GMM-based Bussgang Estimator

We start by deriving the GMM-based estimator, which parameterizes a
componentwise Bussgang estimator. Generally, a GMM is a probability density
function (PDF) of the form

$\displaystyle
f_{{\bm{h}}}^{(K)}({\bm{h}})=\sum_{k=1}^{K}\pi_{k}\mathcal{N}_{\mathbb{C}}({\bm{h}};{\bm{\mu}}_{{\bm{h}}|k},{\bm{C}}_{{\bm{h}}|k})$
(14)

where $K$ is the number of mixture components, and
$\\{\pi_{k},{\bm{\mu}}_{{\bm{h}}|k},{\bm{C}}_{{\bm{h}}|k}\\}_{k=1}^{K}$ is the
set of parameters of the GMM, namely the mixing coefficients, the means, and
the covariances of the Gaussian components. The parameters of the GMM are
fitted via the EM algorithm for a given training dataset $\mathcal{H}$ of
channel samples, cf. [47, Ch. 9]. An essential property of GMMs is that for a
given data sample, the responsibility of each component can be computed as,
cf. [47, Ch. 9],
$p(k|{\bm{h}})\propto\pi_{k}\mathcal{N}_{\mathbb{C}}({\bm{h}};{\bm{\mu}}_{{\bm{h}}|k},{\bm{C}}_{{\bm{h}}|k})$.

To ensure the validity of the Bussgang theorem for each component, we enforce
the component means to be zero, i.e., ${\bm{\mu}}_{{\bm{h}}|k}={\bm{0}}$,
$\forall k\in\\{1,\dots,K\\}$. To reflect this constraint in the fitting
process, the component means are set to zero in every M-step of the EM
algorithm. Naturally, the zero-mean constraint diminishes the capability of
the GMM to a certain extent concerning its ability to approximate the true
underlying distribution. Nevertheless, since a feasible wireless channel
distribution is considered to be zero-mean with a decreasing probability
density towards higher amplitudes, cf., e.g., [38], the loss of accuracy of
the model can be considered to be small. Moreover, restricting the component
means prevents overfitting and allows to model high-dimensional data [48].

After the GMM is trained, it is used for channel estimation similar to [29]
but with multiple adaptions in order to take the quantization effect into
account. We first note that if the channel distribution is modeled as a zero-
mean GMM, also the distribution of the unquantized receive signal ${\bm{y}}$
follows a zero-mean GMM with covariances
${\bm{C}}_{{\bm{y}}|k}={\bm{A}}{\bm{C}}_{k}{\bm{A}}^{\operatorname{H}}+\sigma^{2}\operatorname{\mathbf{I}}$,
$\forall k\in\\{1,\dots,K\\}$, cf. (1). Consequently, for each component, we
can apply the Bussgang decomposition to find a statistically equivalent model:

$\displaystyle{\bm{r}}|k=({\bm{B}}_{k}{\bm{A}}{\bm{h}}+{\bm{q}}_{k})|k,$ (15)

where ${\bm{B}}_{k}$ is the Bussgang gain of component $k$, and
${\bm{q}}_{k}={\bm{B}}_{k}{\bm{n}}+{\bm{\eta}}$, cf. (7). The componentwise
Bussgang gain is computed via (8) or (9) by plugging in the $k$th covariance
${\bm{C}}_{{\bm{y}}|k}$ for the multi- or one-bit case, respectively. The
evaluation of the discrete distribution $p({\bm{r}}|k)$ is intractable in
general since the cardinality of its discrete support increases exponentially
in the number of dimensions [14]. Thus, in order to evaluate the
responsibility for a given pilot signal, we assume that the quantized receive
signal follows a zero-mean GMM distribution with the same second-order
moments; this assumption effectively resembles approximate inference [47, Ch.
10]. The covariance matrix of component $k$, named ${\bm{C}}_{{\bm{r}}|k}$, is
thereby computed via (11) or (13) for the one- or multi-bit quantization case,
respectively, by plugging in the component’s unquantized covariance
${\bm{C}}_{{\bm{y}}|k}$. Thereby, in the case of a covariance matrix
${\bm{C}}_{{\bm{h}}|k}$ with a non-constant diagonal, the scaling parameter
$\rho_{k}$ in (13) for the $k$th component is approximated via
$\rho_{k}=\min(\frac{1}{N}\sum_{i=1}^{N}[{\bm{B}}_{k}]_{i,i},1)$, which
ensures that the resulting matrix is PSD. This results in the following
responsibility evaluation of the quantized receive signal:

$\displaystyle
p(k|{\bm{r}})\approx\frac{\pi_{k}\mathcal{N}_{\mathbb{C}}({\bm{r}};{\bm{0}},{\bm{C}}_{{\bm{r}}|k})}{\sum_{i=1}^{K}\pi_{i}\mathcal{N}_{\mathbb{C}}({\bm{r}};{\bm{0}},{\bm{C}}_{{\bm{r}}|i})}.$
(16)

The final channel estimate is computed via the convex combination of the
componentwise Bussgang estimators, similar to the high-resolution case [29],
parameterized by the GMM covariances, which yields

$\displaystyle\hat{{\bm{h}}}_{\text{BGMM}}^{(K)}({\bm{r}})=\sum_{k=1}^{K}p(k|{\bm{r}}){\bm{C}}_{{\bm{h}}|k}{\bm{A}}^{\operatorname{H}}{\bm{B}}_{k}^{\operatorname{H}}{\bm{C}}_{{\bm{r}}|k}^{-1}{\bm{r}}.$
(17)

As shown in [31], it is possible to enforce different structural constraints
for the GMM’s covariances, such as a circulant (“GMM circ”) or Toeplitz (“GMM
toep”) structure. The covariance matrix of the $k$th GMM component is thereby
constrained to be of the form
${\bm{C}}_{{\bm{h}}|k}={\bm{Q}}^{\operatorname{H}}\operatorname{diag}({\bm{c}}_{{\bm{h}}|k}){\bm{Q}}$
where ${\bm{Q}}$ is an (oversampled) discrete Fourier transform (DFT) matrix
and $[{\bm{c}}_{{\bm{h}}|k}]_{i}\in\mathbb{R}_{+}$. These structural
constraints reflect the typical array geometries of a BS, e.g., a ULA, and
result in a reduced number of parameters and a lower online complexity of the
estimator due to the usage of fast Fourier transforms (FFTs) [31]. The
structural constraints are, without limitation, also applicable in coarsely
quantized systems. The necessary memory overhead and computational complexity
of the resulting estimators are discussed in more detail in Section VII.

###  IV-B MFA-based Bussgang Estimator

A related concept to the GMM is the MFA model, which, in addition to a
discrete latent variable $k$ which describes the mixture component, also
contains a continuous latent variable ${\bm{z}}\in\mathbb{C}^{L}$ of lower
dimension, i.e., $L<N$ holds [49, Ch. 12], [32]. This effectively models the
data on a piecewise linear subspace. After integrating out the continuous
latent variable ${\bm{z}}\sim\mathcal{N}({\bm{0}},\operatorname{\mathbf{I}})$,
the PDF of the MFA model is a special form of a GMM with low-rank plus
diagonal-constrained covariances of the form

$\displaystyle
f^{(K,L)}({\bm{h}})=\sum_{k=1}^{K}\pi_{k}\mathcal{N}_{\mathbb{C}}({\bm{h}};{\bm{\mu}}_{{\bm{h}}|k},{\bm{W}}_{{\bm{h}}|k}{\bm{W}}_{{\bm{h}}|k}^{\operatorname{H}}+{\bm{\Psi}}_{{\bm{h}}|k}),$
(18)

where ${\bm{W}}_{{\bm{h}}|k}\in\mathbb{C}^{N\times L}$ is the factor loading
matrix and ${\bm{\Psi}}_{{\bm{h}}|k}\in\mathbb{C}^{N\times N}$ is a diagonal
matrix. In order to fit the parameters
$\\{\pi_{k},{\bm{\mu}}_{{\bm{h}}|k},{\bm{W}}_{{\bm{h}}|k},{\bm{\Psi}}_{{\bm{h}}|k}\\}_{k=1}^{K}$
of the MFA model for a given dataset $\mathcal{H}$ of channel realizations, an
EM algorithm can be used [49, Ch. 12]. After training, by defining
${\bm{C}}_{{\bm{h}}|k}={\bm{W}}_{{\bm{h}}|k}{\bm{W}}_{{\bm{h}}|k}^{\operatorname{H}}+{\bm{\Psi}}_{{\bm{h}}|k}$,
the model can be effectively treated as a GMM. A zero-mean MFA model with
${\bm{\mu}}_{{\bm{h}}|k}={\bm{0}}$ $\forall k\in\\{1,\dots,K\\}$ can be
similarly enforced as in the GMM case. Similar to [32], we set
${\bm{\Psi}}_{{\bm{h}}|k}=\psi_{{\bm{h}}|k}\operatorname{\mathbf{I}}$,
$\forall k\in\\{1,\dots,K\\}$.

The main advantage of the MFA model in the context of channel estimation, in
contrast to the GMM, lies in the reduced number of parameters due to the low-
dimensional latent space, which mitigates overfitting effects during training.
Thus, it is a more robust model for lower numbers of training data, as
demonstrated in [32]. Since the MFA also parameterizes a conditionally
Gaussian distribution, the componentwise Bussgang estimator is equivalently
applicable as for the GMM case, see Section IV-A, which yields the MFA-
parameterized Bussgang channel estimator $\hat{{\bm{h}}}_{\text{BMFA}}$ by
substituting the respective covariances in (17).

###  IV-C VAE-based Bussgang Estimator

Figure 1: Proposed adapted VAE architecture with the encoder, latent space,
and decoder together with the parameterized distributions and the
reparameterization trick.

The VAE was introduced in [50] and has attracted a lot of interest in the area
of generative modeling due to its strong performance, which builds on the
basis of NNs that are used for the encoder and decoder of the VAE, cf. Fig. 1.
In [33, 34], the VAE was successfully utilized for channel estimation in high-
resolution systems. In contrast to the GMM and MFA, the VAE comprises a
continuous and nonlinear latent space, encoded by the low-dimensional latent
vector ${\bm{z}}\in\mathbb{R}^{L}$. The most common design choice is a
Gaussian model for the latent vector, i.e.,
${\bm{z}}\sim\mathcal{N}({\bm{0}},\operatorname{\mathbf{I}})$. Since the
variational inference task is no longer tractable by the classical EM
algorithm, NNs in combination with the reparameterization trick are used to
train the VAE [51, 50]. In this work, we choose the parameterized encoder and
decoder distributions as

$\displaystyle q_{\phi}({\bm{z}}|{\bm{r}})$
$\displaystyle=\mathcal{N}({\bm{z}};{\bm{\mu}}_{\phi}({\bm{r}}),\operatorname{diag}({\bm{\sigma}}_{\phi}^{2}({\bm{r}}))),$
(19) $\displaystyle p_{\theta}({\bm{h}}|{\bm{z}})$
$\displaystyle=\mathcal{N}_{\mathbb{C}}({\bm{h}};{\bm{0}},\operatorname{\mathbf{F}}^{\operatorname{H}}\operatorname{diag}({\bm{c}}_{\theta}({\bm{z}}))\operatorname{\mathbf{F}}),~{}{\bm{c}}_{\theta}\in\mathbb{R}_{+}^{N},$
(20)

respectively. The matrix $\operatorname{\mathbf{F}}$ is a DFT matrix such that
the parameterized channel covariance matrix of the VAE is a circulant matrix,
cf. [34]. This choice results in a reduced number of parameters and is
justified by the array geometry of the BS, similar to the GMM case, cf.
Section IV-A.

For every training data point ${\bm{h}}_{t}\in\mathcal{H}$, the VAE computes
the ELBO on the log-likelihood [50]

$\displaystyle\log p_{\theta}({\bm{h}}_{t})\geq{\operatorname{E}}_{\phi}[\log
p_{\theta}({\bm{h}}_{t}|{\bm{z}})]-{\operatorname{D}}_{\text{KL}}(q_{\phi}({\bm{z}}|{\bm{h}}_{t})\hskip
1.0pt||\hskip 1.0ptp({\bm{z}})).$

By plugging in (19) and (20) and ignoring the constant terms (they do not
influence the optimization), the ELBO is utilized as the loss function
$L_{\phi,\theta}$ of the VAE, which reads as, cf. [34],

$\displaystyle
L_{\phi,\theta}({\bm{h}}_{t})=\sum_{n=1}^{N}-\log[{\bm{c}}_{\theta}]_{n}-{\bm{h}}^{\operatorname{H}}\operatorname{\mathbf{F}}^{\operatorname{H}}\operatorname{diag}({\bm{c}}_{\theta})^{-1}\operatorname{\mathbf{F}}{\bm{h}}$
(21)
$\displaystyle+\sum_{\ell=1}^{L}\log[{\bm{\sigma}}_{\phi}]_{\ell}-\frac{1}{2}[{\bm{\mu}}_{\phi}]_{\ell}^{2}-\frac{1}{2}[{\bm{\sigma}}_{\phi}]_{\ell}^{2}.$

As demonstrated in [33, 34], the VAE can be utilized to parameterize a channel
estimator by means of the conditional Gaussianity at the output of the decoder
of the VAE in combination with linear MMSE filters. After the training, we
utilize the VAE to parameterize a channel covariance matrix for each quantized
receive pilot ${\bm{r}}$ by forwarding the pilot through the encoder and then
using the latent mean ${\bm{\mu}}_{\phi}({\bm{r}})$ as input to the decoder
which yields the channel covariance matrix
${\bm{C}}_{\theta}=\operatorname{\mathbf{F}}^{\operatorname{H}}\operatorname{diag}({\bm{c}}_{\theta})\operatorname{\mathbf{F}}$.
The VAE-parameterized Bussgang estimator then reads as

$\hat{{\bm{h}}}_{\text{BVAE}}({\bm{r}})=\operatorname{\mathbf{F}}^{\operatorname{H}}\operatorname{diag}({\bm{c}}_{\theta})\operatorname{\mathbf{F}}{\bm{A}}^{\operatorname{H}}{\bm{B}}_{{\bm{z}}}^{\operatorname{H}}{\bm{C}}_{{\bm{r}}|{\bm{z}}}^{-1}{\bm{r}},$
(22)

where ${\bm{B}}_{{\bm{z}}}$ is computed by plugging
${\bm{C}}_{{\bm{y}}|{\bm{z}}}={\bm{A}}{\bm{C}}_{\theta}{\bm{A}}^{\operatorname{H}}+\sigma^{2}\operatorname{\mathbf{I}}$
into (8) or (9). Similarly, the quantized receive covariance
${\bm{C}}_{{\bm{r}}|{\bm{z}}}$ is computed by plugging
${\bm{C}}_{{\bm{y}}|{\bm{z}}}$ into (11) or (13).

For the encoder and decoder, we use a four-layer feed-forward NN with
rectified linear unit (ReLU) activation functions, respectively, for which we
stack the real and imaginary parts of the pilot observation at the input to
the encoder. In the case of multiple pilot observations, we add a
convolutional NN (CNN) with $P/2$ layers and ReLU activation functions before
the encoder, which perform $1\times 1$ convolutions in order to always have a
$2N$-dimensional input at the encoder. This modification drastically reduces
the number of parameters and simplifies a forward pass through the VAE. The
complete VAE architecture is detailed in Fig. 1.

##  V Learning from Quantized Data

As discussed above, the availability of a large training dataset $\mathcal{H}$
consisting of representative ground-truth channel samples for a whole BS cell
(radio propagation scenario) is questionable in practical communication
systems. Although different approaches exist, such as ray tracing, to generate
a training dataset that mimics the underlying channel distribution, it is
unclear whether they can sufficiently capture the characteristics of a real
communication scenario. A different idea is to train the respective models
directly on pilot observations and mitigate imperfections, e.g., additive
noise or sparsely allocated pilot positions, through model-based training
adaptations. This has already been shown to work well for high-resolution
scenarios [34, 52]. However, learning a generative model from quantized data
poses a significant challenge due to the pronounced nonlinear distortion
resulting from low-resolution quantization.

###  V-A Covariance Recovery for GMM Approximation

Recovering the unquantized covariance matrix of the input to a one-bit
quantizer solely from quantized data has gained a lot of interest very
recently [53, 54, 55]. Since the amplitude information is lost in the case of
a zero-threshold one-bit quantization, only the normalized correlation matrix
can be obtained via the inverse arcsine law [44]. To resolve this issue, non-
zero or time-varying quantizer thresholds were considered in order to be able
to estimate the variances of the input signal and, thus, the whole covariance
matrix [53, 54, 55]. The work in [56] validates that the covariance recovery
technique from [53] can be used in combination with the Bussgang estimator in
order to perform channel estimation. However, the considered quantizer designs
with non-zero thresholds are challenging to implement in communication systems
and may require more sophisticated analog and digital signal processing,
ultimately resulting in performance losses. In contrast, considering multi-bit
quantization of the input signal, coarse amplitude information is preserved
because of the multi-level quantization, even with a fixed zero-threshold.
However, up to now, there is no covariance recovery algorithm proposed for
this case.

We derive a novel low-complexity covariance recovery algorithm for the multi-
bit case by splitting the task into estimating the correlation matrix and the
variances independently. Let us define the following problem statement where
we reuse the notation from above for simplicity. Consider a dataset
$\mathcal{R}=\\{{\bm{r}}_{t}\\}_{t=1}^{T}$ of $T$ samples of the form
${\bm{r}}_{t}=Q_{B}({\bm{y}}_{t})$ where
${\bm{y}}_{t}\sim\mathcal{N}_{\mathbb{C}}({\bm{0}},{\bm{C}}_{{\bm{y}}})$ and
$B>1$. The task is to recover ${\bm{C}}_{{\bm{y}}}$ from the $T$ quantized
samples. Since no closed-form solution for the correlation matrix
${\bm{R}}_{{\bm{y}}}=\operatorname{diag}({\bm{C}}_{{\bm{y}}})^{-\frac{1}{2}}{\bm{C}}_{{\bm{y}}}\operatorname{diag}({\bm{C}}_{{\bm{y}}})^{-\frac{1}{2}}$
in the case of multi-bit quantization exists, we simply discard the samples’
amplitude information, effectively treating them as one-bit quantization data.
Because of that, the closed-form expression for estimating the unquantized
correlation matrix ${\bm{R}}_{{\bm{y}}}$ by means of the one-bit sample
covariance matrix
$\hat{{\bm{C}}}_{\text{1bit}}=\frac{1}{T}\sum_{t=1}^{T}Q_{1}({\bm{r}}_{t})Q_{1}({\bm{r}}_{t})^{\operatorname{H}}$
can be obtained via the inverse arcsine law:

$\displaystyle\hat{{\bm{R}}}_{{\bm{y}}}=\sin\left(\frac{\pi}{2}\Re(\hat{{\bm{C}}}_{\text{1bit}})\right)+{\operatorname{j}}\sin\left(\frac{\pi}{2}\Im(\hat{{\bm{C}}}_{\text{1bit}})\right).$
(23)

Unfortunately, although having a closed-form solution, the resulting
correlation estimate is not necessarily PSD [53, 57]. However, this can be
resolved by a projection onto the set of PSD matrices, as discussed later.

Since the quantization acts elementwise on the real and imaginary part
independently, it is sufficient to derive the variance estimation for a real-
valued scalar $y\sim\mathcal{N}(0,\xi^{2})$ and $r=Q_{B}(y)\in\mathcal{R}$ for
ease of notation. We note that the amplitude of $y$ follows the half-normal
distribution. The corresponding CDF of the half-normal distribution is given
by ${\operatorname{P}}(|y|\leq\tau)=\operatorname{erf}(\tau/\sqrt{2\xi^{2}})$.
Because the CDF is fully parameterized by the input signal’s variance
$\xi^{2}$, one can utilize the coarse amplitude information after the
quantizer for its estimation. By defining the positive quantization thresholds
as $\tilde{\tau}_{i}<\infty$, $i\in\\{1,\dots,2^{B-1}-1\\}$, i.e.,
$\tilde{\tau}_{i}=\tau_{i+2^{B-1}}$, one can estimate the probability of
observing a sample with an amplitude of at most $\tilde{\tau}_{i}$ by
$\hat{{\operatorname{P}}}(|y|\leq\tilde{\tau}_{i})=\frac{1}{T}\sum_{i=1}^{T}\chi(|r_{t}|\leq\tilde{\tau}_{i})$.

For circularly symmetric Gaussian distributed complex-valued input and
multiple quantization thresholds, i.e., $B>2$, an overdetermined system of
equations can be constructed using the different quantizer thresholds for both
the real- and imaginary parts, yielding $2^{B}-2$ equations. The subtraction
of two comes from the fact that the last quantization regions up to infinity
are uninformative since, in this case, the (sample) probability is always one.
In summary, the equation system accounting for the real part of the input is
of the form

$\displaystyle\operatorname{erf}\left(\frac{\tilde{\tau}_{i}}{\sqrt{2\xi^{2}}}\right)=\frac{1}{T}\sum_{t=1}^{T}\chi\left(|\Re(r_{t})|\leq\tilde{\tau}_{i}\right)$
(24)

with $i\in\\{1,\dots,2^{B-1}-1\\}$. The remaining half of the equation system
is built similarly by replacing the real with the imaginary part. Note that
the equation system only depends on the unknown variance parameter $\xi^{2}$.
Geometrically, we aim to interpolate the sample probabilities belonging to the
different thresholds by a Gaussian CDF curve in a LS sense with the adjustable
variance parameter $\xi^{2}$. Since the derivative of the CDF is trivially
given by the Gaussian PDF, a simple Gauss-Newton approach can be utilized for
solving the nonlinear LS problem. For a robust initial starting point
$\xi_{0}^{2}$, a solution to the equation with the quantizer’s largest
$\tilde{\tau}_{i}$ is used.

In the multi-dimensional case, if the input signal’s variance is assumed to be
different for each antenna, the nonlinear LS problem can be solved for each
dimension independently, yielding an estimate of
$\operatorname{diag}({{\bm{C}}}_{{\bm{y}}})$. Otherwise, the equation systems
for each dimension can be combined to yield a more accurate estimate of the
single variance parameter. Note that in the complex-valued case, the estimated
variance has to be scaled by a factor of two to account for the sum of the
real and imaginary parts. Finally, the full covariance matrix estimate is
computed as

$\hat{{\bm{C}}}_{{\bm{y}}}=\operatorname{diag}(\hat{{\bm{C}}}_{{\bm{y}}})^{\frac{1}{2}}\hat{{\bm{R}}}_{{\bm{y}}}\operatorname{diag}(\hat{{\bm{C}}}_{{\bm{y}}})^{\frac{1}{2}}.$
(25)

The derived covariance recovery scheme is now used in order to fit the GMM’s
covariances by only using quantized data. For simplicity, we assume that the
training data stems from single snapshot observations, i.e.,
${\bm{r}}=Q_{B}({\bm{h}}+{\bm{n}})$ with $P=1$, which can always be enforced
by pre-processing. In each iteration of the EM algorithm, the M-step is
adapted by using the proposed covariance recovery algorithm for estimating the
unquantized covariance matrix due to the Gaussianity of each GMM component.
The necessary change to the purely Gaussian setting from before is that the
responsibilities, computed in each E-step, are used in order to weight the
sample probability for component $k$ accordingly as

$\hat{{\operatorname{P}}}(|\Re([{\bm{y}}]_{n})|\leq\tilde{\tau}_{i}|k)=\frac{1}{N_{k}}\sum_{t=1}^{T}p(k|{\bm{r}}_{t})\chi(|\Re([{\bm{r}}_{t}]_{n})|\leq\tilde{\tau}_{i}),$
(26)

$\forall i\in\\{1,\dots,2^{B-1}-1\\}$, where
$N_{k}=\sum_{t=1}^{T}p(k|{\bm{r}}_{t})$. Note that $\Re([{\bm{r}}_{t}]_{n})$
is replaced by $\Im([{\bm{r}}_{t}]_{n})$ for the second half of the equation
system. Since the quantizer’s input signal is also distorted with AWGN, the
M-step adaptation from [52, Th. 1] for noisy data is used in addition by means
of subtracting the noise covariance and afterward projecting to the set of PSD
matrices by performing an eigenvalue decomposition (EVD) and truncating the
negative eigenvalues. This also accounts for the possibly non-PSD correlation
estimate from (23). After estimating the channel covariance
$\hat{{\bm{C}}}_{{\bm{h}}|k}$ of component $k$ in this way, one first
determines $\hat{{\bm{C}}}_{{\bm{y}}|k}$ to eventually construct the
covariance of the quantized observation $\hat{{\bm{C}}}_{{\bm{r}}|k}$ by using
one of the approximations given in (12) or (13). Since the complexity is not
crucial in the offline learning, we utilize the accurate formula for the
variance (12). The necessary adaptations in the M-step are concisely
summarized in Algorithm 1. The so-found covariance matrix
$\hat{{\bm{C}}}_{{\bm{r}}|k}$ is afterward used to compute the
responsibilities in the E-step, similar to (16).

By the law of large numbers, the variance estimate based on the sample
probability is a consistent estimator and, together with the closed-form
solution for the correlation matrix, yields a consistent covariance estimator.
This is a very convenient property since pilot observations are received in
large numbers during regular operation at the BS, and thus, a large training
dataset $\mathcal{R}$ can be constructed.

Algorithm 1 Adapted M-step for quantized training data.

0: $\mathcal{R}=\\{{\bm{r}}_{t}\\}_{t=1}^{T}$, $\sigma^{2}$,
$\\{p(k|{\bm{r}}_{t})\\}_{t=1}^{T}$

1: for $k=1$ to $K$ do

2: $N_{k}=\sum_{t=1}^{T}p(k|{\bm{r}}_{t})$

3:
$\hat{{\bm{C}}}_{{\bm{r}}|k}^{\text{1bit}}=\frac{1}{N_{k}}\sum_{t=1}^{T}p(k|{\bm{r}}_{t})Q_{1}({\bm{r}}_{t})Q_{1}({\bm{r}}_{t})^{\operatorname{H}}$

4:
$\hat{{\bm{R}}}_{{\bm{y}}|k}=\sin\left(\frac{\pi}{2}\Re(\hat{{\bm{C}}}_{{\bm{r}}|k}^{\text{1bit}})\right)+{\operatorname{j}}\sin\left(\frac{\pi}{2}\Im(\hat{{\bm{C}}}_{{\bm{r}}|k}^{\text{1bit}})\right)$

5: for $n=1$ to $N$ do

6: Construct equation system for
$\operatorname{diag}(\hat{{\bm{C}}}_{{\bm{y}}|k})$ via (26).

7: Solve equation system via Gauss-Newton.

8: end

9:
$\hat{{\bm{C}}}_{{\bm{y}}|k}=\operatorname{diag}(\hat{{\bm{C}}}_{{\bm{y}}|k})^{\frac{1}{2}}\hat{{\bm{R}}}_{{\bm{y}}|k}\operatorname{diag}(\hat{{\bm{C}}}_{{\bm{y}}|k})^{\frac{1}{2}}$

10:
${\bm{V}}\operatorname{diag}(\hat{{\bm{c}}}_{{\bm{h}}|k}){\bm{V}}^{\operatorname{H}}=\operatorname{EVD}(\hat{{\bm{C}}}_{{\bm{y}}|k}-\sigma^{2}\operatorname{\mathbf{I}})$

11:
$\hat{{\bm{c}}}_{{\bm{h}}|k}^{\text{PSD}}=\max({\bm{0}},\hat{{\bm{c}}}_{{\bm{h}}|k})$
(elementwise max)

12:
$\hat{{\bm{C}}}_{{\bm{h}}|k}={\bm{V}}\operatorname{diag}(\hat{{\bm{c}}}_{{\bm{h}}|k}^{\text{PSD}}){\bm{V}}^{\operatorname{H}}$

13:
$\hat{{\bm{C}}}_{{\bm{y}}|k}=\hat{{\bm{C}}}_{{\bm{h}}|k}+\sigma^{2}\operatorname{\mathbf{I}}$

14: Compute $\operatorname{diag}(\hat{{\bm{C}}}_{{\bm{r}}|k})$ from
$\hat{{\bm{C}}}_{{\bm{y}}|k}$ via (12).

15: Compute ${\bm{B}}_{k}$ from $\hat{{\bm{C}}}_{{\bm{y}}|k}$ via (8).

16:
$\operatorname{nondiag}(\hat{{\bm{C}}}_{{\bm{r}}|k})=\operatorname{nondiag}({\bm{B}}_{k}\hat{{\bm{C}}}_{{\bm{y}}|k}{\bm{B}}_{k}^{\operatorname{H}})$

17: end

###  V-B Loss Function Adaptation for the VAE

Similar to the GMM, the VAE model can be trained directly from noisy pilot
observations by properly adapting the ELBO loss function [33, 34]. The idea,
thereby, is to modify the parameterized channel’s covariance matrix, which
stems from the output of the decoder such that it represents the covariance
matrix of the pilot observations in the loss function (21). Consequently, only
the channel covariance matrix is learned by the VAE through gradient updates.
For the case of coarsely quantized training samples from $\mathcal{R}$, the
expression of the quantized covariance matrix with respect to the channel
covariance matrix is not in closed-form and thus not differentiable, cf. (12).
However, one can use the approximation from (13) as shown in the following.

Since the training is done in the Fourier-transformed domain in order to
parameterize a circulant covariance matrix
${\bm{C}}_{\theta}=\operatorname{\mathbf{F}}^{\operatorname{H}}\operatorname{diag}({\bm{c}}_{\theta})\operatorname{\mathbf{F}}$,
we first simplify the expression for the diagonal term
$\operatorname{diag}({\bm{C}}_{\theta})=\operatorname{diag}(\operatorname{\mathbf{F}}^{\operatorname{H}}\operatorname{diag}({\bm{c}}_{\theta})\operatorname{\mathbf{F}})=\frac{1}{N}\sum_{n=1}^{N}[{\bm{c}}_{\theta}]_{n}\operatorname{\mathbf{I}}$,
which can be simply shown by writing
$\operatorname{diag}({\bm{C}}_{\theta})=\frac{1}{N}\operatorname{tr}({\bm{C}}_{\theta})\operatorname{\mathbf{I}}$
and utilizing the properties of the trace together with the unitary DFT matrix
$\operatorname{\mathbf{F}}$. Using this, we can now write
$\operatorname{diag}({\bm{C}}_{{\bm{y}},\theta})=\frac{1}{N}\sum_{n=1}^{N}([{\bm{c}}_{\theta}]_{n}+\sigma^{2})\operatorname{\mathbf{I}}$
when assuming single snapshot pilot observations as
${\bm{r}}=Q_{B}({\bm{h}}+{\bm{n}})$, such that
${\bm{C}}_{{\bm{r}},\theta}=\operatorname{\mathbf{F}}^{\operatorname{H}}\operatorname{diag}({\bm{c}}_{{\bm{r}},\theta})\operatorname{\mathbf{F}}$
where ${\bm{c}}_{{\bm{r}},\theta}$ equals

$(1-\rho_{\theta})^{2}({\bm{c}}_{\theta}+\sigma^{2}{\bm{1}})+\rho_{\theta}(1-\rho_{\theta})\frac{1}{N}\sum_{n=1}^{N}([{\bm{c}}_{\theta}]_{n}+\sigma^{2}){\bm{1}}$
(27)

and where $\rho_{\theta}$ is computed via (8) by plugging in
$\operatorname{diag}({\bm{C}}_{{\bm{y}},\theta})$. Thus, the VAE model can be
learned from quantized data by replacing ${\bm{c}}_{\theta}$ with
${\bm{c}}_{{\bm{r}},\theta}$ from (27) and ${\bm{h}}$ with ${\bm{r}}$ in the
loss function (21).

##  VI Baseline Channel Estimators

We compare the proposed parameterized generative modeling-aided channel
estimators with state-of-the-art baseline channel estimators for coarse
quantization systems. First, for the 3GPP channel model from Section II-3, we
have genie access to the true underlying channel covariance matrix
${\bm{C}}_{{\bm{h}}|{\bm{\delta}}}$ from (6) in the simulation. This allows us
to evaluate the genie-aided Bussgang estimator $\hat{{\bm{h}}}_{\text{Buss-
genie}}={\bm{C}}_{{\bm{h}}|{\bm{\delta}}}{\bm{A}}^{\operatorname{H}}{\bm{B}}_{{\bm{\delta}}}^{\operatorname{H}}{\bm{C}}_{{\bm{r}}|{\bm{\delta}}}$
where ${\bm{B}}_{{\bm{\delta}}}$ and ${\bm{C}}_{{\bm{r}}|{\bm{\delta}}}$ are
found by plugging ${\bm{C}}_{{\bm{h}}|{\bm{\delta}}}$ into (8) or (9) and (13)
or (11) for the multi-bit or one-bit case, respectively. Note that this
estimator is not feasible in practice but only serves as a lower bound on the
performance of the Bussgang estimator, which is the best linear estimator. The
corresponding curves are labeled as “Buss-genie”.

A practicably feasible approach that is primarily used in the literature is to
use the sample covariance matrix
$\hat{{\bm{C}}}_{{\bm{h}}}=\frac{1}{T}\sum_{t=1}^{T}{\bm{h}}_{t}{\bm{h}}_{t}^{\operatorname{H}}$
in combination with the Bussgang estimator (10), labeled as “Buss-Scov”. Note
that for this case, a training dataset of ground-truth channels $\mathcal{H}$
is necessary.

A simple baseline is the LS estimate based on the Bussgang decomposition (7),
i.e.,
$\hat{{\bm{h}}}_{\text{BLS}}={\bm{A}}^{\dagger}{\bm{B}}^{\dagger}{\bm{r}}$,
labeled as “BLS”, where ${\bm{A}}^{\dagger}$ is the pseudo-inverse of
${\bm{A}}$. For computing the Bussgang gain (8) or (9), we use the sample
covariance matrix $\hat{{\bm{C}}}_{{\bm{h}}}$ from above.

In [20], a CS-based channel estimator is proposed, which is a combination of
the EM algorithm for approximating the channel PDF in the sparse angular
domain and the GAMP algorithm to solve the sparse recovery problem. Note that
in this case, an EM algorithm is deployed online for each transmission link
and pilot observation, which is fundamentally different from the GMM approach,
which exploits the EM solely in the offline phase. We applied the EM-GM-GAMP
algorithm to estimate the channel parameters $\hat{{\bm{x}}}$ in the angular
domain, such that the final channel estimate is computed to
$\hat{{\bm{h}}}_{\text{EM-GM-GAMP}}=\operatorname{\mathbf{F}}\hat{{\bm{x}}}$,
labeled as “EM-GM-GAMP”.

We also evaluate a deep learning-based estimator, similar to [26], where a
three-layered feed-forward NN is trained to directly map the pilot observation
to a channel estimate. The ReLU function is used as the activation function in
all layers except the output layer. To achieve a fair comparison to the
proposed approaches, a single network is trained for the whole SNR range.
Similar to [26], the best performance was achieved by a drastic increase of
the neurons in the hidden layers. We, therefore, set the number of neurons in
both hidden layers to $2N^{2}$. The corresponding curves are labeled as “DNN”.

##  VII Memory and Complexity Analysis

TABLE I: Computational complexity and number of parameters of the discussed
estimators with example numbers for the case of $K=N=64$, $L=16$, and $P=4$.

Name |  Model Parameters (real-valued) |  Example (rounded) |  Online Complexity
---|---|---|---
GMM full |  $K\left(N^{2}+1\right)-1$ |  $2.62\cdot 10^{5}$ |  $\mathcal{O}(KPN^{2})$ (parallel in $K$)
GMM toep |  $K\left(4N+1\right)-1$ |  $1.64\cdot 10^{4}$ |  $\mathcal{O}(KPN^{2})$ (parallel in $K$)
GMM circ |  $K\left(N+1\right)-1$ |  $4.16\cdot 10^{3}$ |  $\mathcal{O}(KPN\log(PN))$ (parallel in $K$)
MFA |  $K\left(2LN+2\right)-1$ |  $1.31\cdot 10^{5}$ |  $\mathcal{O}(KPN^{2})$ (parallel in $K$)
VAE |  $6N^{2}+10NL+\frac{1}{2}P^{2}+\frac{1}{2}P+\mathcal{O}(N)$ |  $3.53\cdot 10^{4}$ |  $\mathcal{O}(N^{2}+PN\log(PN))$
DNN [26] |  $4N^{4}+(4P+4)N^{3}+4N^{2}+2N$ |  $7.24\cdot 10^{7}$ |  $\mathcal{O}(N^{3})$
BLMMSE |  $N^{2}$ |  $5.0\cdot 10^{3}$ |  $\mathcal{O}(PN^{2})$
EM-GM-GAMP [20] |  - |  - |  $\mathcal{O}(PN\log(PN))$
BLS |  - |  - |  $\mathcal{O}(PN)$

The offline memory requirements of data-based techniques and the algorithmic
online complexity are key features for channel estimation in real-time
systems. The number of parameters of the (structured) zero-mean GMM and the
MFA model is determined by the $K$ covariances and the number of mixing
coefficients. The corresponding linear MMSE filters for each component and SNR
value are fixed after the offline training, which means that they can be pre-
computed. Since this can be similarly done for the evaluation of the
responsibilities in (16), the overall online complexity is determined by
matrix-vector products for each component [29]. Notably, the computation of
the $K$ filters/responsibilities is trivially parallelizable, which is of
great importance in practical systems. For the case of circulant-structured
GMM covariances, cf. Section IV-A, the complexity reduces due to the usage of
FFTs [29].

For the VAE approach, the number of parameters and the complexity for a
forward pass through the network depend on the network architecture, cf.
Section IV-C. The resulting filter is computable by means of FFTs since a
circulant covariance matrix is parameterized, similar to the circulant GMM
case. We further note that the approaches that learn from quantized data, cf.
Section V, are only adapted in the training procedure and thus have the same
memory overhead and online complexity as the models learned with perfect CSI.

Table I summarizes the memory overhead and computational online complexity of
all proposed approaches as well as the baseline methods. It can be seen that
the proposed approaches vary in the number of parameters and online complexity
to allow for a smooth trade-off with respect to the desirable performance and
practical system requirements. Of particular importance is the comparison to
the deep NN approach, which is adapted from [26] and directly provides a
channel estimate at the output, i.e., it does not parameterize an analytical
estimator. It becomes apparent that the proposed approaches exhibit a much
lower number of parameters as well as a lower online complexity compared to
the deep NN approach. The main reason for this is the drastic increase of
neurons in the hidden layers in the NN approach, cf. [26]; in contrast, the
proposed models are comprised of a latent space which enforces a compression,
and thus, a reduced memory and complexity overhead. As shown in the following
numerical results, the estimation performance of the proposed approaches is,
in most cases, even better, although having fewer parameters and reduced
online complexity.

##  VIII Achievable Rate Lower Bound

The achievable rate is of great interest in quantized systems [22, 3]. We
evaluate a lower bound on the corresponding achievable rate of a respective
data transmission system that is taking the CSI mismatch into account. To this
end, after estimating the channel with the pilot transmission in (1), the data
symbol $s$ is transmitted over the same channel, i.e.,
${\bm{r}}=Q_{B}({\bm{h}}s+{\bm{n}})={\bm{B}}{\bm{h}}s+{\bm{q}}$; in the second
equation, the linearized model with Bussgang’s decomposition is used where
${\bm{q}}={\bm{B}}{\bm{n}}+{\bm{\eta}}$. We make the worst-case assumption
that the aggregated noise is Gaussian, i.e.,
${\bm{q}}\sim\mathcal{N}_{\mathbb{C}}({\bm{0}},{\bm{C}}_{{\bm{q}}}={\bm{C}}_{{\bm{r}}}-{\bm{B}}{\bm{C}}_{{\bm{h}}}{\bm{B}}^{\operatorname{H}})$,
cf. [58]. Furthermore, the BS is assumed to perform maximum-ratio combining
(MRC) with the normalized filter
${\bm{g}}_{\text{MRC}}^{\operatorname{H}}=\hat{{\bm{h}}}^{\operatorname{H}}/\|\hat{{\bm{h}}}\|_{2}^{2}$.
Note that the variance of the data symbol $s$ is assumed to be one without
loss of generality. We further assume that the SNR is the same during pilot
and data transmission. Thus, we can evaluate the use-and-then-forget (UatF)
bound as a lower bound on the achievable rate, cf. [22, Lemma 1], as

$R_{\text{UatF}}=\log_{2}\left(1+\frac{\big{|}{\operatorname{E}}[{\bm{g}}_{\text{MRC}}^{\operatorname{H}}{\bm{B}}{\bm{h}}]\big{|}^{2}}{\operatorname{var}[{\bm{g}}_{\text{MRC}}^{\operatorname{H}}{\bm{B}}{\bm{h}}]+{\bm{g}}_{\text{MRC}}^{\operatorname{H}}{\bm{C}}_{{\bm{q}}}{\bm{g}}_{\text{MRC}}}\right)$
(28)

with Monte Carlo simulations.

##  IX Simulation Results

###  IX-A Covariance Recovery

Figure 2: Performance evaluation for covariance estimation with $N=64$
dimensions and $100$ Monte Carlo iterations using covariances obtained from
the 3GPP model (6). Figure 3: MSE performance for the 3GPP channel model (cf.
Section II-3) with one (top) and three (bottom) propagation clusters for
$B\in\\{1,2,3\\}$ quantization bits, $N=64$ antennas, and $P=1$ pilot
observation (single snapshot). Figure 4: Evaluation of the achievable rate
lower bound (28) for the 3GPP channel model (cf. Section II-3) with one (left)
and three (right) propagation clusters for $B=\\{1,2,3,\infty\\}$ quantization
bits, $N=64$ antennas, and $P=1$ pilot observation. Figure 5: MSE (left) and
achievable rate lower bound (right) performance for the QuaDRiGa LOS channel
model (cf. Section II-3) for $N=64$ antennas and $P=1$ pilot.

Before investigating the channel estimation performance, we evaluate the
proposed covariance recovery algorithm from Section V-A in a purely Gaussian
setting without AWGN, comparing it to reasonable baselines. By assuming genie-
knowledge of the unquantized samples, we can evaluate the unquantized sample
covariance matrix, i.e.,
$\hat{{\bm{C}}}_{\text{unquant}}=\frac{1}{T}\sum_{t=1}^{T}{\bm{h}}_{t}{\bm{h}}_{t}^{\operatorname{H}}$.
Note that this baseline requires perfect CSI and thus only serves as a
baseline. A feasible approach is to neglect the quantization effect and
evaluate the quantized sample covariance matrix
$\hat{{\bm{C}}}_{\text{quant}}=\frac{1}{T}\sum_{t=1}^{T}{\bm{r}}_{t}{\bm{r}}_{t}^{\operatorname{H}}$,
where ${\bm{r}}_{t}=Q_{B}({\bm{h}}_{t})$. This baseline becomes more accurate
with more quantization bits $B$ but introduces a systematic error due to the
coarse quantization.

We construct $100$ random covariance matrices
${\bm{C}}_{{\bm{h}}|{\bm{\delta}}}$ from (6) and draw a fixed number of
samples from each covariance matrix as
$\\{{\bm{h}}_{t}\sim\mathcal{N}_{\mathbb{C}}({\bm{0}},{\bm{C}}_{{\bm{h}}|{\bm{\delta}}})\\}_{t=1}^{T}$.
Afterward, the normalized MSE
${\operatorname{E}}[\|{\bm{C}}_{{\bm{h}}|{\bm{\delta}}}-\hat{{\bm{C}}}\|_{F}^{2}]/{\operatorname{E}}[\|{\bm{C}}_{{\bm{h}}|{\bm{\delta}}}\|_{F}^{2}]$
is computed by using those $100$ covariance realizations.

The left plot in Fig 2 shows the normalized MSE versus different numbers $T$
of samples for different numbers $B$ of quantization bits. It can be seen that
the proposed covariance recovery scheme performs equally well for all
quantization levels since it yields a consistent estimator, i.e., the
estimation error steadily decreases for a larger number of samples. This is
similar to the unquantized sample covariance matrix (“Scov $\infty$-bit”) but
with a more or less constant offset, which is mainly caused by the correlation
estimate that is unchanged for varying bits $B$, cf. (23). In contrast, the
quantized sample covariance matrix (“Scov $B$-bit”) is a biased estimator and
shows a relatively high error floor, which decreases for a higher number of
quantization bits, as excepted. We note that the consistency and unbiasedness
of the proposed estimator is a key characteristic since, for the training with
quantized pilot observations, it can be expected that a large dataset
$\mathcal{R}$ can be acquired cheaply during regular operation of the BS; this
is in contrast to a dataset $\mathcal{H}$ consisting of ground-truth channels,
which either requires costly measurement campaigns or intricate modeling of
the underlying propagation environment.

In the right plot in Fig. 2, we evaluate the necessary number of iterations of
the Gauss-Newton algorithm for solving the nonlinear LS problem until
convergence, i.e., until the absolute change of the estimated variances is
smaller than $10^{-5}$. It can be seen that in all cases, only a few
iterations are necessary for the convergence; the number of iterations also
decreases for a higher number of data samples $T$ and for less quantization
bits (due to an increasing number of equations for larger $B$), which makes
the variance estimation very fast. In combination with the closed-form
solution for the correlation matrix (23), the covariance estimator exhibits
low complexity and is applicable for any given number of quantization bits.

###  IX-B Channel Estimation

This section provides numerical results to evaluate the proposed channel
estimators, cf. Section IV and Section V, against the discussed state-oft-the-
art baselines from Section VI. In all simulations, we have fixed the number of
training samples for both $\mathcal{H}$ and $\mathcal{R}$, cf. Section II-4,
to $T=100{,}000$. The normalized MSE
${\operatorname{E}}[\|{\bm{h}}-\hat{{\bm{h}}}\|_{2}^{2}]/{\operatorname{E}}[\|{\bm{h}}\|_{2}^{2}]$
and the achievable rate lower bound (28) are computed by means of $10{,}000$
channel samples, which are not part of the training dataset. If not otherwise
stated, the number of components for the GMM/MFA is $K=64$, the latent
dimension for the MFA/VAE is $L=N/4$, and the data-aided approaches are
trained with the channel dataset $\mathcal{H}$. For both the VAE and the DNN
approach, a single NN architecture is trained for the whole SNR range of
$[-10,20]$dB. Since a low pilot overhead is considered to be a key aspect in
practical systems [2], we, therefore, especially focus on the single snapshot
scenario in this paper.

In Fig. 3, we evaluate the MSE performance of the proposed channel estimators
in comparison to the baseline methods over the SNR for the 3GPP channel model
from Section II-3 with one (top row) and three (bottom row) propagation
clusters for $B\in\\{1,2,3\\}$ quantization bits, $N=64$ BS antennas, and
$P=1$ pilot. In all cases, the approaches “BLS”, “Buss-Scov”, and “EM-GM-GAMP”
are outperformed with a considerable performance gap over the whole SNR range
by the proposed approaches. This is due to the fact that the Bussgang theorem
does not hold for the non-Gaussian distributed channels, which is assumed by
“Buss-Scov”; moreover, the channels are generally not perfectly sparse in the
angular domain (leakage effect), which substantially impacts the CS approach
“EM-GM-GAMP”. Interestingly, for the case of one propagation cluster, the GMM-
based approach is close to the “Buss-genie” approach, which is the Bussgang
estimator with utopian knowledge of the true channel covariance matrix for a
single snapshot; this underlines the powerful estimation abilities of the GMM.
For the considered case of a ULA, the Toeplitz-structured GMM version is
almost on par with the full GMM approach, whereas the circulant-structured
approach exhibits a small performance gap. The observation of having
increasing MSE for higher SNR values beyond a certain SNR level in some cases
is due to stochastic resonance, which is a well-known effect in quantized
systems [59]; thereby, the effect can vary for different estimators, depending
on the parameterization. The “DNN” approach also shows good estimation results
for the different scenarios but is consistently outperformed by at least one
of the proposed estimators, although having a much larger number of parameters
and a higher online complexity, cf. Table I. Overall, the simulation results
in Fig. 3 demonstrate the great potential of the proposed class of
parameterized estimators based on Gaussian latent models in combination with
the Bussgang estimator.

Fig. 4 assesses the achievable rate lower bound from (28) for the 3GPP channel
model with one (left) and three (right) propagation clusters for $N=64$
antennas, $P=1$ pilot observation, and different numbers of quantization bits.
For the case of one propagation cluster, the achievable rate lower bound of
the GMM and VAE approach is almost on par with the “Buss-genie” approach.
Moreover, a substantial gap to the achievable rate lower bound of the “Buss-
Scov” approach is apparent for all considered numbers of quantization bits.
This behavior similarly translates to the case of three propagation clusters
but with an overall reduced gap to the baseline approach. This result
indicates that the better estimation performance of the proposed estimators
can be effectively converted to a higher data rate or to a lower resolution
while preserving the same throughput as the baseline approach “Buss-Scov”.

Figure 6: MSE performance for the QuaDRiGa LOS channel model (cf. Section
II-3) for $B=1$. Left: $P=1$ and $\text{SNR}=10$dB; right: $N=64$ and
$\text{SNR}=5$dB. Figure 7: MSE performance for the QuaDRiGa LOS and mixed
LOS/NLOS channel model (cf. Section II-3) for $B=1$ (left) and $B=3$ (right),
$N=64$, and $P=1$. Figure 8: MSE performance evaluation of the models trained
with quantized data for the 3GPP channel model (cf. Section II-3) for $P=1$
pilot and $N=64$ antennas.

In Fig. 5 (left), the MSE performance is compared for different numbers $B$ of
quantization bits, now for the QuaDRiGa LOS channel model, cf. Section II-3.
Once again, the approaches “BLS”, “Buss-Scov”, and “EM-GM-GAMP” are
outperformed over the whole range of quantization bits. Interestingly, the
“DNN” approach is comparably good for $B=2$ and $B=3$ but, in turn, suffers in
performance in the extreme cases of $B=1$ and infinite resolution, which
indicates better robustness of the proposed approaches in comparison.
Interestingly, the MFA estimator is ranked among the best estimators in this
case, which highlights the fact that estimators’ performances may vary
slightly for different channel models; however, it can also be seen that the
overall performance of the proposed class of estimators is stable and robust
with respect to a different channel model. In Fig. 5 (right), the
corresponding achievable rate lower bound from (28) is evaluated. It can be
seen that the better estimation qualities of the proposed approaches in terms
of the MSE directly translate to a higher achievable rate guarantee, which is
approximately only 1–2 bits/s/Hz below the achievable rate lower bound with
perfect CSI knowledge at the receiver.

The left plot in Fig. 6 examines the estimation quality over different numbers
$N$ of antennas at the BS for $B=1$ bit, $P=1$ pilot, and an SNR of $10$dB for
the QuaDRiGa LOS channel model, cf. Section II-3. In contrast to the baseline
approaches “BLS”, “Buss-Scov”, and “EM-GM-GAMP”, the estimation performance of
the proposed estimators significantly increases for a higher number of
antennas, which is particularly important in massive MIMO systems. The GMM
approach performs best for all considered antenna numbers, whereas the VAE
approach is especially strong in the high number of antennas case; this is
reasoned by the circulant parameterization of the covariances by the VAE,
which only holds asymptotically for high numbers of antennas. A similar
behavior is observed for the “GMM-circ” estimator. The “DNN” approach, which
has a quartic scaling of the number of parameters in the number of antennas,
cf. Table I, is outperformed over the whole range.

The right plot in Fig. 6 shows the MSE performance for an increasing number of
pilot observations by utilizing the pilot design from Section II-1 for $N=64$
antennas and a fixed SNR of $5$dB. The proposed estimators outperform all
baseline approaches over the whole range of pilot observations, including the
“DNN” estimator. Especially the MFA and VAE models, which are comprised of a
nonlinear latent space, perform well in the high number of pilots regime.
Although we focus our analysis primarily on the single-snapshot case, we see
that also with an increasing number of pilots, the proposed estimators perform
very well.

Next, the number $K$ of GMM components that are necessary to achieve a certain
performance is discussed. In the left plot of Fig. 7, the MSE over the number
of GMM components is evaluated for $B=1$ bit, $N=64$ antennas, $P=1$ pilot,
and for varying SNRs for the QuaDRiGa LOS as well as mixed LOS/NLOS channel
model, cf. Section II-3. It can be expected that the superposition of many
sub-paths, as it is the case in the mixed LOS/NLOS scenario, results in a less
structured wireless channel, and thus, less structural information can be
inferred as prior knowledge by the data-aided models. Therefore, it can be
observed that the overall performance is worse for the mixed scenario.
However, in both scenarios, the increase of GMM components continuously
enhances the estimation performance, with a greater improvement in the pure
LOS case. This points towards applications in mmWave communications, where the
high frequency in combination with smaller BS cells results in high LOS
probabilities.

The right plot of Fig. 7 analyzes the same setup but now for $B=3$
quantization bits. In this case, we evaluate the approximation quality of
computing ${\bm{C}}_{{\bm{r}}}$ for a given ${\bm{C}}_{{\bm{y}}}$ via (13) by
comparing it with the estimator that computes the exact variances via (12) and
otherwise uses the same approximation for the off-diagonals, labeled “GMM-
ex.”. As expected, the approximation is highly accurate in the low to medium
SNR region, which is the considered operating range of low-resolution systems.
In high SNR, the approximation is less accurate and shows saturation effects
when increasing the number of GMM components. However, the overall
approximation loss is small, and it still results in a high estimation
accuracy as compared to the baseline approaches. Besides that, the estimation
performance generally steadily increases for a higher number of GMM components
with an overall saturation for high numbers $K$ of components.

Fig. 8 evaluates the proposed training adaptations for the GMM and the VAE, as
detailed in Section V, in order to learn from noisy and quantized pilot
observations $\mathcal{R}$, cf. Section II-4, without having ground-truth
channel samples during training. We refer to the adapted GMM, cf. Algorithm 1,
as “GMM $\mathcal{R}$”, and the GMM learned with ground-truth channels as “GMM
$\mathcal{H}$”. The VAEs are denoted likewise. Note that a main difference
between “GMM $\mathcal{R}$” and “VAE $\mathcal{R}$ is that the training of the
adapted GMM is performed for a fixed SNR, but afterward, the model can be
utilized for the whole SNR range, whereas the VAE is trained directly for the
whole SNR range in order to generalize properly. In order to have a meaningful
comparison, the “GMM $\mathcal{R}$” is trained for each SNR point, whereas the
“VAE $\mathcal{R}$” is trained for the whole SNR range (no performance gain
was seen in the simulation results for an SNR-dependent training).

In Fig. 8 (a) and (b), the case of $N=64$, $P=1$, $B=2$, and the 3GPP channel
model with one and three propagation clusters are considered, respectively,
whereas in Fig. 8 (c) and (d) the same setup with $B=3$ bits is investigated.
Astonishingly, through the model-based adaptations, the models trained with
the dataset $\mathcal{R}$ consisting of coarsely quantized and noisy data
samples are almost on par with their counterparts that are trained with
ground-truth noise-free channel samples from $\mathcal{H}$. Overall, the
estimation quality seems to be most accurate in the low to medium SNR range,
where the approximation in (13) and thus the training adaptations are highly
accurate. This implies that the generally costly dataset $\mathcal{H}$ can be
replaced by $\mathcal{R}$ with almost no performance loss in this regime. In
the high SNR regime, the performance loss of “GMM $\mathcal{R}$” tends to
increase, which is a consequence of the generally indefinite closed-form
correlation estimate (23). After the projection onto the set of PSD matrices
by truncating the negative eigenvalues, cf. Algorithm 1, the resulting
covariance estimate is missing the corresponding eigenvectors, which has a
higher impact on the performance in the high SNR regime. This correlates with
the observation that the adaptations are working best in cases with fewer
multi-path components.

##  X Conclusion

In this work, we presented a novel and promising approach for channel
estimation in coarse quantization systems by utilizing Gaussian latent models
such as GMMs, MFAs, and VAEs. These models successfully learn the unknown and
complex channel distributions present in radio propagation scenarios and,
afterward, utilize this valuable prior information to enable the development
of tractable parameterized linear MMSE estimators based on a conditional
Bussgang decomposition. We have shown that all of the presented estimators
perform well for various channel and system parameters with only minor
differences. This allows for selecting the preferred model using the discussed
memory and complexity overhead.

In addition, we derived model-based training adaptations, i.e., a covariance
recovery algorithm for the GMM and a loss function adaptation for the VAE, in
order to learn these models directly from quantized training data, with only
marginal performance losses. Extensive simulations verified a superior
performance over classical and deep learning-based approaches in terms of MSE
and achievable rate metrics.

The analysis of multi-user systems, especially in the presence of coarsely
quantized training data, is left as future work.

##  XI Acknowledgment

The authors gratefully acknowledge valuable discussions with and input from
Dr.-Ing. Michael Koller in the early stages of this work.

## References

  * [1] J. Liu, Z. Luo, and X. Xiong, “Low-Resolution ADCs for Wireless Communication: A Comprehensive Survey,” _IEEE Access_ , vol. 7, pp. 91 291–91 324, 2019.
  * [2] J. Zhang, L. Dai, X. Li, Y. Liu, and L. Hanzo, “On Low-Resolution ADCs in Practical 5G Millimeter-Wave Massive MIMO Systems,” _IEEE Commun. Magazine_ , vol. 56, no. 7, pp. 205–211, 2018.
  * [3] S. Jacobsson, G. Durisi, M. Coldrey, U. Gustavsson, and C. Studer, “Throughput Analysis of Massive MIMO Uplink With Low-Resolution ADCs,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 6, pp. 4038–4051, 2017.
  * [4] A. Mezghani and J. Nossek, “Capacity Lower Bound of MIMO Channels with Output Quantization and Correlated Noise,” in _Int. Symp. Inf. Theory_ , 2012.
  * [5] A. Mezghani, “Information-theoretic Analysis and Signal Processing Techniques for Quantized MIMO Communications,” Ph.D. dissertation, Dept. Elect. Eng., Tech. Univ. Munich, Munich, Germany, 2016.
  * [6] J. Singh, O. Dabeer, and U. Madhow, “On the Limits of Communication with Low-Precision Analog-to-Digital Conversion at the Receiver,” _IEEE Trans. Commun._ , vol. 57, no. 12, pp. 3629–3639, 2009.
  * [7] J. A. Nossek and M. T. Ivrlač, “Capacity and Coding for Quantized MIMO Systems,” in _Proc. Int. Conf. Wireless Commun. Mobile Comput._ , 2006, p. 1387–1392.
  * [8] S. Wang, Y. Li, and J. Wang, “Multiuser Detection in Massive Spatial Modulation MIMO With Low-Resolution ADCs,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 4, pp. 2156–2168, 2015.
  * [9] A. Mezghani, F. Antreich, and J. A. Nossek, “Multiple Parameter Estimation with Quantized Channel Output,” in _Int. ITG Workshop Smart Antennas (WSA)_ , 2010, pp. 143–150.
  * [10] R. P. David and J. Cal-Braz, “Feedback-Controlled Channel Estimation with Low-resolution ADCs in Multiuser MIMO Systems,” in _IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP)_ , 2019, pp. 4674–4678.
  * [11] F. Liu, X. Shang, Y. Cheng, and G. Zhang, “Computationally Efficient Maximum Likelihood Channel Estimation for Coarsely Quantized Massive MIMO Systems,” _IEEE Wireless Commun. Lett._ , vol. 26, no. 2, pp. 444–448, 2022.
  * [12] H. Lee, Y.-S. Jeon, H. Do, and N. Lee, “Concatenated MMSE Estimation for Quantized OFDM Systems,” in _IEEE Int. Conf. Commun._ , 2019.
  * [13] O. Dabeer and U. Madhow, “Channel Estimation with Low-Precision Analog-to-Digital Conversion,” in _IEEE Int. Conf. Commun._ , 2010.
  * [14] B. Fesl, M. Koller, and W. Utschick, “On the Mean Square Error Optimal Estimator in One-Bit Quantized Systems,” _IEEE Trans. Signal Process._ , vol. 71, pp. 1968–1980, 2023.
  * [15] S. S. Thoota and C. R. Murthy, “Variational Bayes’ Joint Channel Estimation and Soft Symbol Decoding for Uplink Massive MIMO Systems With Low Resolution ADCs,” _IEEE Trans. Commun._ , vol. 69, no. 5, pp. 3467–3481, 2021.
  * [16] P. Sun, Z. Wang, R. W. Heath, and P. Schniter, “Joint Channel-Estimation/Decoding With Frequency-Selective Channels and Few-Bit ADCs,” _IEEE Trans. Signal Process._ , vol. 67, no. 4, pp. 899–914, 2019\.
  * [17] C.-K. Wen, C.-J. Wang, S. Jin, K.-K. Wong, and P. Ting, “Bayes-Optimal Joint Channel-and-Data Estimation for Massive MIMO With Low-Precision ADCs,” _IEEE Trans. Signal Process._ , vol. 64, no. 10, pp. 2541–2556, 2016.
  * [18] L. Jacques, K. Degraux, and C. D. Vleeschouwer, “Quantized Iterative Hard Thresholding: Bridging 1-Bit and High-Resolution Quantized Compressed Sensing,” in _Proc. Int. Conf. Sampling Theory Appl. (SampTA)_ , 2013, p. 105–108.
  * [19] A. Kaushik, E. Vlachos, J. Thompson, and A. Perelli, “Efficient Channel Estimation in Millimeter Wave Hybrid MIMO Systems with Low Resolution ADCs,” in _26th Eur. Signal Process. Conf. (EUSIPCO)_ , 2018, pp. 1825–1829.
  * [20] J. Mo, P. Schniter, and R. W. Heath, “Channel Estimation in Broadband Millimeter Wave MIMO Systems With Few-Bit ADCs,” _IEEE Trans. Signal Process._ , vol. 66, no. 5, pp. 1141–1154, 2018.
  * [21] J. J. Bussgang, “Crosscorrelation Functions of Amplitude-Distorted Gaussian Signals,” MIT Res. Lab. Electron., Tech. Rep. 216, 1952.
  * [22] Y. Li, C. Tao, G. Seco-Granados, A. Mezghani, A. L. Swindlehurst, and L. Liu, “Channel Estimation and Performance Analysis of One-Bit Massive MIMO Systems,” _IEEE Trans. Signal Process._ , vol. 65, no. 15, pp. 4075–4089, 2017.
  * [23] S. Rao, G. Seco-Granados, H. Pirzadeh, J. A. Nossek, and A. L. Swindlehurst, “Massive MIMO Channel Estimation With Low-Resolution Spatial Sigma-Delta ADCs,” _IEEE Access_ , vol. 9, pp. 109 320–109 334, 2021.
  * [24] E. Balevi, A. Doshi, A. Jalal, A. Dimakis, and J. G. Andrews, “High Dimensional Channel Estimation Using Deep Generative Networks,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 1, pp. 18–30, 2021.
  * [25] B. Fesl, M. Koller, N. Turan, and W. Utschick, “Learning a Low-Complexity Channel Estimator for One-Bit Quantization,” in _54th Asilomar Conf. Signals, Syst., Comput._ , 2020, pp. 393–397.
  * [26] Y. Zhang, M. Alrabeiah, and A. Alkhateeb, “Deep Learning for Massive MIMO With 1-Bit ADCs: When More Antennas Need Fewer Pilots,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 8, pp. 1273–1277, 2020.
  * [27] Y. Dong, H. Wang, and Y.-D. Yao, “Channel Estimation for One-Bit Multiuser Massive MIMO Using Conditional GAN,” _IEEE Commun. Lett._ , vol. 25, no. 3, pp. 854–858, 2021.
  * [28] N. Turan, B. Fesl, M. Koller, M. Joham, and W. Utschick, “A Versatile Low-Complexity Feedback Scheme for FDD Systems via Generative Modeling,” 2023, arXiv preprint: 2304.14373.
  * [29] M. Koller, B. Fesl, N. Turan, and W. Utschick, “An Asymptotically MSE-Optimal Estimator Based on Gaussian Mixture Models,” _IEEE Trans. Signal Process._ , vol. 70, pp. 4109–4123, 2022.
  * [30] N. Turan, B. Fesl, M. Grundei, M. Koller, and W. Utschick, “Evaluation of a Gaussian Mixture Model-based Channel Estimator using Measurement Data,” in _Int. Symp. Wireless Commun. Syst. (ISWCS)_ , 2022.
  * [31] B. Fesl, M. Joham, S. Hu, M. Koller, N. Turan, and W. Utschick, “Channel Estimation based on Gaussian Mixture Models with Structured Covariances,” in _56th Asilomar Conf. Signals, Syst., Comput._ , 2022, pp. 533–537.
  * [32] B. Fesl, N. Turan, and W. Utschick, “Low-Rank Structured MMSE Channel Estimation with Mixtures of Factor Analyzers,” in _57th Asilomar Conf. Signals, Syst., Comput._ , 2023, to be published.
  * [33] M. Baur, B. Fesl, and W. Utschick, “Leveraging Variational Autoencoders for Parameterized MMSE Channel Estimation,” 2023, arXiv preprint: 2307.05352.
  * [34] M. Baur, B. Fesl, M. Koller, and W. Utschick, “Variational Autoencoder Leveraged MMSE Channel Estimation,” in _56th Asilomar Conf. Signals, Syst., Comput._ , 2022, pp. 527–532.
  * [35] J. Max, “Quantizing for Minimum Distortion,” _IRE Trans. Inf. Theory_ , vol. 6, no. 1, pp. 7–12, 1960.
  * [36] L. Liu, S. Xue, Y. Ma, N. Yi, and R. Tafazolli, “On the Design of Quantization Functions for Uplink Massive MIMO with Low-Resolution ADCs,” in _IEEE 93rd Veh. Technol. Conf._ , 2021, pp. 1–5.
  * [37] R. Gray and D. Neuhoff, “Quantization,” _IEEE Trans. Inf. Theory_ , vol. 44, no. 6, pp. 2325–2383, 1998.
  * [38] 3GPP, “Spatial Channel Model for Multiple Input Multiple Output (MIMO) Simulations,” 3rd Generation Partnership Project (3GPP), Tech. Rep. 25.996 (V16.0.0), Jul. 2020.
  * [39] D. Neumann, T. Wiese, and W. Utschick, “Learning the MMSE Channel Estimator,” _IEEE Trans. Signal Process._ , vol. 66, no. 11, pp. 2905–2917, Jun. 2018.
  * [40] S. Jaeckel, L. Raschkowski, K. Börner, and L. Thiele, “QuaDRiGa: A 3-D Multi-Cell Channel Model With Time Evolution for Enabling Virtual Field Trials,” _IEEE Trans. Antennas Propag._ , vol. 62, no. 6, pp. 3242–3256, 2014.
  * [41] S. Jaeckel, “Quasi-Deterministic Channel Modeling and Experimental Validation in Cooperative and Massive MIMO Deployment Topologies,” Ph.D. dissertation, Dept. Elect. Eng., Tech. Univ. Ilmenau, Ilmenau, Germany, 2017.
  * [42] A. Papoulis and S. U. Pillai, _Probability, Random Variables and Stochastic Processes_. McGraw-Hill Education, 2002.
  * [43] S. Jacobsson, G. Durisi, M. Coldrey, T. Goldstein, and C. Studer, “Quantized Precoding for Massive MU-MIMO,” _IEEE Trans. Commun._ , vol. 65, no. 11, pp. 4670–4684, 2017.
  * [44] G. Jacovitti and A. Neri, “Estimation of the Autocorrelation Function of Complex Gaussian Stationary Processes by Amplitude Clipped Signals,” _IEEE Trans. Inf. Theory_ , vol. 40, no. 1, pp. 239–245, 1994.
  * [45] W. J. Cody, “Rational Chebyshev Approximations for the Error Function,” _Mathematics of Computation_ , vol. 23, no. 107, pp. 631–637, 1969.
  * [46] P. Borjesson and C.-E. Sundberg, “Simple Approximations of the Error Function Q(x) for Communications Applications,” _IEEE Trans. Commun._ , vol. 27, no. 3, pp. 639–643, 1979.
  * [47] C. M. Bishop, _Pattern Recognition and Machine Learning (Information Science and Statistics)_. Springer New York, 2006.
  * [48] W. Pan and X. Shen, “Penalized Model-Based Clustering with Application to Variable Selection,” _J. Mach. Learn. Research_ , vol. 8, no. 41, pp. 1145–1164, 2007.
  * [49] K. P. Murphy, _Machine Learning: A Probabilistic Perspective_. The MIT Press, 2012.
  * [50] D. P. Kingma and M. Welling, “Auto-Encoding Variational Bayes,” in _Proc. 2nd Int. Conf. Learn. Represent._ , 2014.
  * [51] D. J. Rezende, S. Mohamed, and D. Wierstra, “Stochastic Backpropagation and Approximate Inference in Deep Generative Models,” in _Proc. 31st Int. Conf. Mach. Learn._ , 2014.
  * [52] B. Fesl, N. Turan, M. Joham, and W. Utschick, “Learning a Gaussian Mixture Model from Imperfect Training Data for Robust Channel Estimation,” _IEEE Wireless Commun. Lett._ , 2023.
  * [53] S. Dirksen, J. Maly, and H. Rauhut, “Covariance Estimation under One-bit Quantization,” _Proc. Appl. Math. Mech._ , vol. 21, no. 1, 2021.
  * [54] Y.-H. Xiao, L. Huang, D. Ramírez, C. Qian, and H. C. So, “One-Bit Covariance Reconstruction with Non-zero Thresholds: Algorithm and Performance Analysis,” 2023, arXiv preprint: 2303.16455.
  * [55] A. Eamaz, F. Yeganegi, and M. Soltanalian, “Covariance Recovery for One-Bit Sampled Non-Stationary Signals With Time-Varying Sampling Thresholds,” _IEEE Trans. Signal Process._ , vol. 70, pp. 5222–5236, 2022.
  * [56] T. Yang, J. Maly, S. Dirksen, and G. Caire, “Plug-in Channel Estimation with Dithered Quantized Signals in Spatially Non-Stationary Massive MIMO Systems,” 2023, arXiv preprint: 2301.04641.
  * [57] Y. Xu and E. W. Cheney, “Strictly Positive Definite Functions on Spheres,” _Proc. Amer. Math. Soc._ , vol. 116, no. 4, pp. 977–981, 1992.
  * [58] B. Hassibi and B. Hochwald, “How Much Training is Needed in Multiple-Antenna Wireless Links?” _IEEE Trans. Inf. Theory_ , vol. 49, no. 4, pp. 951–963, 2003.
  * [59] M. D. McDonnell, N. G. Stocks, C. E. M. Pearce, and D. Abbott, _Stochastic Resonance: From Suprathreshold Stochastic Resonance to Stochastic Signal Quantization_. Cambridge University Press, 2008.