A Graph-Based Algorithm for Robust Sequential Localization Exploiting Multipath for Obstructed-LOS-Bias Mitigation

This paper presents a factor graph formulation and particle-based sum-product algorithm (SPA) for robust sequential localization in multipath-prone environments. The proposed algorithm jointly performs data association, sequential estimation of a mobile agent position, and adapts all relevant model parameters. We derive a novel non-uniform false alarm (FA) model that captures the delay and amplitude statistics of the multipath radio channel. This model enables the algorithm to indirectly exploit position-related information contained in the multipath components (MPCs) for the estimation of the agent position without using any prior information such as floorplan information or training data. Using simulated and real measurements in different channel conditions, we demonstrate that the algorithm can provide high-accuracy position estimates even in fully obstructed line-of-sight (OLOS) situations and show that the performance of our algorithm constantly attains the posterior Cramér-Rao lower bound (P-CRLB), facilitating the additional information contained in the presented FA model. The algorithm is shown to provide robust estimates in both, dense multipath channels as well as channels showing specular, resolved MPCs, significantly outperforming state-of-the-art radio-based localization methods.


I. INTRODUCTION
Localization of mobile agents using radio signals in environments such as indoor or urban territories is still a challenging task [1]- [4].These environments are characterized by strong multipath propagation and frequent obstructed line-of-sight (OLOS) situations, which can prevent the correct extraction of the line-of-sight (LOS) component (see Fig. 1).Radio channels resulting from multipath propagation are commonly represented as a superposition of a finite number of specular multipath components (MPCs) [5]- [8].However, cluttered environments with closely-spaced reflecting objects or with diffuse scatters (such as walls covered by shelves or irregular object shapes), along with the finite bandwidth of the measurement equipment, cause dense multipath propagation, which cannot be resolved into specular MPCs anymore [5], [9]- [11].
There exist many safety-and security-critical applications, such as autonomous driving [12], medical services [13], or The financial support by the Christian Doppler Research Association, the Austrian Federal Ministry for Digital and Economic Affairs and the National Foundation for Research, Technology and Development is gratefully acknowledged.

LOS partial OLOS full OLOS multipath
anchor agent track Fig. 1.A mobile agent is moving alongside the anchors on an example trajectory.Due to an obstacle, the LOS to all anchors is not always available.There occur partial and full OLOS situations.Multipath propagation may occur, but there is no prior information about the surrounding environment.keyless entry systems [14], where robustness of the position estimate 1 is of critical importance.

A. State-of-the-Art Methods
New localization and tracking approaches within the context of 5G localization [15] that take advantage of large measurement apertures as ultra-wideband (UWB) systems [6], [16] or mmWave systems [17] seek to mitigate the effect of multipath propagation [18] (commonly referred to as "NLOS propagation") and OLOS situations [8], [19], or even take advantage of MPCs by exploiting inherent position information, turning multipath from impairment to an asset [2], [20]- [23].Prominent examples of such approaches are multipath-based methods that estimate MPCs associate them to virtual anchors representing the locations of the mirror images of an anchor on reflecting surfaces [24].The locations of virtual anchors are assumed to be known a priori [25] or estimated jointly with the position of agents using multipath-based SLAM (MP-SLAM) [21], [22], [26].Jointly estimating the positions of virtual anchors and agents allows MP-SLAM to provide high-accuracy position estimates, even in OLOS situations, or to localize the agent with only a single anchor [27].However, it requires specular, resolved MPCs, which are consistent with the virtual anchor model [28].Other methods exploit cooperation among individual agents [4], [27], [29], [30], or perform robust signal processing against multipath propagation and clutter measurements in general.The latter comprise heuristics [6], [31], machine learning-based approaches [19], [32]- [34] as well as Bayesian methods [35]- [37], and hybrids thereof [38]- [40].Heuristic methods, such as searching for the first amplitude to exceed a threshold value, are fast and easily implementable but suffer from low accuracy as well as a high probability of outage in low signal-to-noise-ratio (SNR) regions [6].In recent years, machine learning methods have grown increasingly popular.
Early approaches [19], [33] extract specific features from the radio channel applying model-agnostic supervised regression methods on these features.While these approaches potentially provide high accuracy estimates at low computational demand (after training), they suffer from their dependence on a large representative measurement database and can fail in scenarios that are not sufficiently represented by the training data.This is why recent algorithms facilitate deep learning and autoencoding based methods to directly operate on the received radio signal and reduce the dependence on training data [34], [41], [42].
Multipath-based localization [21], [22], [26], [36], [37], [43], [44], multiobject-tracking [45]- [47], and parametric channel tracking [48] are applications that pose common challenges, such as uncertainties beyond Gaussian noise, like missed detections and clutter, an uncertain origin of measurements, and unknown and time-varying number of objects to be localized and tracked.These challenges can be well addressed by Bayesian inference leveraging graphical models to perform joint detection and estimation.Since the measurement models of these applications are nonlinear, most methods typically rely on sampling techniques such as recursive Monte Carlo sampling or particle filtering, or use linearized Gaussian models [49], [50].Similarly, the probabilistic data association (PDA) algorithm [45], [51] represents a low-complexity Bayesian method for robust localization and tracking with extension to multiple-sensors PDA [52] and amplitude-information probabilistic data association (AIPDA) [44], [53].All these methods can be categorized as "two-step approaches", in the sense that they do not operate on the received sampled radio signal, but use extracted measurements provided by a preprocessing step, providing a high level of flexibility and a significant reduction of computational complexity.In contrast, "direct positioning approaches" such as [20], [54], [55] directly exploit the received sampled signal, which can lead to a better detectability of low-SNR features, yet, they are computationally very demanding.

B. Contributions
In this paper, we propose a particle-based sum-product algorithm (SPA) that sequentially estimates the position of a mobile agent by utilizing the position-related information contained in the LOS component as well as in MPCs 2 .The proposed algorithm jointly performs probabilistic data association and estimation of the mobile agent state [22], [46] together with all relevant model parameters, employing the SPA on a factor graph [56].Similar to other two-step approaches, it uses signal component delays and amplitudes estimated by a snapshot-based parametric channel estimation and detection algorithm (CEDA) as measurements.The proposed algorithm adapts in an online manner the time-varying component SNR [44] as well as the detection probability of the LOS [43], [57].To this end, we propose a novel detection probability model that allows for both an exhaustive representation of the detection space and a smooth estimate of the SNR.The algorithm exploits a novel non-uniform "false alarm (FA) model" 3 .Additionally, the model couples MPC measurements to the LOS measurement by a jointly inferred bias state.This enables the algorithm to utilize the positionrelated information contained in the MPCs without inferring specific map information, which can increase the accuracy and robustness of the agent's position estimate in challenging environments, characterized by strong multipath propagation and temporary OLOS situations.The proposed algorithm is able to operate without any prior information (no floorplan information or training data are needed).It is demonstrated to provide robust estimates for specular, resolved multipath as well as dense, non-resolvable multipath, while offering sub-second runtime 4 even in environments characterized by strong multipath propagation and, thus, a high number of measurements.The contributions of this paper are as follows.
joint distribution of delay and amplitude measurements instead of using heuristical models, (ii) sequentially inferring all parameters of the NLOS model together with the agent instead of using predetermined constants, (iii) improving the convergence behavior using a modified, "decoupled" SPA (see Sec. VI-A), (iv) demonstrating the performance of the proposed algorithm using simulated radio signals as well as real radio measurements obtained by (v) applying a CEDA, (vi) comparing to the MP-SLAM algorithm [22], [44] and (vii) providing the P-CRLB as a performance benchmark.

II. NOTATIONS AND DEFINITIONS
Column vectors and matrices are denoted by boldface lowercase and uppercase letters.Random variables (RVs) are displayed in san serif, upright font, e.g., x and x and their realizations in serif, italic font, e.g. , x and x; x denotes the true value of x.The same notation applies for stochastic processes x(t) and their realizations x(t).f (x) and p(x) denote, respectively, the probability density function (PDF) or probability mass function (PMF) of a continuous or discrete RV x. (•) T , (•) * , and (•) H denote matrix transpose, complex conjugation and Hermitian transpose, respectively.∥•∥ is the Euclidean norm.|•| represents the cardinality of a set.diag{x} denotes a diagonal matrix with entries in x.I [•] is an identity matrix of dimension given in the subscript.[X] n,n denotes the nth diagonal entry of X.Furthermore, 1 A (x) denotes the indicator function that is 1 A (x) = 1 if x ∈ A and 0 otherwise, for A being an arbitrary set and R + is the set of positive real numbers.We predefine the following PDFs with respect to x: The truncated Gaussian PDF is with mean µ, standard deviation σ, truncation threshold λ and Q(•) denoting the Q-function [61].Accordingly, the Gaussian PDF is f N (x; µ, σ) = f TN (x; µ, σ, -∞).The truncated Rician PDF is [62, Ch. 1.6.7] ) with non-centrality parameter u, scale parameter s and truncation threshold λ.I 0 (•) is the 0th-order modified first-kind Bessel function and Q 1 (•, •) denotes the Marcum Q-function [61].The truncated Rayleigh PDF is [62, Ch. 1.6.7] with scale parameter s and truncation threshold λ.This formula corresponds to the so-called Swirling I model [62].Finally, we define the uniform PDF f

III. RADIO SIGNAL MODEL
At each discrete time n, the mobile agent at position pn transmits a signal s(t) and each anchor j ∈ {1, ..., J} at anchor position p Ay ] T acts as a receiver.The complex baseband signal received at the jth anchor is modeled as The first and second term describe the LOS component and the sum of K(j) n specular MPCs with their corresponding complex amplitudes α(j) n,k ∈ C and delays τ (j) n,k ∈ R + , respectively.The delays are related to respective distances via τ (j) n,k = d(j) n,k /c with c being the speed of light.The third term represents an additive white Gaussian noise (AWGN) process w A ∥. We assume time synchronization between all anchors and the mobile agent 5 .However, our algorithm can be extended to an unsynchronized system along the lines of [2], [21], [64].
The signal r (j) 4) is uniformly sampled with sampling frequency f s at corresponding sampling interval T s = 1/f s and N s samples are collected, yielding a duration T = N s T s .By stacking the samples, we obtain the discrete time signal vector where is the stacked signal vector containing the samples of the transmit signal s(t).The measurement noise vector w n ∈ C Ns×1 is a zero-mean, circularly-symmetric complex Gaussian random vector with covariance matrix σ(j)2 I Ns and noise variance σ(j)2 = Ñ (j) 0 /T s .The MPCs arise from reflection or scattering by unknown objects, since we assume that no map information is available.

(j)
Dn (τ ) is the DPS.Using (7), the SNR of the LOS component is defined as6 SNR (j) n = .

A. Delay Power Spectrum (DPS) Model
We choose to model the DPS S (j) which is a double exponential function with Ω(j) n being the DPS power.The rise time γrn and fall time γ(j) fn are shape parameters.The distance difference ∆ b (τ ; •) is given by n where b(j) n is the NLOS bias, which denotes the difference between the LOS distance d fn γrn ] T collects the NLOS shape parameters for each time n and anchor j.Experimental evidence motivates this model: The DPS typically exhibits an exponentially decaying tail [5], [9] and a smooth onset [9], [68].In particular, when the LOS power is excluded, as is done in (6).Note that γrn is mainly determined by the signal bandwidth and onsetdensity of MPCs.For homogeneous deployment environments the on-set density is well modeled as being invariant.Therefore, γrn is assumed to be the same for all anchors.
For inference, we also define the normalized DPS SD (d, pn , n and the densemultipath-to-noise ratio (DNR) ω(j) / σ(j) , where the DNR denotes the square-root power ratio between the dense multipath component and AWGN.
The proposed algorithm utilizes the position information contained in S (j) Dn (τ ) to improve the position estimate without explicitly exploiting map information.

B. Parametric Channel Estimation
By applying a suitable snapshot-based channel estimation and detection algorithm (CEDA) [5], [69]- [71] to the observed discrete signal vector r n into individual, decorrelated components according to (5), reducing the number of dimensions (as M (j) n is usually much smaller than N s ).It thus can be said to compress the information contained in r See the supplementary material [72, Sec.V] for further details.The stacked vector ] T is used as noisy measurement by the proposed algorithm.

IV. SYSTEM MODEL
We consider a mobile agent to be moving along an unknown trajectory as depicted in Fig. 1.The current state of the agent is described by the state vector T , which is composed of the agent's position p n = [p xn p yn ] T and velocity v n = [v xn v yn ] T .We also introduce the following additional state variables, which represent all RV inferred along with x n : First, we define the augmented agent state xn = [x T n γ rn ] T , which collects all RVs that are common for all anchors.Second, we define the anchor state y fn ] T collecting all continuous RVs, which are modeled separately for each anchor.Third, there are two discrete RVs q (j) n and a (j) n , which denote the LOS probability and association variable, respectively, and are modeled separately for all anchors.For the sake of clarity, all RVs constituting the system model are summarized and described in Table I.
At each time n and for each anchor j the CEDA provides the currently observed measurement vector z (j) n , with fixed M (j) n , according to Sec.III-B.Before the measurements are observed, they are random and represented by the vector z un,m ] T .In line with Sec.III-B we define the nested random vectors z ] T .Also the number of measurements M (j) n is a RV.The vector containing all measurement numbers is defined as n,m either originates from the LOS or is due to an MPC.It is also possible that a measurement z (j) n,m did not originate from any physical component, but from FAs of the CEDA.The presented model only distinguishes between "LOS measurements" originating from the LOS and "NLOS measurements", i.e., measurements due to MPCs or FAs.

A. LOS Measurement Model
The LOS likelihood function (LHF) of an individual distance measurement z (j) dn,m is given by with mean d LOS (p n ) and variance σ 2 d (u The variance is determined based on the Fisher information given by σ 2 d (u , where β bw is the root mean squared bandwidth [1], [2] and u (j) n is the normalized amplitude at anchor j.The LOS LHF of the normalized amplitude measurement z (j) un,m is modeled as 7 [44], [48] with f TRice (•) being a truncated Rician PDF (2).γ is the detection threshold of the CEDA, which is a constant to be chosen.As for the distance LHF, the scale parameter is determined based on the Fisher information given as σ 2 u (u 1/(4N s ).Note that this expression reduces to 1/2 if the AWGN noise variance σ (j)2 is assumed to be known or N s to grow indefinitely (see [48] for a detailed derivation).Note that for (10) the Marcum-Q function in (2) represents the detection probability p D (u

B. NLOS Measurement Model
The NLOS LHF of an individual normalized amplitude measurement z un,m is given as un,m ; s u (z is the NLOS scale function.We used ζ Sn ] T for notational brevity.See the supplementary material [72, Sec.I] for details about the derivation of (11).The shape of (11) with respect to z where ) dd is the normalization constant ensuring integration to 1.The exponential term in (13) corresponds to the probability that at time n for anchor j a NLOS measurement at distance z (j) dn,m is generated.The shape of ( 13) with respect to z (j) dn,m for different values of γ is shown in Fig. 2b.Note that (13) approaches a uniform PDF when γ or ω (j) n approach zero.The presented NLOS measurement model is valid independently of the DPS model chosen in (8).However, ( 8) is a reasonable choice as it is physically motivated [9] and is of moderate computational complexity.

C. Data Association Model
At each time n and for each anchor j, the measurements, i.e., the components of z (j) n are subject to data association uncertainty.Thus, it is not known which measurement z (j) n,m originated from the LOS, or which one is due to an "NLOS measurement", i.e., measurements originating from MPCs or FAs.Based on the concept of PDA [51], we define the association variable a Assuming the number of NLOS measurements to follow a uniform distribution (so called "non-parametric model"), the joint PMF of a n can be shown to be [51] p(a (j)  n , M (j) where p E (u n ) is the probability that there is a LOS measurement for the current set of measurements defined in Sec.IV-D and M max is an irrelevant constant.Incorporating a (j) n into the model, we define the overall distance LHF as 7 The presented model describes the distribution of the amplitude estimates of a complex baseband signal in AWGN obtained using maximum likelihood estimation and generalized likelihood ratio test detection [53], [61], [73].
where we used ζ T n ] T for brevity.The shape of ( 16) is depicted in Fig. 3a.Further, the overall amplitude LHF is given by f (z (j)  un,m |z which is shown in Fig. 3b.Using the common assumption of the measurements to be independent for different values of m [46], the joint LHF for all measurements per anchor j and time n is En ) f (z

D. LOS Existence Probability Model
We model the LOS existence probability given in (15) as n ) is the probability that at time step n and anchor j the agent generates a radio signal component whose amplitude is high enough so that it leads to an LOS measurement.It is modeled by the counter probability of a Rician cumulative distribution function (CDF) given as by assuming that the proposed algorithm is applied after a generalized likelihood ratio test detector.q n is the probability of the event that the LOS is not obstructed, which is referred to as LOS probability in the following, and acts as a prior probability to the detection event.According to [43], [57], [74], we model q (j) n as discrete RV that takes its values from a finite set Q = {λ 1 , ... , λ Q }, where λ i ∈ (0, 1].The LOS probabilities for different sensors j are assumed to be independent.The proposed LOS existence probability model p E (u n ) correctly incorporates the detection process into the system model via p D (u (j) n ) excluding a detection of measurements with z (j) un,m below γ and it can cope with amplitude model mismatch by correcting the amplitude-related probability of detection with q (j) n .With respect to implementation (see Sec. VI-B1) this means that our model allows for smooth sequential inference of slow amplitude variations (e.g., due to path loss) via p D (u n ensures a complete representation of the probability space, covering rapid amplitude variations (e.g., due to OLOS).

E. State Transition model
We model the evolution of xn and y (j) n and q (j) n over time n as independent first-order Markov processes, which are defined by the joint state transition PDF being the elements of the transition matrix.

V. PROBLEM FORMULATION AND FACTOR GRAPH
In this section we formulate the sequential estimation problem of interest and present the joint posterior and the factor graph underlying the proposed algorithm.

A. Problem Statement
The problem considered is the sequential estimation of the agent state x n .This is done in a Bayesian sense by calculating the minimum mean-squared error (MMSE) [73] of the augmented agent state with In order to obtain (21), (22), and ( 23), the respective marginal posterior PDFs need to be calculated.Since direct marginalization of the joint posterior PDF is computationally infeasible [46], we perform message passing by means of the SPA rules on the factor graph that represents a factorized version of the joint posterior of our statistical model discussed in Sec.IV.

B. Joint Posterior and Factor Graph
For each n, let y n = [y We now assume that the measurements z are observed and thus fixed.Applying Bayes' rule as well as some commonly used independence assumptions [25], [46] the joint posterior for all states up to time n and all J anchors can be derived up to a constant factor as f ( x,a,y,q,M |z) ∝ f (z| x,a,y,q) f ( x,a,y,q) = f (z| x,a,y,q) f (a|y,q) f (x) p(q) f (y) where we introduced the state-transition functions n−1 ), and Ψ(q n−1 ).We also introduced the pseudo likelihood function ḡ(z neglecting the constant terms in (15).Note that M vanishes in (24) as it is fixed and thus constant, being implicitly defined by the measurements z.Furthermore note that unlike [22], [44], [46]- [48], [51]- [53] in our model the NLOS LHFs ( 11) and ( 13) are both functions of RVs and, thus, cannot be neglected.The joint posterior PDF in (24) is represented by the factor graph shown in Fig. 4.

A. Marginal Posterior and Sum-Product Algorithm (SPA)
The marginal posterior can be calculated efficiently by passing messages on the factor graph according to the SPA [56].For the proposed algorithm, we specify not to send messages backward in time.This makes the factor graph in Fig. 4 an acyclic graph.For acyclic graphs the SPA yields exact results for the marginal posteriors [56].At time n, the following calculations are performed for all J anchors.The prediction messages are given as n−1 ) fy (y n−1 ) pq (q y 0 Fig. 4. Factor graph representing the factorization of the joint posterior PDF in (24) as well as the respective messages according to the SPA (see Sec. VI-A).
The following short notations are used: ηn ≜ η( xn), ϕ where fx ( xn−1 ), fy (y n−1 ) and pq (q n−1 ) are messages of the previous time n − 1.The measurement update messages are given by ḡ(z (j) n ; p n ,y (j) n ,a (j) n ,q (j) n ) dy (j) n (29) Finally, we calculate the marginal posteriors as We additionally compare the performance of the above optimum SPA to that of a suboptimal message passing algorithm, which we refer to as "decoupled SPA".Inspired by [43], we replace (30) by χ (j) ( xn ) = η( xn ) neglecting the mutual dependency of the uncertainties of individual anchor states y (j) n .We demonstrate this modified algorithm to lead to improved numerical stability for a low number of particles.
Hence, the particle-based implementation discussed in section VI-B1 addresses the decoupled variant of the presented SPA.

B. Implementation Aspects
1) Particle-based Implementation: Since the integrals involved in the calculations of the messages and beliefs ( 26)-( 32) cannot be obtained analytically, we use a computationally efficient sequential particle-based message passing implementation that provides approximate computation.Our implementation uses a "stacked state" [75], comprising the augmented agent state as well as the anchor states of all anchors J = {1, ..., J}.
i) Prediction: The beliefs fx ( xn−1 ) and fy (y n−1 ) for all j ∈ J calculated at the previous time step n − 1, respectively, are represented by I particles and corresponding weights, i.e., {x n−1 } I i=1 and {y for all j ∈ J .Weighted particles {x n } I i=1 and {y i=1 for all j ∈ J , representing the messages η( xn ) and ϕ(y 26) and ( 27) are determined as follows: For each particle x[i] n−1 and y n−1 with i ∈ {1, . . ., I}, one particle x′[i] n and y ) and f (y n−1 ) for all j ∈ J .ii) Measurement Update: The non-normalized weights representing the messages ξ (j) ( xn ) and ν(y 29) and ( 31) are calculated by An approximation of the message β(q 32) is given as iii) Belief Calculation and State Estimation: The above approximate messages are further used for calculating the non-normalized weights corresponding to the beliefs fx ( xn ) and fy (y n ) for all j ∈ J as ŵx and ŵy After normalization, i.e., wx n and wy , an approximation of the MMSE state estimates xMMSE n and ŷ(j)MMSE n in ( 21), ( 22) and ( 23) is given as . To avoid particle degeneracy [49], a resampling step 8 is performed as a preparation for the next time step n + 1 leading to equally weighted particles {x n−1 = 1/I} I i=1 and {y n−1 = 1/I} I i=1 for all j ∈ J representing the beliefs fx ( xn ) and fy (y The resulting problem complexity scales only linearly in the number of particles I and in the number of measurements M (j) n .For computational efficiency of the particle-based implementation the LOS LHF of the normalized amplitude measurement (10) is approximated by a truncated Gaussian PDF, i.e.,

2) Initial State Distributions:
We assume the initial state distributions to factorize as fx ( x0 ) = fp (p 0 ) fv (v 0 ) fγr (γ r0 ) and fy (y ).We propose to initialize the NLOS shape parameters as fγr (γ r0 ) = 0 , Q) taking all values of Q with equal probability.We assume the velocity vector v 0 to be zero mean, Gaussian, with covariance matrix σ 2 v I 2 and σ v = 6 m/s, as we do not know in which direction we are moving.
The remainder of the states are initialized heuristically, by assuming an initial measurement vector z 0 containing M (j) 0 measurements to be available.The normalized amplitude PDFs are initialized as f (j) u0 (u u0,max , 0.05 z (j) u0,max , γ) where z (j) u0,max is the maximum normalized amplitude measurement in z where the proposal distribution f (p init ) is drawn uniformly on two-dimensional discs around each anchor j, which are bounded by the maximum possible distance d max and a sample is drawn from each of the J discs with equal probability.The DNR PDFs are initialized as fω (ω init is determined as described in the supplementary material [72, Sec.II]. 3) Normalization of the NLOS Distance Likelihood: As discussed in Sec.IV-B, the NLOS LHF in (13) must be normalized by n ) cannot be determined analytically and, being a function of p n and ζ (j) n , it needs to be calculated for each individual particle (see Sec. VI-B1).Thus, we need an efficient numerical approximation.For details see the supplementary material [72, Sec.III].

VII. RESULTS
We validate the proposed model and analyze the performance gain caused by the features of the proposed algorithm using both synthetic data obtained using numerical simulation and real radio measurements.The performance is compared with the P-CRLB and that of the AIPDA.For synthetic measurements with geometry related multipath 9 , we also compare to the MP-SLAM algorithm presented in [22], [44].

A. Common Analysis Setup
The following setup and parameters are commonly used for all analyses presented.
The PDF of the joint agent state xn is factorized as yn , where the noise vector ε , σ γr = 0.5 γMMSE rn−1 .While σ a is set according to the maximum agent acceleration [62], for the STV of all other parameters we use values relative to the root mean squared error (RMSE) estimate of the previous time step n − 1 as a heuristic.Note that this choice allows no tuning of the STV to be required for all experiments presented, even though the propagation environments are considerably different.We used 5 • 10 4 particles before the first resampling operation and 5000 particles for inference during the track.We set the detection threshold as low as γ = 1.77 (5 dB) for all simulations, which allows the algorithm to facilitate low-energy MPCs (this choice is further discussed in Sec.VII-B).The set of possible LOS probabilities is chosen as Q = {0.01,0.33, 0.66, 1}.The state transition matrix For all other tuples {i, k}, [Q] i,k = 0 in order to encourage high LOS probabilities [57].For the numerical approximation of Q 0 (p n , ζ  to Table .II.It shows the algorithm variants implemented and the corresponding features that are enabled for an algorithm (x) or not ( ).When "q (j) n tracking" is deactivated, we set q (j) n = 0.999 for all n, j.When we use "decoupled SPA", the suboptimal message passing scheme presented in Sec.VI-A is used.Not applying the "non-uniform f NL " means ( 12) is replaced by s 2 u ≜ 1/2, and for AL4 ′ and AL5 ′ we use 5 • 10 4 particles instead of 5000.Note that AL1 represents a multisensor variant of the conventional AIPDA.The MP-SLAM algorithm is implemented according to [22], [44] using the measurements z (j) m,n , i.e., distance and amplitude measurements, as an input.For consistency, the state transition PDFs and initial state distributions of the agent state and normalized amplitude state are set as described in Sec.VII-A and VI-B2.For convergence, we had to use 5 • 10 4 particles and an anchor driving noise of σ An = 0.02 m, other parameters are P s = 0.999, µ n,1 = 0.05.The mean number of false alarms was approximated as µ FA = N s e −γ 2 (see [44] for definitions).For stability we increased the delay measurement variances of all virtual anchors (not the physical anchors) by a factor of 2 with respect to the Fisher information-based value.
As a performance benchmark, we provide the CRLB on the position error variance considering all visible LOS measurements of a single time step n, which we refer to as the snapshot-based positioning CRLB (SP-CRLB) [7], [59], [76], [77].Furthermore, we provide the corresponding posterior Cramér-Rao lower bound (P-CRLB) [60] that additionally  considers the dynamic model of the agent state and the "P-CRLB-LOS", which is the P-CRLB assuming the LOS component to all anchors is always available and, thus, provides a lower bound for the proposed estimator.See the supplementary material [72, Sec.VI] for further details.

B. Analysis on Synthetic Measurements
For the synthetic setup, we investigate the scenario shown in Fig. 5.The agent moves along a trajectory, with two distinct direction changes, where the agent velocity is set to vary around a magnitude of 0.8 m/s.It is observed at N = 190 discrete time steps n ∈ {1, ... , N } at a constant observation rate of ∆T = 100 ms, resulting in a continuous observation time of 19 s.We simulate three anchors, A1-A3, which are placed in close vicinity to each other.The limited directional diversity of the anchors (corresponding to a poor geometric dilution of precision (GDOP) [78]), poses a challenging setup for delay measurement-based position estimation.Note that the environment setup shown in Fig. 5, i.e., walls and resulting obstructions, are only used in Sec.VII-B2.For all synthetic radio measurements involving the proposed CEDA (see [72,Sec.V]), we choose the transmitted complex baseband signal s(t) to be of root-raised-cosine shape with a roll-off factor of 0.6 and a duration of 2 ns (bandwidth of 500 MHz).The signal is critically sampled, i.e., T s = 1.25 ns, with a total number of N s = 161 samples, amounting to a maximum distance of d max = 60 m.
1) Synthetic Measurements with Stochastic Multipath: In this section we present results using synthetic measurements generated by simulating the MPCs as zero mean stochastic process.More specifically, we compare results obtained by simulating the radio signal according to (6) and applying the CEDA to results obtained using fully synthetic measurements, which are generated according to Sec.IV without involving the CEDA.For fully synthetic measurements the average number of NLOS measurements per time n and anchor j prior to the simulated detection process was approximated as N s .Detection further reduces the prior number of NLOS We start by validating the system model presented in  Sec.IV.For this experiment the relatively defined STV are set with respect to the true values instead of the RMSE values, given as fn , σ γr = 0.5 γrn .Fig. 6, 7, 8a and 8b show the results of the performed numerical simulations.Fig. 6 shows MMSE estimates of all state variables as a function of time t ′ and compares to the respective true values.The MMSE estimates are determined according to ( 21)-( 23) using both fully synthetic measurements and CEDA-based measurements.Fig. 7 compares distance-model-agnostic, bin-based estimates (BBEs) of scale parameter and relative measurement frequency with the presented model functions, i.e., with the NLOS scale function (12) and the NLOS distance LHF (13).Each of the functions is determined both ways, using the MMSE estimates of ζ n \ã 1 n ,n ∈ {180, ... , 200}}.For details about the BBEs see the supplementary material [72, Sec.IV].This analysis is complemented by Figs.8a and 8b which show the position RMSE e RMSE n in two ways.First, as a function of the discrete observation time n and, second, as the cumulative frequency of the RMSE evaluated over the whole time span.Fig. 6 demonstrates that using CEDA-based measurements the MMSE estimates of the parameters of the NLOS LHF (i.e., the MMSE estimates corresponding to ζ (1)  n ) are slightly biased, in particular the DNR estimate ω(1)MMSE n .This effect is a consequence of the asymptotic bandwidth assumption used in the derivation of the NLOS likelihood model (see [72,Sec. I]).However, as in Fig. 6 the model functions parameterized with the MMSE values accurately fit the BBEs, the MMSE estimate of the agent position pMMSE n in Fig. 6 remains unbiased and, thus, the positioning performance in Figs.8a and 8b using "CEDA-based measurements" is identical to the performance using "fully synthetic measurements" up to random deviations.
In addition, in Figs.8a and 8b we compare to fully synthetic measurements with (i) known initial state distributions, slightly lowering the RMSE around n = 0, and (ii) assuming the parameters of ζ (j)  n to be known constants, leading to a significant increase of performance at the end of the full OLOS situation as the bias information does not vanish over discrete time n.With CEDA-based measurements we also compare to results where (i) we calculate the relatively defined STV using the RMSE values of the respective last time step n − 1 according to Sec.VII-A and where (ii) we use a uniform delay intensity function f NL (z (j) dn,m ) = 1/d max showing no significant degradation of performance.The latter result suggests that for low values of γ, the information provided in (13) is insignificant (c.f.Fig. 2b).Therefore, in what follows, we keep the uniform delay intensity function leading to a considerable reduction of runtime since Q 0 (p n , ζ (j) n ) does not need to be calculated (see also Sec.VI-B3 and Sec.VII-D).Next, we investigate the influence of the individual features of our algorithm as described in Sec.VII-A and Table II.Figs.8c and 8d show the RMSE of this experiment as a function of t ′ as well as the cumulative frequency of the RMSE.The RMSE of the multi-sensor AIPDA (AL1) mostly attains the P-CRLB during LOS and partial OLOS situations.A reason for that is that the angle, which the remaining anchors A1 and A3 span with respect to the agent is sufficiently large to provide a reasonable position estimate.However, AL1 shows a slightly increased RMSE around n = 80 due to the agent direction change and significantly deviates from the very beginning of the full OLOS situation, losing the track in every single realization.Comparing the curves of AL2-AL5, one can conclude that every single algorithm feature investigated lowers the RMSE significantly when activated., (c) amplitude domain (û The RMSE of the proposed algorithm AL5 constantly attains the P-CRLB, which indicates no lost track, even falling below the P-CRLB in full OLOS situations.This is possible as it leverages the additional position information contained inside the MPCs via the non-uniform NLOS LHF, which is not considered by the P-CRLB model.In contrast, AL2 loses a large percentage of tracks after the full OLOS situation, because NLOS measurements significantly contribute to the LOS based position hypotheses due to the insufficient representation of the existence probability by the amplitude state particles (see Sec. IV-D).While AL3 constantly attains the P-CRLB during the LOS situation as well the partial OLOS situation, it loses the track for every realization in full OLOS.After a short amount of time in which AL3 can maintain the agent position through the agent state transition model and the decreasing LOS probability, it identifies MPCs as the LOS component due to their coherent appearance and large amplitude, which is not covered by the uniform NLOS model, and loses the track.AL4 shows a seemingly random performance degradation, which is due to the insufficient representation of the high dimensional joint state by the particle filter and some resulting lost tracks, which AL5 overcomes by decoupling the anchor states (see Sec. VI).However, the discrepancy between AL4 and AL5 can be dissolved by using a sufficiently high number of particles (see AL4 ′ and AL5 ′ ), at the cost of significantly increasing the runtime (see Sec. VII-D).
2) Synthetic Measurements with Geometry-related Multipath: In this section, we discuss results using synthetic measurements based on the simple floorplan shown in Fig. 5.The measurements are obtained by simulating a radio signal according to (5), consisting of the LOS component and specular MPCs, and using the proposed CEDA.The MPC delays are calculated out of the floorplan (i.e.W1-W5) using the mirror images (virtual anchors) up to the third order [24].The SNR of the LOS component as well as the MPCs [48] are set to 20 dB at a distance of 1 m and are assumed to follow free-space path loss.The SNR of the individual MPCs are additionally attenuated by 3 dB after each reflection (e.g., 6 dB for a second-order reflection).As depicted in Fig. 5, for this experiment the anchors are obstructed by an obstacle (W5), which leads to partial and full OLOS situations in the center of the investigated trajectory.Figs. 9, 10a, and 10b show results of the performed numerical simulation.Fig. 9 provides a graphical representation of the measurement space, showing a single measurement realization z together with the corresponding MMSE estimates of the proposed algorithm (AL5).The MMSE estimates are determined according to ( 21)- (23).In particular, Fig. 9a shows that (i) the MMSE estimate of the LOS delay d(j)MMSE ) remains stable over the whole OLOS situation and that (ii) the maximum of the NLOS LHF follows the first MPC available.We determine the shape of the NLOS LHF using the respective MMSE estimates of all RVs on which (11) depends.Fig. 9c shows that the DNR estimate ω(j)MMSE n accurately represents the dynamic behavior of the multipath energy, deceasing rapidly when the strongest, first MPC is covered, while the SNR estimate û(j)MMSE n remains stable.For visualization, Fig. 9b shows the NLOS scale function at time n = 133 parametrized with the respective MMSE estimates of all NLOS function parameters.Fig. 9d shows the LOS existence probability q (1) n well representing the OLOS situation.Figs.10a and 10b show the RMSE as a function of the discrete observation time n as well as the cumulative frequency of the RMSE.Again, we investigate the influence of the individual features of our algorithm according to Sec.VII-A and Table II.Comparing the presented curves, we again observe AL5 to significantly outperform AL1-AL5, with the qualitative performance differences being almost identical to those of Sec.VII-B1.The only significant dissimilarity is the seemingly smaller deviation between AL4 and AL5.This is because AL4 does not lose  any tracks during initialization, as the average energy and distance to the LOS component of the measurements of the first time step n = 0 are significantly lower in this scenario, leading to a better coverage of the state space by the particle filter.Thus, we only observe a slightly more unstable local behavior of AL4.The MP-SLAM algorithm (AL6) achieves a significantly reduced RMSE during the first part of the OLOS situation, due to geometric information provided by the specular MPCs, outperforming the proposed method (AL5).However, the investigated scenario is geometrically ambiguous as there is little directional change in the agent movement [28].Also there are many low-SNR components, which disappear and reappear, due to the obstacle (W5).This is why AL6 follows ambiguous paths for many realizations (i.e., it loses the track), leading to a significantly reduced performance after the full OLOS situation.We additionally added AL6 * , which represents the numerical results after removing 20.6% (103 realizations) of diverged tracks.This result demonstrates the dramatically increased accuracy that can be obtained using MP-SLAM.

C. Performance for Real Radio Measurements
For further validation of the proposed algorithm, we use real radio measurements collected in a laboratory hall of NXP Semiconductors, Gratkorn, Austria.The hall, shown in Fig. 11a, features a wide, open space and includes a demonstration car (Lancia Thema 2011), furniture, and metallic surfaces, thereby representing a typical multipath-prone industrial environment.An agent is assumed to move along a pseudo-random trajectory (selected out of a grid of agent positions), obtained in a static measurement setup.We selected N = 195 measurements, assuming an observation rate of ∆T = 170 ms.The agent velocity is set to vary around a magnitude of 0.35 m/s.This leads to a corresponding continu-ous observation time 33.15 s.At each selected position, a radio signal was transmitted from the assumed agent position, which was received by 4 anchors.The agent was represented by a polystyrene build, while the anchor antennas were mounted on the demonstration car.The agent as well as the anchors were equipped with a dipole antenna with an approximately uniform radiation pattern in the azimuth plane and zeros in the floor and ceiling directions.The radio signal was recorded by an M-sequence correlative channel sounder with frequency range 3 − 10 GHz.Within the measured band, the actual signal band was selected by a filter with root-raised-cosine impulse response s(t), with a roll-off factor of 0.6, a twosided 3-dB bandwidth of B = 499.2MHz and a center frequency of 7.9872 GHz (corresponding to channel 9 of IEEE 802.15.4a), and critically sampled with T s = 1/(1.6B).We used N s = 161 samples, amounting to a maximum distance of d max = 60 m for the CEDA.We created two full OLOS situations at n ∈ [80, 92] and n ∈ [159, 170] using an obstacle consisting of a metal plate covered with attenuators as shown in Fig. 11b.A floor plan showing the track, the environment (i.e, the car, other reflecting objects and walls), the antenna positions, and OLOS conditions with respect to all antennas is shown in Fig. 11c.The metal surface of the car strongly reflected the radio signal, leading to a radiation pattern of 270 • for A1 and A2 and 180 • for A3 and A4.Thus, during large parts of the trajectory the LOS of 2 or 3 out of 4 anchors is not available.Moreover, the pulse reflected by the car surface strongly interferes with the LOS pulse, leading to significant fluctuations of the amplitudes.In addition, this leads to the channel estimator being prone to produce a high SNR component just after the LOS component.As this violates our signal model, we processed the CEDA measurements attenuating all components, where z  and A2) are visible at the track starting point, the position estimate obtained by trilateration is ambiguous.In the scenario presented, the relative antenna position with respect to the car can be assumed to be known.Thus, for this experiment, we used the antenna pattern as prior information for initialization of the position state.For the numerical evaluation presented, we added AWGN to the real radio signal obtained.We set ∥r (j) raw ∥ 2 /σ (j)2 = 20 dB, where ∥r (j) raw ∥ 2 is the average energy of the real measured signal per anchor j.Figs.10c and 10d show the RMSE as a function of the discrete observation time n as well as the cumulative frequency of the RMSE.Again, we analyze the influence of the individual features of our algorithm according to Sec.VII-A and Table II and observe AL5 to significantly outperform the other algorithm variants.Different to Sec.VII-B2 all presented algorithms fail to reach the P-CRLB over parts of the track.The exact consistency in progression of the RMSE curves suggests unmodeled effects (e.g.diffraction at the vehicle body) as well as inaccuracies in the reference as a probable reason.

D. Runtime
Table III shows the average runtime of the proposed algorithm (A5) and compares it to the runtime of the multi-sensor AIPDA (AL1) and that of the MP-SLAM algorithm (AL6).All runtimes are estimated using Matlab implementations executed on an AMD Ryzen Threadripper 1900X 8-Core Processor with up to 4 GHz for all scenarios investigated.We also show the average number of measurements (over all anchors and time steps) M mean , the number of anchors J and the number of particles, which determine the algorithm complexity per time step.The runtime of our algorithm (AL5) is of the same order of magnitude than that of the multi-sensor AIPDA (AL1), which is in the range of tens of milliseconds for all scenarios investigated.In contrast, the runtime of the MP-SLAM algorithm (AL6) is significantly higher, since it requires joint data association between all map features [22] and a higher number of particles for numerical stability.

VIII. CONCLUSION
We have presented a particle-based sum-product algorithm (SPA) that sequentially estimates the position of a mobile agent using range and amplitude measurements provided by a snapshot-based channel estimation and detection algorithm (CEDA).We introduced a novel non-line-of-sight (NLOS) model that is adapted to the delay power spectrum (DPS) of the multipath radio channel.We analyzed the performance of the proposed algorithm using both numerically simulated and real measurements in different channel conditions and showed that the additional information provided by the NLOS model can support the estimation of the agent position.Our algorithm significantly outperformed the conventional AIPDA filter and consistently attained the P-CRLB in partial OLOS situations (i.e., no lost tracks).While multipath-based SLAM (MP-SLAM) can naturally outperform our algorithm in channels showing resolved, specular MPCs, we demonstrate the proposed algorithm to offer a significantly smaller number of lost tracks at reduced execution time in a geometrically ambiguous scenario.A possible direction for future research includes extending the model to multiple biases with respect to several MPCs using joint probabilistic data association and dynamic MPC initialization [46], [48] or to several MPC clusters by using data association with extended objects [47].
While the effect of estimating σ becomes negligible 2 for a large number of samples N s , which is true for wideband ranging applications in general [9], estimating τ leads to a small but visible bias in the scale parameter [10].Following the steps of [10], this bias, which is usually ignored in amplitude models [5], [11]- [13], can be represented by replacing the truncated Rayleigh distribution of (5) with a truncated Rician distribution [1, Eq. ( 2)] with non-centrality parameter of √ 0.5.However, different from [10] the expected power of the stochastic process in (1) is not a constant with respect to the elements of r N and, thus, neither is the corresponding scale parameter (6) with respect to τ .To take this behavior into account, we model the expected amplitude of the stochastic process to be constant in the local environment, i.e., we modify the non-centrality parameter of √ 0.5 by the ratio of the current scale parameter to the scale parameter for AWGN, given as s N (τ, σ)/ √ 0.5.We get f (u N ; H 1 , τ, σ) = f TRice (u N ; s N (τ, σ), s N (τ, σ), γ) .
We validate the above model by a numerical simulation study.The results are provided in Fig. 1.We show (6) and (8) together with estimates obtained by applying the proposed CEDA (see Sec. V) to simulated radio signals.The simulated radio signals are obtained using numerical simulation according to (1), i.e., they consist of a stochastic process only.The DPS parameters are set to constant values, chosen in line with [1, Sec.VII-B1] for simulation.We used γ = 0 dB in order accept all CEDA estimates.For simplicity, the distance to the LOS component was assumed to be equal to zero.For this experiment, we use two versions of the CEDA that solve the optimization problem in Sec.V, (23) in two ways: First, we use the normal variant as suggested in Sec.V involving continuous unconstrained optimization [14] (BBE Rayleigh opt, BBE Rice opt.) and, second, we use grid-based optimization only, where the grid values correspond to the sampling grid of the radio signal (BBE Rayleigh grid).Out of the estimates provided by the CEDA we then estimate the scale parameter by grouping the estimates into delay/distance bins and applying, respectively, the maximum likelihood estimator for a truncated Rayleigh distribution (BBE Rayleigh opt, BBE Rayleigh grid), or truncated Rice distribution (BBE Rice opt.); 2 For unknown σ, the statistic of two times the squared normalized amplitude 2 u 2 N is described by a Fisher distribution [8,Ch. 15.10.3] with numerator degrees of freedom equal to 2 and denominator degrees of freedom equal to 2 Ns.For large Ns the statistic of 2 u 2 N can be well approximated by a χ 2 distribution [7, Ch. 2.2] and, therefore, the statistic of u N by the Rayleigh distribution described in the main text.see Sec.IV for details abut the bin-based estimation process.
We simulated 600 signals amounting to approximately 500 samples per estimation bin (empirically determined value) at the signal parameters and detection threshold configured.
Analyzing Fig. 1 one can observe that the simplified, asymptotic model ( 8) significantly underestimates the spread parameter as it neglects correlations in C S (τ ) occurring due to finite signal bandwidth.The correlation-aware model in (5) accurately represents the spread parameter for grid-based estimates of τ (BBE Rayleigh grid).However, continuous optimization of τ (BBE Rayleigh opt) leads to an offset, which can be considered using the Rician model ( 9) instead (BBE Rice opt.).However, as we assumed the noise level to be constant in the local environment, (9) cannot represent the influence of finite signal bandwidth with respect to estimation of τ .This effect is visible at lag 0, due to the rapid change of variance in the adjacent region.
The above analysis showed that the approximate model consisting of ( 5) and ( 8) is insufficient as a model for generating measurements.However, in [1, Sec.VII-B1] we demonstrate it to suffice as a model for inference: We obtain no loss in performance3 of the proposed algorithm when comparing results with data generated according to (8) to results using the stochastic radio signal model [1, Eq. ( 6)] and the proposed CEDA (see Sec. V) for generating measurements 4 .In contrast, the runtime of the overall algorithm using the approximate model is orders of magnitude lower than using the full model consisting of ( 9) and ( 6), especially since the latter requires numerical approximation of the convolutions in (7).See [1, Sec.VII-B1] for further details.

II. DENSE-MULTIPATH-TO-NOISE RATIO (DNR) INITIALIZATION
First we determine the ML estimator for a set of i.i.d.samples X = {x 1 , ..., x |X | } following a truncated Rayleigh distribution, i.e. we solve arg max s x∈X f TRayl (x; s, λ).This the additional information leveraged using the proposed NLOS model.However, in contrast to mapping approaches [24]- [26], which can facilitate multipath information via estimated map features (virtual anchors [27], [28]), our model just allows to mitigate the NLOS bias between MPC-related distance measurements and the LOS component distance.Thus, a strict lower bound can be obtained by assuming the LOS component to be available at all times n, i.e., setting 1 V (j) n ≜ 1 in (26).We refer to the corresponding P-CRLB as P-CRLB-LOS in the main text.

n
(t) with double-sided power spectral density Ñ (j) 0 /2.The LOS distance is geometrically related to the agent position via d(j) n,0 ≜ d (j) LOS ( pn ) with d (j)

n
, one obtains, at each time n and anchor j, a number of M (j) n measurements denoted by z (j) n,m with m ∈ M (j) n = {1, ... , M (j) n }.Each z (j) n,m = [z (j) dn,m z (j) un,m ] T contains a distance measurement z (j) dn,m ∈ [0, d max ], with maximum distance d max = c T , and a normalized amplitude measurement z (j) un,m .The CEDA decomposes the discrete signal vector r (j) is shown in Fig.2a.The NLOS LHF of the distance measurement z (j) dn,m is given by

Fig. 3 .
Fig. 3. Graphical representation of the (a) the overall distance LHF f (z (j) dn,m |ζ (j) En ) and (b) overall amplitude LHF f (z (j) un,m |z (j) dn,m ,ζ (j) En ), all at fixed ζ (j) En in line with Fig. 2. In (b) we also fix z (j) dn,m to z d1 = 0 m or z d2 = 8 m. with f ( xn | xn−1 ) and f (y (j) n |y (j) n−1 ) being the respective state transition PDFs.For the discrete RV q (j) n the first-order Markov process model results in a conventional Markov chain, with [Q (j) ] i,k = p(q

n
) as discussed in Sec.VI-B3, we used K T = 30.The results are shown in terms of the RMSE of the estimated agent position e RMSE n = E{∥ pMMSE n − p n ∥ 2 }, evaluated using a numerical simulation with 500 realizations.For each of the scenarios investigated, we consistently analyze the influence of the individual features of our algorithm according

Fig. 6 .
Fig.6.MMSE estimates (determined using fully synthetic measurements and CEDA-based measurements) and true values of all state variables v.s.discrete time n for the experiment described in Sec.VII-B1.We show the mean MMSE estimate and the corresponding range from minimum to maximum value.Anchor state variables are only shown for anchor A1.Different shades of gray represent different numbers of anchors in OLOS according to Fig.5.

Fig. 7 .
Fig. 7. Bin-based estimates (BBEs) of (a) the Rayleigh scale parameter of the amplitude measurements and (b) the relative frequency of the distance measurements compared to (a) the NLOS scale function and (b) the NLOS distance LHF.All values are shown as a function of the difference of the distance measurement and the corresponding LOS component distance, given as z (1) dn,m − d (1) LOS ( p200 ) for anchor A1.

Fig. 8 .
Fig. 8. Performance of different algorithm variants in terms of the RMSE of the estimated agent position (a), (c) as a function of the discrete observation time n, and (b), (d) as the cumulative frequency in inverse logarithmic scale, determined from numerical simulation of stochastic multipath according to Sec.VII-B1.Different shades of gray represent different numbers of anchors in OLOS according to Fig. 5.

( 1 )
200 of the last time step, given as ζ(1)MMSE 200 and using the respective true values used for simulation ζ(1) .The BBEs are determined using all NLOS measurements (the LOS measurements are removed) of the last 20 time steps, given as{z

2 uFig. 9 .
Fig. 9.A single measurement realization and the respective MMSE estimates using the proposed algorithm (AL5) for geometry-related synthetic measurements (see Sec. VII-B2).We show MMSE estimates in (a) distance domain d(1)MMSE n,0 MMSE2 n ), (d) LOS probability domain q(1)MMSE n , and (b) amplitude as a function of distance domain for n ∈ {130, ... , 136}.The measurements in (b) correspond to the upper and in (c) to the left axis.Note that in (a) and (c) only measurements for every third time n are shown (i,e., n ∈ {1, 4, 7, ...}).

Fig. 10 .
Fig. 10.Performance of all algorithm variants of Table II in terms of the RMSE of the estimated agent position (a), (c) as a function of the discrete observation time n, and (b), (d) as the cumulative frequency in inverse logarithmic scale, determined from numerical simulation of specular MPCs according to Sec.VII-B2 in (a), (b) and using real radio measurements according to Sec.VII-C in (c), (d).Different shades of gray represent different numbers of anchors in OLOS according to Fig. 5.

Fig. 11 .
Fig. 11.Measurement setup for real radio-signal experiments described in Sec.VII-C.We show pictures of (a) the overall scenario and (b) the OLOS setup used, as well as (c) the abstracted floorplan and trajectory.

Fig. 1 .
Fig. 1.Bin-based estimates (BBEs) of the scale parameter of the respective truncated Rayleigh or Rician distribution v.s. the NLOS scale function [1, Eq. (12)] as a function of the difference of a distance measurement to the LOS component distance.The non-parametric estimates are determined as shown in Sec.IV.

TABLE I SUMMARY
AND DESCRIPTION OF ALL UNOBSERVED RVS OF THE SYSTEM MODEL.
, where the agent motion, i.e. the state transition PDF f (x n |x n−1 ) of the agent state x n , is described by a linear, constant velocity and stochastic acceleration model[62, p. 273], given as x n = A x n−1 +B w n , with the acceleration process w n being i.i.d.across n, zero mean, and Gaussian with covariance matrix σ 2 a I 2 , σ a is the acceleration standard deviation, and A ∈ R 4x4 and B ∈ R 4x2 are defined according to[62, p. 273], with sampling period ∆T .The state transition of the rise distance γ rn , i.e., the state transition PDF f (γ rn |γ rn−1 ), is γ rn = γ rn−1 + ε γrn , where the noise ε γrn is i.i.d.across n, zero mean, Gaussian, with variance σ 2 γr .Similarly, the state transition model of the joint anchor state y n , i.e. the state transition PDF f (y n |y n−1 ), is chosen as y