# Bayesian Learning for Pilot Decontamination in Cell-Free Massive MIMO

Christian Forsch, Zilu Zhao, Dirk Slock, Laura Cottatellucci

► To cite this version:

Christian Forsch, Zilu Zhao, Dirk Slock, Laura Cottatellucci. Bayesian Learning for Pilot Decontamination in Cell-Free Massive MIMO. WSA 2025, 28th International Workshop on Smart Antennas, Sep 2025, Erlangen, Germany. pp.19-25, 10.1109/WSA65299.2025.11202775. hal-05509066

HAL Id: hal-05509066

https://hal.science/hal-05509066v1

Submitted on 13 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Bayesian Learning for Pilot Decontamination in
Cell-Free Massive MIMO

Christian Forsch*, Zilu Zhao†, Dirk Slock†, and Laura Cottatellucci*

*Institute for Digital Communications, Friedrich-Alexander-Universität Erlangen-Nürnberg, Erlangen, Germany

†Communication Systems Department, EURECOM, Sophia Antipolis, France

Email: {christian.forsch, laura.cottatellucci}@fau.de, {zilu.zhao, dirk.slock}@eurecom.fr

***Abstract***—Pilot contamination (PC) arises when the pilot sequences assigned to user equipments (UEs) are not mutually orthogonal, eventually due to their reuse. In this work, we propose a novel expectation propagation (EP)-based joint channel estimation and data detection (JCD) algorithm specifically designed to mitigate the effects of PC in the uplink of cell-free massive multiple-input multiple-output (CF-MaMIMO) systems. This modified bilinear-EP algorithm is distributed, scalable, demonstrates strong robustness to PC, and outperforms state-of-the-art Bayesian learning algorithms. Through a comprehensive performance evaluation, we assess the performance of Bayesian learning algorithms for different pilot sequences and observe that the use of non-orthogonal pilots can lead to better performance compared to shared orthogonal sequences. Motivated by this analysis, we introduce a new metric to quantify PC at the UE level. We show that the performance of the considered algorithms degrades monotonically with respect to this metric, providing a valuable theoretical and practical tool for understanding and managing PC via iterative JCD algorithms.

***Index Terms***—Cell-free massive MIMO, pilot contamination, joint channel estimation and data detection, expectation propagation, non-orthogonal pilot sequences.

# I. INTRODUCTION

Distributed multiple-input multiple-output (MIMO) communications, particularly in the form of cell-free massive MIMO (CF-MaMIMO) networks, are expected to play a key role in advancing next-generation mobile communication systems by enabling high-rate and energy-efficient communication everywhere in the coverage area [1]–[3]. In this network architecture, a large number of spatially distributed access points (APs) are communicating with a smaller number of user equipments (UEs). One major challenge in practical CF-MaMIMO networks is pilot contamination (PC), which arises from the use of non-orthogonal pilot sequences for channel estimation and deteriorates the overall system performance. In real-world networks, the use of non-orthogonal pilot sequences is necessary due to the potentially large number of UEs in the network. Ensuring orthogonality would require prohibitively long pilot sequences, reducing spectral efficiency and throughput per user. Furthermore, unlike centralized MaMIMO, channel hardening and favorable propagation typically do not hold in CF-MaMIMO [4]–[7], which precludes the use of existing pilot decontamination methods proposed for centralized MaMIMO, e.g., [8]–[12]. Even though PC has been shown not to represent a fundamental limitation for centralized and CF-MaMIMO systems [13], [14], it still

remains a significant practical challenge, especially in scalable CF-MaMIMO systems with a limited number of APs and AP antennas. This motivates the development of new effective, efficient, and distributed pilot decontamination schemes capable of leveraging other degrees of freedom beyond those considered in [13] and [14], which rely primarily on pilot symbols and channel statistics. These additional methods are presented below.

An efficient approach to mitigate PC is through optimized pilot assignment schemes because strong PC arises when UEs in close proximity utilize the same pilot sequence. In [1], the authors propose a greedy pilot assignment scheme based on minimizing PC for the user with the lowest rate. Here, the PC for a given UE is quantified by the average interference power after projecting the received pilot signal onto the pilot sequence used by the given UE. In [15], a K-means and a user-group-based pilot assignment scheme are presented where the distances between all UEs and APs and their serving relationships are utilized to assign the pilot sequences. An alternative strategy to alleviate PC leverages joint channel estimation and data detection (JCD). The authors in [16] developed a JCD scheme based on forward backward splitting which exploits the sparsity of CF-MaMIMO channels and employs non-orthogonal pilot sequences. Bayesian learning methods have also been explored in the literature for JCD. In [17], a distributed semi-blind JCD message-passing algorithm has been presented for CF-MaMIMO networks which applies expectation propagation (EP) on a factor graph. A similar approach for grant-free CF-MaMIMO has been presented in [18]. The authors in [19] proposed a semi-blind JCD algorithm based on Bethe free energy optimization which combines variational Bayes (VB) and EP referred shortly to as VB-EP throughout this paper.

In this work, we propose a novel JCD algorithm based on EP, specifically designed to enhance robustness against PC in CF-MaMIMO systems. The algorithm builds upon the bilinear-EP algorithm presented in [17] and incorporates a modified scheduling and message passing for bilinear structure to effectively exploit the inherent structure of the received data signals and suppress the impact of PC. We evaluate the proposed method against state-of-the-art Bayesian learning algorithms and demonstrate superior detection and estimation performance, particularly under severe PC. This thorough analysis of PC was not conducted in prior work. Our anal-

<footnote>This work was funded by the DFG (German Research Foundation) for FAU and by the ANR for EURECOM – Project CellFree6G CO 1311/1-1, Project ID 491320625.</footnote>
ysis further considers both contamination caused by non-orthogonal sequences and by the reuse of identical orthogonal sequences, revealing that Bayesian methods exhibit greater robustness when non-orthogonal pilots with low correlation are used. This insight motivates the introduction of a novel metric tailored to quantify PC and assess its impact on iterative JCD algorithms.

**Notation:** $(\cdot)^T$ and $(\cdot)^H$ are the transpose and the conjugate transpose operator, respectively. $\delta(\cdot)$ and $\mathbb{1}(\cdot)$ denote the Dirac delta and the indicator function, respectively. $\mathcal{CN}(x|\mu, C)$ represents the circularly-symmetric multivariate complex Gaussian distribution of a complex-valued vector $x$ with mean $\mu$ and covariance matrix $C$. $\pi(x)$ denotes the categorical distribution of a discrete random variable $x$. The notation $x \sim p$ indicates that the random variable $x$ follows the distribution $p$. The message sent from the factor node $\Psi_\alpha$ to the variable node $x_\beta$ in a factor graph is denoted as $m_{\Psi_\alpha; x_\beta}$ and consists of parameters of a distribution in the exponential family which are denoted with the same subscript of the message, e.g., mean $\mu_{\Psi_\alpha; x_\beta}$ and covariance matrix $C_{\Psi_\alpha; x_\beta}$ for a Gaussian distribution or probability values $\pi_{\Psi_\alpha; x_\beta}$ for a categorical distribution. The same holds for variable-to-factor messages $m_{x_\beta; \Psi_\alpha}$.

## II. SYSTEM MODEL

We consider the uplink of a CF-MaMIMO network consisting of $L$ geographically distributed $N$-antenna APs and $K$ synchronized single-antenna UEs. All APs are connected to a central processing unit (CPU) via fronthaul links to share information. During the channel coherence time of $T$ channel uses, the received signal $\mathbf{Y} \in \mathbb{C}^{LN \times T}$ at all the APs is given by

$$ \mathbf{Y} = \mathbf{H}\mathbf{X} + \mathbf{N} \quad (1) $$

where $\mathbf{H} = [\mathbf{H}_1^T \cdots \mathbf{H}_L^T]^T \in \mathbb{C}^{LN \times K}$ is the channel matrix and $\mathbf{H}_l = [\mathbf{h}_{l,1} \cdots \mathbf{h}_{l,K}] \in \mathbb{C}^{N \times K}$ denotes the channel between AP $l$ and all UEs; $\mathbf{X} = [\mathbf{x}_1 \cdots \mathbf{x}_K]^T \in \mathbb{C}^{K \times T}$ is the transmit symbol matrix and $\mathbf{x}_k \in \mathbb{C}^{T \times 1}$ represents the transmit sequence of UE $k$; and $\mathbf{N} \in \mathbb{C}^{LN \times T}$ is the matrix of additive white Gaussian noise (AWGN) with independent and identically distributed (i.i.d.) elements $n \sim \mathcal{CN}(n|0, \sigma_n^2)$. The channels are assumed to be constant during the channel coherence time. We assume block Rayleigh fading channels between UE $k$ and AP $l$, i.e., $\mathbf{h}_{l,k} \sim p_{h_{l,k}}(\mathbf{h}_{l,k}) = \mathcal{CN}(\mathbf{h}_{l,k}|0, \Xi_{l,k})$ where $\Xi_{l,k}$ is the channel correlation matrix with the large-scale fading coefficient (LSFC) $\xi_{l,k} = \frac{1}{N}\text{tr}\{\Xi_{l,k}\}$. The transmit symbol matrix consists of a pilot part $\mathbf{X}^p \in \mathbb{C}^{K \times T_p}$ and a data part $\mathbf{X}^d \in \mathcal{X}^{K \times T_d}$, i.e., $\mathbf{X} = [\mathbf{X}^p \mathbf{X}^d]$, with $T_p + T_d = T$ and $\mathcal{X}$ being the transmit data symbol constellation of cardinality $M = |\mathcal{X}|$. The average transmit symbol power is given by $\sigma_x^2 = E\{|x_{kt}|^2\}$. A similar decomposition in pilot and data part applies to the receive matrix, i.e., $\mathbf{Y} = [\mathbf{Y}^p \mathbf{Y}^d]$ with received pilots $\mathbf{Y}^p \in \mathbb{C}^{L \times T_p}$ and received data $\mathbf{Y}^d \in \mathbb{C}^{L \times T_d}$. Furthermore, we assume that the pilot length is smaller than the number of UEs, $T_p < K$, since in practice the number of UEs can be very large and, thus, it is not practical to assign

orthogonal pilot sequences to the UEs. This gives rise to the PC effect.

## III. PROBLEM FORMULATION

Due to PC, channel estimation based solely on pilot sequences typically presents severely degraded performance. An effective strategy to mitigate the impact of PC is to leverage not only the pilot sequences but also the detected data symbols to iteratively refine both channel estimation and data detection. In this context, the receiver's task is to jointly estimate the channel matrix $\mathbf{H}$ and the user data matrix $\mathbf{X}^d$. The maximum a posteriori (MAP) estimator is given by

$$ (\hat{\mathbf{H}}, \hat{\mathbf{X}}^d) = \arg\max_{\mathbf{H}, \mathbf{X}^d} p_{\text{APP}}(\mathbf{H}, \mathbf{X}^d), \quad (2) $$

with the a posteriori probability (APP) distribution $p_{\text{APP}}(\mathbf{H}, \mathbf{X}^d)$ which can be factorized by applying Bayes' theorem,

$$ \begin{aligned} p_{\text{APP}}(\mathbf{H}, \mathbf{X}^d) = & p_{H,X^d|Y,X^p}(\mathbf{H}, \mathbf{X}^d|\mathbf{Y}, \mathbf{X}^p) \\ & \propto p_{Y|H,X}(\mathbf{Y}|\mathbf{H}, \mathbf{X}) \cdot p_H(\mathbf{H}) \cdot p_X(\mathbf{X}). \end{aligned} \quad (3) $$

Solving the inference problem in (2) is computationally intractable due to the high dimensionality of the involved variables. To address this challenge, we employ Bayesian learning techniques to obtain tractable approximations of the MAP estimates.

## IV. BILINEAR-EP ALGORITHM

In this section, we propose a novel JCD algorithm with enhanced performance under PC. It builds upon the bilinear-EP algorithm introduced in [17] for distributed semi-blind JCD in CF-MaMIMO systems. Unlike the baseline, our approach incorporates knowledge and observations of the pilot sequences in the message-passing procedure. In the following, we present the underlying factorization and message-passing steps, and refer to [17] for general details on bilinear-EP and message derivations.

### A. Factor Graph Representation

In order to solve the MAP problem in (2), the auxiliary variables $\mathbf{z}_{l,kt} := \mathbf{h}_{l,k}x_{kt} \forall l, k, t$ are introduced and stored in the array $\mathbf{Z}$. The APP distribution with respect to the channel, data, and auxiliary variables can be factorized as follows,

$$ p_{\text{APP}}(\mathbf{H}, \mathbf{X}^d, \mathbf{Z}) \propto \prod_{l=1}^{L} \prod_{k=1}^{K} \prod_{t=1}^{T} \left[ p(\mathbf{y}_{l,t} | \mathbf{z}_{l,1t}, \dots, \mathbf{z}_{l,Kt}) \right. \\ \left. \qquad \cdot p(\mathbf{z}_{l,kt} | \mathbf{h}_{l,k}, x_{kt}) \cdot \tilde{p}_{h_{l,k}}(\mathbf{h}_{l,k}) \cdot p_x(x_{kt}) \right], \quad (4) $$

where the independence of channel vectors for different APs and UEs as well as the independence of data symbols for different UEs and time indices is exploited. The probability distribution $\tilde{p}_{h_{l,k}}(\mathbf{h}_{l,k})$ represents modified channel information compared to the prior information $p_{h_{l,k}}(\mathbf{h}_{l,k})$. How to obtain this modified information will be explained in the following. The factor graph representing the APP distribution (4) is shown in Fig. 1. It comprises variable nodes (circles)
<figure><img src="image_2.png" /><figcaption>Fig. 1. Factor graph for bilinear-EP with $\mathcal{T} := T_p + 1$. The numbered red dashed arrows show the message update scheduling according to Algorithm 1.</figcaption></figure>

and factor nodes (rectangles), organized according to their implementation at the CPU and the APs. Each factor node corresponds to one of the following probability distributions,

$$ \Psi_{y_{l,t}} := p(\mathbf{y}_{l,t} | \mathbf{z}_{l,1t}, \dots, \mathbf{z}_{l,Kt}) = \mathcal{CN} \left( \mathbf{y}_{l,t} \mid \sum_{k=1}^{K} \mathbf{z}_{l,kt}, \sigma_n^2 \mathbf{I}_N \right), $$

$$ \Psi_{z_{l,kt}} := p(\mathbf{z}_{l,kt} | \mathbf{h}_{l,k}, x_{kt}) = \delta(\mathbf{z}_{l,kt} - \mathbf{h}_{l,k} x_{kt}), $$

$$ \Psi_{h_{l,k}} := \tilde{p}_{h_{l,k}}(\mathbf{h}_{l,k}) = \mathcal{CN}(\mathbf{h}_{l,k} | \tilde{\boldsymbol{\mu}}_{h_{l,k}}, \tilde{\mathbf{C}}_{h_{l,k}}), $$

$$ \Psi_{x_{kt}} := p_x(x_{kt}) = \begin{cases} 1 & \text{for } t \le T_p \\ \frac{1}{M} \mathbb{1}_{x_{kt} \in \mathcal{X}} & \text{for } t > T_p \end{cases}. $$

Here, $\tilde{p}_{h_{l,k}}(\mathbf{h}_{l,k})$ represents the modified prior channel distribution obtained solely from processing the pilot sequences. In this work, we employ the minimum mean squared error (MMSE) estimator which can be applied independently for each AP without loss of optimality. The resulting estimates are denoted by $\tilde{\boldsymbol{\mu}}_{h_{l,k}}$ with corresponding error covariance matrix $\tilde{\mathbf{C}}_{h_{l,k}}$ given by [20]

$$ \tilde{\boldsymbol{\mu}}_{h_{l,k}} = [\mathbf{F}\mathrm{vec}\{\mathbf{Y}_l^p\}]_{i_k, 1}, \qquad (5) $$

$$ \tilde{\mathbf{C}}_{h_{l,k}} = [(\mathbf{I}_{KN} - \mathbf{F}\tilde{\mathbf{X}}^p)\mathbf{\Xi}_l]_{i_k, i_k}, \qquad (6) $$

where $\mathbf{F} = \mathbf{\Xi}_l \tilde{\mathbf{X}}_l^p \mathbf{H} (\tilde{\mathbf{X}}_l^p \mathbf{\Xi}_l \tilde{\mathbf{X}}_l^p + \sigma_n^2 \mathbf{I}_{N T_p})^{-1}$, $\tilde{\mathbf{X}}_l^p = \mathbf{X}_l^p \mathbf{T} \otimes \mathbf{I}_N$, and $\mathbf{\Xi}_l = \mathrm{blkdiag}(\mathbf{\Xi}_{l,1}, \dots, \mathbf{\Xi}_{l,K})$. The operator $\mathrm{vec}\{\cdot\}$ denotes vectorization, $[\cdot]_{i,j}$ extracts elements located at rows and columns specified by the indices $i$ and $j$, respectively, of the matrix in the square brackets, $\mathbf{i}_k = (k-1)N+1: kN$ represents the index range corresponding to UE $k$, and the function $\mathrm{blkdiag}(\cdot)$ constructs a block diagonal matrix from its arguments.

The key modification of the proposed approach compared to the bilinear-EP algorithm in [17] lies in the inclusion of pilot symbols within the factor graph and, hence, within the message-passing procedure. The corresponding message update rules are detailed in the following section.

## B. Message Updates and Scheduling

The bilinear-EP algorithm models the variables $x_{kt}$ with categorical distributions, while $\mathbf{z}_{l,kt}$ and $\mathbf{h}_{l,k}$ are modeled as multivariate complex Gaussian distributions. Accordingly, EP message-passing update rules are applied to the factor graph in Fig. 1. We present the message scheduling and the final update rules. Detailed derivations can be found in [17]. Note that the mean vector $\boldsymbol{\mu}$ and the covariance matrix $\mathbf{C}$ of a Gaussian random variable are equivalently represented by the natural parameters $\gamma = \mathbf{C}^{-1}\boldsymbol{\mu}$ and $\mathbf{\Lambda} = \mathbf{C}^{-1}$. In the following, we will interchangeably use both these representations without explicitly stating the transformation, i.e., if $\boldsymbol{\mu}_{\Psi_{\alpha};\mathbf{x}_{\beta}}$ and $\mathbf{C}_{\Psi_{\alpha};\mathbf{x}_{\beta}}$ are computed, then $\gamma_{\Psi_{\alpha};\mathbf{x}_{\beta}}$ and $\mathbf{\Lambda}_{\Psi_{\alpha};\mathbf{x}_{\beta}}$ are automatically given and vice versa.

The message initialization is performed as follows: parameters describing the messages $m_{\Psi_{h_{l,k}};h_{l,k}}$, $m_{\mathbf{h}_{l,k};\Psi_{z_{l,kt}}}$, and $m_{\Psi_{z_{l,kt}};\mathbf{z}_{l,kt}}$ $\forall k, l, t$ are initialized per Algorithm 1. All other messages are initialized in an uninformative way.

The messages $m_{\Psi_{y_{l,t}};\mathbf{z}_{l,kt}}$ $\forall k, l, t$ are updated first, in which each AP performs interference cancellation on the received signal using the current knowledge of the auxiliary variables $\mathbf{z}_{l,kt}$,

$$ \boldsymbol{\mu}_{\Psi_{y_{l,t}}; \mathbf{z}_{l,kt}} = \mathbf{y}_{l,t} - \sum_{k' \neq k} \boldsymbol{\mu}_{\Psi_{z_{l,k't}}; \mathbf{z}_{l,k't}}, \qquad (7) $$

$$ \mathbf{C}_{\Psi_{y_{l,t}}; \mathbf{z}_{l,kt}} = \sigma_n^2 \mathbf{I}_N + \sum_{k' \neq k} \mathbf{C}_{\Psi_{z_{l,k't}}; \mathbf{z}_{l,k't}}. \qquad (8) $$

The updated information on the variables $\mathbf{z}_{l,kt}$ at each AP $l$ is used to refine the local beliefs on the data symbols $x_{kt}$, which are subsequently shared with the CPU. This is done by updating the message $m_{\Psi_{z_{l,kt}};x_{kt}}$ $\forall k, l, t > T_p$,

$$ \pi_{\Psi_{z_{l,kt}}; x_{kt}}(x_{kt}) \propto \theta(x_{kt}), \qquad (9) $$

with

$$ \theta(x_{kt}) = \mathcal{CN}(0 | \boldsymbol{\mu}_{\Psi_{y_{l,t}}; \mathbf{z}_{l,kt}} - \boldsymbol{\mu}_{\mathbf{h}_{l,k}; \Psi_{z_{l,kt}}}, x_{kt}, \mathbf{C}_{\Psi_{y_{l,t}}; \mathbf{z}_{l,kt}} + \mathbf{C}_{\mathbf{h}_{l,k}; \Psi_{z_{l,kt}}}, |x_{kt}|^2). $$

Next, the messages $m_{x_{kt};\Psi_{z_{l,kt}}}$ $\forall k, l, t > T_p$ are updated by aggregating the data symbol beliefs from all APs at the CPU, which then sends the following refined beliefs back to the APs,

$$ \pi_{x_{kt}; \Psi_{z_{l,kt}}}(x_{kt}) \propto \prod_{l' \neq l} \pi_{\Psi_{z_{l',kt}}; x_{kt}}(x_{kt}). \qquad (10) $$

The updated beliefs on the symbols $x_{kt}$ at each AP are then used to refine the channel estimates $\mathbf{h}_{l,k}$. This refinement
is achieved through the update of the message $m_{\Psi_{z_{l,kt}};\mathbf{h}_{l,k}}$ $\forall k, l, t$, given by

$$ \mathbf{\Lambda}_{\Psi_{z_{l,kt}};\mathbf{h}_{l,k}} = \hat{\mathbf{\Lambda}}_{1_{l,kt}} - \mathbf{\Lambda}_{\mathbf{h}_{l,k};\Psi_{z_{l,kt}}}, \quad (11) $$

$$ \gamma_{\Psi_{z_{l,kt}};\mathbf{h}_{l,k}} = \hat{\gamma}_{1_{l,kt}} - \gamma_{\mathbf{h}_{l,k};\Psi_{z_{l,kt}}}, \quad (12) $$

with $\hat{\mu}_{1_{l,kt}} = \frac{\tilde{\mu}_{1_{l,kt}}(x_{kt}^p)}{x_{kt}^p}$, $\hat{\mathbf{C}}_{1_{l,kt}} = \frac{\tilde{\mathbf{C}}_{1_{l,kt}}(x_{kt}^p)}{|x_{kt}^p|^2}$ for $t \le T_p$ and

$$ \hat{\mu}_{1_{l,kt}} = \frac{1}{Z_{l,kt}} \sum_{x_{kt} \in \mathcal{X}} \pi_{x_{kt};\Psi_{z_{l,kt}}}(x_{kt}) \cdot \frac{\theta(x_{kt})}{x_{kt}} \cdot \tilde{\mu}_{l,kt}(x_{kt}), $$

$$ \hat{\mathbf{C}}_{1_{l,kt}} = \frac{1}{Z_{l,kt}} \sum_{x_{kt} \in \mathcal{X}} \pi_{x_{kt};\Psi_{z_{l,kt}}}(x_{kt}) \cdot \frac{\theta(x_{kt})}{|x_{kt}|^2} \cdot (\tilde{\mathbf{C}}_{l,kt}(x_{kt}) + \tilde{\mu}_{l,kt}(x_{kt}) \cdot \tilde{\mu}_{l,kt}^{\mathrm{H}}(x_{kt})) - \hat{\mu}_{l,kt} \hat{\mu}_{l,kt}^{\mathrm{H}}, $$

for $t > T_p$ with

$$ Z_{l,kt} = \sum_{x_{kt} \in \mathcal{X}} \pi_{x_{kt};\Psi_{z_{l,kt}}}(x_{kt}) \cdot \theta(x_{kt}), $$

$$ \tilde{\mathbf{\Lambda}}_{l,kt}(x_{kt}) = \mathbf{\Lambda}_{\Psi_{y_{l,t}};\mathbf{z}_{l,kt}} + \mathbf{\Lambda}_{\mathbf{h}_{l,k};\Psi_{z_{l,kt}}}|x_{kt}|^{-2}, $$

$$ \tilde{\gamma}_{l,kt}(x_{kt}) = \gamma_{\Psi_{y_{l,t}};\mathbf{z}_{l,kt}} + \gamma_{\mathbf{h}_{l,k};\Psi_{z_{l,kt}}}|\frac{x_{kt}}{|x_{kt}|^2}|. $$

Then, the messages $m_{\mathbf{h}_{l,k};\Psi_{z_{l,kt}}}$ $\forall k, l, t$ are updated yielding new estimates of $\mathbf{h}_{l,k}$ by combining the information acquired across all time slots with the prior channel information,

$$ \mathbf{\Lambda}_{\mathbf{h}_{l,k};\Psi_{z_{l,kt}}} = \mathbf{\Lambda}_{\Psi_{h_{l,k}};\mathbf{h}_{l,k}} + \sum_{t' \neq t} \mathbf{\Lambda}_{\Psi_{z_{l,kt}'};\mathbf{h}_{l,k}}, \quad (13) $$

$$ \gamma_{\mathbf{h}_{l,k};\Psi_{z_{l,kt}}} = \gamma_{\Psi_{h_{l,k}};\mathbf{h}_{l,k}} + \sum_{t' \neq t} \gamma_{\Psi_{z_{l,kt}'};\mathbf{h}_{l,k}}. \quad (14) $$

The messages $m_{\Psi_{z_{l,kt}};\mathbf{z}_{l,kt}}$ $\forall k, l, t$ are updated last in an EP iteration, generating refined estimates of the variables $\mathbf{z}_{l,kt}$ which are then utilized for interference cancellation in the next iteration,

$$ \mathbf{\Lambda}_{\Psi_{z_{l,kt}};\mathbf{z}_{l,kt}} = \hat{\mathbf{\Lambda}}_{2_{l,kt}} - \mathbf{\Lambda}_{\Psi_{y_{l,t}};\mathbf{z}_{l,kt}}, \quad (15) $$

$$ \gamma_{\Psi_{z_{l,kt}};\mathbf{z}_{l,kt}} = \hat{\gamma}_{2_{l,kt}} - \gamma_{\Psi_{y_{l,t}};\mathbf{z}_{l,kt}}, \quad (16) $$

with $\hat{\mu}_{2_{l,kt}} = \tilde{\mu}_{l,kt}(x_{kt}^p)$, $\hat{\mathbf{C}}_{2_{l,kt}} = \tilde{\mathbf{C}}_{l,kt}(x_{kt}^p)$ for $t \le T_p$ and

$$ \hat{\mu}_{2_{l,kt}} = \frac{1}{Z_{l,kt}} \sum_{x_{kt} \in \mathcal{X}} \pi_{x_{kt};\Psi_{z_{l,kt}}}(x_{kt}) \cdot \theta(x_{kt}) \cdot \tilde{\mu}_{l,kt}(x_{kt}), $$

$$ \hat{\mathbf{C}}_{2_{l,kt}} = \frac{1}{Z_{l,kt}} \sum_{x_{kt} \in \mathcal{X}} \pi_{x_{kt};\Psi_{z_{l,kt}}}(x_{kt}) \cdot \theta(x_{kt}) \cdot (\tilde{\mathbf{C}}_{l,kt}(x_{kt}) + \tilde{\mu}_{l,kt}(x_{kt}) \cdot \tilde{\mu}_{l,kt}^{\mathrm{H}}(x_{kt})) - \hat{\mu}_{1_{l,kt}} \hat{\mu}_{1_{l,kt}}^{\mathrm{H}}, $$

for $t > T_p$ with $Z_{l,kt}$, $\tilde{\mathbf{\Lambda}}_{l,kt}(x_{kt})$, and $\tilde{\gamma}_{l,kt}(x_{kt})$ as given before.

The modified bilinear-EP algorithm is summarized in Algorithm 1. Compared to the EP algorithm presented in [17], it contains additional and augmented message updates. To be more precise, the updates of the messages $m_{\Psi_{y_{l,t}};\mathbf{z}_{l,kt}}$, $m_{\Psi_{z_{l,kt}};\mathbf{h}_{l,k}}$, $m_{\mathbf{h}_{l,k};\Psi_{z_{l,kt}}}$, and $m_{\Psi_{z_{l,kt}};\mathbf{z}_{l,kt}}$ are now also considered for $t \le T_p$ i.e., for the pilot part as well. Furthermore, the update of the message $m_{\mathbf{h}_{l,k};\Psi_{z_{l,kt}}}$ is enhanced by taking into account the additional information from $t \le T_p$.

## Algorithm 1 Modified Bilinear-EP Algorithm

**Input:** Pilot matrix $\mathbf{X}^p$, transmit power $\sigma_x^2$, received signal $\mathbf{Y}$, noise variance $\sigma_n^2$, prior distribution $\tilde{p}_{h_{l,k}}(\mathbf{h}_{l,k}) \equiv (\tilde{\mu}_{h_{l,k}}, \tilde{\mathbf{C}}_{h_{l,k}})$.

**Output:** Estimated channels $\hat{\mathbf{h}}_{l,k}$ and data $\hat{x}_{kt}^d$.

1: $\forall k, l, t$: Initialize all messages uninformatively except

$$ \mu_{\Psi_{h_{l,k}};\mathbf{h}_{l,k}} = \mu_{\mathbf{h}_{l,k};\Psi_{z_{l,kt}}} = \tilde{\mu}_{h_{l,k}}, $$

$$ C_{\Psi_{h_{l,k}};\mathbf{h}_{l,k}} = C_{\mathbf{h}_{l,k};\Psi_{z_{l,kt}}} = \tilde{C}_{h_{l,k}}, $$

$$ \mu_{\Psi_{z_{l,kt}};\mathbf{z}_{l,kt}} = \begin{cases} \tilde{\mu}_{h_{l,k}} x_{kt}^p & \text{for } t \le T_p \\ 0 & \text{for } t > T_p, \end{cases} $$

$$ C_{\Psi_{z_{l,kt}};\mathbf{z}_{l,kt}} = \begin{cases} \tilde{C}_{h_{l,k}} |x_{kt}^p|^2 & \text{for } t \le T_p \\ (\tilde{C}_{h_{l,k}} + \tilde{\mu}_{h_{l,k}} \tilde{\mu}_{h_{l,k}}^{\mathrm{H}}) \sigma_x^2 & \text{for } t > T_p. \end{cases} $$

2: **for** $i=1$ to $i_{\max}$ **do**

3: $\forall k, l, t$: Update $m_{\Psi_{y_{l,t}};\mathbf{z}_{l,kt}}$ via (7), (8).

4: $\forall k, l, t > T_p$: Update $m_{\Psi_{z_{l,kt}};x_{kt}}$ via (9).

5: $\forall k, l, t > T_p$: Update $m_{x_{kt};\Psi_{z_{l,kt}}}$ via (10).

6: $\forall k, l, t$: Update $m_{\Psi_{z_{l,kt}};\mathbf{h}_{l,k}}$ via (11), (12).

7: $\forall k, l, t$: Update $m_{\mathbf{h}_{l,k};\Psi_{z_{l,kt}}}$ via (13), (14).

8: $\forall k, l, t$: Update $m_{\Psi_{z_{l,kt}};\mathbf{z}_{l,kt}}$ via (15), (16).

9: **return** $\hat{\mathbf{h}}_{l,k}$ calculated via (17) $\forall k, l$.

10: **return** $\hat{x}_{kt}^d$ calculated via (18) $\forall k, t$.

After performing the final EP iteration, the channel and data estimates are computed as follows,

$$ \hat{\mathbf{h}}_{l,k} = \hat{\mathbf{\Lambda}}_{\mathbf{h}_{l,k}}^{-1} \hat{\gamma}_{h_{l,k}}, \quad (17) $$

$$ \hat{x}_{kt}^d = \arg \max_{x_{kt}^d \in \mathcal{X}} \hat{p}_{x_{kt}}(x_{kt}^d), \quad (18) $$

with

$$ \hat{\mathbf{\Lambda}}_{\mathbf{h}_{l,k}} = \mathbf{\Lambda}_{\Psi_{h_{l,k}};\mathbf{h}_{l,k}} + \sum_{t=1}^{T} \mathbf{\Lambda}_{\Psi_{z_{l,kt}};\mathbf{h}_{l,k}}, \quad (19) $$

$$ \hat{\gamma}_{\mathbf{h}_{l,k}} = \gamma_{\Psi_{h_{l,k}};\mathbf{h}_{l,k}} + \sum_{t=1}^{T} \gamma_{\Psi_{z_{l,kt}};\mathbf{h}_{l,k}}, \quad (20) $$

and the approximated posterior data distribution

$$ \hat{p}_{x_{kt}}(x_{kt}^d) \propto \prod_{l=1}^{L} \pi_{\Psi_{z_{l,kt}};x_{kt}}(x_{kt}^d). \quad (21) $$

We note that damping is applied to factor-to-variable messages using a damping parameter $\eta \in [0, 1]$ to improve the stability of the bilinear-EP algorithm [17], i.e., each updated parameter is computed as a convex combination of its previous and newly computed values. Furthermore, the parameters of the messages $m_{\Psi_{z_{l,kt}};\mathbf{h}_{l,k}}$ and $m_{\Psi_{z_{l,kt}};\mathbf{z}_{l,kt}}$ in line 6 and 8 of Algorithm 1 are updated only if the corresponding covariance/precision matrices obtained by (11) and (15), respectively, are symmetric positive definite. Otherwise, the parameters from the previous iteration are retained.

## V. QUANTIFYING PILOT CONTAMINATION

In this section, we introduce a metric to quantify PC in CF-MaMIMO networks. The level of PC is influenced by
the choice of the pilot matrix $X^p$ and the resulting correlations between pilot sequences. Mutual coherence, which measures the similarity between pilot sequences, was used in [21] for pilot design to mitigate PC. However, due to the distributed architecture of CF-MaMIMO, user signals can often be separated spatially, especially when the users are far apart, resulting in negligible interference. Therefore, the spatial power profiles of all UEs, captured by the LSFCs $\xi_{l,k}$, are also critical for characterizing PC. This motivates us to develop a new PC metric which is particularly suited for JCD in distributed systems. Inspired by the normalized mean squared error (NMSE) of the pilot-based MMSE channel estimator, we define the PC metric $c_k$ for UE $k$ as

$$c_k = \min_l \frac{\left[ \left( \operatorname{diag}(\xi_{l,1}, \dots, \xi_{l,K})^{-1} + \mathbf{X}^p \mathbf{X}^p \mathbf{H} \sigma_n^{-2} \right)^{-1} \right]_{k,k}}{\xi_{l,k}}, \quad (22)$$

where $\operatorname{diag}(\cdot)$ denotes a diagonal matrix with its inputs on the main diagonal. The rationale for taking the minimum value over all APs is that a strong, low-contamination link to any single AP is sufficient for reliable channel estimation and successful data detection. This reliable information can be used to cancel the interference caused by the corresponding UE and, hence, iteratively remove PC. This concept was formalized in [22] where sufficient and necessary conditions for semi-blind identifiability were established.

## VI. NUMERICAL RESULTS

In this section, we present Monte Carlo simulation results for the modified bilinear-EP algorithm and several state-of-the-art benchmark algorithms. We consider a network spanning an area of $400 \times 400 \text{ m}^2$ comprising $L = 16$ single-antenna APs, i.e., $N = 1$, placed on a regular grid at coordinates $\{(50 + i \cdot 100, 50 + j \cdot 100) \text{ m} | i, j \in \{0, 1, 2, 3\}\}$ and placed at a height of 10 m. A total of $K = 8$ UEs are placed uniformly at random ground locations and transmit $T_p = 4$ pilot symbols and $T_d \in \{10, 30\}$-quadrature amplitude modulation (QAM) uncoded data symbols. The values of $K$ and $T_p$ are chosen such that PC occurs and the complexity of the simulations is not too high. In practice, longer channel coherence times allow for an increase in the number of pilot symbols $T_p$ but also for an increase in the number of UEs $K$, especially when the number of AP antennas $LN$ is increased as well, such that PC still occurs in these more practical cases and needs to be mitigated. We consider two different choices of pilot sequences, referred to as Hadamard and discrete Fourier transform (DFT) pilots. For Hadamard pilots, $T_p$ orthogonal Hadamard pilot sequences are considered and shared among the $K$ UEs. For DFT pilots, the pilot matrix $X^p$ is obtained by truncating a $K \times K$ DFT matrix to the first $T_p$ columns, resulting in non-orthogonal sequences. Hence, in our simulations, the set of Hadamard pilots consists of four orthogonal pilot sequences, each of which is shared between two users. In contrast, the set of DFT pilots consists of eight unique but non-orthogonal pilot sequences. The receiver noise power at each AP is set to $\sigma_n^2 = -96 \text{ dBm}$. The LSFCs

<figure><img src="image_3.png" /><figcaption>Fig. 2. CDF of $c_k$ for different pilot sequences.</figcaption></figure>

are obtained using the 3GPP urban microcell model which incorporates correlated shadow fading [23].

The first set of results pertains to the PC metric introduced in (22) and evaluated over $10^5$ indenpedent large-scale fading realizations. The cumulative distribution function (CDF) of the PC metric $c_k$ is illustrated in Fig. 2. It can be observed that the non-orthogonal DFT pilots result in lower PC than the orthogonal Hadamard pilots which are reused among the UEs.

In the following, we show the channel estimation and data detection performance in terms of the NMSE of the channel estimates and the symbol error rate (SER). We compare the proposed modified bilinear-EP algorithm with the bilinear-EP algorithm in [17], the VB-EP algorithm presented in [19], and different MMSE estimators. For channel estimation, we consider the pilot-based MMSE estimator and the genie-aided MMSE estimator with perfect knowledge of the transmitted data symbols. For symbol detection, we employ the centralized MMSE MIMO detector with channel state information (CSI) obtained by the pilot-based MMSE estimator and with perfect CSI. Note that the MMSE detector with pilot-based CSI serves as an upper performance bound for any linear detection scheme with imperfect CSI, including those proposed in [14] and [23]. Both of the bilinear-EP algorithms execute 20 iterations for JCD, whereas the VB-EP algorithm runs for 40 iterations. All iterative algorithms use a damping parameter of $\eta = 0.5$. In the following figures, solid and dashed lines distinguish between different pilot sequences or data lengths, while colors and markers indicate the respective algorithms.

The next set of results present the NMSE and the SER as functions of the UE transmit power $\sigma_x^2$. The NMSE of the channel matrix estimate $\hat{\mathbf{H}}$ is defined as NMSE := $E\{\frac{\|\mathbf{H}-\hat{\mathbf{H}}\|^2}{\|\mathbf{H}\|^2}\}$, and the SER is obtained by averaging across all UEs, i.e., SER := $E\{\sum_k \sum_t \mathbb{1}_{x_{kt}^d \neq \hat{x}_{kt}^d} / (KT_d)\}$. The expectation operator in both definitions is computed with respect to the channel realizations. In our simulation results, the performance is obtained by averaging $10^4$ block transmissions where each block transmission corresponds to an independent realization of the UE positions. The SER results are shown in Fig. 3. It can be observed that the proposed bilinear-EP algorithm outperforms the linear MMSE
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. SER versus transmit power.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. NMSE versus transmit power.</figcaption></figure>

detector as well as the bilinear-EP algorithm in [17] and the
VB-EP algorithm in [19]. Systems adopting non-orthogonal
DFT pilots show significantly better performance compared
to those using orthogonal Hadamard pilots. Additionally, the
performance gain from increasing the number of data symbols
from T<sub>d</sub> = 10 to T<sub>d</sub> = 30 is more pronounced for DFT
pilots. The VB-EP algorithm, being designed for orthogonal
pilots, is not applicable to systems adopting the proposed
DFT pilots. The proposed bilinear-EP algorithm offers greater
flexibility in the design of pilot sequences, which represents an
additional advantage. Fig. 4 illustrates the NMSE performance
as a function of the transmit power for the same settings
discussed above. The proposed bilinear-EP algorithm, when

<figure><img src="image_7.png" /><figcaption>Fig. 5. NMSE versus iterations for T<sub>d</sub> = 10 and σ<sub>x</sub><sup>2</sup> = 16 dBm.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. CDF of SER<sub>k</sub> for σ<sub>x</sub><sup>2</sup> = 16 dBm.</figcaption></figure>

used with non-orthogonal DFT pilot sequences, closely ap-
proaches the performance of the genie-aided MMSE estimator,
a bound which is not attained with orthogonal Hadamard
pilot sequences. Fig. 5 shows the convergence behavior of
the iterative algorithms. The bilinear-EP algorithm converges
significantly faster than the VB-EP algorithm. As the bilinear-
EP algorithm is executed for only 20 iterations, the result for
iteration number 20 is extended to all following iterations to
enable direct comparison with the benchmark algorithms.

For the following set of results, the UE transmit power
is fixed to σ<sub>x</sub><sup>2</sup> = 16 dBm and the performance is as-
sessed in terms of the CDFs of the NMSE and the SER
per user, i.e., NMSE<sub>k</sub> := E{||h<sub>k</sub> - h&#773;<sub>k</sub>||<sup>2</sup> / ||h<sub>k</sub>||<sup>2</sup>} and SER<sub>k</sub> :=
E{∑<sub>t</sub> 1<sub>x<sub>kt</sub><sup>d</sup> ≠ x&#773;<sub>kt</sub><sup>d</sup></sub> / T<sub>d</sub>}. Here, 1000 independent realizations of
the UE positions are considered which results in 1000 · K =
8000 data points for the CDF. For each UE positioning real-
ization, the performance is averaged over 1000 independent
block transmissions, accounting for small-scale fading and
noise realizations. The corresponding results are illustrated in
Figs. 6 and 7. The CDFs show that the proposed modified
bilinear-EP algorithm outperforms the benchmark schemes.
Furthermore, the algorithms applied to systems utilizing non-
orthogonal DFT pilots show greater performance gains from
increased data lengths, especially when considering the 95%-.
likely performance.

For the final set of results, the exact same setup as before
<figure><img src="image_9.png" /><figcaption>Fig. 7. CDF of NMSE<sub>k</sub> for σ<sub>x</sub><sup>2</sup> = 16 dBm.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. SER<sub>k</sub> versus PC metric for T<sub>d</sub> = 30 and σ<sub>x</sub><sup>2</sup> = 16 dBm.</figcaption></figure>

is used with the UE-based performance metrics now evaluated as a function of the PC metric $c_k$. The results are averaged over UEs experiencing a similar level of PC quantified by $c_k$. Fig. 8 presents the corresponding results for the SER. It can be observed that as $c_k$ increases, the performance degrades which validates $c_k$ as an appropriate PC metric. Combined with the observations from Fig. 2, it explains the superior average performance of DFT pilots over Hadamard pilots in previous results. Furthermore, the proposed algorithm consistently achieves the best performance for a given level of PC. These results further show that the JCD schemes applied to systems using DFT pilots are more effective at mitigating PC for a given level of PC $c_k$.

## VII. CONCLUSION

In this work, we proposed a novel JCD algorithm based on EP, designed to improve robustness against PC in CF-MaMIMO systems. The algorithm extended the bilinear-EP method in [17] and achieved improved channel estimation and data detection performance, particularly under severe PC and for systems using non-orthogonal pilot sequences. It consistently outperformed optimal linear detectors in [14] and [23] and state-of-the-art JCD algorithms. We also compared systems employing orthogonal and non-orthogonal pilots and showed that non-orthogonal sequences provided significant performance gains. Finally, we introduced a new metric to

quantify the impact of PC on iterative JCD algorithms and
demonstrated its relevance and consistency.

## REFERENCES

[1] H. Q. Ngo, A. Ashikhmin, H. Yang, E. G. Larsson, and T. L. Marzetta, "Cell-free massive MIMO versus small cells," *IEEE Trans. Wireless Commun.*, vol. 16, no. 3, pp. 1834–1850, 2017.

[2] H. Q. Ngo, L.-N. Tran, T. Q. Duong, M. Matthaiou, and E. G. Larsson, "On the total energy efficiency of cell-free massive MIMO," *IEEE Trans. Green Commun. Netw.*, vol. 2, no. 1, pp. 25–39, 2018.

[3] M. Mohammadi, Z. Mobini, H. Quoc Ngo, and M. Matthaiou, "Next-generation multiple access with cell-free massive MIMO," *Proc. IEEE*, vol. 112, no. 9, pp. 1372–1420, 2024.

[4] H. Yin, D. Gesbert, and L. Cottatellucci, "Dealing with interference in distributed large-scale MIMO systems: A statistical approach," *IEEE J. Sel. Topics Signal Process.*, vol. 8, no. 5, pp. 942–953, 2014.

[5] Z. Chen and E. Björnson, "Channel hardening and favorable propagation in cell-free massive MIMO with stochastic geometry," *IEEE Trans. Commun.*, vol. 66, no. 11, pp. 5205–5219, 2018.

[6] R. Gholami, L. Cottatellucci, and D. Slock, "Favorable propagation and linear multiuser detection for distributed antenna systems," in *Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP)*, 2020.

[7] ——, "Channel models, favorable propagation and multistage linear detection in cell-free massive MIMO," in *Proc. IEEE Int. Symp. Inf. Theory (ISIT)*, 2020.

[8] H. Q. Ngo and E. G. Larsson, "EVD-based channel estimation in multicell multiuser MIMO systems with very large antenna arrays," in *Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP)*, 2012.

[9] H. Yin, D. Gesbert, M. Filippou, and Y. Liu, "A coordinated approach to channel estimation in large-scale multiple-antenna systems," *IEEE J. Sel. Areas Commun.*, vol. 31, no. 2, pp. 264–273, 2013.

[10] L. Cottatellucci, R. R. Müller, and M. Vehkapera, "Analysis of pilot decontamination based on power control," in *Proc. IEEE 77th Veh. Technol. Conf. (VTC-Spring)*, 2013.

[11] R. R. Müller, L. Cottatellucci, and M. Vehkapera, "Blind pilot decontamination," *IEEE J. Sel. Topics Signal Process.*, vol. 8, no. 5, pp. 773–786, 2014.

[12] H. Yin, L. Cottatellucci, D. Gesbert, R. R. Müller, and G. He, "Robust pilot decontamination based on joint angle and power domain discrimination," *IEEE Trans. Signal Process.*, vol. 64, no. 11, pp. 2990–3003, 2016.

[13] E. Björnson, J. Hoydis, and L. Sanguinetti, "Massive MIMO has unlimited capacity," *IEEE Trans. Wireless Commun.*, vol. 17, no. 1, pp. 574–590, 2018.

[14] A. Á. Polegre, L. Sanguinetti, and A. G. Armada, "Pilot decontamination processing in cell-free massive MIMO," *IEEE Commun. Lett.*, vol. 25, no. 12, pp. 3990–3994, 2021.

[15] S. Chen, J. Zhang, E. Björnson, J. Zhang, and B. Ai, "Structured massive access for scalable cell-free massive MIMO systems," *IEEE J. Sel. Areas Commun.*, vol. 39, no. 4, pp. 1086–1100, 2021.

[16] H. Song, T. Goldstein, X. You, C. Zhang, O. Tirkkonen, and C. Studer, "Joint channel estimation and data detection in cell-free massive MU-MIMO systems," *IEEE Trans. Wireless Commun.*, vol. 21, no. 6, pp. 4068–4084, 2022.

[17] A. Karataev, C. Forsch, and L. Cottatellucci, "Bilinear expectation propagation for distributed semi-blind joint channel estimation and data detection in cell-free massive MIMO," *IEEE Open J. Signal Process.*, vol. 5, pp. 284–293, 2024.

[18] C. Forsch, A. Karataev, and L. Cottatellucci, "Distributed joint user activity detection, channel estimation, and data detection via expectation propagation in cell-free massive MIMO," in *Proc. IEEE 25th Int. Workshop Signal Process. Advances Wireless Commun. (SPAWC)*, 2024, pp. 531–535.

[19] Z. Zhao and D. Slock, "Decentralized message-passing for semi-blind channel estimation in cell-free systems based on Bethe free energy optimization," in *Proc. 58th Asilomar Conf. Signals, Syst., and Comput.*, 2024, pp. 1443–1447.

[20] S. M. Kay, *Fundamentals of Statistical Signal Processing: Estimation Theory*. Prentice Hall PTR, 1993.

[21] H. Iimori, T. Takahashi, K. Ishibashi, G. T. F. de Abreu, and W. Yu, "Grant-free access via bilinear inference for cell-free MIMO with low-coherence pilots," *IEEE Trans. Wireless Commun.*, vol. 20, no. 11, pp. 7694–7710, 2021.

[22] R. Gholami, L. Cottatellucci, and D. Slock, "Tackling pilot contamination in cell-free massive MIMO by joint channel estimation and linear multi-user detection," in *Proc. IEEE Int. Symp. Inf. Theory (ISIT)*, 2021.

[23] E. Björnson and L. Sanguinetti, "Making cell-free massive MIMO competitive with MMSE processing and centralized implementation," *IEEE Trans. Wireless Commun.*, vol. 19, no. 1, pp. 77–90, 2020.