BEACHES: Beamspace Channel Estimation for Multi-Antenna mmWave Systems and Beyond

TL;DR

BEACHES is a low-complexity, beamspace-based channel estimation algorithm for multi-antenna mmWave systems that denoises sparse channel information leveraging directional wave propagation, with an efficient tuning stage for minimal mean-square error.

Contribution

It introduces BEACHES, a novel beamspace channel estimation method that efficiently denoises sparse channels and includes a tuning stage to optimize accuracy in large-antenna systems.

Findings

Effective in LoS and non-LoS scenarios

Reduces complexity compared to traditional methods

Achieves minimal mean-square error in large-antenna limit

Abstract

Massive multi-antenna millimeter wave (mmWave) and terahertz wireless systems promise high-bandwidth communication to multiple user equipments in the same time-frequency resource. The high path loss of wave propagation at such frequencies and the fine-grained nature of beamforming with massive antenna arrays necessitates accurate channel estimation to fully exploit the advantages of such systems. In this paper, we propose BEAmspace CHannel EStimation (BEACHES), a low-complexity channel estimation algorithm for multi-antenna mmWave systems and beyond. BEACHES leverages the fact that wave propagation at high frequencies is directional, which enables us to denoise the (approximately) sparse channel state information in the beamspace domain. To avoid tedious parameter selection, BEACHES includes a computationally-efficient tuning stage that provably minimizes the mean-square error of the channel estimate in the large-antenna limit. To demonstrate the efficacy of BEACHES, we provide simulation results for line-of-sight (LoS) and non-LoS mmWave channel models.

Tables

Table 1. TABLE I: MATLAB runtimes in milliseconds (and normalized runtimes).

Scenario	BEACHES	NOMP	ANM
LoS	1.64 (1$\times$)	199.9 (120$\times$)	47 968 (29 000$\times$)
non-LoS	1.45 (1$\times$)	2 204 (1 500$\times$)	83 750 (58 000$\times$)

Equations

$$\displaystyle\mathbf{h}=\sum_{\ell=0}^{L-1}\alpha_{\ell}\mathbf{a}(\Omega_{\ell}),\,\,\mathbf{a}(\Omega)\!=\!\big{[}e^{j0\Omega},e^{j1\Omega},\ldots,e^{j(B-1)\Omega}\big{]}^{\textnormal{T}}\!\!.$$

$${\hat{\mathbf{h}}}^{\star }=\underset{{\hat{\mathbf{h}}}^{\prime }\in {\mathbb{C}}^B}{\mathrm{argmin}}\frac{1}{2}\Vert \hat{\mathbf{y}}-{\hat{\mathbf{h}}}^{\prime }{\Vert }_2^2+\tau \Vert {\hat{\mathbf{h}}}^{\prime }{\Vert }_1.$$

$$[\eta (\hat{\mathbf{y}},\tau ){]}_b=\frac{{\hat{y}}_b}{|{\hat{y}}_b|}\max \{|{\hat{y}}_b|-\tau ,0\},b=1,\dots ,B,$$

$$\text{MSE}=\frac{1}{B}E\left[\Vert \hat{}{\mathbf{h}}^{\star }-\hat{\mathbf{h}}{\Vert }_2^2\right],$$

$$\text{SURE}=$$

$$+\frac{1}{B}{E}_0\sum _{b=1}^B\left(\frac{\partial [{\mu }_{\mathcal{R}}(\hat{\mathbf{y}}){]}_b}{\partial [{\hat{\mathbf{y}}}_{\mathcal{R}}{]}_b}+\frac{\partial [{\mu }_{\mathcal{I}}(\hat{\mathbf{y}}){]}_b}{\partial [{\hat{\mathbf{y}}}_{\mathcal{I}}{]}_b}-2\right)$$

$${\text{SURE}}_{\tau }=$$

$$-\frac{E_0}{B}\tau \sum _{b:|{\hat{y}}_b|>\tau }\frac{1}{|{\hat{y}}_b|}-2\frac{E_0}{B}\sum _{b:|{\hat{y}}_b|<\tau }1,$$

$$\underset{B\to \infty }{\lim }{\text{SURE}}_{\tau }=\text{MSE}.$$

$${\text{SURE}}_{\tau }=$$

$$-\frac{E_0}{B}\tau \sum _{b=k}^{B+1}({\hat{y}}_b^s{)}^{-1}-2\frac{E_0}{B}(k-1),$$

$$\displaystyle\tau^{\star}_{k}=\left\{\begin{array}[]{ll}\,\tau_{k}^{Q}&\hat{y}^{s}_{k-1}<\tau_{k}^{Q}<\hat{y}^{s}_{k}\\ \,\hat{y}^{s}_{k-1}&\tau_{k}^{Q}<\hat{y}^{s}_{k-1}\\ \,\hat{y}^{s}_{k}&\tau_{k}^{Q}>\hat{y}^{s}_{k},\end{array}\right.$$

$$\text{MSE}=E\left[\frac{1}{B}\Vert {\hat{\mathbf{h}}}^{\star }-\hat{\mathbf{h}}{\Vert }_2^2\right]=E\left[\frac{1}{B}\Vert \mu (\hat{\mathbf{y}})-\hat{\mathbf{h}}{\Vert }_2^2\right],$$

$$\text{MSE}=$$

$$=$$

$$+E\left[\frac{2}{B}{\left[g(\hat{\mathbf{y}}{)}^{\text{H}}(\hat{\mathbf{y}}-\hat{\mathbf{h}})\right]}_{\mathcal{R}}\right].$$

$$\frac{2}{B}E\left[{\left[g(\hat{\mathbf{y}}{)}^{\text{H}}(\hat{\mathbf{y}}-\hat{\mathbf{h}})\right]}_{\mathcal{R}}\right]=$$

$$+\frac{2}{B}E\left[g_{\mathcal{I}}(\hat{\mathbf{y}}{)}^{\text{T}}({\hat{\mathbf{y}}}_{\mathcal{I}}-{\hat{\mathbf{h}}}_{\mathcal{I}})\right].$$

$$\frac{2}{B}$$

$$\stackrel{\text{(a)}}{=}$$

$$\exp \left(-\frac{\Vert {\hat{\mathbf{y}}}_{\mathcal{R}}-{\hat{\mathbf{h}}}_{\mathcal{R}}{\Vert }^2}{2\frac{E_0}{2}}\right)\frac{E_0}{2}\frac{\partial [g_{\mathcal{R}}(\hat{\mathbf{y}}){]}_b}{\partial [{\hat{\mathbf{y}}}_{\mathcal{R}}{]}_b}d{\hat{\mathbf{y}}}_{\mathcal{R}}d{\hat{\mathbf{y}}}_{\mathcal{I}}$$

$$=$$

$$\frac{2}{B}E\left[g_{\mathcal{I}}(\hat{\mathbf{y}}{)}^{\text{T}}({\hat{\mathbf{y}}}_{\mathcal{I}}-{\hat{\mathbf{h}}}_{\mathcal{I}})\right]=\frac{E_0}{B}E\left[{\sum }_{b=1}^B\left(\frac{\partial [g_{\mathcal{I}}(\hat{\mathbf{y}}){]}_b}{\partial [{\hat{\mathbf{y}}}_{\mathcal{I}}{]}_b}-1\right)\right].$$

$$\text{MSE}=$$

$$+\frac{E_0}{B}E\left[{\sum }_{b=1}^B\left(\frac{\partial [{\mu }_{\mathcal{R}}(\hat{\mathbf{y}}){]}_b}{\partial [{\hat{\mathbf{y}}}_{\mathcal{R}}{]}_b}+\frac{\partial [{\mu }_{\mathcal{I}}(\hat{\mathbf{y}}){]}_b}{\partial [{\hat{\mathbf{y}}}_{\mathcal{I}}{]}_b}-2\right)\right].$$

$$\text{SURE}=$$

$$+\frac{E_0}{B}{\sum }_{b=1}^B\left(\frac{\partial [{\mu }_{\mathcal{R}}(\hat{\mathbf{y}}){]}_b}{\partial [{\hat{\mathbf{y}}}_{\mathcal{R}}{]}_b}+\frac{\partial [{\mu }_{\mathcal{I}}(\hat{\mathbf{y}}){]}_b}{\partial [{\hat{\mathbf{y}}}_{\mathcal{I}}{]}_b}-2\right),$$

$$\frac{\partial [{\eta }_{\mathcal{R}}(\hat{\mathbf{y}},\tau ){]}_b}{\partial [{\hat{\mathbf{y}}}_{\mathcal{R}}{]}_b}=\frac{\partial [{\eta }_{\mathcal{I}}(\hat{\mathbf{y}},\tau ){]}_b}{\partial [{\hat{\mathbf{y}}}_{\mathcal{I}}{]}_b}=0.$$

$$\frac{\partial [{\eta }_{\mathcal{R}}(\hat{\mathbf{y}},\tau ){]}_b}{\partial [{\hat{\mathbf{y}}}_{\mathcal{R}}{]}_b}=$$

$$=$$

$$\frac{\partial [{\eta }_{\mathcal{I}}(\hat{\mathbf{y}},\tau ){]}_b}{\partial [{\hat{\mathbf{y}}}_{\mathcal{I}}{]}_b}=$$

$${\text{SURE}}_{\tau }=$$

$$\displaystyle\,\textstyle+\frac{E_{0}}{B}\sum_{b:|\hat{y}_{b}|>\tau}\Big{(}2-\tau\frac{1}{\sqrt{[\hat{\mathbf{y}}_{\mathcal{R}}]_{b}^{2}+[\hat{\mathbf{y}}_{\mathcal{I}}]_{b}^{2}}}-2\Big{)}$$

$$+\frac{E_0}{B}{\sum }_{b:|{\hat{y}}_b|<\tau }\left(0-2\right).$$

$$\underset{B\to \infty }{\lim }\frac{1}{B}{\sum }_{b=1}^B\gamma (\eta ({\hat{y}}_b,\tau ),{\hat{h}}_b)$$

$$=E\left[\gamma (\eta (H+\sqrt{E_0}Z,\tau ),H)\right].$$

$$\underset{B\to \infty }{\lim }\frac{1}{B}{\sum }_{b=1}^B|\eta ({\hat{y}}_b,\tau )-{\hat{y}}_b{|}^2=E_{{\hat{y}}_{\tilde{b}}}\left[|\eta ({\hat{y}}_{\tilde{b}},\tau )-{\hat{y}}_{\tilde{b}}{|}^2\right],$$

$$\underset{B\to \infty }{\lim }\frac{1}{B}\Vert \eta ({\hat{y}}_b,\tau )-{\hat{y}}_b{\Vert }_2^2=E_{\hat{\mathbf{y}}}\left[\frac{1}{B}\Vert \eta (\hat{\mathbf{y}},\tau )-\hat{\mathbf{y}}{\Vert }_2^2\right].$$

Full text

I Introduction

Massive multiuser (MU) multiple-input multiple-output (MIMO) [1] and millimeter-wave (mmWave) communication [2, 3] are among the key technologies of next-generation wireless systems. The high path loss of wave propagation at mmWave or terahertz (THz) frequencies and the fact that massive MU-MIMO enables fine-grained beamforming, requires the basestations (BSs) to acquire accurate channel state information (CSI) [4, 5]. In addition, the trend towards low-precision data converters in all-digital massive MU-MIMO BSs to reduce power, interconnect bandwidth, and costs [6] renders accurate channel estimation increasingly important.

At mmWave or terahertz frequencies, wave propagation is highly directional and real-world channels typically comprise only a small number of dominant propagation paths [2, 3]. These unique properties enable the deployment of channel estimation algorithms that effectively suppress noise [7, 8, 9]. As a consequence, compressive sensing (CS)-based methods have been proposed for mmWave channel estimation in [10, 11]. Most of such methods use a discretization procedure of the number of propagation paths that can be resolved in the beamspace domain [12], resulting in the well-known basis mismatch problem [13]. Methods that perform off-grid CS, such as atomic norm minimization (ANM) [14] or Newtonized orthogonal matching pursuit (NOMP) [15], avoid this basis mismatch problem. These methods, however, exhibit excessively high complexity, which renders the design of corresponding hardware designs challenging.

I-A Contributions

We propose a new channel estimation algorithm for massive MU-MIMO mmWave/terahertz communication systems that relies on Stein’s unbiased risk estimator (SURE). Our algorithm, called BEAmspace CHannel EStimation (BEACHES), exploits sparsity of mmWave/terahertz channels in the beamspace domain and adaptively denoises the channel vectors at low complexity. We prove that BEACHES minimizes the mean square error (MSE) of the channel estimate in the large-antenna limit. We evaluate BEACHES for LoS and non-LoS mmWave channels and demonstrate that it performs on par with ANM and NOMP but at orders-of-magnitude lower complexity.

I-B Notation

Lowercase and uppercase boldface letters designate column vectors and matrices, respectively. For a vector $\mathbf{a}$, the $k$th entry is $[\mathbf{a}]_{k}=a_{k}$; the real and imaginary parts are $[\mathbf{a}]_{\mathcal{R}}=\mathbf{a}_{\mathcal{R}}$ and $[\mathbf{a}]_{\mathcal{I}}=\mathbf{a}_{\mathcal{I}}$, respectively. The transpose and conjugate transpose of matrix $\mathbf{A}$ are $\mathbf{A}^{\textnormal{T}}$ and $\mathbf{A}^{\textnormal{H}}$, respectively. A complex Gaussian vector $\mathbf{a}$ with mean vector $\mathbf{m}$ and covariance matrix $\mathbf{K}$ is written as $\mathbf{a}\sim\mathcal{C}\mathcal{N}(\mathbf{m},\mathbf{K})$ and its probability density function (PDF) as $f^{\mathcal{C}\mathcal{N}}(\mathbf{a};\mathbf{m},\mathbf{K})$. A real Gaussian vector $\mathbf{a}$ with mean vector $\mathbf{m}$ and covariance matrix $\mathbf{K}$ is written as $\mathbf{a}\sim\mathcal{N}(\mathbf{m},\mathbf{K})$ and its PDF as $f^{\mathcal{N}}(\mathbf{a};\mathbf{m},\mathbf{K})$. The expectation operator is $\operatorname{\mathbb{E}}\mathopen{}\left[\cdot\right]$.

II System Model and Beamspace Representation

II-A System Model

We consider an all-digital mmWave/THz massive MU-MIMO uplink system as illustrated in \freffig:system_overview. The BS is equipped with a $B$-antenna uniform linear array (ULA) and communicates with $U$ single-antenna user equipments (UEs) in the same time-frequency resource. For simplicity, we focus on pilot-based channel estimation for flat-fading channels, where the BS estimates the $B$-dimensional complex channel vector $\mathbf{h}\in\mathbb{C}^{B}$ for each UE. By assuming that (i) wave propagation is predominantly directional [16, 4], and (ii) the distance between UEs and BS is sufficiently large, the channel vectors in the antenna domain can be modeled as follows [17]:

[TABLE]

Here, $L$ refers to the total number of paths arriving at the antenna array (including a potential LoS path), $\alpha_{\ell}\in\mathbb{C}$ is the complex-valued channel gain of the $\ell$th path, and $\mathbf{a}(\Omega_{\ell})$ represents a complex-valued sinusoid containing the relative phases between BS antennas, where $\Omega_{\ell}\in[0,2\pi)$ is determined by the incident angle of the $\ell$th path to the antenna array. We model the estimated channel vector in the antenna domain as $\mathbf{y}=\mathbf{h}+\mathbf{e}$, where $\mathbf{e}\sim\mathcal{C}\mathcal{N}(\mathbf{0}_{B\times 1},E_{0}\mathbf{I}_{B})$ represents channel estimation error with variance $E_{0}$ per complex entry.

II-B Beamspace Channel Vector Denoising

The channel vectors $\mathbf{h}$ as modeled in \frefeq:planewavemodel are a superposition of $L$ complex-valued sinusoids. Hence, it is useful to transform the vector $\mathbf{h}$ into the discrete Fourier transform (DFT) domain, also known as the beamspace domain, $\hat{\mathbf{h}}=\mathbf{F}\mathbf{h}$, where $\mathbf{F}$ is the $B\times B$ unitary DFT matrix. In the beamspace domain, each entry of $\hat{}\mathbf{h}$ is associated to a specific incident angle with respect to the BS antenna array [12]. If the number of paths $L$ is smaller than the number of BS antennas $B$, then the beamspace channel vector $\hat{}\mathbf{h}$ will be (approximately) sparse [8]. This key property enables the use of denoising algorithms. More specifically, by transforming $\mathbf{y}$ into the beamspace domain $\hat{\mathbf{y}}=\mathbf{F}\mathbf{y}=\hat{\mathbf{h}}+\hat{\mathbf{e}}$, where $\hat{\mathbf{e}}=\mathbf{F}\mathbf{e}$ has the same statistics as $\mathbf{e}$, one can suppress noise while preserving the strong beamspace components. Prominent methods for beamspace denoising are ANM [14] and NOMP [15], which require high complexity.

III BEACHES: BEAmspace CHannel EStimation

III-A Channel Vector Denoising via Soft-Thresholding

A widely-used sparsity-based denoising method is the least absolute shrinkage and selection operator (LASSO) [18, 19], which corresponds to the following optimization problem:

[TABLE]

Here, $\tau\in\mathbb{R}_{+}$ is a suitably-chosen denoising parameter. The solution to \frefeq:LASSO in the complex case is the well-known soft-thresholding operator $\eta(\hat{\mathbf{y}},\tau)$ defined as [20, App. A]

[TABLE]

where we define $y/|y|=0$ for $y=0$. For sparsity-based denoising via soft-thresholding, the performance strongly depends on the choice of the denoising parameter $\tau$ [18, 21]. In wireless systems, it is particularly important to design robust methods to select this parameter, as many factors such as the propagation conditions, the number of arriving paths, and the signal power can vary widely over time.

III-B Computing the Optimal Denoising Parameter

In what follows, we are interested in the optimal parameter $\tau^{\star}$ that minimizes the estimation MSE defined as

[TABLE]

where $\hat{}\mathbf{h}^{\star}=\eta(\hat{\mathbf{y}},\tau^{\star})$ is the associated denoised beamspace vector. Determining the optimal parameter $\tau^{\star}$ requires knowledge of the noiseless beamspace vector $\hat{\mathbf{h}}$, which is unknown. To avoid the need for knowing the ground truth $\hat{\mathbf{h}}$, we propose to use Stein’s unbiased risk estimate (SURE) [18] as a proxy for the MSE function. The following result provides SURE in the complex domain and shows that it is an unbiased estimator for the MSE. The proof is given in \frefapp:MSEtoSURE.

Theorem 1**.**

Let $\hat{}\mathbf{h}\in\mathbb{C}^{B}$ be an unknown vector and $\hat{\mathbf{y}}\in\mathbb{C}^{B}$ a noisy observation vector distributed as $\hat{}\mathbf{y}\sim\mathcal{C}\mathcal{N}(\hat{}\mathbf{h},E_{0}\mathbf{I}_{B})$. Let $\mu(\hat{\mathbf{y}})$ be an estimator of $\hat{}\mathbf{h}$ from $\hat{\mathbf{y}}$ that is weakly differentiable and operates element-wise on vectors. Then,

[TABLE]

is an unbiased estimate of the MSE, i.e., $\operatorname{\mathbb{E}}\mathopen{}\left[\textit{SURE}\right]=\textit{MSE}$.

The following theorem shows that SURE for the soft-thresholding operator $\eta(\hat{\mathbf{y}},\tau)$ converges to the MSE in the large-antenna limit $B\to\infty$. The proof is given in \frefapp:shrinkageSURE.

Theorem 2**.**

For the soft-thresholding function $\mu(\hat{\mathbf{y}})=\eta(\hat{\mathbf{y}},\tau)$ in \frefeq:softthresholdingfunction, SURE in \frefeq:complexSURE is given by111As discussed in \frefapp:shrinkageSURE, the value of SURE is undefined for $\tau=\hat{y}_{b}$, $b=1,\ldots,B$, due to the non-differentiability of the function $\eta$.

[TABLE]

which, in the limit $B\to\infty$ converges to the MSE, i.e.,

[TABLE]

SURE in \frefeq:shrinkageSURE is independent of the true beamspace channel vector $\hat{\mathbf{h}}$. The expression \frefeq:shrinkageSURE only depends on the magnitudes of the observed beamspace channel vector $\hat{\mathbf{y}}$, the channel estimation error variance $E_{0}$, the number of BS antennas $B$, and the denoising parameter $\tau$. Thanks to \frefeq:SURE_convergence and the fact that $B$ is large in massive MU-MIMO systems, we can use SURE as a surrogate to minimize the MSE and to determine the optimal denoising parameter. While no closed-form expression for the minimum of \frefeq:shrinkageSURE is known, reference [21] proposes a bisection procedure to approximate the optimal value of a similar expression for sparse recovery. We next propose an efficient algorithm that computes the optimal parameter $\tau^{\star}$ using a deterministic procedure with complexity $O(B\log(B))$.

III-C The BEACHES Algorithm

Reference [18] outlines an efficient procedure to minimize SURE for real-valued wavelet denoising. We propose a similar strategy to minimize \frefeq:shrinkageSURE for the complex-valued case with soft-thresholding. Instead of continuously sweeping the denoising parameter $\tau$ in the interval $[0,\infty)$, we first sort the absolute values of the vector $\hat{}\mathbf{y}$ in ascending order which we call $\hat{\mathbf{y}}^{s}$. We then search for the optimal denoising parameter $\tau$ only between each pair of consecutive elements of the sorted vector, i.e., $\tau\in\left(\hat{y}^{s}_{k-1},\hat{y}^{s}_{k}\right)$ for $k=1,\ldots,B+1$. In each such interval, SURE is a quadratic function in $\tau$ given by

[TABLE]

where we define $\hat{y}^{s}_{0}=0$ and $\hat{y}^{s}_{B+1}=\infty$ to account for the first $\left(0,\hat{y}^{s}_{1}\right)$ and last interval $\left(\hat{y}^{s}_{B},\infty\right)$. For each $k\in\{1,\ldots,B+1\}$, we compute the optimal value of $\tau$ that minimizes SURE in each interval $\tau\in\left(\hat{y}^{s}_{k-1},\hat{y}^{s}_{k}\right)$. Since there is a discontinuity in the SURE expression when progressing from one interval to the next, the minimal value in each interval is either the minimum of the quadratic function \frefeq:QuadraticSURE or one of the boundaries of the interval $\hat{y}^{s}_{k-1}$ and $\hat{y}^{s}_{k}$.222Note that SURE is not defined for $\tau=\hat{y}^{s}_{k-1}$ and $\tau=\hat{y}^{s}_{k}$. Instead, we compute $\textit{SURE}_{\tau}$ for two values arbitrarily close to these boundaries, i.e., $\tau=\hat{y}^{s}_{k-1}+\epsilon$ and $\tau=\hat{y}^{s}_{k}-\epsilon$, where $\epsilon>0$ is small compared to $\tau$. The minimum value of \frefeq:QuadraticSURE is given by $\tau_{k}^{Q}=\textstyle\frac{E_{0}}{2(B-k+1)}{\sum_{b=k}^{B+1}(\hat{y}^{s}_{b})^{-1}}$. Since the function $\text{SURE}_{\tau}$ is convex, we can identify the optimal parameter $\tau$ in the interval $\left(\hat{y}^{s}_{k-1},\hat{y}^{s}_{k}\right)$ indexed by $k$ by knowing the value of $\tau_{k}^{Q}$ with respect to the interval boundaries $\hat{y}^{s}_{k-1}$ and $\hat{y}^{s}_{k}$. Put simply, the optimal denoising parameter $\tau^{\star}_{k}$ in each interval $k=1,\ldots,B+1$ is given by

[TABLE]

or simply $\tau^{\star}_{k}=\max\{\hat{y}^{s}_{k-1},\min\{\hat{y}^{s}_{k},\tau_{k}^{Q}\}\}$. By knowing the optimal value of $\tau$ in each interval, we only need to find the minimal value of $\text{SURE}_{\tau^{\star}_{k}}$ for $k=1,\ldots,B+1$. It is now key to realize that we do not need to recalculate SURE in \frefeq:QuadraticSURE from scratch while searching through $k=1,\ldots,B+1$. Instead, we can sequentially update the quantities $S=\sum_{b=0}^{k-1}{(\hat{y}^{s}_{b})^{2}}$ and $V=\sum_{b=k}^{B+1}{(\hat{y}^{s}_{b})^{-1}}$, thanks to sorting the magnitudes of the vector $\hat{\mathbf{y}}$. The resulting procedure, called BEACHES, is summarized in \frefalg:BEACHES. The computational complexity of BEACHES is only ${O}(B\log(B))$, which is caused by the FFT, sorting, and IFFT operations—the remaining computations in the for-loop (lines 5–13) are simple scalar operations.

IV Performance and Runtime of BEACHES

IV-A Bit Error-Rate Performance

To assess the performance of BEACHES, we consider an all-digital massive MU-MIMO system in which $U=16$ UEs communicate with a $B=256$ antenna BS. We focus on the situation in which the UEs first send orthogonal pilots, which are used to acquire maximum-likelihood channel estimates $\mathbf{y}_{u}$ for each UE $u=1,\ldots,U$. The channel matrices are generated for both LoS and non-LoS conditions using the QuaDRiGa mmMAGIC UMi model [22] at a carrier frequency of 60 GHz with a ULA using $\lambda/2$ antenna spacing. The UEs are placed randomly within a $120^{\circ}$ circular sector with minimum and maximum distance of $10$ and $110$ meters from the BS antenna array, respectively. We enforce UE separation of at least $1^{\circ}$ (relative to the BS antenna array) and assume optimal UE power control. We then use different channel vector denoising methods, including (i) ANM-based denoising, where we use the debiased output of the code provided in [14], (ii) NOMP with a (manually tuned) false alarm rate of $P_{\text{fa}}=0.5$ using the code provided in [15], and (iii) “perfect CSI,” which is a baseline that uses the noiseless channel vectors. Finally, we transmit 16-QAM symbols and perform linear minimum MSE (L-MMSE) equalization with the denoised matrix to detect the transmitted bits. The resulting uncoded bit error-rate (BER) is used to assess the performance of various denoising methods.

Figure 2 shows that channel vector denoising in the beamspace domain provides $2$ dB to $3$ dB SNR performance improvements at $\textit{BER}=10^{-2}$ compared to conventional ML channel estimation. The achieved performance gains are more pronounced under LoS conditions. Quite surprisingly, we observe that BEACHES performs on par to ANM and NOMP. This observation indicates that off-the-grid denoising methods, such as ANM and NOMP, do not provide a critical performance advantage over BEACHES (in terms of BER).

IV-B Runtime Comparison

While the BER performance of BEACHES is comparable to ANM and NOMP, it exhibits significantly lower complexity. To support this claim, we measured their MATLAB runtimes in milliseconds on an Intel core i5-7400 CPU with 16 GB RAM at an SNR of $5$ dB. \freftbl:runtimes demonstrates that the runtime of BEACHES is orders of magnitude lower than that of NOMP (up to $1\,500\times$) and ANM (up to $58\,000\times$), while the speedup is more pronounced for the non-LoS scenario.

V Conclusions

We have proposed a new channel denoising algorithm for massive MU-MIMO mmWave and terahertz communication systems called BEAmspace CHannel EStimation (BEACHES). BEACHES exploits sparsity of mmWave/terahertz channels in the beamspace domain to perform adaptive soft-thresholding via Stein’s unbiased risk estimate (SURE). We have shown that BEACHES minimizes the mean square error in the large-antenna limit and performs on par with sophisticated channel estimation algorithms for realistic LoS and non-LoS channel models but at orders-of-magnitude lower complexity. There are many avenues for future work. An extension of BEACHES to systems with low-precision quantizers and single-carrier transmission is an open research problem.

Appendix A Proof of \frefthm:MSEappox

The MSE for $\hat{\mathbf{h}}^{\star}=\mu(\hat{\mathbf{y}})$ is defined as

[TABLE]

where we decompose the complex-valued vector $\hat{\mathbf{y}}$ into the real part $\hat{\mathbf{y}}_{\mathcal{R}}\sim\mathcal{N}(\mathbf{h}_{\mathcal{R}},\textstyle\frac{E_{0}}{2}\mathbf{I}_{B})$ and imaginary part $\hat{\mathbf{y}}_{\mathcal{I}}\sim\mathcal{N}(\mathbf{h}_{\mathcal{I}},\textstyle\frac{E_{0}}{2}\mathbf{I}_{B})$. Note that expectation is with respect to the noisy observation $\hat{\mathbf{y}}$. Define $g(\hat{\mathbf{y}})=\mu(\hat{\mathbf{y}})-\hat{\mathbf{y}}$. Hence,

[TABLE]

The last term can be simplified as

[TABLE]

We can now expand $\frac{2}{B}\operatorname{\mathbb{E}}\mathopen{}\left[{g_{\mathcal{R}}(\hat{\mathbf{y}})^{\textnormal{T}}(\hat{\mathbf{y}}_{\mathcal{R}}-\hat{\mathbf{h}}_{\mathcal{R}})}\right]$ which yields

[TABLE]

where (a) follows from integration by parts. Similarly, we have

[TABLE]

Recall that $g(\hat{\mathbf{y}})=\mu(\hat{\mathbf{y}})-\hat{\mathbf{y}}$ and replace \frefeq:ReRe and \frefeq:ImIm in the original MSE expression in \frefeq:MSE_complex. This leads to

[TABLE]

The second term in the MSE expression above equals $E_{0}$. For the first and third term we remove their expectations to arrive at the following SURE expression:

[TABLE]

which establishes the fact that $\operatorname{\mathbb{E}}\mathopen{}\left[\textit{SURE}\right]=\textit{MSE}$.

Appendix B Proof of \frefthm:SURE_shrinkage

SURE in \frefeq:complexSURE for $\mu(\hat{\mathbf{y}})=\eta(\hat{\mathbf{y}},\tau)$ is derived as follows. The only unknowns in the expression of SURE are its derivative of real and imaginary parts. For $|\hat{y}_{b}|<\tau$, we have

[TABLE]

For $|\hat{y}_{b}|>\tau$, we have

[TABLE]

and

[TABLE]

Note that at $|\hat{y}_{b}|=\tau$, there is a discontinuity and thus the derivative and consequently SURE are not defined for this value. The complex-valued SURE expression reduces to

[TABLE]

We now prove the convergence of SURE in \frefeq:SURE_convergence. In [23, Lemma 4.3.], the authors prove convergence of SURE to MSE in the real domain for the soft-thresholding function. We follow the same procedure for the complex domain. Using [20, Thm. III.15 & III.16], we have that for any pseudo-Lipschitz function $\gamma:\mathbb{C}\to\mathbb{R}$ the following equality holds:

[TABLE]

Here, $Z\sim\mathcal{C}\mathcal{N}(0,1)$ and $H$ is a random variable with the sparse distribution of the channel vector in the beamspace domain $\hat{h}_{b}$. Using \frefeq:lipschitz_stuff, we have the following result

[TABLE]

where, $\hat{y}_{\tilde{b}}$ is any element of the random vector $\hat{\mathbf{y}}$. The expression above can be rewritten as

[TABLE]

Now, since $\frac{\partial{[\eta_{\mathcal{R}}(\hat{\mathbf{y}},\tau)]_{b}}}{\partial[\hat{\mathbf{y}}_{\mathcal{R}}]_{b}}+\frac{\partial{[\eta_{\mathcal{I}}(\hat{\mathbf{y}},\tau)]_{b}}}{\partial[\hat{\mathbf{y}}_{\mathcal{I}}]_{b}}$ is bounded, it is pseudo-Lipschitz and, hence, we can use \frefeq:lipschitz_stuff to obtain the following convergence result:

[TABLE]

By summing \frefeq:term1 and \frefeq:term3 and $\operatorname{\mathbb{E}}\mathopen{}\left[\frac{1}{B}\|\hat{\mathbf{y}}-\hat{\mathbf{h}}\|_{2}^{2}\right]=E_{0}$, we have established that $\lim_{B\to\infty}\text{SURE}_{\tau}=\text{MSE}$.