I Introduction
Multipath channel estimation using WiFi frequency domain channel state information (CSI) has found numerous applications in indoor localization and ranging in recent years [1, 2, 3, 4]. Such applications rely on the channel delay profile to estimate the distance between two devices or detect the motion of an object. The accuracy of such schemes is limited by the resolution with which they can estimate the delay components of the channel. This resolution is limited by the inverse of the overall bandwidth of the WiFi device. For example, using CSI measurements over a bandwidth of MHz, one can localize channel delay components separated by ns, whereas when the bandwidth increases to GHz, one can improve the resolution to ns [3]. Multiplying the delay resolution with the speed of the electromagnetic wave, results in a ranging error of up to m in the former versus up to only cm in the latter case. Therefore, a large channel bandwidth is crucial for any application demanding precise localization. ^{†}^{†}This work was carried out within the “Secure FogConnection Layer for IoT Applications (SecureFog)” project, funded by the Federal Ministry of Education and Research (BMBF), Germany, under grant number FKZ: 16KIS0776.
Since single WiFi bands usually span over MHz, using each of them separately results in poor localization. However, the total bandwidth allocated to a WiFi device, which consists of all the used WiFi bands, may scale up to several hundred MHz or even multiGHz, making WiFi CSI measurements over the wide bandwidth a useful resource for decimeterlevel localization tasks. The process of exploiting such measurements over multiple bands is referred to as multiband splicing [3]. Despite the potential of multiband splicing, channel impulse response (CIR) estimation using raw WiFi CSI measurements is challenging, mainly due to imperfections of the underlying hardware. Perband CSI measurements are performed independently based on sequential frequency hopping, which introduces random and independent phase offsets to the measurements at each band. Furthermore, time synchronization errors introduce a linear phase rotation to the CSI, which is independent from one band to the other [2].
The existence of such distortions makes a straightforward CIR estimation via multiband splicing impossible. A few methods have been recently proposed to solve this problem, by removing the phase distortions using signal processing tools. In [2], the Chronos system was proposed, which uses the CSI measured only on a single, specially chosen subcarrier per band in order to recover the time of flight (ToF), which is the delay corresponding to the first path, assuming a sparse multipath channel. Using a single subcarrier per band may lead to a poor estimation of the ToF, since for channels with a large delay spread, it results in an insufficient sampling rate and an aliased estimation of the CIR. In addition, the CSI measurements over the remaining subcarriers of a band can provide valuable information, especially in low signaltonoise ratio (SNR) scenarios, but are neglected. The problem of spectrum splicing with the application of accurate ToF estimation was also treated in [3] and [4]
, where heuristic methods were proposed with overall performances inferior to the method of
[2], which we take as the state of the art for comparison to our work.In this paper we develop a method for multiband splicing with phasedistorted and noisy WiFi CSI measurements. This method operates only on the magnitude of the channel frequency samples over multiple WiFi bands to estimate the CIR, a technique which is known as phase retrieval (PR) in the literature [5]. As is well known, applying PR to frequency samples results in a fundamental ambiguity in the estimated delay domain signal, in terms of constant shifts and conjugatereflections with respect to the origin [6]. We resolve this ambiguity using a handshaking process between the transmitter and the receiver and exploiting the channel reciprocity [7]. We show that our method can recover the channel delay components and their coefficients with a high resolution for various SNR values, which results in significant ranging accuracy.
Notation: for a positive integer we denote by the set
. For a vector
, denotes elementwise absolute values. For a function, we define the Fourier transform as
. For a set of integers and a vector , denotes the elements of in .Ii System Setup
Assume that the propagation environment between a WiFi transmitter and a receiver consists of scatterers, each imposing a path gain and delay , where denotes the maximum delay spread. Using OFDM signaling, the WiFi device transmits pilots to the receiver over frequency bands, each consisting of equispaced subcarriers.^{1}^{1}1Note that the bands need not be adjacent, as is sometimes the case where the WiFi device uses multiple bands around GHz and GHz [2]. The receiver collects the CSI of each band separately, by downconversion through removing the carrier frequency. The baseband representation of the CIR corresponding to band is thereby given as where denotes the carrier frequency of the th band and denotes Dirac’s delta. Since only a few propagation paths contribute to the channel, i.e. is a small number, the CIR admits a sparse representation over the delay domain, which will be exploited in our proposed method. The channel frequency response (CFR) samples over the set of subcarriers ^{2}^{2}2
is an odd integer.
, with being the carrier spacing, are denoted by the vector where(1) 
where .
Unfortunately, these (raw) CSI measurements are subject to several phase and magnitude distortions as well as additive noise due to hardware imperfections. We denote the CSI of band by the vector where each element can be expressed as [2, 3, 8]
(2) 
where is an affine phase distortion term with denoting the bandspecific time offset due to the packet detection delay (PDD) and denoting the constant phase offset. Also , by assuming that the support^{3}^{3}3The support of a function denotes the set . of the effective CIR (including the shift caused by the PDD) is contained within the OFDM cyclic prefix, is the additive white Gaussian noise (AWGN) vector, and is the power control gain. For simplicity, we set in what follows, since its effect can be subsumed to the noise level. We refer the interested reader to references [2, 3, 4] for an extensive account on the physical source of these distortions and focus on the mathematical presentation of the problem.
Iia Denoising and ZeroSubcarrier CSI Estimation
Before proceeding with the statement of the problem, we should address an issue with the raw WiFi measurements. The WiFi devices do not transmit on the zerosubcarrier of each band, i.e. on the frequencies , since measurements on these subcarriers overlap with the hardware DC offsets that are hard to remove [2]. This means that in each band , the measurement is not observed. However, the CSI on the zero subcarrier is very useful, since the phase distortion term is free of PDD and reduces to a constant offset, i.e.,
. This property will be exploited in the final step of our proposed method. Here, we use the measured CSI over the remaining subcarriers of the band, to interpolate the missing CSI measurement on the zerosubcarrier.
The authors of [2] used a cubic spline interpolator for this task. Motivated by the fact that the CIR is sparse over the delay domain (), we propose a different interpolation/estimation scheme. Our estimator is based on norm minimization, which exploits channel sparsity and will additionally denoise the raw CSI measurements. Recalling (2), one can write where and . Since the signal component of is a linear combination of a small number () of the vectors , it is natural to use a sparse recovery method for both estimating and denoising the observed values of . Here we use the basis pursuit denoising (BPDN) method for this purpose. Define the oversampled DFT matrix where and . This matrix serves as a dictionay over a discretized delay grid with delays over . The BPDN program, thereby, can be formulated as follows [9]
(3) 
where denotes the norm, is an appropriate regularization scalar and . In practice can be any subcarrier index set for which the measurements are available (i.e., pilot symbols are effectively transmitted). Solving (3) gives us an estimate of the sparse coefficients vector , using which we can not only estimate the zero subcarrier CSI as , but also denoise the observed vector . We denote the overall denoised CSI vector over band as .
The above program is convex and is solved for each band separately, using any of the numerous convex solvers. For a fast implementation we use the method proposed in [10], known as SpaRSA. Now, by having (denoised) measurements over all subcarriers in all bands, we are ready to state the problem of CIR estimation.
IiB Problem Statement
Let denote the vector containing the CFR samples over all bands. Using (1) one can easily see that contains the frequency samples of the CIR over the ordered set of subcarrier frequencies . In the same way we define and as the denoised measurements vector and the noise vector, respectively. It follows that the CSI measurement vector can be written as where, recalling (2), is an unknown diagonal matrix, with unit modulus phasedistortion entries. The main problem addressed in this paper is formally stated as follows.
Problem 1
Given the distorted CSI measurements vector , estimate the sparse CIR .
Since the samples are only distorted in phase, and not in magnitude, our idea in this paper is to estimate the CIR by applying a PR algorithm to the magnitude of the CSI, i.e. . An important feature of formulating and solving the problem as such is that, unlike previous methods in the literature, we do not perform any timeconsuming preprocessing to remove the effect of phase distortions and we are able, in principle, to estimate the CIR as soon as a single snapshot of CSI measurement across all bands is collected. This is crucial in delaysensitive applications.
Iii CIR Estimation via Phase Retrieval
Using the formulation developed in the previous section, define
(4) 
where denotes the vector of magnitudes and . It is easy to see that represents the samples of the Fourier transform of the CIR autocorrelation function
(5) 
for , where denotes convolution. In order to recover from we first recover and then estimate the set of delays and their corresponding coefficients (up to a phase constant) from the estimated autocorrelation.
Iiia Sparse Recovery of the CIR Autocorrelation
From (5) it follows that the autocorrelation is conjugate symmetric with respect to . We can write where , and for some . Furthermore, is a sparse function, i.e. assuming to be small, it consists of a few number of delta functions over . Using (4) we can write
(6)  
where , , and denotes the all ones vector of size . In order to recover , we define a dictionary of the vectors on a dense grid of points over and apply sparse recovery methods. Let with be a dense grid over . Define the dictionary where and . We can approximate the signal part of with a linear combination of the columns of and therefore we have where is a sparse vector, where and denote the real and imaginary parts of the coefficients approximating the frequency samples vector of the autocorrelation, respectively. In order to estimate the sparse vector given , we solve the following BPDN problem,
(7) 
where is an appropriate regularization scalar which controls the balance between sparsity of the solution and the approximation accuracy. This parameter depends on the noise level and can be chosen empirically. As explained in Section IIA, the BPDN problem is a convex program and can be efficiently solved using any convex solver. Similar to problem (3), we solve (7) using the SpaRSA algorithm.
Once we solve (7), we have an estimate of a discrete approximation of . Using the notation , the coefficients of this discrete approximation over are given by for and , whereas the coefficients over are simply given by .
This discretized approximation of the autocorrelation usually contains spurious nonzero elements, since the true support entries do not necessarily lie on the grid and also due to noise effects. This means that the support set of , i.e. generally has more than elements. However, by knowing the number of paths contributing to the CIR,^{4}^{4}4The number of paths corresponds to the slowvarying geometry of the environment and can be estimated offline during a preprocessing step. we can refine the support of as follows. Let denote the set of points in whose corresponding coefficients in are large enough, where is a small, predefined threshold. The points in are clusters around the true support points . Knowing the value
, we apply the kmeans clustering method
[11] and obtain the center points of these clusters, denoting them by . Eventually, the estimated autocorrelation support is given by . Now, the coefficients corresponding to points in the estimated support set can be obtained by solving a leastsquares problem(8) 
where . In what follows we refer to the estimated coefficient corresponding to , as . Note that the coefficient corresponding to will be due to conjugate symmetry. This completes the estimation of the autocorrelation function .
IiiB CIR Recovery
In order to recover the CIR from its autocorrelation estimate , we rely on a method developed in [12] which performs support and magnitude recovery in succession. First, note that the support of the autocorrelation is the difference set of the support of the CIR , i.e. (Minkowski difference). Also is a symmetric set containing the point . Now we have an estimate of and want to recover . A crucial point here is that, recovering from involves a fundamental ambiguity even when : for any constant the generic sets and have the same difference set as the . Thus, we can hope to recover only up to constant shifts and reflections with respect to 0. These uncertainties will be resolved in a later step of our algorithm and in this section, we only focus on estimating up to the mentioned ambiguities. We further assume that there exists no collision in the difference set, i.e. there are no two distinct pair of points and for which , which is almost surely guaranteed given that the support elements are driven uniformly at random over the delay domain. The support recovery algorithm is summarized in Algorithm 1 and the interested reader is referred to [12] for further details.
Once we obtain the CIR support estimate , we proceed with estimating the CIR coefficients . Estimation of the coefficients from the autocorrelation would be up to a common constant phase, since it is easy to show that the CIR coefficients and , result in the same autocorrelation. For any two elements we can find their corresponding differences in and furthermore, the corresponding estimated coefficients from the discussion of Section IIIA. From (5), we know that and since is an estimate of we expect . Now construct the matrix such that
(9) 
Note that for . Define the sum of the elements of as . When , using these equations we can obtain as

(10) 
which gives us an estimate of the amplitude of the CIR coefficients.^{5}^{5}5The solution for can be trivially computed and is not discussed here due to a lack of space. Let denote the phase of , i.e., . We estimate the CIR coefficients by treating as realvalued which is of no consequences due to the fundamental phase ambiguity. This vector of coefficients, denoted by is estimated as follows. For any , we simply find such that . The corresponding estimated autocorrelation coefficient for is given by . By dividing to we obtain .
IiiC Resolving Ambiguities
At this point, we have estimated the CIR up to the ambiguity arising from the shift and reflection of the support elements. To resolve this ambiguity, we rely on an observation made in [2]. As discussed in Section IIA, the zero subcarrier in band is free from the PDP phase error and contains only the constant phase error term . This constant phase stems from the PLL phase offset and has the same absolute value but opposite signs on the transmitter and the receiver (see [2] Eqs. (11) and (12)), i.e.,
(11)  
where and denote CSI on the zerosubcarrier of band at the transmitter and the receiver and and denote their corresponding noise terms, respectively. The CFR sample is the same on both ends due to channel reciprocity [7]. During frequency hopping the transmitter and the receiver send packets to each other and they can exchange their CSI measurements. We use this process to exchange interpolated zero subcarrier CSI and by multiplying them at either the receiver or transmitter WiFi device, we get
(12) 
where denotes the noise and cross terms. The vector gives us the extra information with which we resolve the ambiguity about the shift and reflection in the support set . This is done by constructing a simple hypothesis test. Let the elements of be ordered as . In addition, we know that in an indoor environment the first delay component is within a certain limits. Let denote this limit. Our first hypothesis is that the estimated CIR is only a shifted version of the true CIR. Let the tentative CIR representing this hypothesis be denoted by where by construction we made the first estimated delay component to appear at zero, i.e. (see Algorithm 1, line 7). Also define the conjugatereflected and shifted version of to be which represents the CIR of the second hypothesis . Define the Fourier transform of the shifted CIRs corresponding to and by a value as and , respectively. The squared samples of these functions represent the squared CFR of the tentative solutions (up to a phase constant) and we will realize which hypothesis is true and which shift value gives the true CIR by comparing these samples with the ones in . Let and denote the squared frequency samples of the tentative solutions for a shift . For the first hypothesis , i.e. the CIR being only shifted, form the cost function . In a similar way, for the second hypothesis , i.e. the CIR being conjugatereflected define the cost function as . We select over if
and over otherwise. In addition, the optimal value of the shift parameter is given by where is the winning hypothesis. After obtaining , the eventual set of estimated delays is given by . This completes the CIR estimation process and denotes the estimated ToF.
Iv Simulation Results
We compare our proposed method to the Chronos system developed in [2]. The method used in this system first interpolates the CSI on the zero subcarrier in each WiFi band using a cubic spline interpolator. Then, the zerosubcarrier CSI is exchanged between the WiFi device and the client to get the vector in (12). This vector represents the frequency samples of the convolution of the CIR with itself. A sparse recovery method is applied to to get a discretized estimate of the convolution. The first nonzero component in the estimated vector is located at twice the ToF and this gives a way to estimate the ToF and the eventual localization.
To run the simulations we consider a CIR with
paths whose delays are randomly located over the delay domain, such that no elements in the generated delay difference set, corresponding to different delay values, are the same. To each path a complex Gaussian coefficient is assigned. The coefficient assigned to a path with a larger delay has a smaller variance, since paths with larger delay usually have smaller power. We consider
adjacent WiFi bands, each with subcarriers, with a carrier spacing of kHz. We set the parameters in programs (3) and (7) as and , respectively. Appropriate regularization scalars and are computed as a function of the noise power using a training set prior to the experiments. Given an estimate for the first delay component for either of the methods as , the ranging error is computed as (in meters), where denotes the speed of light. We run the experiment for MonteCarlo trials. The ranging error CDF is illustrated in Fig. 1 for two SNR values of and dB. As we can see, our proposed method outperforms Chronos in about of the trials, and has a ranging error less than cm in about of the trials, indicating a cmlevel localization accuracy. Notice that even when SNR=10 dB, the PRbased method performs better than Chronos fed with CSI measurements with SNR=20 dB in about of the times, showing the noise resilience of our proposed method. In a small fraction of the experiments our method has a larger ranging error and performs worse compared to Chronos. These outlying results occur due to the errors in the identification of the difference set of the delays.V Conclusion
We proposed a novel method for indoor localization via channel impulse response estimation using WiFi multiband CSI measurements. Our method was based on the phase retrieval of CSI magnitudes, overcoming the inherent signal phase distortions caused by the WiFi hardware. We empirically showed that this method improves the channel delay profile resolution, which in turns results in a more accurate localization of the clients in the environment in both high and lowSNR regimes.
References
 [1] M. Kotaru, K. Joshi, D. Bharadia, and S. Katti, “Spotfi: Decimeter level localization using WiFi,” in ACM SIGCOMM Computer Communication Review, vol. 45, no. 4. ACM, 2015, pp. 269–282.
 [2] D. Vasisht, S. Kumar, and D. Katabi, “Decimeterlevel localization with a single WiFi access point.” in NSDI, vol. 16, 2016, pp. 165–178.
 [3] Y. Xie, Z. Li, and M. Li, “Precise power delay profiling with commodity WiFi,” IEEE Trans. on Mobile Computing, 2018.
 [4] Y. Zhuo, H. Zhu, H. Xue, and S. Chang, “Perceiving accurate CSI phases with commodity WiFi devices,” in INFOCOM 2017IEEE Conference on Computer Communications, IEEE. IEEE, 2017, pp. 1–9.
 [5] E. J. Candes, T. Strohmer, and V. Voroninski, “Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming,” Communications on Pure and Applied Mathematics, vol. 66, no. 8, pp. 1241–1274, 2013.
 [6] K. Jaganathan, S. Oymak, and B. Hassibi, “Sparse phase retrieval: Uniqueness guarantees and recovery algorithms.” IEEE Trans. Signal Process., vol. 65, no. 9, pp. 2402–2410, 2017.
 [7] M. Guillaud, D. T. Slock, and R. Knopp, “A practical method for wireless channel reciprocity exploitation through relative calibration.” in ISSPA, 2005, pp. 403–406.
 [8] M. Speth, S. A. Fechtel, G. Fock, H. Meyr et al., “Optimum receiver design for wireless broadband systems using OFDM: Part I,” IEEE Trans. Commun., vol. 47, no. 11, pp. 1668–1677, 1999.
 [9] Y. C. Eldar and G. Kutyniok, Compressed sensing: theory and applications. Cambridge University Press, 2012.
 [10] S. J. Wright, R. D. Nowak, and M. A. Figueiredo, “Sparse reconstruction by separable approximation,” IEEE Trans. Signal Process., vol. 57, no. 7, pp. 2479–2493, 2009.
 [11] C. M. Bishop, Pattern recognition and machine learning. Springer, 2006.
 [12] G. Baechler, M. Kreković, J. Ranieri, A. Chebira, Y. M. Lu, and M. Vetterli, “Super resolution phase retrieval for sparse signals,” arXiv preprint arXiv:1808.01961, 2018.
Comments
There are no comments yet.