Kalman filter with a linear state model for PDR+WLAN positioning and its application to assisting a particle filter

Indoor positioning based on wireless local area network (WLAN) signals is often enhanced using pedestrian dead reckoning (PDR) based on an inertial measurement unit. The state evolution model in PDR is usually nonlinear. We present a new linear state evolution model for PDR. In simulated-data and real-data tests of tightly coupled WLAN-PDR positioning, the positioning accuracy with this linear model is better than with the traditional models when the initial heading is not known, which is a common situation. The proposed method is computationally light and is also suitable for smoothing. Furthermore, we present modifications to WLAN positioning based on Gaussian coverage areas and show how a Kalman filter using the proposed model can be used for integrity monitoring and (re)initialization of a particle filter.


Introduction
Wireless local area network (WLAN) access points (APs) are numerous and ubiquitous in most indoor environments. Although WLAN is meant for data transfer, the WLAN signals can also be used for user localization. Because the WLAN APs are not meant for positioning, they do not usually send information on their own location to clients. WLAN positioning therefore makes use of a 'radio map' that describes certain features of the WLAN signals at given locations. A radio map is created and updated using data collected on-site called fingerprints (FPs). A FP is a report that contains at least the receiver location and the IDs and received signal strength (RSS) values of APs within reception range. A radio map is constructed off-line from the collected FPs. The accuracy of positioning based on WLAN signals depends on the model and environment. In small scale (a few buildings), it is possible to achieve positioning results of the order of a couple of meters [1]. For large-scale positioning (a city or larger), when the size of the database is a limiting *Correspondence: matti.raitoharju@tut.fi Department of Automation Science and Engineering, Tampere University of Technology, P.O.Box 692, FI-33101, Tampere, Finland factor, the positioning accuracy can be of the order of tens of meters [2].
Pedestrian dead reckoning (PDR) uses an inertial measurement unit (IMU) to detect when a user takes footsteps and how the direction changes between footsteps. The IMU has three-axis accelerometers and gyroscopes. The user heading change is computed by projecting the gyroscope measurements to the horizontal plane, which is estimated from the accelerometer [3]. It is also possible to use a compass for heading determination [4].
The footstep length may also be estimated from the IMU data. If the sensor is mounted on the foot, it is possible to detect when the foot is motionless and then integrate the footstep length from the IMU data [5]. If the IMU is handheld, this method cannot be used, but the footstep length can be inferred from the IMU data also by other methods, see for example [6,7]. A PDR system can greatly improve the positioning locally as the position estimate may be updated at every footstep, but because the errors accumulate over time, PDR is often combined with other sensors that can, at least occasionally, provide absolute position information.
In this paper, we investigate models for fusing heading change information from IMU with WLAN measurements. We propose a linear state model for the state evolution, whereas in the literature, the state model used with PDR system is usually nonlinear [7][8][9][10] or the heading information is used as nonlinear measurements [4]. In general, there is no closed form optimal algorithm for nonlinear estimation. In this paper, we use the Kalman filter (KF), the extended Kalman filter (EKF), and the unscented Kalman filter (UKF), which are computationally light algorithms, for evaluation of the models.
Particle filters (PFs) are algorithms that can be used for nonlinear estimation. We show how the proposed model with a linear KF can be used for initialization and integrity monitoring of a PF. Section 2 contains the filtering and smoothing algorithms that are used to estimate the user's kinematic state. In Section 3, we present the WLAN model that is used for positioning. The evaluated PDR models are presented in Section 4. In Section 5, we discuss the possibilities of using the proposed model in a particle filter, and in Section 6, we evaluate the performance of different models with real and simulated data. Section 7 concludes the paper.
This paper is an extended version of a paper [11] presented at the Conference on Design and Architectures for Signal and Image Processing (DASIP 2013). In general, the notations in this paper are changed to be more consistent. We added UKF to Section 2. Section 3 has a new subsection containing modifications to the WLAN coverage area positioning. The first model presented in Section 4 has cross terms in the covariance matrix of the state transition matrix to make it mathematically sound. Section 5 is new. In Section 6, tests are new or done with more test tracks and the UKF is also used.

Filtering algorithms
The KF is an algorithm for estimating the state of the system given a time series of measurements in the case of linear state evolution and measurement models. If the measurement and state transition errors are independent, white, and Gaussian processes, the KF computes the posterior mean that is the optimal estimator in minimum mean square error sense (i.e., the KF is the exact Bayesian filter) [12].
The state model is as follows: where x t is the state at the tth time instant, F t is the state transition matrix, and Q t is a zero mean Gaussian noise term with covariance matrix Q t . The KF prediction is computed from the previous posterior N(x t−1|t−1 , P t−1|t−1 ) via the update formulas: that give the new prior N(x t−1|t , P t−1|t ).
The prior is updated using the linear measurement model: where matrix H t defines the relationship between the state and the measurement vector y t , and R t is a zero mean Gaussian noise term with covariance matrix R t . In the update part of the KF, the prior is updated via: whereȳ t is the realized measurement value, S t is the innovation covariance matrix, K t is the Kalman gain,x t|t is the posterior mean, and P t|t the posterior covariance matrix. If there are no measurements, the latest prior is also the new posterior, and if there are several independent measurements in one time instance, the updates can be applied consecutively.
In more general situations, the state transition model is nonlinear, of the form: where f t is a nonlinear state transition function and ε Q t does not necessarily have zero mean. In the EKF [13], the mean propagation formula (Equation 2) has to be replaced with: where ε Q t is the mean of the noise term. The covariance propagation formula (Equation 3) changes to: wherê Similar approximations can also be done for the measurement update equations, but in this paper, we use only linear measurements.
The EKF requires analytical computation of the partial derivatives and approximates the measurement function in a point. One commonly used algorithm for computing the update that does not need analytical differentiation and approximates the state transition function in larger area is the UKF [14]. In the UKF, a set of sigma points is chosen so that they have the same mean and covariance as the original distribution. The sigma points are propagated through nonlinear functions, and the estimate is updated using the transformed sigma points. The transformation is called unscented transformation. The UKF can be used also for measurement update, but because we use only linear measurements, we present only the state propagation of the UKF.
When computing the state propagation, the UKF uses augmented state: A popular choice of sigma points is as follows: where d is the dimension of the augmented state, ξ UKF is an algorithm parameter, and √ P is a matrix square root for which: This matrix square root can be computed, for example, with Cholesky decomposition. Parameter ξ UKF is defined as: where parameters α UKF and κ UKF define how much the sigma points are spread. The prior mean and covariance computed using sigma points are as follows: where the sigma points are computed from the posterior of the previous time step using Equation 15 and the weights are as follows: The parameter β UKF is related to the distribution of the state. In the case of Gaussian distribution, β UKF = 2 is optimal [15].
The Rauch-Tung-Striebel smoother [16] can be used to enhance the state estimates when measurements of 'future' time instants are available, for example, when plotting the track over a given time interval. The recursive smoothing equations are as follows: for t = m−1, m−2, . . . , 1, where m is the final time index.

Coverage area positioning
In its simplest form, probabilistic coverage area (CA) positioning is a method for radio map construction in which the reception area of each AP is modeled as a bivariate normal distribution. Because the radio map does not contain any raw RSS data, computational, memory, and communication complexity are much lower compared to those of conventional FP positioning methods [17]. The algorithm is explained in [18] and here we only briefly outline it.
The coverage area estimate is as follows: where μ i,n i and i,n i are the mean and covariance of the CA estimate of ith AP based on n i FPs, z i,j is the location of the jth FP where AP i was received and B is the prior covariance. When a new measurement is received, these parameters can be updated by: In positioning, the CA model from every received AP could be applied in separate KF updates. To reduce computational cost, it is equivalent to first compute the static position estimate from these CA models using: where i ∈ FP is used to indicate AP:s received in the current FP and the second index n i is not shown because in positioning, the latest CA estimate is always used. When using CA models in KF update (Equations 5 to 8), the linear and Gaussian measurement model is as follows: where ε R t ∼ N(0,ˆ ). Here, it is assumed that the position variables are the first components of the state.
In [19], it was shown that using separate CA models for strong RSS values improves the positioning accuracy. That is, in addition to CA models generated using all FPs, one generates CA models using only FPs with strong RSS values. In the positioning stage, the models are chosen according to signal strengths. In [17], the two-level CA positioning was compared with other WLAN positioning methods. Of the compared methods, there was no significantly faster or more accurate method than the two-level CA positioning.

Modifications to CA positioning
In this section, we present three different heuristics to mitigate some shortcomings of CA positioning. In order to preserve the assumption of independent measurements, the proposed methods do not use the prior information provided by the Kalman filter. The heuristics are applied to measurements in the order they are presented.

Minimum size
When there are many FPs, the effect of the prior B used in Equation 27 diminishes. As a result, if many FPs are collected in a small area, the covariance matrix becomes too small. In our experience, this is a problem especially with CAs that model weak RSS in multilevel CA models. In [11], the prior B for weak CAs was set to a very high value 1, 000 2 m 2 I to compensate. Here, we propose to impose a lower limit on the variance of the Gaussian CA model.
The variances along the main axes can be computed using eigendecomposition: Diagonal matrix consists of the variances corresponding to the main axes (columns of V ) on its diagonal. We propose to impose a minimum variance σ 2 min for main axes by thresholding:

Outlier removal
Because WLAN APs are freely deployable, the APs may be moved. If the database is not updated continually, the CA models of moved APs provide false information. For this kind of situations, we use outlier removal that is based on squared Mahalanobis distance to detect outliers.
In the outlier removal process, the mean estimate (Equation 31) and each CA model are compared separately. If the mean is outside the 95% ellipse, we discard the CA model. This is done by checking if: where χ −2 2 is the inverse cumulative density function of the chi-squared distribution with 2 degrees of freedom.
We propose that the CAs are removed one by one, removing the CA with largest d i until all are close enough to the mean. After removing a CA, the estimate (Equations 30 to 31) can be updated by: where j is the ID of the least probable CA model and variables with subscript + are the updated variables. If the reduction leads to a case where there are only two CAs left and they are distant to each other according to Equation 35, the WLAN measurements are not to be used at all.
There are examples in the statistics literature where the Mahalanobis distance does not work very well [20], but we have found it satisfactory in the tests reported here. For more advanced treatment, the effect of outliers can be mitigated by doing the positioning using distributions with heavier tails like Student's t-distribution. This could be done by adapting the algorithm presented in [21] to CA positioning.

MIMO compensation
Multiple-input multiple-output (MIMO) WLAN devices host several APs on one physical device. Such APs have naturally very similar CAs. If these CAs are used as independent measurements, the estimate covariance becomes too optimistic and the estimate mean gets biased towards the MIMO CA mean. The same problem arises also with non-MIMO APs that have similar coverage areas. This problem could be handled during the construction of the radio map by taking the correlations of received APs into account. This would however make the size of the radio map grow because there would be a need to store the cross correlation terms between every AP pair received simultaneously.
We propose the following heuristic method for compensating the overlapping MIMO CAs in estimation. The method is inspired by [22], where dissimilarity of two components of a Gaussian mixture is measured by: We use a slightly modified exponential form of Equation 38 given by: as the measure of dissimilarity between CAs. This measure is used to scale covariance of CAs such that if there is any number of identical CAs that are not similar with other CAs, they have the same effect as if there was only one CA, and that if two identical APs are further than 1 · σ from each other, they are considered independent. These properties are achieved by using MIMO-compensated CA variances MIMO in the position estimation (Equations 30 and 31), given by: Proofs that this has the properties stated above are given in the Appendix.
In Figure 1, the outlier detection followed by MIMO compensation is applied in a case of five APs received. One CA is an outlier, and there is a cluster of three similar CAs. The figure shows one-sigma ellipses of CAs and different user location estimates.

Pedestrian dead reckoning
We consider positioning systems in which PDR measurements are used in estimation in the state transition models and the WLAN measurements are then used in the update to fuse the measurements together in a Kalman filter. In most of the literature, the state model used in PDR is nonlinear [7][8][9][10]. For comparison purposes, we consider two traditional models. In the first one, the state contains the user location r and the direction of the movement θ: where θ t is the change of heading obtained from gyroscopes and s t is the footstep length estimated from accelerometer data. These are modeled as Gaussian random variables with means θ t and s t , which are the measured values, and variances σ 2 θ and σ 2 s . In EKF, the predicted mean is as follows: the linearized state transition matrix is as follows: and the linearized state transition noise covariance is as follows: where σ 2 is the variance of the variable in the subscript.
In the UKF prediction, the mean of the augmented state is as follows: and the covariance is as follows: In our second traditional model, the footstep length is also an estimated variable: where s is the estimated footstep length change, which is modeled as a zero mean Gaussian noise term with variance σ 2 s .
For this model, the EKF prior mean is as follows: the linearized state model is as follows: and because the state transition noises are additive, we have: In the UKF prediction, the mean of the augmented state is as follows: and the covariance is as follows: (52) Figure 2 shows position estimates after taking one footstep from a known position using this model. The figure shows how the EKF linearization does not take the 'curvature' of the true model into account and how UKF takes it into account by increasing the covariance.
We propose the linear state model: where v is the footstep vector estimate and v is the associated state transition noise. For this model, the predicted mean and covariance are computed using Equations 2 and 3. To make state transition function linear and noise additive and Gaussian, v is modeled as a Gaussian that is independent of state and measurement value θ t . The process noise covariance matrix is as follows: To compare the amount of state transition noise in traditional model (Equation 47) and linear model (Equation 53), we investigate the amount of noise propagated from step variables to position variables in one step. For simplicity, we use here independent step variables. In linear model, the covariance of the state variables increases every step by: and in the EKF, the position covariance is increased by: From this, we can see that when: the amount of propagated error is the same for both models. Because the noise of footstep variables is modeled additive and independent of the state, the increase of position variable variance caused by variance of one footstep is same using the different models when: Compared to the traditional models, the proposed model has the benefit that the Kalman filter can be used to compute the optimal estimate for this model when the measurements are linear and Gaussian. In this paper, we use only linear WLAN location measurements presented in the previous section and the information from PDR is used only in the state transition models. A drawback of this model is that it does not have different variances for the heading and footstep length. Constructing a noise model that has different variances for heading change and footstep length is possible by computing θ and s from step vector and adding noise similar to Equation 56 to the step vector. This would, however, make the noise nonlinearly dependent on the state estimate and linearizations would be required. In Section 6.1, we test the performance of linear model in a situation, where σ 2 θ s 2 = σ 2 s and show that the linear model can perform well also in that situation.
Another drawback is that the measured footstep length cannot be used in the linear model. If it was used as a measurement, the model would not be linear. The first presented nonlinear state model uses all the information, and the second presented nonlinear model uses the same amount of information as the proposed linear model. We will test also this in Section 6 and show that even though the proposed method cannot use all the available information, it can outperform the nonlinear model with footstep length when the initial heading is unknown.

Integration with a particle filter
Besides footstep length, other information that cannot be used with a linear KF includes the use of floorplan information or nonlinear measurement models in estimation. This information can be such that it cannot be used well in EKF, e.g., floorplan. PFs can be used in these situations [23][24][25]. Particle filters estimate the probability distribution using a set of weighted point masses called particles. The number of particles required for estimation is case dependent. The larger the uncertainties in the proposal distribution compared to measurement variance, the more particles are needed. Particle filters are usually computationally significantly more expensive than KFs, and so, running a KF alongside a PF does not significantly increase the computational complexity. In [25], the WLAN+PDR KF presented in [11] was compared with a PF that used wall information and a nonlinear state model. The PF with wall information reduced the mean error of KF from 4.8 to 2.0 m, while PF without wall information had slightly worse accuracy than KF. In the following, we present algorithms in a general form so they can be applied to many known PF variants [26].

Number of particles in initialization
In this section, we consider the situation where the prior covariance is so large that the PF would need an excessively large number of particles to adequately sample the region. We propose to use the KF until the estimate uncertainty is small enough for initialization of a PF that has a given number of particles.
We assume that the PF used is initialized by sampling from the prior distribution. Because the number of particles required for good estimation depends on many factors, there cannot be a simple algorithm for determining the number of particles. We assume that it is known how many particles n 1 is enough with a prior covariance P 1 and then compute the number of particles n 2 that is required with covariance matrix P 2 .
In general, a four-dimensional vector x is inside an ellipsoid that contains fraction p of the probability mass of distribution N(x 1 , P 1 ) when: where χ −2 4 is the inverse cumulative density function of the chi-squared distribution with 4 degrees of freedom. The volume of this four-dimensional ellipsoid is as follows: and the expected number of particles inside it is pn 1 , so the expected particle density is pn 1 V 1 . We propose to choose the number of particles n 2 in a such way that ellipsoid: has the same expected particle density as ellipsoid of Equation 59. Thus, we get rule:

Integrity monitoring
There is a risk that the PF estimate diverges far from the true location. Divergence may be caused by too few particles, error in the floorplan, etc. The linear estimate can be used to detect whether the PF has diverged using the algorithm proposed in [25], as follows. We assume that the particle filter operates in a state space that contains at least the same information as the linear KF. If the particle filter state has different parametrization, e.g., footstep length and heading instead of step vector, the particle states are converted to the same state space that the linear model uses. First, compute standardized deviances d i of all positively weighted particles: where c i is the location of ith particle and P and x are the parameters of the KF estimate. Then, check if the smallest deviance is below a set threshold. If it is not, then particles are reinitialized. Because also the KF estimate may diverge, the resampling can be done so that some of the original particles are retained and the rest are sampled from the KF estimate.
6 Performance evaluation

Synthetic data
In simulations, we tested how the linear model performs against the traditional models. If not otherwise specified, the true track is generated by sampling the initial state from prior and then sampling subsequent true states using the model given in Equation 47. In tests with simulated data, the WLAN measurements are independent and the improvements proposed in Section 3.2 are not needed.
Each simulated test track consists of 50 footsteps, and at every time step, there is a 10% probability of receiving a location measurement with variance 10 2 m 2 I. First, we test the effect of the initial heading variance on the positioning accuracy of the traditional models and comparing it with estimate obtained with linear model that does not use initial heading information. The traditional models require an initial mean for the heading variable, while the proposed linear model can be initialized with a zero mean step vector. It would also be possible to initialize the linear model with accurate heading and footstep length information, but when the initial heading variance is larger, it is not a straightforward task to select initial covariance for the linear model that is similar to the one used with nonlinear model, but that avoids the problems with linearizations. In the first simulation, the effect of the heading error of the initial mean is tested. The standard deviation of the initial heading is varied from 0°to 180°. The standard deviation of footstep length is set to σ s = 0.01m, and the standard deviation of heading change is σ θ = 0.01 0.7 rad. The initial footstep length is set to 0.7m for the true state and for the traditional model. The linear model is initialized to have a zero footstep vector with initial standard deviation σ v 0 = 1m. If the footstep length does not change much during the track, the linear model, where σ v = 0.01m, has the same amount of propagated error in position as the linearized traditional model. All models are initialized to have correct location mean and the initial variance for location dimensions is 10 2 m 2 . The standard deviation for initial footstep length is set to σ s 0 = 0.2 m. For the nonlinear model with given footstep length, we use accurate measurements and σ s = 0. The footstep length variable is not random in this situation, and in the UKF implementation, it is not used in the augmented state. Figure 3 shows the mean of the final position error of 500 simulated tracks for different standard deviations of the initial heading. The traditional models outperform the linear model when the standard deviation is small. This is expected because the traditional models use prior information of the footstep length and initial heading, unlike the linear model. The models with known footstep lengths are more accurate when the standard deviation is small. When the initial heading variance becomes larger, the linear model outperforms the traditional models. This is caused by the increase of the linearization errors in the traditional models.
In the second simulated test, the standard deviation of the initial heading is fixed to 90°, but the initial angle error of the traditional models is varied. Otherwise, the simulation parameters are the same as in the first test setup. Figure 4 shows how the traditional model is advantageous only if the heading variance is small; small realized error is not enough. This means that to gain the benefit of the known initial heading, the initial heading variance should also be set correctly. Here, the models that use known footstep length have worse accuracy than the models that estimate the footstep length.
Compared to simulation results presented in Figure 2 of [11], the main differences in this test are as follows: in this test, the footstep length variance is smaller, and thus, the traditional model does not end up into situation with negative footstep length that often. In [11], there was error in simulation variance: for location variables, it was 10m 2 instead of 10 2 m 2 . When using smaller initial variance, the accuracy of the linear model is on the same level as the accuracy of the traditional models. Next, we investigate how the different methods perform when the heading error (i.e., gyroscope accuracy) is changed. The traditional models use the correct σ θ and σ s = 0.01 m, and the linear model is tested with σ v = σ s and σ v = 0.7σ θ . In this test, we use initial state noise parameters σ r 0 = 0.2 m, σ θ = 0.2 0.7 rad, σ s 0 = 0.2 m, and σ v = 0.2 m. Figure 5 shows the results for this test. The linear model and traditional models that do not use the footstep length have similar accuracy when the heading change error is small, while the models that use footstep length are more accurate. When the noise level increases the linear model that uses the footstep length's variance as the σ v becomes less accurate than the other methods, but the linear model that uses step vector variance computed from the heading change variance with expected footstep length (0.7 m) according to Equation 58 is either on the same level or better than traditional models without known footstep length. Figure 3 of [11] shows a similar test's results. In [11], the initial variance is large compared to this article's test. Figure 6 shows how the 1 · σ ellipse of footstep noise is propagated to state noise in one footstep with different models, when the standard deviation for heading change noise is 20°. The linear model that uses variance of the footstep length has a small covariance compared to others, which causes high position errors. On other hand, the linear model with noise level computed from heading noise using Equation 58 has the largest noise covariance. It seems that this does not cause as much problems as having a too small error estimate. Similarly, when the footstep length is known, the EKF linearization does not take the curvature of the true distribution of the nonlinear model  into account and has bad accuracy, while the UKF has a better variance approximation, which results in better accuracy.

True data
In the true-data test, the algorithm is tested in two floors of a building in the Tampere University of Technology. The radio map of both floors contains approximately 300 WLAN APs and 5,500 FPs. The measurements were collected using a XSENS MTi IMU and Acer Iconia tablet. Both devices were carried in hand while doing measurements. Because XSENS was handheld, the errors in heading are caused, in addition to IMU noise, by the small changes in hand orientation compared to the movement direction. We used four test tracks. WLAN scans were done on average every 3.2 s.
For the WLAN positioning, we use the two-level coverage area method presented in Section 3 with all modifications from Section 3.2, if not otherwise stated. We generate two coverage areas: a weak CA that is constructed using all FPs and a strong CA that is constructed with only FPs that have signal strength ≥ −60 dBm. The prior for weak CA is B = 100 2 m 2 I and for strong CA B = 20 2 m 2 I. Minimum variances for main axes were set to 40 2 m 2 for weak CAs and to 5 2 m 2 for strong CAs. Most of the weak CAs had smaller variances than this limit, and they were extended using Equation 34. According to [19],   the two-level normal CA models have around 10% larger errors than traditional location fingerprinting, but require only storage of ten parameters for each AP. Figure 7 shows filtered and smoothed routes for a test track. The initial heading for traditional models is 90°e rror with variance of (90°) 2 , and the initial footstep length estimate is 0.7 m with variance (0.2 m) 2 . The traditional method is tested also with footstep lengths estimated from the sensor data with σ s = 0.1 m. For the linear model, the initial footstep vector is 0 with variance 1 m 2 I. The initial position variance is large, i.e., the prior is almost uninformative, which causes all models to have the first estimate at the location of the first WLAN estimate, which happens to be inside a wrong corridor.
In filtered routes, the traditional models start almost equally, and in the end, the traditional model that estimates the footstep length is close to the proposed linear method. The smoothed routes have some difference near the beginning of the route. The traditional models are less accurate than the linear model. This is caused by the wrong linearization of the F matrix in the beginning of the track. The smoothing was done using these F matrices in Equation 23; better results might be obtained with  more complex methods such as the unscented Rauch-Tung-Striebel smoother [27]. Some error statistics are given in Table 1. In the table, θ 0 is the initial error on the heading. 'Static' denotes the WLAN-only position estimates without filtering. When there are only footstep measurements without WLAN measurements, the 'Static' uses the last position computed from WLAN measurements as the position estimate. In the end of the table, there are results for linear model and static estimates when the WLAN measurement interval is increased so that there is always at least 15 s between consecutive measurements. A long WLAN scan interval may be used to conserve energy.
All the filtering methods result in similar accuracy with 0°and 45°initial heading errors. When the initial heading is 90 or more degrees off, neither of the traditional methods improves "Static" when the route is filtered. The UKF is more accurate than the EKF in filtering in almost every test situation, but in smoothing, the results are opposite. This indicates that the UKF linearizations of this model with uncertain prior do not work well with the Rauch-Tung-Striebel smoother. The longer WLAN scan intervals make the static estimate 3.0 m worse, filtered estimate 2.0 m worse, and the smoothed estimate 1.1 m worse.
Footstep length approximations given by different models are shown in Figure 8. The footstep lengths are com-puted in the case shown in Figure 7 (θ 0 = 90°). Here, we see that the footstep length estimated with the traditional model (s) and the proposed linear model (||v||) are rather different in the beginning of the track, but become similar towards the end of the track. The estimated footstep lengths are in general shorter than lengths computed from the sensor data.
In Table 2, the effect of CA positioning improvements presented in Subsection 3.2 is shown. The 50% column shows the 50% consistency value of the estimate, i.e., how often the position estimate mean and true position are closer to each other than 50% of the probability mass. The 95% column is defined analogously. Consistencies are computed taking only position variables into account. The closer the values are to 50% and 95%, the better. Good consistency is important for the integration with the PF because the initialization and integrity checks use the uncertainty of the KF estimate.
The position estimate with and without MIMO compensation are almost the same, but the 50% and 95% columns are much worse without the MIMO compensation. The outlier detection improves all results slightly, but not much. This is probably because the positioning data had only a few CAs that were considered outliers. The use of minimum sizes for coverage areas improved both the positioning accuracy and consistency. Even though the estimate with all WLAN improvements has the best consistency values, the error estimates are still too small. This is probably because the radio map is collected inside the building and the CA positioning can be interpreted as computation of a weighted mean. Because there are no APs that have CA center outside the building, the estimates near the edges of the building get biased towards the inside of the building. Figure 9 contains a smoothed track showing this effect. Position estimates that have the position error within 95% error ellipse are shown with a narrow blue line, and the estimates where error is outside 95% error ellipse are shown with a thicker red line. From the figure, it is evident that most of the red lines are close to edges of the building.
In Figure 10, 50% ellipses are plotted on different parts of a track. The ellipse labels show the relative number of particles required for PF initialization (Equation 63). In the figure, most residuals (red lines) are slightly out from the 50% ellipses as expected from the consistency data in Table 2, but none of them is completely wrong. In the beginning of the track, the required number of particles is 3· 10 6 even though the 50% ellipse covers less than half of the building. The next plotted ellipse on the track has 350 particles and has similar size to the remaining ellipses, which all have fewer than 200 particles. This is because the footstep vector part of state has still more uncertainty than in the rest of the track.

Conclusions
We proposed in this paper a novel linear model for PDR in indoor personal positioning and compared it to models that are common in the literature. The evaluation shows that although the model is simpler than the traditional methods, it performs well and is especially suited for situations where the initial heading and position are not known. The model assumes that uncertainty of state transition noise is equal in every direction, but our simulations show that the model is not very sensitive to different footstep length and heading variances. Because the proposed model is linear, it can be smoothed optimally with the Rauch-Tung-Striebel smoother.
In this paper, we used linear Gaussian coverage area models to do positioning with WLAN and proposed improvements for the positioning method. The results show that the proposed improvements reduced positioning errors from 7.9 to 7.0 m and significantly improved the estimate's consistency.
In addition to stand-alone position estimation, the proposed model with a KF can be used to initialize and monitor the integrity of PF. Because the linear KF is a computationally light algorithm compared to a PF, running it together with a PF does not significantly increase computational burden.

Proofs concerning MIMO compensation Proof that two identical CAs with 1 · σ distance are considered independent
Let L be a matrix square root of covariance matrix, i.e., LL T = P. Let means μ i and μ j be 1·σ from each other, i.e., Let unitary matrix U have L −1 (μ i − μ j ) as its first column. Now, Ue 1 = L −1 (μ i − μ j ), where e 1 is the first column of an identity matrix.
Using these the nominator of Equation 39 becomes: For such cases, W i,j = 2 and the MIMO correction is zero.

Proof that a set of identical CAs act as one
For identical CAs, the MIMO correction is 1. If there are n identical CAs, each covariance is multiplied by n. The resulting compensated covariance matrix of n measurements is as follows: which is the same as the covariance of one CA.

Competing interests
The authors declare that they have no competing interests.