Outdoor localization, thanks to GPS technology, has revolutionized navigation-based applications running on e.g. automotive GPS enabled devices and smart phones.  Applications range from guiding drivers to their destination to providing a point-by-point direction to the closest cinema or coffee shop.  The success of GPS can be attributed to the reliability, availability of line-of-sight (LOS) and practical accuracy that the orbiting satellites covering the globe provide. 
Recently as the number of smart devices and mobile users has increased significantly, the potential for new indoor/urban location-based services and challenges has emerged.  The lucrative business opportunities of location-enabled services are not limited to outdoors.  In fact the potential of indoor location-enabled services has been projected by a number of sources as an untapped industry.  The range of potential indoor applications touches every aspect of our lives: from tracking children in day-care centers, tracking elderly in nursing homes, to tracking inventories in warehouses, tracking medical devices in hospitals and tracking personnel in emergency/first responder applications (e.g. firefighters).  The major challenges facing this emerging industry are non-existent GPS coverage, multipath propagation and the non-line-of-sight (NLOS) problem. 
Position estimation is typically achieved through ranging (distance estimation or angle estimation) to different base stations (BSs) with known coordinates and solving a set of non-linear equations or through geometric relations.  Geolocation technologies that depend on techniques such time of arrival (TOA), time difference of arrival (TDOA) and Angle of Arrival (AOA) require a clear LOS between the BSs (the orbiting satellite in GPS) and the mobile device (MD). 
In low system bandwidths, multipath propagation can induce tens of meters of range errors; while NLOS propagation causes significant errors/biases (tens/hundreds of meters) to distance/angle estimation which affects the localization accuracy directly.  Thus to enable accurate and robust localization in indoor/urban environments it is important to address and mitigate the biases that affect the range/distance estimation. 
Existing NLOS bias mitigation literature can be generally sub-divided into two main groups.  The first is NLOS identification and mitigation where the range measurements to different BSs are classified into LOS/NLOS and mitigated for.  The second is NLOS bias tracking and correction which generally assume dynamic movement of the MD. 
In the former approach the NLOS measurements are either ignored or are integrated in a weighted location optimization approach such as weighted least squares (WLS) or weighted constrained optimization algorithms [1]. 
In [2] a bias tracking algorithm is proposed where the biases are incorporated into a Kalman Filter (KF) formulation but this approach requires an a priori knowledge of the bias covariance matrix; which in realistic propagation conditions is difficult to obtain due to rapid fluctuations of the bias resulting a major weakness. 
In [3] a statistical bias correction technique that works on the range measurements prior to the KF tracking stage is proposed.  It is based on statistical processing of a record of measurements taken over a window but it relies on the statistical estimate of the NLOS measurement ratio present in the record.  As a result it suffers from several flaws.  First it assumes that there is zero bias in LOS environments; which is not true.  Second it assumes that the statistics of the biases to be time-invariant within a short period which might not be the case in highly dynamic scenarios.  Third, the approach requires a priori statistical characterization of the biases which is difficult to maintain in a highly changing dynamic environments such as indoors. 
In [4], an extended KF (EKF) tracking is proposed with the state matrix augmented by the unknown bias vectors.  The approach requires prefect knowledge of LOS/NLOS identification which is highly unrealistic; given that NLOS identification is a difficult prospect on its own and perfect identification is rarely, if ever, achieved.  In addition the incorporation of the biases into the state vector requires estimating the bias covariance in real-time which is difficult in highly dynamic environments where the bias fluctuates significantly. 
Another approach based on an improved KF is described for mobile tracking in NLOS environments where a smoothing stage is introduced to suppress NLOS errors in addition to integrating velocity and heading from motion sensors to improve the performance [5].  In addition to the zero bias assumption in LOS condition, the technique requires identification of LOS and NLOS through running mean/variance. 
A NLOS bias correction through a weighting mechanism in wireless position systems has been proposed in [6]. 
In [7] a constrained optimization approach is used to estimate the biases which are assumed to affect NLOS only.  The approach requires NLOS identification and the optimization algorithm is computationally expensive where it is incorporated with a tracking algorithm to estimate the position. 
In [8] a NLOS mitigation approach is proposed that is based on constrained optimization (quadratic programming) to estimate and mitigate the impact of the biases.  In additional to the computational complexity, the results of the simulations show that there is limited accuracy gain using this method. 
In [9] a KF and a sliding window are used to identify NLOS conditions and smooth range measurements to improve range estimation accuracy. 
In [10] statistical characterization of the NLOS biases (variance) is used to improve the accuracy of Kalman filter based algorithms. 
What is common from the literature discussed above is that the techniques are usually incorporated with a KF, assume zero bias in LOS and require some form of LOS/NLOS identification in order to implement a statistical approach to bias correction.  The zero bias in LOS and dependence on LOS/NLOS identification affects the practicality and robustness of the proposed techniques. 
TOA-Based Localization and Tracking
TOA and Distance Estimation
Central to any localization and tracking system are the measurements that infer the geometrical relationships between the set of base stations (GPS satellites in GPS Systems, WiFi Access Points (AP) in WLAN Systems) and the mobile device (MD).  Depending on the system used, geometrical information can be obtained through time measurements, angle measurements or both.  The most popular approach is time-based, where the propagation time, τ, between a base station (BS) and a MD is translated to distance through the speed of signal propagation equation
                    τ        =                  d          c                                    (        1        )            where d is the distance between the BS and the MD and c is the speed of signal propagation (usually assumed the speed of light—c=3e8 m/s).  When the MD is stationary it is possible to estimate the position through popular multilateriation techniques and the TOA (range/distance) from the MD to the BSs forms a set of intersecting circles; where the location of the intersection is the location estimate (x, y). 
FIG.  1 illustrates an example of TOA multilateration with 3 BSs 1.  In ideal scenarios and with accurate range measurements, the intersection of the circles provides the exact location of the MD 2.  In reality range measurements are corrupted with noise and biases, which expands the intersection point to an intersection region, thereby increasing the location uncertainty. 
The accuracy of the position estimation will depend directly on the accuracy of the range measurements and also the relative geometry of the BSs 1 and the MD 2.  If the range measurements are accurate and unbiased (zero mean and small measurement variance), then very accurate position estimation is possible.  There are different techniques to estimate the TOA, but the simplest utilizes a two-way ranging mechanism as illustrated in FIG.  2. 
Essentially Device 1 (the BS) initiates the two-way ranging by sending a ranging packet (signal) to Device 2 (the MD) and noting the time as tTX1.  Device 2 receives the signal at tRX2 and prepares its own ranging signal (after a processing delay) and sends out a response ranging signal at time tTX2.  Finally Device 1 receives the response at tRX1 Given that Device 2 shares the time stamp information tRX2 and tTX2 with Device 1 it is now possible to estimate the propagation delay (distance) between the two devices byτ=[(tRX2−tTX1)+(tRX1−tTX2)]/2  (2)
The main challenge using this approach is the delay inherent in software timestamping which has been reported to have a time resolution of 1-5 μs (0.3-1.5 km) [11].  Recently different IEEE standards such as IEEE 802.15.4, IEEE802.1as and IEEE 802.11v have been proposed to provide accurate and precise TOA-based ranging through physical layer amendments and hardware timestamping [12], [13].  The systematic errors due to time-stamping can be modeled as Gaussian and thus the range (distance) estimate can be modeled as{circumflex over (d)}=c{circumflex over (τ)}=d+n  (3)where n is a zero-mean Gaussian random variable with variance σ2; where the variance will depend on the time-stamping method used.  For hardware time-stamping, nanosecond resolution/accuracy can be achieved compared to microsecond resolution in software time-stamping.  Throughout the rest of the report we assume that one of the IEEE [12-17] standards of ranging/synchronization using hardware time stamping is used. The Kalman Filter
In practice the MD 2 is usually moving (along with the person) and the static estimation of the MD location becomes a dynamic state estimation problem.  For ideal zero-mean Gaussian conditions, an optimum solution to the problem is the Kalman Filter (KF) which estimates and tracks the location of the MD.  Typically the KF relies on two main equations: the state equation and the measurement equation.  The state equation describes the recursive relationship between the current state estimate and the next state estimate.  The measurement equation describes the relationship between the available measurement (range/distance measurements) and the state.  Depending on the motion model assumed, the state vector usually contains the coordinates of the MD (x, y) and the velocity (vx, vy).  A popular motion model (which covers the most generic case) is the constant velocity model; where it is assumed that the MD moves in a constant velocity and the state equation is then given by [14]x(k+1)=F(k)x(k)+G(k)w(k)  (4)where x(k)=x(k|k)=[x(k), vx(k), y(k), vy(k)]T is the state vector at the current sample time tk,
      F    ⁡          (      k      )        =      [                            1                                      T            s                                    0                          0                                      0                          1                          0                          0                                      0                          0                          1                                      T            s                                                0                          0                          0                          1                      ]  is the state transition matrix, and
      G    ⁡          (      k      )        =      [                                                      T              s              2                        /            2                                    0                                                  T            s                                    0                                      0                                                    T              s              2                        /            2                                                0                                      T            s                                ]  is the noise gain matrix and w(k)=[wx (k),wy(k)]T is the process noise (acceleration noise—jitter in the constant velocity movements) where wx(k) and wy(k) are Gaussian distributed with zero mean and variance σ2 and Ts=tk−tk-1 is the sample time. 
The measurement equation which relates the range measurements to the state is given byz(k)=H(k)x(k)+n(k)  (5)where H(k) is the measurement matrix and n(k)=[n1(k) n2(k) . . .  nM(k)]T is the measurement noise where nm(k) is i.i.d.  zero-mean Gaussian with variance σn2, mϵ[1, M] and M is the number of BSs covering the MD.  The measurement equation in (5) assumes that the relationship between the state and the measurements is linear.  For the localization problem this is not the case since the measurement equation is non-linear due to the non-linear relationship between the distance and the state or
                              d          ⁡                      (            k            )                          =                  [                                                                                                                                        (                                                                              x                            1                                                    -                                                      x                            ⁡                                                          (                              k                              )                                                                                                      )                                            2                                        +                                                                  (                                                                              y                            1                                                    -                                                      y                            ⁡                                                          (                              k                              )                                                                                                      )                                            2                                                                                                                                                                                                              (                                                                              x                            2                                                    -                                                      x                            ⁡                                                          (                              k                              )                                                                                                      )                                            2                                        +                                                                  (                                                                              y                            2                                                    -                                                      y                            ⁡                                                          (                              k                              )                                                                                                      )                                            2                                                                                                                          ⋮                                                                                                                                                        (                                                                              x                            M                                                    -                                                      x                            ⁡                                                          (                              k                              )                                                                                                      )                                            2                                        +                                                                  (                                                                              y                            M                                                    -                                                      y                            ⁡                                                          (                              k                              )                                                                                                      )                                            2                                                                                                    ]                                    (        6        )            where (xm, ym), are the x- and y-coordinates of the mth BS.  Thus in order to use the KF, the measurement equation has to be linearized around the current state estimate {circumflex over (x)}(k) [14].  When either one or both the state and the measurement equations are non-linear then a linearized version of the KF, known as the extended KF (EKF), can be implemented.  The linearization is achieved through Taylor series expansion and usually the first order term is retained, while the higher order terms are neglected [14].  The linearized measurement matrix is then given by
                                          H            ⁡                          (              k              )                                =                      [                                                                                                      ∂                                                                        d                          1                                                ⁡                                                  (                          k                          )                                                                                      dx                                                                    0                                                                                            ∂                                                                        d                          1                                                ⁡                                                  (                          k                          )                                                                                      dy                                                                    0                                                                                                                        ∂                                                                        d                          2                                                ⁡                                                  (                          k                          )                                                                                      dx                                                                    0                                                                                            ∂                                                                        d                          2                                                ⁡                                                  (                          k                          )                                                                                      dy                                                                    0                                                                              ⋮                                                  ⋮                                                  ⋮                                                  ⋮                                                                                                                        ∂                                                                        d                          M                                                ⁡                                                  (                          k                          )                                                                                      dx                                                                    0                                                                                            ∂                                                                        d                          M                                                ⁡                                                  (                          k                          )                                                                                      dy                                                                    0                                                      ]                          ⁢                                  ⁢                              where            ⁢                                                  ⁢                                          ∂                                                      d                    m                                    ⁡                                      (                    k                    )                                                              dx                                =                                                                      x                  m                                -                                  x                  ⁡                                      (                    k                    )                                                                                                                                          (                                                                        x                          m                                                -                                                  x                          ⁡                                                      (                            k                            )                                                                                              )                                        2                                    +                                                            (                                                                        x                          m                                                -                                                  x                          ⁡                                                      (                            k                            )                                                                                              )                                        2                                                                        ⁢                                                  ⁢            and                          ⁢                                  ⁢                                            ∂                                                d                  m                                ⁡                                  (                  k                  )                                                      dy                    =                                                                      y                  m                                -                                  y                  ⁡                                      (                    k                    )                                                                                                                                          (                                                                        y                          m                                                -                                                  y                          ⁡                                                      (                            k                            )                                                                                              )                                        2                                    +                                                            (                                                                        y                          m                                                -                                                  y                          ⁡                                                      (                            k                            )                                                                                              )                                        2                                                                        .                                               (        7        )            
The main equation in the EKF is the state update equation [14]{circumflex over (x)}(k+1|k+1)={circumflex over (x)}(k+1|k)+W(k+1)v(k+1)  (8)where {circumflex over (x)}(k+1|k+1) is the estimated state, {circumflex over (x)}(k+1|k) is the predicted state, W(k+1) is the filter gain and v(k+1) is the measurement residual (also referred to as the measurement innovation) [14].  The predicted state is obtained from the state equation (4) in the absence of process noise (representing ideal state transitions) and it is given by{circumflex over (x)}(k+1|k)=F(k){circumflex over (x)}(k).   (9)
The measurement residual is the difference between the actual measurements and the predicted measurements and it is given byv(k+1)=z(k+1)−{circumflex over (z)}(k+1|k)  (10)where {circumflex over (z)}(k+1|k) is the predicted range measurement or
                                          z            ^                    ⁡                      (                                          k                +                1                            |              k                        )                          =                              [                                                                                                                                                        (                                                                                    x                              1                                                        -                                                                                          x                                ^                                                            ⁡                                                              (                                                                                                      k                                    +                                    1                                                                    |                                  k                                                                )                                                                                                              )                                                2                                            +                                                                        (                                                                                    y                              1                                                        -                                                                                          y                                ^                                                            ⁡                                                              (                                                                                                      k                                    +                                    1                                                                    |                                  k                                                                )                                                                                                              )                                                2                                                                                                                                                                                                                                      (                                                                                    x                              2                                                        -                                                                                          x                                ^                                                            ⁡                                                              (                                                                                                      k                                    +                                    1                                                                    |                                  k                                                                )                                                                                                              )                                                2                                            +                                                                        (                                                                                    y                              2                                                        -                                                                                          y                                ^                                                            ⁡                                                              (                                                                                                      k                                    +                                    1                                                                    |                                  k                                                                )                                                                                                              )                                                2                                                                                                                                          ⋮                                                                                                                                                                          (                                                                                    x                              M                                                        -                                                                                          x                                ^                                                            ⁡                                                              (                                                                                                      k                                    +                                    1                                                                    |                                  k                                                                )                                                                                                              )                                                2                                            +                                                                        (                                                                                    y                              M                                                        -                                                                                          y                                ^                                                            ⁡                                                              (                                                                                                      k                                    +                                    1                                                                    |                                  k                                                                )                                                                                                              )                                                2                                                                                                                  ]                    .                                     (        11        )            
The filter gain is then [14]W(k+1)=P(k+1|k)H(k+1)TS(k+1)−1  (12)which is a function of the linearized measurement matrix H(k), the state prediction covariance P (k+1|k) and the residual covariance S(k+1).  The state prediction covariance is given by
                                                                        P                ⁡                                  (                                                            k                      +                      1                                        |                    k                                    )                                            =                            ⁢                              cov                ⁡                                  (                                                            x                      ^                                        ⁡                                          (                                                                        k                          +                          1                                                |                        k                                            )                                                        )                                                                                                        =                            ⁢                                                                    F                    ⁡                                          (                      k                      )                                                        ⁢                                      P                    ⁡                                          (                                              k                        |                        k                                            )                                                        ⁢                                      F                    ⁡                                          (                      k                      )                                                                      +                                                      G                    ⁡                                          (                      k                      )                                                        ⁢                                      Q                    ⁡                                          (                      k                      )                                                        ⁢                                                            G                      ⁡                                              (                        k                        )                                                              T                                                                                                          (        13        )            where P (k|k)=cov({circumflex over (x)}(k|k)) is the state covariance and Q(k) is the covariance of the process noise and it is given by Q(k)=cov(w(k)).  The residual covariance is given byS(k+1)=R(k+1)+H(k+1)P(k+1|k)H(k+1)T  (14)where R(k)=cov(n(k)) is the measurement covariance. Localization and Tracking in Harsh Multipath Environments
Ranging and localization in urban and indoor environments face two major challenges that affect the performance.  The first is multipath and the second is the non-line-of-sight (NLOS) problem.  Both have received considerable attention recently in research and different approaches/algorithms have been proposed to enable more accurate and robust localization systems.  In this section the problems are discussed and an explanation is given of how these challenges introduce biases into the distance/range estimation which ultimately degrade the localization and tracking accuracy. 
The Multipath Problem
In a multipath propagation environment, the transmitted signal undergoes reflections, attenuations and diffractions prior to arriving at the receiver.  At the receiver, replicas of the transmitted signal arrive attenuated, phase shifted and time-delayed.  The multipath signal is a combination of those multiple signal arrivals given by
                              r          ⁡                      (            t            )                          =                              ∑                          l              =              1                                      L              p                                ⁢                                          ⁢                                    α              l                        ⁢                          e                              j                ⁢                                                                  ⁢                                  ϕ                  i                                                      ⁢                          s              ⁡                              (                                  t                  -                                      τ                    l                                                  )                                                                        (        15        )            where s(t) is the transmitted signal waveform, r(t) is the received waveform, where Lp is the number of MPCs, and αl, ϕl and τl are amplitude, phase and propagation delay of the signal traveling the ith path, respectively.  The received signal is essentially the transmitted signal convolved with the channel impulse response (CIR) or r(t)=s(t)*h(t).  Then it follows that the CIR is given by,
                              h          ⁡                      (            τ            )                          =                              ∑                          l              =              1                                      L              p                                ⁢                                          ⁢                                    α              l                        ⁢                          e                              j                ⁢                                                                  ⁢                                  ϕ                  i                                                      ⁢                          δ              ⁡                              (                                  τ                  -                                      τ                    l                                                  )                                                                        (        16        )            where δ(□) is the Dirac delta function. 
In LOS conditions, multipath signals (especially in dense cluttered environment) tend to arrive fairly close to the direct path.  If the inter-arrival time between the multipath components is much smaller than the time-domain resolution of the system (low bandwidth systems) then at the receiver the multiple signals will combine to create a cluster.  The TOA estimate (from the receiver's point of view) will then be the peak of the cluster.  FIG.  3 illustrates an example of a channel impulse response which shows how the multipath bias can corrupt TOA-based range estimates.  It is clear in this case that the actual TOA is not equal to the estimated TOA.  This difference in estimation is called the “multipath bias”. 
For higher system bandwidths, the multipath bias in LOS environments is usually smaller due to the higher time-resolution.  In LOS the distance estimate can be modeled by{circumflex over (d)}LOS=d+bLOS(ω,d)±nLOS  (17)where bLOS=bm (ω, d) is a bias induced by the multipath and it is a function of the system bandwidth and the distance between BS and MD [15].  nLOS is a zero-mean additive measurement noise.  For large bandwidths (such as in Ultra Wideband systems), the biases range from centimeters to a meter and are usually positive [15].  Both exponential and lognormal distributions have been reported in literature as possible fit to the data [15].  In this report we model the biases with an exponential distribution due to the simplicity of the model and its ability to simulate different range error behaviors.  The exponential model is given byfb(b|LOS)=bmin+λLOSe−λLOSb  (18)where λLOS is the rate parameter and bmin is a constant bias that can model the minimum experienced multipath bias in different indoor environments.  Clearly bmin in open areas is smaller than cluttered environments.  Also the minimum bias experienced in low bandwidth systems (WiFi) is much larger than that experienced in large bandwidth systems (UWB).The NLOS Problem
In NLOS conditions there is an obstruction in the path of the transmitter and receiver.  Depending on the type of obstruction and the relative distances of the transmitter/receiver to the obstruction, the channel behavior and ranging bias can vary significantly.  There are two specific NLOS cases or conditions that occur in typical obstructed environments.  The first is when the DP signal is attenuated but detected (albeit with weak SNR).  This situation can arise naturally when the transmitter and receiver are separated by “light” obstructions such as a glass door, wooden door or sheet-rock thin walls.  Indeed in this scenario TOA estimates can be obtained with good accuracy due to the detection of the DP signal.  The second NLOS case can occur when there is a “heavy” obstruction between the transmitter and receiver which attenuates the DP severely making it difficult for the receiver to detect it.  The first non-Direct path (NDP) component is then used for TOA estimation.  This results in a significant bias that corrupts the TOA estimation and ultimately the position estimate.  FIG.  4 provides an illustration of the two specific conditions occurring in NLOS environments.  FIG.  4(a) shows “light” NLOS/NLOS-DP in which the direct path is attenuated but can be detected.  FIG.  4(b) shows severe NLOS/NLOS-NDP in which the direct path is not detected at all. 
In NLOS the estimated distance can be modeled as{circumflex over (d)}NLOS=d+bNLOS+nNLOS  (19a)bNLOS=bm+bpd+bB  (19b)where bm is the bias induced by the multipath, bpd is the bias caused by the propagation delay through obstacles other than free-space and bB is a deterministic but spatially random (due to the unknown nature of the obstacle) additive bias that models the amount of bias induced due to the obstruction of the DP. 
It is reasonable to model the range measurement biases similar to those in LOS but with some variation on the parameters.  For example if we would assume that NLOS biases follow an exponential distribution then it can be valid given that there is a minimum bias which shifts the exponential distribution.  Thus a valid model in this case is given byfb(b|NLOS)=bmin+λNLOSe−λNLOSb  (20)where bmin>0 and λNLOS<λLOS.  The rate parameter λ will control the magnitude of ranging biases and subsequently the probability of occurrence.  bmin can control the minimum bias that can occur where this addition is very useful as it can relate the range biases to a specific scenario. Tracking Performance of EKF with Biased Range Measurements
In order to appreciate the impact of the biases on the tracking performance the present inventors tested the EKF under LOS and NLOS using different models, where the modeling parameters have been adopted from experimental findings reported in [15].  FIG.  5 illustrates the Cumulative Distribution Function (CDF) plots of the biases in LOS, NLOS and severe NLOS for two system bandwidths (WiFi (FIG.  5(a)) and UWB (FIG.  5(b)). 
The exponential modeling parameters of the CDF plots are summarized in Table 1.  Note that severe NLOS refers to NLOS but with a higher bmin; indicating that the MD is constantly ranging behind “heavy” obstructions. 
TABLE 1Range bias exponential modeling parameters for 2 system bandwidthsSystem BandwidthChannel Conditionbmin (m)λWiFi - 20 MHzLOS1⅓NLOS1⅛Severe NLOS5⅛UWB - 3 GHzLOS02NLOS0⅓Severe NLOS2⅓
Next the impact of the biases in each condition/system on the EKF tracking performance is illustrated. 
FIG.  6 illustrates the simulated tracking results for the different channel conditions using a WiFi system with 20 MHz bandwidth.  In the simulation 4 BSs 1 are placed in an indoor environment surrounding the MD track 3 (actual track).  FIG.  6(a) shows the tracking (predicted track 4) in LOS conditions, FIG.  6(b) shows the tracking the NLOS condition and FIG.  6(c) shows the tracking in severe NLOS conditions.  For each condition, all range measurements from the MD to the BSs 1 experience error from the respective statistical distribution in FIG.  5(a) and Table 1. 
For WiFi, range measurements in LOS experience large mulitpath biases because of the low time-resolution (low system bandwidth), see FIG.  5(a).  The tracking performance of WiFi in LOS is impacted significantly.  For NLOS condition, the biases are so large that the estimated track 4 is significantly biased.  This indicates that if the MD receives range information from the BS biased around 5-15 meters with high probability, then it is extremely difficult to estimate and track the position accurately. 
FIG.  7 illustrates the simulated tracking results for different conditions using a UWB system with 3 GHz bandwidth.  As with FIG.  6, FIG.  7(a) is LOS, FIG.  7(b) is NLOS and FIG.  7(c) is severe NLOS.  For each condition, all range measurements from the MD to the BSs 1 experience error from the respective statistical distribution in FIG.  5(b) and Table 1. 
The performance of UWB systems in LOS improves significantly due to the high time-resolution capabilities thus combating multipath; which is a well known attribute of the technology [15]: see FIG.  7(a).  Thus the predicted path 4 and the actual path 3 of the MD are almost identical and can hardly be distinguished. 
The UWB performance in NLOS conditions degrades due to higher biases, but is better than the localization performance in WiFi systems (compare FIGS.  7(b) & (c) with FIGS.  6(b) and (c)).  This is mainly due to the fact that UWB has higher time resolution enabling more accurate estimation of first arrival of non-direct path components. 
The present invention aims to provide methods and systems which provide improved bias estimation, tracking and range/distance measurement correction as well as improved localization, particularly in NLOS/severe NLOS conditions.