WLAN fingerprinting based indoor positioning in the presence of censored and dropped data / Manh Kha Hoang ; Erster Gutachter: Prof. Dr.-Ing. Reinhold Häb-Umbach, zweiter Gutachter: Prof. Dr.-Ing. Peter A. Höher

WLAN Fingerprinting based Indoor Positioning in the Presence of Censored and Dropped Data Von der Fakultät für Elektrotechnik, Informatik und Mathematik der Universität Paderborn zur Erlangung des akademischen Grades Doktor der Ingenieurwissenschaften (Dr.-Ing.) genehmigte Dissertation von M.Sc. Manh Kha Hoang Erster Gutachter: Prof. Dr.-Ing. Reinhold Häb-Umbach Zweiter Gutachter: Prof. Dr.-Ing. Peter A. Höher Tag der mündlichen Prüfung: 04. März 2016 Paderborn 2016 Diss. EIM-E/321 Acknowledgments First of all, I would like to give a great thanks to Prof. Dr.-Ing. Reinhold Häb-Umbach for being my research supervisor. I would have never been able to complete my research work and my dissertation without his continuous support, understanding and patience. He has kindly motivated me for the new challenges and patiently guided me to overcome the difficulties. I feel very lucky to have met Prof. Häb-Umbach and worked with him. I would also like to thank Prof. Dr.-Ing. Peter A. Höher who have evaluated my dissertation and given me the valuable responses to improve its quality. Moreover, I would like to thank Vietnamese government for granting the scholarship for my first three years of studying in Germany. Without it, I would have never had such a good chance to commence this PhD degree in Germany. I also owe thanks to the Department of Communications Engineering, University of Paderborn, for providing me the living expense for my last two years of my research, in particular, once again thanks to Prof. Häb-Umbach. I really appreciate all the support from all staff members of the World University Service who have introduced me to Prof. Häb-Umbach, helped me to arrange the accomodation when I first came to Germany and continuously encouraged me during my research. Special thanks to my colleagues and students at the Department of Communications Engineering who have always supported me to proceed my research and improve my dissertation. In particular, I wish to thank Dr.-Ing Jörg Schmalenströer for all of his support, valuable technical ideas and comments, without that it would be very difficult for me to complete my work. Last but not least, I wish to thank my wife, Thi Hien Trang Vu, and my little daughter, Tra My Hoang, for without their encouragement, support and patience, this work would have not been completed. In addition, I would like to thank my parents who always beside, continuously support and encourage me for not only these years but also whole my life. 3 Contents 1. Introduction 1 2. Fundamentals of Indoor Positioning and State of Research 5 2.1. Proximity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2. Lateration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3. Angulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4. Fingerprinting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.5. Dead Reckoning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.6. Comparison of Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3. Scientific Objectives 18 4. Parameter Estimation of Censored and Dropped Gaussian Data 20 4.1. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.2. Parameter Estimation using EM Algorithm . . . . . . . . . . . . . . . . . 22 4.2.1. Introduction to the EM Algorithm . . . . . . . . . . . . . . . . . . 22 4.2.2. EM algorithm for Censored Gaussian Data . . . . . . . . . . . . . 23 4.2.3. EM algorithm for Censored and Dropped Gaussian Data . . . . . . 30 4.3. Optimal Classification Rule for Censored and Dropped Gaussian Data . . . 33 5. Smartphone Adaptation within the MLLR Framework 35 5.1. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 5.2. Model Adaptation in the Presence of Censored and Dropped Data . . . . . 36 5.2.1. Mean Adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 5.2.2. Variance Adaptation . . . . . . . . . . . . . . . . . . . . . . . . . 40 5.3. Regression Classes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 6. Hidden Markov Model for Indoor User Tracking 43 6.1. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 6.2. Hidden Markov Model for Indoor User Tracking . . . . . . . . . . . . . . 43 6.3. Forward Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 6.4. Movement Vector Estimation . . . . . . . . . . . . . . . . . . . . . . . . . 46 6.4.1. Step Detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 6.4.2. Movement Heading Estimation . . . . . . . . . . . . . . . . . . . 50 6.4.3. Movement Vector Calculation . . . . . . . . . . . . . . . . . . . . 51 6.5. Introduction of Pseudo States . . . . . . . . . . . . . . . . . . . . . . . . . 52 i Contents ii 7. Experimental Results on Indoor Positioning 55 7.1. Parameter Estimation and Classification . . . . . . . . . . . . . . . . . . . 55 7.1.1. EM algorithm for censored data . . . . . . . . . . . . . . . . . . . 55 7.1.2. EM algorithm for censored and dropped data . . . . . . . . . . . . 58 7.2. Smartphone Adaptation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 7.2.1. Classification on Artificial Data . . . . . . . . . . . . . . . . . . . 59 7.2.2. Classification on Field Data . . . . . . . . . . . . . . . . . . . . . 60 7.3. HMM for Indoor User Tracking . . . . . . . . . . . . . . . . . . . . . . . 62 7.3.1. Artificial Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 7.3.2. Field Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 8. Server Based Indoor Navigation System 65 8.1. Overview of Server Architecture . . . . . . . . . . . . . . . . . . . . . . . 68 8.2. Database . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 8.2.1. Map Tile Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 8.2.2. RSSI data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 8.3. Shared Memory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 8.4. Overview of Smartphone Application . . . . . . . . . . . . . . . . . . . . 74 8.5. System Operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 8.5.1. Data Gathering . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 8.5.2. Localization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 8.5.3. Navigation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 8.5.4. Features Under Development . . . . . . . . . . . . . . . . . . . . . 83 8.6. Communication Sercurity . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 9. Conclusions 88 A. Appendix 91 A.1. Derivation of EM Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 91 A.1.1. Computation of I0. . . . . . . . . . . . . . . . . . . . . . . . . . 91 A.1.2. Computation of I1. . . . . . . . . . . . . . . . . . . . . . . . . . 91 A.1.3. Computation of I2. . . . . . . . . . . . . . . . . . . . . . . . . . 92 A.1.4. Computation of W Entries . . . . . . . . . . . . . . . . . . . . . . 94 A.1.5. Computation of Information Matrix I. . . . . . . . . . . . . . . . 95 A.2. Parameters of Kalman Filter . . . . . . . . . . . . . . . . . . . . . . . . . 99 A.3. Android Smartphone Application . . . . . . . . . . . . . . . . . . . . . . . 99 A.3.1. Manifest File . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 A.4. Requests and Responses of the Communication in Indoor Navigation System 100 A.4.1. Set FringerPrint Request and Response . . . . . . . . . . . . . . . 100 List of abbreviations 102 Notations and Symbols 104 List of figures 108 List of tables 111 Contents iii References 112 List of own publications 119 “In the age of automation the ability to navigate persons and devices in indoor environments has become increasingly important for a rising number of applications.” Rainer Mautz [1] “Despite its current limitations, indoor navigation’s huge potential economic and sociological capabilities are pushing forward the research, development, implementation, and sale of low-cost systems. In the near future, this will change the way we interact with our surroundings, with many advantages for the various stakeholders involved in any indoor business.” Andrea Bottino, Giovanni Malnati and Paolo Montuschi [2] 1. Introduction Positioning and navigation have played important roles in many aspects of human civilization for thousands of years. Demand of this service is increasing steadily in many aspects called location-based services (LBS) or location-aware systems [3, 4] such as transportation, security, social networking, marketing, and so on. While positioning or localization is the process to determine the coordinates of the target objects, navigation is the process to estimate the route from one location to another. At the beginning, positioning and navigation techniques were investigated for outdoor environment to support the transportation. Outdoor positioning and navigation were done by a combination of celestial measurements and land marks, which helped people to control the vessel over a very large distance. With the development of technology, satellite was invented and subsequently, satellite based navigation systems such as Global Positioning System (GPS) by United States and Global Navigation Satellite System (GLONASS) by Russia were developed. At the moment, there are two other systems being developed namely Galileo by European Union and the BeiDou Navigation Satellite System by China. Among them, GPS [5, 6] is the most popular and successful navigation system. GPS was first developed for military purpose by the U.S. Army and it was made available for civilian use in the 1990s. With GPS, position estimate in three dimentions is obtained by applying a circular lateration technique, which relies on range measurements. Nowadays, GPS is used to support the transportation of most of the vehicles such as airplanes, cars and ships. Moreover, GPS can now be used on most smart devices (smartphones and tablet computers) to support the daily activities with lower positioning accuracy than the traditional GPS devices since lower-quality GPS chipsets are used on portable devices due to price limitation. In the ideal case of line-of- sight condition, i.e., no or almost no obstacles between the device and the satellites, GPS can locate a receiver with an accuracy of a few meters. However, in an urban area, the accuracy degrades dramatically due to non-line-of-sight problems. Beside the successful GPS, information from some other sources, e.g., Global System for Mobile Communications (GSM), are also used for the outdoor positioning and navigation purposes. GSM based positioning was developed to satisfy the Enhanced 911 (E-911) mandate from the US Federal Communications Commission which requires cellular providers to track the location of their subscribers to within 50 m for over 67% of the time. The GSM positioning system may either solely employ the GSM information [7, 8, 9] or combine it with other sources of information, e.g., GPS, inertial sensors and WiFi, resulting in hybrid systems [10, 11] in order to produce a more accurate position estimation. In an indoor environment, GPS signals are normally blocked or unreliable due to the attenuation of signal through roofs or walls. As a result, the positioning accuracy is poor. Therefore, developing a reliable indoor positioning and navigation system is of a particular interest. This topic has been attracted the consideration of many researchers over the last de- 1 Fundamentals of Indoor Positioning and State of Research 8 where x=x y(2.4) A=   x2−x1y2−y1 . . .. . . xN−x1yN−y1   (2.5) b=1 2   (x2 2+y2 2)−(x2 1+y2 1)−(r2 2−r2 1) . . . (x2 N+y2 N)−(x2 1+y2 1)−(r2 N−r2 1)   (2.6) The LS solution for the above system of equations can be obtained as follows: x= (ATA)−1ATb.(2.7) BS1 BS2 BS3 M r1 r2 r3 Figure 2.2.: Cicular lateration based positioning: the estimated position of the mobile user Mis determined based on the estimated distances between the mobile user and the reference points. At least three reference points are needed to calculate the user position. •Hyperbolic lateration: These methods use range differences between the mobile object and any pair of reference points to form a system of equations of hyperbolas (Eq. (2.8)) which are then used to compute the position of the receiver: dij =ri−rj =p(x−xi)2+ (y−yi)2−q(x−xj)2+ (y−yj)2,(2.8) here, dij denotes the range difference between i-th and j-th reference points, ∀i, j;i6= j. The solution can be easily obtained by applying the LS method in a similar fashion as discussed in the circular lateration based method. The above system of equations can Fundamentals of Indoor Positioning and State of Research 9 be solved as follows [3]: Since di1=ri−r1,(2.9) we have (r1+di1)2=r2 i ⇔r2 i−r2 1=d2 i1+ 2r1di1(2.10) then (x−xi)2+ (y−yi)2−(x−x1)2−(y−y1)2=r2 i−r2 1 =d2 i1+ 2r1di1 ⇔x2 i+y2 i−(x2 1+y2 1)−2x(xi−x1)−2y(yi−y1) = d2 i1+ 2r1di1 x(xi−x1) + y(yi−y1) + di1r1=1 2(x2 i+y2 i)−(x2 1+y2 1)−d2 i1; (2.11) where i= 1,···, N. The above system of equations in Eq. 2.11 can be written in matrix form: Ax =b,(2.12) where x=  x y r1 (2.13) A=     x2−x1y2−y1d21 . . .. . . . . . xN−x1yN−y1dN1     (2.14) b=1 2   (x2 2+y2 2)−(x2 1+y2 1)−d2 21 . . . (x2 N+y2 N)−(x2 1+y2 1)−d2 N1   (2.15) The LS solution for the above system of equations is then obtained as: x= (ATA)−1ATb.(2.16) Many indoor localization methods employing lateration techniques have been proposed, see the descriptions in [46, 47, 48] and the references therein. The positioning accuracy depends on how accurate the estimate of the distance between transmitter and receiver is. Lateration based techniques can be divided into 2categories: time based lateration and signal strength based lateration. Fundamentals of Indoor Positioning and State of Research 10 Time based lateration techniques employ the TOA or TDOA information to compute the distances. In [49], a time based lateration positioning system has been developed. The object to be located is the mobile tag transmitting ultrasonic pulses in a roughly hemispherical pattern around the top of it, once requested from a master station via a radio signal. A central positioning unit polls the network of the ultrasonic receivers mounted on the ceiling of the room to retrieve the time of flight (TOF) of the ultrasound emitted by the mobile tag (if detected at the receiver). The TOF information is measured as the time difference between the moment the master sends out the request to mobile tags for emitting ultrasonic pulse and the first signal peak detected on the receiver. For time synchronization, the master station sends the reset signal to all receivers simultaneously with the ultrasound emitting request. Distances between mobile tag and receivers are then computed and the position of the mobile tag is determined by applying the LS approach. The reported experimental results are very impressive with a positioning error in the order of centimeters. In addition to time based lateration, signal strength based lateration techniques are also employed in indoor positioning. These techniques employ the dependence of signal strength on propagation distance. Therefore, the accuracy of these techniques mostly depends on how accurate the estimated path loss model is. In [18], a modified path loss model for indoor environment was developed, where the attenuation of signal penetrating through walls is considered. The parameters of the path loss model are trained empirically. Once the distances between transmitters and receiver are obtained, the final estimate of receiver location is obtained by various methods, e.g., by solving the system of circular or hypecbolic equations using the LS approach as summarized in [46]. In [47], methods for improving positioning results using RSSI based lateration techniques were explored. Instead of using the theoretical path loss model, the authors proposed two approaches, regression based and correlation based. In the regression based method, polynomial regression was used to model the relationship between the RSSI and distance, and the coefficients of the polynomial are estimated by employing LS approximation. The correlation based method utilizes the fact that the signal propagation from close-by locations to an AP is highly correlated as they face the similar propagation environment. This method recursively performs localization in gradually reduced local areas to approach the true location. In each iteration, a certain number of RSSI readings which are the closest training points to the previous estimated position are kept and used to fit the propagation model. Experimental results on both artificial data and real data showed the considerable improvements compared to the original RSSI based lateration methods. An overview of lateration based techniques for Ultra-Wideband signals is given in [48] where positioning schemes for time based lateration and RSSI based lateration are discussed. In general, time based lateration positioning techniques are more precise than RSSI based lateration techniques, however, they require additional hardware. 2.3. Angulation Angulation based positioning techniques employ the angle of arrival of a wireless signal to determine the position of a mobile object. Theoretically, these techniques are able to compute the position of the receiver when there are at least two reference points, and if the position Fundamentals of Indoor Positioning and State of Research 11 of the mobile object does not lie on the line connecting the reference points. BS1BS2 M α1α2 Figure 2.3.: Angulation based positioning: the estimated position of the mobile user Mis determined based on the estimated angle αibetween the mobile user and the reference points. At least two reference points are needed to calculate the user position with the condition that the user position does not lie on the line connecting the reference points. Fig. 2.3 illustrates the angulation based positioning techniques. In this system, the connected line between two references can be considered as an internal reference. The angle between the transmitter and the receiver αican be determined by: tan αi=y−yi x−xi ⇔xsin αi−ycos αi=xisin αi−yicos αi;i= 1,···, N. (2.17) However, these techniques also suffer from the errors caused by the NLOS problem. If there are more than two reference points, the system of equations might not produce a unique solution. To solve the system of equations, an approximate solution is needed. Again, the LS method is the most common solution for solving the system of equations which is defined by Eq. (2.17). Eq. (2.17) can be written in matrix form as follows Ax =b,(2.18) where x=x y(2.19) A=  −sin α1cos α1 . . .. . . −sin αNcos αN   (2.20) b=   y1cos α1−x1sin α1 . . . yNcos αN−xNsin αN   (2.21) Fundamentals of Indoor Positioning and State of Research 12 The LS solution for the above system of equations can be obtained as follows: x= (ATA)−1ATb.(2.22) An indoor satellite positioning system has been proposed in [38], however, different from the satellite based GPS. This system employs the infrared signals instead of eletromagnetic wave and angulation technique instead of lateration technique. The system consists of 3 infrared light sources, called indoor satellites, mounted at fixed positions as emitters, and the infrared incident angle sensors at the mobile objects measuring the incident angle from each emitter. The position of the mobile objects are then obtained by applying the LS approach. In [50], angulation based techniques have been employed to estimate the current position of an object. ML, LS, Total Least Squares (TLS) and Weighted LS (WLS) algorithms were applied to solve the system of equations. To improve the performance of the positioning system, a method based on Weighted TLS (WTLS) was derived. The optimistic results with WTLS based approach were presented with simulation data which approach the Cramer-Rao Lower Bound (CRLB). Another example of an AOA based localization system can be found in [51]. The system consists of a set of passive thermal infrared sensors installed in the room edges to detect the thermal radiation of the human skin. The thermal sensors are thermopile-arrays, where each contains a number of pixels and each pixel has a field of view. The heat source position is determined via the principle of AOA by computing the intersection point created by the directions of the pixels with the highest outputs. 2.4. Fingerprinting Fingerprinting based positioning techniques are the methods to estimate the position of an object which rely on training data from a set of reference points (anchor points) with known locations. Fig. 2.4 illustrates such a fingerprinting based system. BS1 BS2 BS3BS4 M Reference Point Figure 2.4.: Fingerprinting based positioning: the filled circles indicates the reference (anchor) points, while the triangles show base station locations. The user position is the position of the anchor point whose training data best match the online measurement Fingerprinting based methods can be well formulated as a machine learning and pattern recognition approach. These techniques generally consist of two phases: training phase and classification phase (also refered as offline phase and online phase). In the training phase, Fundamentals of Indoor Positioning and State of Research 13 the training data, i.e., RSSI, are collected at the anchor points and used to build the database which is often called radio map which describes the RSSI-position relationship. The radio map can be defined as follows R={(ℓk,Fk)|k= 1,···, K}(2.23) where ℓkis the k-th reference position, Kis the number of reference positions in the deployment area, Fkcan be either the set of raw measurements, i.e., Fk=Xk=xk,1,···,xk,N , where xk,n = [xk,n,1,···xk,n,NAP ]Tis the n-th measurement vector which contains the RS- SIs from NAP base stations, or it contains the class conditional probability density functions (PDFs), i.e., Fk=p(x|ℓk), estimated from the measurements at the k-th position. During the classification phase, the online measurements are compared against the training data at every anchor point. The position of the anchor point whose training data best match the online data can be considered as the estimated position of the receiver. As discussed in chapter 1, there are two most common approaches for calculating the similarity between training data and online data: deterministic approaches, e.g., the k-nearest neighbor rule, and probabilitic approaches, e.g., the Bayesian classification rule. Though radio fingerprinting based positioning techniques can be applied to both indoor and outdoor environments, it seems that these techniques are more suitable for indoor since the distinction of radio signal strengths observed at different positions indoor is much higher than outdoor due to the density of obstructions in indoor environments. In the very well-known fingerprinting based positioning system Radar [18, 21], the knearest neighbor method is employed. The position of the receiver is computed by averaging the coordinates of the kanchor points which have highest similarity between the training data and the online observation. The similarity is measured by computing the Euclidian distance between the training data and online observation in signal strength space as follows D(o,xk,n) = v u u t NAP X i=1 (oi−xk,n,i)2,∀k, ∀n(2.24) where o= [o1,···, oNAP ]Tand Ddenote the online observation and the distance, respectively. The j-th position is considered to be the nearest neighbor, i.e., the training data that best match the online observation, if there is a training measurement xj,m which satisfies D(o,xj,m)≤ D(o,xk,n),∀k6=j, ∀n(2.25) In [22], an extensive analysis of the fingerprinting based positioning system that employs the Euclidian distance is presented. The effect of the number of access points, the number of training samples per position, the density of the anchor points, etc. on the performance of the positioning system are analyzed. The analysis results provide a guideline on choosing parameters to design and deploy an indoor positioning system. In probabilitic approaches, in the training phase the statistical parameters of the PDF of the signal strength are trained. During the classification phase, the similarity between training data and online data is computed by calculating the likelihood of observing the online data given the PDF of signal strength at anchor points [24, 25, 8], or some other criteria to measure the similarity such as the Bhattacharyya distance [52], then again a k-nearest neighbor algorithm can be applied to compute the final estimated position of the receiver. Fundamentals of Indoor Positioning and State of Research 14 There are two methods to estimate the PDF of the training data, parametric and nonparametric density estimation techniques [3]. Parametric estimation methods assume a model, e.g., a Gaussian, for the density function and aim to estimate the parameters of the model [25, 8], i.e., Fk=θk= (µk,Σk). On the other hand, non-parametric methods do not assume any prior model but estimate the class conditional PDF by using the histogram method [24, 52] or Kernel density estimators [23, 53]. After constructing the class conditional PDF, the final position estimate can be obtained by computing the likelihood or the posterior using the Bayes’ rule. The position which achieves the highest likelihood or posterior is then considered as the position of the mobile receiver. Maximum likelihood estimator aims to find the position that maximizes the probability of observing the online measurement o, given the class conditional PDFs at reference points, as follows ˆ ℓ= argmax ℓk p(o|ℓk).(2.26) A maximum a posteriori estimator finds the position which has the highest posterior probability as ˆ ℓ= argmax ℓk P(ℓk|o).(2.27) Applying Bayes’ rule, P(ℓk|o)can be obtained as follows P(ℓk|o) = p(o|ℓk)P(ℓk) p(o) =p(o|ℓk)P(ℓk) PK k′=1 p(o|ℓk′)P(ℓk′),(2.28) assuming that the prior P(ℓk)is given. In [25], two different types of class conditional PDF are considered, one is the parametric model (Gaussian) and the other is the non-parametric model. In the reported results, the method employing the parametric model outperforms the other. This is because, as stated, the parametric technique smooths the distribution shape to account for missing signal strength values in the training phase, (due to the finite number of training samples) which avoids obtaining a zero probability for any signal strength value that was not observed in the training phase. 2.5. Dead Reckoning Dead reckoning (DR) techniques employ the information from a system of inertial sensors, i.e., acceleration and gyroscope sensors, in order to estimate the position of an object based on the previous estimated position. DR techniques had been first developed for the aviation and marine industry, nowadays they are also used in robotics [54, 55, 56], indoor positioning [57, 58, 59, 32], and many more systems. Fig. 2.5 illustrates the positioning procedure of DR techniques where the current position is estimated based on the previous estimate. vtis the movement vector computed from displacement information and movement heading estimation. Fundamentals of Indoor Positioning and State of Research 15 With the development of Microelectromechanical systems (MEMS) technology, ordinary smartphones now often have a built-in inertial measurement unit (IMU) which allows to employ the dead reckoning techniques for indoor pedestrian tracking without additional hardware. However, they have lower accuracy in comparison with the IMU system used on airplanes or in space science. As discussed in chapter 1, these techniques are often developed in a combination with other methods using other sources of information, such as WiFi signal or GPS, which are used to regularly reset the error of DR techniques which accumulates over time and distance. Displacement estimation of the indoor user in DR based indoor positioning is done via step detection and step length estimation where the observed acceleration data are used. Various methods have been proposed for step detection procedure such as peak detection [57], zero-crossing detection as mentioned in [59] or autocorrelation approach [32]. Step length estimation based on acceleration data is a tricky task since the accuracy of built-in sensors on smartphones are often very low. As a result, most of the step length estimation methods are based on a calibration phase to estimate some heuristic factors, as summarized in [60]. A Kalman filter is often used for movement heading estimation which utilizes the information from gyroscope and magnetic sensors [59]. v1 v2v3 Figure 2.5.: Dead reckoning: the user position is estimated based on the estimated movement vectors assuming that the starting point is given. 2.6. Comparison of Techniques All the above mentioned radio signal based techniques for positioning have their own advantages and disadvantages and their suitability depends on the deployment environment. •Proximity based techniques can be applied for any low cost positioning system for either outdoor or indoor (no requirement of additional hardware) which does not require a high positioning accuracy. The positioning error depends on the coverage range of the base stations, e.g., the proximity positioning systems using GSM base stations have error in the order of hundreds of meters. An indoor positioning system may use these techniques employing the wireless local area network (WLAN) system. However the positioning error is about the coverage range of an AP, i.e., roughly 100 m which is unacceptable in the indoor environment. Moreover, in indoor environments, the contour of signal strength radiated from an AP is obviously not symmetric due to the presence of obstructions. This leads to misclassification results, because the receiver may receive a higher RSSI from a farther AP because of line of sight to this AP while it observes a lower RSSI from a closer AP due to absence of line of sight. Positioning using iBeacon advertising packets is another approach which aims to produce a room precise localization system. However, the higher the accuracy requirements, the more hardware (iBeacon transmitters) needs to be installed. Fundamentals of Indoor Positioning and State of Research 16 •Lateration and angulation based techniques are able produce highly precise positioning results in free space environment where signals can propagate by the direct path from the transmitter to the receiver and almost no multipath problem occurs. Unfortunately, the indoor environment is definitely not free space. It has a lot of obstructions such as walls, furniture or movement of people that makes the accuracy of the positioning using these techniques very limited. Because of the special requirements of the environment, these techniques may only be suitable for positioning in a single room. Another limitation of these techniques is the requirement of additional specialized hardware, i.e., antenna, microphone arrays, etc., and deep hardware access which seems to be not possible for a positioning system employing the existing WLAN infrastructure and ordinary mobile devices, i.e., smartphones. Clock synchronization is another challenge for lateration based techniques which determines the accuracy of positioning results. •Fingerprinting based methods as discussed above consist of two phases, training phase and classification phase. The first and the foremost disadvantage of these techniques is the necessity of a training phase which is very time consuming and need to be done very carefully because the classification results are strongly influenced by the quality of the training data. Another disadvantage of these techniques is that the training data need to be updated regularly in order to adapt to the changes of the infrastructure and environment in order to maintain the accuracy of the positioning results. Despite these disadvantages, fingerprinting based techniques are still the most promising methods for indoor environments for three main reasons: First, they can produce a reasonable positioning result without any additional special hardware. Second, these techniques are suitable to be employed in an already installed WLAN system, and WLAN is available in many places. Third, systems using these techniques can be applied in a large area and have a good scalability. •Dead reckoning techniques are able to produce precise positioning results given the initial position without any knowledge of the infrastructure in the coverage area, however, only over a short time and movement distance. The high accuracy is not kept long because of the accumulation of error over time and distance of movement. The errors are caused, for example, by the drift or bias of sensor systems, and also the noise in the environments, for example the measurement from a magnetic sensor is strongly affected by metallic materials arround the sensor and not only the magnetic field of the earth. From the analysis of the the advantages and disadvantages of the mentioned radio signal based techniques, for developing a real indoor positioning system based on the regular smartphones and the already installed WLAN system, the combination of fingerprinting based and dead reckoning techniques appears to be the best solution because of the following reasons: •Sufficient accuracy: The positioning accuracy meets the requirement of indoor positioning for the common purposes such as user tracking and security, in large indoor areas, i.e., universities, museums, hospitals or stores. •Cost efficiency: no additional hardware is needed to build up a positioning system. Fundamentals of Indoor Positioning and State of Research 17 •Deployment simplicity: WLAN is today available in almost every indoor environment. Consequently, the positioning system can be deployed in many indoor areas easily. •Robustness: The combination of WiFi fingerprinting based and DR techniques makes the system robust, i.e., the mobile devices can do self localization based on DR techniques when no WiFi signal is available. Parameter Estimation of Censored and Dropped Gaussian Data 24 used with the same meaning. The goal is to estimate the parameters θ={µ, σ2}of the underlying Gaussian. E-Step Employing the EM algorithm we identify yand xto be the complete and the observed data, respectively. Thus the expected log-likelihood of the complete data is given by Q(θ;θ(κ)) = Eln (pY(y;θ)) |x;θ(κ)(4.4) = N X n=1 Z∞ −∞ ln (pY(yn;θ)) pyn|xn;θ(κ)dyn(4.5) =: N X n=1 fnθ;θ(κ),(4.6) where fnθ;θ(κ)=R∞ −∞ ln (pY(yn;θ)) pyn|xn;θ(κ)dynand where κis the iteration index. The term p(yn|xn;θ(κ))can be determined as follows: p(yn|xn;θ(κ)) = (δ(yn−xn), if xn> c p(xn|yn;θ(κ))p(yn;θ(κ)) p(xn;θ(κ)), if xn=c(4.7) For the case xn=cwe know that yn≤c.p(yn|xn;θ(κ))can thus be calculated as follows p(yn|xn;θ(κ)) = p(xn|yn;θ(κ))p(yn;θ(κ)) p(xn;θ(κ)) =p(xn|yn;θ(κ))p(yn;θ(κ)) Rc −∞ p(xn|yn;θ(κ))p(yn;θ(κ))dyn =δ(xn−c)p(yn;θ(κ)) Rc −∞ δ(xn−c)p(yn;θ(κ))dyn =Nyn;θ(κ) I0(θ(κ)) Here we have used the notation Ij(θ(κ)) = Zc −∞ yjNy;θ(κ)dy(4.8) with j= 0.I0(θ(κ))can be calulated as follows (see Appendix A.1.1 for the detailed computation): I0(θ(κ)) = 1 2erfc −c−µ(κ) √2σ(κ).(4.9) Then we obtain p(yn|xn;θ(κ)) = (δ(yn−xn), if xn> c N(yn;θ(κ)) I0(θ(κ)), if xn=c.(4.10) Parameter Estimation of Censored and Dropped Gaussian Data 25 Introducing the binary random variable Znwith realization zn, where zn= 0 and zn= 1 indicate that the n-th measurement is not censored or censored, respectively, the summand in Eq. (4.6) can be written as fnθ;θ(κ)=zn I0(θ(κ))Zc −∞ ln (N(yn;θ)) Nyn;θ(κ)dyn + (1 −zn) ln (N(xn;θ)) .(4.11) here, the fact that yn=xnwhen zn= 0 is exploited. M-Step The parameter re-estimation formulas are obtained by computing the derivatives of Eq. (4.11) w.r.t. the elements of θand set them to zero: ∂ ∂µfn(θ;θκ) = zn I0(θ(κ))Zc −∞ yn−µ σ2Nyn;θ(κ)dyn + (1 −zn)xn−µ σ2 =zn σ2 I1θ(κ)−µI0θ(κ) I0(θ(κ))+ (xn−µ)1−zn σ2(4.12) ∂ ∂σfn(θ;θκ) = zn I0(θ(κ))Zc −∞ −1 σ+(yn−µ)2 σ3Nyn;θ(κ)dyn + (1 −zn)−1 σ+(xn−µ)2 σ3 =zn σ1 σ2I2(θ(κ)) I0(θ(κ))−2µI1(θ(κ)) I0(θ(κ))+µ2−1 +1−zn σ(xn−µ)2 σ2−1,(4.13) Here, we have used the notation as defined in Eq. (4.8). Setting the sumations of the above derivatives to zeros, re-estimation formulas are readily obtained: N X n=1 ∂ ∂µfn(θ;θκ)µ=µκ+1 = 0 : µ(κ+1) =1 N I1(θ(κ)) I0(θ(κ)) N X i=1 zn+1 N N X i=1 (1 −zn)xn(4.14) N X n=1 ∂ ∂σfn(θ;θκ)σ2=(σ2)κ+1 = 0 : σ2(κ+1) =I2(θ(κ)) I0(θ(κ))−2µ(κ)I1(θ(κ)) I0(θ(κ))+µ2(κ)1 N N X i=1 zn +1 N N X i=1 (1 −zn)xn−µ(κ)2.(4.15) Parameter Estimation of Censored and Dropped Gaussian Data 26 The computation of I1(θ(κ))and I2(θ(κ))are given in Appendix A.1.2 and Appendix A.1.3, respectively. As can be seen in Eq. (4.14) and Eq. (4.15), the unobservable data, i.e., censored data, also contribute to the estimates beside the observable ones. Intuitively, in case of absence of censored data, i.e., zn= 0 for all measurements, the re-estimation formulas reduce to the typical ML estimation formulas for mean µand variance σ2. To have a better intuition of the re-estimation formulas, the situation after convergence of the estimates is considered. Let µ(κ+1) ≈µ(κ)=: ˆµand (σ2)(κ+1) ≈(σ2)(κ)=: ˆσ2, using this in Eq. (4.14) and Eq. (4.15) and solving for the estimates, we arrive at ˆµ=M N 1 M M X i=1 xn+1−M NRc −∞ ypY(y;ˆ θ)dy Rc −∞ pY(y;ˆ θ)dy,(4.16) ˆσ2=M N 1 M M X i=1 (xn−ˆµ)2+1−M NRc −∞(y−ˆµ)2pY(y;ˆ θ)dy Rc −∞ pY(y;ˆ θ)dy,(4.17) where we assumed without loss of generality that the first M=Pi(1 −zn)observations are the uncensored ones. These expressions lend themselves to the following interpretation: Mean and variance estimates are the weighted average between their ML estimates which are computed from the observed data and the mean and variance of the assumed truncated Gaussian of the unobservable parts. The weights are the relative frequencies of the uncensored and censored measurements, respectively. Properties of the Estimates In the following, the main properties of the proposed estimator are investigated. As will be shown, the proposed EM algorithm delivers virtually biasfree and efficient estimates. Unbiasness and Convergence In order to study the convergence properties the expected values of the difference between the estimates, Eqs. (4.14) and (4.15), and the true values of the parameters are computed. First note that Znis a Bernoulli random variable with P(Zn= 1) = Zc −∞ pY(y;θ)dy=I0(θ) and P(Zn= 0) = 1 −I0(θ). Thus: E [Zn] = E[Z2 n] = I0(θ). Taking the expectation of mean estimate Eq. (4.14) we have Eµ(κ+1)=1 N N X i=1 E[(1 −Zn)Xn] + 1 N N X i=1 EI1(θ(κ)) I0(θ(κ))Zn.(4.18) It has to be noted that the estimate θ(κ)depends on the data and is thus random. Thus the expectation operator also applies to this quantity. The first expectation can be evaluated as follows: E[(1 −Zn)Xn] = 1 X zn=0 Z∞ −∞ (1 −zn)xnP(zn|xn)pY(xn)dxn =Z∞ c xnN(xn;θ)dxn=µ−I1(θ).(4.19) Parameter Estimation of Censored and Dropped Gaussian Data 27 Note that here we have used the notation pY(xn)because x=yin case of not censored. Here we have employed the fact that P(zn= 0|xn) = 1,if xn> c 0,else .(4.20) Using further E hI1(θ(κ)) I0(θ(κ))zni≈EhI1(θ(κ)) I0(θ(κ))iE[Zn]and subtracting µfrom either side of Eq. (4.18) we obtain E˜µ(κ+1)=−I1(θ) + I0(θ)EI1(θ(κ)) I0(θ(κ)).(4.21) where ˜µ(κ+1) =µ(κ+1) −µ. In order to compute the remaining expectation, a Taylor series expansion around the true parameter values θ= (µ, σ2)is applied and truncated after the linear term: EI1(θ(κ)) I0(θ(κ))≈I1(θ) I0(θ)+E˜µ(κ)∂ ∂µ I1(θ) I0(θ)+Eh˜σ2(κ)i∂ ∂σ2 I1(θ) I0(θ).(4.22) Expectation of variance estimate Eq. (4.15) can be calulated in a similar procedure as follows: Eσ2(κ+1)=E"1 N N X n=1 (1 −zn)(xn−µ)2# +E"1 N N X n=1 zn I2(θ(κ))−2I1(θ(κ))µ+I0(θ(κ))µ2 I0(θ(κ))# =1 N N X n=1 E(1 −zn)(xn−µ)2 +1 N N X n=1 Ezn I2(θ(κ))−2I1(θ(κ))µ+I0(θ(κ))µ2 I0(θ(κ))(4.23) The first expectation can be evaluated as follows: E(1 −zn)(xn−µ)2= 1 X zn=0 Z∞ −∞ (1 −zn)(xn−µ)2P(zn|xn)pY(xn)dxn =Z∞ c (xn−µ)2N(xn;θ)dxn =σ2−(I2(θ)−2µI1(θ) + µ2I0(θ)).(4.24) Using further E hzn I2(θ(κ))−2I1(θ(κ))µ+I0(θ(κ))µ2 I0(θ(κ))i≈EhI2(θ(κ))−2I1(θ(κ))µ+I0(θ(κ))µ2 I0(θ(κ))iE[Zn]and subtracting σ2from either side of Eq. (4.23) we obtain Eh˜σ2(κ+1)i=−I2(θ)−2µI1(θ) + µ2I0(θ) +I0(θ)EI2(θ(κ))−2I1(θ(κ))µ+I0(θ(κ))µ2 I0(θ(κ))(4.25) Parameter Estimation of Censored and Dropped Gaussian Data 28 where (˜σ2)(κ+1) = (σ2)(κ+1) −σ2. In order to compute the remaining expectation, a Taylor series expansion around the true parameter values θ= (µ, σ2)is again applied and truncated after the linear term: EI2(θ(κ))−2I1(θ(κ))µ+I0(θ(κ))µ2 I0(θ(κ))≈I2(θ)−2I1(θ)µ+I0(θ)µ2 I0(θ) +E˜µ(κ)∂ ∂µ I2(θ)−2I1(θ)µ+I0(θ)µ2 I0(θ) +Eh˜σ2(κ)i∂ ∂σ2 I2(θ)−2I1(θ)µ+I0(θ)µ2 I0(θ). (4.26) Writing Eq. (4.21) and Eq. (4.25) in matrix form, we arrive at E˜µ(κ+1) Eh(˜σ2)(κ+1)i!≈WµWµσ Wσµ Wσ E˜µ(κ) Eh(˜σ2)(κ)i! =WµWµσ Wσµ Wσκ+1 E˜µ(0) Eh(˜σ2)(0)i!,(4.27) where Wµ=I0(θ)∂ ∂µ I1(µ) I0(µ), Wµσ =I0(θ)∂ ∂σ2 I1(θ) I0(θ), Wσ=I0(θ)∂ ∂σ2 I2(θ)−2I1(θ)µ+I0(θ)µ2 I0(θ), Wσµ =I0(θ)∂ ∂µ I2(θ)−2I1(θ)µ+I0(θ)µ2 I0(θ). The Wentries can be readily computed by using the derivatives given in Appendix A.1.4. Investigating all the possible values of the magnitudes of the eigenvalues of the matrix with the Wentries by changing the clipping threshold cfrom −∞ to ∞, we observed that they are always less than one (approaching one for c→ ∞, i.e., no observable data). We thus conclude that the estimates are biasfree because the expectation of the error decreases exponentially over iterations. Since, however, in practice µhas to be replaced by its estimate in the estimation of the variance, we exhibit the same bias as in ordinary ML estimation of the variance. In the application considered here, the bias of the ML estimate of the variance can be neglected due to the large number of samples. Further, an estimate of the convergence speed of the EM algorithms can also be obtained from Eq. (4.27), see Fig. 4.3. It can be seen that the number of iterations quickly rises once more than 50% of the data are clipped. Precision In order to evaluate the precision of our proposed estimator, we calculated the Cramer-Rao Lower Bound (CRLB) for the estimator. CRLB of the estimation error variance of the mean and variance estimates can be obtained by inverting the information matrix Iwhich contains the expected values of the second-order partial derivatives of the log-likelihood function. From the derivation of the EM algorithm it can be seen that the measurements consist of two types of data, the number Mof noncensored observations and the observations themselves. While the former type is binomially distributed, the latter type is drawn from a truncated Parameter Estimation of Censored and Dropped Gaussian Data 29 c−µ σ No. of EM iterations −2−1.5−1−0.500.51 1.52 0 100 200 300 Figure 4.3.: Theoretical number of EM iterations κrequired to reduce the estimation error to 10−4 of its initial value as a function of (c−µ)/σ. Initial values µ(0),σ2(0) have been set to the ML estimates of µ, σ2computed from the unclipped observations only. Gaussian. It is noted that the draws from the Gaussian are independent of the binomial random variable. Following [65] the likelihood is thus given by, where I0(θ)is the probability of not observing a sample and 1 √2πσ2exp −1 2xj−µ σ2!is the probability of observing xj>−100 dBm p(x;θ) = N! M!(N−M)!IN−M 0(θ)1 √2πσ2M exp −1 2 M X j=1 xj−µ σ2!,(4.28) and the log-likelihood function is obtained as follows: ln p(x;θ) = ln N! M!(N−M)!+ (N−M) ln (I0(θ)) −M 2ln 2πσ2−1 2 M X j=1 xj−µ σ2 .(4.29) The information matrix Iis defined as: I= Eh−∂2 ∂µ2ln p(x;θ)iEh−∂2 ∂µ∂σ2ln p(x;θ)i Eh−∂2 ∂µ∂σ2ln p(x;θ)iEh−∂2 ∂σ2∂σ2ln p(x;θ)i .(4.30) The entries of matrix Ican be calculated as follows (see Appendix A.1.5 for the detailed computation): E−∂2 ∂µ2ln p(x;θ)=N σ2I2 1−I2I0 I0 + 1(4.31) E−∂2 ∂µ∂σ2ln p(x;θ)=N σ2 − ∂I1 ∂σ2I0−I1∂I0 ∂σ2 I0!(4.32) E−∂2 ∂σ2∂σ2ln p(x;θ)=N σ2 1 2σ2−1 2σ2 ∂I2 ∂σ2I0−I2∂I0 ∂σ2 I0−2µ ∂I1 ∂σ2I0−I1∂I0 ∂σ2 I0!! (4.33) where the computation of the derivative terms are given in Appendix A.1.4. Parameter Estimation of Censored and Dropped Gaussian Data 30 CRB (σ2) CRB (µ) Simulation (σ2) Simulation (µ) log10(MSE) c−µ σ −2−1.5−1−0.500.51 1.52 −2 −1 0 1 2 3 Figure 4.4.: Comparison of CRLB for mean and variance with MSE obtained from simulation for σ2= 25 and N= 1000. Computing the I−1, we obtain CRLB on estimation error variance of µand σ2: Var(ˆµ)≥I−1(1,1) (4.34) Var(ˆσ2)≥I−1(2,2) (4.35) Fig. 4.4 compares the CRLB with the mean squared error (MSE) of the proposed estimators for mean µand variance σ2obtained from a simulation. It can be seen that the estimator practically achieves the bound. The differences are so small that they are not visible in the graph. We therefore conclude that the estimator is efficient. Furthermore, for the limiting case of completely uncensored data the well-known results for ML parameter estimation from a normal population are obtained: MSE(µ) = σ2/N; MSE(σ2) = σ4(2N−1)/N2. 4.2.3. EM algorithm for Censored and Dropped Gaussian Data In this subsection, the proposed EM algorithm for parameter estimation of censored Gaussian data is extended to be able to cope not only with the censored data but also the dropped data. Notations The measurement model for censored and dropped data is shown in Fig. 4.5. Here, the hidden variables dn,n= 1,...N, indicate whether an observation is dropped (dn= 1) or not (dn= 0), where P(dn= 1) is the dropping rate, in the following, P(dn= 1) is denoted by π. The variables are gathered in the binary vector d= [d1,...,dN], where Nis the number of measurements. Furthermore, define y=y1, ..., yN, where the ynare i.i.d. with Gaussian probability density function (PDF) pY(yn) = N(yn;µ, σ2)if dn= 0, and yn=cif dn= 1, where we set cto the smallest measurable RSSI value (e.g., c=−100 dBm). It should be Parameter Estimation of Censored and Dropped Gaussian Data 31 noted that yis not the complete data anymore since the data are possibly dropped. Observable data are the censored, possibly dropped data x=x1, ..., xN, where xn= max(yn, c), with the censoring threshold cas above. xn yn max(yn, c) ∼ N(µ, σ2) c dn dn= 0 dn= 1 Dropping Censoring Figure 4.5.: Measurement model. The difference between censoring and dropping can be expressed as follows: Even if the parameters of the Gaussian are such that practically no censoring occurs (because µ≫c), dropping will still possibly occur. The goal is to estimate the parameters θ={µ, σ2, π}. This will be achieved by the Expectation Maximization (EM) algorithm, where the complete data are {y,d}and the observable are {x}. E-Step It is noted that xndoes not convey more information about θthan yn. So, the complete data are {y,d}rather than {x,y,d}. The expected log-likelihood of the complete data, given the observed one, is given by: Q(θ;θ(κ)) = Eln (p(y,d;θ)) |x;θ(κ) = N X n=1 1 X dn=0 Z∞ −∞ ln (p(yn, dn;θ)) p(yn, dn|xn;θ(κ))dyn,(4.36) where κis the iteration index. p(yn, dn|xn)can be written as p(yn, dn|xn;θ(κ)) = p(yn|dn, xn;θ(κ))P(dn|xn;θ(κ)).(4.37) For calculating p(yn|dn, xn;θ(κ)), two cases are considered: •For data that are not dropped (dn= 0), see Eq. (4.10), p(yn|dn= 0, xn;θ(κ)) = (δ(yn−xn),if xn> c N(yn;θ(κ)) I0(θ(κ)),if xn=c(4.38) •For data that are dropped (dn= 1) p(yn|dn= 1, xn;θ(κ)) = 0,if xn> c δ(yn−c),if xn=c.(4.39) Parameter Estimation of Censored and Dropped Gaussian Data 32 Furthermore, the posterior of the hidden variable dncan be computed using Bayes’ rule P(dn|xn;θ(κ)) = p(xn|dn;θ(κ))P(dn) P1 dn=0 p(xn|dn;θ(κ))P(dn).(4.40) Since there are two variables, four cases can be discerned. Using the notation βn(d, z) := P(dn|zn;θ(κ)), where zn= 1 indicates that xn=cand zn= 0 that xn> c, we obtain: βn(1,1) = P(dn= 1|xn=c;θ(κ)) =p(xn=c|dn= 1)P(dn= 1) p(xn=c|dn= 0)P(dn= 0) + p(xn=c|dn= 1)P(dn= 1) =π(κ) I0(θ(κ))(1 −π(κ)) + π(κ),(4.41) where the result from the previous section is employed: p(xn=c|dn= 0) = I0(θ(κ))is the probability that a measurement is censored, i.e., unobservable, given that it is not dropped. In addition, p(xn=c|dn= 1) is always equal to 1since if the measurement is dropped, it is assigned the value of c.βn(1,1) is the probability that the measurement is dropped given it is unobservable. The probability that the measurement is not dropped given it is unobservable, βn(0,1), is then easily obtained since βn(0,1) + βn(1,1) = 1. Furthermore, it is obvious that if a measurement is observable, it cannot be dropped, so βn(1,0) = 0, and thus βn(0,0) = 1. We further note that p(yn, dn;θ) = P(dn;θ)p(yn|dn;θ) =πδ(yn−c),if dn= 1 (1 −π)N(yn;θ),if dn= 0 (4.42) Using all the above in (4.36) we arrive at Q(θ;θ(κ)) = N X n=1 znZc −∞ ln ((1 −π)N(yn;θ)) N(yn;θ(κ)) I0(θ(κ))βn(0,1)dyn + (1 −zn)Z∞ c ln ((1 −π)N(yn;θ)) δ(yn−xn)βn(0,0)dyn +znZc −∞ ln (πδ(yn−c)) δ(yn−c)βn(1,1)dyn = N X n=1 znln(1 −π)βn(0,1) +zn I0(θ(κ))Zc −∞ ln (N(yn;θ)) N(yn;θ(κ))dynβn(0,1) + (1 −zn) ln(1 −π) + (1 −zn) ln (N(xn;θ)) +znln(π)βn(1,1).(4.43) Parameter Estimation of Censored and Dropped Gaussian Data 33 M-Step Computing the derivative of the auxiliary function Eq. (4.36) the following iterative parameter estimation formulas can be readily derived: µ(κ+1) =PN n=1(1 −zn)xn+I1(θ(κ)) I0(θ(κ))PN n=1 znβn(0,1) N−PN n=1 znβn(1,1) (4.44) σ2(κ+1) =PN n=1 znI2(θ(κ)) I0(θ(κ))−2µI1(θ(κ)) I0(θ(κ))+µ2βn(0,1) N−PN n=1 znβn(1,1) +PN n=1 (1 −zn)(xn−µ)2 N−PN n=1 znβn(1,1) (4.45) π(κ+1) =PN n=1 znβn(1,1) N.(4.46) These formulas reduce to the re-estimation formulas for censored data presented in the previous section if the dropping rate is set to π= 0. Then βn(1,1) = 0 and βn(0,1) = 1. As can be seen in Eqs. (4.44), (4.45) and (4.46), not only observable and censored data, but also dropped data contribute to the estimates. 4.3. Optimal Classification Rule for Censored and Dropped Gaussian Data Above ML parameter estimation is done for all possible user locations ℓkand it is done for the measurements of each AP seperately, assuming that data from different APs are independent. The final parameter estimates are denoted by ˆ µk= [ˆµk,1,···,ˆµk,NAP ]Tand ˆ Σk=diag ˆσ2 k,1,···,ˆσ2 k,NAP , where NAP is number of observable APs. Indoor localization can be formulated as a classification problem, where the classes are the positions from which RSSI measurements are taken during the offline training phase. For each position ℓkthe parameters of a Gaussian class-conditional density pY(x|ℓk)of RS- SI measurements are estimated using the EM algorithm of the last sections. During online classification, to estimate the user’s location, an optimal classification rule is developed. Let x=x1,···, xNAP be the vector of the online measurement. First the posterior is calculated as follows p(ℓk|x) = QNAP i=1 p(xi|ℓk)P(ℓk) PK k′=1 QNAP i=1 p(xi|ℓk′)P(ℓk′)(4.47) where Kis the number of offline training locations and xiis the RSSI of i-th AP. P(ℓk)is the prior on the position. Here we have employed the assumption of the independence of the RSSIs of different APs. It is noted that the test data are also subject to censoring. The likelihood p(xi|ℓk)in Eq. (4.47) can be calculated as follows for censored Gaussian data p(xi|ℓk) = Z∞ −∞ p(xi|yi)pY(yi|ℓk)dyi,(4.48) Smartphone Adaptation within the MLLR Framework 40 w(c(k)) i(κ+1) can be obtained as follows w(c(k)) i(κ+1) = R X r=1 N X n=1 γkr(n) σ2 kr,i zn,i I1(λ(κ) kr,i) I0(λ(κ) kr,i)βkr,i(0,1) + (1 −zn,i)xn,i!ξT kr R X r=1 N X n=1 γkr(n) σ2 kr,i (1 −zn,iβkr,i(1,1)) ξkrξT kr!−1 (5.20) 5.2.2. Variance Adaptation Variance adaptation is performed after mean adaptation. It is noted that we assume the RSSI measurements from different APs are independent. Let Σkbe the diagonal covariance matrix of the Gaussian random vector modelling the RSSI readings of the APs at position ℓk. The adapted variance can be calculated as follows [70]: ˆ Σk=BT kH(c(k))Bk(5.21) where H(c(k)) is the transform to be estimated and Bkis the Choleski factor of Σk. Since Σk are diagonal covariance matrices, the re-estimation formula Eq. (5.21) simplifies to ˆσ2 k,i = h(c(k)) iσ2 k,i;i= 1,...,NAP . To estimate H(c(k)), we again employ the EM algorithm. The expected log-likelihood of the complete data, given the observed, has the same form as in mean adaptation. However, here λk,i ={πk,i, h(c(k)) i, θk,i}with h(c(k)) iis the i-th component of the main diagonal of H(c(k)). Using all in the expected log-likelihood function Eq. (5.12), we arrive at: Q(λ;λ(κ)) = K X k=1 N X n=1 γk(n) NAP X i=1 (zn,iβk,i(0,1) ln(1 −πk,i) +zn,iβk,i(0,1) I0(λ(κ) k,i )Zc −∞ ln  1 q2πh(c(k)) iσ2 k,i exp −(yn,i −ˆµk,i)2 2h(c(k)) iσ2 k,i ! N(yn,i;λ(κ) k,i )dyn,i + (1 −zn,i) ln(1 −πk,i) + (1 −zn,i) ln  1 q2πh(c(k)) iσ2 k,i exp −(xn,i −ˆµk,i)2 2h(c(k)) iσ2 k,i !  +zn,iβk,i(1,1) ln(πk,i))(5.22) The re-estimation fomula for h(c(k)) iis readily obtained by computing the derivative of Smartphone Adaptation within the MLLR Framework 41 Eq. (5.22) w.r.t. h(c(k)) iand setting it to zero: h(c(k)) i(κ+1) =1 σ2 k,i PN n=1 γk(n) (1 −zn,iβk,i(1,1)) N X n=1 γk(n) ((1 −zn,i)(xn,i −ˆµk,i)2βk,i(0,0) +zn,i I2(λ(κ) k,i ) I0(λ(κ) k,i )−2I1(λ(κ) k,i ) I0(λ(κ) k,i )ˆµk,i + ˆµ2 k,i!βk,i(0,1)).(5.23) Here ˆµk,i are the new updated means which were discussed in the previous subsection. As can be seen in the derivation, beside estimating the adaptation matrices for means Eq. (5.18) and variances Eq. (5.23), the dropping rates of the adaptation data are simultaneously estimated. Eq. (5.18) and Eq. (5.23) show that not only the observable measurements, but also the unobservable ones contribute to the estimate of the adaptation matrices. If neither censoring nor dropping occurs, i.e., πk,i = 0 and zn,i = 0, Eq. (5.18) and Eq. (5.23) reduce to the formulas for mean and variance adaptation presented in [70]. 5.3. Regression Classes In the previous section, the method for aligning training model with adaptation data in the presence of clipped and dropped data has been presented for the general case of multiple transformation matrices. If the relationship between two sets of training data and adaptation data follows only one linear rule as mentioned in [61], only one set of adaptation matrices, i.e., mean adaptation matrix and variance adaptation matrix, for all states is needed. However, according to our data preliminary study, the relationship of RSSI data measured by any pair of devices does not follow only one linear rule. Two possible main factors which influence the relationship between two sets of data are signal strength range and signal frequency. This is obvious because of the fact that the sensitivities of radio signal sensors depend on the strength and the frequency of the measured radio signal. To relax the linearity assumption, a clustering approach has been applied to seperate the user positions in multiple regression classes that share the same linear relationship, however different from the other clusters. However, a critical issue is how to define the regression classes, within which the parameters of all states, i.e., fingerprint positions, change in the same manner. Since no prior knowledge of which models should share the same transformation rule is available, a criterion for finding regression classes must be found. Various criteria could be used for clustering, for example, the positions with similar RSSI PDFs could share the same linear rule. The similarity of the PDFs can be measured by computing the likelihood of observing the training data collected at one position given the training models of the others. Using these criteria, means and variances of the training models are employed. However, from our observation, the states, i.e., positions, which transform in the same manner have similar RSSI mean values, irrespective of the actual identity of the APs. As a result, those positions should be assigned into the same cluster. We thus apply a k-means clustering on the mean vectors µk Smartphone Adaptation within the MLLR Framework 42 of the probability density function of all location ℓkto obtain Cclusters. For each cluster, transformation matrices are estimated as presented in the previous section. 6. Hidden Markov Model for Indoor User Tracking 6.1. Motivation In chapter 4, we have discussed the method to estimate the parameters of censored and dropped Gaussian data and the optimal classification rule for localizing indoor users. The proposed algorithms belong to the most common approach of indoor positioning which is the fingerprinting based approach. Fingerprinting techniques are able to produce a stable positioning accuracy over time and movement distance of the user. However, there is another approach which can also be applied successfully in indoor positioning, namely the Dead Reckoning technique. This technique uses data from inertial sensors to locate the user position. However, as discussed in chapter 2, the Dead Reckoning technique is often applied in a combination with other positioning techniques since it can only produce precise position estimates in a short period of time and small distance. For long term use and large distance, the position error is unreasonably large because the errors are accumulated over time and distance. As discussed, the stability of fingerprinting based indoor positioning techniques make them become the common candidates to accompany with DR techniques in indoor positioning. Various combinations have been proposed to keep the advantages and to reduce the disadvantages of those two approaches. Sensor fusion can be achieved by a Kalman filter [31], particle filter [32] or with the use of a Hidden Markov Model (HMM) [29, 30]. Among those possible solutions, we decided to employ the HMM in our system because it gives us the possibility to employ the knowledge of possible walking path which may improve the positioning accuracy. 6.2. Hidden Markov Model for Indoor User Tracking In this section, the method to employ a HMM to fuse RSSI measurements and step detection information for position estimation is presented. First, it should be recalled that the standard HMM can be depicted as in Fig. 6.1. Here, stis the hidden state variable and xtis the observation at time t. It is noted that the hidden states cannot be observed, however, the most likely sequence of states can be inferred from the sequence of observations generated by HMM. The HMM is determined by three sets of parameters: the emission probabilities, the state transition probabilities and the initial state probabilities. Fig. 6.2 shows the unfolding of a state diagram over time, called trellis diagram. Here it is assumed that the HMM has K= 4 different states, i.e., K= 4 different values 43 Hidden Markov Model for Indoor User Tracking 44 st−1stst+1 xt−1xtxt+1 Figure 6.1.: Standard Hidden Markov Model: stis the hidden state variable and xtis the observation at time t of st:st∈ {1,2,3,4}. Each column in the graph (corresponding to a time instant) contains nodes representing the states of the HMM. Each node in the graph has connections to at least one node at an earlier and one node at a later time. The connections represent the possible transitions between states. The transitions from one state at one time instant to the possible states at the later time instant have probabilities which sum up to 1, i.e., PK j=1 Pi,j = 1, where Pi,j =P(st+1=j|st=i), assuming the transition probabilities are independent of time, given any i= 1,...,K. Without any prior knowledge, the transition probabilities from one state to the next can be assumed to follow a uniform distribution. 1 1 1 2 2 2 3 3 3 4 4 4 P1,1 P1,2 P1,3 Figure 6.2.: Trellis diagram with four different states st∈ {1,2,3,4}and Pi,j are the transition probabilities from state iat one time instant to state jat the later time instant The remaining parameters of a HMM are the emission probability distributions. For each HMM state, the emission probability distribution p(x|st)gives the likelihood of observing x at state st. Those emission PDFs are supposed to be different amongst states. For positioning purpose, we identify stwith the position of the user at time t: if st=kthen the user is at the location ℓkat time t, where ℓkis a two-dimensional vector describing the user’s location. In WiFi fingerprinting based indoor positioning system, the RSSI training models and the online RSSI measurements can be considered as emission probabilities and observation sequence, respectively. However, with step detection information, there is another set of observations, i.e., move- Hidden Markov Model for Indoor User Tracking 45 ment observations. To incorporate the step detection information, the emission probability distributions for this kind of observation need to be defined. Intuitively step detection information gives us the information about moving from one state to another state, i.e., the state transitions, instead of the information about a single state. So the solution here is to associate an observation with a state transition rather than a state. This ends up with a modification to the HMM to combine RSSI and step detection observations as depicted in Fig. 6.3. In the following, this model is called modified HMM. st−1stst+1 xt−1xtxt+1 vtvt+1 Figure 6.3.: Modified HMM for position estimation based on the fusion of RSSI and movement vector observations xand v, respectively. Here, vtis the two-dimensional movement vector which denotes the traversed route from st−1to st.vtis obtained from inertial sensor measurements. Its computation will be discussed in section 6.4. In the modified HMM, RSSI measurements and movement vectors are taken as observations attached to the hidden states and the state transitions, respectively. The transition probabilities between the states are chosen to reflect which positions are accessible from a given state within one measurement interval. Only transitions to those positions which are accessible within one measurement interval will be given a probability greater than zero. With the Samsung Android smartphones we have at hand, the measurement interval is roughly 1,5 s which is the required time to update a new WiFi scan of the devices. 6.3. Forward Algorithm With the employment of a HMM, the estimation of the user position can then be carried out either by the Forward algorithm, the Viterbi algorithm or the Forward-Backward algorithm. While Forward algorithm computes the probability of being in a certain state by gathering the probabilities over all possible predecessor states, the Viterbi algorithm considers only the most probable predecessor. The Forward-Backward algorithm is only of academic interest because the induced latency is not acceptable for an online positioning system. In the following, the Forward algorithm which aims to estimate the most probable state at time instant t given the history of observations was chosen to decode the user position. Let x1:t= [x1,...,xt]be the sequence of RSSI measurements up to time t. The step detection information is gathered in the sequence v1:t= [v1,...,vt]. Our goal is to compute P(st=k|v1:t,x1:t), i.e., the probability of being in the k-th state for all possible user positions ℓk, k = 1,...,K, given all RSSI values and step detection vectors measured sofar. Using Bayes’ rule, the probability can be expressed as follows: P(st=k|v1:t,x1:t) = p(st=k, v1:t,x1:t) p(v1:t,x1:t) ∝p(st=k, v1:t,x1:t) =: αt(k),(6.1) Hidden Markov Model for Indoor User Tracking 46 where the so-called Forward variable αt(k)is the probability of being at time tin state k, while having observed the sequence of x1:tand v1:t. Using marginal probability and chain rule, the forward variable can be written as follows αt(k) = X i p(st=k, st−1=i, v1:t,x1:t) =X i p(vt|st=k, st−1=i, v1:t−1,x1:t−1,xt) ·p(xt|st=k, st−1=i, v1:t−1,x1:t−1) ·P(st=k|st−1=i, v1:t−1,x1:t−1) ·p(st−1=i, v1:t−1,x1:t−1).(6.2) Applying the properties of the HMM, which are depicted in the graphical model of Fig. 6.3, i.e., xtis independent of all other variables if stis given, stis independent of x1,...,xt−1 and v1,...,vt−1if st−1is given, and vtis independent of everything if st−1and stare given, and assuming the step detection and RSSI information to be statistically independent of each other given the user location, the above formula simplifies to αt(k) = X i p(vt|st=k, st−1=i)·p(xt|st=k) ·P(st=k|st−1=i)·p(st−1=i, v1:t−1,x1:t−1) |{z } =αt−1(i) (6.3) which is a recursion of the forward variable. The final location estimate ˆ ℓis obtained by the weighted average over the set Pof the most likely positions: ˆ ℓ=1 P k∈P αt(k)X k∈P αt(k)·ℓk.(6.4) Equation (6.3) shows how the different knowledge sources are combined. The transition probabilities P(st=k|st−1=i)are nonzero only for those locations ℓkthat can be reached from position ℓiwithin one time step. The choice of the transition probabilities thus encodes our knowledge about the floor plan. The term p(xt|st=k)is the likelihood of the RSSI measurement xt, assuming that the user’s position is ℓk. The computation of emission PDF estimation and likelihood calculation were discussed in Chapter 4. The movement information gathered from the step detection is captured by the term p(vt|st=k, st−1=i), i.e., the likelihood of observing the movement vector vtwhen moving from position ℓito ℓk. It is noted that this derivation requires that RSSI measurements and the movement information are obtained at the same rate. In fact, the rate of user’s steps is often higher than RSSI measurements. To synchronize these 2sources of data, movement vector vtis the accumulated results of the step information during 1RSSI measurement interval. 6.4. Movement Vector Estimation In this section, the method to infer the step detection information from data measured by acceleration sensor, gyroscope sensor and magnetic sensor is summarized. More detailed information about movement vector estimation can be found in [71]. Hidden Markov Model for Indoor User Tracking 47 The system to estimate the movement vectors is depicted in Fig. 6.4. a m g gyr kak ka′k − G ka′kLP Rdψ b ψ Steps ψ Step length Lstep v1:t k·k Lowpass filter Orientation Angular velocity Kalman filter Movement calculation Step detection Figure 6.4.: Step detection and position estimation system overview Here, ais the 3-dimentional acceleration vector from the accelerometer of the smartphone, G= 9,81 m/s2is the gravity constant, which will be used to perform the step detection procedure. m,gand gyr are the magnetometer data, gravity information and gyroscope data, respectively, will be used to estimate the movement heading of the user. The gravity information vector gis the 3-dimentional vector which contains the force of gravity along three axes of the smartphone. Ris the rotation matrix, ψis the angle between movement heading and the magnetic north direction, and dψ denotes the changes of ψover time. The detail information about these data and how to retrieve them by an Android application can be found on the website for android developers [72]. Movement vector vtis computed by accumulating the estimated movement of all detected steps within one RSSI measurement interval. 6.4.1. Step Detection The step detection procedure is based on the measured acceleration data from the smartphone. The obtained samples are stored in 3-dimensional vectors a= [ax, ay, az]Twhere x, y, z are the axes in the coordinate system. This is defined in the relation with the screen of the smartphone as shown in Fig. 6.5(a) (a) Device coordinates system (b) World coordinate system Figure 6.5.: Device coordinate system (a) and world coordinate system (b) Hidden Markov Model for Indoor User Tracking 48 The characteristics of the user movement is captured in the acceleration data as depicted in Fig. 6.6. As can be seen in the example the acceleration component along the z-axis clearly shows the up and down motion caused by walking. It is noted that the measured data are always influenced by the force of gravity. Therefore, to obtain the real acceleration data, the contribution of the force of gravity must be eliminated. One can use this azfor step detection purpose if that is the stable feature of the obtained data. However, this is just an example of acceleration data when the smartphone was held in the way that its screen is in parallel with the earth surface. If the smartphone, however, is held in an arbitrary attitude, the observed data will not follow the example data in the figure any more. Because of that reason, instead of using data from one specific component, the Euclidean norm of the acceleration data from all three components is used to perform step detection, see Eq. (6.5). Fig. 6.7 shows the calculated norm values after substracting the force of gravity constant G=9.81. kak=qa2 x+a2 y+a2 z(6.5) az Samples ay Samples ax Samples Raw Acceleration Data [m/s2] 128 192 256 320 128 192 256 320 128 192 256 320 −10 0 10 20 30 −20 0 20 −20 0 20 Figure 6.6.: Example of raw acceleration data if the smartphone is held with display in parallel to earth surface when walking To reduce the variation caused by sensor errors, a 5-Hz lowpass filter is applied, as is suggested in [57]. Fig. 6.8 shows the output signal of the lowpass filter. Here, a Butterworth low-pass of order 20, with a cutoff frequency of 0.2·fn= 5Hz, where fncorresponds to a half of the sampling rate (50Hz/2) is used. In the literature, two common approaches which have been employed in many research to detect steps using acceleration data are peak detection and zero crossing count. In [57] a peak detection with a minimum threshold is used to count the steps. This has the disadvantage that sometimes the maximum is split into two local maxima and thus an additional step is detected. The other method for detecting steps calculates the zero crossing rate of the Hidden Markov Model for Indoor User Tracking 49 ka′k=kak−G Samples 128 192 256 320 −10 −5 0 5 10 Figure 6.7.: Norm of acceleration data after substracting the force of gravity constant G=9.81 acceleration data [59], which in case of noisy data may also result in additionally detected steps. The proposed approach combines the strengths of peak detection and zero crossing count approaches by counting the crossing rate of acceleration data over a certain threshold. Our method is then able to cope with multiple local maxima and noisy data in the measured acceleration data. Fig. 6.9 illustrates the threshold-crossing approach. The step counter is only increased if the threshold (magenta line) is larger than the safety threshold (red line) to avoid an erroneous counting when the user is not moving but there is still some fluctuation of the observed acceleration data. The threshold is not a fixed value through the whole process, instead it is calculated for every single data block. Here a block size of 64 samples is chosen regarding the sampling rate of acceleration sensor and processing speed. First, for each data block, the peak mean is determined by calculating the mean of all the detected peaks. Second, the threshold is computed by multiplying the peak mean with an experimentally determined weighting factor. In our experiment, the weighting factor of 0.6produces the best results. Recently, another method, namely autocorrelation method, for detecting steps has been proposed in [32]. This method employs the fact that the acceleration data exhibits a very repetitive pattern. The advantage of the autocorrelation method compared to peak detection and zero crossing count methods is that it is able to eliminate the hand gestures when the user is not moving, transition from sitting to standing and vice versa, and so on. According to our experimental results, the autocorrelation method outperforms the other two traditional methods. Therefore, at the current state of our project, this method is employed in our system for detecting steps. However, this section aims to present our proposed approach to detect steps which has be- Experimental Results on Indoor Positioning 56 Table 7.1.: Mean and standard deviation of APs at 2 positions for artifical data experiment AP index AP1 AP2 AP3 AP4 AP5 µ1,i -102 -103 -97 -89 -95 µ2,i -105 -100 -99 -86 -101 σ1,i 4.8 4.9 5.0 5.2 5.1 σ2,i 5.0 4.8 4.8 5.4 5.0 Table 7.2.: Classification error rate on artificial censored data Method Error rate (%) Plain trng + recog 30.7 EM trng + plain recog 26.9 EM trng + censored recog 22.5 3-strongest APs 35.1 1-nearest neighbor 36.8 •Plain training (trng) + recognition (recog): ML parameter estimation is carried out assuming normally distributed, uncensored data. Also recognition is performed disregarding any censoring. •EM trng + plain recog: ML parameter estimation in training accounts for the censored data using the proposed EM algorithm, while the presence of censored data is still disregarded in recognition. •EM trng + censored recog: Training with the proposed EM algorithm and recognition employing eq. (4.52). •3-strongest APs: Select three strongest APs of each location in the training phase, then apply EM trng + censored recog. •1-nearest neighbor classification rule. Table 7.2 clearly shows the superiority of the schemes which are aware of the censoring. Considering the presence of censored data in training improved the error rate from 30.7% to 26.9%, and a further improvement to 22.5% is obtained by accounting for censored data also in recognition. We can also see the important role of weak APs for the recognition accuracy: using only the three strongest APs raises the error rate to 35.1%.1-nearest neighbor performs worst with the error rate of 36.8%. Field Data In order to evaluate the performance of the proposed EM algorithm on real world data, real WiFi RSSI measurements were gathered on a floor of an office building consisting of 10 office rooms and a long aisle having an overall size of 12 m by 30 m (see Fig. 7.1). RSSI values were taken at 25 different positions, roughly evenly distributed, resulting in an average distance of 2,7 m between two locations. Two measurement campaigns were carried out using a smartphone, with 100 measurements taken per position per campaign. Data of the Experimental Results on Indoor Positioning 57 first measurement campaign served the purpose of training and the second one was for testing purposes. For the training data set, the percentage of uncensored observations, averaged over all APs which were observable at each location, was found to be 36.7%. Figure 7.1.: Floor plan of the area where field data has been conducted. To compare the proposed approach to a state-of-the-art system, the algorithm from [52] was implemented. There, for each location ℓkthe probability distributions of the 10 strongest APs are determined during the training phase and compared to those of the online measurements employing the Bhattacharyya coefficient. A 3-nearest neighbor rule is then applied to decide on the user location. Furthermore, we compared with the well-known system, RADAR [18], where classification is performed with a 3-nearest neighbor rule, employing the Euclidian distance. When applying the proposed algorithm, 5online measurements were used for each estimate, the final location estimate was then computed using the 3most likely positions, see Eq. (4.53). It is noted that 5online measurements were also employed in the implementation of the algorithm of [52]. Fig. 7.2 shows the cumulative distribution function (CDF) of the error as a function of the distance for each method. The CDF is defined as the probability that the positioning error ǫ is lower than a certain distance d: CDFǫ(d) = P(ǫ≤d)d≥0.(7.1) The results in Fig. 7.2 show that the proposed method outperforms the other two, especially for the 40% error quantile. Note, also, that the computational cost of the proposed method during the online phase is smaller than those of computing the Bhattacharyya distances between probability distributions [52] or the nearest-neighbor based [18] methods. After [18] After [52] Proposed algorithm Distance d[m] CDFǫ(d) 0 2 46 8 10 0 0.2 0.4 0.6 0.8 1 Figure 7.2.: CDF of the positioning error for different systems. Experimental Results on Indoor Positioning 58 Table 7.3.: Classification error rate on artificial censored and dropped data Method Error rate (%) EM trng + censored recog 29.7 Adv. EM trng + censored & dropped recog 25.7 7.1.2. EM algorithm for censored and dropped data In the following, experimental results showing the effectiveness of EM algorithm when being aware of dropped data in addition to censored data are presented. For convenience, let us call the EM algorithm for parameter estimation of censored and dropped data as the advanced EM algorithm. Artificial Data The artificial data were generated according to the parameters which were used in Section 7.1.1, however, the means of all APs at all positions are increased by 10 dBm to show more impact of dropped data on the estimated parameters and consequently on classification results. The generated Gaussian data were first censored and then dropped with the dropping rate of 20%. Since the effectiveness of EM algorithm for censored data has been proved in the previous sections, this section aims to compare the classification results between two proposed EM algorithms to show the importance of the awareness of dropped data besides censored data. Table 7.3 shows the classification results on censored and dropped data using the proposed EM algorithms. It can be seen that the parameter estimation formulas derived in Section 4.2.3 are able to cope with an unknown drop-out rate, which leads to the improvement in the classification error rate from 29.7% to 25.7% if indeed drop-outs occur. The reason is if all unobservable data are considered as censored data, the parameters estimated by the EM algorithm which is aware of censored data only, are inaccurate, i.e., the estimated means are biased to the left of the true means and the estimated variances are higher than the true variances. These inaccurate estimated parameters cause the degradation of the classification performance. Field Data To demonstrate the effectiveness of the advanced EM algorithm, an experiment on real field data has been done. This experiment was done by using the same data sets as used in the experiment of the EM algorithm for censored Gaussian data. Fig. 7.3 shows the experimental results on real field data using the two proposed EM algorithms. As can be seen, being aware of dropped data improves the positioning accuracy, though the improvement is moderate. The reasons might be either the amount of training data, 100 samples per position, are not sufficient for the advanced EM algorithm since it tries to estimate more parameters, or the true dropping rates are low. Since the data sets were gathered in a real indoor environment, true dropping rates are unknown. As our study on the results, the estimated means and variances from advanced EM are better than the EM algorithm awaring of censored data only, and the estimated dropping rate averaged over all APs and all locations was approximately Experimental Results on Indoor Positioning 59 0.37, we concluded that the limited improvement of positioning accuracy is because of the not very accurate estimated dropping rates. As presented in section 4.3, dropping rate plays an important role in likelihood calculation for the unobservable data, so the inaccurate estimated dropping rates limit the improvement in classification results. In the bad case of low accuracy of the estimated dropping rates, the positioning accuracy might even decrease. EM Censored Adv. EM CDFǫ(d) Distance d[m] 01 2 3 4 56 7 8910 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 7.3.: Comparing the experimental results on real data using 2proposed EM algorithms 7.2. Smartphone Adaptation This section investigates the impact of the proposed adaptation approach on the accuracy of fingerprinting based indoor positioning. Some variations of the adaptation approach within the MLLR framework such as mean and variance adaptation, mean adaptation only with full or diagonal adaptation matrix have been implemented. In addition, the impact of the number of adaptation data and of the number of regression classes on the performance of adaptation, and consequently on the positioning accuracy, are also investigated. 7.2.1. Classification on Artificial Data In this set of experiments, artificial data were generated to assess the impact of the drop-out rate and the amount of adaptation data on the classification performance. 1000 samples of training data were generated for each of K= 6 positions and for NAP = 2 access points. A single affine transformation of the means of the training data is used for each location according to ˆ µk=Aµk+b, where A= [0.9,0; 0,0.9] and b= [−3,2]T. Adaptation and test data were then generated by sampling from the transformed Gaussians with means ˆ µkand variances as those of the training data. While the amount of adaptation data was gradually increased from 1% to 5, 10, 50 and 100% of the number of training samples, the number of test samples was fixed at 100 observations per position. The value of the drop-out rate is set to 0% (no drop-outs, π= 0) and 20% (π= 0.2). Although the effectiveness of being aware of dropped data on classification results was demonstrated in section 7.1.2, for the convenience of evaluating the effectiveness of the adaptation approach, classification results employing advanced EM with the current set of parameters are still presented. Table 7.4 discusses the impact of the drop-out rate without Experimental Results on Indoor Positioning 60 Table 7.4.: Classification results (% positions correctly classified) when training and test data are generated from the same parameter set Aware of dropping yes no π= 0 86 86 π= 0.275 64 Table 7.5.: Adaptation performance (% positions correctly classified) on artifical data as a function of amount of adaptation data No. Pos with adapt. data π0% 1% 5% 10% 50% 100% 3 0 68 81 85 85 86 86 0.2 56 70 73 74 75 75 6 0 68 85 86 86 86 86 0.2 56 72 75 75 75 75 adaptation. Here, the ’adaptation’ data were used for a completely new training from scratch for each of the 6positions. 1000 samples of adaptation data were used in this experiment. As can be seen, if no dropping occurs, the two EM algorithms produce the same classification result. Once the dropping occurs, the advanced EM algorithm outperforms the other. Table 7.5 shows the effectiveness of applying the proposed adaptation algorithm. As can be seen, the results when only 3locations have adaptation data is comparable to the case when adaptation data are available for all 6locations. The reason is even if there are only 3locations with adaptation data, still the PDFs of all locations are well adapted, since they are from a single regression class. However, it is also noticeable that when the number of adaptation data per position is small, i.e., less than 10% of the training data, classification results using 6positions with adaptation data are about 1to 4% better than those when adaptation data are only available for 3positions. The column with 0% of adaptation data shows the results when no adaptation is performed. Comparing with the results of Table 7.4 it is clear that the classification results after adaptation approach those of retraining, except when there is only 1% of adaptation data, even if adaptation data are only available for 3out of 6positions. For the case of only 1% of adaptation data available, there is still a dramatic improvement in classification results when adaptation is employed compared to the case if no adaptation is performed. 7.2.2. Classification on Field Data To examine the effectiveness of the adaptation approach on real world data, measurements were collected on 3floors of an office building with roughly 30 rooms (lecture halls, office and laboratory rooms), where each floor has an overall size of 35 m by 35 m, RSSI values were taken at 60 different positions with an average distance of 5.0m between 2positions. 200 measurements were taken per position with 2different smartphones at each position. Data of the first smartphone were used to estimate the training models while data from the second smartphone were divided into 2sets at each position. The first was used for adaptation (0, 5, 25, or 75 samples) or retraining (150 samples) from scratch and the second was for testing (50 samples). For the measured data, the estimated dropping rate averaged over all Experimental Results on Indoor Positioning 61 Table 7.6.: RMS positioning error (in [m]) as a function of the amount of positions having adaptation data and the amount of adatation data Condition Adaptation Method No. Pos with Amount of µ&σ2,µonly, µonly, adapt. data adapt. data Afull Afull Adiag 0 6.15 15 5 5.08 4.76 3.93 25 4.60 4.41 3.93 75 4.62 4.52 3.93 30 5 3.86 3.96 3.84 25 3.82 3.96 3.85 75 3.76 3.91 3.85 60 5 3.74 3.93 3.84 25 3.67 3.87 3.83 75 3.65 3.87 3.84 retrain 2.43 APs and all positions was approximately 0.3. For the adaptation procedure, the estimated training means were sorted at each position in descending order, and only the 8strongest APs were used to estimate the adaptation matrices since the contribution of the remaining APs to the likelihood was negligible. The adaptation matrices are then used to calculate the adapted parameters of the 8strongest APs at all positions being in the same regression class. Table 7.6 shows the dependency of the positioning accuracy on the number of locations for which adaptation data are available and the amount of available adaptation data at each position. The results in Table 7.6 are the average of the root mean square (RMS) positioning error of 50 experiments. In each experiment, test data, adaptation data and the positions with adaptation data are randomly selected. As expected, the more positions there are with adaptation data and the more adaptation samples per position, the better the positioning results. For all considered adaptation methods, improvements in positioning accuracy are obtained, even when very few adaptation data are available. However, mean and variance adaptation with a completely filled matrix Adelivers only the best results if there are sufficiently many adaptation data, while the use of only mean adaptation with a diagonal Ais superior if only few adaptation data are available, since fewer parameters need to be estimated. From the presented results, indoor positioning system developers could have an idea of which and how parameters can be efficiently adapted given the available adaptation data. However, the positioning accuracy when all 60 positions have adaptation data is still well below the accuracy achievable, when training and test data are collected by the same smartphone, which is 2.43 m. The reason could be one transformation rule could not describe well the relationship between training and adaptation data for all positions which actually follows a nonlinear rule as discussed in chapter 5. To relax the linearity assumption between training and adaptation data, multiple regression classes are employed. For estimating the regression classes, the means of the 8strongest APs at each position are sorted in descending order, the resulting 8-dimensional vectors are Experimental Results on Indoor Positioning 62 Table 7.7.: Effect of number of clusters on RMS positioning error No. of clusters 1 2 3 RMS pos. error [m] 3.76 3.48 3.25 clustered using k-means, as discussed in section 5.3. Table 7.7 shows the positioning results when doing clustering and estimating different adaptation matrices for each cluster, assuming that in each cluster 50% of positions have adaptation data with 75 adaptation samples per position. It has to be noted that the optimal number of regression classes depends on the amount of available adaptation data. The more adaptation data the more transformation matrices can be reliably estimated. In our setup the best results were achieved with 3clusters, which led to a reduction of the RMS positioning error from 3.76 m to 3.25 m. 7.3. HMM for Indoor User Tracking This section evaluates the effectiveness of the proposed algorithms discussed in Chapter 6 for an indoor positioning problem both on artificially generated data and real field data, especially the impact of the introduction of pseudo states on the classification/positioning results. 7.3.1. Artificial Data Fig. 7.4 shows the HMM states with the locations of the regular and pseudo states marked with red and green circles, respectively. 530 535 540 545 550 555 560 565 985 990 995 1000 1005 1010 1015 x[m] y[m] Figure 7.4.: The modified HMM states, i.e., the allowable user positions, and the transitions between them. Red and green circles indicate regular and pseudo states, respectively. The RSSI measurements of 15 randomly placed APs for training the Gaussian emission PDF for the regular HMM states are generated artificially as follows: The signal strength fol- Experimental Results on Indoor Positioning 63 Table 7.8.: Mean positioning error on artificial data Method Mean error [m] RSSI only [74] 1.74 RSSI + step det. [73] 1.37 RSSI + step det. + pseudo states [75] 1.02 lows a large-scale log-normal fading model with an additional zero mean Gaussian random variable with standard deviation σL= 5 to model small-scale fading. At each position, we generated a set of 300 RSSI measurements as the training data, then estimated the parameters of RSSI distribution of each AP at each position as described in Chapter 4. The step detection information is modeled as a bivariate Gaussian with the mean µi,j = ℓj−ℓiand a diagonal covariance matrix Σvwith entries Σv[1,1] = Σv[2,2] = 0,25 m2. This value has been determined in offline experiments. In the experiment, pseudo states were introduced in the way that the Euclidean distance between neighboring states is not more than at most 0,75 m (which is close to the measured average step length within the experimental data), while the distance between two neighboring regular states is about 3-5m except for some special areas such as the stairs. The total number of pseudo states is 125, which has to be compared to the number of 81 regular states. In the experiments we assume that a user cannot move faster than 3 m/s. User movement is simulated by a random walk on the HMM graph. Since movement vectors and RSSI measurements are generated every 1,5 s, only a limited number Uof HMM states can be reached from any given state st−1=i(on the average U=15 states) and the corresponding transition probability P(st=j|st−1=i)is set to 1 U. For all states outside this neighborhood the corresponding transition probabilities are set to zero. Table 7.8 presents the mean positioning error in meters, averaged over 100 experiments, where each experiment corresponds to a different random walk on the HMM grid of length of about 200 m. We compare the performance of the proposed algorithm with our earlier work of [73], which also fused RSSI and step detection information, however without the introduction of pseudo states. Using pseudo states improves the mean positioning error from 1,37 m to 1,02 m. If step detection information is neglected and user positioning relies only on RSSI information, a mean positioning error of 1,74 m is obtained. 7.3.2. Field Data The proposed approach is also evaluated with the field data recorded in the office building depicted in Fig. 7.1. In the training phase, 100 RSSI measurements per position were collected and the RSSI distribution were estimated as described in Chapter 4. In the online testing phase, the smartphone user randomly went through the whole floor area to collect the test data, i.e., WiFi data and inertial sensor data. Two different trajectories were recorded, each trajectory consists of approximately 140 test potitions. The position estimate was performed every 1,5 s using the proposed approach. The step detection information is modeled in the same fashion as in the experiments using artificial data. According to the obtained data and experimental results, it seemed that the step detection information is more reliable than the RSSI information, as a consequence, a heuristic Experimental Results on Indoor Positioning 64 weighting factor λwas introduced in the calculation of the forward variable as follows αt(j) = [p(ot|st=j)]λX i [p(vt|st=j, st−1=i)](1−λ)P(st=j|st−1=i)·αt−1(i).(7.2) For the determination of λa jackknife procedure was employed: The data of 1out of the 2 trajectories was used for the estimation of λ, whereas tests were conducted on the held-out data. This was repeated 2times, every time, one trajectory was used to estimate λ. With our collected data, the estimated value was always λ≈0.003. The very small value of λ can be explained as follows: since the likelihood of RSSI data is much smaller than the likelihood of step detection information, this value of λwould help to avoid the problem that the contribution of the likelihood of step detection information to the calculated forward variable is dominated by the likelihood of RSSI data. The proposed method was compared to our ealier work: first, using RSSI only for user positioning [74], and second, using HMM model as described in Chapter 6, however without the introduction of pseudo states [73]. For both trajectories the experimental results showed that the new approach outperforms the others, especially for the 90% error quantile, in terms of the CDF of the positioning error Eq. (7.1). Although the test area is limited, the experimental results in Fig. 7.5 indicate that the proposed approach is significantly better than the other approaches. RSSI only RSSI + Step Det. RSSI + Step Det. + Pseudo States CDFǫ(d) Distance d[m] 0 1 2345 6 78 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 7.5.: CDF of the positioning error for different systems. Average over 2test trajectories. 8. Server Based Indoor Navigation System Navigation Server Clients Figure 8.1.: Server based Indoor Navigation System. Together with developing theoretical algorithms for the enhancement of positioning accuracy, a real indoor navigation system has been developed for over a period of three years with the contributions from many employees and students at the Department of Communications Engineering, University of Paderborn. The system is the result of our attempt to bring the theoretical research results into reality. Before building the indoor navigation system, one needs to answer the question: Which criteria must the system meet? To our understanding, the most important criteria to evaluate a real system are costs, ease of deployment, stability, security, scalability and maintainability. The better the criteria are fullfilled, the better the system is. System architecture has the most important impact on those criteria. As a consequence, to achieve the required criteria the best, different system architectures have been analysed. There exist two common architectures of positioning systems: First, smartphone based system in which all data are stored and processes are performed on the smartphone itself, e.g., the solutions from WIFARER [77], 65 Server Based Indoor Navigation System 72 markerid lon lat room layer ref markerlocations connectionid markerid linkedmarker distance highway access markerconnections Table 8.1.: Map information tables in the database: storage of the geographical information of fingerprint positions, and the information about the direct connection between any pair of positions. we decided to combine all the seperate OSM map files into one global OSM map file. The navigation database is then built up by importing the information from the global map file to the PostGIS database. The rendered map tiles are stored on the server and delivered to the clients for map displaying purposes once the server receives map requests from clients. Map information stored in the navigation database is the information of fingerprinting positions and the connections among those. This information is gathered from the standard tables by a C++ program and stored in two new tables (see Table 8.1 for an illustration). The bold fields are the so-called unique keys or master keys of the tables. Table markerlocations contains all the information about any fingerprinting position such as longitude, latitude, room number, layer (floor level) and ref (corridor, toilet, ect.). These fingerprinting positions are the positions where training data will be collected. It is very convenient to distribute the training positions during the map editing procedure since one can choose the critical and interesting points on the maps to mark them as fingerprinting points. It is also easy to manage the density of the fingerprint grid. Table markerconnections contains the information about the direct connection between any pair of positions, e.g., ID of current considered position (markerid), ID of the connected position (linkedmarker), spatial distance between these positions (distance), type of connection (highway), i.e., footway or bicycle, and accessibility information (access) to indicate whether the connection is blocked or if there is free access. The information in Table markerconnections is needed for the route estimation procedure and the employment of the HMM which was discussed in Chapter 6. 8.2.2. RSSI data RSSI data are divided into two types: raw data (RSSI measurements) and processed data (the estimated training models). Table 8.2 shows the tables in the database which store the measurement data. There are several reasons for storing the raw data. For example, raw data can be used for future research, i.e., if new algorithms are used to estimate the training model from measured data, or for comparison purposes, i.e., evaluating the positioning accuracy of other algorithms on the same data set. By storing the raw data in the way as described in Table 8.2, it is possible to re-generate the original measurements. Tables scanresults, measurements and hardwareinformation are used together to store the RSSI measurements, the position, ID and orientation of the device, and the time stamp of the measu- Server Based Indoor Navigation System 73 scanid measurementid mac rssi scanresults measurementid lon lat layer hwid orientation measurementtime measurements hwid modelid uniqueid hardwareinformation Table 8.2.: Measurement tables in the database: contain the RSSI measurements, the position, ID and orientation of the device which has been used to collect data, and the time stamp of the measurements parameterid markerid macid numberofmeasurement numberofobservation mean variance rate sumofobservations sumofsquare updatedtime parameters macid bssid macaddresses Table 8.3.: Training model information tables in the database: contain the estimated parameters and the sufficient statistics information. rements. The information about the smartphone hardware is stored to support the adaptation procedure as described in chapter 5. In addition to storing the raw data, the RSSI database also contains the tables to store the information of the estimated training models as well as the sufficient statistics information as depicted in Table 8.3. Table parameters stores the estimated parameters as well as the sufficient statistics information. The estimated parameters are the means, variances and rate (dropping rate) of the observed RSSI data of the AP, specified by macid, at a location, specified by markerid. Sufficient statistics parameters such as number of measurements, number of observations, sum of observations and sum of observations squared, are used for incrementally updating the training model to reduce the estimation time using the EM algorithm. Obviously, estimation time is not a serious problem when the amount of training data is small. However, training data are collected over the course of time and model estimation from scratch would take much more time than applying the incremental update method. For convenience of AP searching, the table macaddresses is created to store the MAC addresses of all the observed access points. Server Based Indoor Navigation System 74 8.3. Shared Memory As mentioned above, shared memory is used to represent the database in order to speed up the data reading procedure. The structure of the data stored in shared memory is similar to the structure used in the PostGIS database as presented in the previous section. Seperate shared memory areas are created to store the information of fingerprinting positions, connections, MAC addresses and trained models. In addition, to speed up the localization procedure, other areas are created to store the information of the pre-computed positions at which a given MAC address was observed in the training data. The reason is that in fingerprinting based positioning, online measured data are compared with the training data in the database to come up with the decision of the user position estimate. Obviously the bigger the database becomes, the more computation time is needed to complete the comparison procedure. To get rid of this problem, instead of computing the similarity between the online measurement with the whole database, the localization module has to compute it at the possible positions only, where at least one of the APs, which is present in the online measurement, was observed in the training data. 8.4. Overview of Smartphone Application This section presents an overview of the smartphone application. Although the developed application is able to support many tasks, e.g., data gathering, localization, navigation, social features (peer group finding) and 3-D representation of the routing path, in the following, only the structure which is related to the main features, which are localization and navigation, is presented. Fig. 8.5 gives an overview of the reduced version of the smartphone application architecture as a class diagram. It should be noted that in Fig. 8.5, the details of each class, i.e., variables and methods, are not presented since this would extremely extend the size of the graph. In the following, the role of each class in the application and the relation of classes will be summarized. As can be seen in Fig. 8.5, besides the regular classes, there are several special types of classes in an Android application such as activity, service, application which operate differently [72] as summarized below: •Activity: An activity is an application component that provides a single window with a user interface that interacts with the users in order to do an action. An application may consist of several activities to handle different tasks. An activitiy can start another activity. Once a new activity starts, the previous activity is stopped. The system preserves the previous activity in a “last in, first out” stack. Therefore, when the user is done with the current activity and presses the “Back” button of the smartphone, the previous activity will resume. •Service: A service is an application component that can operate in the background and does not provide a user interface. A service can be started by another application component, i.e., activity, and keeps running in the background even if the user switches to another application. Several types of service implementation can be used to manage Server Based Indoor Navigation System 75 Figure 8.5.: Class diagram of Java application architecture the lifecycle of a service that the service may stop itself when it completes the task or its binding component stops or an activity stops it. •Application: An application is a base class to maintain global application state. In our application, the main activity is “IndoorNavMap” which is launched when the user starts the application. All other activities, i.e., WiFi scanning, setting, etc., can be activated from this activity, see the manifest file of the application in Appendix A.3.1 for information of all activities. Fig 8.6 shows a snapshot of the application when it is launched by the users. In the following discussion, the relation between the classes in Fig. 8.5 and the role of each class are summarized: •“IndoorNavMap” uses the setting information and initialized information stored in “MapScrollApplication” and the appropriate methods implemented in “CellMapSurfa- Server Based Indoor Navigation System 76 Figure 8.6.: Java application showing floor plan information and user position (red dot). ceView” to display map and other information such as current user position, navigation path, etc. . •“MapScrollApplication”, the application class, is implemented to maintain all the global states and the data which are needed for map displaying, offline localization, navigation, and so on. For map drawing, positioning and navigation, “MapScrollApplication” stores the metadata of map data, i.e., objects of “Map” and “MapTile” classes, and the information for localization and navigation purposes, i.e., “BinaryModel“, ”MacObject“, ”NodeObject“ and ”RoutingPathElement“, as descibed in section 8.2. •“CellMapSurfaceView” extends the “SurfaceView” class provided by Android APIs, and contains all methods to manage map drawing and user interaction. •“IServerConnection” is the interface declaring all the methods which are implemented in “HttpServerConnection” to support the communication procedure between the application and the server. These methods can be called by “IndoorNavMap” for data downloading or by “LocalizationService” in case of doing localization in online mode. •As all the messages are in the pre-defined XML format, methods for writing/parsing data to/from XML format messages are implemented in “XmlWriter” and “XmlReader” for communicating with the server. •“LocalizationService” is the service for handling the tasks related to the localization procedure. This service is started by “IndoorNavMap” and calls the methods implemented in “LocalizationManager” to perform localization. This service is stopped automatically by the system once the user terminates the application. •“IndoorNavMapBroadcastReveiver” and “LocalizationServiceBroadcastReceiver” are implemented to support the data exchange between “LocalizationService” and “IndoorNavMap”. These two classes extend the “BroadcastReceiver” provided by the Android APIs. •“WiFiScanner” is responsible for WiFi data acquisition for localization procedure. Server Based Indoor Navigation System 77 •“RouteCalculate” contains all the methods related to navigation path calculation. These methods are activated by “IndoorNavMap” via “MapScrollApplication” where the information of the calculated route is stored. By doing this, if the user switches to any other activities or applications, once user gets back to the “IndoorNavMap” activity, the navigation path will be displayed without re-calculation. More details about the navigation procedure will be presented in section 8.5.3. 8.5. System Operation This section presents the operation of the main features of our indoor positioning system from a user’s view, as well as the flow sequence of the smartphone application for each feature. The developed indoor navigation system is able to support several services such as gathering of training data, localization, navigation, some initial versions of social features such as peer group finding, and 3-D representation of the routing path. The main features, i.e., gathering of training data, localization, and navigation, are discussed in detail, a short introduction about the other features is given at the end of the section. 8.5.1. Data Gathering This feature is responsible for gathering training data. To do this, the smartphone application performs WiFi scans and writes the measured data to a “SetFingerprintRequest”. It should be noted that on the smartphones, at the moment, this feature is only available when the application runs in the development mode in order to avoid problems in case someone tries to harm the database. In future, for extending the system towards online learning, this feature could be performed in the background of the smartphone application during the usage by users. The development mode is activated by the developers by choosing the option “Activate developer view”, as shown in Fig. 8.7. Figure 8.7.: Java application showing setting options. The development view of the smartphone application is illustrated in Fig 8.8. In the figure, the red circles indicate the positions at which RSSI training data have already been collected. Server Based Indoor Navigation System 78 By doing so, the developers can then choose the fingerprint positions which have not been trained to collect RSSI data. The red “+” sign on the map shows the position where WiFi data are going to be collected. The developers can modify the position of the red “+” sign by moving the map. To collect the RSSI data, since each position needs a sufficient amount of measurements to accurately estimate training models, an activity named “WLANDataAssemblerActivity” was developed to support this requirement. Figure 8.8.: Java application showing the development mode: the red circles are the positions where training data were already collected, the red “+” sign on the map shows the position where WiFi data are going to be collected. Once the developer starts the data gathering procedure, a window appears which allows the developer to verify the information of the current position and enter additional information, i.e., name of the position, amount of scans and scan interval, see Fig. 8.9(a). Depending on the sampling rate of the WiFi sensor of the test smartphone, the scan interval need to be large enough to avoid duplication of the measurement data, We observed the sufficient scan interval of about 1,5 s for Samsung smartphones, while for a Sony device it is roughly 5 s. Once the information is validated and entered, the measurement process will be started by simply clicking the “Start scan” button. Fig. 8.9(b) shows the screenshot when the smartphone performs a WiFi scan, where the number of measurements and scanned data are shown. The measured data are then written into an XML file which is stored on the local storage of smartphone once the scanning procedure is completed. The XML file contains the measured RSSIs, MAC addresses of the observed APs, the smartphone information, as well as the information of the location where the measurements are taken. It should be noted that each XML file contains the training data at one position only. After the measurement campaign, all XML files are imported into the database on the server using a C++ module named setfingerprint. This module parses the request and accesses the PostGIS database to insert the measured data. For each file, setfingerprint produces a response message to inform whether the data are successfully inserted to the database or not. Appendix A.4.1 shows an example of the “SetFingerprintRequest” and response messages. Our system is also able to support the online mode of the measurement campaign, i.e., the measured data are sent directly to the server after each measurement via a wireless network, Server Based Indoor Navigation System 79 (a) Setting for WiFi scan (b) WiFi scanning Figure 8.9.: WiFi training data gathering if we use the setfingerprint module as a FastCGI module. This option provides the possibility of gathering data during the usage of users for an online learning approach. However, at the moment, no algorithm has been developed to measure the reliability of a measurement. Therefore, we temporarily disable this option to avoid accidental measurements containing wrong information which may harm the database. The indoor environment and WLAN network are subject to change over time, this requires regular database maintenance (update) to keep the positioning accuracy. It is infeasible to retrain the database after every change of the environment or the network infrastructure. The reason is collecting training data is very time consuming, especially for the large deployment area. Therefore, a method for automatically updating database, i.e., the online learning approach, could be a solution. The sequence diagram of Fig. 8.10 shows how the smartphone application performs WiFi scanning in offline mode, i.e., the scanned data are written in an XML file and stored in the local storage of the smartphone. 8.5.2. Localization For localization, users can choose either the offline localization mode or the online localization mode in the settings of the smartphone application . In offline mode, position estimation is performed locally on the smartphone using the database which is stored on its local storage without any need of internet connection. In online localization mode, which requires network connection for data exchange, position estimation is performed on the server. Figure 8.11 shows the options for localization, i.e., operation mode and the data source. The sequence diagrams shown in Fig. 8.12 and Fig. 8.13 present the communication among the components of the smartphone application to perform localization in offline mode and online mode using WiFi information, respectively. To perform the localization using WiFi data, smartphones periodically log the WiFi data which are the MAC addresses of the observed APs and their measured RSSIs. An estimator is then used to process the scanned data to determine the user positions In offline localiza- Server Based Indoor Navigation System 80 Figure 8.10.: Sequence diagram showing the data gathering procedure. tion, the estimator is implemented as a method in the smartphone Java application which uses the local database to compare with observed data. Several issues may arise in offline localization, for example, the local database or the positioning algorithms might be out of date, as discussed in the beginning of this chapter. To solve this problem, the local RSSI database can be reloaded from server regularly, i.e., once a week, automatically or manually if network connection is available. The problem with positioning algorithms is unsolvable unless the user updates the application. Within this work, an automatic update procedure for smartphone application has not been developed yet. For online localization, scanned data are sent to the server via the WLAN connection in an XML format. On the server, a FastCGI module named “positionestimate” was developed to parse the request and perform the localization procedure. The position estimation algorithm was implemented using the classification rule which was discussed in Chapter 4, the same algorithm is used for offline localization. This module uses the training models stored in the shared memory to carry out the localization procedure. Sensor fusion for the improvement in positioning accuracy as discussed in Chapter 6, is not implemented on the server side, but on the client side. The reason is that the server must be a stateless server since there would be an explosion of the amount of HMM states which need to be stored if many clients request the localization service at the same time using RSSI Server Based Indoor Navigation System 81 (a) Setting options (b) Localization modes Figure 8.11.: Setting screen showing options for localization. and step detection information. This would dramatically degrade the response time of the server. Once the position estimate is obtained, a localization response is sent from the server to the client which contains the information of the estimated user position. The client parses the response and shows the user location on the smartphone screen, see Fig. 8.6 where the red dot represents the current user position. 8.5.3. Navigation For navigation purpose, the user can input the information of the source and destination positions, i.e., room numbers, and the options for the expected route, i.e., using elevators or stairs. Fig. 8.14 shows an example of a smartphone screen when the user starts the routing function of the smartphone application. Navigation, similar to localization, can be performed in either offline mode or online mode. In online mode, a navigation request with the information inserted by the user is generated by the smartphone application and sent to the remote server. On the server, a FastCGI module named “routeestimate” was developed which is able to parse the request and perform the route calculation. Data for the route calculation process are obtained by reading the shared memory. To perform the path search, several path finding algorithms have been considered such as the well known Dijkstra’s algorithm, the A* algorithm (a variant of Dijkstra’s algorithm) and the jump point algorithm. In navigation, the path cost is simply the geographical distance of any pair of connected positions. While Dijkstra’s algorithm examines all nodes to find the shortest path between the source and the destination, the A* algorithm is trying to examine the nodes which are potentially on the shortest path first to optimize the computational demand. This is the reason why the A* is also called goal-oriented Dijkstra’s algorithm. However, the A* algorithm needs heuristic weights to guarantee the solution is the shortest path. The jump point algorithm tries to optimize A* in case of uniform-cost grids, which is normally not suitable for an indoor environment since distributing the indoor fingerprinting 9. Conclusions Within this thesis, the techniques to improve the accuracy of WiFi fingerprinting based indoor positioning are presented. The accuracy of indoor positioning employing WLAN information can be enhanced by a statistical approach which is able to account for the variation of the measured RSSIs in the indoor environment. The positioning accuracy depends on two factors: training models and classification rule. As our careful study on WiFi data in indoor environment showed, RSSI measurements suffer from two problems namely censoring and dropping. As discussed in chapter 2 and chapter 4, these two problems of the indoor WiFi data have not been addressed in any previous research. Therefore, within this work, methods for estimating the parameters of the training model and classification in the presence of censored and dropped data were proposed. For parameter estimation, an EM algorithm was proposed in chapter 4 which efficiently copes with the censoring and dropping problem. The proposed EM algorithm for censored Gaussian data was proved to be a virtually biasfree and efficient estimator. Improvements in positioning accuracy are demonstrated both on artificially generated data and in real field data experiments compared to some other approaches. The experiments presented in chapter 7 showed the superiority of the statistical approach compared to a deterministic approach. It is noted that the computational demand for positioning using a parametric statistical approach is less than the other mentioned approaches. Another problem that has been considered in this work is the mismatch between the measured data of the training device and the test devices as discussed in chapter 5 which leads to a serious reduction of the positioning accuracy. In the literature, we found only one solution which tried to address this problem using Least Squares approach where the authors assumed a linear relationship between the RSSI readings of different devices. However, what we observed is that this relation is not linear which renders the LS approach inappropriate. An effective method to cope with this problem must be developed to make WiFi signal based indoor positioning realistic. Therefore, we proposed a method for aligning the training model with the properties of the test devices while relaxing the linearity assumption in chapter 5. The proposed aligning method called “smartphone adaptation” was developed within the MLLR framework which is a very well known and successful technique for speaker adaptation in automatic speech recognition. It has to be noted that the proposed adaptation method is able to cope with censored and dropped data in the adaptation data, resulting in reliable adapted models. During the adaptation procedure, the dropping rate of the adaptation data is also estimated which will be used in the classification procedure. Experimental results presented in chapter 7 demonstrated the effectiveness of doing adaptation in indoor positioning. Assuming the RSSI readings from different devices follow one linear relationship, applying the proposed approach showed a big improvement in positioning accuracy, however, still well far below the achievable accuracy. Employing clustering approach before 88 Conclusions 89 doing adaptation relaxed the linearity assumption resulting in assuming piecewise linearity. As a result, better position accuracy was obtained. Improvements in positioning accuracy can be obtained by fusing the information from different sources such as WiFi information and inertial sensor information. These two kinds of information are obtainable on most modern smartphones. Inertial navigation which utilizes the data from the built-in sensors of the smartphones is able to produce precise position estimation in a short term (time and distance) only. Unfortunately the precise positioning results cannot be maintained for a longer period of time due to error accumulation over time and distance, resulting in unreliable position estimates. WiFi fingerprinting based positioning, on the other hand, can provide a stable positioning accuracy without the error accumulation problem. As discussed in chapter 2, these two positioning techniques can be combined in order to produce better positioning results compared to using any individual approach alone. Therefore, chapter 6 presented a modified HMM for data fusion of WiFi data and inertial sensor data and utilizing the possible walking path of the user to come up with position estimates. More accurate positioning results were obtained by employing HMM in comparison with using WiFi information or inertial sensor information alone. Furthermore, a method to reduce the quantization error caused by the coarse grid of trained positions was proposed in chapter 6. By introducing pseudo states to the HMM inbetween the regular states and synthesizing the emission PDFs of the pseudo states from those of neighboring regular states, we obtained a dense grid of states without additional training effort. Experimental results showed that employing the extended HMM improved the positioning accuracy. In addition to the theoretical research, a real server based indoor positioning and navigation system was developed as presented in chapter 8. At the moment the system is able to provide the localization and navigation services for the students of the University of Paderborn. The system employs the lighttpd web server and the FastCGI protocol which allows to handle thousands of connections in parallel on the server. In addition, the FastCGI modules, which are written in C++, run in seperate processes which make the system stable, and further, easy to scale and easy to maintain. A Java application for Android smartphones was developed which is able to run as a standalone positioning system or as a client in the system. As discussed in chapter 8, developing two possible operation modes for the smartphone application solves the problem of loss of internet connection, unsynchronization of map data and fingerprinting training model, and so on. As a result, the system seems to satisfy the requirements of low cost, high robustness and simplicity in deployment, scaling and maintenance. Outlook An interesting direction for future research is how to improve/update the radio map with the measured data reported by the users during the online phase. To do that, a possible solution is to develop a semi-supervised online learning system where the already built database is continuously updated during the system operation using the measurements reported by users. This kind of system can automatically handle the changes of the indoor environment or WLAN infrastructure without re-training the database from scratch after a certain amount of time. This also helps to improve the accuracy of the training model since more training data are available. The challenge is how to determine the position where the measurement is Conclusions 90 taken. The REDPIN system believes in the information that is reported from any users which does not seem to be fault tolerant since the database can be corrupted on purpose or unintentionally through wrong user positions. Simultaneous Localization And Mapping (SLAM) is a common approach in the robotics community which mainly relies on inertial sensor information to track the position of a mobile robot and build the map simultaneously. This approach can be used to determine the position at which a specific RSSI measurement is collected. However, as discussed before, this technique suffers from the error accumulation over time and distance and, as a consequence, estimated positions are not reliable. Map matching can be formulated as a fingerprinting based localization problem where the user position is assigned to a route segment based on the online observation and training data. Therefore, the combination of SLAM and map matching approach would be a relevent solution to estimate the position of the user where RSSI measurement is collected since the accumulated error in inertial navigation are corrected globally by RSSI data and floor plan information. Our proposed EM algorithms presented in chapter 4 are efficient for estimating the parameters of censored and dropped Gaussian data. However, we employed the empirical fixed clipping threshold for data measured by all devices. This might degrade the parameter estimation performance if the clipping threshold of different devices are not identical. Therefore, estimating clipping threshold from training data should be done before parameter estimation procedure to ensure the precision of the estimated parameters. As a result, a method for clipping threshold estimation is needed. Employing a Gaussian mixture model estimation to estimate the RSSI distribution instead of the assumption of a single Gaussian would be an interesting try. This method might help to improve the precision of signal strength distribution estimation, since it is able to capture all the modes in the training data distribution. However, since censoring and dropping are severe in WiFi data, methods for dealing with censored and dropped data during the parameter estimation procedure would need to be taken into consideration. A. Appendix A.1. Derivation of EM Algorithm A.1.1. Computation of I0 I0(θ(κ)) = Zc −∞ Ny;θ(κ)dy =Zc −∞ 1 √2πσ(κ)exp −(y−µ(κ))2 2(σ2)(κ)dy(A.1) Let t=y−µ(κ) √2σ(κ)and change the variable of the integral, we arrive at: I0(θ(κ)) = 1 √πZc−µ(κ) √2σ(κ) −∞ exp(−t2)dt(A.2) Since exp(−t2)is an even function, the limits of the integral in Eq. (A.2) can be modified by changing their signs and swapping lower and upper limits, then I0(θ(κ))can be easily obtained by using complementary error function I0(θ(κ)) = 1 2erfc −c−µ(κ) √2σ(κ)(A.3) A.1.2. Computation of I1 I1(θ(κ)) = Zc −∞ yNy;θ(κ)dy =Zc −∞ y1 √2πσ(κ)exp −(y−µ(κ))2 2(σ2)(κ)dy(A.4) Employing the integration by parts rule, let (u=y dv=1 √2πσ(κ)exp −(y−µ(κ))2 2(σ2)(κ)dy⇒(du=dy v=1 2erf y−µ(κ) √2σ(κ),(A.5) then we arrive at I1(θ(κ)) = y1 2erf y−µ(κ) √2σ(κ) c −∞ −1 2Zc −∞ erf y−µ(κ) √2σ(κ)dy |{z } =A1 (A.6) 91 Appendix 92 Again integration by parts rule is used for computing the remaining integal A1as follows, to simplify the integral, let t=y−µ(κ) √2σ(κ),A1can be written as A1=√2σ(κ)Zc−µ(κ) √2σ(κ) −∞ erf (t)dt(A.7) Let u=erf (t) dv=dt⇒du=2 √πexp (−t2)dt v=t,(A.8) then we arrive at A1=√2σ(κ) terf (t) c−µ(κ) √2σ(κ) −∞ −1 √πZc−µ(κ) √2σ(κ) −∞ 2texp −t2dt  =√2σ(κ)terf (t) + 1 √πexp −t2 c−µ(κ) √2σ(κ) −∞ (A.9) Using A1in Eq. (A.6), I1(θ(κ))is readily obtained: I1(θ(κ)) = µ(κ)I0(θ(κ))−1 √2πσ(κ)exp −c−µ(κ) √2σ(κ)2!(A.10) A.1.3. Computation of I2 I2(θ(κ)) = Zc −∞ y2Ny;θ(κ)dy =Zc −∞ y21 √2πσ(κ)exp −(y−µ(κ))2 2(σ2)(κ)dy(A.11) Employing the integration by parts rule, let (u=y2 dv=1 √2πσ(κ)exp −(y−µ(κ))2 2(σ2)(κ)dy⇒(du= 2ydy v=1 2erf y−µ(κ) √2σ(κ),(A.12) then we arrive at I2(θ(κ)) = y21 2erf y−µ(κ) √2σ c −∞ −Zc −∞ yerf y−µ(κ) √2σ(κ)dy |{z } A2 (A.13) To compute A2, integration by parts is again applied, let (u1=y dv1=erf y−µ(κ) √2σ(κ)dy⇒(du1=dy v1=Rerf y−µ(κ) √2σ(κ)dy,(A.14) Appendix 93 For calculating v1, let t=y−µ(κ) √2σ(κ) v1=√2σ(κ)Zerf (t)dt(A.15) and applying integration by parts u2=erf (t) dv2=dt⇒du2=2 √πexp (−t2)dt v2=t,(A.16) we arrive at v1=√2σ(κ)t.erf (t)−1 √πZ2texp −t2dt =√2σ(κ)t.erf (t) + 1 √πexp −t2.(A.17) Using t=y−µ(κ) √2σ(κ)in v1then we obtain v1=√2σ(κ)y−µ(κ) √2σ(κ)erf y−µ(κ) √2σ(κ)+1 √πexp −(y−µ(κ) √2σ(κ))2 = (y−µ(κ))erf y−µ(κ) √2σ(κ)+σ(κ)r2 πexp −y−µ(κ) √2σ(κ)2!.(A.18) Now, A2can be calculated as follows A2=u1v1 c −∞ −Zc −∞ v1du1 =y((y−µ(κ))erf y−µ(κ) √2σ(κ)+σ(κ)r2 πexp −y−µ(κ) √2σ(κ)2!) c −∞ −Zc −∞ ((y−µ(κ))erf y−µ(κ) √2σ(κ)+σ(κ)r2 πexp −y−µ(κ) √2σ(κ)2!)dy =y((y−µ(κ))erf y−µ(κ) √2σ(κ)+σ(κ)r2 πexp −y−µ(κ) √2σ(κ)2!) c −∞ −Zc −∞ y.erf y−µ(κ) √2σ(κ)dy |{z } A2 +µ(κ)Zc −∞ erf y−µ(κ) √2σ(κ)dy |{z } v1c −∞ −2σ2(κ)Zc −∞ 1 √2πσ(κ)exp −y−µ(κ) √2σ(κ)2!dy |{z } I0(θ(κ)) (A.19) Appendix 94 Simply using the limits, A2is readily obtained: A2=1 2y2erf y−µ(κ) √2σ(κ) c −∞ −1 2(µ(κ)2 +σ2(κ)).erf c−µ(κ) √2σ(κ)+ 1 −c+µ(κ)σ(κ)r2 πexp −c−µ(κ) √2σ(κ)2!) =1 2y2erf y−µ(κ) √2σ(κ) c −∞ −(µ(κ)2 +σ2(κ))I0(θ(κ)) +c+µ(κ).σ(κ)1 √2πexp −c−µ(κ) √2σ(κ)2!(A.20) Plugging A2into Eq. (A.13), the calculation of I2ends up with I2(θ(κ)) = σ2(κ)+µ2(κ)I0(θ(κ))−1 √2πσ(κ)(µ(κ)+c) exp −c−µ(κ) √2σ(κ)2!(A.21) A.1.4. Computation of W Entries ∂ ∂µI0(θ) = ∂ ∂µ Zc −∞ p(y;θ)dy =1 σ2Zc −∞ (y−µ)p(y;θ)dy =1 σ2(I1(θ)−µI0(θ)) (A.22) ∂ ∂µI1(θ) = ∂ ∂µ Zc −∞ y.p(y;θ)dy =1 σ2Zc −∞ y(y−µ)p(y;θ)dy =1 σ2(I2(θ)−µI1(θ)) (A.23) ∂ ∂µI2(θ) = ∂ ∂µ Zc −∞ y2.p(y;θ)dy =1 σ2Zc −∞ y2(y−µ)p(y;θ)dy =1 σ2(I3(θ)−µI2(θ)) (A.24) ∂ ∂σ2I0(θ) = ∂ ∂σ2Zc −∞ p(y;θ)dy =−1 2σ2I0(θ) + 1 2σ4I2(θ)−2I1(θ)µ+I0(θ)µ2(A.25) Appendix 95 ∂ ∂σ2I1(θ) = ∂ ∂σ2Zc −∞ yp(y;θ)dy =−1 2σ2I1(θ) + 1 2σ4I3(θ)−2I2(θ)µ+I1(θ)µ2(A.26) ∂ ∂σ2I2(θ) = ∂ ∂σ2Zc −∞ y2p(y;θ)dy =−1 2σ2I2(θ) + 1 2σ4I4(θ)−2I3(θ)µ+I2(θ)µ2(A.27) where I3(θ)and I4(θ)can be computed by using integration by parts as the computation of I2(θ), however lengthier computation I3(θ) = Zc −∞ y3p(y;θ)dy =µ3σ2+µ2I0(θ)−1 √2πσ[2σ2+µ2+µc +c2] exp −c−µ(κ) √2σ(κ)2!(A.28) I4(θ) = Zc −∞ y4p(y;θ)dy = (3σ4+ 6µ2σ2+µ4)I0(θ) −1 √2πσ[σ2(5µ+ 3c) + µ3+µ2c+µc2+c3] exp −c−µ(κ) √2σ(κ)2!(A.29) A.1.5. Computation of Information Matrix I For convenience, the log-likelihood function and the information matrix are recalled as follows ln p(x;θ) = ln N! M!(N−M)!+ (N−M) ln (I0(θ)) −M 2ln 2πσ2−1 2 M X j=1 xj−µ σ2 .(A.30) I= Eh−∂2 ∂µ2ln p(x;θ)iEh−∂2 ∂µ∂σ2ln p(x;θ)i Eh−∂2 ∂µ∂σ2ln p(x;θ)iEh−∂2 ∂σ2∂σ2ln p(x;θ)i .(A.31) In the following, to shorten the notation, the parameters of the Ijfunctions are removed. Here the derivative computations of the Ijwhich were given in Eq. (A.22) to (A.27) are employed. Since I0is the probability that a measurement is censored, in the following derivation, the expected number of uncensored measurements is determined by E [M] = N(1 −I0), where Nis the total number of measurements. Appendix 96 •Eh−∂2 ∂µ2ln p(x;θ)i: First derivative of log-likelihood function w.r.t. µ ∂ ∂µ ln p(x;θ) = (N−M)1 I0 ∂ ∂µI0−1 2 M X j=1 −2xj−µ σ2 = (N−M)1 I0 1 σ2I1−µI0+1 σ2 M X j=1 (xj−µ)(A.32) Second derivative of log-likelihood function w.r.t. µ: ∂2 ∂µ2ln p(x;θ) = ∂ ∂µ (N−M)1 I0 1 σ2I1−µI0+1 σ2 M X j=1 (xj−µ)! =N−M σ2 ∂ ∂µ I1 I0−µ−M σ2 =N−M σ2 ∂I1 ∂µ I0−I1∂I0 ∂µ I2 0−1!−M σ2 =N−M σ21 σ2(I2−µI1)I0−I11 σ2(I1−µI0) I2 0−1−M σ2 =N−M σ4I2I0−I2 1 I2 0−σ2−M σ2 =N−M σ4I2I0−I2 1 I2 0−N σ2(A.33) Expectation of the second derivative of log-likelihood function w.r.t. µ: E−∂2 ∂µ2ln p(x;θ)=E−N−M σ4I2I0−I2 1 I2 0+N σ2 =−E[N−M]1 σ4I2I0−I2 1 I2 0+N σ2 =−NI0 1 σ4I2I0−I2 1 I2 0+N σ2 =N σ2I2 1−I2I0 I0 + 1(A.34) Appendix 97 •Eh−∂2 ∂µ∂σ2ln p(x;θ)i: ∂2 ∂µ∂σ2ln p(x;θ) = ∂ ∂σ2 (N−M)1 I0 1 σ2I1−µI0+1 σ2 M X j=1 (xj−µ)! = (N−M)∂ ∂σ21 σ2I1 I0−µ−1 σ4 M X j=1 (xj−µ) = (N−M)(−1 σ4I1 I0−µ+1 σ2 ∂I1 ∂σ2I0−I1∂I0 ∂σ2 I2 0!) −1 σ4 M X j=1 (xj−µ)·(A.35) Expectation of the second derivative of log-likelihood function w.r.t. µand σ2: E−∂2 ∂µ∂σ ln p(x;θ)=NI0(1 σ4I1 I0−µ−1 σ2 ∂I1 ∂σ2I0−I1∂I0 ∂σ2 I2 0!) +1 σ4E"M X j=1 (xj−µ)#,(A.36) where the remaining expectation can be computed as follows E"M X j=1 (xj−µ)#=E[M] (E[x]−µ) =N(1 −I0)1 1−I0Z∞ c yp(y;θ)dy −µ =N(1 −I0)1 1−I0 (µ−I1)−µ =N(µI0−I1).(A.37) Using this in Eq. (A.36) we arrive at E−∂2 ∂µ∂σ2ln p(x;θ)=N σ2− ∂I1 ∂σ2I0−I1∂I0 ∂σ2 I0(A.38) •E−∂2 ∂σ2∂σ2ln p(x;θ): First derivative of log-likelihood function w.r.t. σ2 ∂ ∂σ2ln p(x;θ) = (N−M)1 I0 ∂ ∂σ2I0−M 2 1 σ2−1 2−1 σ4M X j=1 (xj−µ)2 = (N−M)−1 2σ2+1 2σ4I2 I0−2µI1 I0 +µ2 −M 2σ2+1 2σ4 M X j=1 (xj−µ)2(A.39) Notations and Symbols General Notations and Functions E[·]. . . . . . . . . . . . . . . . Expectation (·)T. . . . . . . . . . . . . . . . Transpose Q(·). . . . . . . . . . . . . . . Auxiliary function (·)−1. . . . . . . . . . . . . . . Inversion ln(·). . . . . . . . . . . . . . . Natural logarithm function p(x|ℓk). . . . . . . . . . . . . Class conditional probability density function δ(·). . . . . . . . . . . . . . . . Diract delta function erf (·). . . . . . . . . . . . . . Error function erfc (·). . . . . . . . . . . . . Complementary error function N(θ). . . . . . . . . . . . . . Gaussian distribution parameterized by θ p(x;θ). . . . . . . . . . . . . Probability density function parameterized by θ P. . . . . . . . . . . . . . . . . Summation of a sequence Q. . . . . . . . . . . . . . . . . Products of a sequence (.)! . . . . . . . . . . . . . . . . . Factorial of a non negative integer Fundamentals of Indoor Positioning and State of Research R. . . . . . . . . . . . . . . . . . Mathematical expression of radio map in fingerprinting based techniques ℓk. . . . . . . . . . . . . . . . . . Vector consists of coordinates of the k-th position Fk. . . . . . . . . . . . . . . . . Fingerprint at the k-th position Xk. . . . . . . . . . . . . . . . . Set of measurement vectors at the k-th position xk,n ................ Then-th measurement vector at the k-th position D(o,xk,n). . . . . . . . . . Euclidian distance between online sample oand training sample xk,n θk. . . . . . . . . . . . . . . . . . Set of the parameters of the Gaussian describing the signal strength distribution at the k-th position µk. . . . . . . . . . . . . . . . . Mean vector of the Gaussian describing the signal strength distribution at the k-th position Σk. . . . . . . . . . . . . . . . . Covariance matrix of the Gaussian describing the signal strength distribution at the k-th position 104 Notations and Symbols 105 Parameter Estimation of Censored and Dropped Gaussian Data c. . . . . . . . . . . . . . . . . . . Clipping threshold Θ. . . . . . . . . . . . . . . . . . Set of parameters of a GMM πk. . . . . . . . . . . . . . . . . Mixing weight of the k-th component of a GMM θk. . . . . . . . . . . . . . . . . . Set of the parameters of the k-th component of a GMM X. . . . . . . . . . . . . . . . . . Set of observable data Z. . . . . . . . . . . . . . . . . . Set of hidden variables y. . . . . . . . . . . . . . . . . . Set of unobservable, non-censored, possibly dropped data yn................. Then-th measurement, scalar value x. . . . . . . . . . . . . . . . . . Set of observable data θ. . . . . . . . . . . . . . . . . . . Set of the parameters of a univariate single Gaussian µ. . . . . . . . . . . . . . . . . . Mean of a univariate single Gaussian σ. . . . . . . . . . . . . . . . . . Variance of a univariate single Gaussian θ(κ). . . . . . . . . . . . . . . . Set of the estimated parameters of a univariate single Gaussian after the κ-th iteration of an EM algorithm Ijθ(κ)............ Thej-th moment of the truncated part of a censored Gaussian, defined in Eq. (4.8) zn. . . . . . . . . . . . . . . . . . Realization of the binary random variable Znindicate whether the n-th measurement is censored (zn= 1) or not (zn= 0) N. . . . . . . . . . . . . . . . . . Number of measurements M. . . . . . . . . . . . . . . . . Number of observable measurements ˜µ(κ). . . . . . . . . . . . . . . . Difference between the estimated mean after the κ-th EM iteration and the true mean (˜σ2)(κ). . . . . . . . . . . . . Difference between the estimated variance after the κ-th EM iteration and the true variance I. . . . . . . . . . . . . . . . . . . Information matrix dn. . . . . . . . . . . . . . . . . Hidden variables indicate whether the n-th measurement is dropped (dn= 1) or not (dn= 0) π. . . . . . . . . . . . . . . . . . Dropping rate, defined as π=P(dn= 1) βn(d, z). . . . . . . . . . . . Probability that dn= 0 or dn= 1 given zn= 0 or zn= 1 Smartphone Adaptation within MLLR Framework µk. . . . . . . . . . . . . . . . . Mean vector consists of means of all APs at the k-th position ξk. . . . . . . . . . . . . . . . . Extended mean vector consists of means of all APs at the k-th position and a offset term factor ˆ µk. . . . . . . . . . . . . . . . . Adapted mean vector at the k-th position NAP . . . . . . . . . . . . . . . Total number of observable APs c(k). . . . . . . . . . . . . . . . Regression class c(k)where the Gaussian describing k-th position belongs to C. . . . . . . . . . . . . . . . . . Total number of regression classes W(c(k)) . . . . . . . . . . . . . Mean transformation matrix to be applied to the regression class c(k) γk(n). . . . . . . . . . . . . . The posterior probability that RSSI measurement vector xnis from position ℓk Notations and Symbols 106 Y. . . . . . . . . . . . . . . . . . Set of measurement vectors yn................. Then-th measurement vector consists of RSSI from NAP APs yn,i . . . . . . . . . . . . . . . . The non-censored, possibly dropped measured data of the n-th measurement from the i-th AP D. . . . . . . . . . . . . . . . . . Set of hidden variable vectors dn................. Then-th hidden variable vector consists of random variables, each indicates whether the measured data from corresponding AP is dropped or not dn,i . . . . . . . . . . . . . . . . Random variable indicates whether the n-th measurement from the i-th AP is dropped or not X. . . . . . . . . . . . . . . . . . Set of observable, censored, possibly dropped measurement vectors xn................. Then-th measurement vector consists of observable, censored, possibly dropped measured data from NAP APs xn,i . . . . . . . . . . . . . . . . The censored, possibly dropped measured data of the n-th measurement from the i-AP w(c(k)) i. . . . . . . . . . . . . . The row vector which is the i-th row of W(c(k)) λ. . . . . . . . . . . . . . . . . . Set of parameters, i.e., dropping rate, mean and variance adaptation matrices, to be estimated for adaptation purpose λk,i . . . . . . . . . . . . . . . . Set of parameters, i.e., dropping rate, mean and variance adaptation matrices, to be estimated of the i-th AP at the k-th position πk,i . . . . . . . . . . . . . . . . Dropping rate of the adaptation data of the i-th AP at the k-th position µk,i . . . . . . . . . . . . . . . . Mean of the Gaussian describing the RSSI distribution of the i-th AP at the k-th position σk,i . . . . . . . . . . . . . . . . Variance of the Gaussian describing the RSSI distribution of the i-th AP at the k-th position zn,i . . . . . . . . . . . . . . . . . Random variable indicates whether the n-th measurement from the i-th AP is observed or not βk,n,i(d, z). . . . . . . . . . Probability that dn,i = 0 or dn,i = 1 given zn,i = 0 or zn,i = 1 at the k-th position ˆ Σk. . . . . . . . . . . . . . . . . Adapted covariance matrix of the Gaussian describing the RSSI readings of the APs at the k-th position Bk. . . . . . . . . . . . . . . . . The Choleski factor of Σk H(c(k)) . . . . . . . . . . . . . Variance transformation matrix to be applied to the regression class c(k) h(c(k)) i.............. Thei-th component of the main diagonal of H(c(k)) Hidden Markov Model for Indoor User Tracking st. . . . . . . . . . . . . . . . . . Hidden state variable at time t xt. . . . . . . . . . . . . . . . . . RSSI observation at time t Pi,j . . . . . . . . . . . . . . . . Transition probability from the i-th to the j-th state vt. . . . . . . . . . . . . . . . . . Two-dimensional movement vector computed from inertial sensor data Notations and Symbols 107 αt(k). . . . . . . . . . . . . . . Forward variable: The probability that the user is at time tin state kgiven the observed sequence of o1:tand v1:t a.................. 3-dimentional acceleration data vector ax. . . . . . . . . . . . . . . . . Acceleration data along x-axis ay. . . . . . . . . . . . . . . . . Acceleration data along y-axis az. . . . . . . . . . . . . . . . . . Acceleration data along z-axis G. . . . . . . . . . . . . . . . . . Gravity constant m. . . . . . . . . . . . . . . . . Magnetor data vector g. . . . . . . . . . . . . . . . . . Gravity data vector gyr . . . . . . . . . . . . . . . . Gyroscope data vector R. . . . . . . . . . . . . . . . . . Rotation matrix ψ. . . . . . . . . . . . . . . . . . Azimuth angle α(n). . . . . . . . . . . . . . . State vector that contains the absolute Azimuth angle and its derivative F(n). . . . . . . . . . . . . . . Transition matrix which indicates a transition of the system from time nto n+ 1 in Kalman filter ν1(n). . . . . . . . . . . . . . White Gaussian system noise in Kalman filter ν2(n). . . . . . . . . . . . . . White Gaussian measurement noise in Kalman filter G(n). . . . . . . . . . . . . . . Matrix that describes the influences of the system noise on the state vector in Kalman filter Q1. . . . . . . . . . . . . . . . . Covariance matrix of the white Gaussian system noise in Kalman filter Q2. . . . . . . . . . . . . . . . . Covariance matrix of the white Gaussian measurement noise in Kalman filter H(n). . . . . . . . . . . . . . . Measurement matrix that describes the influences of the measurement on the state vector in Kalman filter Lstep . . . . . . . . . . . . . . . Step length Experimental Results on Indoor Positioning λ. . . . . . . . . . . . . . . . . . Weighting factor between inertial sensor information and WiFi information in the calculation of forward variable List of figures 2.1. Proximity-based positioning: the estimated position of the mobile user Miis the position of the base station BSjif Midetects only signal from BSj, or the observed signal strength of BSjis the strongest among the detected singals 6 2.2. Cicular lateration based positioning: the estimated position of the mobile user Mis determined based on the estimated distances between the mobile user and the reference points. At least three reference points are needed to calculate the user position. . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3. Angulation based positioning: the estimated position of the mobile user M is determined based on the estimated angle αibetween the mobile user and the reference points. At least two reference points are needed to calculate the user position with the condition that the user position does not lie on the line connecting the reference points. . . . . . . . . . . . . . . . . . . . . . . . 11 2.4. Fingerprinting based positioning: the filled circles indicates the reference (anchor) points, while the triangles show base station locations. The user position is the position of the anchor point whose training data best match the online measurement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.5. Dead reckoning: the user position is estimated based on the estimated movement vectors assuming that the starting point is given. . . . . . . . . . . . . 15 4.1. Histogram of real field data illustrates censoring and dropping problem of WiFi data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.2. Censoring problem: left figure: PDF from where yis drawn; right figure: PDF from where the observation xis drawn. . . . . . . . . . . . . . . . . . 23 4.3. Theoretical number of EM iterations κrequired to reduce the estimation error to 10−4of its initial value as a function of (c−µ)/σ. Initial values µ(0),(σ2)(0) have been set to the ML estimates of µ, σ2computed from the unclipped observations only. . . . . . . . . . . . . . . . . . . . . . . . . . 29 4.4. Comparison of CRLB for mean and variance with MSE obtained from simulation for σ2= 25 and N= 1000. . . . . . . . . . . . . . . . . . . . . . . 30 4.5. Measurement model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 6.1. Standard Hidden Markov Model: stis the hidden state variable and xtis the observation at time t. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 6.2. Trellis diagram with four different states st∈ {1,2,3,4}and Pi,j are the transition probabilities from state iat one time instant to state jat the later time instant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 108 List of figures 109 6.3. Modified HMM for position estimation based on the fusion of RSSI and movement vector observations xand v, respectively. . . . . . . . . . . . . 45 6.4. Step detection and position estimation system overview . . . . . . . . . . . 47 6.5. Device coordinate system (a) and world coordinate system (b) . . . . . . . 47 6.6. Example of raw acceleration data if the smartphone is held with display in parallel to earth surface when walking . . . . . . . . . . . . . . . . . . . . 48 6.7. Norm of acceleration data after substracting the force of gravity constant G=9.81 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 6.8. ka′kafter lowpass filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 6.9. Step detection procedure: In each data block, the number of steps is determined based on the crossing rate of the acceleration data over a threshold (magenta lines). The safety threshold is the minimum threshold value that the amplitude of noisy acceleration data cannot exceed. . . . . . . . . . . . 51 6.10. Introduction of pseudo states: The red circles are the regular states with training data, the green circle is the pseudo state introduced inbetween the regular states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 7.1. Floor plan of the area where field data has been conducted. . . . . . . . . . 57 7.2. CDF of the positioning error for different systems. . . . . . . . . . . . . . . 57 7.3. Comparing the experimental results on real data using 2proposed EM algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 7.4. The modified HMM states, i.e., the allowable user positions, and the transitions between them. Red and green circles indicate regular and pseudo states, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 7.5. CDF of the positioning error for different systems. Average over 2test trajectories. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 8.1. Server based Indoor Navigation System. . . . . . . . . . . . . . . . . . . . 65 8.2. Server architecture of the indoor navigation system. . . . . . . . . . . . . . 69 8.3. Examples of original map and edited map . . . . . . . . . . . . . . . . . . 71 8.4. Map editing with JOSM . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 8.5. Class diagram of Java application architecture . . . . . . . . . . . . . . . . 75 8.6. Java application showing floor plan information and user position (red dot). 76 8.7. Java application showing setting options. . . . . . . . . . . . . . . . . . . . 77 8.8. Java application showing the development mode: the red circles are the positions where training data were already collected, the red “+” sign on the map shows the position where WiFi data are going to be collected. . . . . . . . . 78 8.9. WiFi training data gathering . . . . . . . . . . . . . . . . . . . . . . . . . 79 8.10. Sequence diagram showing the data gathering procedure. . . . . . . . . . . 80 8.11. Setting screen showing options for localization. . . . . . . . . . . . . . . . 81 8.12. Sequence diagram presents the offline localization procedure. . . . . . . . . 82 8.13. Sequence diagram presents the online localization procedure. . . . . . . . . 83 8.14. Java application showing the routing options. . . . . . . . . . . . . . . . . 84 8.15. Sequence flow presents the offline navigation procedure. . . . . . . . . . . 85 8.16. Java application showing the indoor routing path. . . . . . . . . . . . . . . 86 8.17. Java application showing the indoor and outdoor routing path. . . . . . . . 86 List of figures 110 8.18. Peer group finding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 8.19. 3-D representation of the routing path . . . . . . . . . . . . . . . . . . . . 87 List of tables 7.1. Mean and standard deviation of APs at 2 positions for artifical data experiment 56 7.2. Classification error rate on artificial censored data . . . . . . . . . . . . . . 56 7.3. Classification error rate on artificial censored and dropped data . . . . . . . 58 7.4. Classification results (% positions correctly classified) when training and test data are generated from the same parameter set . . . . . . . . . . . . . . . 60 7.5. Adaptation performance (% positions correctly classified) on artifical data as a function of amount of adaptation data . . . . . . . . . . . . . . . . . . . 60 7.6. RMS positioning error (in [m]) as a function of the amount of positions having adaptation data and the amount of adatation data . . . . . . . . . . . . 61 7.7. Effect of number of clusters on RMS positioning error . . . . . . . . . . . 62 7.8. Mean positioning error on artificial data . . . . . . . . . . . . . . . . . . . 63 8.1. Map information tables in the database: storage of the geographical information of fingerprint positions, and the information about the direct connection between any pair of positions. . . . . . . . . . . . . . . . . . . . . . . . . 72 8.2. Measurement tables in the database: contain the RSSI measurements, the position, ID and orientation of the device which has been used to collect data, and the time stamp of the measurements . . . . . . . . . . . . . . . . 73 8.3. Training model information tables in the database: contain the estimated parameters and the sufficient statistics information. . . . . . . . . . . . . . . 73 111 References [1] R. Mautz, “Indoor Positioning Technologies,” in Habilitation Thesis submitted to ETH Zurich, Zurich, February 2012. [2] “Indoor Positioning and Navigation,” http://www.computer.org/portal/web/computingnow/archive/june2013. [3] A. Kushki, K. N. Plataniotis, and A. N. Venetsapopoulos, WLAN Positioning Systems, Cambridge, New York, United States, 2012. [4] M. D. Rodríguez, J. Favela, E. A. Martínez, and M. A. Muñoz, “Location-aware access to hospital information and services,” IEEE Transactions on Information Technology in Biomedicine, vol. 8, no. 4, pp. 448–455, 2004. [5] “GPS website,” http://www.gps.gov/. [6] P. Enge and P. Misra, “Special Issue on Global Positioning System,” in Proceedings of the IEEE, January 1999, vol. 87, pp. 3–15. [7] A. Varshavsky, M. Y. Chen, E. Lara, J. Froehlich, D. Haehnel, J. Hightower, A. LaMarca, F. Potter, T. Sohn, K. Tang, and I. Smith, “Are GSM phones THE solution for localization,” in Proceedings of the 7th IEEE Workshop on Mobile Computing Systems and Applications, Orcas Island, WA, August 2006. [8] R. Haeb-Umbach and S. Peschke, “A Novel Similarity Measure for positioning Cellular Phones by a Comparison With a Database of Signal Power Levels,” IEEE Transactions on Vehicular Technology, pp. 368–372, 2007. [9] M. Ibrahim and M. Youssef, “CellSense: A Probabilistic RSSI-based GSM Positioning System,” in Proceedings of the GLOBECOM, Miami, PL, USA, December 2010. [10] S. Peschke, M. Bevermeier, and R. Haeb-Umbach, “A GPS positioning approach exploiting GSM velocity estimates,” in Proceedings of the 6th Workshop on Positioning Navigation and Communication (WPNC 2009), Dresden, Germany, 2009. [11] S. Panzieri, F. Pascucci, and G. Ulivi, “An outdoor navigation system using GPS and inertial platform,” IEEE/ASME Transactions on Mechatronics, vol. 7, no. 2, pp. 134– 142, 2002. [12] H. Liu, H. Darabi, P. Banerjee, and J. Liu, “Survey of wireless indoor positioning techniques and systems,” IEEE Transactions on Systems, Man, and Cybernetics, Part C: Applications and Reviews, vol. 6, pp. 1067–1080, 2007. 112 References 113 [13] H. Cho, H. Jang, and Y. Baek, “Practical localization system for consumer devices using zigbee networks,” IEEE Transactions on Consumer Electronics, vol. 56, no. 3, pp. 1562–1569, 2010. [14] J. Yoon, J. Kim, W. Lee, and D. Eom, “A TDoA-Based Localization Using Precise Time-Synchronization,” in Proceedings of the 14th International Conference on Advanced Communication Technology (ICACT), PyeongChang, February 2012. [15] X. Li and K. Pahlavan, “Super-resolution toa estimation with diversity for indoor geolocation,” IEEE Transactions on Wireless Communications, vol. 3, no. 1, pp. 224–234, 2004. [16] A. Ali and A.S. Omar, “Time of Arrival Estimation for WLAN Indoor Positioning Systems using Matrix Pencil Super Resolution Algorithm ,” in Proceedings of the 2th Workshop on Positioning, Navigation and Communication (WPNC), Hannover, Germany, March 2005. [17] M. Li and Y. Lu, “Angle-Of-Arrival Estimation for Localization and Communication in Wireless Networks,” in Proceedings of the 16th European Signal Processing Conference (EUSIPCO), Lausanne, Switzerland, August 2008. [18] P. Bahl and V.N. Padmanabhan, “RADAR: An In-Building RF-Based User Location and Tracking System,” in INFOCOM 2000. Nineteenth Annual Joint Conference of the IEEE Computer and Communications Societies. IEEE, 2000, vol. 2, pp. 775–784. [19] H. Nurminen, J. Talvitie, S. Ali-Löytty, P. Müller, E. Lohan, R. Piché, and M. Renfors, “Statistical path loss parameter estimation and positioning using rss measurements,” Journal of Global Positioning Systems, vol. 12, no. 1, pp. 13–27, 2013. [20] A. W. S. Au, C. Feng, S. Valaee, S. Reyes, S. Sorour, D. Gold, K. Gordon, and M. Eizenman, “Indoor tracking and navigation using received signal strength and compressive sensing on a mobile device,” IEEE Transactions on Mobile Computing, vol. 12, no. 10, pp. 2050–2062, 2013. [21] P. Bahl and V.N. Padmanabhan, “Enhancements to the RADAR User Location and Tracking System,” in Technical Report MSR-TR-2000-12. Microsoft Research, February 2000. [22] K. Kaemarungsi and P. Krishnamurth, “Modeling of indoor positioning systems based on location fingerprinting,” in Proceedings of the INFOCOM, Hong Kong, March 2004. [23] T. Roos, P. Myllymaki, H. Tirri, P. Misikangas, and J. Sievanen, “A probabilistic approach to wlan user location estimation,” International Journal of Wireless Information Networks, vol. 9, no. 3, pp. 155–164, 2002. [24] A. Agiwal, P. Khandpur, and H. Saran, “Locator: location estimation system for wireless LANs,” in Proceedings of the 2nd ACM International workshop on Wireless mobile applications and services on WLAN hostpots. ACM, 2004, pp. 102–109.