The Time-of-Arrival Offset Estimation in Neural Network Atomic Denoising in Wireless Location

Abstract

With the increasing demand for wireless location services, it is of great interest to reduce the deployment cost of positioning systems. For this reason, indoor positioning based on WiFi has attracted great attention. Compared with the received signal strength indicator (RSSI), channel state information (CSI) captures the radio propagation environment more accurately. However, it is necessary to take signal bandwidth, interferences, noises, and other factors into account for accurate CSI-based positioning. In this paper, we propose a novel dictionary filtering method that uses the direct weight determination method of a neural network to denoise the dictionary and uses compressive sensing (CS) to extract the channel impulse response (CIR). A high-precision time-of-arrival (TOA) is then estimated by peak search. A median value filtering algorithm is used to locate target devices based on the time-difference-of-arrival (TDOA) technique. We demonstrate the superior performance of the proposed scheme experimentally, using data collected with a WiFi positioning testbed. Compared with the fingerprint location method, the proposed location method does not require a site survey in advance and therefore enables a fast system deployment.

1. Introduction

Due to the emergence of 5G and 6G technology, people have higher and higher requirements for indoor positioning [1,2]. According to Market&Markets research data, the global market size related to indoor location is expected to reach USD 18.74 billion by 2025 [3]. The huge market demand has become a key driver of location-based applications in the future. However, despite the growing demand in a large market, it is not easy to provide a viable indoor positioning solution in many cases. One of the key challenges facing distance-based positioning systems is that the reflection and refraction of wireless signals in indoor environments may lead to a serious decrease in positioning accuracy.

At present, in addition to satellite-based positioning technology [1,2], various positioning technologies are available, such as Bluetooth [4,5], radar [6,7], radio frequency identification [8], ultra-wideband [9], geomagnetic field [10], visible light [11], thermal infrared [12], sound [13], 5G cellular networks [3], and so on. Most of these positioning technologies, however, require additional hardware anchor points and therefore incur additional costs. With the widespread deployment of WiFi systems, it is advantageous to exploit the ubiquitous WiFi [14,15,16] for positioning, which provides broader convenience for the promotion of technology and business cases.

There are currently two CSI-based positioning technologies, which are based on fingerprint and triangulation, respectively. The fingerprint-based localization method needs to collect CSI data at various locations of the deployment site in advance. The FILA approach proposed by Kaishun et al. [17] achieves a 50 percent accuracy at less than 0.6 m in the laboratory. The C-Map proposed by WEN et al. [18] achieves an average positioning error in the comprehensive environment of 0.97 m. CHAUR-HEH et al. [19] adopted a convolutional neural network to localization in a room, which performs the best, resulting in a maximal localization error of 0.92 m with a probability of 99.97%. YUAN et al. In [20], a multi-view discriminant learning approach was developed for indoor localization that exploits both the amplitude and the phase information of CSI to create feature images for each location, and the minimum distance errors for the laboratory and corridor experiments were 0.205 m and 0.109 m, respectively. In [21], the authors proposed to transform the measured data from CSI into images, extract their features, and use deep learning networks to estimate localization, and achieved an accuracy of more than 90% in laboratories. In [22], the authors proposed to transform the measured data of CSI into images and use the classification ability of the convolutional neural network (CNN) for localization, and 70% of the test cases had a localization error under 1.5 m. The fingerprinting-based approach achieves a relatively high positioning accuracy, but requires the deployment site to be surveyed in advance and whenever the environment is changed, which increases the operating costs of commercial applications.

Data fusion-based positioning can achieve a higher positioning accuracy. Zhao et al. [23] proposed a data fusion method of fingerprinting using RSS and CSI data from single access points, which can achieve a positioning accuracy of 1.79 m in a typical laboratory. Li et al. [24] proposed an enhanced particle-filter-based positioning method that combines CSI information and inertial sensor information to achieve an average accuracy of 1.3 m. However, such techniques require more data sources and therefore are not favorable for commercial applications.

Triangulation-based positioning uses angle and/or distance measurements to locate targets. Unlike fingerprinting-based approaches, it does not require data collection in advance, which brings great convenience to commercial deployment. The systems only needs to be updated when the anchor points are moved, which involves surveying the coordinates of anchors only. In [25], Manikanta et al. used multiple measurements of the angle-of-arrival to achieve positioning on standard WiFi equipment, and 60% of the localization error reached approximately one-meter in indoor office deployment. In [14], WiFi devices were located by using the round-trip delay and angle of arrival (AOA) measurements, which require changes to the firmware of target devices. The multiple signal classification (MUSIC) algorithm proposed by Schmidt et al. [26] had a good effect in terms of resolution, estimation accuracy, and stability under the condition of multiple antennas. ArrayTrack [27] proposed by Jie et al. is similar to [26], requiring eight antennas.

This paper proposes a novel algorithm to estimate the TOA from CSI, which enables an accurate localization of standard WiFi devices. The main contributions of this work are as follows:

• A novel algorithm based on the CS is proposed to analyze the CSI and extract the CIR.
• A dictionary denoising method based on the direct weight determination of a neural network is derived for denoising CIR signals in compressed sensing. The sparsity of a dictionary is used as the sparse matrix of compressed sensing to estimate the TOA, which improves the accuracy of estimation.
• The proposed method is validated experimentally and shown to outperform conventional approaches based on MUSIC or estimating signal parameters via rotational invariance techniques (ESPRIT).

The remainder of this paper is organized as follows. The system model, including system architecture and signal structure, is presented in Section 2. The problem-solving process is presented in Section 3. We present the process of solving TOA by the signal reconstruction proposed. Experimental validation is provided in Section 4, followed by concluding remarks in Section 5.

2. System Model

In this section, we briefly introduce the localization testbed that we used, and then introduce the experiment data collected with the system.

2.1. System Architecture

In this paper, as shown in Figure 1, we used the WiFi-based wireless ad hoc system for positioning (WiFi-WASP) [2,28,29] developed by the Commonwealth Scientific and Industrial Research Organization (CSIRO) of Australia for our experiments. We installed six custom-built WiFi sniffers as anchors for localization in a conference hall and placed an access point (AP) approximately in the middle of the conference hall. The sniffer measures the time stamp of the communication between the sniffers and the AP to estimate the clock offset and offset of the system clock. The sniffers also measure the time stamp of the communication between the target and the AP to locate the target. A computer was placed in the conference hall to collect time stamps collected by sniffers and estimate the location of target devices by using the TDOA. The system is passive and does not interfere with existing standard WiFi systems.

In WiFi-WASP systems, the sniffer clock is post-synchronized by estimating and compensating for clock skew and clock offset to estimate the arrival time of packets based on a common reference clock [30]. This indicates that the locations of the sniffers and the target device is ${\mathit{s}}_j={[x_j,y_j,z_j]}^T(j=1,\dots ,N)$, where $\mathit{c}v={[x,y,z]}^T$ is the coordinate value of the device and N is the number of sniffers. Each sniffer in the system has a local clock that is not synchronized with each other or with any external references. Assuming that the target device transmits a packet at $t_{Tx}$, the arrival time of the corresponding radio signal measured by each sniffer is

$$r_j=(1+{\alpha }_j)(t_{Tx}+\frac{d_j}{c}+{\beta }_j+\Delta t_j+m_j),j=1,\dots ,J,$$

(1)

where ${\alpha }_j$ and ${\beta }_j$ denote the clock skew and clock offset of sniffer j, $d_j\triangleq |{\mathit{s}}_j-\mathit{c}v|$ denotes the distance between the target device and sniffer j, c is the speed of light, $\delta t_j$ is the hardware delay (e.g., delay caused by radio frequency (RF) circuitry), and $m_j$ is the time measurement error. During system operation, the measured packet arrival time is corrected according to the relevant AGC settings. Since the values of $\alpha$ and $\Delta t$ are small, in order to directly subtract excessive hardware delay from $r_j$, Equation (1) can be rewritten as

$$r_j=\Delta t_j+(1+{\alpha }_j)(t_{Tx}+\frac{d_j}{c}+{\beta }_j+m_j),j=1,\dots ,J,$$

(2)

By estimating the clock skew ${\alpha }_j$, clock offset ${\beta }_j$, and the hardware delay $\Delta t_j$ of sniffer j, as well as the measurement value of $r_j$, the measurement value $t{r}_j$ of the synchronous arrival time of sniffer j (less than 0.5 s) in a small time can be estimated, which can be expressed as

$$t{r}_j=\frac{r_j-\Delta t_j}{1+{\alpha }_j}-{\beta }_j-m_j,$$

(3)

Suppose a standard WiFi device transmits two packets in succession at $t_n$ and $t_{n+1}$, respectively. ${\alpha }_j$ can be written approximately as

$${\alpha }_j-{\alpha }_k\approx \frac{r_j^{n+1}-r_j^n}{r_k^{n+1}-r_k^n}-1$$

(4)

where j and k are different sniffers. By assigning one of the sniffers as the reference clock (i.e., ${\alpha }_j$ = 0), the clock skews of other sniffers can be obtained from (4).

Since the locations of the sniffers are also known, the distance between the transmitter and the sniffer can be readily obtained. Thus, ${\beta }_j$ can be written approximately as

$${\beta }_j-{\beta }_k\approx \frac{r_j}{1+{\widehat{\alpha }}_j}-\frac{r_k}{1+{\widehat{\alpha }}_k}-\frac{d_j}{c}+\frac{d_k}{c}$$

(5)

where j and k represent different sniffers and $\widehat{\alpha }$ represents the clock skew estimated by (4).

2.2. Signal Structure

The WiFi-WASP system equipment used the orthogonal frequency-division multiplexing (OFDM) modulation scheme. Figure 2 illustrates the system model for obtaining the TOA calculation considered in this paper from OFDM [31]. At the receiving end, the analog data collected by the antenna went through the processes of the analog/digital converter, the serial transmission data were then converted to parallel data, circular prefixes were removed, discrete Fourier transform was conducted [32], and, finally, the frequency domain CSIs were obtained for localization. The CSI in the frequency domain is given by

$$H\left(f\right)=\sum _{k=0}^{L-1}{\alpha }_{k}exp(-j2\pi (f_0+n△f){\tau }_k),k\in [1,L],$$

(6)

where $f_0$ is the central frequency and $△f$ represents the frequency space between adjacent subcarriers. ${\alpha }_k$ and ${\tau }_k$ denote the complex Gaussian channel coefficient and the time delay of the k-th signal propagation path, respectively. L is the number of multipath components. Since the system adopts an 80 MHz bandwidth, there were 256 subcarriers in total, and the subcarrier space was 80 MHz/256 = 312.5 khz. The number of CSI values in the frequency domain was 245 because 11 subcarriers were unused. f denotes the center frequency of each subcarrier. We denote $\tau =\{{\tau }_0,{\tau }_1,\dots ,{\tau }_{L-1}\}$ as the set of true L time delays.

Multipath propagation characterized by the channel frequency response (CFR) in the frequency domain signals includes the amplitude–frequency response and phase–frequency response, while, in wireless channels, amplitude attenuation, phase offset, and time delay can be characterized by CIR. In the CS method, assuming that the sparse matrix is an identity matrix and there is only a rotation factor in the sparse matrix, the multiplication of the sparse matrix by the time domain signal is equivalent to performing the Fourier transform. The physical meaning of the sparse matrix is to transform the signal to be sampled into another domain, and, in this domain, the signal is sparse. Thus, we used the CS method to convert the CFR into the CIR in the time domain to estimate the TOA, and we obtained

$$h\left(\tau \right)=\sum _{k=0}^{L-1}{\alpha }_k\delta (\tau -{\tau }_k),k\in [1,L],$$

(7)

where $a_k$ and ${\tau }_k$ denote the amplitude of the multipath component and the complex attenuation and propagation delay of the kth path, respectively. $\delta \left(\tau \right)$ is the Dirac delta function and L is the number of multipath components, while $1\le k\le L$ are in ascending order. In a line-of-sight (LOS) environment, the first peak of the impulse response is known to be a good choice for the TOA of the direct path.

3. Proposed Approach

In this section, the direct weight determination method based on a neural network is introduced to filter dictionary atoms, and the filtering effect is proven by using the quintic polynomial. In order to increase the resolution of TOA, an augmented matrix was added for compressed sensing, and the sparsity and coherence of the proposed method are proven.

Figure 3 illustrates the steps taken to locate a target based on CSI collected by the considered WiFi positioning testbed. A Hamming window was first applied to the CSI to prevent frequency leakage. In order to estimate high-precision TOA, we used neural network weight determination method to filter dictionary atoms. The phase space reconstruction theory was used to reconstruct the input signal to address the problem of insufficient input training signals. The sparsity of dictionary was used as the sparse matrix of compressed sensing, and the resolution of TOA was improved by adding and expanding matrix. However, in a non-LOS or multipath environment, the first arrival path tends to decay significantly and its signal strength may be lower than the reflection path. A peak filtering algorithm was employed to improve the accuracy of estimated TOA. The system clocks among anchors were then synchronized and used to correct the measured TOA values. Finally, the target was located based on the differences between the TOA measured by different anchors.

3.1. Frequency Domain Windowing

In the process of digital signal processing, we can only transform the limited length of signal data, so it is necessary to carry out signal truncation. Even if it is a periodic signal, if the truncation time is not an integer multiple of the period, then the intercepted signal will have frequency leakage. To minimize this leakage error, we need to use a weighting function, also known as a window function. There are many kinds of window functions, such as rectangular window, Hann window, and Hamming window. Rectangle window only requires an accurate reading of the main lobe frequency, regardless of amplitude accuracy. Hann window is superior to rectangular window, but the main lobe of Hann window is widened, which is equivalent to the analysis bandwidth broadening, and the frequency resolution decreases. If the test signal has multiple frequency components, the spectrum performance is very complex, and the purpose of the test is more concerned with the frequency points rather than the size of the energy, so Hann window should be selected. Both Hamming window and Hann window are cosine windows, also known as improved ascending cosine windows, but the weighting coefficient is different, which makes the sidelobe smaller, and the sidelobe attenuation speed of the Hamming window is slower than that of Hann window. Due to the fact that the amplitude–frequency characteristic of Hamming window is that the side lobe attenuation is large [33], the peak attenuation of the main lobe and the first side lobe can reach 40 dB. Thus, we have

$$w\left(l\right)={\alpha }_0-(1-{\alpha }_0)·cos\left(\frac{2\pi l}{L-1}\right),0\le l\le L-1,$$

(8)

where ${\alpha }_0$ = 25/46, and this value is intended to generate zero-crossing at the frequency $5\pi /(L-1)$. Since the sample point is 256, L = 256. Therefore, the CFR data after applying the Hamming window is expressed as

$$F=H·\mathbf{w},$$

(9)

where $H={[h_1,h_2,\cdots ,h_l,\cdots ,h_L]}^T\in \mathbb{C}$ is the frequency domain CSI collected by WiFi-WASP system, $\mathbf{w}=[w_0,w_1,\cdots ,w_l,\cdots ,w_{L-1}]\in \mathbb{R}$ is the Hamming window, and $F={[s_1,s_2,\cdots ,s_k,\cdots ,s_K]}^T\in \mathbb{C}$.

3.2. Dictionary Atomic Filtering

In order to solve the problem of insufficient data of neural network training, the phase space reconstruction theory is proposed to reconstruct signal space. A direct weight determination method based on neural network is proposed to realize dictionary atomic filtering, and its effectiveness is verified.

3.2.1. The Phase Space Reconstructs the Input Signal

The theory of phase space reconstruction holds that the prime mover system model can be reconstructed from an observation of the system. The evolution of any component in the system is determined by the other components interacting with it, so the information of these related components implies that, in the development process of any component, a spatial system model of the original signal can be reconstructed with an observation of the system.

Due to the fact that each set of signals sampled is different, there are too little data to train the dictionary. According to the phase space reconstruction theory [34,35], for the noise signal $F\in {\mathbb{C}}^{K\times 1}$ sampled from WiFi-WASP system and through the Hamming window, it can be embedded into the attritor orbit matrix $\mathit{F}$ in dimension K so as to solve the problem of insufficient training data $\mathit{F}$.

$$\mathit{F}=\left[\begin{array}{cccc}{s}_1,& s_2,& \cdots ,& s_{K-N+1}\\ s_2,& s_3,& \cdots ,& s_{K-N+2}\\ ⋮& ⋮& \ddots & ⋮\\ s_N,& s_{N+1},& \cdots ,& s_K\end{array}\right],$$

(10)

where $\mathit{F}\in {\mathbb{C}}^{N\times K}$.

3.2.2. Dictionary Learning Algorithms

The dictionary learning model is widely concerned and has been successfully applied in various fields of signal processing, such as signal processing, image processing, image fusion, video key frame extraction, and so on. In the field of signal processing, sparse representation of signals is a simple and effective signal coding method, which describes signals by using linear combinations of as few basis vectors as possible in the basis space. The span of basis vectors is also called a dictionary, and the basis vectors in the corresponding space are atoms.

The dictionary learning algorithm is also very suitable for signal recovery and reconstruction and image noise and blur removal applications. For the dictionary learning algorithm of signal denoising, the trained data can be divided into two categories [36]. First, the method that uses data without noise for dictionary learning is called the external prior dictionary learning algorithm. Second, the method of using noisy data for dictionary learning is called the internal prior dictionary learning algorithm. The external prior dictionary learning algorithm is not adaptive to signal denoising, and it may not be able to reconstruct some fine-scale noise signals well. Compared with the external prior dictionary learning algorithm, the internal prior dictionary learning algorithm adopts the signal data containing noise as the training data; although it has good adaptability, the prior learning is greatly affected by noise. The denoising effect of using such a dictionary to denoise signals is inevitably lower than that of using a fixed atom dictionary. Therefore, reducing the atomic noise of a learning dictionary will definitely improve the denoising effect of signals. The main idea of dictionary learning is to use the dictionary matrix $\Psi \in {\mathbb{R}}^{N*K}$, which contains K atoms ${\mathit{\psi }}_{\mathit{k}}$ and a sparse linear representation of the original sample $\mathit{Y}\in {\mathbb{R}}^{N*L}$, where m represents the number of samples, n represents the attributes of samples, and, in the ideal case, $\mathit{Y}=\Psi \mathit{X}$. This can be transformed into an optimization problem.

$$\begin{array}{cc}\hfill & \underset{\Psi ,\mathit{X}}{min}{\parallel \mathit{Y}-\Psi \mathit{X}\parallel }_F^2,\hfill \\ \hfill & s.t.\forall i,\parallel {\mathit{X}}_i{\parallel }_0\le T_0,\hfill \end{array}$$

(11)

where $\mathit{X}\in {\mathbb{R}}^{K\times L}$ is the sparse decomposition coefficient of the signal, and ${\mathit{X}}_i$ is the row vector in the matrix. $\parallel {\mathit{X}}_i{\parallel }_0$ represents the zero-order norm, which represents the number of non-zero numbers in a vector. $T_0$ represents the maximum sparsity of $\parallel {\mathit{X}}_i{\parallel }_0$. There are two optimization variables $\Psi$ and $\mathit{X}$ to solve this optimization problem. It is common to fix one optimization variable, optimize the other variable, and so on. The sparse matrix $\mathit{X}$ in Formula (11) is solved by the least absolute shrinkage and selection operator classical algorithm.

3.2.3. Direct Weight Determination of Neural Networks

Neural networks with their massively parallel processing, distributed storage, and elastic topology, and other significant characteristics of highly redundant and nonlinear arithmetic, such as the control and nonlinear signal processing, pattern recognition, and robots, have been widely used in such fields as [37]. According to the model of the artificial neural network, the structure can be divided into a feedforward network (also known as a multilayer perceptron network) and the feedback network (also known as the Hopfield net). Among neural network models, the back propagation (BP) neural network model is one of the most widely used neural network models at present. Its core is the error back propagation algorithm, and its model is shown in Figure 4.

The model shows a three-layer BP neural network model that includes an input layer, hidden layer, and output layer. The input layer contains N neurons, so the corresponding input ${\mathit{I}}_i=[x_1,x_2,\cdots ,x_n,\cdots ,x_N]$ is an N dimension vector. There are L neurons in the hidden layer and the excitation function is denoted as $f_1$. In addition, the expected output ${\mathit{O}}_i$ is an $\mathit{I}$ dimensional vector and ${\omega }_{nl}$ is the connection weight between the n neuron of the input layer and the L neuron of the hidden layer, whereas ${\varpi }_{lj}$ is used to represent the connection weight between the l neuron of the hidden layer and the j neuron of the output layer. ${\vartheta }_l$ and ${\theta }_j$ are used to represent the threshold value of the neuron at the hidden layer and the neuron at the output layer, respectively.

${\nu }_1$ denotes the first l hidden layer neurons of the output, and $y_j$ represents the actual output of the j neuron in the output layer, ${\mathit{O}}_i=[O_{i1},O_{i2},\cdots ,O_{ij},\cdots ,O_{IJ}]$ for the sample output. According to the direction of negative gradient descent, the BP algorithm iteratively adjusts the weights and thresholds of the network to achieve the reduction in the error function value, and its weight ${\omega }_{nl}$ and ${\varpi }_{lj}$ correction iterative formula can be described as

$$\left\{\begin{array}{cc}\hfill & {\omega }_{nl}⟵{\omega }_{nl}+\Delta {\omega }_{nl}\hfill \\ \hfill & {\varpi }_{lj}⟵{\varpi }_{lj}+\Delta {\varpi }_{lj},\hfill \end{array}\right.$$

(12)

where

$$\left\{\begin{array}{cc}\hfill \Delta {\varpi }_{lj}& =-\eta \frac{\partial E}{\partial {\varpi }_{lj}}=-\eta \frac{\partial E}{\partial y_j}·\frac{\partial y_j}{{\varpi }_{lj}}=\eta (O_{ij}-y_j)·f_2^{\prime }·v_l\hfill \\ \hfill \Delta {\omega }_{nl}& =-\eta \frac{\partial E}{\partial {\omega }_{lj}}=-\eta \frac{\partial E}{\partial y_j}·\frac{\partial y_j}{\partial v_l}·\frac{\partial v_l}{\partial {\omega }_l}\hfill \\ \hfill & =\eta \sum _{j=1}^J(O_{ij}-y_j)·f_2^{\prime }·{\varpi }_{lj}·f_1^{\prime }·x_n.\hfill \end{array}\right.$$

(13)

In the above formula, $v_l$ represents the output of the lth neuron, $y_j$ represents the output of the jth neuron, and E represents the neural network learning error function corresponding to the jth sample pair under the incremental processing mode. Their derivation formulas are as follows:

$$\left\{\begin{array}{cc}\hfill & v_l=f_1\left(\sum _{n=1}^N{\omega }_{nl}{x}_n-{\vartheta }_l\right),l=1,2,\cdots ,L\hfill \\ \hfill & y_j=f_2\left(\sum _{l=1}^L{\varpi }_{lj}{v}_l-{\theta }_j\right),j=1,2,\cdots ,J\hfill \\ \hfill & E=\frac{1}{2}\sum _{j=1}^J{(y_j-O_{ij})}^2,\hfill \end{array}\right.$$

(14)

The BP network has some problems, such as a slow convergence, local minimization, and fixed learning rate. To solve these problems, Zhang yu-nong et al. proposed an improved power-excited forward neural network model by using polynomial interpolation and approximation theory [38]. Different from the traditional idea of obtaining neural network weights through lengthy BP iterative training, the direct weight determination method can directly obtain the weights, greatly shortening the determination time of network weights. It overcomes the inherent defect that the traditional iterative training algorithm is susceptible to a dynamic learning process and that it is difficult to determine the optimal topology of the neural network. This model can still be regarded as a BP neural network model, as shown in Figure 5. In this figure, the number of hidden layer neurons is n, and the connection weight between hidden layer neurons and output layer neurons is denoted as ${\varpi }_j,j=0,1\cdots ,n-1$, but the relation between the input xand the output y of the neural network can be obtained as $y={\varpi }_0{x}^0+{\varpi }_1{x}^1+{\varpi }_2{x}^2+\cdots +{\varpi }_{n-1}{x}^{n-1}$.

In the figure, the sample set $(x_i,y_i),i=1,2,\cdots ,n$ is used to train the power-excited forward neural network, and the learning error function E is defined as $E=\frac{1}{2}{\displaystyle \sum _{i=1}^n}(y_i-{\displaystyle \sum _{p=0}^{n-1}}{\varpi }_p{x}_i^p)$. Therefore, the weight iterative formula can be expressed as:

$$\begin{array}{cc}\hfill {\varpi }_j(k+1)& ={\varpi }_j\left(k\right)-\eta \frac{\partial E}{\partial {\varpi }_j}\hfill \\ \hfill & ={\varpi }_j\left(k\right)-\eta \sum _{i=1}^m{x}_i^j\left[\left(\sum _{p=0}^{n-1}{\varpi }_p{x}_i^p-y_i\right)\right],\hfill \end{array}$$

(15)

where $j=[1,2,\cdots ,n-1]$. Let $\varpi ={[{\varpi }_0,{\varpi }_1,\cdots ,{\varpi }_{n-1}]}^T\in {\mathbb{R}}^n$ and $y={[y_0,y_1,\cdots ,y_{n-1}]}^T\in {\mathbb{R}}^m$, respectively, where the superscript T expresses the matrix transpose operation. Let $\mathit{X}$ for

$$\mathit{X}=\left(\begin{array}{cc}\hfill x_1^0& x_1^1{x}_1^2\cdots x_1^{n-1}\hfill \\ \hfill x_2^0& x_2^1{x}_2^2\cdots x_2^{n-1}\hfill \\ \hfill & \ddots ⋮\hfill \\ \hfill x_m^0& x_m^1{x}_m^2\cdots x_m^{n-1}\hfill \end{array}\right)\in {\mathbb{R}}^{m\times n}.$$

(16)

Thus, Formula (15) can be expressed in matrix vector form as follows:

$$\begin{array}{c}\hfill \mathit{\varpi }(k+1)=\mathit{\varpi }\left(k\right)-\eta {\mathit{X}}^T[\mathit{X}\mathit{\varpi }\left(k\right)-y].\end{array}$$

(17)

For Formula (15), when the network training reaches the steady state, we have

$$\begin{array}{c}\hfill \underset{k\to +\infty }{lim}\mathit{\varpi }(k+1)=\underset{k\to +\infty }{lim}\mathit{\varpi }\left(k\right)=\mathit{\varpi },\end{array}$$

(18)

Then, there is $-\eta {\mathit{X}}^T[\mathit{X}\mathit{w}-\mathit{y}]=0$. The optimal weight of the power-excited neural network can be directly determined as $\mathit{\varpi }$ by using the matrix pseudo-inverse idea as follows:

$$\begin{array}{c}\hfill \mathit{\varpi }\left(k\right)={\left({\mathit{X}}^T\mathit{X}\right)}^{-1}{\mathit{X}}^T\mathit{y}.\end{array}$$

(19)

It can be seen from the above that the direct weight determination method can automatically, quickly, effectively, and accurately determine the optimal weight, so as to achieve network performance optimization.

3.2.4. Dictionary Atomic Denoising Verification

Regardless of using the recursive least squares dictionary learning algorithm (RLS-DLA), K-singular value decomposition (K-SVD)K-SVD, or other dictionary learning algorithms, as long as the training data contain noise, the resulting dictionary is bound to contain noise. In order to prove that the dictionary atom trained by RLS-DLA contains noise and verify the effect of dictionary learning by the direct weight determination method, we set the real quintic polynomial phase signal as the input signal, which can be expressed as

$$\begin{array}{c}\hfill y\left(n\right)=A\left(n\right)sin\left(\sum _{i=0}^5{a}_i{(\Delta n)}^i\right)+v\left(n\right),n\in [0,N],\end{array}$$

(20)

where $a\left(n\right)=1$. $v\left(n\right)$ is Gaussian white noise, and the signal-to-noise ratio was set to 10 dB. In order to avoid the ambiguity of $a_i$, it follows $|a_i|\le \frac{\pi }{m!n^{m-1}{\Delta }^m}$, so we have $a=(a_5,a_4,a_3,a_2,a_1,a_0)=(4e-10,4e-9,2.75e-6,1.5e-3,\pi /8,\pi /3)$.

According to the algorithm idea of a neural network, a dictionary learning algorithm for denoising polynomial phase signals based on the direct weight determination method is proposed. $\mathit{D}\in {\mathbb{R}}^{N\times L}$ initializes the dictionary. Each column in the dictionary represents an atom, and L represents the number of atoms in the dictionary. Since a redundant dictionary is used, $L\gg N$. $\mathit{Y}\in {\mathbb{R}}^{N\times K}$ represents the training set, where $K\gg L$; thus, the coefficient matrix $\mathit{W}\in {\mathbb{R}}^{L\times K}$. Since the denoising signal itself is used as training data, the observation length of the signal needs to be considered. In the processing of polynomial phase signals, the signal length is not very long, and the number of trained signals K in the training signal set $\mathit{Y}$ is much larger than the number of atoms in the redundant dictionary. Therefore, if the length of the object signal $y\left(i\right)(i=1,2,\cdots ,M)$ is not long enough, in order to satisfy $K\gg L$, the training signal set $\mathit{Y}$ is reconstructed by referring to the phase space reconstruction theory, which can be expressed as

$$\mathit{Y}=\left[\begin{array}{cccc}{y}_1,& y_2,& \cdots ,& y_{N-K+1}\\ y_2,& y_3,& \cdots ,& y_{N-K+2}\\ ⋮& ⋮& \ddots & ⋮\\ y_K,& y_{K+1},& \cdots ,& y_N\end{array}\right].$$

(21)

At the same time, any observed signal in segment K is taken as the initial dictionary $\mathit{D}$, and the ith atom ${\mathit{d}}_i={[d_{i1},d_{i2},\cdots ,d_{iN}]}^T$ of dictionary $\mathit{D}$ is obtained. The input excitation matrix $\mathit{X}$ can be expressed as

$$\mathit{X}=\left[\begin{array}{cccc}{x}_1^0,& x_1^1,& \cdots ,& x_1^{n-1}\\ x_2^0,& x_2^1,& \cdots ,& x_2^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ x_N^0,& x_N^1,& \cdots ,& x_N^{n-1}\end{array}\right].$$

(22)

According to the neural network weight determination method, $w={\left({\mathit{X}}^T\mathit{X}\right)}^{-1}{\mathit{X}}^T{\mathit{d}}_i$. Therefore, the predicted atomic value after filtering can be obtained as ${\mathit{d}}_i^{*}=\mathit{X}\mathit{w}$, and the filtered dictionary can be expressed as ${\mathit{D}}^{*}=[{\mathit{d}}_1^{*},{\mathit{d}}_2^{*},\cdots ,{\mathit{d}}_i^{*}\cdots ,{\mathit{d}}_L^{*}]$.

As shown in Figure 6, it can be seen from the figure that the real quintic polynomial phase signal with additive Gaussian white noise is a smooth curve after learning by the direct weight determination method. The dictionary atom after RLS-DLA learning obviously contains additive white Gaussian noise. Therefore, the dictionary composed of atoms learned by the direct weight determination method is denoised by sparse representation, and the denoising effect is very obvious.

3.3. Sparse Blind the TOA Offset Estimation

In this section, sparse dictionary learning is used to represent signals sparsely. Meanwhile, in order to improve the resolution of TOA, the signal aug-expanding matrix is added and its sparsity and coherence are proven.

3.3.1. Construction of Sparse Basis Matrix

Sparsity is the premise of the compressed sensing algorithm and the necessary condition for signals. However, most signals do not meet the conditions of sparsity or compressibility. In this case, signals need to be transformed into sparse signals in a certain transform domain, and the transformed signals become compressible signals. The signal $\Re \left(\mathit{F}\right)={[{\varsigma }_1,{\varsigma }_2,\cdots ,{\varsigma }_l,\cdots ,{\varsigma }_L]}^T\in \mathbb{R}$ sampled by us was a discrete signal in the frequency domain with dimension L. The dictionary $\Psi$ estimated by the direct weight determination method was used to estimate the sparse representation coefficient $\alpha$ as the sparse matrix. The signal $\Re \left(\mathit{F}\right)$ can be expressed as

$$\Re \left(\mathit{F}\right)=\Psi \mathit{\alpha },$$

(23)

where $\Psi \in {\mathbb{R}}^{L\times L}$, and $L=256$ is the number of dimensions of $\mathit{F}$ of the sampled signal. $\mathit{\alpha }={[{\alpha }_1,{\alpha }_1,\cdots ,{\alpha }_l,\cdots ,{\alpha }_L]}^T$ is the coefficient vector, which is another representation of $\mathit{F}$ and is sparse on ${\Psi }^{L\times L}$.

3.3.2. Construction of Measurement Matrix

In the whole process of compressed sensing, the design of the measurement matrix is a key step. The properties of the measurement matrix are related to whether the compression can be achieved and whether the signal can be reconstructed accurately. How can an appropriate observation matrix be designed, which can not only achieve the purpose of compressed sampling but also ensure that the signal can be reconstructed without distortion? The restricted isometry property (RIP) can solve this problem in theory [39]. The random Gaussian measurement matrix is the most widely used in compressed sensing. We constructed a matrix $\Phi$ with the size of $K*L$, and each element in the matrix independently followed a Gaussian distribution with mean value $\upsilon =0$ and variance of $\sigma =\sqrt{1/K}$, which can be expressed as follows

$${\Phi }_{k,l}\sim N(\upsilon ,\sigma ),l\in [1,L],k\in [1,K],$$

(24)

where L is the dimensions of the sampled data, and N is the normal distribution function. K is the number of CIR dimensions to be restored. The random Gaussian measurement matrix has a strong randomness. It can be proven that, when the measurement number $K\ge cMlog(L/M)$ of the random Gaussian measurement matrix, RIP conditions will be met with a great probability. In general, there is $M<<K\le L$, and c is a tiny constant. In the compressed sensing process, a random Gaussian measurement matrix is widely used mainly because it is not related to most orthogonal bases or orthogonal dictionaries, and the number of measurements required for accurate reconstruction is relatively small.

3.3.3. CIR Resolution Augmentation Matrix

In order to improve the time domain resolution and recover CIR to estimate a more accurate TOA, we constructed a matrix to extend the arrival delay $\tau$ [40], which can be expressed as

$${\mathbf{S}}_{nk}=sinc(\pi \left({\mathbf{A}}_{nk}\right)/T),n\in [1,N],k\in [1,K],$$

(25)

where $sinc$ is a filter constructed by the $sinc$ function, often used in anti-aliasing techniques. T is the sampling period. N is the resolution of CIR of data to be restored, and N is the control parameter of the resolution of data to be restored, which can be used to improve the resolution of data to be restored, and $N>K$. ${\mathbf{A}}^{N\times K}\in \mathbb{R}$ can be represented as

$$\begin{array}{c}\hfill \mathbf{A}=\left[\begin{array}{ccccc}& 0,& {\nu }_{1,1},& \cdots ,& {\nu }_{1,k-1}\\ & -{\nu }_{2,2},& -{\nu }_{2,1},& \cdots ,& {\nu }_{2,k-3}\\ & -{\nu }_{3,4},& -{\nu }_{3,3},& \cdots ,& {\nu }_{3,k-5}\\ & \ddots & ⋮\\ & -{\nu }_{n,2*(n-1)},& {\nu }_{n,2*(n-1)-1},& \cdots ,& {\nu }_{n,k-2*(n-1)-1}\end{array}\right],\end{array}$$

(26)

where $\nu$ is the sampling interval and ${\mathbf{S}}^{N\times K}\in \mathbb{R}$.

3.3.4. Sparse Recovery

The $K*L$ dimension matrix formed by the product of measurement matrix ${\Phi }^{K\times L}$ and transformation basis ${\Psi }^{L\times L}$ is called the perception matrix, and observes the signal and obtains the observation vector $\mathbf{y}={[y_1,y_2,\cdots ,y_n,\cdots ,y_N]}^T\in \mathbb{R}$. The observed signal can be expressed as

$$\mathbf{y}=\Phi \Re \left(\mathit{F}\right)=\Phi \Psi \mathit{\alpha }.$$

(27)

Meanwhile, in order to improve the resolution, we constructed a matrix ${\mathbf{S}}^{N\times K}$ to extend the arrival delay, so (27) can be expressed as

$$\mathbf{y}=\mathbf{S}\Phi \Re \left(\mathit{F}\right)=\mathbf{S}\Phi \Psi \mathit{\alpha }.$$

(28)

In the context of CS, the CIR reconstruction problem is now formulated as

$$\begin{array}{cc}\hfill & \underset{\alpha }{min}{\parallel \mathit{\alpha }\parallel }_0,\hfill \\ \hfill & {s.t.\parallel \mathbf{y}-\mathbf{S}\Phi \Psi \mathit{\alpha }\parallel }_2^2\le \epsilon ,\hfill \end{array}$$

(29)

where, for some small $\epsilon$, ${\parallel ·\parallel }_0$ denotes $l_0$-norm, which counts the number of non-zero elements, and ${\parallel ·\parallel }_2$ stands for $l_2$-norm. Note that $\epsilon$ approximates ${\sigma }^2$ if we have enough data. However, $l_0$-norm minimization is a complex combinatorial problem [41]. In order to improve the computational efficiency, it is usually extended to $l_1$-norm minimization, because it can be solved by convex optimization [42,43]. Compared with the inverse Fourier transform, this CS-based approach enhances sparsity and effectively minimizes the pseudo-sidelobe.

3.3.5. CS with Complex-Valued Targets

The CS estimates a sparse signal vector in a probabilistic manner. It maximizes the posterior probability of the sparse signal vector at a given observed value. In the posterior probability maximization, the CS theory is applicable to the processing of real data, while CSI data are complex. In order to apply this theory to TOA estimation, observation data should be realizable in order to build an optimization model of a real number field [44]. For the application of CS in complex valued CIR estimation, Equation (28) is first extended to the equivalent real-valued form by

$$\begin{array}{cc}\hfill \tilde{H}& =\left[\begin{array}{c}\Re (\Phi \Psi \mathbf{S})-\Im (\Phi \Psi \mathbf{S})\\ \Im (\Phi \Psi \mathbf{S})\Re (\Phi \Psi \mathbf{S})\end{array}\right]*\left[\begin{array}{c}\Re \left(\mathit{\alpha }\right)\\ \Im \left(\mathit{\alpha }\right)\end{array}\right]\hfill \\ \hfill & =\tilde{\mathbf{S}}\tilde{\Phi }\tilde{\Psi }\tilde{\mathit{\alpha }},\hfill \end{array}$$

(30)

where $\tilde{H}\in {\mathbb{R}}^{2N\times 1}$, $\tilde{\Phi }\in {\mathbb{R}}^{2K\times 2L}$, $\tilde{\Psi }\in {\mathbb{R}}^{2L\times 2L}$, $\tilde{\mathbf{S}}\in {\mathbb{R}}^{2N\times 2K}$, and $\tilde{\mathit{\alpha }}\in {\mathbb{R}}^{2L\times 1}$. This expansion increases the dimension of the observation matrix and separates the reconstruction of real and imaginary values of a channel magnitude in each delay tap.

3.3.6. Sparse and Coherence Proof

In the whole process of compressed sensing, the design of a measurement matrix is a key step. The properties of the measurement matrix are related to whether the compression can be achieved and whether the signal can be reconstructed accurately. For a given measured value y, the solution of Formula (28) is an underdetermined problem. Candes et al. proposed a solution: as long as the measurement matrix conforms to the RIP, a definite solution can be obtained; that is, the existence, accuracy, and uniqueness of the original signal reconstruction can be guaranteed. RIP is defined as:

$$\begin{array}{c}\hfill {(1-\delta )\parallel x\parallel }_2^2\le {\parallel x\parallel }_2^2\le (1+\delta ){\parallel x\parallel }_2^2,\end{array}$$

(31)

where $\delta$ is a constant and $\delta \in (0,1)$. ${\parallel ·\parallel }_2$ stands for ${\ell }_2$-norm. In this paper, the measurement matrix is a random Gaussian measurement matrix, and there is a high probability that RIP conditions are met.

According to Baraniuk, the equivalent condition of RIP is that the observation matrix $\Phi$ and the sparse representation basis $\Psi$ are irrelevant. However, since we added the $\mathbf{S}$ matrix to improve the CIR resolution, the equivalent condition of RIP is therefore irrelevant between $\Phi$, $\Psi$ and $\mathbf{S}$. In practical application, we can use an equivalent condition of criterion, namely incoherence, to guide the design of the measurement matrix. Incoherency means that the row vector ${\Phi }_i$ of the matrix observation system $\Phi$ cannot be sparsely represented by the column vector ${\Psi }_i$ in the basis $\Psi$, nor can it be represented by the column vectors of the ${\mathbf{S}}_i$ matrix. It can be measured by coherence and is expressed as

$$\begin{array}{c}\hfill \mu (\Phi ,\Psi ,\mathbf{S})=\underset{1\le k\le K,1\le l\le L,1\le n\le N}{max}\mid \langle {\Phi }_k,{\Psi }_l,{\mathbf{S}}_n\rangle \mid .\end{array}$$

(32)

In this paper, ${\Phi }^{K\times L}$ is a sensing matrix of random normal distribution with mean $\upsilon =0$ and variance $\sigma =\sqrt{1/K}$. According to the normal distribution function, the ratio within $3\sqrt{1/K}$ is approximately 99.73%, so

$$0\le |{\Phi }_{kl}^{K\times L}{|}_1\lesssim 3\sqrt{1/K},k\in [1,K],l\in [1,L].$$

(33)

The sparse matrix $\Psi$ estimated by dictionary learning is derived from collected CSI data. The indoor location signal model of the wiFI-WASP system used to collect CSI data is suitable for the Saleh–Valenzuela propagation model. The Saleh–Valenzuela propagation model is an indoor radio channel statistical multipath model that fits well with our measurements, through which, the received signal rays arrive in clusters. These rays have an independent uniform phase, as well as independent Rayleigh amplitudes whose variances decay exponentially with the cluster and ray delay. Light in clusters and within clusters forms Poisson arrival processes with different but fixed rates. According to the model, the received signal amplitude follows a normal distribution, which can be expressed as

$$a_k=N(0,{\sigma }_k^2)+jN(0,{\sigma }_k^2),$$

(34)

where k indicates the number of paths and $N(0,{\sigma }_k^2)$ shows a value that is in accordance with the normal distribution. Equation (34) creates the Rayleigh fading. The variance ${\sigma }_k^2$ is the average power of the k path, so the strength of the paths within the clusters is given as [45,46]:

$${\sigma }^2\left(k\right)=e^{-T_l/\Gamma }{e}^{{-\tau }_{il}/\gamma },$$

(35)

where $T_l$ is the arrival time of the lth cluster and ${\tau }_{il}$ is the arrival time of the ith path in the lth cluster. We set the time constants $T_l$ and ${\tau }_{il}$ or the inter and intra cluster as 300 ns and 5 ns, respectively, as Poisson distributions. $\Gamma$ is a constant of cluster arrival decay time, $\gamma$ indicates a constant of ray arrival decay time, $\Gamma$ = 60 ns, and $\gamma$ = 20 ns. According to the normal distribution function, the ratio within $3\sqrt{1/{\sigma }^2\left(k\right)}$ is approximately 99.73%. Thus, it can be seen that the range of $l_1$ norm of any element in $\Psi$ can be expressed as

$$0\le |{\Psi }_{uv}^{L\times L}{|}_1\lesssim 3\sigma .$$

(36)

The $sinc$ function of CIR resolution enhancement matrix ${\mathbf{S}}^{N\times K}$ can be expressed as

$${\mathbf{S}}_{nk}=sinc\frac{\pi {\mathbf{A}}_{nk}}{T}=\frac{sin(\left(\pi {\mathbf{A}}_{nk}\right)/T)}{\pi {\mathbf{A}}_{nk}/T},n\in [1,N],k\in [1,K].$$

(37)

Obviously, the value range of any element in matrix $\mathbf{S}$ is

$$0\le |{\mathbf{S}}_{nk}{|}_1\le 1.$$

(38)

Therefore, the coherence of $\Phi$, $\Psi$, and $\mathbf{S}$ when multiplied by any rows and columns is in the range of

$$\begin{array}{c}\hfill 0<\mu (\Phi ,\Psi ,\mathbf{S})\lesssim \frac{9L\sigma }{\sqrt{K}},\end{array}$$

(39)

where the coherence of $\mu$ is closely related to the number of rows and columns of $\Phi$ and $\Psi$, but not to the number of rows and columns of $\mathbf{S}$ and $\sigma \ll 1$. Since the elements in the $\mathbf{S}$ matrix are generally less than 1, it is still helpful for incoherence $\mu$.

4. Experimental Results

The performance of the proposed algorithm was evaluated under outdoor line-of-sight (LOS) conditions and indoor conditions, respectively. In each experiment, an 802.11ac wireless LAN was built, and six sniffers were deployed, running in 149 channels with a bandwidth of 80 MHz. One target device was used for location verification and one PC was used to collect data from the sniffer and target device online to estimate the location of the target device.

In practical applications with LOS propagation, usually the first peak of the impulse response is the TOA of the direct path. However, in a non-LOS or multipath environment, the first-arrival path usually attenuates significantly and may even have a lower signal intensity than the reflection path. In addition, multiple reflected signals that overlap the direct path signal distort the first peak of the combined signal [10]. Existing ToA estimation algorithms are based on thresholding the first peak, using a prior model [3], or the MUSIC/ESPRIT algorithm. The performance of the proposed approach was compared against MUSIC and ESPRIT.

4.1. Outdoor Test

The system was tested under outdoor LOS conditions to evaluate its performance. We deployed six sniffers around a WiFi network for the outdoor LOS tests with one AP and one WiFi dongle. The topology of the system is shown in Figure 7. The size of the experimental area was 900 m${}^2$. The target moved across the test area. Figure 5 shows the positioning results for the outdoor test. It can be seen that the positioning results are consistent with the actual locations.

Figure 8 shows the positioning accuracy of our system when the TOA is estimated using the CS and MUSIC and ESPRIT algorithm. It can be seen that the 72-percentile positioning error is 0.5 m, and 1 m in 90% with the CS. Compared with MUSIC and ESPRIT’s TOA estimation methods, the accuracy of CS’s TOA estimation method is much higher.

4.2. Indoor Test

Localization accuracy depends on the multipath environment, the materials used in the walls, the presence of metal objects, the density of sniffer deployment, and many other factors, so indoor localization testing is an important step in validating wireless positioning methods. The assembly hall covers an area of around 70 m${}^2$. We deployed six sniffers and a standard Wi-Fi device in the conference hall and placed them at known locations to collect data from 17 locations of target devices inside the conference hall. Figure 9 shows the real and estimated location of the target. It can be seen that, except for the poor reflection effect of the seat in the Y-axis direction, the estimated location is consistent with the actual location, which can meet the needs of many indoor positions.

Figure 10 shows the positioning effect when using CS, MUSIC, and ESPRIT algorithms to estimate the TOA. We observed that CS achieves a localization error of approximately the 48th percentile at 0.5 m and 1 m, respectively. Compared with MUSIC and ESPRIT’s TOA estimation method, CS’s TOA estimation method has a higher accuracy. The SpotFi positioning system is an indoor localization method proposed by m. Kotaru et al. The SpotFi localization system has three antennas and is deployed in the indoor office with 5–6 AP. The localization error of approximately 18% is 0.5 m and approximately 60% is 1 m [25]. The ArrayTrack localization system is an indoor localization method proposed by Jie Xiong et al. In the same experimental environment as SpotFi, the localization error of approximately 30 % is 1 m.

5. Conclusions

With the increasing demand for positioning, it is difficult for existing positioning methods to meet the convenience of deployment equipment without satellite positioning. The proposed TOA estimation method is expected to better adapt to the business requirements of equipment deployment. In this paper, a method based on a neural network to determine the direct weight of dictionary filtering was designed, and the filtering effect is proven by quintic polynomial. The learned dictionary was used as the sparse matrix of compressed sensing, the augmentation and expansion matrix was added to increase the CIR resolution after sparse recovery, and the sparsity of the dictionary, measurement matrix, augmentation, and expansion matrix was proven. The proposed algorithm was validated by experiments on the indoor and outdoor positioning system. The results show that the proposed algorithm has a high positioning accuracy. Compared with the fingerprint positioning method, the positioning method adopted in this paper has the advantage of not needing to sample the data in the positioning environment in advance, so it has a broad application prospect.