Decimeter Level Indoor Localization Using WiFi Channel State Information

—Indoor localization using WiFi signal parameters is challenging, with encouraging decimeter localization results available with enough line-of-sight coverage and hardware infrastructure. This paper proposes a new 2-dimensional multiple packets based matrix pencil (2D M-MP) method to estimate the Angle of Arrival (AoA) and Time of Flight (ToF) based on WiFi channel state information (CSI). Compared with the conventional parameter estimation algorithms, this method has two advantages. First, 2D M-MP method uses the discrete Fourier transform (DFT) to convert the complex computation into real computation to reduce the computational complexity signiﬁcantly without losing accuracy. Second, it accumulates multiple CSI packets to improve the parameter estimation accuracy effectively, especially at low values of signal-to-noise-ratio (SNR) environment. To verify the practicability of our proposed 2D M-MP method, we set up a localization system in an actual scenario using commodity WiFi cards which demonstrates that the performance of 2D M-MP method is better than conventional parameter estimation algorithms and can achieve a localization accuracy of 42 cm in indoor hall deployment.


I. INTRODUCTION
L OCATION information plays an important role in many applications, such as mobile navigation, smart home, and augmented reality, and accurate localization has always been a hot research field. In recent years, indoor localization based on WiFi signals has been greatly developed due to channel state information (CSI), which describes the amplitude and phase of WiFi signals and can be used to estimate signal parameters, such as the Angle of Arrival (AoA) [1], [5], Time of Flight (ToF) [2]- [4], Angle of Departure (AoD) [16], [17], and the Doppler frequency shift [8], [19], [20] to achieve high precision localization and tracking [7], [8]. At present, many researches focus on the estimation of WiFi signal parameters. SpotFi [1] achieves joint estimation of AoA and ToF using an improved 2D Multiple Signal Classification (MUSIC) algorithm and attains a median localization accuracy of 40 cm. ToneTrack [2] uses 1D MUSIC algorithm to estimate ToF and a novel spectrum identification scheme is proposed to retrieve useful information from a ToF profile to achieve a median localization accuracy of 90 cm. Through constructing an array manifold matrix and designing a spatial smoothing method, 2D AoA and ToF are estimated in [5]. Multi-dimensional parameters, such as AoA, ToF, and Doppler frequency shift are jointly estimated in [9], which achieves passive localization and behavior perception. However, conventional parameter estimation algorithms, such as MUSIC algorithm [21], Space-alternating Generalized Expectationmaximization (SAGE) algorithm [22], Orthogonal Matching Pursuit (OMP) algorithm [23], have some limitations in some ways. MUSIC algorithm relies on a computationally expensive grid search for parameter estimation. The parameter estimation accuracy of the SAGE algorithm depends on the correct number of multipath inputs, and errors in estimating the number of multipath at low signal-to-noise-ratio (SNR) will lead to a local optimal problem for parameter estimation. The accuracy of the parameter estimation of OMP algorithm is related to the SNR of the received signals, if the values of SNR are low, the matching degree between the received signals and the dictionary matrix will be limited, and the performance of OMP algorithm will decrease.
Indoor localization based on CSI parameter estimation is faced with the following challenges. First, the indoor environment is complex, the reflection of the wall and the walking of pedestrians will lead to the low values of SNR, some conventional algorithms can no longer effectively resolve the signal parameters under the low values of SNR environment. Second, indoor localization and tracking have high requirements for real-time performance. However, most of the stateof-art systems have high complexity on parameter estimation and are far from meeting its requirements. This paper proposes to use 2-dimensional multiple packets based matrix pencil (2D M-MP) method to solve the above two challenges to achieve AoA and ToF estimation jointly. 2D M-MP method process the multiple CSI packets to construct the 2D Hankel complex matrix and the discrete Fourier transform (DFT) is used to convert the complex computation into real computation to reduce the computational complexity significantly without losing accuracy. Meanwhile, 2D M-MP method accumulates multiple CSI packets to estimate signal parameters, which can get better parameter estimation accuracy compared with only using a single CSI packet. The figure posted in the abstract depicts the system diagram, the CSI obtained from APs in the experimental deployment is used in 2D M-MP method to obtain AoAs and ToFs, the clustering method [1] is applied to get AoAs of direct path, and weighted least square method (WLS) [31] is used to estimate the target locations at last. Our experimental results show the performance and robustness of this method in three different scenarios.
To summarize, our main contributions of this paper are as follows: • We propose a novel 2D M-MP method to jointly estimate AoA and ToF based on WiFi channel state information (CSI), which can accumulate multiple CSI packets to improve the parameter estimation accuracy effectively, especially at low values of SNR.
• The low complexity of 2D M-MP method provides a prospect for the future real-time localization and tracking platform.
• Our 2D M-MP method is implemented and evaluated with the commodity WiFi platform. Extensive experiments demonstrate the high accuracy and significant performance improvement over the conventional parameter estimation algorithms.
The rest of this paper is structured as follows. The related work is provided in Section II, and the signal model is presented in Section III. The detailed process of 2D M-MP method and the localization scheme are presented in Section IV and V. The experimental setup is discussed in Section VI. We conclude this paper in section VII.

II. RELATED WORK
Various indoor localization systems based on signal parameters have been proposed in recent years, and we divide them into two categories: Received Signal Strength Indication (RSSI) -based, CSI-based techniques.
RSSI-based: RSSI is widely used for indoor localization because it can be easily obtained in the WiFi commercial network card and Bluetooth. A common method is to establish a signal propagation model for ranging and use triangulation [24], [25] for locating the target. However, RSSI values have fluctuations because there are multipath propagations in an indoor environment, and these systems can achieve a median localization accuracy of 2-5 m. The other method is RSSI fingerprint-based [26], which is composed of two phases: offline phase and online phase. In the offline phase, the RSSI values of the reference positions are collected to establish a mapping relationship between reference positions and RSSI values as a fingerprint database. In the online phase, the RSSI obtained from the measurement is used to match the RSSI values of the fingerprint database for locating the target. However, this method will consume a lot of time establishing the fingerprint database. Moreover, these systems are difficult to deploy and are vulnerable to environmental changes.
ToF-based localization method obtains the distance between the transmitter and the receiver by multiplying ToF with the light speed and uses triangulation for locating the target [2], [3], [18], [29]. To our knowledge, most research works [27], [28] use software-defined radio (SDR) to obtain ToF, which is an open-source and fine-tuned platform, but it is expensive to purchase and cannot be expanded to commercialization. For common WiFi commercial network card, the real ToF is difficult to obtain due to the insufficiency of synchronization between transmitter and receiver. To solve this problem, Chronos [29] used the zero subcarriers (at different frequency bands) to measure ToF, which dodged all the deteriorations in CSI phase. Splicer [3] adopted the probability distribution of phase errors to eliminate interferences. However, the two systems need to use frequency hopping technology and modify the communication protocol. SAIL [18] was another system that can localize a target with a single access point using round trip ToF. However, it relied on external IMU sensors on the target. The ToF-based localization method can achieve high accuracy by using real ToF since the transmitter and receiver modules are clock-synchronous. RFind [14] enabled emulating a large virtual bandwidth to estimate ToF and achieved centimeter-scale localization accuracy.
With the development of MIMO technology, antenna arraybased techniques that use multiple antennas at the access point (AP) have attracted wide attention. Specifically, AoA can be estimated by using a corresponding phase shift introduced across the multiple antennas in the array, and triangulation can be applied to localize [1], [5], [15]- [17]. The localization accuracy is related to the accuracy of AoAs estimation and the number of AP. For AoAs estimation of multipath signals, the number of antennas at the AP should be greater than the number of signal paths, otherwise, the accuracy of AoAs estimation will seriously deteriorate. Arraytrack [15] configured multiple antennas at the APs to receive the signals and applied MUSIC algorithm to estimate the AoA of the signals, but it needed to modified the hardware equipment. SpotFi [1] can accurately estimate AoA of the direct signal by using CSI, and an improved 2D MUSIC algorithm was applied to overcome the limits of antenna number and coherent signals to some extent. Generally speaking, accurate localization requires information sharing and coordination of several APs.
Multi-dimensional signal parameters, such as AoA, ToF, and AoD can be applied to realize joint localization. Theoretically, it can achieve better localization accuracy than using the single signal parameter. M3 [16] presented a super-resolution algorithm SAGE to jointly estimate channel parameters including AoA, AoD, and relative ToF, and centimeter-level localization accuracy was achieved. mD-Track [9] incorporated information from as many signal dimensions as possible to advance the accuracy of passive wireless sensing in a multipath environment and achieved passive multi-target localization and motion tracking. MonoLoco [17] had designed a new model for the multipath channel that enables super-resolution methods to estimate AoA, AoD, and the relative ToF between paths, and it leveraged multipath reflections to estimate the location of a target and a reflector with respect to the receiver.

III. SIGNAL MODEL
Considering an AP with M identical antennas, the antennas are arranged in a straight line with uniform spacing in space, and the uniform spacing is a half wavelength. At the same time, the transmitter uses a single antenna to transmit Orthogonal Frequency Division Multiplexing (OFDM) signals, within the bandwidth, the N subcarrier frequencies {f 0 , . . . , f N −1 } are equally spaced by ∆f .
For I CSI packets, the i-th (1 ≤ i ≤ I) CSI measurement can be expressed as is the CSI measurement at the m-th antenna of AP and the n-th subcarrier frequency of i-th CSI packet, for L multipath propagation paths, h i m,n can be denoted as Here, d and f are the equal spacing between consecutive antennas and signal frequency, respectively, α i l and c are the signal (complex) amplitude of the l-th path and the light speed, respectively, θ i l and τ i l are the angle of arrival and propagation delay of the l-th path, respectively, ω i m,n represents the noise. For simplicity, we have . The 2D M-MP method proposed in this paper will focus on the estimation of u l and v l , which mean the CSI phase difference between two adjacent antennas in space domain and the CSI phase difference between two adjacent subcarriers in time domain, respectively. we can write where θ l and τ l present the AoA and ToF of the lth (1 ≤ l ≤ L) path in multiple CSI packet data scenery.

A. Multiple CSI Packet Data Matrix Construction
As part of our overall method, we first need to construct multiple CSI packet data matrix. By using the principle of original 2D MP method [12], a data matrix V i can be constructed as where dimensional Hankel matrix, K and P are the pencil parameter, which can be adjusted to filter the noise effectively. According to [11], K and P should be satisfied as In general, K and P can be chosen between M/3 and M/2, N/3 and N/2 to increase the estimation accuracy. The original MP method only considers using a single CSI packet. For I CSI packets, we can get data matrices V 1 , . . . , V i , . . . , V I and the multipath enhanced data matrix V E can be constructed as Theoretically, the i-th CSI packet data matrix V i can also be written as where where in (8) and (14), denotes the Khatri-Rao product, and the enhanced data matrix V E can be also written as

B. DFT Transformation
The original MP method directly calculates the complex data matrix to obtain the required parameters. In M-MP method, we involve only real-valued computation from start to finish after introducing the DFT matrix [34]. If W K is a DFT matrix of K × K size, and the k-th row of the conjugate centro-symmetrized version W K H will be The row vector w k H represents a DFT beam steered at a spatial frequency µ=k ×2π/K. The W P is also a DFT matrix and the p-th row of W P H can be expressed in the same way as w k H . Now we introduce two K × L and P × L beamspace manifold matrices whose l-th columns can be written as where Comparing b k u i l with b k+1 u i l , the numerator of b k+1 u i l is observed to be the negative of that of b k u i l . Thus, trigonometric manipulations enable the relationship between b k u i l and b k+1 u i l can be described as Compiling all K equations in vector form yields an invariance relationship for the beamspace manifold b K u i l as follows where With L multiple propagation paths, the beamspace manifold relation in (21) can further translate into the matrix relation as where The relationship between b p v i l and b p+1 v i l is also similar to b k u i l and b k+1 u i l , and we can also get where the differences between Γ 3 , Γ 4 and Γ 1 , Γ 2 are that K is replaced by P and the u of Ω i u is replaced by v. Next, we define Γ u1 = Γ 1 ⊗ I P and Γ u2 = Γ 2 ⊗ I P , where ⊗ denotes the Kronecker product, I P is P × P identity matrix, and we can get According to (24) and (27) in the above equation and we get Similarly, we can also get where Γ v1 = I K ⊗ Γ 3 and Γ v2 = I K ⊗ Γ 4 , and I K is K × K identity matrix. Now we use the DFT matrix W K H ⊗ W P H to multiply data matrix V E from left to obtain a matrix V as where the superscript H denotes the conjugate transpose, and U s wr contains the left singular vectors span the signal subspace corresponding to L largest singular values. Now we get the matrix U s wr = TU s wr , where T is defined as T = [r(1), r(1 + P ), . . . , r(1 + (K − 1)P ), r(2), r(2 + P ), . . . , r(2 + (K − 1)P ), . . . , r(P ), r(P + P ), . . . , r(P + (K − 1)P )] T where r (p + (k − 1) P ) is a KP size row vector and the value of the l-th column is 1 and others are zeros. We know U s wr and B i K B i P share the same column space, according to (28) and (29), we can get \\ AoAs and ToFs parameter calculation; where Ψ u and Ψ v are both L×L matrices and the eigenvalue of them can be defined as λ l and ζ l (1 ≤ l ≤ L), corresponding to tan u l 2 and tan v l 2 , respectively, and the θ l and τ l can be calculated as The pairing algorithm [10], [12] can be used to pair the obtained {θ 1 , . . . , θ L } and {τ 1 , . . . , τ L } to find the appropriate group for them, then complete the two-dimensional AoA and ToF estimation. It can be realized the multipath channel poles of the time and space dimensions have the same generalized eigenvectors Ω when they are decomposed using the generalized eigenvalue decomposition. Based on this property, according to Equation (33), we can get where Ψ u and Ψ v are both L×L matrices and the eigenvalue of them can be defined as λ l and ζ l (1 ≤ l ≤ L). In this way, the poles of two dimensions are in a grouped form, and the required parameters can be calculated as where the obtained parameters are in pair form {θ 1 , τ 1 ; . . . ; τ L , θ L } without the need to solve an eigenvalue problem for each dimension. We summarize 2D M-MP method in Algorithm 1. The input multiple CSI packet matrix H 1 , . . . , H i , . . . , H I is used to construct Hankel block matrix V E , and then the DFT matrix W K H ⊗ W P H is constructed to get the matrix V by left multiplication of V E . Next, an SVD is used on [Re {V} , Im {V}] to obtain U s wr and [U s wr ] to achieve the calculation of the two-dimensional parameters θ l and τ l .

V. CLUSTERING AND LOCALIZING THE TARGETS
The AoA and ToF parameters can be obtained by using 2D M-MP method, considering ToF is not available due to the phase errors, we only focus on AoA for localization. The environment of the indoor scenario is complex, and the walls, ceilings and desks will cause signal reflection, generally, the APs will receive 2-4 multipath signals with relatively strong power, which makes it difficult for estimating the AoA of the direct path. Here, the clustering method [1] is used to identify AoA of the line of sight (LoS) path. As Fig. 1 shows, the red dotted line corresponds to the LoS path, whose power is more stable and the fluctuation of direct signal scatters is smaller than other reflected signal scatters. Next, as Fig. 2 shows, we use AoAs from two APs for localization by meeting at a point of two straight lines or use AoAs from three APs for localization by using the WLS method [31] to solute the fuzzy region.

VI. EXPERIMENTAL SETUP
As Fig. 3 shows, we use ASUS RT-AC86U as the AP with four antennas, which is equipped with the Nexmon CSI Extractor Tool [13] to obtain CSI information and supports 80 MHz bandwidth at 5.805 GHz, and the distance between two adjacent antennas is equal to 2.58 cm (half a wavelength). Honor X1 is used as the transmitter with one antenna whose packet transmission rate is set to 50 Hz. We operate the transmitter and a PC equipped with a WiFi network card in the same local area network and control the PC to request the data from the transmitter through the ping command. Meanwhile, the APs operate in monitor mode, and the CSI data of 256 subcarriers of 4 antennas can be obtained by using Nexmon CSI extractor tool. In TABLE I, our experimental parameters were listed in detail. We use different algorithms to estimate the parameters based on the CSI data obtained and then locate the target locations. All the data processing uses MATLAB R2016a software in our experiments.  We conduct the experiments in three different scenarios, as Fig. 4 shows, the first is an open roof (18m × 20m) with almost no reflection, the second is an indoor hall (8m × 15m) with less reflection, and the third is a small laboratory (7.5m × 12.5m) with strong multipath reflection. The blue APs and the red points are the receivers and the transmitter target positions, respectively, and three experiments test 26, 56, and 40 target locations respectively. We collect two sets of CSI data in each group of target position, which contains about 300 CSI packets in total. The height of the APs and transmitters is identical and ground truth location and orientation were confirmed by using laser rangefinders, building protractors, and geometric relationships between transmitter and receiver. In our experiments, the proposed 2D M-MP method is compared to some previous works, which are based on the 2D MUSIC algorithm [1], [9], [15], original 2D MP algorithm [10], [12] and 2D OMP algorithm [30], [32]- [33] for evaluating the parameter estimation accuracy and localization accuracy. Due to the lack of synchronization between transmitter and receiver, ToF is not available, so we only evaluate AoA accuracy and use triangulation for localization. The comparison is conducted from three aspects. First, we analyze the AoA accuracy and localization accuracy in each experimental scenario. Second, the overall AoA accuracy and localization accuracy are evaluated between different algorithms under different scenarios. Third, we compare the computational complexity of different algorithms through the complex multiplication numbers and the program execution time. Everything in each experiment is the same except for the different algorithms.

B. Experiment 2
As Fig. 4 shows, we repeat the same experiment as before in the indoor hall. Fig. 6 shows the AoA estimation accuracy and localization accuracy in the indoor hall. Considering the estimation error under the ratio of 0.667, we can see the AoA estimation accuracy of 2D M-MP method, 2D MP method, 2D MUSIC algorithm, and 2D OMP algorithm can reach about 4.56 • , 6.84 • , 7.25 • , and 8.93 • , respectively, and the localization accuracy of them can reach about 0.42 m, 0.58 m, 0.76 m, and 1.09 m, respectively. Compared with the previous scenario, the indoor hall is less spacious and the reflected signals are more, and the AoA estimation accuracy and localization accuracy are reduced to some extent. Meanwhile, we can see that 2D M-MP still has the best performance.

C. Experiment 3
As Fig. 4 shows, we deploy two APs and test 40 target positions between the desks in the small laboratory. Fig. 7 shows our AoA estimation accuracy and localization accuracy  Considering the environment of the small laboratory is most complex, only two APs are used for localization which has led to greater fluctuation and variance of the localization error than before. So the CDF in Fig. 7(right) is not smooth enough compared to the two localization error CDF before, the AoA estimation accuracy and localization accuracy are further reduced than before. Meanwhile, the 2D M-MP method still shows the best performance of them.

D. Overall Evaluation
Parameter estimation error and localization error are related to many factors, such as the multipath environment, the number of AP, the wall materials, and so on. To intuitively compare the parameter estimation accuracy and localization accuracy of the previous three experiments, we show their estimation error of each algorithm under the ratio of 0.667 from the previous CDF in Fig. 8. We can see that no matter which algorithm is adopted, the AoA estimation accuracy and localization accuracy of the open roof are the best, followed by the indoor hall and the small laboratory, meanwhile, 2D M-MP method always show the best performance, and we analyze

E. Computational Complexity Evaluation
In this part, the complexity of 2D M-MP method is compared with 2D MP method, 2D MUSIC algorithm and 2D OMP algorithm. For the conventional 2D MP method, the complexity mainly focuses on an SVD of a data matrix, and the data matrix is complex of size KP × (M − K + 1)(N − P + 1) and its SVD requires 11(KP ) 3 + 4(KP ) 2 (M − K + 1) (N − P + 1) complex multiplications. For the 2D MUSIC algorithm, the complexity mainly focuses on eigenvalue decomposition and seeking the spectrum peak. For 2D OMP algorithm, its complexity mainly lies in the construction and search of the dictionary matrix. For 2D M-MP method in this paper, the complexity analysis can be divided into two sub-steps. In the first sub-step, a matrix W K H ⊗ W P H of dimension KP × KP is multiplied with the data matrix V E from left. In the second sub-step, an SVD of a real matrix of size KP × I(M − K + 1)(N − P + 1) is computed. TABLE II shows the complexity comparison of the four algorithms in detail, we list the main steps with the high complexity of different algorithms and the computational burden of other steps is negligible and therefore not discussed here. Based on our experiments, we quantify the complexity for comparing intuitively and the meanings and values of variables are shown in TABLE I. We can see that the complexity of 2D M-MP method is the lowest among them. Theoretically, the more CSI packets accumulated, the higher accuracy of parameter estimation will obtain, while it will increase the computational complexity. But the introduction of DFT can reduce the complexity to some extent. In reality, the CSI packet number will be adjusted according to the actual demand. According to the complexity analysis before, we know that the complexity of 2D M-MP method is the lowest. To intuitively observe the real-time performance of this method, we compare the program execution time of the four algorithms. Note that except for using different algorithms, others are the same, including using the same PC, processing the same CSI data, etc. To make a more effective comparison, we run the program of each algorithm several times and record its average running time. As shown in Fig. 9, we can see the program execution time of 2D M-MP method, 2D MP method, 2D MUSIC algorithm and 2D OMP algorithm can reach about 0.074s, 0.362s, 5.367s and 4.281s, respectively, and 2D M-MP method is the lowest. To a certain extent, this method breaks through the problem of the poor real-time performance of the conventional algorithm, and has important significance for the modern real-time commercial platform.

VII. CONCLUSION
In this paper, we propose 2D M-MP method to jointly estimate AoA and ToF parameters based on WiFi signals, which accumulates multiple CSI packets to improve the parameter estimation accuracy effectively. At the same time, based on the low complexity of the method, it is of great significance to the application of real-time localization and tracking platform in the future. Our parameter estimation and localization experiments in three typical environments show that the performance and robustness of 2D M-MP method are better than state-of-art systems.