# On the Optimal MMSE Channel Estimation for One-Bit Quantized MIMO Systems

Minhua Ding, Italo Atzeni,
Antti Tölli, and A. Lee Swindlehurst M. Ding, I. Atzeni, and A. Tölli are with
the Centre for Wireless Communications, University of Oulu, Finland (e-mails:
{minhua.ding, italo.atzeni, antti.tolli}@oulu.fi). A. L. Swindlehurst is with
the Center for Pervasive Communications and Computing, University of
California, Irvine, CA, USA (e-mail: swindle@uci.edu). This work was supported
by the Research Council of Finland (318927 6G Flagship, 336449 Profi6, 348396
HIGH-6G, and 357504 EETCAMD) and by the U. S. National Science Foundation
(CCF-2225575)

###### Abstract

This paper focuses on the minimum mean squared error (MMSE) channel estimator
for multiple-input multiple-output (MIMO) systems with one-bit quantization at
the receiver side. Despite its optimality and significance in estimation
theory, the MMSE channel estimator has not been fully investigated in this
context due to its general non-linearity and computational complexity.
Instead, the typically suboptimal Bussgang linear MMSE (BLMMSE) estimator has
been widely adopted. In this work, we develop a new framework to compute the
MMSE channel estimator that hinges on computation of the orthant probability
of the multivariate normal distribution. Based on this framework, we determine
a necessary and sufficient condition for the BLMMSE channel estimator to be
optimal and equivalent to the MMSE estimator. Under the assumption of specific
channel correlation or pilot symbols, we further utilize the framework to
derive analytical expressions for the MMSE channel estimator that are
particularly convenient for computation when certain system dimensions become
large, thereby enabling a comparison between the BLMMSE and MMSE channel
estimators in these cases.

###### Index Terms:

Channel estimation, massive MIMO, minimum mean squared error (MMSE), one-bit
quantization, orthant probability.

##  I Introduction

Due to its potential to significantly enhance spectral efficiency and
reliability [1, 2, 3], state-of-the-art wireless networks employ multiple-
input multiple-output (MIMO) transceivers with large numbers of antennas, and
such systems are envisioned to assume an even more prominent role in future
wireless systems [4]. However, deploying transceivers with many antennas
motivates study of the use and impact of low-cost, energy-efficient hardware
on system performance metrics.

One of the main challenges in implementing a fully digital high-resolution
massive MIMO system is the prohibitively high power consumption of the analog-
to-digital/digital-to-analog converters (ADCs/DACs), which scales
exponentially with the number of resolution bits [5, 6, 7, 8]. Against this
backdrop, there is a large body of literature investigating massive MIMO
systems with low-resolution ADCs [6, 9, 7, 8, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19]. In particular, one-bit quantization for massive MIMO systems has
received special attention due to its energy efficiency and low complexity.
The capacity of MIMO channels with one-bit ADCs was investigated in [6]
assuming perfect channel state information (CSI) at both the transmitter and
the receiver, whereas capacity bounds for one-bit MIMO Gaussian channels with
analog combining were studied in [12]. Assuming perfect CSI at the receiver,
one-bit MIMO-OFDM detection based on inexact variants of the classical
expectation-maximization method was considered in [13].

Most massive MIMO systems require CSI to facilitate beamforming design.
Channel estimation based on one-bit quantized measurements is thus of great
importance [7, 14, 15, 16, 17, 19]. In this respect, an approach based on
maximum likelihood (ML) was used for channel estimation in connection with
near-ML data detection in [14], while the Cramér-Rao performance bounds for
channel estimation in one-bit quantized massive MIMO were studied in [15].
Channel estimation from one-bit quantized measurements with an unknown
threshold was considered in [16]. The Bussgang linear minimum mean squared
error (BLMMSE) channel estimator proposed in [7] elegantly used the second-
order statistics of one-bit quantized signals [20, 21], and its linear
structure inspired subsequent signal processing methods applied to, e.g., data
detection [17, 18] and precoding design [22, 23].

Despite its optimality in terms of mean squared error (MSE) and its
significance in estimation theory, the minimum MSE (MMSE) channel estimator,
also known as the conditional mean estimator (CME) [19, 24], [25, p. 17,
(2.36)], has not been well investigated for one-bit quantized MIMO systems.
The authors of [19] established the equivalence between the BLMMSE and the
optimal MMSE channel estimates in single-input single-output (SISO) systems
with one pilot symbol, as well as in noiseless SISO systems with multiple
(specifically chosen) pilot symbols. However, some questions remain unanswered
regarding the CME for one-bit quantized MIMO systems. Furthermore, explicit
computation of the MMSE estimate still appears elusive due to the required
multi-dimensional integration involving a multivariate normal (MVN)
distribution [19, 24].

In this work, we aim at providing further results on the optimal MMSE channel
estimator for MIMO systems with one-bit ADCs. We consider a point-to-point
MIMO system to focus on the fundamental estimation problem, although our
results can be framed for a multi-user uplink system with single-antenna users
and a multi-antenna base station (BS). Despite the complexity of the
estimation problem, we obtain several new results, which are summarized as
follows.

  1. 1.

We first develop a general framework for computing the MMSE channel estimate,
which relies on computation of the _orthant probability_ of an MVN
distribution. The Gaussian orthant probability has been investigated in the
literature [26, 27, 28, 29, 30, 31] and used in various engineering fields
[32, 33]. Therefore, our framework connects MMSE estimation to a wider class
of problems.

  2. 2.

Based on the above framework, we determine a necessary and sufficient
condition for the equivalence of the BLMMSE and MMSE channel estimates in one-
bit quantized MIMO systems. Moreover, we examine specific cases where this
equivalence holds.

  3. 3.

The simplest case where the BLMMSE channel estimate becomes suboptimal is
identified to correspond to a single-input multiple-output (SIMO) system with
real-valued full channel correlation among three receive antennas. For this
case, we present a closed-form expression of the optimal (non-linear) MMSE
estimate.

  4. 4.

We derive an expression for the MMSE channel estimate for SIMO channels with a
single pilot symbol assuming real-valued equal correlation among all the
receive antennas. In addition, we derive an expression for the optimal MMSE
channel estimate for SISO channels or spatially white SIMO channels with
multiple real-valued pilot symbols. These formulas only require one-
dimensional integration and can thus be efficiently computed using common
software when either the number of receive antennas or the number of pilot
symbols becomes large.

_Outline._ The rest of the paper is organized as follows. Section II describes
the system model and the goal of our investigation. Section III develops the
framework for computing the MMSE channel estimate in one-bit quantized MIMO
systems. After providing background on the Gaussian orthant probability in
Section IV, the condition for the optimality of the BLMMSE channel estimate is
presented in Section V. Various results on the MMSE channel estimate are
discussed in Section VI, whereas simulation results are provided in Section
VII. Finally, Section VIII concludes the paper.

_Notation._ Lowercase letters, bold lowercase letters, and bold uppercase
letters denote scalars, vectors, and matrices, respectively. The quantity
$x_{k}$ represents the $k$-th element of vector $\mathbf{x}$ and
$[\mathbf{X}]_{ik}$ denotes the $(i,k)$-th element of matrix $\mathbf{X}$. We
let $(\cdot)^{*},(\cdot)^{\mathrm{T}}$, and $(\cdot)^{\mathrm{H}}$ denote
complex conjugate, transpose, and Hermitian transpose, respectively. The
$n\times n$ identity matrix is written as $\mathbf{I}_{n}$. The operator
$\otimes$ represents the Kronecker product, $|\mathbf{X}|$ denotes the
determinant of a square matrix $\mathbf{X}$, $\textrm{vec}(\cdot)$ denotes
vectorization, $\|\mathbf{x}\|=\sqrt{\mathbf{x}^{\mathrm{H}}\mathbf{x}}$,
$|x|$ denotes the modulus of $x$, and $j=\sqrt{-1}$. The sets of real-valued
non-negative, real-valued, and complex-valued $n$-dimensional vectors are
denoted by $\mathbb{R}_{+}^{n}$, $\mathbb{R}^{n}$, and $\mathbb{C}^{n}$,
respectively, whereas the sets of real- and complex-valued $m\times n$
matrices are denoted by $\mathbb{R}^{m\times n}$ and $\mathbb{C}^{m\times n}$,
respectively. Real-valued and circularly symmetric complex-valued MVN
distributions with mean $\bm{\mu}$ and covariance matrix $\bm{\Sigma}$ are
respectively indicated by $\mathcal{N}(\bm{\mu},\bm{\Sigma})$ and
$\mathcal{CN}(\bm{\mu},\bm{\Sigma})$. The operators $\mathsf{E}(\cdot)$ and
$\Pr(\cdot)$ represent statistical expectation and probability, respectively.
In addition, $[z_{k}]_{k=1}^{m}$ represents the vector with elements
$z_{1},\ldots,z_{m}$, $\mathtt{diag}{{(}}\mathbf{z})$ is a diagonal matrix
with the elements of $\mathbf{z}$ on its diagonal,
$\widetilde{\mathtt{diag}{{}}}(\mathbf{X})$ is the diagonal matrix formed
using the diagonal elements of the arbitrary square matrix $\mathbf{X}$, and
$\mathtt{DIAG}{{(}}\mathbf{Y}_{1},\ldots,\mathbf{Y}_{n})$ is a block-diagonal
matrix with $\mathbf{Y}_{1},\ldots,\mathbf{Y}_{n}$ representing the blocks
along the diagonal.

##  II System Model and Problem Statement

###  II-A System Model

Consider a general MIMO system with $N_{T}$ transmit antennas and $N_{R}$
receive antennas. Assume slow flat Rayleigh fading, and let
$\mathbf{H}\in\mathbb{C}^{N_{R}\times N_{T}}$ denote the channel matrix
between the transmitter and the receiver.

For channel estimation, $\tau$ pilot symbols per antenna are transmitted,
which are collectively denoted by
$\mathbf{S}^{\mathrm{T}}=[\mathbf{s}_{1}\;\ldots\;\mathbf{s}_{\tau}]\in\mathbb{C}^{N_{T}\times\tau}$.
We define the received signal (at the input of the ADCs) as

$\displaystyle\mathbf{B}=\mathbf{H}\mathbf{S}^{\mathrm{T}}+\mathbf{N}=[\mathbf{b}_{1}\;\ldots\;\mathbf{b}_{\tau}]\in\mathbb{C}^{N_{R}\times\tau},$
(1)

where
$\mathbf{N}=[\mathbf{n}_{1}\;\ldots\;\mathbf{n}_{\tau}]\in\mathbb{C}^{N_{R}\times\tau}$
represents the noise matrix. Furthermore, we introduce

$\displaystyle\mathbf{b}=\textrm{vec}(\mathbf{B})=(\mathbf{S}\otimes\mathbf{I}_{N_{R}})\mathbf{h}+\mathbf{n}=\mathbf{A}\mathbf{h}+\mathbf{n}\in\mathbb{C}^{\tau
N_{R}},$ (2)

with $\mathbf{h}=\textrm{vec}(\mathbf{H})$ $\in\mathbb{C}^{N_{T}N_{R}}$,
$\mathbf{n}=\textrm{vec}(\mathbf{N})$ $\in\mathbb{C}^{\tau N_{R}}$, and

$\displaystyle\mathbf{A}=\mathbf{S}\otimes\mathbf{I}_{N_{R}}\in\mathbb{C}^{(\tau
N_{R})\times(N_{T}N_{R})}.$ (3)

Based on the channel assumption,
$\mathbf{h}\sim\mathcal{CN}(\mathbf{0},\bm{\Sigma})$, where $\bm{\Sigma}$
represents the $(N_{T}N_{R})\times(N_{T}N_{R})$ channel covariance matrix, and
its joint probability density function (pdf) is given by

$\displaystyle
p(\mathbf{h})=\frac{1}{\pi^{N_{T}N_{R}}|\bm{\Sigma}|}\exp\left\\{-\mathbf{h}^{\mathrm{H}}\bm{\Sigma}^{-1}\mathbf{h}\right\\}.$
(4)

The additive Gaussian noise samples are assumed to be spatially and temporally
white, i.e., $\mathbf{n}\sim\mathcal{CN}(\mathbf{0},\sigma^{2}\mathbf{I}_{\tau
N_{R}})$.

We denote the element-wise memoryless one-bit quantization operation as

$\displaystyle\mathcal{Q}_{\rm
1bit}(\cdot)=\textrm{sgn}\left(\Re(\cdot)\right)+j\textrm{sgn}\left(\Im(\cdot)\right),$
(5)

where $\textrm{sgn}(\cdot)$ denotes the sign function and $\Re(\cdot)$ (resp.
$\Im(\cdot)$) represents the real (resp. imaginary) part. The quantized
received signal (at the output of the ADCs) is given by

$\displaystyle\mathbf{r}=\mathcal{Q}_{\rm 1bit}(\mathbf{b})\in\mathbb{C}^{\tau
N_{R}}.$ (6)

###  II-B Problem Statement

Our goal is two-fold.

  1. 1.

Based on the quantized output $\mathbf{r}$, the MMSE estimate of $\mathbf{h}$
is given by $\hat{\mathbf{h}}_{\rm
MMSE}=\mathsf{E}(\mathbf{h}|\mathbf{r})=\int_{\mathbb{C}^{N_{T}N_{R}}}\mathbf{h}p(\mathbf{h}|\mathbf{r}){\rm
d}\mathbf{h}$, with
$p(\mathbf{h}|\mathbf{r})=\Pr(\mathbf{r}|\mathbf{h})p(\mathbf{h})/\Pr(\mathbf{r})$.
Therefore,

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\frac{\displaystyle\int_{\mathbb{C}^{N_{T}N_{R}}}\mathbf{h}\Pr(\mathbf{r}|\mathbf{h})p(\mathbf{h}){\rm
d}\mathbf{h}}{\Pr(\mathbf{r})},$ (7)

with

$\displaystyle\Pr(\mathbf{r})=\int_{\mathbb{C}^{N_{T}N_{R}}}\Pr(\mathbf{r}|\mathbf{h})p(\mathbf{h}){\rm
d}\mathbf{h}.$ (8)

Since there are no explicit details in the literature on how to compute (7),
we first develop a novel framework for its computation, and then we determine
a necessary and sufficient condition under which the MMSE estimate is linear,
i.e., a necessary and sufficient condition for the BLMMSE channel estimator to
be MSE-optimal.

  2. 2.

We further identify the simplest case where the MMSE estimate admits a closed-
form solution. Despite the generic complexity of computing (7) [19], under
specific assumptions on the channel correlation or pilot symbols, we aim at
obtaining computationally efficient analytical expressions for
$\hat{\mathbf{h}}_{\rm MMSE}$ that facilitate insight into its performance
when certain system dimensions, i.e., the number of receive antennas or the
number of pilot symbols, become large.

##  III A Framework for One-Bit MMSE Channel Estimation

Our goal in this section is to examine the structure of the MMSE channel
estimator with one-bit observations and reformulate the corresponding
conditional mean expression into a form that involves the Gaussian orthant
probability. This probability will then be exploited in the following section
to find a condition that can be used to determine the optimality of the BLMMSE
estimator.

Recall (2) and (5)–(6), and let
$\mathbf{r}=\mathbf{r}_{\mathsf{R}}+j\mathbf{r}_{\mathsf{I}}$, with

$\displaystyle\mathbf{r}_{\mathsf{R}}=\Re(\mathbf{r})=\textrm{sgn}\left(\Re\left(\mathbf{b}\right)\right),\quad\mathbf{r}_{\mathsf{I}}=\Im(\mathbf{r})=\textrm{sgn}\left(\Im\left(\mathbf{b}\right)\right).$
(9)

Then, $r_{\mathsf{R},k}=\Re(r_{k})$ and $r_{\mathsf{I},k}=\Im(r_{k})$, with
$r_{\mathsf{R},k},r_{\mathsf{I},k}\in\\{\pm 1\\}$, for $k=1,\ldots,\tau
N_{R}$. Correspondingly, we have

$\displaystyle r_{k}$
$\displaystyle=r_{\mathsf{R},k}+jr_{\mathsf{I},k}=\textrm{sgn}(\Re(b_{k}))+j\textrm{sgn}(\Im(b_{k}))$
(10)
$\displaystyle=\textrm{sgn}(\Re(\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h})+\Re(n_{k}))+j\textrm{sgn}(\Im(\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h})+\Im(n_{k})),$
(11)

where $\mathbf{a}_{k}^{\mathrm{T}}$ denotes the $k$-th row of $\mathbf{A}$,
$\Re(n_{k})$ and $\Im(n_{k})$ are independent and identically distributed
(i.i.d.) $\mathcal{N}(0,\frac{\sigma^{2}}{2})$ random variables, and

$\displaystyle(\Re(\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h}))^{2}+(\Im(\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h}))^{2}=\mathbf{h}^{\mathrm{H}}\mathbf{a}_{k}\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h},$
(12)

for $k=1,\ldots,\tau N_{R}$. It can be readily shown that

$\displaystyle\Pr(r_{\mathsf{R},k}|\mathbf{h})$
$\displaystyle=\int\limits_{0}^{\infty}\frac{\exp\left\\{-\frac{1}{\sigma^{2}}\left(x_{k}-r_{\mathsf{R},k}\Re(\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h})\right)^{2}\right\\}}{\sqrt{\pi\sigma^{2}}}{\rm
d}x_{k},$ (13) $\displaystyle\Pr(r_{\mathsf{I},k}|\mathbf{h})$
$\displaystyle=\int\limits_{0}^{\infty}\frac{\exp\left\\{-\frac{1}{\sigma^{2}}\left(y_{k}-r_{\mathsf{I},k}\Im(\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h})\right)^{2}\right\\}}{\sqrt{\pi\sigma^{2}}}{\rm
d}y_{k}.$ (14)

Based on the distribution of the noise, we have

$\displaystyle\Pr(r_{k}|\mathbf{h})=\Pr(r_{\mathsf{R},k}|\mathbf{h})\Pr(r_{\mathsf{I},k}|\mathbf{h}),$
(15)

$k=1,\ldots,\tau N_{R}$, and thus

$\displaystyle\Pr(\mathbf{r}|\mathbf{h})=\prod_{k=1}^{\tau
N_{R}}\Pr(r_{k}|\mathbf{h})$ (16)
$\displaystyle=\displaystyle\int\limits_{\mathbb{R}_{+}^{\tau
N_{R}}}\int\limits_{\mathbb{R}_{+}^{\tau N_{R}}}\frac{e^{-\overset{\tau
N_{R}}{\underset{k=1}{\sum}}\frac{\left(x_{k}-r_{\mathsf{R},k}\Re(\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h})\right)^{2}+\left(y_{k}-r_{\mathsf{I},k}\Im(\mathbf{a}_{k}^{\mathrm{T}}\mathbf{h})\right)^{2}}{\sigma^{2}}}}{(\pi\sigma^{2})^{\tau
N_{R}}}{\rm d}\mathbf{x}{\rm d}\mathbf{y},$ (17)

with

$\displaystyle\mathbf{x}=[x_{1}\;\ldots\;x_{\tau
N_{R}}]^{\mathrm{T}},\quad\mathbf{y}=[y_{1}\;\ldots\;y_{\tau
N_{R}}]^{\mathrm{T}}.$ (18)

Define

$\displaystyle\bm{\Lambda}_{\mathsf{R}}=\mathtt{diag}{{(}}\mathbf{r}_{\mathsf{R}}),\quad\bm{\Lambda}_{\mathsf{I}}=\mathtt{diag}{{(}}\mathbf{r}_{\mathsf{I}}),$
(19)

and note that
$\bm{\Lambda}_{\mathsf{R}}\bm{\Lambda}_{\mathsf{R}}=\bm{\Lambda}_{\mathsf{I}}\bm{\Lambda}_{\mathsf{I}}=\mathbf{I}_{\tau
N_{R}}$. Then, (17) is equivalent to

$\displaystyle\Pr(\mathbf{r}|\mathbf{h})=\int_{\mathbb{R}_{+}^{\tau
N_{R}}}\int_{\mathbb{R}_{+}^{\tau
N_{R}}}\frac{\exp\left\\{-\frac{\|\bm{\Lambda}_{\mathsf{R}}\mathbf{x}+j\bm{\Lambda}_{\mathsf{I}}\mathbf{y}-\mathbf{A}\mathbf{h}\|^{2}}{\sigma^{2}}\right\\}}{(\pi\sigma^{2})^{\tau
N_{R}}}{\rm d}\mathbf{x}{\rm d}\mathbf{y}.$ (20)

The covariance matrix of $\mathbf{b}$ in (2) is given by

$\displaystyle\bm{\Omega}_{\rm
b}=\mathbf{A}\bm{\Sigma}\mathbf{A}^{\mathrm{H}}+\sigma^{2}\mathbf{I}_{\tau
N_{R}}\in\mathbb{C}^{(\tau N_{R})\times(\tau N_{R})}.$ (21)

For subsequent derivation, let
$\bm{\Theta}=\left(\bm{\Sigma}^{-1}+\frac{1}{\sigma^{2}}\mathbf{A}^{\mathrm{H}}\mathbf{A}\right)^{-1}$
and let $\bm{\mu}_{\rm
h}=\frac{1}{\sigma^{2}}\bm{\Theta}\mathbf{A}^{\mathrm{H}}\left(\bm{\Lambda}_{\mathsf{R}}\mathbf{x}+j\bm{\Lambda}_{\mathsf{I}}\mathbf{y}\right)$.
Using the matrix inversion identity [34], we have
$\frac{1}{\sigma^{2}}\bm{\Theta}\mathbf{A}^{\mathrm{H}}=\bm{\Sigma}\mathbf{A}^{\mathrm{H}}\bm{\Omega}_{\rm
b}^{-1}$ and thus, $\bm{\mu}_{\rm
h}=\bm{\Sigma}\mathbf{A}^{\mathrm{H}}\bm{\Omega}_{\rm
b}^{-1}\left(\bm{\Lambda}_{\mathsf{R}}\mathbf{x}+j\bm{\Lambda}_{\mathsf{I}}\mathbf{y}\right)$.
Based on (4) and (20), it can be shown that (8) is given by

$\displaystyle\Pr(\mathbf{r})$
$\displaystyle=\displaystyle\int_{\mathbb{C}^{N_{T}N_{R}}}\Pr(\mathbf{r}|\mathbf{h})\left\\{\frac{\exp\left\\{-\mathbf{h}^{\mathrm{H}}\bm{\Sigma}^{-1}\mathbf{h}\right\\}}{\pi^{N_{T}N_{R}}|\bm{\Sigma}|}\right\\}{\rm
d}\mathbf{h}$ (22) $\displaystyle=\displaystyle\int_{\mathbb{R}_{+}^{\tau
N_{R}}}\int_{\mathbb{R}_{+}^{\tau
N_{R}}}\frac{e^{-\frac{1}{\sigma^{2}}\|\bm{\Lambda}_{\mathsf{R}}\mathbf{x}+j\bm{\Lambda}_{\mathsf{I}}\mathbf{y}\|^{2}}e^{\bm{\mu}_{\rm
h}^{\mathrm{H}}\bm{\Theta}^{-1}\bm{\mu}_{\rm h}}}{\pi^{\tau
N_{R}}|\bm{\Omega}_{\rm b}|}{\rm d}\mathbf{x}{\rm d}\mathbf{y},$ (23)

where we have used the fact that $\mathcal{CN}(\bm{\mu}_{\rm h},\bm{\Theta})$
integrates to one over its entire support. Further algebraic manipulations
yield

$\displaystyle\Pr(\mathbf{r})$
$\displaystyle=\displaystyle\int\limits_{\mathbb{R}_{+}^{\tau
N_{R}}}\int\limits_{\mathbb{R}_{+}^{\tau
N_{R}}}\frac{e^{-(\mathbf{x}^{\mathrm{T}}\bm{\Lambda}_{\mathsf{R}}\mathbf{D}_{\mathsf{R}}\bm{\Lambda}_{\mathsf{R}}\mathbf{x}+\mathbf{y}^{\mathrm{T}}\bm{\Lambda}_{\mathsf{I}}\mathbf{D}_{\mathsf{R}}\bm{\Lambda}_{\mathsf{I}}\mathbf{y}-2\mathbf{x}^{\mathrm{T}}\bm{\Lambda}_{\mathsf{R}}\mathbf{D}_{\mathsf{I}}\bm{\Lambda}_{\mathsf{I}}\mathbf{y})}}{\pi^{\tau
N_{R}}|\bm{\Omega}_{\rm b}|}{\rm d}\mathbf{x}{\rm d}\mathbf{y},$ (24)

with

$\displaystyle\mathbf{D}_{\mathsf{R}}$
$\displaystyle=\Re\left(\bm{\Omega}_{\rm
b}^{-1}\right)=\mathbf{D}_{\mathsf{R}}^{\mathrm{T}}\in\mathbb{R}^{(\tau
N_{R})\times(\tau N_{R})},$ (25) $\displaystyle\mathbf{D}_{\mathsf{I}}$
$\displaystyle=\Im\left(\bm{\Omega}_{\rm
b}^{-1}\right)=-\mathbf{D}_{\mathsf{I}}^{\mathrm{T}}\in\mathbb{R}^{(\tau
N_{R})\times(\tau N_{R})}.$ (26)

At this point, (24) can be compactly expressed as

$\displaystyle\Pr(\mathbf{r})=\frac{1}{\pi^{\tau N_{R}}|\bm{\Omega}_{\rm
b}|}\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}\exp\left\\{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}\right\\}{\rm
d}\mathbf{z},$ (27)

with

$\displaystyle\mathbf{C}$
$\displaystyle=\begin{bmatrix}\bm{\Lambda}_{\mathsf{R}}\mathbf{D}_{\mathsf{R}}\bm{\Lambda}_{\mathsf{R}}&\bm{\Lambda}_{\mathsf{R}}\mathbf{D}_{\mathsf{I}}^{\mathrm{T}}\bm{\Lambda}_{\mathsf{I}}\\\
\bm{\Lambda}_{\mathsf{I}}\mathbf{D}_{\mathsf{I}}\bm{\Lambda}_{\mathsf{R}}&\bm{\Lambda}_{\mathsf{I}}\mathbf{D}_{\mathsf{R}}\bm{\Lambda}_{\mathsf{I}}\end{bmatrix}\in\mathbb{R}^{(2\tau
N_{R})\times(2\tau N_{R})},$ (28) $\displaystyle\mathbf{z}$
$\displaystyle=\left[\mathbf{x}^{\mathrm{T}}\;\;\mathbf{y}^{\mathrm{T}}\right]^{\mathrm{T}}\in\mathbb{R}^{2\tau
N_{R}}.$ (29)

Similarly, the numerator of (7) is given by

$\displaystyle\displaystyle\int_{\mathbb{C}^{N_{T}N_{R}}}\mathbf{h}\Pr(\mathbf{r}|\mathbf{h})p(\mathbf{h}){\rm
d}\mathbf{h}$
$\displaystyle=\frac{\bm{\Sigma}\mathbf{A}^{\mathrm{H}}\bm{\Omega}_{\rm
b}^{-1}\left[\bm{\Lambda}_{\mathsf{R}}\;\;j\bm{\Lambda}_{\mathsf{I}}\right]}{\pi^{\tau
N_{R}}|\bm{\Omega}_{\rm b}|}$ $\displaystyle\phantom{=}\
\cdot\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}\mathbf{z}e^{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}}{\rm
d}\mathbf{z}.$ (30)

Based on (7), (27), and (30), we have

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\bm{\Sigma}\mathbf{A}^{\mathrm{H}}\bm{\Omega}_{\rm
b}^{-1}\left[\bm{\Lambda}_{\mathsf{R}}\;\;j\bm{\Lambda}_{\mathsf{I}}\right]\left(\frac{\displaystyle\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}\mathbf{z}e^{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}}{\rm
d}\mathbf{z}}{\displaystyle\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}e^{-{\mathbf{z}}^{\mathrm{T}}\mathbf{C}{\mathbf{z}}}{\rm
d}{\mathbf{z}}}\right).$ (31)

Note that the quotient of the integrals in (31) denotes the mean of a
truncated MVN distribution with $\mathbb{R}_{+}^{2\tau N_{R}}$ as its support.
Moreover, since

$\displaystyle\frac{{\rm d}}{{\rm
d}\mathbf{z}}\exp\left\\{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}\right\\}=-2\mathbf{Cz}\exp\left\\{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}\right\\},$
(32)

we have

$\displaystyle\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}\mathbf{z}\exp\left\\{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}\right\\}{\rm
d}\mathbf{z}=\frac{1}{2}\mathbf{C}^{-1}\bm{\alpha},$ (33)

where the $k$-th element of $\bm{\alpha}\in\mathbb{R}^{2\tau N_{R}}$ is given
by

$\displaystyle\alpha_{k}=\int_{\mathbb{R}_{+}^{2\tau
N_{R}-1}}\exp\left\\{-\mathbf{z}_{-k}^{\mathrm{T}}\mathbf{C}_{-k}\mathbf{z}_{-k}\right\\}{\rm
d}\mathbf{z}_{-k},$ (34)

with $\mathbf{z}_{-k}=[z_{1}\;\ldots\;z_{k-1}\;z_{k+1}\;\ldots\;z_{2\tau
N_{R}}]^{\mathrm{T}}\in\mathbb{R}^{2\tau N_{R}-1}$ and where
$\mathbf{C}_{-k}\in\mathbb{R}^{(2\tau N_{R}-1)\times(2\tau N_{R}-1)}$ is
obtained by deleting the $k$-th row and column of $\mathbf{C}$. Consequently,
(31) is equivalent to

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\frac{\bm{\Sigma}\mathbf{A}^{\mathrm{H}}\bm{\Omega}_{\rm
b}^{-1}\left[\bm{\Lambda}_{\mathsf{R}}\;\;j\bm{\Lambda}_{\mathsf{I}}\right]\mathbf{C}^{-1}\bm{\alpha}}{\displaystyle
2\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}\exp\left\\{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}\right\\}{\rm
d}\mathbf{z}}.$ (35)

######  Remark 1.

In the expression for $\hat{\mathbf{h}}_{\rm MMSE}$ used in [19], the
quantized outputs (i.e., (9) or (19)) affect the integration limits. In
contrast, in (31) and (35), they appear outside the integrals and in the
integrands through $\mathbf{C}$ and $\mathbf{C}_{-k}$, for $k=1,\ldots,2\tau
N_{R}$, which facilitates linking $\hat{\mathbf{h}}_{\rm MMSE}$ to the
Gaussian orthant probability, as will be shown later.

######  Remark 2.

The diagonal elements of $\mathbf{C}$ are not affected by the quantized
output. Each off-diagonal element of $\mathbf{C}$ is affected by exactly two
different elements from
$[\mathbf{r}_{\mathsf{R}}^{\mathrm{T}}\;\;\mathbf{r}_{\mathsf{I}}^{\mathrm{T}}]^{\mathrm{T}}$.
The same comments also apply to $\mathbf{C}_{-k}$, for $k=1,\ldots,2\tau
N_{R}$. These observations are instrumental for deriving simplified but
equivalent expressions for $\hat{\mathbf{h}}_{\rm MMSE}$ under specific
assumptions on the channel correlation or pilot symbols (see Theorems 3 and 4
in Section VI).

In (35), the expression for $\hat{\mathbf{h}}_{\rm MMSE}$ is directly related
to the orthant probability of the real-valued MVN distribution defined below
[26, 32, 33]:

$\displaystyle\mathcal{P}(\bm{\Psi})=\frac{1}{(2\pi)^{\frac{L}{2}}\left|\bm{\Psi}\right|^{\frac{1}{2}}}\displaystyle\int_{\mathbb{R}_{+}^{L}}\exp\left\\{-\frac{1}{2}\mathbf{u}^{\mathrm{T}}\bm{\Psi}^{-1}\mathbf{u}\right\\}{\rm
d}\mathbf{u},$ (36)

where $\bm{\Psi}$ denotes the $L\times L$ covariance matrix in
$\mathcal{N}(\mathbf{0},\bm{\Psi})$. Using the notation in (36), (34) is
equivalently expressed as

$\displaystyle\alpha_{k}=\left(2\pi\right)^{\frac{2\tau
N_{R}-1}{2}}\left|\frac{1}{2}\mathbf{C}_{-k}^{-1}\right|^{\frac{1}{2}}\mathcal{P}\left(\frac{1}{2}\mathbf{C}_{-k}^{-1}\right),$
(37)

and the integral in the denominator of (35) is given by

$\displaystyle\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}e^{-\mathbf{z}^{\mathrm{T}}\mathbf{C}\mathbf{z}}{\rm
d}\mathbf{z}=\left(2\pi\right)^{\tau
N_{R}}\left|\frac{1}{2}\mathbf{C}^{-1}\right|^{\frac{1}{2}}\mathcal{P}\left(\frac{1}{2}\mathbf{C}^{-1}\right).$
(38)

It can also be shown that, for $k=1,\ldots,2\tau N_{R}$,

$\displaystyle|\mathbf{C}_{-k}^{-1}|/|\mathbf{C}^{-1}|=|\mathbf{C}|/|\mathbf{C}_{-k}|=1/[\mathbf{C}^{-1}]_{kk}.$
(39)

Based on (37)–(39), (35) is shown to be equivalent to

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\frac{\bm{\Sigma}\mathbf{A}^{\mathrm{H}}\bm{\Omega}_{\rm
b}^{-1}\left[\bm{\Lambda}_{\mathsf{R}}\;\;j\bm{\Lambda}_{\mathsf{I}}\right]\mathbf{C}^{-1}\left(\widetilde{\mathtt{diag}{{}}}(\mathbf{C}^{-1})\right)^{-\frac{1}{2}}\mathbf{g}}{2\sqrt{\pi}\mathcal{P}\left(\frac{1}{2}\mathbf{C}^{-1}\right)},$
(40)

where the $k$-th element of $\mathbf{g}$ is given by

$\displaystyle
g_{k}=\mathcal{P}\left(\frac{1}{2}\mathbf{C}_{-k}^{-1}\right)\quad(k=1,\ldots,2\tau
N_{R}).$ (41)

######  Remark 3.

From (3), (28), (31), (35), and (40)–(41), the channel covariance matrix
$\bm{\Sigma}$ as well as the pilot matrix $\mathbf{S}^{\mathrm{T}}$ determines
the computational complexity of $\hat{\mathbf{h}}_{\rm MMSE}$ through
$\bm{\Omega}_{\rm b}$ that is embedded in $\mathbf{C}$ (see (28)). The number
of transmit antennas $N_{T}$ is not involved in computing the orthant
probabilities.

Since the Gaussian orthant probability plays a central role in the expressions
for $\hat{\mathbf{h}}_{\rm MMSE}$ in (40)–(41), we will study the orthant
probability more deeply in the next section.

##  IV Orthant Probability of a Real-Valued MVN Distribution

Referring back to (36), we define

$\displaystyle\bm{\Psi}_{\rm
std}=\left[\widetilde{\mathtt{diag}{{}}}(\bm{\Psi})\right]^{-\frac{1}{2}}\bm{\Psi}\left[\widetilde{\mathtt{diag}{{}}}(\bm{\Psi})\right]^{-\frac{1}{2}}.$
(42)

Here, $\bm{\Psi}_{\rm std}$ denotes the covariance matrix of the standardized
version of $L$ zero-mean random variables whose original covariance matrix is
$\bm{\Psi}$ in (36) [35]. The diagonal elements of $\bm{\Psi}_{\rm std}$ are
all equal to one, while the off-diagonal elements represent the correlation
coefficients between the random variables defined by $\bm{\Psi}$. Let

$\displaystyle\left[\bm{\Psi}_{\rm
std}\right]_{ik}=\psi_{ik}\in(-1,1)\quad(i\neq k,\;i,k=1,\ldots,L).$ (43)

Through a change of variables, one can readily show that [27]

$\displaystyle\mathcal{P}(\bm{\Psi}_{\rm std})=\mathcal{P}(\bm{\Psi}).$ (44)

In [26, 27, 30], it was shown that, for $L=2$,

$\displaystyle\mathcal{P}(\bm{\Psi})=\frac{1}{4}+\frac{\arcsin(\psi_{12})}{2\pi},$
(45)

and, for $L=3$,

$\displaystyle\mathcal{P}(\bm{\Psi})=\frac{1}{8}+\frac{\sum_{i=1}^{2}\sum_{k=i+1}^{3}\arcsin(\psi_{ik})}{4\pi}.$
(46)

Furthermore, for $L=4$, an explicit expression for $\mathcal{P}(\bm{\Psi})$ is
provided in Appendix A based on [27], which requires finite integration with
one variable. There is no general closed-form expression for
$\mathcal{P}(\bm{\Psi})$ when $L\geq 4$ [27, 29, 28]. In [27, Eq. (6)], an
expression for the orthant probability was derived based on characteristic
functions, where the number of terms to be evaluated grows exponentially with
$L$ and the complexity of evaluating some individual integral terms also grows
with $L$. Alternatively, one can use Monte Carlo methods to compute (36),
e.g., the algorithm in [31] that was used in [19]. However, Monte Carlo
methods also incur a high computational load when $L$ becomes large, which, in
turn, generally implies a high computational complexity for determining
$\hat{\mathbf{h}}_{\rm MMSE}$ with large values of $N_{R}$ and/or $\tau$.

######  Lemma 1.

Let $\mathbf{q}=[q_{1}\;\;q_{2}]^{\mathrm{T}}$, with $q_{1},q_{2}\in\\{\pm
1\\}$. Let $p(\mathbf{u})$, the pdf of the real-valued random vector
$\mathbf{u}=[u_{1}\;\;u_{2}]^{\mathrm{T}}$, be given by
$\mathcal{N}(\mathbf{0},\bm{\Phi})$, with $[\bm{\Phi}]_{ii}=\sigma_{ui}^{2}$,
$i=1,2$, and
$[\bm{\Phi}]_{12}=[\bm{\Phi}]_{21}=q_{1}q_{2}\phi_{12}\sigma_{u1}\sigma_{u2}$,
$\phi_{12}\in(-1,1)$. Let
$\bm{\Lambda}_{\mathbf{q}}=\mathtt{diag}{{(}}\mathbf{q})$. Then, the vector

$\displaystyle\bm{\Lambda}_{\mathbf{q}}\displaystyle\int_{\mathbb{R}_{+}^{2}}\mathbf{u}p(\mathbf{u}){\rm
d}\mathbf{u}\Bigg{/}\displaystyle\int_{\mathbb{R}_{+}^{2}}p(\mathbf{u}){\rm
d}\mathbf{u}$ (47)

is linear in $\mathbf{q}$.

###### Proof:

Based on (44) and (45),

$\displaystyle\displaystyle\int_{\mathbb{R}_{+}^{2}}p(\mathbf{u}){\rm
d}\mathbf{u}=\mathcal{P}(\bm{\Phi})=\frac{1}{4}+\frac{q_{1}q_{2}\arcsin(\phi_{12})}{2\pi}.$
(48)

Based on the reduction in (32)–(34), we obtain

$\displaystyle\displaystyle\int_{\mathbb{R}_{+}^{2}}\mathbf{u}p(\mathbf{u}){\rm
d}\mathbf{u}=\frac{1+q_{1}q_{2}\phi_{12}}{2\sqrt{2\pi}}\left[\sigma_{u1}\;\;\sigma_{u2}\right]^{\mathrm{T}}.$
(49)

Inserting (48)–(49) into (47) and capitalizing on the fact that
$q_{1},q_{2}\in\\{\pm 1\\}$, after some algebraic manipulations we obtain

$\displaystyle\left.\bm{\Lambda}_{\mathbf{q}}\displaystyle\int_{\mathbb{R}_{+}^{2}}\mathbf{u}p(\mathbf{u}){\rm
d}\mathbf{u}\middle/\displaystyle\int_{\mathbb{R}_{+}^{2}}p(\mathbf{u}){\rm
d}\mathbf{u}\right.$
$\displaystyle=\frac{\begin{bmatrix}\sigma_{u1}\left(1-\frac{\phi_{12}\arcsin(\phi_{12})}{\pi/2}\right)&\sigma_{u1}\left(\phi_{12}-\frac{\arcsin(\phi_{12})}{\pi/2}\right)\\\
\sigma_{u2}\left(\phi_{12}-\frac{\arcsin(\phi_{12})}{\pi/2}\right)&\sigma_{u2}\left(1-\frac{\phi_{12}\arcsin(\phi_{12})}{\pi/2}\right)\end{bmatrix}\mathbf{q}}{\sqrt{\frac{\pi}{2}}\left(1-\left(\frac{\arcsin(\phi_{12})}{\pi/2}\right)^{2}\right)}.$
(50)

∎

######  Remark 4.

Relation between Lemma 1 and the optimality of the BLMMSE estimator: Lemma 1
holds based on the use of one-bit symmetric quantization, i.e.,
$q_{1},q_{2}\in\\{\pm 1\\}$. This linearity does not hold generally when
$\mathbf{u}$ and $\mathbf{q}$ have three or more elements, as shown in
Appendix B. These results are pivotal in determining the condition for the
BLMMSE estimator to be MSE-optimal, which will be presented in the next
section.

##  V Condition for the Optimality of the BLMMSE Estimator

The equivalence between the BLMMSE and MMSE channel estimators for one-bit
quantized SISO systems with a single pilot symbol or one-bit quantized
noiseless SISO systems with multiple specifically chosen pilots was shown in
[19]. Here, our objective is to establish a general condition for their
equivalence in one-bit quantized MIMO systems.

Defining

$\displaystyle\mathbf{D}_{\Omega}=\widetilde{\mathtt{diag}{{}}}(\bm{\Omega}_{\rm
b})$ (51)

from (21), the BLMMSE channel estimate for the system setup in Section II-A
can be expressed as [7]

$\displaystyle\hat{\mathbf{h}}_{\rm BLM}$
$\displaystyle=\frac{\sqrt{\pi}}{2}\bm{\Sigma}\mathbf{A}^{\mathrm{H}}\mathbf{D}_{\Omega}^{-\frac{1}{2}}\left(\arcsin\left(\mathbf{D}_{\Omega}^{-\frac{1}{2}}\Re\left(\bm{\Omega}_{\rm
b}\right)\mathbf{D}_{\Omega}^{-\frac{1}{2}}\right)\right.$
$\displaystyle\phantom{=}\
+\left.j\arcsin\left(\mathbf{D}_{\Omega}^{-\frac{1}{2}}\Im\left(\bm{\Omega}_{\rm
b}\right)\mathbf{D}_{\Omega}^{-\frac{1}{2}}\right)\right)^{-1}\mathbf{r}.$
(52)

######  Theorem 1.

Necessary and sufficient condition for the BLMMSE channel estimator to be MSE-
optimal: Assuming Rayleigh flat fading, spatially and temporally white
Gaussian noise, and one-bit ADCs, the BLMMSE channel estimator is MSE-optimal
if and only if the off-diagonal elements of $\mathbf{C}$ in (28) satisfy
$[\mathbf{C}]_{il}\neq 0$ for at most one value of $l$ other than $l=i$, i.e.,
an element $z_{i}$ of the random vector $\mathbf{z}$ defined in (29) is
correlated with at most one other element $z_{l}$, for $l,i=1,\ldots,2\tau
N_{R}$ and $l\neq i$.

###### Proof:

First, we show that when the condition on $\mathbf{C}$ in (28) is satisfied,
the BLMMSE channel estimator is optimal. Note that
$\frac{1}{2}\mathbf{C}^{-1}$ represents the covariance matrix of $\mathbf{z}$
in (29). Based on (28), a suitable permutation matrix $\bm{\Gamma}$ is
introduced to make
$\bm{\Gamma}\mathbf{C}\bm{\Gamma}^{\mathrm{T}}=\mathbf{C}_{\rm perm}$ block-
diagonal with $2\times 2$ blocks, with
$\bm{\Gamma}\bm{\Gamma}^{\mathrm{T}}=\bm{\Gamma}^{\mathrm{T}}\bm{\Gamma}=\mathbf{I}_{2\tau
N_{R}}$. These $2\times 2$ sub-matrices can be diagonal. Let
$\mathbf{M}=\frac{1}{2}\mathbf{C}_{\rm perm}^{-1}$ and note that $\mathbf{M}$
is also block-diagonal with $2\times 2$ blocks. For $i=1,\ldots,\tau N_{R}$,
the $i$-th diagonal block of $\mathbf{M}$, denoted as $\mathbf{M}_{i}$,
represents the cross-covariance between a pair of correlated random variables
if its off-diagonal elements are non-zero. Furthermore,
$\mathcal{P}(\mathbf{M})=\prod_{i=1}^{\tau N_{R}}\mathcal{P}(\mathbf{M}_{i})$.
The pdf of $\mathbf{t}=\bm{\Gamma}\mathbf{z}$, denoted by $p(\mathbf{t})$, is
given by $\mathcal{N}(\mathbf{0},\mathbf{M})$. After applying this change of
variable to (31), we obtain

$\displaystyle\hat{\mathbf{h}}_{\rm MMSE}$
$\displaystyle=\bm{\Sigma}\mathbf{A}^{\mathrm{H}}\bm{\Omega}_{\rm
b}^{-1}\left[\bm{\Lambda}_{\mathsf{R}}\;\;j\bm{\Lambda}_{\mathsf{I}}\right]\bm{\Gamma}^{\mathrm{T}}$
$\displaystyle\phantom{=}\ \cdot\frac{\displaystyle\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}\frac{\mathbf{t}\exp\left\\{-\frac{\mathbf{t}^{\mathrm{T}}\mathbf{M}^{-1}\mathbf{t}}{2}\right\\}}{(2\pi)^{\tau
N_{R}}\left|\mathbf{M}\right|^{\frac{1}{2}}}{\rm
d}\mathbf{t}}{\displaystyle\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}\frac{\exp\left\\{-\frac{\mathbf{t}^{\mathrm{T}}\mathbf{M}^{-1}\mathbf{t}}{2}\right\\}}{(2\pi)^{\tau
N_{R}}\left|\mathbf{M}\right|^{\frac{1}{2}}}{\rm d}\mathbf{t}}$ (53)
$\displaystyle=\bm{\Sigma}\mathbf{A}^{\mathrm{H}}\bm{\Omega}_{\rm
b}^{-1}\left[\mathbf{I}_{\tau N_{R}}\;\;\mathbf{I}_{\tau
N_{R}}\right]\bm{\Gamma}^{\mathrm{T}}$ $\displaystyle\phantom{=}\
\cdot\bm{\Gamma}\;\mathtt{diag}{{\left([\mathbf{r}_{\mathsf{R}}^{\mathrm{T}}\;\;j\mathbf{r}_{\mathsf{I}}^{\mathrm{T}}]^{\mathrm{T}}\right)}}\bm{\Gamma}^{\mathrm{T}}\frac{\displaystyle\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}\mathbf{t}p(\mathbf{t}){\rm d}\mathbf{t}}{\mathcal{P}(\mathbf{M})}.$
(54)

The vector $\mathbf{t}$ defined above can be written as the concatenation of
$\tau N_{R}$ $2$-dimensional vectors, i.e.,

$\displaystyle\mathbf{t}=[\widetilde{\mathbf{t}}_{i}]_{i=1}^{\tau
N_{R}}=[\widetilde{\mathbf{t}}_{1}^{\mathrm{T}}\;\ldots\;\widetilde{\mathbf{t}}_{\tau
N_{R}}^{\mathrm{T}}]^{\mathrm{T}},$ (55)

where $\widetilde{\mathbf{t}}_{i}$ consists of $t_{2i-1}$ and $t_{2i}$ from
$\mathbf{t}$, $i=1,\ldots,\tau N_{R}$. Using the notation in (55), we have

$\displaystyle\frac{\displaystyle\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}\mathbf{t}p(\mathbf{t}){\rm
d}\mathbf{t}}{\mathcal{P}(\mathbf{M})}=\left[\frac{\left(\displaystyle\int_{\mathbb{R}_{+}^{2}}\widetilde{\mathbf{t}}_{k}p(\widetilde{\mathbf{t}}_{k}){\rm
d}\widetilde{\mathbf{t}}_{k}\right)\prod_{\begin{subarray}{c}i=1\\\ i\neq
k\end{subarray}}^{\tau
N_{R}}\mathcal{P}(\mathbf{M}_{i})}{\mathcal{P}(\mathbf{M})}\right]_{k=1}^{\tau
N_{R}}$ (56)
$\displaystyle=\left[\frac{\displaystyle\int_{\mathbb{R}_{+}^{2}}\widetilde{\mathbf{t}}_{1}^{\mathrm{T}}p(\widetilde{\mathbf{t}}_{1}){\rm
d}\widetilde{\mathbf{t}}_{1}}{\mathcal{P}(\mathbf{M}_{1})}\;\ldots\;\frac{\displaystyle\int_{\mathbb{R}_{+}^{2}}\widetilde{\mathbf{t}}_{\tau
N_{R}}^{\mathrm{T}}p(\widetilde{\mathbf{t}}_{\tau N_{R}}){\rm
d}\widetilde{\mathbf{t}}_{\tau N_{R}}}{\mathcal{P}(\mathbf{M}_{\tau
N_{R}})}\right]^{\mathrm{T}}.$ (57)

Since each component of (57) is the ratio of two integrals as in (47), we
conclude based on Lemma 1 and (50) that

$\displaystyle\bm{\Gamma}\;\mathtt{diag}{{\left([\mathbf{r}_{\mathsf{R}}^{\mathrm{T}}\;\;j\mathbf{r}_{\mathsf{I}}^{\mathrm{T}}]^{\mathrm{T}}\right)}}\bm{\Gamma}^{\mathrm{T}}\displaystyle\int_{\mathbb{R}_{+}^{2\tau
N_{R}}}\frac{\mathbf{t}p(\mathbf{t})}{\mathcal{P}(\mathbf{M})}{\rm
d}\mathbf{t}$ (58)

is linear in $\mathbf{r}_{\mathsf{R}}$ and $\mathbf{r}_{\mathsf{I}}$.
Therefore, $\hat{\mathbf{h}}_{\rm MMSE}$ in (54) is linear in $\mathbf{r}$.
Since the BLMMSE estimator has the smallest MSE [7, 24] among all linear
estimators, the BLMMSE estimator must be MSE-optimal.

Conversely, when the BLMMSE estimator is optimal, i.e., when
$\hat{\mathbf{h}}_{\rm MMSE}$ is linear in $\mathbf{r}$, based on Remark 4 and
Appendix B, there cannot be any full $m\times m$ (with $m\geq 3$) sub-matrix
in $\frac{1}{2}\mathbf{C}^{-1}$ with only non-zero elements. The dimension of
any full sub-matrix of $\frac{1}{2}\mathbf{C}^{-1}$ must be at most $2\times
2$. The same can be concluded for $\mathbf{C}$. ∎

The proof of Theorem 1 hinges on Lemma 1 as well as Remark 4 and Appendix B.
Guided by this framework, below we present specific cases where
$\hat{\mathbf{h}}_{\rm MMSE}$ is directly computed to be equivalent to
$\hat{\mathbf{h}}_{\rm BLM}$.

###  V-A Spatially White Channel and Unitary Pilot Matrix

Here, we assume $\bm{\Sigma}=\mathbf{I}_{N_{T}N_{R}}$ and a unitary pilot
matrix with $\tau=N_{T}$. Thus,
$\mathbf{S}\mathbf{S}^{\mathrm{H}}=\eta\mathbf{I}$, with $\eta>0$, and

$\displaystyle\mathbf{D}_{\mathsf{R}}=\frac{1}{\eta+\sigma^{2}}\mathbf{I}_{\tau
N_{R}},\quad\mathbf{D}_{\mathsf{I}}=\mathbf{0}.$ (59)

For this case,
$\mathbf{C}=\mathtt{DIAG}{{(}}\mathbf{D}_{\mathsf{R}},\mathbf{D}_{\mathsf{R}})$
is diagonal, thereby satisfying the condition of Theorem 1.

Due to the fact that $\mathbf{D}_{\mathsf{I}}=\mathbf{0}$, from (29) and (31)
we have

$\displaystyle\hat{\mathbf{h}}_{\rm MMSE}$
$\displaystyle=\bm{\Sigma}\mathbf{A}^{\mathrm{H}}\bm{\Omega}_{\rm
b}^{-1}\left(\frac{\bm{\Lambda}_{\mathsf{R}}\displaystyle\int_{\mathbb{R}_{+}^{\tau
N_{R}}}\mathbf{x}e^{-\mathbf{x}^{\mathrm{T}}\bm{\Lambda}_{\mathsf{R}}\mathbf{D}_{\mathsf{R}}\bm{\Lambda}_{\mathsf{R}}\mathbf{x}}{\rm
d}\mathbf{x}}{\displaystyle\int_{\mathbb{R}_{+}^{\tau
N_{R}}}e^{-\mathbf{x}^{\mathrm{T}}\bm{\Lambda}_{\mathsf{R}}\mathbf{D}_{\mathsf{R}}\bm{\Lambda}_{\mathsf{R}}\mathbf{x}}{\rm
d}\mathbf{x}}\right.$ $\displaystyle\phantom{=}\
+\left.\frac{j\bm{\Lambda}_{\mathsf{I}}\displaystyle\int_{\mathbb{R}_{+}^{\tau
N_{R}}}\mathbf{y}e^{-\mathbf{y}^{\mathrm{T}}\bm{\Lambda}_{\mathsf{I}}\mathbf{D}_{\mathsf{R}}\bm{\Lambda}_{\mathsf{I}}\mathbf{y}}{\rm
d}\mathbf{y}}{\displaystyle\int_{\mathbb{R}_{+}^{\tau
N_{R}}}e^{-\mathbf{y}^{\mathrm{T}}\bm{\Lambda}_{\mathsf{I}}\mathbf{D}_{\mathsf{R}}\bm{\Lambda}_{\mathsf{I}}\mathbf{y}}{\rm
d}\mathbf{y}}\right).$ (60)

Further calculations based on (59)–(60) and (52) lead to

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\frac{\left(\mathbf{S}^{\mathrm{H}}\otimes\mathbf{I}_{N_{R}}\right)\mathbf{r}}{\sqrt{\pi(\eta+\sigma^{2})}}=\hat{\mathbf{h}}_{\rm
BLM}.$ (61)

######  Remark 5.

Based on (61), for a spatially white Rayleigh fading channel with a unitary
pilot matrix, the low-complexity BLMMSE estimator is MSE-optimal and thus its
use is justified for low spatial correlation. A similar conclusion was reached
in [19], but only for the SIMO case with a single pilot symbol.

The above result can be generalized even further, as shown next.

###  V-B Transmit-Only Channel Correlation

Assume $\tau=N_{T}$ and let $\mathbf{H}=\mathbf{H}_{\rm w}\bm{\Sigma}_{\rm
TX}^{\frac{1}{2}}$, where the elements of $\mathbf{H}_{\rm w}$ are i.i.d.
$\mathcal{CN}(0,1)$ random variables. Then,
$\mathbf{h}=\textrm{vec}(\mathbf{H})\sim\mathcal{CN}(\mathbf{0},\bm{\Sigma})$,
with $\bm{\Sigma}=\bm{\Sigma}_{\rm TX}\otimes\mathbf{I}_{N_{R}}$. Denote the
eigenvalue decomposition of $\bm{\Sigma}_{\rm TX}$ as

$\displaystyle\bm{\Sigma}_{\rm
TX}=\mathbf{U}\bm{\Xi}\mathbf{U}^{\mathrm{H}},\quad\bm{\Xi}=\mathtt{diag}{{}}(\left[\xi_{i}\right]_{i=1}^{N_{T}}),$
(62)

with $\xi_{i}>0$, $\forall i$. Setting
$\mathbf{S}=\sqrt{\eta}\mathbf{U}^{\mathrm{H}}$, it can be shown that
$\mathbf{D}_{\mathsf{R}}=\mathtt{diag}{{\left(\left[\frac{1}{\eta\xi_{i}+\sigma^{2}}\right]_{i=1}^{N_{T}}\right)}}\otimes\mathbf{I}_{N_{R}}$
and $\mathbf{D}_{\mathsf{I}}=\mathbf{0}$. Moreover,
$\mathbf{C}=\mathtt{DIAG}{{(}}\mathbf{D}_{\mathsf{R}},\mathbf{D}_{\mathsf{R}})$
is diagonal and satisfies the condition in Theorem 1. Since
$\mathbf{D}_{\mathsf{I}}=\mathbf{0}$, (60) holds again and the integrals
therein are over scalar variables. Finally, it is straightforward to show that

$\displaystyle\hat{\mathbf{h}}_{\rm MMSE}$
$\displaystyle=\left(\left(\mathbf{U}\;\mathtt{diag}{{}}\left(\left[\frac{\xi_{i}\sqrt{\eta}}{\sqrt{\eta\xi_{i}+\sigma^{2}}}\right]_{i=1}^{N_{T}}\right)\right)\otimes\mathbf{I}_{N_{R}}\right)\frac{\mathbf{r}}{\sqrt{\pi}}$
(63) $\displaystyle=\hat{\mathbf{h}}_{\rm BLM}.$ (64)

Clearly, (61) is a special case of (64) with $\bm{\Sigma}_{\rm
TX}=\mathbf{I}_{N_{T}}$. If $\bm{\Sigma}_{\rm TX}$ is circulant, then a
Fourier matrix can be used as pilot matrix.

######  Remark 6.

Assume a standardized $\bm{\Sigma}_{\rm TX}$ as in (42). Based on (64) and [7,
17], it can be shown that, when the pilot symbol power is sufficiently high
or, equivalently, when the noise power is sufficiently low, the MSE per
antenna for both estimators is given by $1-\frac{2}{\pi}$, i.e.,

$\displaystyle\lim_{\sigma^{2}\rightarrow
0}\frac{\mathsf{E}(\|\hat{\mathbf{h}}_{\rm
MMSE}-\mathbf{h}\|^{2})}{N_{T}N_{R}}$
$\displaystyle=\lim_{\sigma^{2}\rightarrow
0}\frac{\mathsf{E}(\|\hat{\mathbf{h}}_{\rm
BLM}-\mathbf{h}\|^{2})}{N_{T}N_{R}}$ (65) $\displaystyle=1-\frac{2}{\pi}.$
(66)

###  V-C Spatially White Channel with $\tau=2$ Real-Valued Pilots

Here, $\bm{\Sigma}=\mathbf{I}_{N_{T}N_{R}}$ and $\tau=2$. Let
$\mathbf{s}_{1},\mathbf{s}_{2}$ be two real-valued $N_{T}$-dimensional pilot
vectors in (1) and
$\mathbf{S}^{\mathrm{T}}=[\mathbf{s}_{1}\;\;\mathbf{s}_{2}]$. Clearly,
$\mathbf{D}_{\mathsf{R}}=\left(\mathbf{S}\mathbf{S}^{\mathrm{T}}+\sigma^{2}\mathbf{I}_{2}\right)^{-1}\otimes\mathbf{I}_{N_{R}}$.
Moreover, $\mathbf{D}_{\mathsf{I}}=\mathbf{0}$ and (60) holds. Defining

$\displaystyle\bm{\Delta}_{\rm R}$
$\displaystyle=\mathtt{diag}{{\left(\left[r_{\mathsf{R},i}r_{\mathsf{R},N_{R}+i}\right]_{i=1}^{N_{R}}\right)}},$
(67) $\displaystyle\bm{\Delta}_{\rm I}$
$\displaystyle=\mathtt{diag}{{\left(\left[r_{\mathsf{I},i}r_{\mathsf{I},N_{R}+i}\right]_{i=1}^{N_{R}}\right)}},$
(68) $\displaystyle\mathbf{C}_{{\rm sub,R}}$
$\displaystyle=\begin{bmatrix}(\|\mathbf{s}_{2}\|^{2}+\sigma^{2})\mathbf{I}_{N_{R}}&-\mathbf{s}_{1}^{\mathrm{T}}\mathbf{s}_{2}\bm{\Delta}_{\rm
R}\\\ -\mathbf{s}_{1}^{\mathrm{T}}\mathbf{s}_{2}\bm{\Delta}_{\rm
R}&(\|\mathbf{s}_{1}\|^{2}+\sigma^{2})\mathbf{I}_{N_{R}}\end{bmatrix},$ (69)

and letting $\mathbf{C}_{{\rm sub,I}}$ be identical to $\mathbf{C}_{{\rm
sub,R}}$ except that all $\bm{\Delta}_{\rm R}$ terms are replaced by
$\bm{\Delta}_{\rm I}$, we obtain

$\displaystyle\mathbf{C}=\frac{\mathtt{DIAG}{{(}}\mathbf{C}_{{\rm
sub,R}},\mathbf{C}_{{\rm
sub,I}})}{\left|\mathbf{S}\mathbf{S}^{\mathrm{T}}+\sigma^{2}\mathbf{I}_{2}\right|}.$
(70)

Since $x_{i}$ and $x_{i+N_{R}}$, as well as $y_{i}$ and $y_{i+N_{R}}$,
$i=1,\ldots,N_{R}$, are pair-wise correlated, Theorem 1 holds. Define

$\displaystyle\mathbf{s}_{1,{\rm
new}}=\frac{\mathbf{s}_{1}}{\sqrt{\|\mathbf{s}_{1}\|^{2}+\sigma^{2}}},\quad\mathbf{s}_{2,{\rm
new}}=\frac{\mathbf{s}_{2}}{\sqrt{\|\mathbf{s}_{2}\|^{2}+\sigma^{2}}},$ (71)

$\gamma_{s}=\mathbf{s}_{1,{\rm new}}^{\mathrm{T}}\mathbf{s}_{2,{\rm new}}$,
and $\mathbf{S}_{\rm w1}^{\mathrm{T}}=\left[\mathbf{s}_{1,{\rm
new}}\;\;\mathbf{s}_{2,{\rm new}}\right]$. Using the method in the proof of
Theorem 1, we obtain

$\displaystyle\hat{\mathbf{h}}_{\rm MMSE}=\hat{\mathbf{h}}_{\rm BLM}$ (72)
$\displaystyle=\left(\left(\mathbf{S}_{\rm
w1}^{\mathrm{T}}\begin{bmatrix}1&\frac{\arcsin\left(\gamma_{s}\right)}{\pi/2}\\\
\frac{\arcsin\left(\gamma_{s}\right)}{\pi/2}&1\end{bmatrix}^{-1}\right)\otimes\mathbf{I}_{N_{R}}\right)\frac{\mathbf{r}}{\sqrt{\pi}}.$
(73)

###  V-D Spatially White Channel and $\tau=2$ Complex-Valued Pilots with
$\pi/2$ Phase Difference

Let $\tau=2$ and $\bm{\Sigma}=\mathbf{I}_{N_{T}N_{R}}$. By assumption,

$\displaystyle\mathbf{S}^{\mathrm{T}}=[\mathbf{s}_{1}\;\;\mathbf{s}_{2}],\quad\frac{\mathbf{s}_{2}}{\|\mathbf{s}_{2}\|}=j\frac{\mathbf{s}_{1}}{\|\mathbf{s}_{1}\|}.$
(74)

Then,

$\displaystyle\mathbf{S}\mathbf{S}^{\mathrm{H}}$
$\displaystyle=\begin{bmatrix}\|\mathbf{s}_{1}\|^{2}&-j\|\mathbf{s}_{1}\|\|\mathbf{s}_{2}\|\\\
j\|\mathbf{s}_{1}\|\|\mathbf{s}_{2}\|&\|\mathbf{s}_{2}\|^{2}\end{bmatrix},$
(75) $\displaystyle\mathbf{D}_{\mathsf{R}}$
$\displaystyle=\frac{\mathtt{diag}{{\left(\left[\|\mathbf{s}_{2}\|^{2}+\sigma^{2}\;\;\|\mathbf{s}_{1}\|^{2}+\sigma^{2}\right]^{\mathrm{T}}\right)}}\otimes\mathbf{I}_{N_{R}}}{\sigma^{2}\left[\sigma^{2}+\textrm{tr}(\mathbf{S}\mathbf{S}^{\mathrm{H}})\right]},$
(76) $\displaystyle\mathbf{D}_{\mathsf{I}}$
$\displaystyle=\frac{\begin{bmatrix}0&\|\mathbf{s}_{1}\|\|\mathbf{s}_{2}\|\\\
-\|\mathbf{s}_{1}\|\|\mathbf{s}_{2}\|&0\end{bmatrix}\otimes\mathbf{I}_{N_{R}}}{\sigma^{2}\left(\sigma^{2}+\textrm{tr}(\mathbf{S}\mathbf{S}^{\mathrm{H}})\right)}.$
(77)

Here, (60) no longer holds. However, by defining

$\displaystyle\bm{\Delta}_{\rm X}$
$\displaystyle=\mathtt{diag}{{\left(\left[r_{\mathsf{R},i}r_{\mathsf{I},N_{R}+i}\right]_{i=1}^{N_{R}}\right)}},$
(78) $\displaystyle\bm{\Delta}_{\rm Y}$
$\displaystyle=\mathtt{diag}{{\left(\left[r_{\mathsf{I},i}r_{\mathsf{R},N_{R}+i}\right]_{i=1}^{N_{R}}\right)}},$
(79) $\displaystyle\mathbf{C}_{{\rm sub,A}}$
$\displaystyle=\mathtt{DIAG}{{}}\left((\|\mathbf{s}_{2}\|^{2}+\sigma^{2})\mathbf{I}_{N_{R}},(\|\mathbf{s}_{1}\|^{2}+\sigma^{2})\mathbf{I}_{N_{R}}\right),$
(80) $\displaystyle\mathbf{C}_{{\rm sub,B}}$
$\displaystyle=\begin{bmatrix}\mathbf{0}&-\|\mathbf{s}_{1}\|\|\mathbf{s}_{2}\|\bm{\Delta}_{\rm
X}\\\ \|\mathbf{s}_{1}\|\|\mathbf{s}_{2}\|\bm{\Delta}_{\rm
Y}&\mathbf{0}\end{bmatrix},$ (81)

we have

$\displaystyle\mathbf{C}=\frac{1}{\sigma^{2}\left(\sigma^{2}+\textrm{tr}(\mathbf{S}\mathbf{S}^{\mathrm{H}})\right)}\begin{bmatrix}\mathbf{C}_{\rm
sub,A}&\mathbf{C}_{\rm sub,B}\\\ \mathbf{C}_{\rm
sub,B}^{\mathrm{T}}&\mathbf{C}_{\rm sub,A}\end{bmatrix},$ (82)

which clearly satisfies the condition in Theorem 1. Reusing the definitions
from (71) in the new context here, denote
$\widetilde{\gamma}_{s}=\|\mathbf{s}_{1,{\rm new}}\|\|\mathbf{s}_{2,{\rm
new}}\|$ and $\mathbf{S}_{\rm w2}^{\mathrm{H}}=\left[\mathbf{s}_{1,{\rm
new}}^{*}\;\;\mathbf{s}_{2,{\rm new}}^{*}\right]$. Then, (31) yields

$\displaystyle\hat{\mathbf{h}}_{\rm MMSE}=\hat{\mathbf{h}}_{\rm BLM}$ (83)
$\displaystyle=\left(\left(\mathbf{S}_{\rm
w2}^{\mathrm{H}}\begin{bmatrix}1&\frac{-j\arcsin\left(\widetilde{\gamma}_{s}\right)}{\pi/2}\\\
\frac{j\arcsin\left(\widetilde{\gamma}_{s}\right)}{\pi/2}&1\end{bmatrix}^{-1}\right)\otimes\mathbf{I}_{N_{R}}\right)\frac{\mathbf{r}}{\sqrt{\pi}}.$
(84)

######  Corollary 1.

With $N_{T}=1$, $\tau=2$, and arbitrary $N_{R}$, $\hat{\mathbf{h}}_{\rm BLM}$
and $\hat{\mathbf{h}}_{\rm MMSE}$ are equivalent for spatially white channels
if QPSK symbols are used as pilot symbols.

###### Proof:

If the phase difference of the two pilot symbol $s_{1}$ and $s_{2}$ is $\pi$,
the result in Section V-C applies. If the phase difference is $\frac{\pi}{2}$,
the result in Section V-D applies. In both cases, the condition in Theorem 1
is satisfied. ∎

###  V-E SIMO with Real-Valued Channel Correlation, $N_{R}=2$, and
$\tau=N_{T}=1$

The non-linearity of the MMSE channel estimate for noiseless correlated SIMO
channels was considered in [19]. Here, we focus on a general noisy correlated
SIMO channel with a single complex-valued pilot symbol (i.e., $\tau=N_{T}=1$)
and the system model (2) becomes

$\displaystyle\mathbf{y}=s\mathbf{h}+\mathbf{n}.$ (85)

Let $\bm{\Sigma}$ be real-valued, e.g., as in the exponential channel model
[36, 37]. Then,

$\displaystyle\mathbf{D}_{\mathsf{R}}^{-1}=\bm{\Omega}_{\rm
b}=|s|^{2}\bm{\Sigma}+\sigma^{2}\mathbf{I}_{N_{R}},\quad\mathbf{D}_{\mathsf{I}}=\mathbf{0},$
(86)

and (60) holds. Furthermore,
$\mathbf{C}=\frac{1}{2}\mathtt{DIAG}{{(}}\mathbf{W}_{\mathsf{R}}^{-1},\mathbf{W}_{\mathsf{I}}^{-1})$,
with

$\displaystyle\mathbf{W}_{\mathsf{R}}=\frac{1}{2}\bm{\Lambda}_{\mathsf{R}}\mathbf{D}_{\mathsf{R}}^{-1}\bm{\Lambda}_{\mathsf{R}},\quad\mathbf{W}_{\mathsf{I}}=\frac{1}{2}\bm{\Lambda}_{\mathsf{I}}\mathbf{D}_{\mathsf{R}}^{-1}\bm{\Lambda}_{\mathsf{I}}.$
(87)

Hence, (60) is equivalent to

$\displaystyle\hat{\mathbf{h}}_{\rm MMSE}$
$\displaystyle=s^{*}\bm{\Sigma}\mathbf{D}_{\mathsf{R}}\left(\frac{\bm{\Lambda}_{\mathsf{R}}\displaystyle\int_{\mathbb{R}_{+}^{N_{R}}}\mathbf{x}e^{-\frac{\mathbf{x}^{\mathrm{T}}\mathbf{W}_{\mathsf{R}}^{-1}\mathbf{x}}{2}}{\rm
d}\mathbf{x}}{\displaystyle\int_{\mathbb{R}_{+}^{N_{R}}}e^{-\frac{\mathbf{x}^{\mathrm{T}}\mathbf{W}_{\mathsf{R}}^{-1}\mathbf{x}}{2}}{\rm
d}\mathbf{x}}\right.$ $\displaystyle\phantom{=}\
+\left.\frac{j\bm{\Lambda}_{\mathsf{I}}\displaystyle\int_{\mathbb{R}_{+}^{N_{R}}}\mathbf{y}e^{-\frac{\mathbf{y}^{\mathrm{T}}\mathbf{W}_{\mathsf{I}}^{-1}\mathbf{y}}{2}}{\rm
d}\mathbf{y}}{\displaystyle\int_{\mathbb{R}_{+}^{N_{R}}}e^{-\frac{\mathbf{y}^{\mathrm{T}}\mathbf{W}_{\mathsf{I}}^{-1}\mathbf{y}}{2}}{\rm
d}\mathbf{y}}\right).$ (88)

When $N_{R}=2$, given the structure of $\mathbf{C}$ above, Theorem 1 holds
trivially. Consistent with this observation, based on Lemma 1, the two ratios
in (88) are linear in $\mathbf{r}_{\mathsf{R}}$ and $\mathbf{r}_{\mathsf{I}}$,
respectively. Assume a standardized $\bm{\Sigma}$ and denote the correlation
coefficient as $\rho_{12}$. Defining
$\beta_{12}=\frac{\rho_{12}|s|^{2}}{|s|^{2}+\sigma^{2}}$, we obtain

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\frac{s^{*}\bm{\Sigma}\begin{bmatrix}1&\frac{\arcsin(\beta_{12})}{\pi/2}\\\
\frac{\arcsin(\beta_{12})}{\pi/2}&1\end{bmatrix}^{-1}\mathbf{r}}{\sqrt{\pi(|s|^{2}+\sigma^{2})}}=\hat{\mathbf{h}}_{\rm
BLM}.$ (89)

##  VI Cases with Nonlinear MMSE Estimators

We now consider two cases where the assumption of Theorem 1 does not hold and
thus $\hat{\mathbf{h}}_{\rm MMSE}$ is non-linear in $\mathbf{r}$. For general
MIMO systems with arbitrary channel correlation and pilots, one can always
resort to our explicit expressions in (40)–(41) and use existing methods
discussed in Section IV (e.g., the algorithm for Monte Carlo simulations in
[31]) to compute the required Gaussian orthant probabilities for one
covariance matrix of dimension $(2\tau N_{R})\times(2\tau N_{R})$ and $2\tau
N_{R}$ covariance matrices of dimension $(2\tau N_{R}-1)\times(2\tau
N_{R}-1)$. Clearly, when either $\tau$ or $N_{R}$ is large, Monte Carlo
methods become computationally cumbersome and cannot guarantee accuracy.
Nevertheless, using our insights above (e.g., Remark 2) and existing results
on Gaussian orthant probabilities, computationally efficient analytical
expressions for $\hat{\mathbf{h}}_{\rm MMSE}$ can be derived for various SIMO
channels under specific assumptions on the channel correlation or pilot
symbols.

###  VI-A SIMO with Real-Valued Correlation, $N_{R}>2$, and $\tau=N_{T}=1$

Consider the same SIMO system described by (85) in Section V-E with
$\tau=N_{T}=1$ and $N_{R}>2$. In this case, (86)–(88) remain applicable.
Without loss of generality, let the channel covariance matrix $\bm{\Sigma}$ be
standardized and let the correlation coefficients be denoted as $\rho_{ik}$,
for $i,k=1,\ldots,N_{R}$ and $i\neq k$, unless otherwise stated.

Let $\mathbf{D}_{\mathsf{R},-k}$ denote the $(N_{R}-1)\times(N_{R}-1)$ matrix
obtained by deleting the $k$-th row and column of $\mathbf{D}_{\mathsf{R}}$.
Likewise, let $\bm{\Lambda}_{\mathsf{R},-k}$ (resp.
$\bm{\Lambda}_{\mathsf{I},-k}$) denote the $(N_{R}-1)\times(N_{R}-1)$ matrix
obtained by deleting the $k$-th row and column of $\bm{\Lambda}_{\mathsf{R}}$
(resp. $\bm{\Lambda}_{\mathsf{I}}$). Moreover, for $k=1,\ldots,N_{R}$, define

$\displaystyle\mathbf{W}_{\mathsf{x},k}$
$\displaystyle=\frac{1}{2}\bm{\Lambda}_{\mathsf{R},-k}\mathbf{D}_{\mathsf{R},-k}^{-1}\bm{\Lambda}_{\mathsf{R},-k},$
(90) $\displaystyle\mathbf{W}_{\mathsf{y},k}$
$\displaystyle=\frac{1}{2}\bm{\Lambda}_{\mathsf{I},-k}\mathbf{D}_{\mathsf{R},-k}^{-1}\bm{\Lambda}_{\mathsf{I},-k}.$
(91)

Using the definitions of $\mathbf{W}_{\mathsf{R}}$ and
$\mathbf{W}_{\mathsf{I}}$ in (87), it is shown in Appendix C that

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\frac{s^{*}\bm{\Sigma}}{2\sqrt{\pi\left(|s|^{2}+\sigma^{2}\right)}}\left(\frac{\bm{\Lambda}_{\mathsf{R}}{\mathbf{g}}_{\rm
R}}{\mathcal{P}(\mathbf{W}_{\mathsf{R}})}+\frac{j\bm{\Lambda}_{\mathsf{I}}{\mathbf{g}}_{\rm
I}}{\mathcal{P}(\mathbf{W}_{\mathsf{I}})}\right),$ (92)

where the $k$-th element of $\mathbf{g}_{\rm R}$ (resp. $\mathbf{g}_{\rm I}$)
is equal to $\mathcal{P}(\mathbf{W}_{\mathsf{x},k})$ (resp.
$\mathcal{P}(\mathbf{W}_{\mathsf{y},k})$). When $N_{R}=2$, (92) reduces to
(89) in Section V-E. When $N_{R}\geq 3$, (92) is the starting point for
further results below.

######  Theorem 2.

Consider the one-bit quantized, spatially correlated Rayleigh fading SIMO
system given by (85) with $N_{T}=\tau=1$ and $N_{R}=3$. Then,

$\hat{\mathbf{h}}_{\rm
MMSE}=\frac{s^{*}\bm{\Sigma}}{2\sqrt{\pi\left(|s|^{2}+\sigma^{2}\right)}}\left(\frac{\mathbf{v}_{\rm
3R}}{P_{\rm 3R}}+j\frac{\mathbf{v}_{\rm 3I}}{P_{\rm 3I}}\right),$ (93) with
$\displaystyle\beta_{ik}$
$\displaystyle=\frac{\rho_{ik}|s|^{2}}{|s|^{2}+\sigma^{2}},\quad(i\neq
k,\;i,k=1,2,3)$ (93a) $\displaystyle\mathbf{v}_{\rm 3R}$
$\displaystyle=\begin{bmatrix}\frac{r_{\mathsf{R},1}}{4}+\frac{r_{\mathsf{R},1}r_{\mathsf{R},2}r_{\mathsf{R},3}}{2\pi}\arcsin\left(\frac{\beta_{23}-\beta_{12}\beta_{13}}{\sqrt{1-\beta_{12}^{2}}\sqrt{1-\beta_{13}^{2}}}\right)\\\
\frac{r_{\mathsf{R},2}}{4}+\frac{r_{\mathsf{R},1}r_{\mathsf{R},2}r_{\mathsf{R},3}}{2\pi}\arcsin\left(\frac{\beta_{13}-\beta_{12}\beta_{23}}{\sqrt{1-\beta_{12}^{2}}\sqrt{1-\beta_{23}^{2}}}\right)\\\
\frac{r_{\mathsf{R},3}}{4}+\frac{r_{\mathsf{R},1}r_{\mathsf{R},2}r_{\mathsf{R},3}}{2\pi}\arcsin\left(\frac{\beta_{12}-\beta_{13}\beta_{23}}{\sqrt{1-\beta_{13}^{2}}\sqrt{1-\beta_{23}^{2}}}\right)\end{bmatrix},$
(93b) $\displaystyle P_{\rm 3R}$
$\displaystyle=\frac{1}{8}+\frac{\sum_{i=1}^{2}\sum_{k=i+1}^{3}r_{\mathsf{R},i}r_{\mathsf{R},k}\arcsin(\beta_{ik})}{4\pi},$
(93c)

and where $\mathbf{v}_{\rm 3I}$ (resp. $P_{\rm 3I}$) is the same as
$\mathbf{v}_{\rm 3R}$ (resp. $P_{\rm 3R}$) except that $r_{\mathsf{R},i}$ in
(93b) (resp. (93c)) is replaced by $r_{\mathsf{I},i}$ for $i=1,2,3$.

###### Proof:

The proof starts from (92) and is based on (44), (45), and (46). The detailed
calculations are similar to those in Appendix B, hence omitted for brevity. ∎

On the other hand, under the same assumptions for Theorem 2, the BLMMSE
channel estimate is given by

$\displaystyle\hat{\mathbf{h}}_{\rm
BLM}=\frac{s^{*}\bm{\Sigma}\begin{bmatrix}1&\frac{\arcsin(\beta_{12})}{\pi/2}&\frac{\arcsin(\beta_{13})}{\pi/2}\\\
\frac{\arcsin(\beta_{12})}{\pi/2}&1&\frac{\arcsin(\beta_{23})}{\pi/2}\\\
\frac{\arcsin(\beta_{13})}{\pi/2}&\frac{\arcsin(\beta_{23})}{\pi/2}&1\end{bmatrix}^{-1}\mathbf{r}}{\sqrt{\pi\left(|s|^{2}+\sigma^{2}\right)}}.$
(94)

In general, (94) is not equivalent to (93). However, we have
$\hat{\mathbf{h}}_{\rm BLM}=\hat{\mathbf{h}}_{\rm MMSE}$ if one of the
following holds: $\rho_{23}=\rho_{13}=0$, $\rho_{13}=\rho_{12}=0$, or
$\rho_{12}=\rho_{23}=0$. This is consistent with Theorem 1.

######  Remark 7.

Theorem 2 is significant that, for one-bit ADCs, it provides the simplest
closed-form expression for $\hat{\mathbf{h}}_{\rm MMSE}$ that demonstrates
non-linear dependence on $\mathbf{r}$.

When $N_{R}=4$, for the same SIMO system given by (85) with a real-valued
channel covariance matrix, $\mathcal{P}(\mathbf{W}_{\mathsf{R}})$ and
$\mathcal{P}(\mathbf{W}_{\mathsf{I}})$ in (92) can be obtained based on (109)
in Appendix A. In addition, for the $3\times 3$ matrices
$\mathbf{W}_{\mathsf{x},k}$ and $\mathbf{W}_{\mathsf{y},k}$,
$\mathcal{P}(\mathbf{W}_{\mathsf{x},k})$ and
$\mathcal{P}(\mathbf{W}_{\mathsf{y},k})$ required in (92) are calculated using
(46), $k=1,\ldots,4$. Indeed, one can obtain a semi-closed-form expression for
(92) with only one-dimensional finite integration from (109). Details are
omitted for brevity.

Next, we consider the SIMO system in (85) with arbitrary $N_{R}$ but with
equal positive channel cross-correlation for all pairs of receive antennas.
Let $Q(\cdot)$ denote the $Q$-function, i.e., the right tail distribution
function of the standard normal distribution $\mathcal{N}(0,1)$.

######  Theorem 3.

Consider the one-bit quantized, spatially correlated Rayleigh fading SIMO
system given by (85) with $N_{T}=\tau=1$. Let the correlation coefficients in
$\bm{\Sigma}$ be given by $\rho_{ik}=\rho$ (with $0\leq\rho\leq 1$), for
$i,k=1,\ldots,N_{R}$ and $i\neq k$. Then,

$\displaystyle\hat{\mathbf{h}}_{\rm
MMSE}=\frac{s^{*}\bm{\Sigma}}{2\sqrt{\pi\left(|s|^{2}+\sigma^{2}\right)}}\left(\frac{\bm{\Lambda}_{\mathsf{R}}\mathbf{v}_{\rm
RE}}{P_{\rm RE}}+\frac{j\bm{\Lambda}_{\mathsf{I}}\mathbf{v}_{\rm IM}}{P_{\rm
IM}}\right),$ (95)

with

$\displaystyle P_{\rm
RE}=\displaystyle\int_{-\infty}^{\infty}\prod_{k=1}^{N_{R}}Q\left(\frac{r_{\mathsf{R},k}\sqrt{\rho}|s|t}{\sqrt{|s|^{2}(1-\rho)+\sigma^{2}}}\right)\frac{e^{-t^{2}/2}}{\sqrt{2\pi}}{\rm
d}t,$ (96)

and where $P_{\rm IM}$ is the same as $P_{\rm RE}$ except that
$r_{\mathsf{R},k}$ in (96) is replaced by $r_{\mathsf{I},k}$, for
$k=1,\ldots,N_{R}$. The $k$-th element of $\mathbf{v}_{\rm RE}$ is given by

$\displaystyle v_{{\rm
RE},k}=\displaystyle\int_{-\infty}^{\infty}\prod_{i=1,i\neq
k}^{N_{R}}Q\left(\frac{r_{\mathsf{R},i}\sqrt{\rho}|s|t}{\sqrt{|s|^{2}+\sigma^{2}}}\right)\frac{e^{-t^{2}/2}}{\sqrt{2\pi}}{\rm
d}t,$ (97)

for $k=1,\ldots,N_{R}$. The $k$-th element of $\mathbf{v}_{\rm IM}$, i.e.,
$v_{{\rm IM},k}$, is the same as $v_{{\rm RE},k}$ except that
$r_{\mathsf{R},i}$ is replaced by $r_{\mathsf{I},i}$ for $i$ in (97) and
$k=1,\ldots,N_{R}$.

###### Proof:

The proof utilizes the key observation made in Remark 2 that each off-diagonal
element of $\mathbf{C}$ in (28) is affected by exactly two different elements
from
$\left[\mathbf{r}_{\mathsf{R}}^{\mathrm{T}}\;\;\mathbf{r}_{\mathsf{I}}^{\mathrm{T}}\right]^{\mathrm{T}}$.
In particular, the off-diagonal elements of $\mathbf{W}_{\mathsf{R}}$ and
$\mathbf{W}_{\mathsf{I}}$ in (87), as well as those of
$\mathbf{W}_{\mathsf{x},k}$ and $\mathbf{W}_{\mathsf{y},k}$ in (90)–(91), for
$k=1,\ldots,N_{R}$, are each affected by one pair of elements in
$\mathbf{r}_{\mathsf{R}}$ and $\mathbf{r}_{\mathsf{I}}$, respectively. Based
on this key fact, the proof further invokes the results in [30, p. 192]. A
detailed proof is provided in Appendix D. ∎

######  Remark 8.

Both (96) and (97)) require integration with only one variable. Therefore,
(95) can be conveniently evaluated using common computing software even for
large $N_{R}$ [30].

###  VI-B SISO and Spatially White SIMO with Real-Valued Pilots, Arbitrary
$\tau$, and $N_{T}=1$

For spatially white SIMO systems, the MMSE channel estimate derived for one
receive antenna holds for all receive antennas. This allows us to focus on
SISO channels with multiple real-valued pilot symbols. Denote the real-valued
pilot vector $\mathbf{s}$ as

$\displaystyle\mathbf{s}=[s_{1}\;\ldots\;s_{\tau}]^{\mathrm{T}}.$ (98)

When $\tau=1$, this scenario becomes a special case of Section V-A with
$\tau=N_{T}=N_{R}=1$. Furthermore, when $\tau=2$, the setup reduces to that in
Section V-C with $N_{T}=N_{R}=1$. Therefore, the non-trivial new case here
corresponds to $\tau>2$, although subsequent results are applicable to
$\tau=1,2$.

When $N_{R}=N_{T}=1$, the system model (2) reduces to
$\mathbf{y}=h\mathbf{s}+\mathbf{n}$, with $h\sim\mathcal{CN}(0,1)$.
Correspondingly, we have

$\displaystyle\mathbf{D}_{\mathsf{R}}=\left(\mathbf{s}\mathbf{s}^{\mathrm{T}}+\sigma^{2}\mathbf{I}_{\tau}\right)^{-1},\quad\mathbf{D}_{\mathsf{I}}=\mathbf{0}.$
(99)

Using the Sherman-Morrison formula [34, p. 124], we have
$\mathbf{D}_{\mathsf{R}}=\frac{1}{\sigma^{2}}\left[\mathbf{I}_{\tau}-\frac{\mathbf{s}\mathbf{s}^{\mathrm{T}}}{\mathbf{s}^{\mathrm{T}}\mathbf{s}+\sigma^{2}}\right]$.
Let $\widetilde{\mathbf{D}}_{\mathsf{R},-k}$ denote the
$(\tau-1)\times(\tau-1)$ matrix obtained by deleting the $k$-th row and column
of $\mathbf{D}_{\mathsf{R}}$ in (99). Then,

$\displaystyle\widetilde{\mathbf{D}}_{\mathsf{R},-k}$
$\displaystyle=\frac{1}{\sigma^{2}}\left(\mathbf{I}_{\tau-1}-\frac{\mathbf{s}_{-k}\mathbf{s}_{-k}^{\mathrm{T}}}{\mathbf{s}^{\mathrm{T}}\mathbf{s}+\sigma^{2}}\right),$
(100) $\displaystyle\widetilde{\mathbf{D}}_{\mathsf{R},-k}^{-1}$
$\displaystyle=\sigma^{2}\left(\mathbf{I}_{\tau-1}+\frac{\mathbf{s}_{-k}\mathbf{s}_{-k}^{\mathrm{T}}}{s_{k}^{2}+\sigma^{2}}\right),$
(101)

where
$\mathbf{s}_{-k}=[s_{1}\;\ldots\;s_{k-1}\;s_{k+1}\;\ldots\;s_{\tau}]^{\mathrm{T}}$.
Define

$\displaystyle\mathbf{W}_{\mathsf{A}}=\frac{1}{2}\bm{\Lambda}_{\mathsf{R}}\mathbf{D}_{\mathsf{R}}^{-1}\bm{\Lambda}_{\mathsf{R}},\quad\mathbf{W}_{\mathsf{B}}=\frac{1}{2}\bm{\Lambda}_{\mathsf{I}}\mathbf{D}_{\mathsf{R}}^{-1}\bm{\Lambda}_{\mathsf{I}},$
(102)

and let $\widetilde{\bm{\Lambda}}_{\mathsf{R},-k}$ (resp.
$\widetilde{\bm{\Lambda}}_{\mathsf{I},-k}$) denote the
$(\tau-1)\times(\tau-1)$ matrix obtained by deleting the $k$-th row and column
of $\bm{\Lambda}_{\mathsf{R}}$ (resp. $\bm{\Lambda}_{\mathsf{I}}$). Lastly,
define

$\displaystyle\mathbf{W}_{\mathsf{A},k}$
$\displaystyle=\frac{1}{2}\widetilde{\bm{\Lambda}}_{\mathsf{R},-k}\widetilde{\mathbf{D}}_{\mathsf{R},-k}^{-1}\widetilde{\bm{\Lambda}}_{\mathsf{R},-k},$
(103) $\displaystyle\mathbf{W}_{\mathsf{B},k}$
$\displaystyle=\frac{1}{2}\widetilde{\bm{\Lambda}}_{\mathsf{I},-k}\widetilde{\mathbf{D}}_{\mathsf{R},-k}^{-1}\widetilde{\bm{\Lambda}}_{\mathsf{I},-k}.$
(104)

######  Theorem 4.

For a Rayleigh fading SISO channel with a real-valued pilot vector and with
one-bit quantization at the receiver, the MMSE channel estimate is given by

$\displaystyle\hat{h}_{\rm
MMSE}=\frac{1}{2\sqrt{\pi}}\sum_{k=1}^{\tau}\frac{s_{k}\left(\frac{r_{\mathsf{R},k}\mathcal{P}(\mathbf{W}_{\mathsf{A},k})}{\mathcal{P}(\mathbf{W}_{\mathsf{A}})}+\frac{jr_{\mathsf{I},k}\mathcal{P}(\mathbf{W}_{\mathsf{B},k})}{\mathcal{P}(\mathbf{W}_{\mathsf{B}})}\right)}{\sqrt{s_{k}^{2}+\sigma^{2}}},$
(105)

with

$\displaystyle\mathcal{P}(\mathbf{W}_{\mathsf{A}})$
$\displaystyle=\displaystyle\int_{-\infty}^{\infty}\prod_{i=1}^{\tau}Q\left(\frac{r_{\mathsf{R},i}s_{i}t}{\sigma}\right)\frac{e^{-\frac{t^{2}}{2}}}{\sqrt{2\pi}}{\rm
d}t,$ (106) $\displaystyle\mathcal{P}(\mathbf{W}_{\mathsf{A},k})$
$\displaystyle=\displaystyle\int_{-\infty}^{\infty}\prod_{i=1,i\neq
k}^{\tau}Q\left(\frac{r_{\mathsf{R},i}s_{i}t}{\sqrt{s_{k}^{2}+\sigma^{2}}}\right)\frac{e^{-\frac{t^{2}}{2}}}{\sqrt{2\pi}}{\rm
d}t,$ (107)

for $k=1,\ldots,\tau$. Note that, $\forall k$,
$\mathcal{P}(\mathbf{W}_{\mathsf{B},k})$ is the same as
$\mathcal{P}(\mathbf{W}_{\mathsf{A},k})$ in (107) except that
$r_{\mathsf{R},i}$ is replaced by $r_{\mathsf{I},i}$, for $i\neq k$.
Similarly, $\mathcal{P}(\mathbf{W}_{\mathsf{B}})$ is the same as
$\mathcal{P}(\mathbf{W}_{\mathsf{A}})$ in (106) except that $r_{\mathsf{R},i}$
is replaced by $r_{\mathsf{I},i}$, $i=1,\ldots,\tau$.

###### Proof:

The main idea is similar to that for the proof of Theorem 3. Details are given
in Appendix E. ∎

######  Remark 9.

Let $N_{T}=N_{R}=1$. As expected, when $\tau=1$, (105) is equivalent to (61),
whereas when $\tau=2$, (105) is the same as (73). In both cases, $\hat{h}_{\rm
MMSE}=\hat{h}_{\rm BLM}$, and the corresponding asymptotic MSE is given by

$\displaystyle\lim\limits_{\sigma^{2}\rightarrow 0}\mathsf{E}(|\hat{h}_{\rm
MMSE}-h|^{2})=\lim\limits_{\sigma^{2}\rightarrow 0}\mathsf{E}(|\hat{h}_{\rm
BLM}-h|^{2})=1-\frac{2}{\pi}.$ (108)

In fact, from both (73) and (84), (108) is the asymptotic MSE per antenna for
both estimators when the noise power is sufficiently low.

##  VII Numerical Results

We now provide simulation results to complement the above analyses. The BLMMSE
channel estimate is obtained using (52) from [7] and the corresponding
theoretical MSE per antenna can be readily calculated as in [17]. For cases
with a single pilot, the signal-to-noise ratio (SNR) is defined as
$\textrm{SNR}=|s|^{2}/\sigma^{2}$. For cases with a $\tau$-dimensional pilot
vector $\mathbf{s}$, $\textrm{SNR}=\|\mathbf{s}\|^{2}/(\tau\sigma^{2})$. When
a pilot matrix is used,
$\textrm{SNR}=\textrm{tr}(\mathbf{S}\mathbf{S}^{\mathrm{H}})/(\tau
N_{T}\sigma^{2})$. The MSE per antenna is averaged over many fading
realizations. Without loss of generality, the channel covariance matrices used
here are standardized.

###  VII-A MISO and MIMO with No Receive Correlation and Unitary Pilot Matrix

Here, we investigate the impact of transmit antenna correlation and the number
of transmit antennas $N_{T}$ for the case considered in Section V-B with
$\hat{\mathbf{h}}_{\rm BLM}=\hat{\mathbf{h}}_{\rm MMSE}$. In the absence of
channel correlation at the receiver for MIMO systems, the MSE is the same for
all receive antennas. Therefore, a MISO setup is considered for simplicity.
For transmit channel correlation in (62), we use the exponential correlation
model [37, 36] and let $[\bm{\Sigma}_{\rm TX}]_{ik}=0.5^{|i-k|}$ or
$[\bm{\Sigma}_{\rm TX}]_{ik}=0.9^{|i-k|}$, for $i,k=1,\ldots,N_{T}$. The
results are shown in Fig. 1. Clearly, increasing both $N_{T}$ and the amount
of transmit correlation improves the MSE performance at low SNR. At high SNR,
the MSE per antenna converges to $1-\frac{2}{\pi}$, confirming the claim in
Remark 6. Similar trends are observed for other correlation models such as
those in [38, 39].

Figure 1: MSE comparison for MISO case with $\tau=N_{T}=16,32$ and exponential
transmit correlation model using a tailored pilot matrix. In this scenario
$\hat{\mathbf{h}}_{\rm MMSE}=\hat{\mathbf{h}}_{\rm BLM}$.

###  VII-B SIMO with $\tau=1$ and $N_{R}=3,4$

Here, we assume $N_{R}=3$ and $N_{R}=4$ corresponding to the discussion in
Section VI-A. Fig. 2 shows the results using the exponential channel
correlation model [36, 37] with $\rho_{ik}=0.5^{|i-k|}$ or
$\rho_{ik}=0.9^{|i-k|}$, $i,k=1,\ldots,N_{R}$. We observe that the difference
in the MSE performance between $\hat{\mathbf{h}}_{\rm BLM}$ and
$\hat{\mathbf{h}}_{\rm MMSE}$ is more pronounced at moderate-to-high SNR when
both the channel correlation coefficient and $N_{R}$ grow large. When the
channel correlation is low, the performance of the two is similar, as expected
according to Theorem 1 (see Remark 5).

Figure 2: MSE comparison for SIMO case with $N_{R}=3,4$ and exponential
channel correlation model.

###  VII-C SIMO with $\tau=1$ and $\rho_{i,k}=\rho$

Figure 3: MSE comparison for SIMO case with $N_{R}=4,16,32$ and the same
channel correlation coefficient $\rho=0.9$ between all antennas. Figure 4: MSE
comparison for SIMO case with $N_{R}=8,16,32,64$ and the same channel
correlation coefficient $\rho=0.9$ between all antennas.

Using Theorem 3, we provide simulation results in Fig. 3 for $N_{R}=4,16,32$
and in Fig. 4 for $N_{R}=8,16,32,64$, with $\rho_{i,k}=\rho=0.9$, for $i\neq
k$. In Fig. 3, the curves for $N_{R}=4$ are obtained using the following three
methods: 1) based on Theorem 3, with the corresponding curves labeled
“Theorem”; 2) based on the semi-closed-form expression that relies on Appendix
A [27], with the corresponding curve labeled “Childs”; and 3) based on (92)
along with the algorithm in [31] to obtain the required orthant probabilities
via Monte Carlo simulations, with the corresponding curve labeled “Genz”. The
results from all three methods coincide, which corroborates the analysis
leading to Theorem 3. However, when $N_{R}$ becomes larger, the semi-closed-
form expressions do not exist and the Monte Carlo method [31] also becomes
computationally more demanding and less accurate. Consequently, we resort to
Theorem 3 to obtain all the curves for the MMSE channel estimate in Fig. 4. We
see in Figs. 2–4 that when the channel correlation is high,
$\hat{\mathbf{h}}_{\rm MMSE}$ clearly outperforms $\hat{\mathbf{h}}_{\rm BLM}$
across a wide range of SNRs, particularly for moderate-to-high SNRs. The per-
antenna MSE of $\hat{\mathbf{h}}_{\rm MMSE}$ improves considerably with
$N_{R}$ at all SNR values, whereas for $\hat{\mathbf{h}}_{\rm BLM}$, the
improvement with $N_{R}$ is predominantly observed for relatively low SNRs.

###  VII-D Spatially White SIMO or SISO with Multiple Real-Valued Pilots

We investigate the effect of using multiple pilot symbols based on Theorem 4.
For the simulations, we use the second column of $\tau\times\tau$ Hadamard
matrices as the pilot vector $\mathbf{s}$. Fig. 5 presents the results for
$\tau=2,16,32$. The curves corresponding to $\tau=2$ corroborate the claim
made in Remark 9 and demonstrate the equivalence of the two estimators as well
as the equivalence of (73) and (105). It is also visible in Figs. 3–4 that the
MSE of both $\hat{\mathbf{h}}_{\rm MMSE}$ and $\hat{\mathbf{h}}_{\rm BLM}$
first improves and then deteriorates with the SNR, since it is noise-limited
at low SNR and quantization-error-limited at high SNR. For $\tau>2$, the two
estimators exhibit the same asymptotic MSE performance at both low and high
SNR. However, significant MSE gains are achievable from non-linear processing
for intermediate SNRs. The simulation results validate the high-SNR asymptotic
analysis for the case with $\tau=2$ in (108). Interestingly, the same behavior
at high SNR is also observed for $\tau>2$.

Figure 5: MSE comparison for SISO or SIMO with no receive correlation,
$\tau=2,16,32$.

##  VIII Conclusions

In this work, we have investigated optimal MMSE channel estimation for MIMO
systems in slow, flat Rayleigh fading channels with one-bit quantization at
the receiver. We have developed an explicit framework for the computation of
the MMSE estimate, which is related to the orthant probability of the MVN
distribution. Leveraging insights gained from this framework, we have fully
unraveled the condition for the BLMMSE channel estimator to be MSE optimal.
Furthermore, we have derived computationally efficient expressions for the
non-linear MMSE estimators under specific assumptions regarding channel
correlation or pilot symbols. Our findings reveal that the MSE performance of
the low-cost BLMMSE estimator closely approximates that of the MMSE estimator
under conditions of low channel correlation, relatively low SNR, or very high
SNR with multiple pilots per transmit antenna. However, with a large number of
highly correlated receive antennas or multiple pilot symbols, the MMSE channel
estimator yields significant gains over the BLMMSE estimator across a wide
range of SNRs. Our results offer unique insights into the one-bit channel
estimation problem and provide guidance for applications when $N_{R}$ or
$\tau$ is large.

##  Appendix A The orthant probability in (36) with $L=4$

Here, we adopt the same notation as in (43). Based on [27], when $L=4$,

$\displaystyle\mathcal{P}(\bm{\Psi})$
$\displaystyle=\frac{1}{16}+\frac{\arcsin(\psi_{34})+\sum_{l=3}^{4}\arcsin(\psi_{2l})}{8\pi}$
$\displaystyle\phantom{=}\
+\frac{\sum_{k=2}^{4}\arcsin(\psi_{1k})}{8\pi}+\frac{\sum_{i=1}^{3}\mathcal{J}_{i}}{4\pi^{2}},$
(109) with $\displaystyle\mathcal{J}_{1}$
$\displaystyle=\displaystyle\int_{0}^{1}\frac{\psi_{12}\arcsin\left(\frac{d_{34,a}-t^{2}d_{34,b}}{\sqrt{\left(\mu_{23,a}-t^{2}\mu_{23,b}\right)\left(\mu_{24,a}-t^{2}\mu_{24,b}\right)}}\right)}{\sqrt{1-t^{2}\psi_{12}^{2}}}{\rm
d}t,$ (109a) $\displaystyle\mathcal{J}_{2}$
$\displaystyle=\displaystyle\int_{0}^{1}\frac{\psi_{13}\arcsin\left(\frac{d_{24,a}-t^{2}d_{24,b}}{\sqrt{\left(\mu_{23,a}-t^{2}\mu_{23,b}\right)\left(\mu_{34,a}-t^{2}\mu_{34,b}\right)}}\right)}{\sqrt{1-t^{2}\psi_{13}^{2}}}{\rm
d}t,$ (109b) $\displaystyle\mathcal{J}_{3}$
$\displaystyle=\displaystyle\int_{0}^{1}\frac{\psi_{14}\arcsin\left(\frac{d_{23,a}-t^{2}d_{23,b}}{\sqrt{\left(\mu_{24,a}-t^{2}\mu_{24,b}\right)\left(\mu_{34,a}-t^{2}\mu_{34,b}\right)}}\right)}{\sqrt{1-t^{2}\psi_{14}^{2}}}{\rm
d}t,$ (109c) and $\displaystyle d_{23,a}$
$\displaystyle=\psi_{23}-\psi_{24}\psi_{34},$ (109d) $\displaystyle d_{24,a}$
$\displaystyle=\psi_{24}-\psi_{23}\psi_{34},$ (109e) $\displaystyle d_{34,a}$
$\displaystyle=\psi_{34}-\psi_{23}\psi_{24},$ (109f) $\displaystyle\mu_{ik,a}$
$\displaystyle=1-\psi_{ik}^{2}\quad(i<k,\;i,k=2,3,4),$ (109g)
$\displaystyle\mu_{ik,b}$
$\displaystyle=\psi_{1i}^{2}+\psi_{1k}^{2}-2\psi_{1i}\psi_{1k}\psi_{ik}\quad(i<k,\;i,k=2,3,4),$
(109h) $\displaystyle d_{23,b}$
$\displaystyle=\psi_{12}\psi_{13}+\psi_{14}^{2}\psi_{23}-\psi_{12}\psi_{14}\psi_{34}-\psi_{13}\psi_{14}\psi_{24},$
(109i) $\displaystyle d_{24,b}$
$\displaystyle=\psi_{12}\psi_{14}+\psi_{13}^{2}\psi_{24}-\psi_{12}\psi_{13}\psi_{34}-\psi_{13}\psi_{14}\psi_{23},$
(109j) $\displaystyle d_{34,b}$
$\displaystyle=\psi_{13}\psi_{14}+\psi_{12}^{2}\psi_{34}-\psi_{12}\psi_{14}\psi_{23}-\psi_{12}\psi_{13}\psi_{24}.$
(109k)

##  Appendix B Computational Considerations for Remark 4

With a slight abuse of notation, let
$\mathbf{q}=[q_{1}\;q_{2}\;q_{3}]^{\mathrm{T}}$, with
$q_{1},q_{2},q_{3}\in\\{\pm 1\\}$. Let the pdf $p(\mathbf{u})$ of the real-
valued Gaussian random vector $\mathbf{u}=[u_{1}\;u_{2}\;u_{3}]^{\mathrm{T}}$
be given by $\mathcal{N}(\mathbf{0},\bm{\Phi})$, with
$[\bm{\Phi}]_{ii}=\sigma_{ui}^{2}$ and
$[\bm{\Phi}]_{ik}=[\bm{\Phi}]_{ki}=q_{i}q_{k}\phi_{ik}\sigma_{ui}\sigma_{uk}$,
$\phi_{ik}\in(-1,1)$, for $i,k=1,2,3$ and $i\neq k$. Let
$\bm{\Lambda}_{\mathbf{q}}=\mathtt{diag}{{(}}\mathbf{q})$. We show through
evaluation that
$\left.\bm{\Lambda}_{\mathbf{q}}\int_{\mathbb{R}_{+}^{3}}\mathbf{u}p(\mathbf{u}){\rm
d}\mathbf{u}\middle/\int_{\mathbb{R}_{+}^{3}}p(\mathbf{u}){\rm
d}\mathbf{u}\right.$ is not in general linear in $\mathbf{q}$. The same claim
holds when there are more than three elements in $\mathbf{q}$ and
$\mathbf{u}$.

Clearly, $\int_{\mathbb{R}_{+}^{3}}p(\mathbf{u}){\rm
d}\mathbf{u}=\mathcal{P}(\bm{\Phi})$. Based on (44) and (46),

$\displaystyle\displaystyle\int_{\mathbb{R}_{+}^{3}}p(\mathbf{u}){\rm
d}\mathbf{u}=\frac{1}{8}+\frac{\sum_{i=1}^{2}\sum_{k=i+1}^{3}q_{i}q_{k}\arcsin(\phi_{ik})}{4\pi}.$
(110)

Using the reduction in (32)–(34),

$\displaystyle\displaystyle\int_{\mathbb{R}_{+}^{3}}\mathbf{u}p(\mathbf{u}){\rm
d}\mathbf{u}=\frac{\bm{\Phi}\left[\displaystyle\int_{\mathbb{R}_{+}^{2}}e^{-\frac{\hat{\mathbf{u}}_{i}^{\mathrm{T}}\bm{\Phi}_{-i}^{-1}\hat{\mathbf{u}}_{i}}{2}}{\rm
d}\hat{\mathbf{u}}_{i}\right]_{i=1}^{3}}{(2\pi)^{\frac{3}{2}}|\bm{\Phi}|^{\frac{1}{2}}},$
(111)

where $\hat{\mathbf{u}}_{i}$ denotes the $2$-dimensional vector derived by
deleting the $i$-th element of $\mathbf{u}$ and $\bm{\Phi}_{-i}^{-1}$ is the
$2\times 2$ matrix obtained by deleting the $i$-th row and column of
$\bm{\Phi}^{-1}$, for $i=1,2,3$. Clearly,

$\displaystyle\displaystyle\int_{\mathbb{R}_{+}^{2}}\exp\left\\{-\frac{1}{2}\hat{\mathbf{u}}_{i}^{\mathrm{T}}\bm{\Phi}_{-i}^{-1}\hat{\mathbf{u}}_{i}\right\\}{\rm
d}\hat{\mathbf{u}}_{i}=2\pi|\bm{\Phi}_{-i}|^{\frac{1}{2}}\mathcal{P}(\bm{\Phi}_{-i}).$
(112)

It can be shown that $|\bm{\Phi}_{-i}|=|\bm{\Phi}|/\sigma_{ui}^{2}$, and from
(45),

$\displaystyle\mathcal{P}(\bm{\Phi}_{-i})=\frac{1}{4}+\frac{q_{k}q_{l}}{2\pi}\arcsin\left(\frac{\phi_{kl}-\phi_{ik}\phi_{il}}{\sqrt{1-\phi_{ik}^{2}}\sqrt{1-\phi_{il}^{2}}}\right),$
(113)

for $i,k,l=1,2,3$, $i\neq k$, $i\neq l$, and $k\neq l$. Define
$\bm{\Lambda}_{\sigma}=\mathtt{diag}{{([\sigma_{u1}\;\;\sigma_{u2}\;\;\sigma_{u3}]^{\mathrm{T}})}}$.
Further capitalizing on the fact that $q_{i}\in\\{\pm 1\\},$ for $i=1,2,3$, we
obtain

$\displaystyle\bm{\Lambda}_{\mathbf{q}}\displaystyle\int_{\mathbb{R}_{+}^{3}}\mathbf{u}p(\mathbf{u}){\rm
d}\mathbf{u}=\frac{\bm{\Lambda}_{\sigma}}{\sqrt{2\pi}}\begin{bmatrix}1&\phi_{12}&\phi_{13}\\\
\phi_{12}&1&\phi_{23}\\\ \phi_{13}&\phi_{23}&1\end{bmatrix}$
$\displaystyle\cdot\begin{bmatrix}\frac{q_{1}}{4}+\frac{q_{1}q_{2}q_{3}}{2\pi}\arcsin\left(\frac{\phi_{23}-\phi_{12}\phi_{13}}{\sqrt{1-\phi_{12}^{2}}\sqrt{1-\phi_{13}^{2}}}\right)\\\
\frac{q_{2}}{4}+\frac{q_{1}q_{2}q_{3}}{2\pi}\arcsin\left(\frac{\phi_{13}-\phi_{12}\phi_{23}}{\sqrt{1-\phi_{12}^{2}}\sqrt{1-\phi_{23}^{2}}}\right)\\\
\frac{q_{3}}{4}+\frac{q_{1}q_{2}q_{3}}{2\pi}\arcsin\left(\frac{\phi_{12}-\phi_{13}\phi_{23}}{\sqrt{1-\phi_{13}^{2}}\sqrt{1-\phi_{23}^{2}}}\right)\end{bmatrix}.$
(114)

From (110) and (B), it is clear that, in general,
$\left.\bm{\Lambda}_{\mathbf{q}}\int_{\mathbb{R}_{+}^{3}}\mathbf{u}p(\mathbf{u}){\rm
d}\mathbf{u}\middle/\int_{\mathbb{R}_{+}^{3}}p(\mathbf{u}){\rm
d}\mathbf{u}\right.$ is not linear in $\mathbf{q}$.

##  Appendix C Derivation of (92)

Recall (85)–(88) and (90)–(91). It is clear that
$\left|\mathbf{W}_{\mathsf{R}}\right|=\left|\mathbf{W}_{\mathsf{I}}\right|=2^{-N_{R}}\left|\bm{\Omega}_{\rm
b}\right|$. Using (32) and (36), and based on (60), we obtain

$\displaystyle\hat{\mathbf{h}}_{\rm MMSE}$
$\displaystyle=s^{*}\bm{\Sigma}\mathbf{D}_{\mathsf{R}}\left(\frac{\bm{\Lambda}_{\mathsf{R}}\displaystyle\int_{\mathbb{R}_{+}^{N_{R}}}\mathbf{x}e^{-\frac{\mathbf{x}^{\mathrm{T}}\mathbf{W}_{\mathsf{R}}^{-1}\mathbf{x}}{2}}{\rm
d}\mathbf{x}}{(2\pi)^{\frac{N_{R}}{2}}\left|\mathbf{W}_{\mathsf{R}}\right|^{\frac{1}{2}}\mathcal{P}(\mathbf{W}_{\mathsf{R}})}\right.$
$\displaystyle\phantom{=}\
+\left.\frac{j\bm{\Lambda}_{\mathsf{I}}\displaystyle\int_{\mathbb{R}_{+}^{N_{R}}}\mathbf{y}e^{-\frac{\mathbf{y}^{\mathrm{T}}\mathbf{W}_{\mathsf{I}}^{-1}\mathbf{y}}{2}}{\rm
d}\mathbf{y}}{(2\pi)^{\frac{N_{R}}{2}}\left|\mathbf{W}_{\mathsf{I}}\right|^{\frac{1}{2}}\mathcal{P}(\mathbf{W}_{\mathsf{I}})}\right)$
(115)
$\displaystyle=\frac{s^{*}\bm{\Sigma}\mathbf{D}_{\mathsf{R}}}{\pi^{\frac{N_{R}}{2}}\left|\bm{\Omega}_{\rm
b}\right|^{\frac{1}{2}}}\left(\frac{\bm{\Lambda}_{\mathsf{R}}\mathbf{W}_{\mathsf{R}}\;\mathbf{g}_{\rm
x}}{\mathcal{P}(\mathbf{W}_{\mathsf{R}})}+\frac{j\bm{\Lambda}_{\mathsf{I}}\mathbf{W}_{\mathsf{I}}\;\mathbf{g}_{\rm
y}}{\mathcal{P}(\mathbf{W}_{\mathsf{I}})}\right),$ (116)

where $\mathbf{g}_{\rm x}$ (resp. $\mathbf{g}_{\rm y}$) is not the same as
$\mathbf{g}_{\rm R}$ (resp. $\mathbf{g}_{\rm I}$) in (92). For
$k=1,\ldots,N_{R}$, the $k$-th element of $\mathbf{g}_{\rm x}$ is given by

$\displaystyle g_{{\rm x},k}$
$\displaystyle=\displaystyle\int_{\mathbb{R}_{+}^{N_{R}-1}}e^{-\frac{1}{2}\mathbf{x}_{-k}^{\mathrm{T}}\mathbf{W}_{\mathsf{x},k}^{-1}\mathbf{x}_{-k}}{\rm
d}\mathbf{x}_{-k}$ (117)
$\displaystyle=(2\pi)^{\frac{N_{R}-1}{2}}\left|\mathbf{W}_{\mathsf{x},k}\right|^{\frac{1}{2}}\mathcal{P}(\mathbf{W}_{\mathsf{x},k}),$
(118)

with
$\mathbf{x}_{-k}=[x_{1}\;\ldots\;x_{k-1}\;x_{k+1}\;\ldots\;x_{N_{R}}]^{\mathrm{T}}$.
Similarly,

$\displaystyle\mathbf{g}_{\rm
y}=\left[(2\pi)^{\frac{N_{R}-1}{2}}\left|\mathbf{W}_{\mathsf{y},k}\right|^{\frac{1}{2}}\mathcal{P}(\mathbf{W}_{\mathsf{y},k})\right]_{k=1}^{N_{R}}.$
(119)

From (39), we obtain
$\left|\mathbf{D}_{\mathsf{R},-k}^{-1}\right|=\frac{|\bm{\Omega}_{\rm
b}|}{|s|^{2}+\sigma^{2}}$, and thus from (90)–(91), for $k=1,\ldots,N_{R}$, we
have

$\displaystyle\left|\mathbf{W}_{\mathsf{x},k}\right|=\left|\mathbf{W}_{\mathsf{y},k}\right|=\frac{2^{-(N_{R}-1)}|\bm{\Omega}_{\rm
b}|}{|s|^{2}+\sigma^{2}}.$ (120)

Based on the above steps, (116) is further simplified to (92).

##  Appendix D Proof of Theorem 3

###### Proof:

Recall (85)–(88) and the starting point (92). Each off-diagonal element of
$\mathbf{W}_{\mathsf{R}}$ (resp. $\mathbf{W}_{\mathsf{I}}$) is affected by two
different elements from $\mathbf{r}_{\mathsf{R}}$ (resp.
$\mathbf{r}_{\mathsf{I}}$). Specifically, the $(i,k)$-th element of the
standardized version of $\mathbf{W}_{\mathsf{R}}$ is given by
$\left(r_{\mathsf{R},i}\sqrt{\frac{|s|^{2}\rho}{|s|^{2}+\sigma^{2}}}\right)\left(r_{\mathsf{R},k}\sqrt{\frac{|s|^{2}\rho}{|s|^{2}+\sigma^{2}}}\right)$,
for $i,k=1,\ldots,N_{R}$ and $i\neq k$. Based on Lemma [30, p. 192] and using
(44), $\mathcal{P}(\mathbf{W}_{\mathsf{R}})$ is given by (96). Similarly,
$\mathcal{P}(\mathbf{W}_{\mathsf{I}})$ is obtained.

By assumption, $\bm{\Sigma}$ is a standardized circulant matrix with all the
off-diagonal elements equal to $\rho$. The inverse of $\bm{\Sigma}$ is also a
circulant matrix [40] given by:

$\displaystyle\bm{\Sigma}^{-1}=\frac{\left[1+(N_{R}-2)\rho\right]\bm{\Sigma}_{\rm
inv,circ}}{1+(N_{R}-2)\rho-(N_{R}-1)\rho^{2}},$ (121)

where $\bm{\Sigma}_{\rm inv,circ}$ is a standardized circulant matrix whose
off-diagonal elements are equal to
$\left(-\frac{\rho}{1+(N_{R}-2)\rho}\right)$. Clearly, here $\bm{\Omega}_{\rm
b}$ and its inverse $\mathbf{D}_{\mathsf{R}}$ are also circulant matrices.
Consider $\mathbf{D}_{\mathsf{R},-k}$, for $k=1,\ldots,N_{R}$. Using (121),
the $(N_{R}-1)\times(N_{R}-1)$ circulant matrix
$\mathbf{D}_{\mathsf{R},-k}^{-1}$ is given by

$\displaystyle\mathbf{D}_{\mathsf{R},-k}^{-1}=\frac{\left((|s|^{2}+\sigma^{2})^{2}-|s|^{4}\rho^{2}\right)\mathbf{D}_{\rm
inv,circ}}{|s|^{2}+\sigma^{2}},$ (122)

where $\mathbf{D}_{\rm inv,circ}$ is a standardized circulant matrix with off-
diagonal elements given by $\frac{\rho|s|^{2}}{[|s|^{2}(1+\rho)]+\sigma^{2}}$.
From (122) and (90), the $(i,l)$-th element of the standardized version of
$\mathbf{W}_{\mathsf{x},k}$ is given by
$\left(r_{\mathsf{R},i}\sqrt{\frac{|s|^{2}\rho}{|s|^{2}(1+\rho)+\sigma^{2}}}\right)\left(r_{\mathsf{R},l}\sqrt{\frac{|s|^{2}\rho}{|s|^{2}(1+\rho)+\sigma^{2}}}\right)$,
for $i,l=1,\ldots,N_{R}$, $i\neq l$, $i\neq k$, and $l\neq k$. Invoke the
results in [30] again to obtain $\mathcal{P}(\mathbf{W}_{\mathsf{x},k})$ as
given by (97), and similarly obtain $\mathcal{P}(\mathbf{W}_{\mathsf{y},k})$,
$k=1,\ldots,N_{R}$. ∎

##  Appendix E Proof of Theorem 4

###### Proof:

Recall (99) with $\mathbf{D}_{\mathsf{I}}=\mathbf{0}$, so that (60) holds.
Also recall (102) and (103)–(104). Following similar steps to those leading to
(116) in Appendix C, we have

$\displaystyle\hat{h}_{\rm
MMSE}=\frac{\mathbf{s}^{\mathrm{T}}\mathbf{D}_{\mathsf{R}}}{(2\pi)^{\frac{\tau}{2}}}\left(\frac{\bm{\Lambda}_{\mathsf{R}}\mathbf{W}_{\mathsf{A}}\mathbf{g}_{\rm
A}}{\left|\mathbf{W}_{\mathsf{A}}\right|^{\frac{1}{2}}\mathcal{P}(\mathbf{W}_{\mathsf{A}})}+\frac{j\bm{\Lambda}_{\mathsf{I}}\mathbf{W}_{\mathsf{B}}\mathbf{g}_{\rm
B}}{\left|\mathbf{W}_{\mathsf{B}}\right|^{\frac{1}{2}}\mathcal{P}(\mathbf{W}_{\mathsf{B}})}\right),$
(123)

where, for $k=1,\ldots,\tau$, the $k$-th element of $\mathbf{g}_{\rm A}$ is
given by

$\displaystyle g_{{\rm A},k}$
$\displaystyle=\displaystyle\int_{\mathbb{R}_{+}^{\tau-1}}e^{-\frac{\widetilde{\mathbf{x}}_{-k}^{\mathrm{T}}\mathbf{W}_{\mathsf{A},k}^{-1}\widetilde{\mathbf{x}}_{-k}}{2}}{\rm
d}\widetilde{\mathbf{x}}_{-k}$ (124)
$\displaystyle=(2\pi)^{\frac{\tau-1}{2}}\left|\mathbf{W}_{\mathsf{A},k}\right|^{\frac{1}{2}}\mathcal{P}(\mathbf{W}_{\mathsf{A},k}),$
(125)

and
$\widetilde{\mathbf{x}}_{-k}=[x_{1}\;\ldots\;x_{k-1}\;x_{k+1}\;\ldots\;x_{\tau}]^{\mathrm{T}}$.
Similarly, we have $\mathbf{g}_{\rm
B}=\left[(2\pi)^{\frac{\tau-1}{2}}\left|\mathbf{W}_{\mathsf{B},k}\right|^{\frac{1}{2}}\mathcal{P}(\mathbf{W}_{\mathsf{B},k})\right]_{k=1}^{\tau}.$
According to [30, p. 192], $\mathcal{P}(\mathbf{W}_{\mathsf{A}})$ is given by
(106), and $\mathcal{P}(\mathbf{W}_{\mathsf{B}})$ is similarly obtained. Using
(101) and (103)–(104), again based on [30, p. 192],
$\mathcal{P}(\mathbf{W}_{\mathsf{A},k})$ is shown to be given by (107) and
$\mathcal{P}(\mathbf{W}_{\mathsf{B},k})$ is similarly derived. From [34, p.
475], for a non-singular $L\times L$ matrix $\mathbf{E}$ and $L$-dimensional
vectors $\mathbf{e}_{a}$ and $\mathbf{e}_{b}$,
$\left|\mathbf{E}+\mathbf{e}_{a}\mathbf{e}_{b}^{\mathrm{T}}\right|=\left|\mathbf{E}\right|\left(1+\mathbf{e}_{b}^{\mathrm{T}}\mathbf{E}^{-1}\mathbf{e}_{a}\right)$.
Apply this identity to obtain
$\left|\mathbf{W}_{\mathsf{A}}\right|=\left|\mathbf{W}_{\mathsf{B}}\right|=\left(\frac{\sigma^{2}}{2}\right)^{\tau}\left(1+\frac{\mathbf{s}^{\mathrm{T}}\mathbf{s}}{\sigma^{2}}\right)$
and
$\left|\mathbf{W}_{\mathsf{A},k}\right|=\left|\mathbf{W}_{\mathsf{B},k}\right|=\left(\frac{\sigma^{2}}{2}\right)^{\tau-1}\left(\frac{\mathbf{s}^{\mathrm{T}}\mathbf{s}+\sigma^{2}}{s_{k}^{2}+\sigma^{2}}\right)$.
Inserting the above results into (123), we obtain (105) after some algebraic
manipulations. ∎

## References

  * [1] E. Larsson, O. Edfors, F. Tufvesson, and T. Marzetta, “Massive MIMO for next generation wireless systems,” _IEEE Commun. Mag._ , vol. 52, no. 2, pp. 186–195, 2014.
  * [2] F. Rusek, et al., “Scaling up MIMO: Opportunities and challenges with very large arrays,” _IEEE Signal Process. Mag._ , vol. 30, no. 1, pp. 40–60, 2013.
  * [3] E. Björnson, J. Hoydis, and L. Sanguinetti, “Massive MIMO networks: Spectral, energy, and hardware efficiency,” _Found. and Trends Signal Process._ , vol. 11, no. 3–4, pp. 154–655, 2017.
  * [4] N. Rajatheva, I. Atzeni, E. Björnson _et al._ , “White paper on broadband connectivity in 6G,” 2020. [Online]. Available: http://jultika.oulu.fi/files/isbn9789526226798.pdf
  * [5] R. Walden, “Analog-to-digital converter survey and analysis,” _IEEE J. Sel. Areas Commun._ , vol. 17, no. 4, pp. 539–550, 1999.
  * [6] J. Mo and R. Heath, Jr., “Capacity analysis of one-bit quantized MIMO systems with transmitter channel state information,” _IEEE Trans. Signal Process._ , vol. 63, no. 20, pp. 5498–5512, 2015.
  * [7] 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.
  * [8] I. Atzeni, A. Tölli, and G. Durisi, “Low-resolution massive MIMO under hardware power consumption constraints,” in _Proc. Asilomar Conf. Signals, Syst., and Comput. (ASILOMAR)_ , 2021.
  * [9] C. Mollén, J. Choi, E. G. Larsson, and R. W. Heath, “Uplink performance of wideband massive MIMO with one-bit ADCs,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 1, pp. 87–100, 2017.
  * [10] 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.
  * [11] J. Mo, P. Schniter, and R. W. Heath, Jr., “Channel estimation in broadband millimeter wave MIMO systems with few-bit ADCs,” _IEEE Trans. Signal Process._ , vol. 66, no. 5, pp. 1141–1154, 2018.
  * [12] N. Bernardo, J. Zhu, Y. C. Eldar, and J. Evans, “Capacity bounds for one-bit MIMO Gaussian channels with analog combining,” _IEEE Trans. Commun._ , vol. 70, no. 11, pp. 7224–7239, 2022.
  * [13] M. Shao and W.-K. Ma, “Divide and conquer: One-bit MIMO-OFDM detection by inexact expectation maximization,” in _Proc. IEEE Int. Conf. Acoust., Speech, and Signal Process. (ICASSP)_ , 2021.
  * [14] J. Choi, J. Mo, and R. W. Heath Jr., “Near maximum-likelihood detector and channel estimator for uplink multiuser massive MIMO systems with one-bit ADCs,” _IEEE Trans. Commun._ , vol. 64, no. 5, pp. 2005–2018, 2016.
  * [15] S. Rao, A. Mezghani, and A. L. Swindlehurst, “Channel estimation in one-bit massive MIMO systems: Angular versus unstructured models,” _IEEE J. Sel. Topics Signal Process._ , vol. 13, no. 5, pp. 1017–1031, 2019.
  * [16] M. S. Stein, S. Bar, J. A. Nossek, and J. Tabrikian, “Performance analysis for channel estimation with 1-bit ADC and unknown quantization threshold,” _IEEE Trans. Signal Process._ , vol. 66, no. 10, pp. 2557–2571, 2018.
  * [17] I. Atzeni and A. Tölli, “Channel estimation and data detection analysis of massive MIMO with 1-bit ADCs,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 6, pp. 3850–3867, 2022.
  * [18] I. Atzeni and A. Tölli, “Uplink data detection analysis of 1-bit quantized massive MIMO,” in _Proc. IEEE Int. Workshop Signal Process. Adv. in Wireless Commun. (SPAWC)_ , 2021.
  * [19] 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.
  * [20] J. J. Bussgang, “Crosscorrelation functions of amplitude-distorted Gaussian signals,” _Res. Lab. Electron., MIT, Cambridge, MA, U.S.A._ , 1952.
  * [21] 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.
  * [22] A. K. Saxena, I. Fijalkow, and A. L. Swindlehurst, “Analysis of one-bit quantized precoding for the multiuser massive MIMO downlink,” _IEEE Trans. Signal Process._ , vol. 65, no. 17, pp. 4624–4634, 2017.
  * [23] S. Jacobsson, G. Durisi, M. Coldrey, and C. Studer, “Linear precoding with low-resolution DACs for massive MU-MIMO-OFDM downlink,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 3, pp. 1595–1609, 2019.
  * [24] S. M. Kay, _Fundamentals of Statistical Signal Processing: Estimation Theory_. Prentice-Hall, 1993.
  * [25] R. E. Curry, _Estimation and Control with Quantized Measurements_. The MIT Press, 1970.
  * [26] I. G. Abrahamson, “Orthant probabilities of the quadrivariate normal distribution,” _Ann. Math. Statist._ , vol. 35, no. 4, pp. 1685–1703, 1964.
  * [27] D. R. Childs, “Reduction of the multivariate normal integral to characteristic form,” _Biometrika_ , vol. 54, no. 1–2, pp. 293–300, 1967.
  * [28] R. L. Plackett, “A reduction formula for normal multivariate integrals,” _Biometrika_ , vol. 41, no. 3–4, pp. 351–360, 1954.
  * [29] H.-J. Sun, “A general reduction method for $n$-variate normal orthant probability,” _Comm. Statist. Theory Methods_ , vol. 17, no. 11, pp. 3913–3921, 1988.
  * [30] Y. L. Tong, _The Multivariate Normal Distribution_. Springer-Verlag, 1990.
  * [31] A. Genz, “Numerical computation of multivariate normal probabilities,” _J. Comput. Graphical Statist._ , vol. 1, no. 2, pp. 141–149, 1992.
  * [32] D. M. Keenan, “A time series analysis of binary data,” _J. Amer. Statist. Assoc._ , vol. 77, no. 380, pp. 816–821, 1982.
  * [33] M. Sinn and K. Keller, “Covariances of zero crossings in Gaussian processes,” _Theory Probab. Appl._ , vol. 55, no. 3, pp. 485–504, 2011.
  * [34] C. D. Meyer, _Matrix Analysis and Applied Linear Algebra_. SIAM, 2000.
  * [35] R. A. Johnson and D. W. Wichern, _Applied Multivariate Statistical Analysis_. Pearson Education Inc., 2007.
  * [36] V. Aalo, “Performance of maximal-ratio diversity systems in a correlated Nakagami-fading environment,” _IEEE Trans. Commun._ , vol. 43, no. 8, pp. 2360–2369, 1995.
  * [37] M. Chiani, M. Z. Win, and A. Zanella, “On the capacity of spatially correlated MIMO Rayleigh-fading channels,” _IEEE Trans. Inf. Theory_ , vol. 49, no. 10, pp. 2363–2371, 2003.
  * [38] A. Abdi and M. Kaveh, “A space-time correlation model for multielement antenna systems in mobile fading channels,” _IEEE J. Sel. Areas Commun._ , vol. 20, no. 3, pp. 550–560, 2002.
  * [39] W. C. Jakes, _Microwave Mobile Communications_. Wiley-IEEE Press, 1974.
  * [40] R. M. Gray, “Toeplitz and Circulant Matrices: A Review,” _Found. and Trends Commun. and Inf. Theory_ , vol. 2, no. 3, pp. 155–239, 2006.