# Multipath Component Power Delay Profile Based Ranging

Fangqing Xiao, Zilu Zhao, Dirk Slock

► To cite this version:

Fangqing Xiao, Zilu Zhao, Dirk Slock. Multipath Component Power Delay Profile Based Ranging. IEEE Journal of Selected Topics in Signal Processing, 2024, 18 (5), pp.950-963. 10.1109/JSTSP.2024.3491580. hal-05045244

HAL Id: hal-05045244

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

Submitted on 24 Apr 2025

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 Attribution 4.0 International License
Multipath Component Power Delay Profile Based
Ranging

Fangqing Xiao, Student Member, IEEE, Zilu Zhao, Student Member, IEEE
and Dirk T. M. Slock, Fellow, IEEE

*Abstract*—Precision ranging technology has become indispensable for ensuring efficient, reliable, and low-latency fifth-generation (5G) networks. In this paper, we propose a novel ranging method which is multipath component (MPC) power delay profile (PDP) based ranging. Whereas the Received Signal Strength (RSS) only summarizes the PDP into a single characteristic, we aim to furthermore exploit the range dependent curvature of the PDP envelope over its delay spread. However, the multipath propagation only allows to sample the PDP envelope at the path delays and suffers from (slow) fading. Hence our approach involves constructing a statistical fading model of the PDP and establishing a relationship between the distribution parameters and the propagation distance. To theoretically validate the feasibility of our proposed method, we adopt the widely accepted Nakagami-m fading model, which renders traditional estimation methods difficult to apply. Therefore we introduce the Expectation Maximization (EM)-Revisited Vector Approximate Message Passing (ReVAMP) algorithm. This algorithm is specifically designed to handle difficulties in parameter estimation for Gaussian linear models (GLMs) with hidden random variables and intractable posterior distributions. Extensive numerical simulation results have been conducted which exhibit preliminary evidence of the effectiveness of our MPCPDP-based ranging method compared to the received signal strength (RSS)-based method. Moreover, the versatility of the EM-ReVAMP algorithm allows for its extension to other statistical fading models beyond the Nakagami-m model with minor modifications, which opens the door to potential improvements based on more accurate statistical fading models. Nevertheless, the applicability of our MPCPDP-based ranging method should be validated in real-world scenarios in future studies.

*Index Terms*—Ranging Estimation, Multipath Component, Power Delay Profile, Expectation Maximization, Revisited Vector Approximate Message Passing Algorithm

I. INTRODUCTION

WITH the advent of 5G communications, ranging technology has become indispensable, playing a pivotal role in ensuring efficient, reliable, and low-latency networks [1]. Its applications extend to both indoor and outdoor positioning services [2], supplying dependable location information for various scenarios, including smart cities [3], intelligent transportation [4], and the Internet of Things (IoT) [5].

Consequently, there is a growing interest in utilizing wireless signals for determining distances accurately. In complex environments, receivers may encounter multiple signals arriving through various paths, including both line-of-sight (LoS) and non-line-of-sight (NLoS) paths, intensifying the multipath effect as depicted in Fig. 1. Numerous studies have focused on mitigating or exploiting multipath components (MPCs) for

ranging and localization [6]–[19]. Some methods [6], [7] aim to mitigate the multipath interference. For example, Dardari et al. [6] describe a maximum-likelihood estimator (MLE) for Time of Arrival (ToA)-based ranging, along with practical low-complexity schemes to reduce the multipath interference. Wymeersch et al. [7] present a machine learning approach to directly mitigate bias in both LoS and NLoS conditions. Other studies leverage MPCs for ranging and localization [8]–[12]. Li et al. [12] propose a belief propagation (BP)-based algorithm for sequential channel estimation and detection (CEDA) of MPC parameters, such as the distance and the angle of arrival (AoA) based on radio signals. Building on CEDA, Venus et al. [10], [11] address obstructed LoS situations to provide high-accuracy position estimation. However, these methods often require additional information like the ToA, the direction of arrival (DoA), and the AoA, with channel fading amplitudes used solely for MPC detection. Additionally, some approaches jointly estimate and exploit scatterer positions, assuming a single-bounce multipath model and a simple topological map [13]–[19]. Despite constraints such as single-bounce models, these methods still necessitate at least three parameters (e.g., ToA, DoA, DoD) per path to provide useful information. Methods based on ToA typically assume very precise ToA estimates, leading to accurate range estimates as well. In summary, the aforementioned methods require additional information such as the DoA or the ToA, which necessitate extra hardware or synchronization. In contrast, received signal strength (RSS)-based ranging [20]–[22] offers

<figure><img src="image_2.png" /><figcaption>Fig. 1. Example of LoS and NLoS links.</figcaption></figure>
an alternative that may avoid these requirements.

As an integral of power delay profile (PDP), however, it was asserted that RSS, as a single number that exploits attenuation as a function of distance, is not an ideal metric in complex environments due to its vulnerability to multipath effects [23]–[25]. To solve this problem, in the state of the art, the PDP has been exploited in [26], [27] to handle multipath and NLoS via PDP fingerprinting, something that has become again popular more recently [28]. However, these methods still need additional information like ToA, etc..

In this paper, we propose a novel MPCPDP-based rang-
ing method requiring no additional hardware or estimation
information, by exploiting the evolution of the attenuation
over the whole delay spread and capturing the attenuation
at a number of delays (distances), which in particular allows
to capture the curvature of the PDP envelope (the distance
dependent attenuation function). Whereas the attenuation is
sensitive to calibration errors in synchronization and Tx/Rx
gain estimation, the curvature is insensitive to such offsets.
As demonstrated in our previous work [29], we confirm
that MPCPDP-based ranging achieves more precise range
estimation compared to RSS-based methods. However, the
PDP envelope is much harder to estimate from a single
channel realization than the single number which is RSS,
being sensitive to not only fast fading but also shadowing,
scatterer spread, and the fact that the PDP in some sense only
samples its envelope at the subset of multipath delays that
are present. To handle this last effect, we exploit the specular
part of the channel and focus on the attenuations of the MPC,
hence MPCPDP. The other challenge is the proper statistical
modeling of the shadowing, in combination with the distance-
dependent attenuation.

One widely accepted statistical model for indoor multi-
path propagation is the Saleh-Valenzuela model [30], which
takes into account factors such as reflection, diffraction and
scattering caused by indoor structures. However, this model
may not accurately represent channel behavior in outdoor or
wide-area environments, as it does not explicitly consider
wide-area path loss, shadowing, and other outdoor-specific
phenomena [31]. To accurately model outdoor or wide-area
channel behavior, different empirical models based on exten-
sive outdoor measurements, such as the Okumura-Hata model
[32] or some 3GPP models [33], are commonly used. Given
the multitude of options, selecting the appropriate statistical
model for validation is paramount. On one hand, we must
ensure its applicability across the majority of cases, while
on the other hand, we need to confirm that the algorithm
employed for this model can be extrapolated to other models.
In comparison to alternative models like Rayleigh, Rician,
or log-normal distributions, the Nakagami-m distribution [34]
demonstrated superior versatility and accuracy in fitting a wide
range of experimental data [35]. This superiority stems from
its capability to accommodate the superposition of primary
and clutter signals resulting from diffuse reflections within a
single path, making it a better fit than the Rayleigh distri-
bution [36]. It is worth noting that Nakagami-m and Rician
distribution models behave approximately equivalently near
their mean value. Therefore, we concluded that the Nakagami-

m decay model is well-suited for validating the feasibility of
our PDP-based ranging approach. However, the Nakagami-m
fading model diverges from the Rayleigh fading model in that
obtaining an analytic form for the likelihood function becomes
impractical due to intractable integrals, as will be elucidated
further in this paper.

Taking into account the uniform distribution of the phase
varying from 0 to 2π for each MPC, we considered how
to establish the relationship between the parameters of the
Nakagami-m distribution and the propagation distance to en-
able distance estimation. The shape parameter *m* of Nakagami-
m distribution is closely associated with the environment,
while the scale parameter *Ω* represents the average attenu-
ation power intensity [37], which is directly linked to the
propagation distance. Consequently, when the received data
contains sufficient information about the attenuation of dif-
ferent paths, we can directly estimate the distance based on
these measurements. While vast studies have explored the use
of MPCs for ranging/localization estimation, to the best of
our knowledge, no prior research has specifically focused on
directly estimating the propagation distance of the LoS path by
assuming that both the LoS and NLoS paths of PDP conform
to specific fading distributions.

With formulating the range estimation problem, we found
that traditional estimation approaches are not available. To
overcome this problem, we proposed the EM-ReVAMP algo-
rithm as an effective solution. The EM algorithm is employed
to handle estimation problem with hidden variables. In cases
where the analytic formula for the posterior probability density
function (pdf) is unavailable within the EM algorithm, we
introduce the ReVAMP inference algorithm to approximate the
posterior distribution [38]. Compared to the original VAMP
[39], which only provides averaged variances, reVAMP can
yield distinct variances. In reVAMP, each marginal extrinsic
distribution is approximated using a complex Gaussian distri-
bution through the process of approximation belief propaga-
tion. The simulations verify the theoretical feasibility of our
PDP-based ranging and validate the effectiveness of our EM-
ReVAMP algorithm.

However, we have not yet had a chance to test the proposed
MCPDP-based ranging on experimental data in various real
environments. This may call into question the applicability
of the Nakagami-m model in all scenarios. Fortunately, the
versatility of the EM-ReVAMP algorithm allows the MPCPDP
based ranging to extend to other statistical channel fading
models beyond the Nakagami-m model with minor modifica-
tions [38]. Even though our MCPDP-based ranging approach
may not be as precise as some state-of-the-art (SoTA) methods
that require more information and hence more procedures or
additional hardware, it outperforms the RSS-based ranging
method. Our simulation results validate this finding, making
our method a low-cost alternative when additional hardware
is not available. The paper makes the following key contribu-
tions:

* Proposal of a novel MPCPDP-based ranging method requiring only the channel impulse response, building statistical attenuation models for each MPC in the PDP.
* Introduction of the EM-ReVAMP algorithm as a practical and robust solution for EM based Maximum Likelihood estimation involving random hidden variables with intractable posteriors, as in the proposed MPCPDP based ranging approach.

* Verification of the superior accuracy and robustness of the proposed MPCPDP-based ranging method and the EM-ReVAMP algorithm with Nakagami-m statistical models through comprehensive simulations.

The organization of the remaining sections is as follows: Section II presents the system model, encompassing the orthogonal frequency-division multiplexing (OFDM) model, Nakagami-m amplitude fading of MPCs, and MPCPDP-based ranging estimation. Sections III and IV delve into detailed explanations of the EM algorithm and the ReVAMP algorithm, respectively. In Section V, we describe the procedure for estimating the LoS distance using the EM-ReVAMP algorithm and derive the Cramér-Rao bound (CRB) for our MPCPDP-based ranging method when the Nakagami-m distribution becomes a Rayleigh distribution with a shape parameter $m = 1$. Subsequently, Section VI showcases the simulation results. Finally, Section VII concludes with our findings and conclusions.

**Notation:** The following notation will be used throughout this paper. Column vectors are denoted by lowercase bold $\mathbf{x}$. Matrices are denoted by uppercase bold $\mathbf{X}$. Scalars are represented without bold, such as $\mathbf{x}$. The $i^{th}$ entry of a vector $\mathbf{x}$ is designated as $\mathbf{x}[i]$ or $x_i$. The element at row $i$ and column $j$ of matrix $\mathbf{X}$ is denoted as $X_{ij}$. The operation diag($\mathbf{X}$) is used to extract the column vector consisting of the main diagonal elements of the matrix $\mathbf{X}$. $p(\mathbf{x})$ denotes the pdf of continuous random vector. $(\cdot)^T$, $(\cdot)^H$ and $(\cdot)^*$ denote matrix transpose, Hermitian transpose and complex conjugation, respectively. $\|\cdot\|$ represents the Euclidean norm and $|\cdot|$ is the absolute value of a real number set or the module value of a complex number set. diag($\mathbf{X}$) denote the diagonal entries of matrix $\mathbf{X}$. $I_N$ is an identity matrix of dimension $N$ and $0_N$ denotes zero vector of dimension $N$. $C$ and $Z^+$ denote the domain of complex numbers and the set of positive integers, respectively.

## II. SYSTEM MODEL

In this section, we will offer a succinct introduction to the system model. We will begin with a brief overview of the OFDM signal model in Subsection II-A. Following that, Subsection II-B will delve into the Nakagami-m fading of MPCs' amplitudes. Lastly, in Subsection II-C, we will present our MPCPDP-based ranging estimation.

### A. OFDM model

The widely preferred modulation technique in communication networks is OFDM, which finds extensive application in 5G-NR [40]. In the OFDM model, the received baseband signal can be mathematically expressed as the convolution of the transmitted OFDM signal, denoted as $s(t)$, and the channel impulse response, denoted as $g(t)$. Additionally, complex

additive white Gaussian noise, represented as $v(t)$, is added to the received signal. This relationship can be represented as:

$$r(t) = s(t) * g(t) + v(t), \quad (1)$$

where * denotes the convolution operation. After the received signal, $r(t)$, is sampled at a rate of $T_s$, time and frequency synchronizations are performed prior to the $N$-point fast Fourier transform (FFT) operation. The output of the FFT, denoted as $\mathbf{y}$, can be written as:

$$\mathbf{y} = \mathbf{X}\mathbf{h} + \mathbf{v} \in \mathcal{C}^{N \times 1}, \quad (2)$$

where $\mathbf{X}$ is an $N \times N$ diagonal matrix containing the transmitted symbols on its diagonal, $\mathbf{h}$ represents the channel frequency response (CFR) as a vector, and $\mathbf{v}$ is a vector of independently and identically distributed (i.i.d.) complex zero-mean Gaussian noise samples with equal variance $\sigma_v^2$.

In the case of a block fading channel that remains constant over the duration of a packet, the channel impulse response (CIR) can be described as follows: [41]

$$g(t) = \sum_{l=0}^{L-1} a_l \delta(t - \kappa_l T_s), \quad (3)$$

where $a_l \in \mathcal{C}$ and $\kappa_l T_s (\kappa_0 < \kappa_1 < ... < \kappa_{L-1}$ and $\kappa_l \in Z^+)$ represent the gain and delay of the $l^{th}$ path, respectively, and $\delta(t)$ denotes the Kronecker delta function. Let

$$\mathbf{h} = [h_0, h_1, \dots, h_{N-1}]^T, \quad (4)$$

be the discrete CFR. Under the assumption that the sampling starts at $t = 0$, the $n$-th element of $\mathbf{h}$ can be written as [42]:

$$h_n = \sum_{l=0}^{L-1} a_l e^{-j\kappa_l\omega} \Big|_{\omega = \frac{2\pi[n]_N}{N}}, \quad (5)$$

where

$$[n]_N = \begin{cases} n, & n \le N/2 - 1, \\ n - N, & n \ge N/2 + 1. \end{cases} \quad (6)$$

Therefore, we can present (5) as

$$\mathbf{h} = \mathbf{T}\mathbf{a} \in \mathcal{C}^{N \times 1} \quad (7)$$

where $\mathbf{a} \in \mathcal{C}^{L \times 1}$ is a vector filled with fading gains and $\mathbf{T} \in \mathcal{C}^{N \times L}$ is a transformation matrix that $T_{kl} = e^{-j\kappa_l\omega} \Big|_{\omega = \frac{2\pi[k]_N}{N}}$.

### B. Nakagami-m amplitude fading of MPCs

As discussed in the previous subsection, the received signal in OFDM can be represented as follows:

$$\mathbf{y} = \mathbf{X}\mathbf{T}\mathbf{a} + \mathbf{v} = \mathbf{H}\mathbf{a} + \mathbf{v}; \quad \mathbf{v} \sim \mathcal{CN}(0, \sigma_v^2\mathbf{I}), \quad (8)$$

where $\mathbf{a} \in \mathcal{C}^{L \times 1}$ denotes the complex attenuation coefficients (amplitude $\mathbf{m}$ and phase $\phi$). For each individual element $a_i = m_i e^{j\phi_i}$ of $\mathbf{a}$, we assume its magnitude $m_i$ with a Nakagami-m distribution and phase $\phi_i$ with a uniform distribution. Therefore, the pdf of magnitude and phase can be expressed as follows:

$$p(m_i|\Omega_i) = \frac{2m^m m_i^{2m-1}}{\Gamma(m) \Omega_i^m} \exp\left[-\frac{m m_i^2}{\Omega_i}\right], m_i > 0, m_i \ge 0.5; \quad (9a)$$

$$p(\phi_i) = \frac{1}{2\pi}, \quad \phi_i \in [0, 2\pi), \quad (9b)$$
where $\Gamma(\cdot)$ denotes the gamma function $m$ is the shape parameter of the Nakagami-m distribution and $\Omega_i$ is the average power intensity of path $i$. The shape parameter $m$ controls the fading characteristics of the distribution. For lower values of $m$, the distribution resembles a Rayleigh distribution with a more rapid decay. As $m$ increases, the distribution becomes more concentrated around its mean, resembling a more concentrated fading behavior. In practice, $m$ is often estimated from channel measurements to accurately model the fading characteristics of the specific wireless channel. Referring to [43], the parameter $\Omega_i$ can be defined as:

$$ \Omega_i(d_0) = P_t G_t G_r \left[ \frac{\lambda}{4\pi(d_0 + c\tau_i)} \right]^n = G_0(d_0 + c\tau_i)^{-n}, \quad (10) $$

in the given equation, several variables are defined as follows: $P_t$ represents the transmitting power, $G_t$ denotes the transmitting antenna amplification, $\lambda$ is the wavelength of the electromagnetic wave, $c$ is the velocity of light, $n$ represents the propagation fading factor influenced by the environment, $d_0$ indicates the LoS distance, and $\tau_i$ indicates the propagation delay between the $i$-th path and the LoS path.

In (10), the term $P_t G_t G_r (\frac{\lambda}{4\pi})^n$ can be considered as a constant, denoted as $G_0$, which combines the effects of transmit power, antenna gains, wavelength, and path loss exponent. The propagation fading factor $n$ plays a crucial role in determining the rate of signal attenuation with distance and can vary depending on the characteristics of the wireless channel and the environment in which the signals propagate. As the propagation distance $d_0 + c\tau_i$ increases, $\Omega_i$ decreases following an inverse power-law relationship $(d_0 + c\tau_i)^{-n}$. This allows us to estimate the specific range $d_0$ based on $\Omega_i$ when $\tau_i$ is known in a given environment.

Using the Jacobi determinant [44], we can obtain the pdf of complex fading coefficient $a_i$ as follows:

$$ p_{a_i}(a_i|\Omega_i(d_0)) = \frac{m^m |a_i|^{2m-2}}{\pi \Gamma(m) \Omega_i^m} \exp \left[ -\frac{m |a_i|^2}{\Omega_i} \right]. \quad (11) $$

For simplicity, we denote $p_{a_i}(a_i|\Omega_i(d_0))$ by $p_{a_i}(a_i|d_0)$. Thus, the pdf of the collection $a$ can be given as:

$$ p_a(\mathbf{a}|d_0) = \prod_{i=0}^{L-1} p_{a_i}(a_i|d_0). \quad (12) $$

Before presenting the specific ranging estimation process, we assume the presence of a LoS path with an unknown distance $d_0$, as well as measurable time delays between NLoS paths and the LoS path. While acknowledging the possibility of measurement and calibration biases, this paper does not focus on their effects. Consequently, we disregard these biases in the subsequent estimation process.

### C. MPCDPP-based Ranging Estimation

Our objective is to estimate $d_0$ directly from $\mathbf{y}$. To achieve this, we will employ the maximum likelihood estimation (MLE) method, which transforms the problem into the following equation:

$$ \hat{d}_0 = \arg \max_{d_0} \ell(d_0; \mathbf{y}) = \arg \max_{d_0} \ln \mathcal{L}(d_0; \mathbf{y}), \quad (13) $$

where $\mathcal{L}(\cdot)$ and $\ell(\cdot)$ represent the likelihood function and log-likelihood function, respectively.

Regarding the optimization problem (13), the likelihood function can be expressed as:

$$ \mathcal{L}(d_0; \mathbf{y}) = p(\mathbf{y}|d_0) = \int p(\mathbf{a}, \mathbf{y}|d_0)d\mathbf{a} = \int p_{\mathbf{y}}(\mathbf{y}|\mathbf{a})p_{\mathbf{a}}(\mathbf{a}|d_0)d\mathbf{a}. \quad (14) $$

The pdf $p(\mathbf{y}|d_0)$ is crucial for estimating the LoS range $d_0$ based on the received signal $\mathbf{y}$ in (13). However, solving the integral problem directly to acquire $p(\mathbf{y}|d_0)$ proves to be intractable, as finding an analytical form poses significant challenges. Furthermore, the latent variable $\mathbf{a}$ is unobserved, and its distribution is unknown before reaching $d_0$. To tackle these challenges, the EM-ReVAMP algorithm is introduced in Sections III and IV.

### III. REVIEW OF EXPECTATION MAXIMIZATION (EM)

As we discussed before, in the linear mixing data model described by (8), we have a known measurement matrix $\mathbf{H} \in \mathcal{C}^{M \times L}$ and an non-identically and independent distributed (n.i.i.d.) prior $p_{\mathbf{a}}(\mathbf{a}|d_0) = \prod_{i=0}^{L-1} p_{a_i}(a_i|d_0)$ for the vector $\mathbf{a}$. Additionally, we consider a zero-mean Gaussian measurement noise $p(\mathbf{v}) = \mathcal{CN}(\mathbf{v}; \mathbf{0}_M, \mathbf{C}_{vv})$ with covariance matrix $\mathbf{C}_{vv} \in \mathcal{R}^{M \times M}$.

To address the optimization problem (13), the Expectation-Maximization (EM) algorithm [45] proves to be a suitable solution. This algorithm is effective for estimation problems involving latent variables, such as $\mathbf{a}$, which are unobserved.

Using minorization maximization (MM) [46], we construct a more easily optimized lower bound of the log-likelihood function and iteratively approximate the optimal parameters by continuously optimizing this lower bound. Assuming at $t$-th iteration that we have the estimated $d_0^{(t)}$, which allows us to write

$$ \begin{aligned} \ell(d_0) - \ell(d_0^{(t)}) &= \ln \int p_{\mathbf{y}}(\mathbf{y}|\mathbf{a})p_{\mathbf{a}}(\mathbf{a}|d_0)d\mathbf{a} - \ln p(\mathbf{y}|d_0^{(t)}) \\ &= \ln \int \frac{p_{\mathbf{y}}(\mathbf{y}|\mathbf{a})p_{\mathbf{a}}(\mathbf{a}|d_0)}{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} p(\mathbf{a}|\mathbf{y}, d_0^{(t)})d\mathbf{a} \\ &\quad - \int p(\mathbf{a}|\mathbf{y}, d_0^{(t)}) \ln p(\mathbf{y}|d_0^{(t)})d\mathbf{a}, \end{aligned} \quad (15) $$

where $p(\mathbf{a}|\mathbf{y}, d_0)$ is the posterior distribution of $\mathbf{a}$ as which can be expressed by Bayes' rule as

$$ p(\mathbf{a}|\mathbf{y}, d_0) = \frac{p_{\mathbf{y}}(\mathbf{y}|\mathbf{a})p_{\mathbf{a}}(\mathbf{a}|d_0)}{p(\mathbf{y}|d_0)} = \frac{p_{\mathbf{y}}(\mathbf{y}|\mathbf{a})p_{\mathbf{a}}(\mathbf{a}|d_0)}{\int p_{\mathbf{y}}(\mathbf{y}|\mathbf{a})p_{\mathbf{a}}(\mathbf{a}|d_0)d\mathbf{a}}. \quad (16) $$

By using the concavity of $\ln(\cdot)$ and Jensen's inequality, (15)
becomes

$$
\begin{align}
\ell(d_0) - \ell(d_0^{(t)}) &\ge \int p(\mathbf{a}|\mathbf{y}, d_0^{(t)}) \ln \frac{p_y(\mathbf{y}|\mathbf{a})p_\mathbf{a}(\mathbf{a}|d_0)}{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} d\mathbf{a} \nonumber \\
&\quad - \int p(\mathbf{a}|\mathbf{y}, d_0^{(t)}) \ln p(\mathbf{y}|d_0^{(t)}) d\mathbf{a} \nonumber \\
&= \int p(\mathbf{a}|\mathbf{y}, d_0^{(t)}) \ln \frac{p_y(\mathbf{y}|\mathbf{a})p_\mathbf{a}(\mathbf{a}|d_0)}{p(\mathbf{y}|d_0^{(t)})p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} d\mathbf{a} \nonumber \\
&= \mathcal{E}_{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} \left[ \ln \frac{p_y(\mathbf{y}|\mathbf{a})p_\mathbf{a}(\mathbf{a}|d_0)}{p(\mathbf{y}|d_0^{(t)})p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} \right]. \tag{17}
\end{align}
$$

The lower bound $\mathcal{B}(d_0, d_0^{(t)})$ can be obtained as

$$
\mathcal{B}(d_0, d_0^{(t)}) = \ell(d_0^{(t)}) + \mathcal{E}_{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} \left[ \ln \frac{p_y(\mathbf{y}|\mathbf{a})p_\mathbf{a}(\mathbf{a}|d_0)}{p(\mathbf{y}|d_0^{(t)})p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} \right]. \quad (18)
$$

The updated $d_0^{(t+1)}$ can be obtained from

$$
\begin{align}
d_0^{(t+1)} &= \arg \max_{d_0} \mathcal{B}(d_0, d_0^{(t)}) \nonumber \\
&= \arg \max_{d_0} \left\{ \ell(d_0^{(t)}) + \mathcal{E}_{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} \left[ \ln \frac{p_y(\mathbf{y}|\mathbf{a})p_\mathbf{a}(\mathbf{a}|d_0)}{p(\mathbf{y}|d_0^{(t)})p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} \right] \right\} \nonumber \\
&= \arg \max_{d_0} \mathcal{E}_{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} [\ln p(\mathbf{y}, \mathbf{a}|d_0)]. \tag{19}
\end{align}
$$

At convergence we get $d_0^{(t)} = d_0^{(t+1)}$ which leads to the inequality becoming an equality:

$$
\begin{align}
\ln \int \frac{p_y(\mathbf{y}|\mathbf{a})p_\mathbf{a}(\mathbf{a}|d_0^{(t+1)})}{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} p(\mathbf{a}|\mathbf{y}, d_0^{(t)}) d\mathbf{a} &= \ln p(\mathbf{y}|d_0^{(t+1)}) \nonumber \\
&= \int p(\mathbf{a}|\mathbf{y}, d_0^{(t)}) \ln \frac{p_y(\mathbf{y}|\mathbf{a})p_\mathbf{a}(\mathbf{a}|d_0^{(t+1)})}{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} d\mathbf{a} \nonumber \\
\iff \ell(d_0^{(t+1)}) - \ell(d_0^{(t)}) &= 0. \tag{20}
\end{align}
$$

This proves that the EM algorithm can converge to the (local)
optimal point. Moreover, the EM iteration can be specified as:

$$
\begin{align}
d_0^{(t+1)} &= \arg \max_{d_0} \mathcal{E}_{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} [\ln p(\mathbf{a}, \mathbf{y}|d_0)] \nonumber \\
&= \arg \max_{d_0} \mathcal{E}_{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} [\ln p_{\mathbf{a}}(\mathbf{a}|d_0) + \ln p_y(\mathbf{y}|\mathbf{a})] \nonumber \\
&= \arg \max_{d_0} \mathcal{E}_{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} [\ln p_{\mathbf{a}}(\mathbf{a}|d_0)]. \tag{21}
\end{align}
$$

When considering the pdf of **a** as described in (11) and (12),
the EM iteration in (21) can be transformed as follows:

$$
\begin{align}
d_0^{(t+1)} &= \arg \max_{d_0} \mathcal{E}_{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} \left[ \sum_{i=0}^{L-1} \left( -\ln \Omega_i(d_0) - \frac{|a_i|^2}{\Omega_i(d_0)} \right) \right] \nonumber \\
&= \arg \min_{d_0} \sum_{i=0}^{L-1} \left[ \ln \Omega_i(d_0) + \frac{\mathcal{E}_{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})} [|a_i|^2]}{\Omega_i(d_0)} \right], \tag{22}
\end{align}
$$

where $\Omega_i(d_0)$ was defined in (10). Here, as it is easy to find
that both $\ln \Omega_i(d_0)$ and $\frac{1}{\Omega_i(d_0)}$ are convex functions w.r.t. $d_0$.
This characteristic ensures that the entire optimization function
is convex, with only one global minimum point for $d_0$.

However, in this scenario, the EM algorithm remains
intractable because obtaining the posterior distribution

$p(\mathbf{a}|\mathbf{y}, d_0^{(t)})$ is challenging due to the integration involved in (16). Therefore, it becomes crucial to develop an algorithm that approximates this posterior distribution with another tractable distribution. To achieve this goal, we propose an algorithm called Revisited Vector Approximate Message Passing (Re-VAMP).

IV. REVISITED VECTOR APPROXIMATE MESSAGE PASSING (REVAMP)

Our objective is to find a distribution $q(\mathbf{a})$ that approximates
the posterior distribution $p(\mathbf{a}|\mathbf{y}, d_0^{(t)})$ with minimal Kullback-
Leibler divergence (KLD) between $p(\mathbf{a}|\mathbf{y}, d_0^{(t)})$ and $q(\mathbf{a})$,
where $d_0^{(t)}$ is given. Mathematically, we have:

$$
\hat{q}(\mathbf{a}) = \arg \min_{q(\mathbf{a})} D_{KL} [p(\mathbf{a}|\mathbf{y}, d_0) \| q(\mathbf{a})]. \quad (23)
$$

As seen in (22), the posterior distribution $p(\mathbf{a}|\mathbf{y}, d_0^{(t)})$ is used to get second order moment of each $a_i$. To accomplish this, we choose $q(\mathbf{a})$ as complex Gaussian distribution $\mathcal{CN}(\mathbf{a}; \mathbf{m}, \mathcal{C}_m)$, which turns to match the first-order and second-order moments of $q(\mathbf{a})$ and $p(\mathbf{a}|\mathbf{y}, d_0^{(t)})$ [47] if (23) is satisfied as follows:

$$
\mathcal{E}_{q(\mathbf{a})}[\mathbf{a}] = \mathcal{E}_{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})}[\mathbf{a}], \quad (24a)
$$

$$
\mathcal{E}_{q(\mathbf{a})}[\mathbf{a}\mathbf{a}^H] = \mathcal{E}_{p(\mathbf{a}|\mathbf{y}, d_0^{(t)})}[\mathbf{a}\mathbf{a}^H]. \quad (24b)
$$

It ensures that even with an approximate distribution $q(\mathbf{a})$, the
updated estimation value of $d_0^{(t+1)}$ in (22) will be the same
as if it were calculated using the true posterior distribution
$p(\mathbf{a}|\mathbf{y}, d_0^{(t)})$. However, solving equation (23) directly is com-
putationally intractable. Therefore, we introduce the ReVAMP
algorithm [38] as a solution.

To begin, we consider the factorization of the joint distri-
bution as follows:

$$
p(\mathbf{a}, \mathbf{y}|d_0) = p_y(\mathbf{y}|\mathbf{a}) \prod_{i=0}^{L-1} p_{a_i}(a_i|d_0). \quad (25)
$$

This factorization can be represented as a factor graph, as
shown in Fig. 2. In this graph, the variable nodes $a_i$, where
$i = 0, \dots, L-1$, are connected to the factor nodes $p_y(\mathbf{y}|\mathbf{a})$
and $p_{a_i}(a_i)$.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Factor graph of ReVAMP</figcaption></figure>
Applying the sum-product rule [39], the message $\mu_{py \to a_i}(a_i)$ passed from left factor node $p_y(y|a)$ to variable node $a_i$ can be expressed as:

$$ \mu_{py \to a_i}(a_i) \propto \int p_y(y|a) \prod_{j \neq i} \mu_{a_j \to p_y}(a_j) da_j, \quad (26) $$

where $\mu_{a_j \to p_y}(a_j)$ represents the message passed from variable node $a_j$ to left factor node $p_y(y|a)$. In addition, the message $\mu_{p_{a_i} \to a_i}(a_i)$ passed from right factor node $p_{a_i}$ to variable node $a_i$ can be represented as:

$$ \mu_{p_{a_i} \to a_i}(a_i) = p_{a_i}(a_i). \quad (27) $$

At variable node $a_i$, the received messages $b_{sp}(a_i) \propto \mu_{py \to a_i}(a_i)\mu_{p_{a_i} \to a_i}(a_i)$ are approximated with a complex Gaussian belief (approximated posterior) $b_{app}(a_i) = q(a_i) = CN(a_i; \hat{a}_i, \tau_{a_i})$ by minimizing the KLD. This approximation is formulated as:

$$ \hat{b}_{app} = \arg \min_{b_{app}} D_{KL}(b_{sp}(a_i) || b_{app}(a_i)). \quad (28) $$

Analogous to the sum-product rule, the messages $\mu_{a_i \to p_y}(a_i)$ shown in (26) can be calculated by:

$$ \mu_{a_i \to p_y}(a_i) = \frac{b_{app}(a_i)}{\mu_{py \to a_i}(a_i)}. \quad (29) $$

In the following, we will provide a detailed derivation of Re-VAMP for the case where the measurement noise is assumed to be complex Gaussian.

### A. Extrinsic to variable nodes

Assuming that at each iteration, the message $\mu_{a_i \to p_y}(a_i)$ passed from each variable node $a_i$ to left factor node $p_y(y|a)$ is redefined as assumed prior $q_i(a_i)$ for all $i = 0, \dots, L-1$, where $q_i(a_i)$ is supposed to be a complex Gaussian distribution. In this case, we can prove by induction that they will remain complex Gaussian. Without loss of generality, let us define $q_i(a_i) = CN(a_i; p_i, \tau_{p_i})$, where $p_i$ and $\tau_{p_i}$ are the mean and variance of assumed prior of each $a_i$, respectively.

In this case, the joint distribution $\prod_{i=0}^{L-1} q_i(a_i)$ is equal to $CN(a; p, D_p)$, where $p = [p_0 \dots p_{L-1}]^T$ and $D_p$ is a diagonal matrix whose $i$th entry is $\tau_{p_i}$. The real posterior $p(a|y)$ is approximated as $q(a) = CN(a; m, C_m) \propto p(y|a) \prod_{i=0}^{L-1} q_i(a_i)$.

In addition, the message $\mu_{py \to a_i}(a_i)$ is recalled as the extrinsic message $q_y(a_i)$ for all $i = 0, \dots, L-1$, where $q_y(a_i)$ is supposed to be a complex Gaussian distribution with mean $r_i$ and variance $\tau_{r_i}$. The extrinsic message for any variable node $a_i$ is obtained by:

$$ q_y(a_i) \propto \frac{\int_{a_{i}} p(y|a) CN(a; p, D_p) da_{i}}{q_i(a_i)} \quad (30) \\ \propto \frac{\int_{a_{i}} CN(a; m, C_m) da_{i}}{CN(a_i; p_i, \tau_{p_i})}, $$

where $a_{i}$ represents a vector that is the same as $a$ except that it excludes the $i$-th entry, with

$$ C_m = (H^H C_{vv}^{-1} H + D_p^{-1})^{-1}, \quad (31a) $$

$$ m = C_m (H^H C_{vv}^{-1} y + D_p^{-1} p). \quad (31b) $$

By following these steps, we can derive the extrinsic messages for the variable nodes, which play a crucial role in the algorithm for handling the complex Gaussian measurement noise. Additionally, we define $\tau_m = [\tau_{m_0} \dots \tau_{m_{L-1}}]^T = \text{diag}(C_m)$. Exploiting the properties of multivariate complex Gaussian distribution and (30), the extrinsic message $q_y(a_i)$ is represented by the complex Gaussian distribution $CN(a_i; r_i, \tau_{r_i})$ with

$$ r_i = \frac{\tau_{p_i} m_i - \tau_{m_i} p_i}{\tau_{p_i} - \tau_{m_i}}, \quad (32a) $$

$$ \tau_{r_i} = \frac{\tau_{m_i} \tau_{p_i}}{\tau_{p_i} - \tau_{m_i}}. \quad (32b) $$

To approximate the belief $q(a_i) = CN(a_i; \hat{a}_i, \tau_{a_i})$ at variable node $a_i$ as a complex Gaussian distribution, we minimize the KLD as

$$ \begin{gathered} \arg \min_{q(a_i)} D_{KL} [p(a_i)CN(a_i; r_i, \tau_{r_i}) || q(a_i)] \\ \Leftrightarrow \arg \min_{\hat{a}_i, \tau_{a_i}} D_{KL} [p(a_i)CN(a_i; r_i, \tau_{r_i}) || CN(a_i; \hat{a}_i, \tau_{a_i})]. \end{gathered} \quad (33) $$

Define the normalization factor as

$$ Z_i(r_i, \tau_i) = \int p(a_i)CN(a_i; r_i, \tau_{r_i})da_i, \quad (34) $$

then we obtain

$$ \hat{a}_i = \frac{\int a_i p_{a_i}(a_i) CN(a_i; r_i, \tau_{r_i}) da_i}{Z_i(r_i, \tau_i)}, \quad (35a) $$

$$ \tau_{a_i} = \frac{\int |a_i - \hat{a}_i|^2 p_{a_i}(a_i) CN(a_i; r_i, \tau_{r_i}) da_i}{Z_i(r_i, \tau_i)}. \quad (35b) $$

It's worth noting that $p_{a_i}(a_i)$ in (35) is not restricted to the Nakagami-m propagation model prior of (11); it can accommodate other priors as well. When transitioning to a different statistical propagation model, minor adjustments of prior should be made to ensure the continued functionality of the reVAMP algorithm.

### B. Passing the Approximation to the Factor Node

The assumed prior $q_i(a_i)$ can be expressed as the quotient of two complex Gaussian pdfs w.r.t. $q(a_i)$ and $q_y(a_i)$. This ensures that the resulting message distribution $q_i(a_i)$ remains complex Gaussian. Specifically, it is defined as:

$$ q_i(a_i) = CN(a_i; p_i, \tau_{p_i}) \propto \frac{CN(a_i; \hat{a}_i, \tau_{a_i})}{CN(a_i; r_i, \tau_{r_i})}. \quad (36) $$

From (36), we can determine $p_i$ and $\tau_{p_i}$ as follows:

$$ p_i = \frac{\tau_{r_i} \hat{a}_i - \tau_{a_i} r_i}{\tau_{a_i} \tau_{r_i}}, \quad (37a) $$

$$ \tau_{p_i} = \frac{\tau_{r_i} \tau_{a_i}}{\tau_{r_i} - \tau_{a_i}}. \quad (37b) $$

It is worth noting that when using the sequential updating method, the complexity of the matrix inverse operation in (31) can be reduced by employing matrix inverse lemma. Let us denote the result of $\tau_{p_i}^{new}$ during the update messages of the $a_i$ and define $\Delta_{p_i} = \frac{\tau_{p_i} - \tau_{p_i}^{new}}{\tau_{p_i} \tau_{p_i}^{new}}$. Moreover, we define $h_C(\cdot)$
as the updating of the $C_m$ with the new value of $\tau_{p_i}^{new}$ as follows:

$$ 
\begin{aligned}
C_m^{new} &= h_C(C_m, e_i, \Delta_{p_i}) = [C_m^{-1} + \Delta_{p_i} e_i e_i^T]^{-1} \\
&= C_m - C_m e_i (1/\Delta_{p_i} + e_i^T C_m e_i)^{-1} e_i^T C_m, \quad (38)
\end{aligned}
$$

where $e_i$ is a unit vector with only the i-th entry set to 1. To handle the cycles, we define $e_0 = e_N$.

The computation for updating $m$ in (31) can also be simplified with the same technique. We define $\Psi_{p_i} = \frac{p_i^{new}}{\tau_{p_i}^{new}} - \frac{p_i}{\tau_{p_i}}$ and denote $h_m(\cdot)$ as its update equation as follows:

$$ 
\begin{aligned}
m^{new} &= h_m(m, C_m, e_i, \Delta_{p_i}, \Psi_{p_i}) \\
&= C_m^{new} (H^T C_{vv}^{-1} y + D_p p + \Psi_{p_i} e_i) \\
&= m + \frac{\Psi_{p_i} - \Delta_{p_i} e_i^T m}{1 + \Delta_{p_i} e_i^T C_m e_i} C_m e_i. \quad (39)
\end{aligned}
$$

In summary, our algorithm iteratively computes messages from factor nodes to variable nodes and subsequently calculates messages from variable nodes back to factor nodes until convergence is achieved. The final approximation for $p(a|y)$ is represented by $q(a) = CN(a; m, C_m)$. Importantly, these update steps can be performed in parallel, yielding a similar algorithm to VAMP but with individual variance updates. By leveraging the matrix inverse lemma, the sequential update method maintains the same complexity as the parallel update method. Algorithm 1 delineates the detailed steps of this process.

### C. Relation to Expectation Propagation (EP)

Algorithm 1 can be regarded as an EP algorithm by approximating the factorization in Equation (25) as follows:

$$ p(a|y) \approx q(a) \propto p(a, y) \approx p_y(y|a) \prod_{i=0}^{L-1} q_i(a_i), \quad (40) $$

where each $q_i(a_i)$ is a Gaussian distribution with mean $p_i$ and variance $\tau_{p_i}$. To further explore the EP connection, let us consider the optimization problem:

$$ 
\begin{aligned}
\arg \min_{q^{new}} D_{KL} \left[ \frac{q(a)}{q_i(a_i)} p(a_i) \middle| q^{new}(a) \right] \\
= \arg \min_{q^{new}} \int_{a_i} \int_{a_{/i}} q(a) da_{/i} \frac{p(a_i)}{q_i(a_i)} \ln \left[ \frac{p(a_i)}{q^{new}(a_i)} \right] da_i. \quad (41)
\end{aligned}
$$

Let us continue by introducing $h(a_i) = \int_{a_{/i}} [q(a) da_{/i}] / q_i(a_i)$ and then we can rewrite the optimization problem (41) as follows:

$$ 
\begin{aligned}
\arg \min_{q^{new}} \int_{a_i} h(a_i) p_{a_i}(a_i) \ln \left[ \frac{h(a_i) p_{a_i}(a_i)}{h(a_i) q^{new}(a_i)} \right] da_i \\
= \arg \min_{q^{new}} D_{KL} [h(a_i) p_{a_i}(a_i) \middle| h(a_i) q^{new}(a_i)]. \quad (42)
\end{aligned}
$$

In Algorithm 1, this marginal extrinsic is represented as a complex Gaussian distribution with mean $r_i$ and variance $\tau_{r_i}$. Lastly, the first equality in (41) holds because during the update for the i-th entry, for all $k \neq i$, the minimum is achieved when $q^{new}(a_k) = q(a_k)$.

### Algorithm 1 ReVAMP (Complex Gaussian measurement noise via sequential updating)

**Ensure:** $m, C_m$

**Require:** $y, H, p_a(a), p(v)$

1. Initialize: $m, C_m, p, \tau_p$

2. repeat

3. repeat [For each $i = 0...L-1$]

4. [Update the extrinsic]

5. $\tau_m = \text{diag}(C_m)$

6. $r_i = \frac{\tau_{p_i} m_i - \tau_{m_i} p_i}{\tau_{p_i} - \tau_{m_i}}$

7. $\tau_{r_i} = \frac{\tau_{m_i} \tau_{p_i}}{\tau_{p_i} - \tau_{m_i}}$

8. [Approximate the marginal posterior]

9. Update $\hat{a}_i$ with (35a)

10. Update $\tau_{a_i}$ with (35b)

11. [Propagate the approximation back]

12. $p_i^{new} = \frac{\tau_{r_i} \hat{a}_i - \tau_{a_i} r_i}{\tau_{r_i} - \tau_{a_i}}$

13. $\tau_{p_i}^{new} = \frac{\tau_{r_i} \hat{a}_i}{\tau_{r_i} - \tau_{a_i}}$

14. $\Delta_{p_i} = \frac{1}{\tau_{p_i}^{new}} - \frac{1}{\tau_{p_i}}$

15. $\Psi_{p_i} = \frac{p_i^{new}}{\tau_{p_i}^{new}} - \frac{p_i}{\tau_{p_i}}$

16. [Update the posterior approximation]

17. $C_m = h_C(C_m, e_i, \Delta_{p_i})$

18. $m = h_m(m, C_m, e_i, \Delta_{p_i}, \Psi_{p_i})$

19. until All i-s have been updated

20. $p = p^{new}$

21. $\tau_p = \tau_p^{new}$

22. until Convergence

### D. Implementation Details

For the practical implementation with finite-dimensional $H$, we suggest incorporating small enhancements to reVAMP, as discussed in Algorithm 1.

Firstly, it is advisable to clip the variances $\tau_{r_i}$ and $\tau_{p_i}^{new}$ within a positive interval $[\gamma_{min}, \gamma_{max}]$. Occasionally, the re-VAMP algorithm may yield negative values for $\tau_{r_i}$ and $\tau_{p_i}^{new}$ if not addressed. In our numerical results presented in Section VI, we utilized $\gamma_{min} = 10^{-10}$ and $\gamma_{max} = 10^{10}$.

Secondly, rather than mandating reVAMP to complete several iterations, we propose stopping the iterations when the normalized difference $||m^{new} - m||/||m^{new}||$ falls below a tolerance threshold $\epsilon$. In Section VI, our numerical results employed $\epsilon = 10^{-4}$.

Lastly, it's important to note that the reVAMP algorithm requires the user to initialize $p, \tau_p, m$ and $C_m$. Generally, all elements of $\tau_p$ and all diagonal elements of $C_m$ must be positive; other initializations do not significantly affect the final result. In our experiments in Section VI, we set all elements in $m, p$, and $\tau_p$ to be 1, and $C_m$ to be an identity matrix.

## V. RANGING ESTIMATION WITH NAKAGAMI-M PRIOR DISTRIBUTION

### A. MPCDPD-based Ranging Method

We propose the EM-ReVAMP algorithm, outlined in Algorithm 2, for estimating $d_0$. This algorithm utilizes ReVAMP
sequentially at each step of the EM algorithm to obtain approximate second-order moments. Specifically, within the ReVAMP part, with given $\Omega_i$ in (10) w.r.t. $\hat{d}_0$ and prior distribution in (11), the marginal posterior approximation involves the calculation of $\hat{a}_i$ and $\tau_{a_i}$ in (35a) and (35b) can be computed as follows:

$$ \hat{a}_i = \frac{m\Omega_i\tau_i \frac{1}{m\tau_{r_i} + \Omega_i} \frac{F_1(m+1; 2; \frac{\Omega_i|r_i|^2}{m\tau_{r_i}^2 + \tau_{r_i}\Omega_i})}{1F_1(m; 1; \frac{\Omega_i|r_i|^2}{m\tau_{r_i}^2 + \tau_{r_i}\Omega_i})}; \quad (43a) $$

$$ \tau_{a_i} = \frac{m\Omega_i\tau_{r_i} \frac{1}{m\tau_{r_i} + \Omega_i} \frac{F_1(m+1; 1; \frac{\Omega_i|r_i|^2}{m\tau_{r_i}^2 + \tau_{r_i}\Omega_i})}{1F_1(m; 1; \frac{\Omega_i|r_i|^2}{m\tau_{r_i}^2 + \tau_{r_i}\Omega_i})} - \hat{a}_i\hat{a}_i^*, \quad (43b) $$

where $F_1(a; b; z)$ represents the confluent hypergeometric function [48], defined by the hypergeometric series:

$$ F_1(a; b; z) = \sum_{k=0}^{\infty} \frac{(a)_k z^k}{(b)_k k!}. \quad (44) $$

The detailed derivation is provided in Appendix.

In Algorithm 2, for accelerating convergence, we typically recommend users to initialize $\hat{d}_0$ based on the actual characteristics of the environment, such as the maximum range or distance resolution. In our subsequent experiments, we initialized it to 0.1 meter. Additionally, we set the EM-loop to terminate when the difference of $d_0$ before and after the iteration is less than 0.1 meters for greater precision and maximum iteration times to be 20. When transitioning from the Nakagami-m fading model to another statistical model, the EM-reVAMP algorithm can be adapted by primarily modifying (12), (13), and (26). These adjustments will enable the algorithm to effectively accommodate the characteristics of the new statistical model while maintaining its functionality.

## B. Theoretical Cramér-Rao Bound (CRB) For $m = 1$

When $m = 1$, the path complex attenuation coefficients $a \in C^{L \times 1}$ that each element $a_i$ ($i = 0, \dots, L-1$) is an i.i.d. complex zero-mean Gaussian random variable, can be expressed as follows:

$$ a \sim CN(0, C_{aa}), \quad C_{aa} = \begin{bmatrix} \Omega_0(d_0) & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \Omega_{L-1}(d_K) \end{bmatrix}. \quad (45) $$

To estimate $d_0$ directly and solely based on $y$ using the maximum likelihood estimator (MLE) of the pdf of $y$ given $\Omega(d_0)$, we proceed as follows:

$$ p(y|\Omega(d_0)) = \pi^{-N}(\det(C_{yy}))^{-1} \exp(-y^H C_{yy}^{-1} y) \quad (46) $$

where

$$ C_{yy} = HC_{aa}H^H + \sigma_v^2 I. \quad (47) $$

To compute the Fisher Information Matrix (FIM) from the pdf $p(y|\Omega(d_0))$, the FIM can be represented as follows:

$$ J_{d_0 d_0} = \text{tr} \left[ C_{yy}^{-1} \frac{\partial C_{yy}}{\partial d_0} C_{yy}^{-1} \frac{\partial C_{yy}}{\partial d_0} \right], \quad (48) $$

## Algorithm 2 EM-ReVAMP

**Ensure:** $\hat{d}_0$

**Require:** $y, H, p_v(v), m, [\tau_0, \dots, \tau_{L-1}], G_0, n$

1: Initialize: $\hat{d}_0$

2: repeat [For $t = 0 \dots L-1$]

3: Initialize: $m, C_m, p, \tau_p$

4: Update $\Omega(\hat{d}_0)$ w.r.t. $G_0, n, \tau_i$ and $\hat{d}_0$ with (10)

5: repeat

6: repeat [For each $i = 0 \dots L-1$]

7: [Update the extrinsic]

8: $\tau_m = \text{diag}(C_m)$

9: $r_i = \frac{\tau_{p_i} m_i - \tau_{m_i} p_i}{\tau_{p_i} - \tau_{m_i}}$

10: $\tau_{r_i} = \frac{\tau_{m_i} \tau_{p_i}}{\tau_{p_i} - \tau_{m_i}}$

11: [Approximate the marginal posterior]

12: Update $\hat{a}_i$ with (43a)

13: Update $\tau_{a_i}$ with (43b)

14: [Propagate the approximation back]

15: $p_i^{\text{new}} = \frac{\tau_{r_i} \hat{a}_i - \tau_{a_i} r_i}{\tau_{r_i} - \tau_{a_i}}$

16: $\tau_{p_i}^{\text{new}} = \frac{\tau_{r_i} a_i}{\tau_{r_i} - \tau_{a_i}}$

17: $\Delta_{p_i} = \frac{1}{\tau_{p_i}^{\text{new}} - \tau_{p_i}}$

18: $\Psi_{p_i} = \frac{p_i^{\text{new}}}{\tau_{p_i}^{\text{new}} - \tau_{p_i}}$

19: [Update the posterior approximation]

20: $C_m = h_C(C_m, e_i, \Delta_{p_i})$

21: $m = h_m(m, C_m, e_i, \Delta_{p_i}, \Psi_{p_i})$

22: until All $i$-s have been updated

23: $p = p^{\text{new}}$

24: $\tau_p = \tau_p^{\text{new}}$

25: until Convergence

26: $\hat{d}_0 = \arg \min_{d_0} \sum_{i=0}^{L-1} \left[ \ln \Omega_i(d_0) + \frac{\tau_{m_i} + |m_i|^2}{\Omega_i(d_0)} \right]$

27: until Convergence

where $\text{tr}$ denotes the trace operator. After some algebraic computations, we obtain the expression in (48)

$$ \frac{\partial C_{yy}}{\partial d_0} = H \frac{\partial C_{aa}}{\partial d_0} H^H, \quad (49a) $$

$$ \frac{\partial C_{aa}}{\partial d_0} = -nG_0 \begin{bmatrix} (d_0 + c\tau_0)^{-n-1} & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & d_0 + c\tau_{K-1}^{-n-1} \end{bmatrix}. \quad (49b) $$

In conclusion, using (48) and (49), the CRB of $d_0$ w.r.t. MPCDP-based ranging can be calculated as follows:

$$ CRB_{d_0} = J_{d_0 d_0}^{-1}. \quad (50) $$

Unfortunately, if $m \neq 1$, the high-dimensional integration required to obtain the likelihood $p(y|d_0)$ is intractable, making it impossible to calculate its CRB.

## VI. SIMULATION RESULTS

This section presents the simulation verification using MATLAB to assess the impact of different parameters. Table I lists the main parameters involved. In general, the primary environmental factors influencing our MPCDP-based ranging method are the number of distinguishable MPCs, the
<table><caption>TABLE I
PARAMETERS SETTING</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>G<sub>0</sub></td><td>1.</td></tr><tr><td>SNR (dB)</td><td>Range from 0 to 20, default 10.</td></tr><tr><td>N</td><td>20.</td></tr><tr><td>L</td><td>Ranging from 2 to 8.</td></tr><tr><td>n</td><td>Ranging from 2 to 4, default 3.</td></tr><tr><td>m</td><td>Ranging from 1 to 10, default 5.</td></tr><tr><td>d<sub>0</sub> (meter)</td><td>20.</td></tr><tr><td>Distance of NLOS path (m)</td><td>Random between 1.1d<sub>0</sub> to 2.0d<sub>0</sub>.</td></tr><tr><td>Test repetitions</td><td>50.</td></tr></tbody></table>

magnitude of Signal-to-Noise Ratio (SNR), the propagation attenuation factor $n$, and the shape parameter of Nakagami-m distribution $m$. In the following subsections, we analyze the effects of these factors on ranging accuracy through simulations. We conduct 10000 times for each scenario and calculate the Root Mean Square Error (RMSE). For $m = 1$, we calculate the square root of CRB (SR-CRB) for our MPCPDP-based ranging method. Moreover, we compare with the SoTA RSS-based ranging method in [49]. In the simulation result figures, we use a solid line to denote the RMSE of the RSS-based ranging method, a dashed line to denote the RMSE of the MPCPDP-based ranging method, and a dotted line to denote the SR-CRB.

### A. Impact of SNR and Number of NLoS Paths

In this set of experiments, we set $n = 3$ and examine the influence of SNR and the number of NLoS paths on the estimation bias. Figs. 3 and 4 illustrate the simulation results for $m = 1$ and $m = 5$, respectively. We observe that varying SNR from 15dB to 40dB does not significantly affect our performance. Moreover, the estimation accuracy of our method gradually improves with an increasing number of NLoS paths which act better than the RSS-based ranging method under the same conditions. Even for $m \neq 1$, it is hard to get its theoretic CRB, comparison to RSS-based ranging method can show our method's high precision.

### B. Impact of the propagation fading factor $n$ and Number of NLoS Paths

This set of experiments investigates the effects of the environmental propagation fading factor $n$ and the number of NLoS paths on the estimation bias. The simulation results, depicted in Figs. 5 and 6 for $m = 1$ and $m = 5$, respectively, indicate that the variation of $n$ within the range of 2 to 4 has a obvious impact on our performance. As $n$ increases, the variance of Nakagami-m actually decreases, thus leading to an increase in estimation accuracy for both the RSS-based ranging method and our method. This is theoretically verified by the SR-CRB when $m = 1$. It is clear to see that our method consistently performs better than the RSS-based ranging method for different $n$.

### C. Impact of the shape parameter $m$ and Number of NLoS Paths

In this simulation, we set $n$ as 3, SNR as 10 dB, $d_0$ as 20 meters, and we examine the influence of the shape param-

<figure><img src="image_4.png" /><figcaption>Fig. 3. The impact of SNR and the number of NLoS paths on $d_0$ estimation with $m = 1$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. The impact of SNR and the number of NLoS paths on $d_0$ estimation with $m = 5$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. The impact of the environment propagation fading factor $n$ and the number of NLoS paths on $d_0$ estimation with $m = 1$</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. The impact of the environment propagation fading factor $n$ and the number of NLoS paths on $d_0$ estimation with $m = 5$</figcaption></figure>

eter of the Nakagami-m distribution, denoted as $m$, and the number of NLoS paths on the estimation bias. Fig. 7 presents the simulation results, indicating that varying $m$ from 1 to 10 significantly affects the performance of MPCPDP-based ranging method. As $m$ increases, the variance of Nakagami-m distribution decreases and the accuracy of estimation method gets higher. Moreover, with the number of NLoS paths increase, the performance of our method also increases but RSS-based ranging generally decreases. Obviously, under the same condition, the method proposed in this paper has the small range error compared to the RSS-based range method. In addition, the range error of the new method increases along with the number of NLoS paths increases.

<figure><img src="image_8.png" /><figcaption>Fig. 7. The impact of the Nakagami-m distribution's shape parameter $m$ and the number of NLoS paths on $d_0$ estimation with $n = 3$</figcaption></figure>

### D. Impact of number of NLoS and the mismatched shape parameter $m$

We investigate the impact of a mismatched shape parameter $m$ in the Nakagami-m distribution, where we set a true value of $m = 5$, SNR = 10 and $n = 3$. Our simulations reveal that the estimation error is relatively small when the deviation of $m$ is not substantial. This observation underscores the robustness of our algorithm under theoretical conditions. However, achieving accurate parameter initialization in practical scenarios is

<figure><img src="image_9.png" /><figcaption>Fig. 8. The impact of the mismatched Nakagami-m distribution's shape parameter $m$ and the number of NLoS paths on $d_0$ estimation with true $m = 5$</figcaption></figure>

challenging, and significant errors in initialization can lead to substantial estimation deviations. One approach to address this challenge is to jointly estimate all parameters using the EM algorithm, which remains a topic for our future research. While theoretically effective of the EM algorithm, this method encounters difficulties due to the complex nonlinear relationships among parameters and non-convex optimization problems. Moreover, increasing the number of parameters to estimate can degrade performance due to constraints imposed by available observed data. Another promising approach that warrants further research attention is minimizing errors resulting from model mismatches.

### E. Simulation Conclusions

Based on experimental simulations with varying SNR values, $n$ and $m$, and the number of NLoS paths, our method has demonstrated strong performance in diverse and complex environments. Notably, the number of NLoS paths has emerged as a crucial parameter, significantly influencing estimation accuracy. As the number of NLoS paths increases, our algorithm's accuracy improves, whereas the performance of the RSS-based method declines. Furthermore, as $m$ or $n$ increases, the variance of the fading channel decreases, thereby enhancing our method's performance. Consequently, in complex environments characterized by significant fading with large $m$ and $n$ and multiple NLoS paths, the MPCPDP-based ranging algorithm proves to be a more effective solution than RSS-based method for estimating the LoS distance.

## VII. CONCLUSIONS

In conclusion, we propose a novel MPCPDP-based ranging method, aimed at exploiting multipath effect. To theoretically validate the feasibility of our approach, we selected the widely adopted Nakagami-m statistical model and established a relationship between distribution parameters and propagation distance. On one hand, this model ensures its applicability across the majority of cases; on the other hand, it exemplifies
many other statistical models that traditional estimation tools
cannot approach. To address the challenges of ranging estima-
tion, which involve hidden random variables and intractable
posterior distributions due to complex statistical models such
as the Nakagami-m model, we introduce the EM-ReVAMP
algorithm. The simulation results convincingly demonstrate
the effectiveness of our approach than the RSS-based ranging
method, providing substantial evidence to support the accuracy
and robustness of our approach. Moreover, the EM-ReVAMP
algorithm can be adapted for other statistical fading models
with minor modifications, assuming that the statistical models
and their parameter initializations are sufficiently accurate.
To further validate the practicality and effectiveness of our
method, our next objective is to collect measurement data from
diverse environments and conduct comprehensive experimen-
tal analysis. This pivotal step will enable us to assess the per-
formance of our method in real-world scenarios. Additionally,
exploring further application scenarios of EM-ReVAMP and
investigating its theoretical performance are essential aspects
that warrant attention. Finally, the robustness and accuracy of
this algorithm need to undergo further testing.

APPENDIX

In this appendix, we will derive the expressions for $\hat{a}_i$ and $\tau_{a_i}$ in (35a) and (35b) based on the pdf $p_{a_i}(a_i|\Omega_i)$. In the EM-ReVAMP algorithm presented in Algorithm 2, we have the following expressions:

$$
\hat{a}_i = \frac{\int a_i p_{a_i}(a_i | \Omega_i) \mathcal{CN}(a_i; r_i, \tau_{r_i}) da_i}{\int p_{a_i}(a_i | \Omega_i) \mathcal{CN}(a_i; r_i, \tau_{r_i}) da_i}; \quad (51a)
$$

$$
\tau_{a_i} = \frac{\int |a_i - \hat{a}_i|^2 p_{a_i}(a_i | \Omega_i) \mathcal{CN}(a_i; r_i, \tau_{r_i}) da_i}{\int p_{a_i}(a_i | \Omega_i) \mathcal{CN}(a_i; r_i, \tau_{r_i}) da_i}, \quad (51b)
$$

where $Z_i$ is defined as:

$$
Z_i = \int p_{a_i}(a_i | \Omega_i) \mathcal{CN}(a_i; r_i, \tau_{r_i}) da_i. \quad (52)
$$

For calculating $Z_i$, we have the integral:

$$
Z_i = \int \frac{m^m |a_i|^{2m-2}}{\pi \Gamma(m) \Omega_i^m} \exp \left[ -\frac{|m| a_i|^2}{\Omega_i} \right] \\
\times \frac{1}{\pi \tau_{r_i}} \exp \left( -\frac{(a_i - r_i)^* (a_i - r_i)}{\tau_{r_i}} \right) da_i. \tag{53}
$$

Inside the integral, we have:

$$
\int_0^{2\pi} \exp \left( \frac{2}{\tau_{r_i}} \left( |r_i| r \cos(\phi - \phi_{r_i}) \right) \right) d\phi = 2\pi J_0 \left( j \frac{2}{\tau_{r_i}} r |r_i| \right), \quad (54)
$$

where $J_0(x)$ is the Bessel function of the first kind. By incorporating the equations and clarifying the expressions, the derivation of $Z_i$ becomes more comprehensible. Moreover, we can simplify the remaining integral using the confluent hypergeometric function:

$$
\begin{align}
& \int_0^{+\infty} x^\mu \exp(-\alpha x^2) J_v(\beta x) dx \nonumber \\
& = \frac{\beta^v \Gamma(0.5(v+\mu+1))}{2^{v+1}\alpha^{0.5(v+\mu+1)}\Gamma(v+1)} {}_1F_1\left(0.5(v+\mu+1); v+1; -\frac{\beta^2}{4\alpha}\right), \tag{55}
\end{align}
$$

where $_1F_1(a; b; z)$ is the confluent hypergeometric function. By utilizing these results, we can evaluate $Z_i$ and proceed with the derivation of $\hat{a}_i$ and $\tau_{a_i}$. Simplifying the integral and applying the properties of the confluent hypergeometric function, we obtain:

$$
\begin{align}
& \int_0^{+\infty} r^{2m-1} \exp\left(-\left(\frac{m}{\Omega_i} + \frac{1}{\tau_{r_i}}\right) r^2\right) J_0\left(j \frac{2}{\tau_{r_i}} r |r_i|\right) dr \nonumber \\
& = \frac{\Gamma(m)}{2 \left(\frac{m}{\Omega_i} + \frac{1}{\tau_{r_i}}\right)^m} {}_1F_1(m; 1; \frac{\Omega_i |r_i|^2}{m \tau_{r_i}^2 + \tau_{r_i} \Omega_i}). \tag{56}
\end{align}
$$

Based on the previous results, we can express the $Z_i$ as
follows:

$$
Z_i = \frac{m^m}{\pi \Omega_i^m \tau_{r_i} \left(\frac{m}{\Omega_i} + \frac{1}{\tau_{r_i}}\right)^m} {}_1F_1(m; 1; \frac{\Omega_i |r_i|^2}{m \tau_{r_i}^2 + \tau_{r_i} \Omega_i}). \quad (57)
$$

Then for $\int a_i p_i(a_i) CN(a_i; r_i, \tau_{r_i}) da_i$, it can be written as:

$$
\begin{align}
& \int a_i p_i(a_i) \mathcal{CN}(a_i; r_i, \tau_{r_i}) da_i \nonumber \\
&= \frac{m^m}{\pi^2 \Gamma(m) \Omega_i^m \tau_{r_i}} \int_0^{+\infty} \int_0^{2\pi} (\cos \phi + j \sin \phi) \nonumber \\
&\qquad \exp \left\{ \frac{2}{\tau_{r_i}} [|r_i| r \cos(\phi - \phi_{r_i})] \right\} d\phi r^{2m-1} \exp \left[ -\left( \frac{m}{\Omega_i} + \frac{1}{\tau_{r_i}} \right) r^2 \right] rdr. \tag{58}
\end{align}
$$

Inside this expression, we have the following intermediate
results:

$$
\begin{align}
& \int_0^{2\pi} (\cos \phi + j \sin \phi) \exp\left(\frac{2}{\tau_{r_i}}\left(|r_i| r \cos(\phi - \phi_{r_i})\right)\right) d\phi \tag{59} \\
& \quad = -j2\pi \exp(j\phi_{r_i}) J_1\left(j\frac{2r|r_i|}{\tau_{r_i}}\right), \nonumber
\end{align}
$$

and

$$
\begin{align}
& \int_0^{+\infty} r^{2m} \exp\left(-\left(\frac{m}{\Omega_i} + \frac{1}{\tau_{r_i}}\right) r^2\right) J_1\left(j \frac{2}{\tau_{r_i}} r |r_i|\right) dr \nonumber \\
&= \frac{j |r_i| \Gamma(m+1)}{2 \tau_{r_i} \left(\frac{m}{\Omega_i} + \frac{1}{\tau_{r_i}}\right)^{(m+1)}} {}_1F_1(m+1; 2; \frac{\Omega_i |r_i|^2}{m \tau_{r_i}^2 + \tau_{r_i} \Omega_i}). \tag{60}
\end{align}
$$

Thus, $\int a_i p_i(a_i) CN(a_i; r_i, \tau_{r_i}) da_i$ can be expressed as follows:

$$
\int a_i p_i(a_i) \mathcal{CN}(a_i; r_i, \tau_{r_i}) da_i = \frac{m^{m+1} |r_i|}{\pi \Omega_i^m \tau_{r_i}^2 \left(\frac{m}{\Omega_i} + \frac{1}{\tau_{r_i}}\right)^{(m+1)}} \\
\times \exp(j\phi_{r_i}) {}_1F_1(m+1; 2; \frac{\Omega_i |r_i|^2}{m\tau_{r_i}^2 + \tau_{r_i}\Omega_i}). \tag{61}
$$

For $\int a_i^* a_i p_i(a_i) CN(a_i; r_i, \tau_{r_i}) da_i$, we can compute it as
follows:

$$
\int a_i^* a_i p_i(a_i) \mathcal{CN}(a_i; r_i, \tau_{r_i}) da_i = \frac{m^m}{\pi^2 \Gamma(m) \Omega_i^m \tau_{r_i}} \int_0^{+\infty} \int_0^{2\pi} \\
\exp\left\{\frac{2}{\tau_{r_i}} \left[|r_i| r \cos(\phi - \phi_{r_i})\right]\right\} d\phi r^{2m} \exp\left[-\left(\frac{m}{\Omega_i} + \frac{1}{\tau_{r_i}}\right) r^2\right] rdr. \quad (62)
$$
Inside this expression, we have the following intermediate results:

$$ \int_{0}^{2\pi} \exp\left(\frac{2}{\tau_{r_i}}\left(|r_i| r \cos(\phi - \phi_{r_i})\right)\right) d\phi = 2\pi J_0\left(j\frac{2}{\tau_{r_i}}r|r_i|\right), \quad (63) $$

and

$$ \int_{0}^{+\infty} r^{2m+1} \exp\left(-\left(\frac{m}{\Omega_i} + \frac{1}{\tau_{r_i}}\right)r^2\right) J_0\left(j\frac{2}{\tau_{r_i}}r|r_i|\right) dr = \frac{\Gamma(m+1)}{2\left(\frac{m}{\Omega_i} + \frac{1}{\tau_{r_i}}\right)^{(m+1)}} {}_1F_1(m+1; 1; \frac{\Omega_i|r_i|^2}{m\tau_{r_i}^2 + \tau_{r_i}\Omega_i}). \quad (64) $$

Therefore, $\int a_i a_i^* p_i(a_i) CN(a_i; r_i, \tau_{r_i}) da_i$ can be given as :

$$ \int a_i a_i^* p_i(a_i) CN(a_i; r_i, \tau_{r_i}) da_i = \frac{m^{(m+1)}}{\pi \Omega_i^m \tau_{r_i} (\frac{m}{\Omega_i} + \frac{1}{\tau_{r_i}})^{(m+1)}} {}_1F_1(m+1; 1; \frac{\Omega_i|r_i|^2}{m\tau_{r_i}^2 + \tau_{r_i}\Omega_i}). \quad (65) $$

Thus, we can obtain the expressions for $\mathcal{E}[a_i]$ and $\mathcal{E}[a_i a_i^*]$ as follows:

$$ \mathcal{E}[a_i] = \frac{m \Omega_i r_i}{m \tau_{r_i} + \Omega_i} \frac{1 F_1(m+1; 2; \frac{\Omega_i |r_i|^2}{m \tau_{r_i}^2 + \tau_{r_i} d_0})}{1 F_1(m; 1; \frac{\Omega_i |r_i|^2}{m \tau_{r_i}^2 + \tau_{r_i} \Omega_i})}; \quad (66a) $$

$$ \mathcal{E}[a_i a_i^*] = \frac{m \Omega_i \tau_{r_i}}{m \tau_{r_i} + \Omega_i} \frac{1 F_1(m+1; 1; \frac{\Omega_i |r_i|^2}{m \tau_{r_i}^2 + \tau_{r_i} \Omega_i})}{1 F_1(m; 1; \frac{\Omega_i |r_i|^2}{m \tau_{r_i}^2 + \tau_{r_i} \Omega_i})}. \quad (66b) $$

Finally, we obtain the expressions for $\hat{a}_i$ and $\tau_{a_i}$ as follows:

$$ \hat{a}_i = \frac{m \Omega_i r_i}{m \tau_{r_i} + \Omega_i} \frac{1 F_1(m+1; 2; \frac{\Omega_i |r_i|^2}{m \tau_{r_i}^2 + \tau_{r_i} \Omega_i})}{1 F_1(m; 1; \frac{\Omega_i |r_i|^2}{m \tau_{r_i}^2 + \tau_{r_i} \Omega_i})}; \quad (67a) $$

$$ \tau_{a_i} = \frac{m \Omega_i \tau_{r_i}}{m \tau_{r_i} + \Omega_i} \frac{1 F_1(m+1; 1; \frac{\Omega_i |r_i|^2}{m \tau_{r_i}^2 + \tau_{r_i} \Omega_i})}{1 F_1(m; 1; \frac{\Omega_i |r_i|^2}{m \tau_{r_i}^2 + \tau_{r_i} \Omega_i})} - \hat{a}_i \hat{a}_i^*, \quad (67b) $$

which are corresponding to (43a) and (43b).

## REFERENCES

[1] N. Shaik and P. K. Malik, “A comprehensive survey 5G wireless communication systems: open issues, research challenges, channel estimation, multi carrier modulation and 5G applications,” *Multimedia Tools and Applications*, vol. 80, no. 19, pp. 28 789–28 827, 2021.

[2] F. Mogyorósi, P. Revisnyei, A. Pašić, Z. Papp, I. Törös, P. Varga, and A. Pašić, “Positioning in 5G and 6G networks—A Survey,” *Sensors*, vol. 22, no. 13, p. 4757, 2022.

[3] M. J. Shehab, I. Kassem, A. A. Kutty, M. Kucukvar, N. Onat, and T. Khattab, “5G networks towards smart and sustainable cities: A review of recent developments, applications and future perspectives,” *IEEE Access*, vol. 10, pp. 2987–3006, 2021.

[4] A. Gohar and G. Nencioni, “The role of 5G technologies in a smart city: The case for intelligent transportation system,” *Sustainability*, vol. 13, no. 9, p. 5188, 2021.

[5] S. Li, L. Da Xu, and S. Zhao, “5G Internet of Things: A survey,” *J. Ind. Inf. Integr.*, vol. 10, pp. 1–9, 2018.

[6] D. Dardari, A. Conti, U. Ferner, A. Giorgetti, and M. Z. Win, “Ranging With Ultrawide Bandwidth Signals in Multipath Environments,” *Proc. IEEE*, vol. 97, no. 2, pp. 404–426, Feb 2009.

[7] H. Wymeersch, S. Marano, W. M. Gifford, and M. Z. Win, “A Machine Learning Approach to Ranging Error Mitigation for UWB Localization,” *IEEE Transactions on Communications*, vol. 60, no. 6, pp. 1719–1728, June 2012.

[8] O. Bialer, D. Raphaeli, and A. J. Weiss, “Maximum-Likelihood Direct Position Estimation in Dense Multipath,” *IEEE Transactions on Vehicular Technology*, vol. 62, no. 5, pp. 2069–2079, Jun 2013.

[9] T. Kazaz, G. J. M. Janssen, J. Romme, and A. J. van der Veen, “Delay Estimation for Ranging and Localization Using Multiband Channel State Information,” *IEEE Trans. Wirel. Commun.*, vol. 21, no. 4, pp. 2591–2607, April 2022.

[10] A. Venus, E. Leitinger, S. Tertinek, and K. Witrisal, “A Graph-based Algorithm for Robust Sequential Localization Exploiting Multipath for Obstructed-LOS-Bias Mitigation,” *IEEE Trans. Wirel. Commun.*.

[11] A. Venus, E. Leitinger, S. Tertinek, F. Meyer, and K. Witrisal, “Graph-based Simultaneous Localization and Bias Tracking for Robust Positioning in Obstructed LOS Situations,” in *56th Asilomar Conf. on Signals, Systems, and Computers*, Pacific Grove, CA, USA, 2022, pp. 777–784.

[12] X. Li, E. Leitinger, A. Venus, and F. Tufvesson, “Sequential Detection and Estimation of Multipath Channel Parameters Using Belief Propagation,” *IEEE Trans. Wirel. Commun.*, vol. 21, no. 10, p. 8385–8402, 2022.

[13] E. Leitinger, P. Meissner, C. Rüdisser, G. Dumphart, and K. Witrisal, “Evaluation of Position-Related Information in Multipath Components for Indoor Positioning,” *IEEE J. Sel. Areas Commun.*, vol. 33, no. 11, pp. 2313–2328, Nov 2015.

[14] C. Gentner, T. Jost, W. Wang, S. Zhang, A. Dammann, and U. C. Fiebig, “Multipath Assisted Positioning with Simultaneous Localization and Mapping,” *IEEE Trans. Wirel. Commun.*, vol. 15, no. 9, pp. 6104–6117, Sept 2016.

[15] E. Leitinger, F. Meyer, F. Hlawatsch, K. Witrisal, F. Tufvesson, and M. Z. Win, “A Belief Propagation Algorithm for Multipath-Based SLAM,” *IEEE Trans. Wirel. Commun.*, vol. 18, no. 12, pp. 5613–5629, Dec 2019.

[16] E. Leitinger, S. Grebien, and K. Witrisal, “Multipath-Based SLAM Exploiting AoA and Amplitude Information,” in *IEEE Int. Conf. on Commu. Workshops (ICC Workshops)*, Shanghai, China, 2019, pp. 1–7.

[17] H. Kim, K. Granström, L. Svensson, S. Kim, and H. Wymeersch, “PMBM-Based SLAM Filters in 5G mmWave Vehicular Networks,” *IEEE Trans. Veh. Technol.*, vol. 71, no. 8, pp. 8646–8661, Aug 2022.

[18] K. Papakonstantinou and D. Slock, “Direct Location Estimation for MIMO Systems in Multipath Environments,” in *IEEE GLOBECOM*, 2008, pp. 1–5.

[19] ——, “Hybrid TOA/AOD/Doppler-Shift localization algorithm for NLOS environments,” in *IEEE 20th Int. Symp. Person. Indoor Mobile Radio Commun. (PIMRC)*, 2009, pp. 1948–1952.

[20] G. Zanca, F. Zorzi, A. Zanella, and M. Zorzi, “Experimental comparison of RSSI-based localization algorithms for indoor wireless sensor networks,” in *Proc. Workshop on Real-world wireless sensor networks*, 2008, pp. 1–5.

[21] H. P. Mistry and N. H. Mistry, “RSSI based localization scheme in wireless sensor networks: A survey,” in *Proc. IEEE 51th Int. Conf. on Advanced Computing and Communication Technologies*, pp. 647–652.

[22] Z. Yang, Z. Zhou, and Y. Liu, “From RSSI to CSI: Indoor localization via channel response,” *ACM Comput. Surv. (CSUR)*, vol. 46, no. 2, pp. 1–32, 2013.

[23] T. He, J. Tan, W. Zhuo, M. Printz, and S.-H. G. Chan, “Tackling Multipath and Biased Training Data for IMU-Assisted BLE Proximity Detection,” in *IEEE Conf. on Comput. Commun. (INFOCOM)*, 2022, pp. 1259–1268.

[24] T. Nowak, M. Hartmann, H.-M. Tröger, L. Patino-Studencki, and J. Thieleke, “Probabilistic Multipath Mitigation in RSSI-Based Direction-of-Arrival Estimation,” in *IEEE Int. Conf. on Commun. Workshops (ICC Workshops)*, 2017, pp. 1024–1029.

[25] T. Aono, K. Higuchi, T. Ohira, B. Komiyama, and H. Sasaoka, “Wireless Secret Key Generation Exploiting Reactance-Domain Scalar Response of Multipath Fading Channels,” *IEEE Trans. Antennas Propag.*, vol. 53, no. 11, pp. 3776–3784, 2005.

[26] M. Triki, D. T. M. Slock, V. Rigal, and P. Francois, “Mobile Terminal Positioning via Power Delay Profile Fingerprinting: Reproducible Validation Simulations,” in *IEEE Veh. Tech. Conf.*, 2006, pp. 1–5.

[27] M. Triki and D. T. Slock, “Mobile Localization for NLOS Propagation,” in *IEEE 18th Int. Symp. Person. Indoor Mobile Radio Commun. (PIMRC)*, 2007, pp. 1–4.

[28] V. B. Vales, T. Dominguez-Bolano, C. J. Escudero, and J. A. Garcia-Naya, “Using the power delay profile to accelerate the training of neural network-based classifiers for the identification of LOS and NLOS UWB propagation conditions,” *IEEE Access*, vol. 8, pp. 220 205–220 214, 2020.

[29] F. Xiao and D. Slock, “A Cramér–Rao Bound for Indoor Power Delay Profile Based Ranging,” 2023.
[30] A. A. Saleh and R. Valenzuela, "A statistical model for indoor multipath propagation," *IEEE J. Sel. Areas Commun.*, vol. 5, no. 2, pp. 128-137, 1987.

[31] A. Meijerink and A. F. Molisch, "On the physical interpretation of the Saleh-Valenzuela model and the definition of its power delay profiles," *IEEE Trans. Antennas Propag.*, vol. 62, no. 9, pp. 4780-4793, 2014.

[32] A. Medeisis and A. Kajackas, "On the use of the universal Okumura-Hata propagation prediction model in rural areas," in *IEEE 51st Veh. Technol. Conf. (VTC)*, vol. 3, 2000, pp. 1815-1818.

[33] 3rd Generation Partnership Project (3GPP), "Technical Specification Group Radio Access Network," 3GPP Technical Specification, 2001, 3GPP TR 25.882 V1.1.0, Release 5. [Online]. Available: https://www.3gpp.org/ftp/tsg_ran/tsg_ran/TSGR_14/Docs/PDFs/RP-010810.pdf

[34] H. Z. Stefanovic, D. N. Milic, D. C. Stefanovic, and S. Z. Milosavljevic, "Some integral characteristics of MRC receiver in Nakagami-m fading environment," *Int. J. Reason.-based Intell. Syst.*, vol. 5, no. 3, pp. 183-188, 2013.

[35] S. M. Abuelenin, "On the similarity between Nakagami-m fading distribution and the Gaussian ensembles of random matrix theory," *arXiv preprint arXiv:1803.08688*, 2018.

[36] H. Stefanovic, A. Savic et al., "Some general characteristics of Nakagami-m distribution," in *Proc. ISICM Conf.*, 2011, pp. 721-731.

[37] N. C. Beaulieu and C. Cheng, "Efficient Nakagami-m fading channel simulation," *2005 IEEE Trans. Veh. Technol.*, vol. 54, no. 2, pp. 413-424.

[38] Z. Zhao, F. Xiao, and D. Slock, "Approximate Message Passing for Not So Large niid Generalized Linear Models," in *2023 IEEE 24th Int. Workshop on Signal Processing Advances in Wire. Commu. (SPAWC)*, pp. 386-390.

[39] S. Rangan, P. Schniter, and A. K. Fletcher, "Vector Approximate Message Passing," 2018.

[40] A. M. Graff, W. N. Blount, P. A. Iannucci, J. G. Andrews, and T. E. Humphreys, "Analysis of OFDM signals for ranging and communications," in *Proc. 34th Int. Tech. Meet. Satellite Div. Inst. Navig. (ION GNSS+ 2021)*, pp. 2910-2924.

[41] A. Goldsmith, *Wireless communications*. Cambridge University Press, 2005.

[42] J. Liu, "Be cautious when using the FIR channel model with the OFDM-based communication systems," *IEEE Trans. Veh. Technol.*, vol. 58, no. 3, pp. 1607-1612, 2008.

[43] V. A. Aalo, C. Mukasa, and G. P. Efhymoglou, "Effect of mobility on the outage and BER performances of digital transmissions over Nakagami-m fading channels," *IEEE Trans. Veh. Technol.*, vol. 65, no. 4, pp. 2715-2721, 2015.

[44] L. Mirsky, *An introduction to linear algebra*. Courier Corporation, 2012.

[45] S. Borman, "The expectation maximization algorithm-a short tutorial," *Submitted for publication*, vol. 41, 2004.

[46] D. R. Hunter and K. Lange, "A Tutorial on MM Algorithms," *Am. Stat.*, vol. 58, no. 1, pp. 30-37, 2004.

[47] T. P. Minka, "Expectation propagation for approximate Bayesian inference," *arXiv preprint arXiv:1301.2294*, 2013.

[48] H. Buchholz, *The confluent hypergeometric function: with special emphasis on its applications*. Springer, 2013, vol. 15.

[49] A. U. Ahmed, R. Arablouei, F. d. Hoog, B. Kusy, N. Bergmann, and R. Jurdak, "Nakagami-m distribution of RSSI in shadowing pathloss model for indoor localization," in *IEEE Sensors Appl. Symp. (SAS)*, 2019, pp. 1-6.

**Zilu Zhao** (Student Member, IEEE) received his B.E. degree in Mechatronics from Tongji University, Shanghai, China, in 2018, and his M.Sc. degree in Electrical Engineering from Friedrich-Alexander-Universität Erlangen-Nürnberg, Germany, in 2022, respectively. He is a coauthor of one IEEE SPAWC'23 Best Student Paper Award at and one IEEE CAMAD'23 Best Paper Award (Runner-up). He is currently pursuing a Ph.D. degree in the Department of Communication Systems at Eurecom, Biot, France. His research interests include statistical inference and cell-free 6G.

**Dirk T. M. Slock** (Fellow, IEEE) received the E.E. degree from Ghent University, Ghent, Belgium, in 1982, the M.S.E.E., M.S. degrees in statistics, and the Ph.D. degree in electrical engineering from Stanford University, Stanford, CA, USA, in 1986, 1989, and 1989, respectively. While with Stanford, he developed new fast recursive least-squares algorithms for adaptive filtering. During 1989-91, he was the Member of the Research Staff with Philips Research Laboratory Belgium. In 1991, he joined EURECOM, where he is currently the Professor and teaches statistical signal processing and signal processing techniques for wireless communications. He invented semi-blind channel estimation, chip equalizer correlator receiver used by 3G HSDPA mobile terminals, spatial multiplexing cyclic delay diversity (MIMO-CDD) now part of LTE, and his work led to the Single Antenna Interference Cancellation (SAIC) integrated in the GSM standard in 2006. Recent keywords are multi-cell multi-user (Massive) MIMO, imperfect CSIT, distributed resource allocation, variational and empirical Bayesian learning techniques, large random matrix analysis, audio source separation, location estimation and exploitation. He has supervised about 40 Ph.D. students, leading to an edited book and more than 500 papers. He cofounded the start-ups SigTone in 2000 (music signal processing products) and Nestwave in 2014 (Ultra Low-Power Indoor and Outdoor Mobile Positioning). He is the coauthor of two IEEE Globecom'98, one IEEE SIU'04, one IEEE SPAWC'05, one IEEE WPNC'16, one IEEE SPAWC'18 and one IEEE SPAWC'23 Best Student Paper Award, and a honorary mention (finalist in Best Student Paper contest) at IEEE SSP'05, IWAENC'06, IEEE Asilomar'06, and IEEE ICASSP'17. He was the Associate Editor for the IEEE-SP Transactions in 1994-96, and IEEE SIGNAL PROCESSING LETTERS in 2009-10. He was the General Chair of the IEEE-SP SPAWC'06 and IWAENC'14 workshops, and EUSIPCO'15. In 1984, he was the recipient of the Fulbright Scholarship for Stanford University, and in 1992, one Best Journal Paper Award from IEEE-SP and one from EURASIP. He was also the recipient of the URSI France medal, in 2018. He is a Fellow of EURASIP.

**Fangqing Xiao** (Student Member, IEEE) received his B.E. degree in Space Science and Technology from XIDIAN University, Xi'an, China, in 2019, and his M.Sc. degree in Automatic System from Polytech Nice Sophia of University Côte d'Azur, Biot, France, in 2022, respectively. He is a coauthor of one IEEE SPAWC'23 Best Student Paper Award and one IEEE CAMAD'23 Best Paper Award (Runner-up). He is currently pursuing a Ph.D. degree in the Department of Communication Systems at Eurecom, Biot, France. His research interests include statistical inference, 5G+ positioning, and integrated sensing and communication.