4.1 Particle filter

Particle filters are a promising method for WiFi RTT positioning, as they are better suited than standard Kalman filters for highly nonlinear problems (Arulampalam et al., Reference Arulampalam, Maskell, Gordon and Clapp2002; Gentner et al., Reference Gentner, Ulmschneider, Kuehner and Dammann2020). WiFi RTT is a nonlinear problem as the RTT measurements are not linear functions of the mobile device position and are non-Gaussian. The latter is mainly a result of NLOS reception errors and multipath effects.

The process for the Particle and Genetic filters follows the diagram shown in Figure 1.

Figure 1. Generic particle and genetic filter process.

An initial position estimate is determined based on a single epoch least squares algorithm.

Initialisation – N _p particles are created using the initial position estimate as the mean. The particles are then randomly distributed around the initial position estimate to account for uncertainty. The standard deviation of the state distribution matches the uncertainty of the initial position estimate. The initial distribution of the particles follows a Gaussian distribution.

Prediction – For positioning in motion, the PDR model described in Section 3 is used. The positions of the particles are adjusted using

(3) $$\;\;X_{k + 1}^n = X_k^n\; + \;cos\left( {\varphi _{k + 1}^n} \right)\;r_{k + 1}^n\;\;\;\;$$

(4) $$\;\;Y_{k + 1}^n = Y_k^n\; + \;sin\left( {\varphi _{k + 1}^n} \right)\;r_{k + 1}^n\;\;\;\;$$

where X are the x coordinates of the particles and Y are the y coordinates of the particles, according to the coordinate frame of the environment, which can vary for indoor environments, k represents the epoch of the particles and n represents the particle index. The displacement, r _k , and orientation, ${\varphi _k}$ , are retrieved as described below and defined in Equations (5) and (6), respectively.

As steps are discrete events, the movement of the device appears to be discrete from the perspective of the filter as the prediction step will update the device’s position at each step. This is an issue as WiFi measurements can be recorded mid-step, which are used to determine the weights of particles. Thus, only updating the particle states per step will not give a precise representation of the device’s state. To resolve this, step-lagged smoothing is used. This method allows the filter to account for the movement of the mobile device during a step by using the time taken between a step along with the step length to determine a velocity of the step. When a new epoch of WiFi RTT ranges is received, the distance travelled is determined using the step velocity and the time taken between the WiFi RTT epochs, allowing the particles access to a more precise representation of the mobile device’s expected position during a step.

The step length and orientations are retrieved as described in Section 3 and then system noise is incorporated into these values, as described as follows for the displacement and heading respectively:

(5) $$r_{k + 1}^n\; = {{N\left( {{q_{k + 1}},\;{\sigma _q}} \right)} \over {{n_e}}}$$

(6) $$\varphi _{k + 1}^n = \;N\left( {{\varphi _{k + 1}},\;\;{\sigma _\varphi }} \right)\;\;$$

where r _k represents the displacement of the device for an epoch, q represents the step length, ${\sigma _q}$ represents the assumed standard deviation of the step length error, ${n_e}$ represents the number of filter epochs that have occurred between a full step and N indicates that each element is sampled from a Gaussian distribution. In Equation 6, ${\varphi _k}$ represents the heading of the mobile device, $\varphi _{k + 1}^n$ represents the heading of the device with Gaussian noise incorporated and ${\sigma _\varphi }\;$ represents the assumed standard deviation of the heading error.

The full motion model algorithm is shown in Figure 2 and occurs during the prediction step of all filters.

Figure 2. Step-lagged PDR motion model.

Update – The weights of each particle are updated by calculating how closely the measured ranges to each landmark match the particles’ distance to the landmarks. The model used for weighting is described as follows:

(7) $${w_{k + 1}} = \;{w_k}\;{1 \over {\sigma _k^i\sqrt {2\pi } }}{e^{ - {1 \over 2}{{( {{{z_k^i - d_k^i} \over {\sigma _k^i}}} )^2}}}}$$

The update step for computing weightings is the step for determining how strongly a particle matches the measured state from state estimates. The Euclidean distance, $d_k^i$ , between each particle and the landmarks are computed, where i is the AP landmark being measured from and k represents the epoch which is then treated as the mean in a Normal distribution alongside a standard deviation. The standard deviation is assumed to be equal across all APs at this stage, which is modified using RSSI-based outlier detection which is described in Section 5. Once the normal distribution is determined, the probability distribution function (PDF) of the distribution at $z_k^i$ , the measurement obtained for the distance between the AP, i, and the mobile device is calculated at epoch k, which gives the particle weight for that landmark. The weights for all landmarks for each epoch are then multiplied together to give a final weight for that particle and this process is repeated for all particles. Following this, the weights are normalised.

Resampling – If Equation (8) evaluates to be true, then Sequential Importance Resampling (SIR) is carried out. In this process, the higher weighted particles are retrieved and replace lower weighted particles whilst resetting all weights.

The effective sample size, N _eff , is determined as shown in Equation (9) so long as Equation (8) is evaluated as true,

(8) $${N_{eff}} \lt {{{N_p}} \over 2}$$

(9) $${N_{eff}} = {1 \over {\sum\nolimits_{n = 1}^{{N_p}} {{{\left( {w_k^n} \right)}^2}} }}$$

where N _p is the number of particles and $w_k^n$ is the normalised weight of the nth particle at epoch k. In most Particle filters, SIR is used (Arulampalam et al., Reference Arulampalam, Maskell, Gordon and Clapp2002), where resampling is triggered because the number of effective particles is too low. SIR essentially replicates higher weighted particles and deletes lower weighted particles. This is done to mitigate particle degeneracy. As the algorithm progresses, there will be particles that have very low weight and contribute nothing to the solution, and are thus wasted computation. The resampling step essentially removes these low weighted particles in favour of the higher weighted particles such that more particles are being used to identify the position of the mobile device.

Estimate – Compute the mean and standard deviation of the states using the particle weights according to the following:

(10) $$\bar x_w^k = \;{{\sum\nolimits_{n = 1}^{{N_p}} {w_k^n} X_k^n} \over {\sum\nolimits_{n = 1}^{{N_p}} {w_k^n} }}$$

(11) $$\bar y_w^k = \;{{\sum\nolimits_{n = 1}^{{N_p}} {w_k^n} Y_k^n} \over {\sum\nolimits_{n = 1}^{{N_p}} {w_k^n} }}$$

where $\bar x_w^k$ is the weighted average of the x coordinate of the position estimate at epoch k and $\bar y_w^k$ is the weighted average of the y coordinate of the position estimate at epoch k.

The Prediction, Update, Resampling and Estimation steps are repeated for each epoch.

4.2 Genetic filter

The Genetic filter (Higuchi, Reference Higuchi1997; Park et al., Reference Park2007) offers an alternative to SIR at the resampling step; instead, the resampling follows a genetic algorithm. The intention of this method is to provide an alternative and potentially more effective way of mitigating particle degeneracy and increase the diversity of the particles during resampling than SIR, such that fewer particles could be used for the same level of accuracy. Mitigating particle degeneracy is important as it reduces the chance of the particles of a Particle filter focusing on false positives too quickly. Essentially, the Genetic filter provides the particles with more opportunity to explore the search space to identify regions which could result in a higher weighting. A Genetic filter has been applied to an ultrawideband (UWB) positioning solution by Zhou et al. (Reference Zhou, Lau, Bai and Moore2021).

The resampling step of a Genetic filter is composed of three steps: classification, crossover and mutation. At the selection/classification step, the algorithm gets a set of strongly weighted and weakly weighted particles. Classification is the process of sorting the particles into mating pools to be used in the crossover step. The particles are sorted into ascending order based on weight. Then, the particles are split into two particle sets, computed by

(12) $$nint\left( p \right)\; \le \;{N_{eff}} \lt nint\left( {p + 1} \right)$$

where $p$ represents the integer that is used as the index to split the ordered particle set, where the set with index 0 to p is the higher weighted set and the set with starting index p+1 to N _p is the lower weighted set.

For the crossover step, to summarise, each strongly weighted particle undergoes an arithmetic crossover with another strongly weighted particle to produce two offspring particles which will replace the two parent particles. Then, all weakly weighted particles undergo an arithmetic crossover with a strongly weighted parent particle and the new particles replace the original particles.

Crossover is the process of taking two random particles from the mating pools and exchanging information between them to construct a new particle (referred to as offspring). The crossover operation conducted for this algorithm is an arithmetic crossover,

(13) $${\rm{\;}}\matrix{ {\omega _k^{off,{\rm{\;}}n} = {\rm{\;}}{\alpha _1}\omega _k^{h,n} + \left( {1 - {\alpha _1}} \right)\omega _k^{h,n + 1}} \cr {{\\{\omega _k^{off,{\rm{\;}}n + 1} = {\rm{\;}}{\alpha _2}\omega _k^{h,n + 1} + \left( {1 - {\alpha _2}} \right)\omega _k^{h,n}}}} \cr }\qquad \qquad {\rm{where}}\,\,n \le p$$

(14) $${\rm{\;}}\matrix{ {{\alpha _1} = {{w_k^{h,n}} \over {\left( {w_k^{h,n} + {\rm{\;}}w_k^{h,n + 1}} \right)}}} \cr {{\alpha _2} = {{w_k^{h,n + 1}} \over {\left( {w_k^{h,n} + {\rm{\;}}w_k^{h,{\rm{\;}}n + 1}} \right)}}} \cr }\qquad \qquad {\rm{where}}\,\,n \le p$$

Equations (13) and (14) describe the arithmetic crossover conducted on the higher weighted set. The two offspring particles are then used for the next epoch. The total number of particles remains the same. Here, $\omega _k^{off}$ represents the states of the offspring particles and $\omega _k^h$ represents the states of the higher weighted particles,

(15) $$\omega _k^{off,n} = {\rm{\;}}\beta \omega _k^{l,{\rm{}}n} + \left( {1 - \beta } \right)\omega _k^{h,0}\qquad \qquad {\rm where\,\,{\it p}\, \le {\it n} + {\rm{ }}1 \le {{\it N}_p}$$

Equation 15 describes the arithmetic crossover conducted for the lower weighted particle set. In this version of the crossover, instead of taking two parents from the same set, a particle from the higher weighted set is crossed over with a random particle from the lower weighted particle set. The parent particle used and the new offspring particle are then passed onto the next generation. Here, $\omega _k^l$ represents the states of the lower weighted particles, where β is a random number between 0 and (N _p − N _eff )/N _p.

For this algorithm, a mutation step was not implemented, as noise is incorporated in the prediction step which achieves the same objective.

4.3 Grid filter

The Grid filter differs from the Genetic filter and Particle filter in the initialisation, prediction and resampling steps. Instead of using particles, the search space is split into grid squares. Each grid square intersection represents a candidate position of the mobile device and fundamentally, instead of particles moving around and each particle’s weight being updated, weights are moved throughout the grid to represent the posterior distribution of the device’s state. This algorithm has also been explored by Zhong & Groves (Reference Zhong and Groves2022). The process follows the following steps, this is also shown in Figure 3.

Figure 3. Grid filter process.

Initialisation – The search area is initialised as a square grid of a given size $\kappa $ (in mm) centred around the initial position estimate, this is then split into two-dimensional (2D) grid squares with grid spacing, $\eta $ (in mm), to form an array of dimensions, $\kappa /\eta $ . Following this, the grid squares are weighted according to their distance to the centre of the grid (the initial position estimate) according to a Gaussian distribution with a standard deviation that matches the sensor noise. For this paper, $\eta $ was set as 50 mm and $\kappa $ was initialised as 2 m.

Prediction – For positioning in motion, the PDR model described in Section 3 is used to retrieve the displacement and heading values. However, the position displacement differs from that of the Particle filters to ensure that the coordinates remain in a grid formation. The displacement of all grid coordinates is determined as described by

$$\;X_{k + 1}^n = X_k^n\; + \;cos\left( {{{\bar \varphi }_k}} \right){\bar r_k}$$

(16) $$\;Y_{k + 1}^n = Y_k^n\; + \;sin\left( {{{\bar \varphi }_k}} \right){\bar r_k}\;\;\;$$

where the mean values of the heading, ${\bar \varphi _k}$ , and displacement, ${\bar r_k}$ , are used to displace the grid.

However, this does not effectively account for sensor noise. Thus, a separate process must be followed to incorporate sensor noise into the grid coordinates and weights, which is described as follows. First, a new set of grid coordinates ‘c grid’ is generated that are equal to the grid coordinates of the main grid. A multivariate Gaussian distribution is generated from the heading uncertainty and step length uncertainty. From this distribution, a new set of weights are generated for the ‘c grid’ by calculating the probability density function of the distribution for the distances between each coordinate and the mean of the coordinates. With the weights of the ‘c grid’ and the weights of the grid coordinates from the previous epoch, the weights are convolved using direct convolution to form a new set of weights for the grid coordinates that incorporate the sensor noise of the PDR.

Weighting – Each grid point is weighted according to its distance from all of the landmarks. The model used for the weighting is the same as the Particle filter and Genetic filter, described in Equation (7). This weighting represents the likelihood that the device is at that grid point. Following this, the weights are normalised.

Using the weights of each grid point, the weighted average position is determined according to Equations (10) and (11). This represents the estimated position of the mobile device at that epoch.

Repeat the Prediction, Weighting and Estimation steps for each epoch.

6.1 Methodology

The experimental layout for static positioning is described in detail by Raja and Groves (Reference Raja and Groves2023). The experiments to test the algorithms in motion involved moving through different routes in a single environment. The environment was irregular in shape and had many walls, both interior and exterior, to create a diverse range of multipath and NLOS scenarios during the routes walked by the pedestrian. The algorithms were tested across three different routes, varying in complexity and NLOS conditions; these are shown in Figure 4. Each route was repeated in both directions, so in total, six scenarios are collected. In Figure 4, the ‘forward’ route starting point is shown with the purple square and the ‘reverse’ route starting point is shown with the red square. The trials involved a pedestrian holding a smartphone and walking on top of fixed step markers placed on the ground along the route of the trial. The markers were approximately 670 mm apart. The location of each marker was measured against the reference origin to generate the step marker’s ground truth coordinates. To align the ground truth data with the measured data, the trials were filmed. The timestamps of the steps from the video were used to determine the expected position of a device at a given time, allowing the algorithms to have an accurate comparison point at each step. The standard deviation of the step length assuming white noise was calibrated as follows. A pedestrian carrying the mobile device walks 20 steps in a straight line 10 times. The real distance travelled is measured and the estimated distance travelled is calculated using the PDR model. The standard deviation of the difference is computed and this is divided by the number of steps to give the standard deviation of a step.

Figure 4. Positioning in motion experimental environments.

6.2 Experimental results

6.2.1 Static positioning

The preliminary results of the static data were presented at the ION GNSS+ 2023 conference by Raja and Groves (Reference Raja and Groves2023), however, the Grid filter algorithms have since been improved and will be discussed here.

The RMSE of the positioning solution for each environment is shown in Table 1. Each row refers to the algorithm being used with and without outlier detection, whilst each column refers to the trial. Each trial is also marked if it had NLOS scenarios. Here, 27 out of 42 or 64.3% of tests produced sub-metre accuracy and 41 out of 42 or 97.6% of trials had an RMSE below 2 metres, indicating there is a strong argument for WiFi RTT being able to produce sub-metre accuracy for positioning, as pitched by Google (2022), as long as the AP biases are calibrated. Environments B, C and D where all APs had a LOS to the mobile device produced sub-metre accuracy for 20 out of 21 trials, indicating that in optimal conditions (with calibration), WiFi RTT could provide a reasonable positioning solution for most pedestrian navigation use cases. In environments A, E and F where there were NLOS signals present, only 7 out of 21 trials achieved sub-metre accuracy. This is to be expected as NLOS and multipath effects vary substantially from environment to environment.

Table 1. Positioning solution RMSE for trials and algorithms configuration in static

Table 2 shows the percentage improvement of each algorithm combination against the baseline least squares algorithm, the mean percentage improvement of each algorithm across all trials is shown on the rightmost column. In environments C and D, all of the filters improved the positioning accuracy over least squares. This is also true for environments A and B. However, in environment F, the Particle filter and Grid filter without outlier detection performs worse than least squares. The Grid filter with outlier detection produced the highest improvement overall of 80% in environment A. The algorithms provided the greatest mean improvement in environment B; this is the simplest environment with the smallest distances between the APs and mobile devices, and no NLOS signals, and would be expected to have the highest accuracy given there are fewer error sources when compared with more complex environments like environments E and F. The best performing algorithm on average across all environments was the Genetic filter with and without outlier detection with 48% and 36% mean percentage improvement over single epoch least squares, respectively.

Table 2. Percentage decrease of RMSE against least squares for each environment in the static case and algorithm configuration

6.2.2 Positioning in motion

In this section, the algorithms will be analysed for all scenarios where the mobile device was in motion. Each scenario going forward and in reverse will be analysed separately. Table 3 summarises the results of the trials such as the minimum positioning solution RMSE at any step in the trial, the maximum positioning solution RMSE at any step in the trial, the average positioning solution RMSE of all steps in the trial and the final positioning estimate RMSE of a trial. The average column is of particular note.

Table 3. RMSE position error statistics for trials in motion and each algorithm combination

As shown in Table 3, the least squares solution positioning accuracy is poorer than the other algorithms. The average positioning RMSE for least squares was 2.79 m when taking an average across all trials. For all trials except Trial 2R, least squares was the worst performing algorithm based on the average RMSE across all steps. Furthermore, for all trials except Trials 2R and 3R, the final position RMSE was highest for least squares. Observing the least squares position estimates in the maps shown in Figures 5, 6, 7, 8, 9 and 10, it is clear that irrespective of the results in Table 3, the position estimates do not correctly follow the pedestrian’s path and are visibly inaccurate, placing the pedestrian in an incorrect room on most occasions, which demonstrates the importance of the PDR model in the positioning algorithm. The following analysis will exclude least squares due to the reasons described above.

Figure 5. Trial 1F position per epoch.

Figure 6. Trial 1R position per epoch.

Figure 7. Trial 2F position per epoch.

Figure 8. Trial 2R position per epoch.

Figure 9. Trial 3F position per epoch.

Figure 10. Trial 3R position per epoch.

In 8 out of 36 trials with filtered solutions, the average position RMSE over all steps of a given trial was sub-metre and in 29 out of 36 trials, the average position RMSE was below two metres. In virtually all trials, the PF and GrF performed similarly to each other, and the GrF performed worse on average in comparison to the other algorithms when looking at average RMSE. All algorithms in Trial 1R had a positioning RMSE greater than 2 m. For Trial 1R, as can be seen in Figure 6, the paths of all filters are translated approximately 2.5 m to the north of the actual path. This positioning error is worse than Trial 1F for all algorithms. The poor accuracy of the algorithms in comparison could be a result of poor initialisation which could be a result of fewer LOS APs during the earlier part of the trial. This discrepancy is potentially caused by the geometry of the APs, as shown in Figure 6, in combination with the pedestrian’s path. Practically, these kinds of poor initialisation errors can be resolved with better PDR algorithms that account for error in the orientation sensor as well as access to more user-specific calibration data for step length. In this trial, the pedestrian blocks AP1 and AP2 with their body. This results in attenuation and, since the range is lower than 4 m, the outlier detection algorithm will not deweight the range estimate and without outlier detection, the algorithms will treat the ranges as reliable, even though they are NLOS due to the pedestrian’s body blocking the signals. This threshold could be calibrated and varied to account for these scenarios, but this has not been explored in this paper.

As shown in Table 4, the initial position RMSE was 3.11 m, which is approximately 5.4 times the initial position RMSE of Trial 1F. This will drastically reduce the accuracy of the results from all algorithms. The Grid filter in Trial 2R also had an average position RMSE greater than 2 m. In Figure 8, this is evident as, compared with the other algorithms, the path of the pedestrian is estimated at a different location. In Table 4, the initial position RMSE was greater than 3 m, contributing to the overall positioning inaccuracy of the trial. For all algorithms, the larger positioning error at the final steps of Trial 2R is due to a poor orientation value at step 7, resulting in an incorrect trajectory of the pedestrian, demonstrating a shortfall in the orientation feature of the PDR model. Poor orientation readings also impacted the performance of the algorithms in 3F and 3R, demonstrating that improving the orientation sensor or orientation processing model could improve the overall positioning accuracy of the algorithms. These errors do not accumulate linearly with step count. The system is mostly impacted by anomalies in the readings, be this from the step length being incorrect or the orientation being incorrect for a given prediction step.

Table 4. Initial position RMSE for each trial in motion

In 33 out of 36 trials, the final position estimate was more accurate than the initial position estimate. This indicates that not only were the algorithms successful at maintaining the path of the pedestrian, they also improved the accuracy of the position estimate compared with the initial position estimate. With outlier detection included, the greatest absolute improvement was made by the PF in Trial 2F with an improvement of 1.32 m and the greatest improvement without outlier detection was made by the PF in Trial 2R with an improvement of 1.23 m. The greatest percentage improvement between the initial and final position estimates without outlier detection was made by the GrF in Trial 2F with an improvement of 287%.

Trials 3F and 3R were the longest trials with a combination of different manoeuvres by the pedestrian. These trials consisted of 35 steps with the pedestrian also walking outside of the building. Figure 9 shows a trial run on the map when the device is moving through the forward scenario (Trial 3F). Figure 10 shows a trial run on the map when the device is moving through the reverse scenario (Trial 3R), beginning outside of the building. In the case of Trial 3F, the GrF final position RMSE is almost double that of the other filters with an RMSE of 2.51 and 2.37 m for the GrF without and with RSSI-based outlier detection, respectively, and from Figure 9, the path is also visibly incorrect, unlike the other filters. However, in Trial 3R, it appears that the GrF better guides the PDR model through the path and most accurately represents the path of the user. The final position estimate for the GrF in this trial was sub-metre both with and without outlier detection.

Table 5 shows the RMSE percentage improvement of the average positioning error for each trial for the algorithms with and without RSSI-based outlier detection. As can be seen in 15 out of 18 trials, the RSSI-based outlier detection improved the positioning solution. Any outlier detection will have a false positive rate, so these results are quite promising as the decrease in accuracy from the expected false positive rate does not result in a drastic change in accuracy. The largest improvement was with the GrF in Trial 2F, where the outlier detection provided an improvement of 25%. The outlier detection provided the greatest improvement for the Particle filter with an average improvement of 14%, whilst the outlier detection provided the least improvement to the GrF with an average of 4%. A large limitation of the RSSI-based outlier detection algorithm is caused by an edge case where the increased WiFi-RTT measurement matches with the RSSI path loss propagation model measured range. This results in the reading not being picked up by the RSSI-based outlier detection model. As a result, the erroneous signal will reduce the accuracy and reliability of the overall positioning solution. Furthermore, the material of the blockage (brick wall, body, wooden door etc.) varying also has an impact on the RSSI property of the signal compared with the WiFi RTT measured range. This property will depend on the attenuating properties of the material.

Table 5. RMSE percentage position error improvement for algorithms and trials in motion