WiFi-Based Indoor Localization via Multi-Band Splicing and Phase Retrieval

TL;DR

This paper introduces a WiFi-based indoor localization method that combines multi-band splicing and phase retrieval to achieve centimeter-level accuracy, overcoming hardware distortions and bandwidth limitations.

Contribution

It proposes a novel multi-band CSI splicing and phase retrieval approach with sparse CIR recovery, significantly improving indoor localization accuracy over existing methods.

Findings

Achieves cm-level localization accuracy in simulations

Outperforms state-of-the-art methods by an order of magnitude

Effective across various SNR regimes

Abstract

We study the problem of indoor localization using commodity WiFi channel state information (CSI) measurements. The accuracy of methods developed to address this problem is limited by the overall bandwidth used by the WiFi device as well as various types of signal distortions imposed by the underlying hardware. In this paper, we propose a localization method that performs channel impulse response (CIR) estimation by splicing measured CSI over multiple WiFi bands. In order to overcome hardware-induced phase distortions, we propose a phase retrieval (PR) scheme that only uses CSI magnitude values to estimate the CIR. To achieve high localization accuracy, the PR scheme involves a sparse recovery step, which exploits the fact that the CIR is sparse over the delay domain, due to the small number of contributing signal paths in an indoor environment. Simulation results indicate that our approach outperforms the state of the art by an order of magnitude (cm-level localization accuracy) for more than 90% of the trials and for various SNR regimes.

Equations

$${\tilde{\mathbf{h}}}_n^{(m)}=\sum _{k=0}^{K-1}{c}_k{e}^{-j2\pi (f_{m,0}+n{f}_s){\tau }_k}=\sum _{k=0}^{K-1}{c}_k{e}^{-j2\pi f_{m,n}{\tau }_k},$$

$${\mathbf{y}}_n^{(m)}=\alpha e^{-j{\varphi }_{m,n}}{\tilde{\mathbf{h}}}_n^{(m)}+{\mathbf{z}}_n^{(m)},$$

$${\tilde{\mathbf{w}}}^{(m)}=\underset{\mathbf{w}\in {\mathbb{C}}^L}{\mathrm{argmin}}\frac{1}{2}\Vert {\mathbf{y}}_J^{(m)}-{\left(\mathbf{D}\mathbf{w}\right)}_J{\Vert }_2^2+\rho \Vert \mathbf{w}{\Vert }_1,$$

$${\mathbf{u}}_i\stackrel{\Delta }{=}|{\mathbf{y}}_i{|}^2=|{\tilde{\mathbf{h}}}_i{|}^2+{\tilde{\mathbf{z}}}_i,$$

$$R(\xi )\stackrel{\Delta }{=}h(\tau )\star h^{\ast }(-\tau ){|}_{\xi }=\sum _{k=0}^{K-1}\sum _{\ell =0}^{K-1}{c}_k{c}_{\ell }^{\ast }\delta (\xi -({\tau }_k-{\tau }_{\ell })),$$

$$\displaystyle\scalebox{0.94}{${\bf u}$}$$

$$\displaystyle\scalebox{0.94}{$=r_{0}\mathbf{1}+2\sum_{s=1}^{\tilde{K}}\left({\rm Re}\left\{r_{s}\right\}{\bf a}_{\text{re}}(\xi_{s})-{\rm Im}\left\{r_{s}\right\}{\bf a}_{\text{im}}(\xi_{s})\right)+\tilde{{\bf z}},$}$$

$${\mathbf{x}}^{\star }=\underset{\mathbf{x}\in {\mathbb{C}}^{2G+1}}{\mathrm{minimize}}\frac{1}{2}\Vert \mathbf{A}\mathbf{x}-\mathbf{u}{\Vert }_2^2+\lambda \Vert \mathbf{x}{\Vert }_1,$$

$${\mathbf{r}}^{\star }=\underset{\mathbf{r}\in {\mathbb{C}}^{2\tilde{K}+1}}{\mathrm{argmin}}\Vert {\mathbf{A}}_{\hat{\mathcal{D}}}\mathbf{r}-\mathbf{u}{\Vert }_2$$

$${\mathbf{C}}_{k,\ell }\stackrel{\Delta }{=}\{\begin{array}{ll}0& k=\ell ,\\ \log |{\hat{r}}_s|=\log |{\hat{c}}_k|+\log |{\hat{c}}_{\ell }|& k\ne \ell .\end{array}$$

$$\log |{\hat{c}}_k|=\frac{1}{(K-2)}\left(\sum _{\ell =1}^K{\mathbf{C}}_{k,\ell }-\frac{\beta }{2(K-1)}\right),$$

$${\mathbf{y}}_{0,tx}^{(m)}=$$

$${\mathbf{y}}_{0,rx}^{(m)}=$$

$${\mathbf{y}}_m^{\prime }:={\mathbf{y}}_{0,tx}^{(m)}{\mathbf{y}}_{0,rx}^{(m)}=({\tilde{\mathbf{h}}}_0^{(m)}{)}^2+{\mathbf{z}}_m^{\prime },$$

$$\underset{{\tau }_{ϵ}\in [0,\bar{\tau }]}{\min }g({\tau }_{ϵ}|H_1)<\underset{{\tau }_{ϵ}\in [0,\bar{\tau }]}{\min }g({\tau }_{ϵ}|H_2),$$

Full text

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 $W=20$ MHz, one can localize channel delay components separated by $\Delta\tau=\frac{1}{W}=50$ ns, whereas when the bandwidth increases to $1$ GHz, one can improve the resolution to $\Delta\tau=1$ ns [3]. Multiplying the delay resolution with the speed of the electromagnetic wave, results in a ranging error of up to $15$ m in the former versus up to only $30$ 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 Fog-Connection 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 $20$ 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 multi-GHz, making WiFi CSI measurements over the wide bandwidth a useful resource for decimeter-level localization tasks. The process of exploiting such measurements over multiple bands is referred to as multi-band splicing [3]. Despite the potential of multi-band splicing, channel impulse response (CIR) estimation using raw WiFi CSI measurements is challenging, mainly due to imperfections of the underlying hardware. Per-band 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 multi-band 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 signal-to-noise 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 multi-band splicing with phase-distorted 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 conjugate-reflections 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 $M$ we denote by $[M]$ the set

$\{0,\ldots,M-1\}$

. For a vector ${\bf x}$, $|{\bf x}|$ denotes elementwise absolute values. For a function $h(\tau)$, we define the Fourier transform as

${\cal F}\{h(\tau)\}|_{f}=\int h(\tau)e^{-j2\pi\tau f}d\tau$

. For a set of integers $J$ and a vector ${\bf x}$, ${\bf x}_{J}$ denotes the elements of ${\bf x}$ in $J$.

II System Setup

Assume that the propagation environment between a WiFi transmitter and a receiver consists of $K$ scatterers, each imposing a path gain $c_{k}\in{\mathbb{C}}$ and delay $\tau_{k}\in\Xi\overset{\Delta}{=}[0,\tau_{\text{max}}]$, $k\in[K]$ where $\tau_{\text{max}}>0$ denotes the maximum delay spread. Using OFDM signaling, the WiFi device transmits pilots to the receiver over $M$ frequency bands, each consisting of $N$ equi-spaced subcarriers.111Note that the bands need not be adjacent, as is sometimes the case where the WiFi device uses multiple bands around $2$ GHz and $5$ GHz [2]. The receiver collects the CSI of each band separately, by down-conversion through removing the carrier frequency. The baseband representation of the CIR corresponding to band $m\in[M]$ is thereby given as $h^{(m)}(\tau)=\sum_{k=0}^{K-1}c_{k}e^{-j2\pi f_{m,0}\tau_{k}}\delta(\tau-\tau_{k}),$ where $f_{m,0}$ denotes the carrier frequency of the $m$-th band and $\delta(\cdot)$ denotes Dirac’s delta. Since only a few propagation paths contribute to the channel, i.e. $K$ 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 $\{n\,f_{s}:n=-\tfrac{N-1}{2},\ldots,\tfrac{N-1}{2}\}$222$N$ is an odd integer., with $f_{s}$ being the carrier spacing, are denoted by the vector $\tilde{{\bf h}}^{(m)}\in{\mathbb{C}}^{N}$ where

[TABLE]

where $f_{m,n}\overset{\Delta}{=}f_{m,0}+nf_{s}$.

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 $m\in[M]$ by the vector ${\bf y}^{(m)}\in{\mathbb{C}}^{N}$ where each element can be expressed as [2, 3, 8]

[TABLE]

where $\phi_{m,n}\overset{\Delta}{=}2\pi nf_{s}\delta_{m}+\psi_{m},$ is an affine phase distortion term with $\delta_{m}$ denoting the band-specific time offset due to the packet detection delay (PDD) and $\psi_{m}$ denoting the constant phase offset. Also ${\bf z}^{(m)}\in{\mathbb{C}}^{N}$, by assuming that the support333The support of a function $f(x)$ denotes the set $\{x:f(x)\neq 0\}$. 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 $\alpha>0$ is the power control gain. For simplicity, we set $\alpha=1$ 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.

II-A Denoising and Zero-Subcarrier 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 zero-subcarrier of each band, i.e. on the frequencies $f_{m,0},~{}m\in[M]$, since measurements on these subcarriers overlap with the hardware DC offsets that are hard to remove [2]. This means that in each band $m$, the measurement ${\bf y}^{(m)}_{0}$ is not observed. However, the CSI on the zero subcarrier is very useful, since the phase distortion term $\phi_{m,0},~{}m\in[M]$ is free of PDD and reduces to a constant offset, i.e., $\phi_{m,0}=\psi_{m}$. 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 zero-subcarrier.

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 ($K\ll N$), we propose a different interpolation/estimation scheme. Our estimator is based on $\ell_{1}$-norm minimization, which exploits channel sparsity and will additionally denoise the raw CSI measurements. Recalling (2), one can write ${\bf y}^{(m)}=\sum_{k=0}^{K-1}w_{k}^{(m)}{\bf b}^{(m)}(\tau_{k})+{\bf z}^{(m)},$ where ${\bf b}^{(m)}_{n}(\tau_{k})=\exp\left(-j2\pi\left(f_{m,n}\tau_{k}+nf_{s}\delta_{m}\right)\right)$ and $w_{k}^{(m)}=c_{k}e^{-j\psi_{m}}$. Since the signal component of ${\bf y}^{(m)}$ is a linear combination of a small number ($K$) of the vectors ${\bf b}^{(m)}(\tau_{k})$, it is natural to use a sparse recovery method for both estimating ${\bf y}_{0}^{(m)}$ and denoising the observed values of ${\bf y}^{(m)}$. Here we use the basis pursuit denoising (BPDN) method for this purpose. Define the oversampled DFT matrix ${\bf D}\in{\mathbb{C}}^{N\times L}$ where $[{\bf D}]_{n,\ell}=\frac{1}{\sqrt{N}}e^{-j2\pi\frac{n\ell}{L}}~{},n\in[N],\,\ell\in[L]$ and $L\gg N$. This matrix serves as a dictionay over a discretized delay grid with $L$ delays over $[0,1]$. The BPDN program, thereby, can be formulated as follows [9]

[TABLE]

where $\|\cdot\|_{1}$ denotes the $\ell_{1}$ norm, $\rho>0$ is an appropriate regularization scalar and $J=\{-\tfrac{N-1}{2},\ldots,\tfrac{N-1}{2}\}\backslash\{0\}$. In practice $J$ 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 $\tilde{{\bf w}}^{(m)}$, using which we can not only estimate the zero subcarrier CSI as $\left({\bf D}\tilde{{\bf w}}^{(m)}\right)_{0}$, but also denoise the observed vector ${\bf y}^{(m)}_{J}$. We denote the overall denoised CSI vector over band $m\in[M]$ as ${\bf y}^{(m)}_{\sf den}={\bf D}\tilde{{\bf w}}^{(m)}$. 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.

II-B Problem Statement

Let $\tilde{{\bf h}}\overset{\Delta}{=}[(\tilde{{\bf h}}^{(0)})^{\sf T},\ldots,(\tilde{{\bf h}}^{(M-1)})^{\sf T}]^{\sf T}\in{\mathbb{C}}^{MN}$ denote the vector containing the CFR samples over all bands. Using (1) one can easily see that $\tilde{{\bf h}}$ contains the frequency samples of the CIR $h(\tau)=\sum_{k=0}^{K-1}c_{k}\delta(\tau-\tau_{k})$ over the ordered set of subcarrier frequencies $\{f_{0,-\tfrac{N-1}{2}},\ldots,f_{M-1,\tfrac{N-1}{2}}\}$. In the same way we define ${\bf y}\overset{\Delta}{=}[(\tilde{{\bf y}}_{\sf den}^{(0)})^{\sf T},\ldots,(\tilde{{\bf y}}_{\sf den}^{(M-1)})^{\sf T}]^{\sf T}\in{\mathbb{C}}^{MN}$ and ${\bf z}\overset{\Delta}{=}[(\tilde{{\bf z}}^{(0)})^{\sf T},\ldots,(\tilde{{\bf z}}^{(M-1)})^{\sf T}]^{\sf T}\in{\mathbb{C}}^{MN}$ as the denoised measurements vector and the noise vector, respectively. It follows that the CSI measurement vector can be written as ${\bf y}=\boldsymbol{\Phi}\tilde{{\bf h}}+{\bf z},$ where, recalling (2), $\boldsymbol{\Phi}\in{\mathbb{C}}^{MN\times MN}$ is an unknown diagonal matrix, with unit modulus phase-distortion entries. The main problem addressed in this paper is formally stated as follows.

Problem** 1**

Given the distorted CSI measurements vector ${\bf y}$, estimate the sparse CIR $h(\tau)=\sum_{k=0}^{K-1}c_{k}\delta(\tau-\tau_{k})$.

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. $\{|{\bf y}_{i}|^{2}\}_{i\in[MN]}$. An important feature of formulating and solving the problem as such is that, unlike previous methods in the literature, we do not perform any time-consuming pre-processing 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 delay-sensitive applications.

III CIR Estimation via Phase Retrieval

Using the formulation developed in the previous section, define

[TABLE]

where ${\bf u}$ denotes the vector of magnitudes and $\tilde{{\bf z}}_{i}\overset{\Delta}{=}\tilde{{\bf h}}_{i}\boldsymbol{\Phi}_{i,i}^{\ast}{\bf z}_{i}+\tilde{{\bf h}}_{i}^{\ast}\boldsymbol{\Phi}_{i,i}{\bf z}_{i}^{\ast}+|{\bf z}_{i}|^{2}$. It is easy to see that $|\tilde{{\bf h}}_{i}|^{2},i\in[MN]$ represents the samples of the Fourier transform of the CIR autocorrelation function

[TABLE]

for $\xi\in[-\tau_{\text{max}},\tau_{\text{max}}]$, where $\star$ denotes convolution. In order to recover $h(\tau)$ from $\{{\bf u}_{i}\}_{i\in[MN]}$ we first recover $R(\xi)$ and then estimate the set of delays $\{\tau_{k}\}_{k\in[K]}$ and their corresponding coefficients $\{|c_{k}|\}_{k\in[K]}$ (up to a phase constant) from the estimated autocorrelation.

III-A Sparse Recovery of the CIR Autocorrelation

From (5) it follows that the autocorrelation $R(\xi)$ is conjugate symmetric with respect to $\xi=0$. We can write $R(\xi)=r_{0}\delta(\xi)+\sum_{s=1}^{\tilde{K}}r_{s}\delta(\tau-\xi_{s})+\sum_{s=1}^{\tilde{K}}r_{s}^{\ast}\delta(\tau+\xi_{s})$ where $\tilde{K}\overset{\Delta}{=}\frac{K(K-1)}{2}$, and $\xi_{s}\overset{\Delta}{=}|\tau_{k}-\tau_{\ell}|\in\Xi$ for some $k,\ell\in[K]$. Furthermore, $R(\xi)$ is a sparse function, i.e. assuming $K$ to be small, it consists of a few number of delta functions over $\Xi\cup\{-\Xi\}$. Using (4) we can write

[TABLE]

where ${\bf a}(\xi_{s})=[e^{-j2\pi f_{0}\xi_{s}},\ldots,e^{-j2\pi f_{MN-1}\xi_{s}}]^{\sf T}$, ${\bf a}_{\text{re}}(\xi_{s})={\rm Re}\{{\bf a}(\xi_{s})\}$, ${\bf a}_{\text{im}}(\xi_{s})={\rm Im}\{{\bf a}(\xi_{s})\}$ and $\mathbf{1}$ denotes the all ones vector of size $MN$. In order to recover $R(\xi)$, we define a dictionary of the vectors ${\bf a}(\xi)$ on a dense grid of points over $\Xi$ and apply sparse recovery methods. Let ${\cal G}=\{\tfrac{\tau_{\text{max}}}{G},\tfrac{2\tau_{\text{max}}}{G},\ldots,\tfrac{G\tau_{\text{max}}}{G}\}$ with $G\gg MN$ be a dense grid over $\Xi$. Define the dictionary ${\bf A}\overset{\Delta}{=}\left[\mathbf{1},\,2{\bf A}_{1}\,,-2{\bf A}_{2}\right]$ where ${\bf A}_{1}={\rm Re}\left\{\left[{\bf a}(\xi_{0}),\ldots,{\bf a}(\xi_{G-1})\right]\right\}\in{\mathbb{R}}^{MN\times G}$ and ${\bf A}_{2}={\rm Im}\left\{\left[{\bf a}(\xi_{0}),\ldots,{\bf a}(\xi_{G-1})\right]\right\}\in{\mathbb{R}}^{MN\times G},~{}\xi_{i}\in{\cal G}$. We can approximate the signal part of ${\bf u}$ with a linear combination of the columns of ${\bf A}$ and therefore we have $|\tilde{{\bf h}}|^{2}\approx{\bf A}{\bf x},$ where ${\bf x}=[x_{0},\,{\bf x}_{\text{re}}^{\sf T},{\bf x}_{\text{im}}^{\sf T}]^{\sf T}$ is a sparse vector, where ${\bf x}_{\text{re}}\in{\mathbb{R}}^{G}$ and ${\bf x}_{\text{im}}\in{\mathbb{R}}^{G}$ 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 ${\bf x}$ given ${\bf u}$, we solve the following BPDN problem,

[TABLE]

where $\lambda>0$ 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 II-A, 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 $R(\xi)$. Using the notation ${\bf x}^{\star}=[x_{0}^{\star},\,({\bf x}_{\text{re}}^{\star})^{\sf T},({\bf x}_{\text{im}}^{\star})^{\sf T}]^{\sf T}$, the coefficients of this discrete approximation over $(0,\tau_{\text{max}}]$ are given by $\hat{R}(\xi_{i})=[{\bf x}_{\text{re}}^{\star}]_{i}+j[{\bf x}_{\text{im}}^{\star}]_{i}$ for $\xi_{i}\in{\cal G},$ and $\hat{R}(\xi_{0})=x^{\star}_{0}$, whereas the coefficients over $[-\tau_{\text{max}},0)$ are simply given by $\hat{R}(-\xi_{i})=\hat{R}(\xi_{i})^{\ast}$.

This discretized approximation of the autocorrelation usually contains spurious non-zero elements, since the true support entries ${\cal D}=\{\xi_{s}\}_{s=1}^{\tilde{K}}$ do not necessarily lie on the grid ${\cal G}$ and also due to noise effects. This means that the support set of $\hat{R}(\xi)$, i.e. $\{\xi:\hat{R}(\xi)\neq 0\}$ generally has more than $2\tilde{K}+1$ elements. However, by knowing the number of paths $K$ contributing to the CIR,444The number of paths corresponds to the slow-varying geometry of the environment and can be estimated offline during a pre-processing step. we can refine the support of $\hat{R}(\xi)$ as follows. Let ${\cal I}=\{\xi_{i}>0:|\hat{R}(\xi)|\geq\varepsilon\}$ denote the set of points in ${\cal G}$ whose corresponding coefficients in $R(\xi)$ are large enough, where $\varepsilon>0$ is a small, predefined threshold. The points in ${\cal I}$ are clusters around the true support points $\{\xi_{s}\}_{s=1}^{\tilde{K}}$. Knowing the value $\tilde{K}$, we apply the k-means clustering method [11] and obtain the center points of these clusters, denoting them by $\{\hat{\xi}_{s}\}_{s=1}^{\tilde{K}}$. Eventually, the estimated autocorrelation support is given by $\hat{{\cal D}}=\{\hat{\xi}_{0}=0\}\cup\{\pm\hat{\xi}_{s}\}_{s=1}^{\tilde{K}}$. Now, the coefficients corresponding to points in the estimated support set can be obtained by solving a least-squares problem

[TABLE]

where ${\bf A}_{\hat{{\cal D}}}=\left[{\bf a}(-\hat{\xi}_{\tilde{K}}),\ldots,{\bf a}(0),\ldots,{\bf a}(\hat{\xi}_{\tilde{K}})\right]$. In what follows we refer to the estimated coefficient corresponding to $\hat{\xi}_{s},~{}s=1,\ldots,\tilde{K}$, as $\hat{r}_{s}$. Note that the coefficient corresponding to $-\hat{\xi}_{s}$ will be $\hat{r}_{s}^{\ast}$ due to conjugate symmetry. This completes the estimation of the autocorrelation function $R(\xi)$.

III-B CIR Recovery

In order to recover the CIR $h(\tau)$ from its autocorrelation estimate $\hat{R}(\xi)$, we rely on a method developed in [12] which performs support and magnitude recovery in succession. First, note that the support of the autocorrelation ${\cal D}=\{\tau_{k}-\tau_{\ell}\}_{k,\ell\in[K]}$ is the difference set of the support of the CIR ${\cal T}=\{\tau_{k}\}_{k\in[K]}$, i.e. ${\cal D}={\cal T}-{\cal T}$ (Minkowski difference). Also ${\cal D}$ is a symmetric set containing the point [math]. Now we have an estimate of ${\cal D}$ and want to recover ${\cal T}$. A crucial point here is that, recovering ${\cal T}$ from $\hat{{\cal D}}$ involves a fundamental ambiguity even when $\hat{{\cal D}}={\cal D}$: for any constant $b\in{\mathbb{R}}$ the generic sets $b+{\cal T}$ and $b-{\cal T}$ have the same difference set as the ${\cal T}$. Thus, we can hope to recover ${\cal T}$ 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 ${\cal T}$ 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 $(\tau_{k},\tau_{\ell})$ and $(\tau_{k^{\prime}},\tau_{\ell^{\prime}})$ for which $\tau_{k}-\tau_{\ell}=\tau_{k^{\prime}}-\tau_{\ell^{\prime}}$, 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 $\hat{{\cal T}}$, we proceed with estimating the CIR coefficients ${\bf c}=[c_{1},\ldots,c_{k}]^{\sf T}$. 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 ${\bf c}$ and $e^{j\phi}{\bf c}$, result in the same autocorrelation. For any two elements $\hat{\tau}_{k},\hat{\tau}_{\ell}\in\hat{{\cal T}}$ we can find their corresponding differences $\pm\hat{\xi}_{s}=\pm|\tau_{k}-\tau_{\ell}|$ in $\hat{{\cal D}}$ and furthermore, the corresponding estimated coefficients $\hat{r}_{s},~{}\hat{r}_{s}^{\ast}$ from the discussion of Section III-A. From (5), we know that $|r_{s}|=|c_{k}c_{\ell}^{\ast}|$ and since $\hat{r}_{s}$ is an estimate of $r_{s}$ we expect $|\hat{r}_{s}|=|\hat{c}_{k}\hat{c}_{\ell}^{\ast}|$. Now construct the matrix ${\bf C}\in{\mathbb{R}}^{K\times K}$ such that

[TABLE]

Note that $\sum_{\ell=1}^{K}{\bf C}_{k,\ell}=(K-2)\log|\hat{c}_{k}|+\sum_{k=1}^{K}\log|\hat{c}_{\ell}|$ for $k=1,\ldots,K$. Define the sum of the elements of ${\bf C}$ as $\beta\overset{\Delta}{=}\sum_{k}\sum_{\ell}{\bf C}_{k,\ell}=2(K-1)\sum_{k=1}^{K}\log|\hat{c}_{\ell}|$. When $K>2$, using these equations we can obtain $|\hat{c}_{k}|,~{}k\in[K]$ as

[TABLE]

which gives us an estimate of the amplitude of the CIR coefficients.555The solution for $K=1,2$ can be trivially computed and is not discussed here due to a lack of space. Let $\theta_{k}$ denote the phase of $\hat{c}_{k}$, i.e., $\hat{c}_{k}=|\hat{c}_{k}|e^{j\theta_{k}}$. We estimate the CIR coefficients by treating $\hat{c}_{0}$ as real-valued which is of no consequences due to the fundamental phase ambiguity. This vector of coefficients, denoted by $\hat{{\bf c}}=[|\hat{c}_{0}|,|\hat{c}_{1}|e^{j\theta_{1}},\ldots,|\hat{c}_{K-1}|e^{j\theta_{K-1}}]^{\sf T},$ is estimated as follows. For any $k=1,\ldots,K-1$, we simply find $\hat{\xi}_{s}$ such that $\hat{\xi}_{s}=\hat{\tau}_{k}-\hat{\tau}_{0}$. The corresponding estimated autocorrelation coefficient for $\hat{\xi}_{i}$ is given by $\hat{r}_{s}\approx\hat{c}_{k}\hat{c}_{0}^{\ast}$. By dividing $\hat{r}_{s}$ to $|\hat{c}_{0}|$ we obtain $|\hat{c}_{k}|e^{j\theta_{k}}$.

III-C 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 II-A, the zero subcarrier in band $m\in[M]$ is free from the PDP phase error and contains only the constant phase error term $\psi_{m}$. 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.,

[TABLE]

where ${\bf y}_{0,tx}^{(m)}$ and ${\bf y}_{0,rx}^{(m)}$ denote CSI on the zero-subcarrier of band $m$ at the transmitter and the receiver and ${\bf z}_{0,tx}^{(m)}$ and ${\bf z}_{0,rx}^{(m)}$ denote their corresponding noise terms, respectively. The CFR sample $\tilde{{\bf h}}^{(m)}_{0}$ 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

[TABLE]

where ${\bf z}_{m}^{\prime}$ denotes the noise and cross terms. The vector ${\bf y}^{\prime}=[{\bf y}_{0}^{\prime},\ldots,{\bf y}^{\prime}_{M-1}]^{\sf T}\in{\mathbb{C}}^{M}$ gives us the extra information with which we resolve the ambiguity about the shift and reflection in the support set $\hat{{\cal T}}$. This is done by constructing a simple hypothesis test. Let the elements of $\hat{{\cal T}}$ be ordered as $\hat{\tau}_{0}(=0)<\ldots<\hat{\tau}_{K-1}$. In addition, we know that in an indoor environment the first delay component is within a certain limits. Let $\bar{\tau}>0$ 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 $H_{1}$ be denoted by $f_{+}(\tau;\tau_{\epsilon})=\sum_{k=0}^{K-1}\hat{c}_{k}\delta(\tau-\hat{\tau}_{k}-\tau_{\epsilon}),$ where by construction we made the first estimated delay component to appear at zero, i.e. $\hat{\tau}_{0}=0$ (see Algorithm 1, line 7). Also define the conjugate-reflected and shifted version of $f_{+}(\tau;\tau_{\epsilon})$ to be $f_{-}(\tau;\tau_{\epsilon})=\sum_{k=0}^{K-1}\hat{c}_{k}^{\ast}\delta(\tau+\hat{\tau}_{k}-\hat{\tau}_{K-1}-\tau_{\epsilon}),$ which represents the CIR of the second hypothesis $H_{2}$. Define the Fourier transform of the shifted CIRs corresponding to $H_{1}$ and $H_{2}$ by a value $\tau_{\epsilon}$ as $F_{+}(f;\tau_{\epsilon})={\cal F}\{f_{+}(\tau;\tau_{\epsilon})\}|_{f}$ and $F_{-}(f;\tau_{\epsilon})={\cal F}\{f_{-}(\tau;\tau_{\epsilon})\}|_{f}$, 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 $\tau_{\epsilon}$ gives the true CIR by comparing these samples with the ones in ${\bf y}^{\prime}$. Let ${\bf p}_{+}(\tau_{\epsilon})=[F_{+}^{2}(f_{0,0};\tau_{\epsilon}),\ldots,F_{+}^{2}(f_{M-1,0};\tau_{\epsilon})]^{\sf T}$ and ${\bf p}_{-}(\tau_{\epsilon})=[F_{-}^{2}(f_{0,0};\tau_{\epsilon}),\ldots,F_{-}^{2}(f_{M-1,0};\tau_{\epsilon})]^{\sf T}$ denote the squared frequency samples of the tentative solutions for a shift $\tau_{\epsilon}$. For the first hypothesis $H_{1}$, i.e. the CIR being only shifted, form the cost function $g(\tau_{\epsilon}|H_{1})=\|{\bf p}_{+}(\tau_{\epsilon})-{\bf y}^{\prime}\|_{2}$. In a similar way, for the second hypothesis $H_{2}$, i.e. the CIR being conjugate-reflected define the cost function as $g(\tau_{\epsilon}|H_{2})=\|{\bf p}_{-}(\tau_{\epsilon})-{\bf y}^{\prime}\|_{2}$. We select $H_{1}$ over $H_{2}$ if

[TABLE]

and $H_{2}$ over $H_{1}$ otherwise. In addition, the optimal value of the shift parameter $\tau_{\epsilon}$ is given by $\tau_{\epsilon}^{\star}=\underset{\tau_{\epsilon}}{{\hbox{arg}}\min}~{}g(\tau_{\epsilon}|H_{i}),$ where $H_{i}$ is the winning hypothesis. After obtaining $\tau_{\epsilon}^{\star}$, the eventual set of estimated delays is given by $\hat{{\cal T}}+\tau_{\epsilon}^{\star}$. This completes the CIR estimation process and $\tau_{\epsilon}^{\star}$ 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 zero-subcarrier CSI is exchanged between the WiFi device and the client to get the vector ${\bf y}^{\prime}$ in (12). This vector represents the frequency samples of the convolution of the CIR with itself. A sparse recovery method is applied to ${\bf y}^{\prime}$ to get a discretized estimate of the convolution. The first non-zero 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 $K=3$ 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 $M=32$ adjacent WiFi bands, each with $N=33$ subcarriers, with a carrier spacing of $f_{s}=312.5$ kHz. We set the parameters in programs (3) and (7) as $L=3N$ and $G=3MN$, respectively. Appropriate regularization scalars $\rho$ and $\lambda$ 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 $\hat{\tau}_{0}$, the ranging error is computed as $e_{d}=|\tau_{0}-\hat{\tau}_{0}|C$ (in meters), where $C$ denotes the speed of light. We run the experiment for $1000$ Monte-Carlo trials. The ranging error CDF is illustrated in Fig. 1 for two SNR values of $10$ and $20$ dB. As we can see, our proposed method outperforms Chronos in about $94\%$ of the trials, and has a ranging error less than $10$ cm in about $90\%$ of the trials, indicating a cm-level localization accuracy. Notice that even when SNR=10 dB, the PR-based method performs better than Chronos fed with CSI measurements with SNR=20 dB in about $92\%$ 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 multi-band 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 low-SNR regimes.