Empirical Fading Model and Bayesian Calibration for Multipath-Enhanced Device-Free Localization

MDFL systems determine presence and location of objects and users not necessarily equipped with localization devices. For localization, multipath-enhanced device-free localization (MDFL) systems exploit user-induced changes in the power of all received signal components, including both line-of-sight and multipath components (MPCs). In this work, we therefore provide a statistical fading model that describes user-induced changes in received power specifically for MPCs. The model is derived and validated empirically using an extensive set of wideband and ultra-wideband measurement data. Since the localization performance of MDFL systems strongly depends on the information about the propagation paths within the wireless network, we further propose a Bayesian calibration approach that estimates the location of the reflection points of MPCs caused by single-bounce reflections. For MPCs caused by single-bounce reflections, the solution space of possible locations of reflection points is constrained to the delay ellipse, which allows the formulation of a computationally efficient one-dimensional estimation problem. Eventually, the problem is solved by sequential Bayesian estimation. The applicability of the proposed approach is demonstrated and evaluated using measurement data. Independent of the underlying measurement system, the Bayesian calibration approach is shown to robustly estimate the locations of the reflection points in different environments. Finally, the localization results of MDFL for an indoor scenario confirm the applicability of the Bayesian calibration approach.


I. INTRODUCTION
Our homes [1], our cities [2], and our industries [3] are becoming smart.Primarily, this trend is driven by the increasing, nigh ubiquitous connectivity.The number of connected devices is steadily growing, and with it the demand for location-based services [4].This demand in location-awareness can be served for users carrying a localization device, e.g., by active radio frequency (RF)-based localization systems [5].However, users are not always equipped with such a device, but the demand in location-awareness still exists.Possible use cases range from monitoring applications for health and elderly care, to home security and intruder detection, to audio applications for HiFi [6].Therefore, alternative passive localization systems are required enabling the localization of passive users that are not equipped with active localization devices.As wireless connectivity increases, RF-based passive localization or sensing The authors are with the Institute of Communications and Navigation, German Aerospace Center (DLR), 82234 Wessling, Germany (e-mail: martin.schmidhammer@dlr.de;christian.gentner@dlr.de;m.walter@dlr.de;stephan.sand@dlr.de;benjamin.siebler@dlr.de;uwe.fiebig@dlr.de).
systems are becoming a viable option.Here, the presence and location of users can be estimated by means of RF sensor networks, exploiting the physical impact of the user on RFsignals [7].We basically distinguish between radar systems using the properties of signals directly scattered off the user, e.g., [8], and device-free localization (DFL) systems using power changes of received signals within a wireless network caused by user-induced diffraction and shadowing effects, e.g., [9].For inferring presence and location of users, current narrowband DFL systems typically measure the received signal strength (RSS) between the network nodes in line-ofsight (LoS) [7], [9]- [14].The location of the user can then be either estimated by computing propagation field images, which is also known as radio tomographic imaging [9], or by modeling the changes in the received power as a function of the user location.The location estimation problem of the latter is then typically solved using Bayesian filtering, as in [14].
In the past, RSS-based DFL systems have attracted considerable attention since narrowband measurements are obtained from many commercially available wireless devices [11].However, recent advances in wireless communications show the trend towards wideband and ultra-wideband (UWB) technologies, which are already on the market.For signals of higher bandwidth capacities, individual multipath components (MPCs) caused by scattering and reflection from the surrounding environment can be resolved and the corresponding signal parameters can be estimated.The work in [15] provides an example for DFL using UWB devices.Here, the wide signal bandwidth is used to isolate the LoS signal for mitigating distortions due to multipath propagation.In [16], we have shown that users also induce fading to the received power of MPCs.Based on these findings, we proposed a novel multipath-enhanced device-free localization (MDFL) approach in [17], since the propagation paths of MPCs inherently differ from those in LoS.We demonstrated that the localization performance improves without additional infrastructure by exploiting MPCs as complementary DFL network links.In our previous work, we assume prior information about the surrounding environment, e.g., provided by a floor plan, for determining the propagation paths during an initial calibration.In this work, we address the problem of determining the propagation paths, even when no prior information about the surrounding environment is available or has altered after calibration.We thus describe a novel Bayesian calibration approach, that estimates the reflection points (RPs) of MPC caused by single-bounce reflections.
For estimating the RPs, the proposed Bayesian calibration approach relates measured changes in the received power of an MPC to the current location of the user in the network.Note that this relation is similarly required for model-based DFL and MDFL approaches.That means, we can apply similar fading models that express the power changes as a function of the user location.For the RSS, i.e., the received power of the direct link between transmitting and receiving node, there is a large number of fading models.We can thereby distinguish between theoretical propagation models based on diffraction theory, e.g., [11], [18], [19], and empirical propagation models, that statistically model perturbations of the received power, e.g., [10], [13], [20].Regarding the received power of MPCs, only few models exist.In [16], we have adapted the diffraction based model of [18] for multipath propagation.Requiring a computational more efficient model, in [17] and [21], we have approximated the fading of MPCs by an exponential model.
In this work, we analyze, adapt, and evaluate the exponential model empirically using an extensive set of wideband and UWB measurement data.
The main contributions of this work are • the empirical derivation and validation of a statistical fading model for multipath components, • the novel Bayesian calibration approach for MDFL that estimates reflection points without a priori information about the surrounding environment, • the demonstration and evaluation of the proposed Bayesian calibration approach using measurement data.
The remainder of this paper is organized as follows.Required preliminaries including network and propagation model are given in Section II.The empirical fading model is derived and evaluated in Section III.The novel Bayesian calibration approach is described in Section IV and evaluated using measurement data in Section V. Finally, Section VI concludes the paper.

A. Network and Propagation Model
For the MDFL system, we consider a network of N Tx transmitting and N Rx receiving nodes.Transmitting and receiving nodes can be collocated or individually placed at known locations r Txi , i ∈ {1, . . ., N Tx }, and r Rxj , j ∈ {1, . . ., N Rx }, respectively.The index set P determines the link configuration of the network, with link (i, j) ∈ P being composed of the i-th transmitting and the j-th receiving node and indexed by l ∈ {1, . . ., |P|}.Further, we model the received signal of link l as a superposition of scaled and delayed replica of a known transmit signal s l (t) [22].Due to reflections and scattering of the surrounding environment, these comprise the LoS component and a finite number of N l MPCs.Thus, we express the received signal as reflecting surface Fig. 1: Propagation path for an MPC due to single-bounce reflection at a reflecting surface and corresponding 2D representation of the delay ellipsoid.
with α l,n (t) as time-variant, complex amplitude, τ l,n as static propagation delay of the n-th MPC, where n = 0 refers to the LoS component, and n l (t) as white circular symmetric normal distributed noise with variance σ 2 y l .The product of the propagation delay τ l,n and the speed of light c defines the geometric length of the propagation path as where for the LoS the length of the propagation path equals the distance between the transmitting and receiving node

B. Virtual Nodes and Reflection Points
For localization, the MDFL system requires information about the geometric propagation paths between the individual network nodes.While for LoS components the propagation paths are inherently defined by the locations of the respective transmitting and receiving nodes, the propagation paths of MPCs still have to be determined.In [17], we have used virtual nodes, i.e., virtual transmitters (VTs) and virtual receivers (VRs), for modeling the propagation paths.Thereby, the virtual nodes are constructed by successive mirroring of the physical nodes on reflecting surfaces.Corresponding VTs and VRs form equidistant pairs of nodes, whose distances equal the length of the propagation delay [17], [23], [24].The paths between these corresponding pairs of nodes intersect at the physical RPs on the reflecting surfaces, as exemplary shown in Fig. 1.Thus, similar to optical ray-tracing, we can geometrically reconstruct the physical propagation paths [23].

C. Geometric Properties for Single-Bounce Reflections
An important special case represents MPCs caused by firstorder reflections or scattering, i.e., single-bounce reflections (SBRs).Compared to MPCs caused by higher order reflections, the received power of an MPC caused by SBR is expected to be less attenuated [25].The comparatively high received power improves the quality of the MPC parameter estimates such as propagation delay and amplitude, cf.(1).
Besides the advantages due to the received power, MPCs caused by SBR comprise also specific geometric properties of the propagation paths, that we will outline in the following.
For SBRs, the set of virtual nodes reduces to one VT and one VR, so that the propagation path of the corresponding MPC can be described by a single RP.As shown in Fig. 1, the respective vector between the physical node, i.e., transmitter or receiver, and the RP is aligned with the vector between the physical and the corresponding virtual node.That means, for the transmitting and receiving node we can equally define the unit vectors and where r VT l,n , r VR l,n , and r RP l,n denote the locations of VT, VR, and RP, respectively, for the n-th MPC of network link l.
As mentioned in Sec.II-B, the distances between the physical transmitting or receiving node and the corresponding virtual node equal the length of the propagation path as By inserting ( 6) into ( 4) and ( 5), we can resolve for the locations of VT and VR as and Besides these unique expressions for the locations of VT and VR, the propagation delay of an MPC caused by SBR allows to constrain the possible locations of the corresponding RP.Specifically, this means that the possible locations of an RP coincide with the locations on the surface of an ellipsoid, each of which equally causes the propagation delay τ l,n [26], [27].For reflections from surfaces orthogonal to the axis between the transmitting and receiving nodes, e.g., vertical walls of buildings, the possible locations are reduced to the circumference of an ellipse, as shown in Fig. 1.Since we expect the majority of MPCs to originate from reflections off surrounding walls, we adopt this two-dimensional (2D) representation of an ellipse.We therefore also assume that all transmitting and receiving nodes are located at the same height, i.e., in the same plane.Note that all vectors and vector operations in the remainder of this paper are considered in 2D.In general, the shape of an ellipse is determined by its principal axes.These include the semi-major and the semiminor axis, which are defined for the n-th MPC of network link l as and with path lengths d l,n and d l as given in ( 2) and (3), see Fig. 1.The center of the ellipse is given by In a local coordinate system with the center of the ellipse r C l at the origin and the semi-major and semi-minor axes along x l -axis and y l -axis, see Fig. 1, the ellipse can be expressed analytically by and, equivalently, by the parametric representation with the ellipse parameter θ ∈ [0, 2π].Note that the local ellipse-centric coordinate system is defined individually for each network link l, but applies to all MPCs and corresponding delay ellipses of that link.Unlike a circle, there is no simple relation between the ellipse parameter with different semiminor and semi-major axes and the arc length [28].In order to describe the distribution of possible locations of an RP along the elliptic arc, we need to define the arc length as a function of the ellipse parameter [28], as with the numerical eccentricity Therewith, we can further calculate the circumference of the ellipse as With E( l,n ) = π/2 0 1 − 2 l,n cos 2 ϑ dϑ as the complete elliptic integral of the second kind, (16) Finally, we want to express the local ellipse-centric coordinates x l in the coordinate system used for the MDFL network.Thus, the local coordinates are first rotated by the rotation matrix R l and translated by the center of the ellipse r C l as Here, we consider the rotation of the local coordinates in the x y-plane with respect to the x-axis about the origin of the coordinate system used for MDFL.Therefore, the rotation matrix is expressed as with rotation angle α l .The angle determines the rotation between the x-axis and the x l -axis of the local coordinate system.With the unit vectors e x in x-direction and e x l = r Rxj − r Txi / r Rxj − r Txi in x l -direction the rotation angle can be calculated by where the direction of the rotation is determined by the sign of the determinant of the unit vectors, i.e., sgn l = sign (det [e x , e x l ]).

III. EMPIRICAL FADING MODEL FOR MPCS
MDFL systems aim to estimate the user location based on user-induced changes in the received signal power.Therefore, the measurement model that relates the measured changes in the received signal power to the user location is essential for any MDFL system.The performance of the localization system depends on the robustness, complexity, accuracy, and validity of the underlying model.In this section, we therefore derive and validate an empirical model describing the userinduced variations in the received power of MPCs, based on measurement data.

A. Data Collection
For investigating the impact of users on the power of MPCs, we have collected an extensive set of wideband and UWB measurement data in different multipath environments, both indoors and outdoors.Accounting for a variety of different influencing parameters, we consider data from three different measurement campaigns.Thereby, the sets contain data from different configurations, different environments, and different types of users.Table I provides a detailed overview of the settings and parameters used during the individual measurement campaigns.The corresponding network configurations and geometric setups are illustrated in Fig. 2. In the following, we briefly describe the individual measurement campaigns.
1) Setup I: The measurement campaign was conducted outdoors on a flight apron providing a largely controlled environment.Using the MEDAV RUSK-DLR wideband channel sounder in 1 × 4 single-input multiple-output (SIMO) mode, we measured repeatedly channel impulse response (CIR) snapshots of four radio channels.The measurement setup consisted of one transmitting and four receiving antennas statically arranged as shown in Fig. 2a.All antennas were placed 1 m above the ground.We used four similar toroidal, omni-directional receive antennas.To avoid an elevation of the noise floor due to limits of the dynamic range, we used a directional transmit antenna in order to reduce the received signal power of the LoS component in each channel [33].The orientation of the transmit antenna is indicated by the light gray radiation pattern in Fig. 2a.For inducing variations on the received signal power, different users were moving individually within the measurement environment.In this campaign, pedestrians, bicyclists, and cars were considered as users.By using a tachymeter in conjunction with a highprecision reflector prism, we recorded the ground truth of the moving users.The reflector prism was attached to a helmet worn by the pedestrians and bicyclists, and was mounted centrally on the roof of the cars, respectively.Fig. 2a illustrates the trajectories of the pedestrian user as an example.The trajectories of the bicyclists and cars were comparable.In total, we collected more than 2 100 000 CIR snapshots for pedestrian users, and more than 1 300 000 each for bicycle and car users.
2) Setup II: The measurement campaign was also conducted outdoors, but at an intersection on the DLR campus.The setup was openly accessible, i.e., the environment was less controlled.In addition, the propagation environment, with nearby buildings, parked cars, and vegetation, resembles that of a typical urban intersection.Similarly to Setup I, the MEDAV RUSK-DLR wideband channel sounder was used in 1 × 4 SIMO mode collecting CIR snapshots of four channels.The measurement setup is shown to scale in Fig. 2b.Again, the antennas were placed 1 m above the ground.Here, however, we used the same toroidal, omni-directional antennas at the transmitter and receivers to account for MPCs from all directions.In this setup, one moving pedestrian is considered, which induces variations to the received signal.Ground truth was recorded using a u-blox F9R global navigation satellite system (GNSS) receiver connected to an antenna attached to the pedestrian's helmet, similar to the prism in Setup I. Using real-time differential GNSS corrections, the receiver computed a GNSS real-time kinematic solution.Due to the fair open-sky environment, a sub-decimeter accuracy can be expected.The trajectory of the pedestrian is also shown in Fig. 2b.In total, more than 182 000 CIR snapshots were recorded. 1) Setup III: The measurement campaign was conducted indoors in a controlled environment.Using an UWB measurement system, based on the commercial off-the-shelf DecaWave DW 1000 chip, we collected the CIR snapshots for seven channels.Thereby, the measurement setup consisted of one transmitting and seven receiving UWB nodes as shown in Fig. 2c.For measuring the CIR, the receiving nodes were individually addressed by the transmitting node in a round-robin manner.All nodes were mounted 1 m above the ground.In this campaign, one moving pedestrian was considered, where the ground truth was recorded using a Vicon high-precision optical motion capture system.Therefore, a reflector was attached on top of the head of the pedestrian, which was tracked by the motion capture system.The Vicon motion capture system is capable of tracking the motion of the reflector with a subcentimeter accuracy.The trajectory is shown in Fig. 2c.In total, we collected more than 2800 CIR snapshots.

B. Parameter Estimation
For each channel measured during the measurement campaigns presented above, we must first determine the specific propagation effects of the static environment.In each campaign, we therefore collected additional CIR samples over a period during which the environment was devoid of any user.For each of these measured CIRs, we can now determine the signal parameter for N l separable MPCs, i.e., the amplitude and delay values as defined in (1), using maximum likelihood estimation.Specifically, we use the space-alternating generalized expectation-maximization (SAGE) algorithm for parameter estimation [35].Note that the measured channels individually addressed, round-robin Transmit antenna 9 dBi directional [30] 8 dBi toroidal, omni-directional [31] 0 dBi omni-dirctional [29], [32] Receive antennas 8 dBi toroidal, omni-directional [31] 8 dBi toroidal, omni-directional [31] 0 dBi omni-dirctional [29], [32] Reference system tachymeter GNSS optical motion capture system Environment outdoors outdoors indoors User types pedestrian pedestrian pedestrian bike car of each campaign correspond to network links of an equivalent MDFL system, which allows the use of the same notation as introduced in Sec.II-A.Subsequently, we average the estimated signal parameters over time and obtain a mean amplitude ᾱl,n and a mean delay value τl,n for each MPC n of channel l.The mean amplitude values allow to calculate the power of the MPCs for the idle channel.In logarithmic domain, the power of the n-th MPC of channel l is which can be used as reference power level for measuring user-induced power changes.
For this purpose, we determine the amplitude values for all MPCs of each channel.Given the mean delay τl,n for the n-th MPC of channel l, we estimate the amplitude using that is, the projection of the residuum signal y res l,n (t) onto the unit transmit signal s l (t) [17], [23].Adjusting the received signal for all MPCs up to the (n − 1)-th, the residuum signal is calculated as Given the amplitude estimate, we can then express the measured power in logarithmic domain as γl,n = 20 log 10 |α l,n | .
By subtracting the reference power level in (21) from the measured power in (24) we can finally express the userinduced power changes as In order to describe the location dependent impact of a user on the power of an MPC, we must next determine the actual propagation path of the MPC.Therefore, for each measurement setup, we have accurately measured the locations of distinct reflective objects in the propagation environment, including walls, lampposts, and fences, as well as vegetation and furniture.Using additional precise floor plan information, we obtain a detailed map of possible locations of RPs.Applying ( 13) and ( 18), we then construct the delay ellipse for τl,n , as shown exemplary in Fig. 1.By calculating the tangent or intersection points of the ellipse within the previously defined map, we determine potential RPs of the MPC.In this way, we determine all uniquely resolvable RPs for the MPCs of each measured channel.As examples, Fig. 9a and Fig. 9b illustrate the resulting propagation paths for Setup I and Setup III, respectively.Please note that while mandatory for the empirical modeling of this section, the manual determination of RPs is extremely cumbersome and can only be performed in postprocessing.Flexible adaptation to a changing environment is thus basically impossible, which further motivates for an automated RP estimation and mapping approach, as presented in Sec.IV.

C. Empirical Exponential Model
After collecting and processing the measurement data, we can now investigate the location dependent impact of a user on the power of an MPC.Therefore, we must first describe the relation between the user location r and the propagation path of an MPC.A common measure that relates the user's location to a signal propagation path is the excess path length [10], [14], [17].For calculating the excess path length of an MPC, we need to geometrically decompose the propagation path by combining successive pairs of virtual and physical nodes (see Sec. II-B).Thus, the excess path length of an MPC caused by SBR is calculated as for the link between physical transmitter and VR, see blue line in Fig. 1, and as for the link between VT and physical receiver, see red line in Fig. 1.Following [17], we approximate the user impact for the two pairs of nodes individually by the empirical exponential model (cf.[10]).Therewith, the user location dependent power changes of an MPC are modeled as with φ l,n and κ l,n being the maximum modeled power change in dB and the spatial decay rate.Note that in practice the model parameter can be determined individually for each link and each MPC.For the following empirical derivation and evaluation of the model parameter, however, we use the measured power changes of all MPCs from all channels to obtain sufficient statistics.Accordingly, Fig. 3a illustrates the power changes measured for all MPCs of each channel, induced by the pedestrian in Setup I. Additionally, the figure highlights the exponential fading model of ( 28) with fitted model parameters.Note that, for illustration purposes, we use the minimum excess path length for the abscissa of Fig. 3a and Fig. 3b.For SBR, the minimum excess path length is defined as with excess path lengths ξ Tx l,n (r) and ξ Rx l,n (r) defined in ( 26) and (27).See Fig. 9a for details on the propagation paths of Setup I used for Fig. 3. Finally, Table II provides a summary of the fitted model parameter for all setups.
Originating from the lognormal shadowing model [36], the measured power in logarithmic domain is very often modeled by a Gaussian random variable, e.g., in [10], [14], [37].Correspondingly, the measured power changes defined in (25) are simply modeled as with the zero-mean Gaussian noise w l,n ∼ N (0, σ2 l,n ).For Setup I, the residual values, i.e., the difference between the fitted exponential fading model (28) and the measured power changes (25), are illustrated in Fig. 3b.Note that unlike the work in [10], we could not statistically confirm that the residuals follow a Gaussian distribution with zero mean, although the overall shape of the residuals would be comparable (cf.Fig. 7 in [10]).Instead, we can clearly observe a location dependent variance as shown in Fig. 3b.In particular, the variance of the residuals increases with decreasing excess path length, i.e., the closer the user is to the propagation path, the higher the variance of the measured MPCs.Accounting for the location dependence of the variance, we modify (30) as with the location dependent and zero-mean Gaussian distributed noise term w l,n (r) ∼ N (0, σ 2 l,n (r)).Consistent with the work in [38] and [18] defining detection and sensitivity areas near LoS paths, we distinguish between areas near and far from the propagation path of an MPC for modeling the variance.Specifically, we use the definition of the excess path length describing an ellipse with foci at (virtual) transmitting and receiving nodes, see (26) and (27).By introducing a threshold ξ th for the excess path length, we therefore achieve an inherent geometric separation, i.e., between locations near the propagation path with excess path lengths below the threshold and locations far from the propagation path with excess path lengths above it.Thus, we model the standard deviation σ l,n (r) piece-wise as where the additive term ∆σ l,n > 0 reflects the increasing variations in power when a user is near a propagation path and implies that σ 1 l,n < σ 2 l,n .As mentioned above, the excess path length geometrically defines an ellipse with the transceiving nodes as foci.Similarly, Fresnel zones are defined as a series of discretized concentric ellipses, also with the transceiving nodes as foci [39], [40].Therefore, the definition of the Fresnel zones is basically equivalent to that of the excess path length, which allows to define a physically motivated threshold ξ th as with wavelength λ and the number of the Fresnel zone n F .The corresponding maximum radius of the Fresnel zone is given by On the basis of the measurement data for pedestrian users, we can observe that the standard deviation remains stable for excess path lengths above a threshold corresponding to the third Fresnel zone, i.e., n F = 3.As an example, Fig. 3b shows the residuals, the threshold, and the course of the standard deviation for Setup I.While the standard deviation of the residuals corresponding to an excess path length below the threshold increases, the standard deviation of the residuals corresponding to an excess path length above the threshold does not substantially change.Since the course of the standard deviation is qualitatively comparable for all measurement setups, the threshold calculated with n F = 3 also applies to Setup II and Setup III.Note that the value of the threshold for Setup III differs from that of Setup I and II due to the different center frequencies of the measurement systems (cf.Table I).Finally, Table III gives an overview of the standard deviation values calculated for the pedestrian users from the measurements of Setups I-III.As expected, the values for the location dependent noise model clearly indicate that the standard deviation is considerably higher when the user is close to the propagation path.This is visually confirmed by the two PDFs of the residuals in Fig. 3c, which are obtained for excess path length above and below the threshold ( 1 and 2 ).The good agreement between the empirical PDFs with the corresponding fitted normal distributions further confirms the approximation of the noise by the piece-wise defined standard deviation in (32).

D. Discussion
The fitted model parameters for all measured setups are provided in Table II and the corresponding values of the additive noise in Table III.In the following, we would like to comment on two observations in more detail.
First, when comparing the model parameters for the different user types of Setup I (cf.Table II), we can see that the values of both the maximum modeled power change φ and the decay rate κ are highest for the car, followed by the bicycle, and lowest for the pedestrian.Accordingly, the values of the model parameters are related to the dimensions of the user, which is consistent to theoretical models based on diffraction theory, e.g., [11], [18], [41] for LoS paths or [16] for MPCs.Since the dimensions of the car are larger than those of the bicycle, which in turn are larger than those of the pedestrian, we can observe that the absolute induced fading on the received power increases with the dimensions of the user.The increase in decay rate can be also explained by the user dimensions, as users are modeled empirically as point masses.Since perturbations in the received power can be measured when the edge of the user approaches the propagation path, the distance of the center point of the user to the propagation path increases with the user's dimensions.This relation is considered by the model with an increasing decay rate.Hence, the threshold on the excess path length (33) needs to be determined individually for different types of users.In this work we have determined the threshold on the excess path length for pedestrian users only.
Second, regarding the maximum modeled power change for pedestrian users (cf.Table II), we observe higher fitted values for Setup II and Setup III compared to Setup I. To explain these values, we refer to the basics of diffraction theory [39].The impact on the received power can be explained by the area of the first Fresnel zone blocked by the user [38], [39].
Due to the ellipsoidal definition of the Fresnel zone with the maximum radius depending on the length of the propagation path (34), the intensity of the induced fading depends where the user moves through the propagation path.As an example, Fig. 4 shows the power changes induced by a pedestrian user versus the user's distance to the respective nodes.It can be seen, that the user induces a stronger attenuation when being close to the nodes, i.e., when the user relatively blocks a larger part of the Fresnel zones.This applies both to the path between transmitting node and RP and between receiving node and RP.Looking at the user trajectory in Setup II, cf.Fig. 2b, we see that for the most part the user circled tightly around the network nodes, resulting in comparatively stronger attenuation, which is reflected in the model parameters.For Setup III, the smaller dimensions of the network geometry lead to comparatively smaller maximum Fresnel radii (34).Therefore, the area blocked by the user is relatively larger resulting in a higher maximum modeled power change.
Referring again to the results in Fig. 4, we can see that the power changes are comparable both qualitatively and quantitatively for both parts of the propagation path.This is an important finding since the user fading is thus independent of which part of the propagation path is affected.That means, it does not matter whether the user crosses the path between transmitting node and RP or between receiving node and RP.This independence of user-induced fading from the respective parts of the propagation path is a prerequisite for the geometric decomposition and thus fundamental to the fading model described in Sec.III-C.Overall, we have therefore introduced a computationally efficient yet accurate exponential fading model for MPCs that additionally accounts for user location dependent noise through a physically motivated spatial segmentation.

IV. BAYESIAN MAPPING OF REFLECTION POINTS
The localization performance of an MDFL system is severely depending on the information about the propagation paths within the network.While the propagation paths for the LoS components are given a priori with the locations of the network nodes, the propagation paths of MPCs still have to be determined.Therefore, in the following we present a sequential Bayesian estimation approach that determines the location of an RP that ultimately defines the propagation path of an MPC.

A. Problem Formulation
For estimating the RPs we can use the geometric properties of MPCs caused by SBR as described in Sec.II-C.That is, the possible locations of an RP are described by a delay ellipse, cf.(13).Thus, we can uniquely describe the location of the RP using the arc length as defined in (14).For the n-th MPC of network link l, we can therewith define a one-dimensional (1D) state at time instant k as with the arc length s l,n (θ k ).Since the RPs can be estimated independently, we omit the indices for MPC and network link in the following description for notational convenience.Considering a static MPC, we expect the location of the RP to be stationary, i.e., the state x k is time-invariant.The transition prior distribution is thus defined as where δ(•) denotes the Dirac delta distribution.For calibrating, i.e., initially estimating the locations of the RPs, the MDFL system continuously measures the variations in the power of the MPC, cf.(25).Assuming a known user location r k , we can model the measured power change equivalently to (31) as assuming zero-mean Gaussian distributed measurement noise w(r k ) ∼ N (0, σ 2 (r k )) with σ(r k ) as defined in (32).Applying (28), we can express the power change depending on the location of the RP as with maximum modeled power change φ and decay rate κ.The user state dependent excess path lengths are defined similarly to ( 26) and ( 27) as and where r Tx and r Rx denote the locations of the transceiving nodes, respectively.The geometric length of the propagation path is determined by the propagation delay τ of the considered MPC as in (2), i.e., d = c τ .The locations of the virtual nodes are calculated according to (7) and ( 8) as and The location of the RP r RP (x k ), expressed in the coordinate system used for the MDFL network, is calculated using (18), i.e., by rotating and translating the location of the RP r RP (θ(x k )), expressed in the local ellipse-centric coordinate system.Finally, we can calculate the ellipse parameter depending on the user state by the inverse function of the arc length defined in (14) as Note that there is no closed form solution for the inverse function of the arc length, i.e., the ellipse parameter needs to be determined numerically [42].

B. Elliptic Probability Density Function
Due to the 1D definition of the state x k , i.e., as elliptical arc length, cf.(35), the corresponding PDF of state x k is circularly distributed on the arc of the delay ellipse.In the following we outline the characteristics of that elliptic PDF.Therefore, Fig. 5 illustrates different examples for the elliptical PDF using the introductory example of an MPC due to SBR (see Fig. 1).As can be seen, the PDF can take four different shapes depending on the user's location, which are explained by the measurement model (37) together with the definitions of the excess path lengths in (39) and (40).Consider the example of a measurement corresponding to an excess path length of ξ Tx (x k , r k ) = ξ Rx (x k , r k ) = 0, i.e., the user is located on or very close to the propagation path, as shown in Fig. 5 for user location 1 (blue area).According to the excess path length, we know that the user needs to be located between the RP and one of the two transceiving nodes.As illustrated for user location 1 , there are exactly two possible locations for the RP.Namely, when the user is located between r Tx and the location of the true RP or, equally likely, when the user is located between r Rx and a second location of the RP.In the same way, we can explain the three remaining shapes of the PDF.Thereby, for non-zero measurements, the most general case is given for user location 3 (green area).Here, a measured power change corresponds to a non-zero excess

Initialize weights
; Normalize weights; Thus, for each part of the propagation path, i.e. between r Tx or r Rx and an RP, two individual RP locations can be determined, resulting in a total of four possible locations.This is also true for user location 2 (orange area), but due to symmetry, for each part of the propagation path, one of the two determined RP locations corresponds to the true RP, resulting in three possible locations for the RP.Finally, a user at location 4 (gray area) does not affect the received power significantly.Due to the shape of the exponential model, no unique locations for the RP can be determined when the measurement is close to zero.But, depending on the user position, it is possible to distinguish between areas on the elliptic arc in which the RP can be located and in which the RP can be excluded.

C. Sequential Bayesian Filtering
Using a sequential Bayesian estimator, the PDF of the state x k is determined by computing the posterior density p(x k |z 1:k ) applying the general Bayesian update recursion [43].While for linear system models, such as linear Gaussian systems, the posterior can be efficiently estimated using Kalman filtering solutions, for nonlinear systems the posterior must be numerically approximated in most cases [44].For the estimation problem outlined in Section IV-A, the nonlinearities are due to the exponential measurement model (37) as well as the state-dependent definition of the excess path lengths (39) and (40).Furthermore, due to the properties of the PDF, as discussed in Section IV-B, the posterior can hardly be assumed to be Gaussian.Possible filter solutions for such nonlinear and non-Gaussian processes that numerically approximate the posterior are given by the particle filter (PF) and the point mass filter (PMF) [43], [45].Given the elliptic PDF of the state x k , a PF can lead to increased complexity, since according [46] sampling from circular distributions can be costly.In contrast to the PF, the PMF does not require any resampling.With respect to complexity, we choose the PMF for estimating the posterior.Therefore, the PMF approximates the posterior distribution with the discrete density where x i represents the i-th grid point of the deterministic grid {x i } Ns i=1 [43].Due to the definition of the state x k as arc length, the grid points can be easily arranged equidistantly around the arc of the ellipse.The grid point spacing is determined by the circumference of the ellipse L, cf.(17), and the number of grid points N s as ∆ x = L/N s .The weights w i k|k are calculated as with the normalization term c k = Ns j=1 w j k|k−1 p(z k |x j , r k ) and the likelihood distribution p(z k |x i , r k ), which expresses the measurement model in (37).The predicted weights are where p(x i |x j ) refers to the transition prior distribution.
Because of the approximation of the posterior density by a discrete density, i.e., by a finite number of stationary grid points, a direct application of the Dirac delta function of ( 36) can lead to estimation problems similar to the problem of loss of diversity [45].Therefore, we approximate the transition prior by an elliptic normal distribution with η as concentration parameter and L as elliptic circumference (17).I 0 (•) is the modified Bessel function of the first kind and order 0. Please refer to the derivation of the elliptic normal distribution in the Appendix.Note that applying the elliptic normal distribution to the discretized grid points around the elliptic arc equals a circular convolution [43], [46].In order to obtain a point estimate xk from the elliptic posterior density of (44), we compute the minimum mean square error (MMSE) estimate.Finally, Algorithm 1 provides a pseudocode for the sequential Bayesian filtering approach described above.

V. EXPERIMENTAL EVALUATION
In this section, we evaluate the above presented Bayesian mapping approach of RPs.Therefore, we demonstrate the applicability of the approach using measurement data.First, for an exemplary MPC for the link between Tx and Rx 1 of Setup II.And second, for estimating all resolvable MPCs of Setup I and Setup III, thus showing the applicability for different environments and different measurement systems.Prior to the subsequent evaluation, we define the evaluation metric.For evaluating the RP location estimation, we need to consider an error measure taking into account periodicity.Analogously to the shortest distance on a circle [46], we define the shortest distance on an ellipse as where L is the circumference of the ellipse, x RP is the arc length of the true RP, and xk is the point estimate of the arc length from the PMF.First, we demonstrate the Bayesian RP mapping approach with measurement data from Setup II.Details on the measurement data are given in Table I.Specifically, we consider the link between Tx and Rx  9) and (10), respectively.Using (17), the circumference is calculated as L = 97.633m.An overview of the considered measurement environment is provided in Fig. 6, including the transmitting and receiving nodes and the aforementioned delay ellipse.The true RP and the propagation path are shown for illustration only.Furthermore, Fig. 6 shows the calibration trajectory of a pedestrian user, here determined using GNSS (cf.Table I).Note that during calibration the user is assumed to provide the position information to the MDFL system.The power changes of the MPC measured while the user was moving are shown in Fig. 7. Thereby, we can clearly observe three time intervals, labeled as t 1 -t 3 , of increased variations of the received power.The user positions that correspond to these intervals are indicated by arrows in Fig. 6.It can be clearly seen, that the time intervals of increased variations in the power coincide with the times when the user passes through the propagation path.Further, we have pointed to a time instant t 0 in both Fig. 6 and Fig. 7.As shown, the user is located on the delay ellipse at t 0 .That means that at this user location, the delay of a directly scattered signal coincides with that of the considered MPC, which explains the variations in the received power around t 0 .These variations are independent of the RP that we want to estimate, and thus would certainly affect the performance of the RP estimation.Since we know the user location during calibration, we can easily mitigate   log(e RP ) [m] Fig. 8: Distance error of estimated RP location over time, for PMF realization with noise model of ( 30), ( ), and with user location dependent noise model ( 32), ( ).The distance error is given in logarithmic domain.
such potentially deceptive measurements when the user is near the delay ellipse.
For estimating the RP, we realize the PMF as introduced in Section IV-C.The grid points of the PMF are uniformly distributed on the elliptical arc with a spacing of ∆ x = 5 cm, i.e., a total of N s = 1953 grid points.Initially, the weights of the PMF are set equal.For the measurement model (37), we choose the parameters φ = −2.5 dB and κ = 0.015 dB.For the weight prediction (46), we use a concentration parameter of η = 1/0.01mm 2 for the elliptic normal distribution.As measurement noise, we consider both the normal distributed (30), with σ = 0.75 dB, as a reference, and the user location dependent normal distributed (32), with σ 1 = 0.5 dB and σ 2 = 1.0 dB.
The resulting distance error over time is shown in Fig. 8. Initially, the weights of the PMF are set equal, i.e., the posterior PDF corresponds to a uniform distribution.That uniform distribution in turn explains the initially large distance error.When the user passes the propagation path for the first time (t 1 ), we can observe that the distance error immediately decreases from the initial value.Here, the PMF realized with the location dependent noise model achieves a higher performance gain when the user passes the propagation path the first time.When the user passes the propagation path the second time (t 2 ), also the PMF realization using the noise model of (30) converges.The results of both PMF realization achieve a similar distance error.After circling around the receive antenna (t 3 ), the distance errors of both realizations further decrease.With a final distance error of e RP,end = 12.8 cm using the noise model of (30), and of e RP,end = 3.5 cm using the user location dependent noise model of (32), both realizations achieve an accurate estimation result.
In order to further demonstrate the applicability of the Bayesian RP mapping approach to measurement data in dif- ferent environments and for different measurement systems, we apply the approach to the resolvable MPCs of Setup I and of Setup III.Therefore, we use the same PMF realization as introduced above.Only differences are the number of grid points, that are adjusted individually for each MPC, due to different ellipse circumferences and a fixed grid point spacing.And second, the concentration parameter of the elliptic normal distribution is adapted for the UWB measurement system.For UWB, we set η = 1/3.24cm 2 accounting for the increased update time T g , see Table I.Further, we choose the normal distributed noise model of (30) for simplicity.The calibration trajectories for both setups are shown in Fig. 9a and Fig. 9b, respectively, along with the resulting RP location estimates.The distance error for each RP is visualized in red.Overall, the estimated RP locations agree very well with the true locations in both setups.In particular, we can not identify any substantially incorrect RP location estimates.The averaged distance errors over all considered MPCs of the respective setups, i.e., ēRP = 0.82 m for Setup I and ēRP = 0.44 m for Setup III, support the visually observed robust estimation performance.Note that the higher distance error of Setup I can be explained by the larger dimensions of the environment of Setup I compared to Setup III, which result in larger absolute errors.In summary, the Bayesian RP mapping approach was thus successfully applied to different measurement systems and different environments.

VI. CONCLUSION
This paper provides a novel Bayesian calibration approach for multipath-enhanced device-free localization (MDFL) systems that determines the propagation paths of multipath components.Therefore, we first present a statistical fading model that describes user induced changes in the received power of a multipath component.The model is thus not only important for calibration, but also an essential building block for MDFL systems that infer the location of a user based on the induced changes in the received power of both multipath components and line-of-sight signal components.Based on an extensive set of wideband and ultra-wideband (UWB) measurement data for both indoor and outdoor environments, we derive and validate the model empirically.As emerges from the measurement data, the presence of users not only attenuates the received power, but also leads to increasing fluctuations.Thus, in addition to the empirical exponential fading model, we propose a location dependent variance model for the measurement noise using an efficient, physically motivated spatial segmentation.
Second, we present the novel Bayesian calibration approach for MDFL that robustly estimates the locations of reflection points from single-bounce reflections.MDFL systems severely depend on the information about the propagation paths within the network.While known for the line-of-sight, the propagation paths have yet to be determined for multipath components, which is ultimately achieved by determining corresponding reflection points.Using the initially derived empirical fading model and given user location during calibration, we can relate measured changes in the received power of the multipath component to the location of a reflection point.Taking advantage of geometrical properties of multipath components caused by single-bounce reflections, we can constrain the possible locations of reflection points to the locations of the corresponding delay ellipse, which allows to formulate a one-dimensional elliptic estimation problem.We propose the point mass filter that approximates the elliptic posterior density by a deterministic grid for efficiently solving the estimation problem.The applicability of the presented approach is demonstrated and evaluated using measurement data for different environments and different measurement systems, including a commercial off-the-shelf UWB system.The Bayesian calibration approach shown to robustly estimate the locations of the reflection points.Independent of the considered environment or measurement system, the estimation approach achieves a sub-meter accuracy.

APPENDIX DERIVATION OF ELLIPTIC NORMAL DISTRIBUTION
The von Mises distribution, also denoted as circular normal distribution, describes a probability distribution around a circle.Its equation is given by [47] p(α|µ, κ) = 1 2πI 0 (κ) exp (κ cos (α − µ)) , α ∈ [0, 2π) , (49) with concentration parameter κ, location parameter µ, and modified Bessel function I 0 (•) of the first kind and order 0. The von Mises probability density function (PDF) is defined on a unit circle.We want to transform the PDF of α to circles of arbitrary length with circumference L c .Therefore, we define s c as the length on the circle corresponding to the angular parameter α, with the variable substitution The Jacobi element, which is needed for the correct normalization factor according to [48] is given by In the following, we want to extend the distribution to an ellipse with circumference L e , cf. (17).With the definition of an arbitrary arc length of the ellipse of we can define, similarly to (50), the variable substitution

Fig. 2 :
Fig. 2: Measurement setups showing transmitter and receiver locations ( ), user trajectories (), and distinct reflection surfaces, e.g.building walls ( ).In (a), the shape of the directional antenna pattern around r Tx is highlighted in light gray.The dotted rectangles in (a) and (b) refer to the map sections used for Fig.4and Fig.6, respectively.

Fig. 3 :
Fig. 3: Measurement data for pedestrian user in Setup I: (a) measured power change values and fitted exponential fading model; (b) corresponding residual power values along with standard deviation (enveloping gray dashed lines), and Fresnel-motivated threshold for excess path length (vertical dashed line); (c) PDFs of residuals for excess path length above and below the threshold ( 1 and 2 ), where dashed lines in black indicate the fitted normal distribution with σ 1 = 0.2813 dB and σ 2 = 0.7957 dB (cf.TableIII).

Fig. 4 :
Fig.4: Exemplary specular reflection for the nodes at r Tx and r Rx 3 of Setup I and averaged power changes induced by the pedestrian user when located on the propagation path depending on the user's distance to r Tx ( 1 ) or r Rx 3 ( 2 ).

Fig. 5 :
Fig. 5: Exemplary network link and propagation path of an MPC due to single-bounce reflection, cf.Fig. 1.Areas of characteristic shapes of the PDF are indicated in color.Corresponding PDFs on the arc of the delay ellipse are provided for user locations 1 -4 .The user locations are indicated by ×, the true RP is indicated by black diamond, and the PDF is highlighted in blue.

Algorithm 1 :
Bayesian Mapping of RP at Time Step k Input: Measured power change: z k User location: r k Grid points of PMF: {x i } Ns i=1 Weights of PMF: {w i k−1|k−1 } Ns i=1 Output: Weights of PMF: {w i k|k } Ns i=1 MMSE estimate: xk 1 if k=0 then 2

Fig. 6 :
Fig. 6: Overview of measurement environment.Grid points located on delay ellipse indicated in blue, true RP indicated by black diamond marker, corresponding propagation path indicated in gray, and trajectory including motion direction highlighted in green.
1 with a distance of d LoS = 31.370m and an MPC with an estimated propagation distance of d = 38.673m.As shown in Section II-C, possible locations for the RP are determined by an ellipse, cf.(13), with principle axis of a = 19.337m and b = 11.308m, see (

Fig. 7 :
Fig.7: Measured power changes in dB for an MPC with propagation distance of 38.673 m.Noticeable variations are highlighted in gray and the corresponding times labeled as t i ,i ∈ {0, 1, 2, 3}.
Fig. 9: Resulting RP location estimates: true propagation paths are shown as gray lines, estimated RP locations are indicated by diamonds, and calibration trajectory is represented by dashed green line for each setup; distance between estimated and true RPs are highlighted by red lines; averaged distance error of RP estimates for (a) Setup I: ēRP = 0.82 m and for (b) Setup III: ēRP = 0.44 m.

TABLE I :
Overview measurement setups and parameters

TABLE II :
Fitted model parameters

TABLE III :
(32)dard deviation of noise models for pedestrian user according to(30)or(32)