A Higher-Order Singular Value Decomposition-Based Target Localization Algorithm for WiFi Array Systems

Abstract

Traditional Angle of Arrival (AoA)-based WiFi array indoor localization algorithms do not fuse Channel State Information (CSI) inter-packet data for estimation, which makes WiFi arrays less effective for localization in complex indoor environments. Most algorithms are overburdened leading to inefficient localization. To address these issues, in this article, an indoor positioning algorithm based on Higher-Order Singular Value Decomposition (HOSVD) is proposed. First, the CSI data are reconstructed as a new measurement matrix by borrowing subcarriers, and a third-order tensor is constructed. Next, tensor compression techniques are used to reduce computational complexity and the signal subspace is obtained by HOSVD. Then, the AoA is obtained by the Reduced Dimension Multiple Signal Classification (RD-MUSIC) method. Finally, the coordinates of the target can be obtained by triangulating the AoAs of the three Access Points (APs). According to the simulation experiments, the AoA can be estimated accurately at a low SNR and with low snapshots. In practical experiments, we can successfully estimate the AoA in complex indoor environments with shorter timelines using HOSVD without modifications to commercial hardware and produce a lower AoA error and localization error rates compared to other algorithms. The effectiveness of our proposed algorithm is proven by simulations and practical experiments.

1. Introduction

The demand for accurate indoor positioning systems has witnessed a significant rise in recent years. While traditional Global Navigation Satellite System (GNSS) technology has revolutionized outdoor navigation, it falls short when it comes to providing precise positioning information within indoor environments. This limitation has spurred the development of WiFi indoor positioning, a technology that utilizes WiFi signals to enable accurate localization and navigation within buildings.

In recent years, with the mass penetration of wireless networks, there have been more opportunities for the much-needed location-based services to be realized [1]. Internet of Things (IoT) technology [2,3] is continuously improving which brings stronger drivers for the development of indoor positioning. In indoor scenarios, CSI defined by the IEEE 802.11 WiFi protocol is a better reflection of CSI [4] than received signal strength information. Therefore, CSI data are beginning to be used by academics for WiFi indoor localization.

The AoA, AoD, and ToF required for indoor positioning can be estimated by traditional algorithms [5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]. The MUSIC algorithm [2] relies on a computationally intensive grid search for parameter estimation. High-accuracy estimates take a long time to obtain. For example, SpotFi [6], a median positioning accuracy of 40cm can be achieved in an office environment using the 2D MUSIC algorithm. However, this method has some disadvantages: the 2D MUSIC algorithm takes a lot of computational time. The estimation performance deteriorates in low snapshot situations and at a low SNR [20,21,22]. Soltanaghaei [8] combines the TDoA estimation and the AoA estimation of multipath CSI signals. A single AP localization is achieved. It relies on a sufficient number of multipaths while ignoring multiple reflected paths. This makes this system ineffective in complex indoor environments. The M3 [11] system was proposed, which uses the super-resolution algorithm SAGE to jointly estimate the AOA, AOD, and relative ToF. The algorithm succeeded in achieving a centimeter-level localization accuracy. However, incorrect multipath numbers can lead to large errors in the parameter estimation results at a low SNR. RoArray [13] transforms the AoA estimation into a sparse recovery problem. This method is able to produce sharp and sparse AoA spectra. Direct path identification is achieved by jointly estimating the ToF and AoA for all paths. However, the long computation time makes it difficult to meet the requirements of real-time localization. The above method has a poor estimation performance and high computational burden at a low SNR.

Some algorithm implementations require modifications to the hardware settings. This makes WiFi array localization difficult to implement on commercial NICs. Chronos [14] is a single AP positioning system. It achieves a decimeter-level location accuracy by combining the ToA and AoA estimations. However, it requires modifications to existing WiFi protocols and increases the bandwidth of the CSI signal through frequency hopping. This means that it is difficult to implement on actual commercial APs. Yang [15] uses a two-dimensional matrix pencil (2D-MMP) approach to achieve a fast estimation of the AoA and ToF with a decimeter-level accuracy. However, the algorithm must meet at least four receiving antennas, which requires modification of the receiving equipment. This is not applicable in the current situation where most commercial network interface cards (NICs) have only three receive antennas. Moreover, the estimation performance is poor at a low SNR.

The above articles all use matrix analysis. Because of the multi-dimensional structural nature of multi-packet CSI data, scholars have started to use tensor analysis to extract multi-dimensional structural information to improve algorithm performance [23]. Zhao et al. [24] proposed a tensor-based localization algorithm that makes full use of the multipath component of the CSI to achieve decimeter-level localization by fusing the arrival time difference and AoA information of each MPC. Gong [25] formulates the location parameter estimation problem as a Tensor Decomposition problem. The received signal is decomposed into three steering vectors, the AoA, AoD, and ToF, by a pre-processing step. Only one AP is used for localization. Both of the above algorithms use one steering vector of the Tensor Decomposition to estimate a parameter. HOSVD [26,27,28,29,30] can make full use of the multi-dimensional structure of the CSI tensor. The HOSVD has higher estimation accuracy at a low SNR with low snapshots; therefore, HOSVD is used to process the tensor.

In this paper, an algorithm for a WiFi array localization system based on HOSVD is presented. This algorithm wants to solve two problems. The first one allows the algorithm to have a better performance for estimation in complex environments. The second one is to solve the problem of estimation efficiency of the algorithm so that the algorithm can have a shorter estimation time. First, CSI data are reconstructed as a new measurement matrix by borrowing subcarriers. The multi-packet data are stored in a tensor to capture the multi-dimensional structure. Then, a tensor that is as small as possible is obtained by tensor compression techniques and preserves most of the valid information of the signal. Next, the signal subspace can be obtained by HOSVD, and the AoA can be obtained by means of the RD-MUSIC method. Finally, the coordinates of the target can be obtained by triangulating the AoAs of the three AP points. We list the main contributions of the proposed method as follows:

(1)

The multi-dimensional structure of the multi-packet CSI data is stored in the tensor, which improves the accuracy of the target estimation. The method implements AoA estimation for WiFi indoor positioning on a tensor domain.

(2)

Reduced computational complexity. To store the same amount of multi-packet CSI data, we constructed a tensor with a smaller dimension than the ones constructed by other algorithms.

(3)

Further reduction in computational complexity. As the tensor stores a large amount of CSI data, tensor compression techniques are used to reduce computational complexity and retain most of the valid information of the signal. With the Lagrange multiplier method, the AoA and ToF are decoupled, which further reduces the dimension of the peak search. Both methods ensure accuracy and reduce the complexity of the calculation.

The notations used in this paper are defined in

Table 1. Notations.

Notations	Definitions
${(·)}^T$	transpose operator
${(·)}^H$	conjugate transpose operator
${(·)}^{†}$	pseudo-inverse
${(·)}^{*}$	conjugate
⊗	Kronecker product operator
⊙	Khatri–Rao product operator
${\mathit{I}}_K$	$K\times K$ identity matrix
${\mathbb{Q}}^{M\times N}$	$M\times N$ dimensional rational numbers matrix set
${\mathbb{C}}^{M\times N}$	$M\times N$ dimensional complex matrix set

.

2. Data Model

The initial matrix of the collected CSI data is first presented in order to better introduce the data processing of the proposed algorithm. Consider an AP with M identical antenna Uniform Linear Array (ULA) arrays spaced at half of a wavelength. The transmitter uses a single antenna to transmit an orthogonal frequency division multiplexing (OFDM) signal with N subcarriers. The n-th subcarrier frequency is $f_N$, and the continuous subcarrier spacing is $\Delta f$. There are L channels in the room, and I CSI packets are collected at a time. The i-th CSI packet can be expressed as

$${\mathit{H}}^i=\left[\begin{array}{ccc}{h}_{0,0}^i& \cdots & h_{0,N-1}^i\\ ⋮& \ddots & ⋮\\ h_{M-1,0}^i& \cdots & h_{M-1,N-1}^i\end{array}\right],$$

(1)

where $h_{M-1,N-1}^i$ is the data received at the M-th antenna and the N-th subcarrier. The measured value of CSI at the AP point at the m-th antenna in the i-th packet of the n-th subcarrier can be expressed as

$$h_{m,n}^i=\sum _{l=1}^L{\alpha }_l^i\times exp\left(-j\frac{2\pi f}{c}mdsin{\theta }_l^i\right)\times exp\left(-j2\pi n\Delta f{\tau }_l^i\right)+w_{m,n}^i,$$

(2)

where f and d are the central frequency and the spacing between consecutive antennas; c and ${\alpha }_l^i$ are the speed of light and the complex fading of the l-th path, respectively. ${\theta }_l^i$ and ${\tau }_l^i$ are the AoA and ToF of the l-th path, and $w_{m,n}^i$ is the Gaussian noise. It was observed through the composition of Equation (2) that there are parameters of interest to us within the CSI matrix, such as the AoA and ToF.

The phase shift AoA function for the l-th propagation path of the i-th packet can be expressed as:

$$\Phi \left({\theta }_l^i\right)=exp\left(-j\frac{2\pi f}{c}mdsin{\theta }_l^i\right),$$

(3)

and the AoA steering vector can be expressed as

$${\mathit{a}}_l^i={\left[\begin{array}{cccc}0\hfill & \Phi \left({\theta }_l^i\right)\hfill & \cdots \hfill & {\Phi }^{M-1}\left({\theta }_l^i\right)\hfill \end{array}\right]}^T$$

(4)

Thus, the array manifold pattern for the AoA of the i-th packet can be expressed as

$${\mathit{A}}^i={\left[\begin{array}{cccc}{\mathit{a}}_1^i\hfill & {\mathit{a}}_2^i\hfill & \cdots \hfill & {\mathit{a}}_L^i\hfill \end{array}\right]}^T$$

(5)

Similarly, the complex exponent of the phase shift introduced between adjacent subcarriers can be expressed as

$$\Omega \left({\tau }_l^i\right)=exp\left(-j2\pi n\Delta f{\tau }_l^i\right)$$

(6)

The steering vector for the ToF of the i-th packet can be expressed as

$${\mathit{b}}_l^i={\left[\begin{array}{cccc}0\hfill & \Omega \left({\tau }_i^i\right)\hfill & \cdots \hfill & {\Omega }^{N-1}\left({\tau }_l^i\right)\hfill \end{array}\right]}^T$$

(7)

In the subsequent process, Equations (4) and (7) will be mathematically computed to solve the problem of receiving antennas resources. The array manifold pattern for the ToF of the i-th packet can be expressed as

$${\mathit{B}}^i={\left[\begin{array}{cccc}{\mathit{b}}_1^i\hfill & {\mathit{b}}_2^i\hfill & \cdots \hfill & {\mathit{b}}_L^i\hfill \end{array}\right]}^T$$

(8)

Then, the signal $\mathit{X}$ of the i-th packet can be expressed as

$${\mathit{X}}^i=({\mathit{A}}^i\odot {\mathit{B}}^i){\mathit{S}}^i+{\mathit{N}}^i,$$

(9)

where ${\mathit{S}}^i$ and ${\mathit{N}}^i$ are the signal matrix and noise matrix of the i-th packet, respectively.

3. Proposed Methods

3.1. Steering Vector Reconstruction

The algorithm proposed in this paper consists of 7 parts, including direction vector reconstruction, tensor construction, tensor compression, estimating the number of sources, obtaining the subspace through HOSVD, RD-MUSIC, and triangular positioning. First of all, the problem of the limited resources of receiving antennas should be solved. The number of multipath signals propagated indoors is generally from 4 to 5; a commercial network card with just 3 antennas is unable to estimate a number of multipath signals that is greater than 3. The signal subspace cannot be estimated from this tensor $\mathcal{X}$ by HOSVD. Therefore, we can use the resource-rich frequency resources to expand the spatial resources because $N\gg M$. We use $n_{tc}$ subcarriers to expand the number of receiving antennas. The expansion process is shown in Figure 1.

As shown in Figure 1, the decomposition of $\mathcal{X}$ can be divided into three steering vectors $\mathit{a}$, $\mathit{b}$, and $\mathit{c}$, with dimensions M, N, and I. However, when M<L, the target’s AoA cannot be estimated. The value of N is much larger than the value of L. Therefore, consider that the number of dimensions in $\mathit{a}$ can be borrowed from the number of dimensions in $\mathit{b}$. The dimension of $\mathcal{X}$ will be transformed into $(M\times n_{tc})\times (N-n_{tc})\times I$. The dimension in $\tilde{\mathit{a}}$ can be larger than the value of L. The mathematical expression for the expansion process of the direction $\tilde{\mathit{a}}$ can be expressed as

$${\tilde{\mathit{a}}}_l^i={\left[\begin{array}{cccc}0\hfill & \Phi \left({\theta }_l^i\right)\hfill & \cdots \hfill & {\Phi }^{M-1}\left({\theta }_l^i\right)\hfill \end{array}\right]}^T\otimes {\left[\begin{array}{cccc}0\hfill & \Omega \left({\tau }_l^i\right)\hfill & \cdots \hfill & {\Omega }^{n_{tc}-1}\left({\tau }_l^i\right)\hfill \end{array}\right]}^T$$

(10)

In this case, there is an $M\times n_{tc}$ array of elements in the $\tilde{\mathit{a}}$ dimension. The ToF steering vector after borrowing the $n_{tc}$ subcarrier can be expressed as

$${\tilde{\mathit{b}}}_l^i={\left[\begin{array}{cccc}0\hfill & \Omega \left({\tau }_l^i\right)\hfill & \cdots \hfill & {\Omega }^{N-1-n_{tc}}\left({\tau }_l^i\right)\hfill \end{array}\right]}^T$$

(11)

There are $N-n_{tc}$ array elements in the $\tilde{\mathit{b}}$ dimension. The problem of receiving antenna resource constraints has been solved by Equations (10) and (11). Then, the steering vector of the i-th packet signal can be expressed as

$${\tilde{\mathit{c}}}_l^i={\tilde{\mathit{a}}}_l^i\odot {\tilde{\mathit{b}}}_l^i$$

(12)

The array manifold of the entire signal of the i-th packet can be expressed as

$${\tilde{\mathit{C}}}^i={\left[\begin{array}{ccc}{\tilde{\mathit{c}}}_1^i\hfill & \cdots \hfill & {\tilde{\mathit{c}}}_L^i\hfill \end{array}\right]}^T$$

(13)

With the above calculation, the received signal of the i-th packet can be converted by Equation (9):

$$\begin{array}{cc}\hfill {\mathit{X}}^i& ={\tilde{\mathit{C}}}^i\times {\mathit{S}}^i+{\mathit{N}}^i={\left(\tilde{\mathit{a}}\right.}^{\mathit{i}}\odot {\tilde{\mathit{b}}}^{\mathit{i}}{\mathit{S}}^{\mathit{i}}+{\mathit{N}}^{\mathit{i}},\hfill \end{array}$$

(14)

where ${\mathit{X}}^i=\left[\begin{array}{ccc}{\mathit{x}}_1^i\hfill & \cdots \hfill & {\mathit{x}}_L^i\hfill \end{array}\right],\mathit{S}={\left[\begin{array}{ccc}{\mathit{s}}_1^i\hfill & \cdots \hfill & {\mathit{s}}_L^i\hfill \end{array}\right]}^T\in {\mathbb{C}}^{L\times I},{\tilde{\mathit{a}}}^i={\left[\begin{array}{ccc}{\tilde{\mathit{a}}}_1^i\hfill & \cdots \hfill & {\tilde{\mathit{a}}}_L^i\hfill \end{array}\right]}^T,{\tilde{\mathit{b}}}^i={\left[\begin{array}{ccc}{\tilde{\mathit{b}}}_1^i\hfill & \cdots \hfill & {\tilde{\mathit{b}}}_L^i\hfill \end{array}\right]}^T,{\mathit{N}}^i=\left[\begin{array}{ccc}{\mathit{n}}_1^i\hfill & \cdots \hfill & {\mathit{n}}_L^i\hfill \end{array}\right]$, and ${\mathit{N}}^i$ is Gaussian white noise.

3.2. Tensor Construction

The second part describes how to perform tensor construction. Firstly, the steering vector reconstruction is first realized by spatial smoothing, and the measurement matrix is obtained. The tensor is constructed by stacking the measurement matrices obtained from different packets, and we find that the tensor constructed in this paper can store the same amount of information with smaller data dimensions. When collecting a set of data packets, we converted them into the tensor we needed using the following method, as shown below.

As shown in Figure 2, when $n_{tc}$ is taken as 3, a part of each row of the computational matrix of a single packet is taken as a subarray for smoothing. Then, the smoothed-out matrix is placed in the same column vector in order. The different column vectors obtained from all packets are put in the same matrix. We can obtain the matrix $\mathit{X}$ of all packet fusion, and then it can be transformed into the required tensor $\mathcal{X}$. The dimension of the tensor $\mathcal{X}$ is 9 × 27 × I. The dimension of the tensor of Tensor Decomposition [17] is 30 × 32 × I. This means that we use a smaller dimensional tensor to store the same amount of CSI data.

3.3. Tensor Compression

The third part is to introduce how to use tensor compression to reduce the computation time. Data processing for a tensor with large dimensions takes a lot of time. Therefore, in order to reduce the complexity of the algorithm, the tensor obtained in the previous part will be compressed into a smaller tensor. Before compressing the tensor, a compression matrix that retains most of the information of the CSI data needs to be obtained. The obtained tensor has a K ($K<I$) rank. The tensor compression process is shown in Figure 3, where $\mathit{W}$ is the first K left singular value vector after the singular value decomposition of the mode-3 expansion of the tensor $\mathcal{X}$.

Firstly, the mode-3 unfolding of the tensor $\mathcal{X}$ can be expressed as

$$\left[{\mathcal{X}}_{\left(3\right)}\right]=\left[\begin{array}{c}{\left[{\mathcal{X}}_{\left(3\right)}\right]}_1\\ {\left[{\mathcal{X}}_{\left(3\right)}\right]}_2\\ ⋮\\ {\left[{\mathcal{X}}_{\left(3\right)}\right]}_I\end{array}\right]=\left(\tilde{\mathit{a}}\odot \tilde{\mathit{b}}\right)S+N_{\left[{\mathcal{X}}_{\left(3\right)}\right]}$$

(15)

Then, we can obtain the compression matrix $\mathit{W}$ by the Tucker 3 decomposition. The tensor after compression can be expressed as

$$\left[{\tilde{\mathcal{X}}}_{\left(3\right)}\right]=\left[{\mathcal{X}}_{\left(3\right)}\right]\ast \mathit{W}=[\tilde{\mathit{a}}\odot \tilde{\mathit{b}}]S\ast \mathit{W}+{\mathit{N}}_{\left[{\mathit{\chi }}_{\left(3\right)}\right]}\ast \mathit{W}=[\tilde{\mathit{a}}\odot \tilde{\mathit{b}}]S^{\prime }+{\mathit{N}}_{\left[{\mathcal{X}}_3\right]}^{\prime }$$

(16)

where $\mathit{W}\in {\mathbb{C}}^{I\ast K}$. Before the use of tensor compression techniques, the signal subspace was obtained with 3 high-complexity mode-n decompositions. After using tensor compression techniques, the signal subspace is obtained with only 1 high-complexity mode-n decomposition and 3 low-complexity mode-n decompositions.

3.4. Estimation of the Number of Sources

The fourth part is the determination of the number of sources. The determination of the exact number of sources is required to accurately estimate the AoA of the target. To improve the accuracy of estimating the number of targets, we use the Minimum Description Length (MDL) [31] to determine the number of targets. The MDL Equation (17) determines the number of incoherent signals required for estimation in an indoor environment. Assuming that the Covariance matrix $X{X}^H$ of Equation (14) is of the order G and the number of incoherent signals is L, L can be determined by minimizing the estimator as follows

$${\widehat{L}}_{MDL}=\underset{L\in \{0,\cdots ,G-1\}}{argmin}MDL\left(L\right)$$

(17)

where

$$\begin{array}{cc}\hfill MDL\left(L\right)& =(G+M-1)\times (G-L)log\left(\frac{\frac{1}{G-L}\times {\sum }_{g=1}^{G-L}\mathit{\lambda }}{{\left({\prod }_{g=1}^{G-L}\mathit{\lambda }\right)}^{\frac{1}{G-L}}}\right)\hfill \\ \hfill & +\frac{L}{2}\times (2G-L+1)\times log(G+M-1),\hfill \end{array}$$

(18)

where $p\in \left\{0,1,\cdots G-1\right\}$, and the integer g that makes the MDL(p) the smallest is the number of estimated incoherent signals. $\mathit{\lambda }$ is the diagonal vector of the eigenvalue matrix of the $X{X}^H$ after EVD, which can be epressed as: $\lambda ={\left[{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_g,\cdots ,{\lambda }_G\right]}^T,{\lambda }_1\ge {\lambda }_2\ge \dots \ge {\lambda }_g\ge \dots \ge {\lambda }_G$.

3.5. Tensor Domain Subspace Estimation Methods

The fifth part is a description of how to find the signal subspace of the tensor using HOSVD. First, we start by introducing the composition of the tensor. The $\tilde{\mathcal{X}}$ decomposition can be expressed as

$$\tilde{\mathcal{X}}=\mathcal{G}\times {}_1{\mathit{U}}_1\times {}_2{\mathit{U}}_2\times {}_3{\mathit{U}}_3,$$

(19)

where $\mathcal{G}$ means the tensor of the core, which can be obtained by HOSVD, ${\mathit{U}}_n\in {\mathbb{C}}^{I_n\times I_n}$$(n=1,2,3)$, and the left singular value can be obtained by expanding the mode n of the tensor, which can be expressed as: ${\left[\tilde{\mathcal{X}}\right]}_n={\mathit{U}}_n{\Lambda }_n{\mathit{V}}_n^H$.

Next, we determine the signaling subspace ${\mathit{U}}_s$ of the tensor. In addition, it is assumed that there are L multipath signals generated in the room, so the rank of $\tilde{\mathcal{X}}$ is L. The tensor-based subspace is obtained by truncating the HOSVD and can be expressed as

$${\tilde{\mathcal{X}}}_t={\mathcal{G}}_t\times {}_1{\mathit{U}}_{t1}\times {}_2{\mathit{U}}_{t2},$$

(20)

where ${\mathcal{G}}_t$ is the truncated kernel tensor, and it can be expressed as ${\mathcal{G}}_t=\mathcal{X}\times {}_1{\mathit{U}}_{t1}^H\times {\mathit{U}}_{t2}^H\times {}_3{\mathit{U}}_{t3}^H$. ${\mathit{U}}_{tn}$ is the eigenvector corresponding to the L maximum eigenvalues of ${\mathit{U}}_n(n=1,2,3)$. Then, substitute $\mathcal{G}$ into Equation (20). Next, the ${\tilde{\mathcal{X}}}_t$ can be expressed as

$${\tilde{\mathcal{X}}}_t=\tilde{\mathcal{X}}\times {}_1\left({\mathit{U}}_{t1}{\mathit{U}}_{t1}^H\right)\times {}_2\left({\mathit{U}}_{t2}{\mathit{U}}_{t2}^H\right)\times {}_3{\mathit{U}}_{t3}^H$$

(21)

The tensor-based subspace can be expressed as

$${\mathit{U}}_s={\left[{\tilde{\mathcal{X}}}_t\right]}_3^T=\left({\mathit{U}}_{t2}{\mathit{U}}_{t2}^H\otimes {\mathit{U}}_{t1}{\mathit{U}}_{t1}^H\right){\left[\mathit{x}\right]}_{\left(3\right)}^T{\mathit{U}}_{t3}^{*},$$

(22)

where ${\left[\tilde{\mathcal{X}}\right]}_{\left(3\right)}={\mathit{U}}_3{\Lambda }_3{\mathit{V}}_3^{\mathrm{H}}$. Similarly, ${\left[\tilde{\mathcal{X}}\right]}_{\left(3\right)}^{\mathrm{T}}\approx {\mathit{V}}_{t3}^{*}{\Lambda }_{t3}{\mathit{U}}_{t3}^{\mathrm{T}}$ is obvious. Then, substitute ${\left[\tilde{\mathcal{X}}\right]}_{\left(3\right)}^{\mathrm{T}}$ into Equation (21), and ${\mathit{U}}_s$ can be expressed as

$$\begin{array}{c}\hfill {\mathit{U}}_s=\left({\mathit{U}}_{t2}{\mathit{U}}_{t2}^H\otimes {\mathit{U}}_{t1}{\mathit{U}}_{t1}^H\right){\mathit{V}}_{t3}^{*}{\Lambda }_{t3}\end{array}$$

(23)

After the orthogonal projection, the noise subspace can be obtained, which can be expressed as

$$\begin{array}{c}\hfill {\mathit{U}}_{\mathrm{noise}}{\mathit{U}}_{\mathrm{noise}}^H={\mathit{I}}_{N{n}_{t}c\left(M-n_{t}c\right)}-{\mathit{U}}_o{\mathit{U}}_o^H\end{array}$$

(24)

Finally, after obtaining the noise subspace, we can construct the spectral peak search function to search the AoA of the sources. The spectrum peak search function is

$$\mathit{f}\left(\theta ,\tau \right)=\frac{1}{{\left[\tilde{\mathit{a}}\left(\theta ,\tau \right)\otimes \tilde{\mathit{b}}\left(\tau \right)\right]}^H\left[{\mathit{I}}_{N{n}_{tc}\left(M-n_{tc}\right)}-{\mathit{U}}_o{\mathit{U}}_o^H\right]\left[\tilde{\mathit{a}}\left(\theta ,\tau \right)\otimes \tilde{\mathit{b}}\left(\tau \right)\right]}$$

(25)

However, this is a two-dimensional spectral peak search function, which requires a lot of computation.

3.6. RD-MUSIC

The sixth part is to introduce how to achieve the estimation of the AoA with RD-MUSIC. Firstly, in order to facilitate the subsequent operation of the AoA and ToF, we set ${\mathit{T}}_1\in {\mathbb{C}}^{d_1\times d_2},{\mathit{T}}_2\in {\mathbb{C}}^{d_2\times d_3},{\mathit{T}}_3\in {\mathbb{C}}^{d_4\times d_5}$, and ${\mathit{T}}_4\in {\mathbb{C}}^{d_5\times d_6}$. We can obtain

$$\left({\mathit{T}}_1{\mathit{T}}_{\mathit{s}}\right)\otimes \left({\mathit{T}}_3{\mathit{T}}_4\right)=\left({\mathit{T}}_1\otimes {\mathit{T}}_3\right)\left({\mathit{T}}_2\otimes {\mathit{T}}_4\right)$$

(26)

$${\mathit{T}}_1\otimes {\mathit{T}}_3=\left({\mathit{T}}_1{\mathit{I}}_{d_2}\right)\otimes \left({\mathit{I}}_{d_4}{\mathit{T}}_3\right)=\left({\mathit{T}}_1\otimes {\mathit{I}}_{{\mathit{d}}_4}\right)\left({\mathit{I}}_{d_2}\otimes {\mathit{T}}_3\right)$$

(27)

Because of $\mathit{a}(\theta ,\tau )\in {\mathbb{C}}^{N{n}_{tc}\times 1}$ and $\mathit{b}\left(\tau \right)\in {\mathbb{C}}^{\left(M-n_{tc}\right)\times 1}$, we can simplify it as follows: $\widehat{\mathit{a}}\left(\theta \right)\in {\mathbb{C}}^{N\times 1}$ and $\widehat{\mathit{b}}\left(\tau \right)\in {\mathbb{C}}^{n_{tc}\left(M-n_{tc}\right)\times 1}$. We can obtain the simplified $\mathit{a}(\theta ,\tau )\otimes \mathit{b}\left(\tau \right)$, and it can be expressed as

$$\begin{array}{c}\hfill \tilde{\mathit{a}}\left(\theta ,\tau \right)\otimes \tilde{\mathit{b}}\left(\tau \right)=\widehat{\mathit{a}}\left(\theta \right)\otimes \widehat{\mathit{b}}\left(\tau \right)=\left[\widehat{\mathit{a}}\left(\mathit{\theta }\right)\otimes {\mathit{I}}_{n_{tc}\left(M-n_{tc}\right)}\right]\widehat{\mathit{b}}\left(\mathit{\tau }\right)\end{array}$$

(28)

Then, the spectrum peak search function can be simplified as

$$\begin{array}{cc}\hfill \mathit{f}(\theta ,\tau )& ={[\widehat{\mathit{a}}\left(\theta \right)\otimes \widehat{\mathit{b}}\left(\tau \right)]}^H{\mathit{U}}_{\mathrm{orth}}[\widehat{\mathit{a}}\left(\theta \right)\otimes \widehat{\mathit{b}}\left(\tau \right)]\hfill \\ \hfill & =\widehat{\mathit{b}}{\left(\tau \right)}^H\left[\widehat{\mathit{a}}\left(\theta \right)\otimes {\mathit{I}}_{n_{tc}\left(M-n_{tc}\right)}\right]{}^H{\mathit{U}}_{\mathrm{orth}}\left[\widehat{\mathit{a}}\left(\theta \right)\otimes {\mathit{I}}_{n_{tc}\left(M-n_{tc}\right)}\right]\widehat{\mathit{b}}\left(\tau \right)\hfill \\ \hfill & =\widehat{\mathit{b}}{\left(\tau \right)}^H\mathit{Q}\left(\theta \right)\widehat{\mathit{b}}\left(\tau \right)\hfill \end{array}$$

(29)

where $\mathit{Q}\left(\theta \right)=\left[\widehat{\mathit{a}}\left(\theta \right)\otimes {\mathit{I}}_{n_{tc}\left(M-n_{tc}\right)}\right]{}^H{\mathit{U}}_{orth}\left[\widehat{\mathit{a}}\left(\theta \right)\otimes {\mathit{I}}_{n_{tc}\left(M-n_{tc}\right)}\right]$. The above equation is a quadratic optimization problem, and we eliminate the mundane solution of $\widehat{\mathit{a}}\left(\theta \right)=0$ with ${\mathit{e}}_1^H\widehat{\mathit{a}}\left(\theta \right)=1$ where ${\mathit{e}}_1={\left[\begin{array}{cccc}1\hfill & 0\hfill & \cdots \hfill & 0\hfill \end{array}\right]}^T\in {\mathbb{C}}^{n_{tc}\left(M-n_{tc}\right)\times 1}$.

Then, in order to realize the reduced dimensional estimation, it is necessary to replace the $\widehat{\mathit{b}}\left(\tau \right)$-relativistic equation with the $\widehat{\mathit{a}}\left(\theta \right)$-relativistic equation. This optimization problem can be reformulated as follows

$$\underset{\theta ,\tau }{min}\widehat{\mathit{b}}{\left(\tau \right)}^H\mathit{Q}\left(\theta \right)\widehat{\mathit{b}}\left(\tau \right),\mathrm{s}.\mathrm{t}.{\mathit{e}}_1^H\widehat{\mathit{a}}\left(\theta \right)=1$$

(30)

The cost function can be expressed as

$$\mathit{L}(\theta ,\tau )=\widehat{\mathit{b}}{\left(\tau \right)}^{\mathit{H}}\mathit{Q}\left(\theta \right)\widehat{\mathit{b}}\left(\tau \right)-\mathit{\lambda }\left({\mathit{e}}_1^{\mathit{H}}{\mathit{a}}_x\left(\mathit{v}\right)-1\right)$$

(31)

where $\lambda$ is constant and derivative of $\widehat{\mathit{b}}\left(\mathit{\tau }\right)$ for Equation (31). Then, we obtain the expressions as follows

$$\frac{\partial }{\partial \widehat{\mathit{b}}\left(\tau \right)}\mathit{L}(\theta ,\tau )=2\mathit{Q}\left(\theta \right)\widehat{\mathit{b}}\left(\tau \right)+\lambda {\mathit{e}}_1=0$$

(32)

According to the above equation, we can obtain $\widehat{\mathit{b}}\left(\tau \right)=\mu {\mathit{Q}}^{-1}\left(\theta \right)e_1$, where $\mu$ is a constant. Because of ${\mathit{e}}_1^H\widehat{\mathit{a}}\left(\mathit{\tau }\right)=1$, we can obtain $\mu =\frac{1}{{\mathit{e}}_1^H\mathit{Q}{\left(\theta \right)}^{-1}{\mathit{e}}_1}$. The expression for $\widehat{\mathit{b}}\left(\mathit{\tau }\right)$ after simplification can be expressed as

$$\widehat{\mathit{b}}\left(\tau \right)=\frac{\mathit{Q}{\left(\theta \right)}^{-1}{\mathit{e}}_1}{{\mathit{e}}_1^H\mathit{Q}{\left(\theta \right)}^{-1}{\mathit{e}}_1}$$

(33)

Finally, the dimensionality reduction can be realized by substituting Equation (33) into Equation (29). The expressions for $\tilde{\theta }$ can be expressed as

$$\begin{array}{cc}\hfill \tilde{\theta }& =argmax\mathit{f}(\theta ,\tau )\hfill \\ \hfill & =argmin\widehat{\mathit{b}}{\left(\tau \right)}^H\mathit{Q}\left(\theta \right)\widehat{\mathit{b}}\left(\tau \right)\hfill \\ \hfill & =argmin{\mathit{e}}_1^{-1}\mathit{Q}\left(\theta \right){\mathit{e}}_1^{-H}\hfill \\ \hfill & =argmax{\mathit{e}}_1^H\mathit{Q}{\left(\theta \right)}^{-1}{\mathit{e}}_1\hfill \end{array}$$

(34)

3.7. Triangular Positioning

After obtaining the AoA of the three AP points, triangulation can be performed. Triangulation is performed as follows in Figure 4.

In Figure 4, the position can be calculated, and the coordinates of the center of mass can be expressed as

$$O\left(x_0,y_0\right)=\left(\frac{1}{3}\left(x_1+x_2+x_3\right),\frac{1}{3}\left(y_1+y_2+y_3\right)\right)$$

(35)

After we obtain the three AoAs relative to the AP point, we can intersect a triangle range from the three AoAs directions based on Figure 4. According to the coordinates of the three points of the triangle, the coordinates of the quality points of this triangle can be found by Equation (35), and the positioning can be completed.

4. Summary of Algorithm Steps

First, we will use a large tensor to store the multi-packet CSI data. Then, smaller tensors will be constructed by tensor compression technology, and the signal subspace will be obtained by HOSVD. Finally, the AoA will be obtained by the MUSIC algorithm. The resulting AoA will be obtained by the MUSIC algorithm. The complete algorithm pseudo-code is shown in Algorithm 1.

Algorithm 1 HOSVD Method

1: Input: Multiple packet CSI matrix: $\left[\begin{array}{cccc}{\mathit{H}}_1& {\mathit{H}}_2& \cdots & {\mathit{H}}_I\end{array}\right]$ 2: Output: AoAs target: ${\theta }_l\left(1\le l\le L\right)$; 3: function HOSVD($\left[\begin{array}{cccc}{\mathit{H}}_1& {\mathit{H}}_2& \cdots & {\mathit{H}}_I\end{array}\right]$);         $\mathit{X}$ is constructed by multi-packet fusion         $\mathit{X}$ = []         for i = 1: I do           $\mathit{X}$ is reconstructed as ${\mathit{X}}_i$ by borrowing subcarriers;           $\mathit{X}=\left[\mathit{X}{\mathit{X}}_i\right]$         end for         Construct the tensor $\mathcal{X}$ with $\mathit{X}$;         for $o=1:M\ast n_{tc}$ do           $\mathcal{X}(o,:,:)$ = $\mathit{X}((o-1)\ast (N-n_{tc})+1:o\ast (N-n_{tc}),:))$         end for         Tensor $\tilde{\mathcal{X}}$ is obtained by tensor compression technique;         $\mathit{W}\leftarrow mode3svd\left(\mathcal{X}\right)$         $\tilde{\mathcal{X}}=\mathcal{X}\ast \mathit{W}$         for $k=1:M\ast n_{tc}$ do            $\tilde{\mathcal{X}}(k,:,:)$ = $\tilde{\mathcal{X}}((k-1)\ast (N-n_{tc})+1:k\ast (N-n_{tc}),:))$         end for         ${\mathit{U}}_{\mathrm{noise}}{\mathit{U}}_{\mathrm{noise}}^H$ is obtained by HOSVD;         ${\mathit{U}}_{noise}{\mathit{U}}_{noise}^H\leftarrow HOSVD\left(\tilde{\mathcal{X}}\right)$         $f\left({\theta }_l,{\tau }_l\right)=\frac{1}{{\left[\mathit{a}\left({\theta }_l,{\tau }_l\right)\otimes \mathit{b}\left({\tau }_l\right)\right]}^H\left[{\mathit{I}}_{n_{tc}\left(M-n_{tc}\right)}-U_o{U}_o^H\right]\left[\mathit{a}\left({\theta }_l,{\tau }_l\right)\otimes \mathit{b}\left({\tau }_l\right)\right]}$         One-dimensional spectrum peak search function is constructed;         ${\theta }_l=argmax{\mathit{e}}_1^H\mathit{Q}{\left({\theta }_l\right)}^{-1}{\mathit{e}}_1$ 4: end function 5: Achieve AoAs estimation.

5. Simulation and Experimentation

In order to verify the proposed method, simulation experiments and practical experiments are performed in this section.

In this experiment, we do not make hardware changes and only use commercial NICs and routers for the experiment. We choose the 5G band that has a bandwidth of 40 MHZ, and this choice will have a better estimation performance than using a 20 MHZ bandwidth.

5.1. Verification of Algorithm Efficiency

The first part is an ablation experiment with tensor compression and RD-MUSIC to verify the effectiveness of both computation methods. The second part is a simulation experiment of the HOSVD algorithm under different conditions. The first aspect is to test the ability of the HOSVD algorithm to estimate different multipath targets. The second aspect is to test the effect of different snapshots on the estimation accuracy of the HOSVD algorithm. The number of times a sensor makes an adoption is called a snapshot. The more samples taken, the better the estimate will be. The third aspect is to test the estimation capability of the HOSVD algorithm with different SNRs. The ratio of signal power to noise power is defined as the SNR, and a ratio higher than 1:1 (SNR > 0dB) means that there is more signal than noise.

According to

Table 2. Simulation and experimental parameters.

Antenna numbers	M	3
Channel number	-	161
Center frequency	-	5.805 GHZ
Packet numbers	I	150
$n_{tc}$	-	3
Antenna distance	d	2.5 cm
Multipath numbers	L	3, 4, 5
Bandwidth	-	40 MHZ
Subcarrier numbers	N	30

, we used the parameters of the chart and conducted simulation experiments using different snapshots. Simulation experiments are performed to estimate the AoA and ToF for each of the five targets. In this simulation experiment, the following four cases are performed separately:

A: Tensor compression and dimensionality RD-MUSIC are not performed.

B: Tensor compression is performed, but no dimensionality RD-MUSIC is performed.

C: Tensor compression is not performed, but dimensionality RD-MUSIC is performed.

D: Both tensor compression and dimensionality RD-MUSIC are performed.

According to

Table 3. Computational time of the algorithm in 4 cases.

Snapshots	A	B	C	D
150	1.23 s	1.09 s	0.39 s	0.22 s
200	1.57 s	1.11 s	0.70 s	0.23 s
300	2.28 s	1.11 s	1.53 s	0.23 s
400	3.35 s	1.14 s	2.74 s	0.24 s
500	7.18 s	1.15 s	4.25 s	0.24 s

, first is the comparison of Case A and Case B. After tensor compression is performed, the computation time does not increase significantly due to the increase in the number of packets. This is because it is a dimensionality reduction process on the data dimension. Then, there is a comparison between Case A and Case B. After RD-MUSIC is performed, the computation time increases dramatically with the increase in data packets. However, with the same number of packets, the computation time of Case C is significantly smaller than that of Case A. This is because RD-MUSIC is making a two-dimensional traversal search reduced to a one-dimensional traversal search. Last is the comparison between Case A and Case D. With both tensor compression and RD-MUSIC, the computation time can be reduced significantly and does not increase with the increase in the number of packets.

According to Table 2, we used the parameters of the chart and conducted simulation experiments using 150 snapshots. The tested AoA and ToF are (−68.4°, 4.0 × ${10}^{-8}$ ns), (−33.1°, 1.8 × ${10}^{-8}$ ns), (10.1°, 2.8 × ${10}^{-8}$ ns), (40.3°, 7.2 × ${10}^{-8}$ ns), and (59.8°, 5.8 × ${10}^{-8}$ ns) with an SNR = −5 dB, and the effect is shown in Figure 5.

As Figure 5 shows, there are three, four, and five obvious and accurate spectral peaks on the spectral peak graph, respectively, when the multipath number L = 3, L = 4, and L = 5. This can prove that the HOSVD method can successfully and accurately estimate the angle when L ≤ 5. We used the parameters of Table 2 and conducted simulation experiments using different snapshots. The tested AoA and ToF are (−68.4°, 4.0 × ${10}^{-8}$ ns), (−33.1°, 1.8 × ${10}^{-8}$ ns), (10.1°, 2.8 × ${10}^{-8}$ ns), (40.3°, 7.2 × ${10}^{-8}$ ns), and (59.8°, 5.8 × ${10}^{-8}$ ns) with an SNR = −5 dB, and the effect is shown in Figure 6.

As Figure 6 shows, only one of the four multipath targets can be estimated when the snapshots = 15. There are four multipath targets that can be estimated completely and accurately when the number of snapshots increases to 150. This can prove that multi-packet fusion can improve the accuracy of estimating the number of targets and the accuracy of the target angle. When the number of snapshots increases from 150 to 300, it can be seen from the peak grap that the increase in the number of snapshots has a positive gain on the estimation accuracy. However, the gain is smaller when the number of snapshots increases to a certain extent. We used the parameters of Table 2 and conducted simulation experiments using 150 snapshots. The tested AoA and ToF are (−68.4°, 4.0 × ${10}^{-8}$ ns), (−33.1°, 1.8 × ${10}^{-8}$ ns), (10.1°, 2.8 × ${10}^{-8}$ ns), (40.3°, 7.2 × ${10}^{-8}$ ns), and (59.8°, 5.8 × ${10}^{-8}$ ns) with different SNRs, and the results are shown in Figure 7.

As Figure 7 shows, with the increase in the SNR, the estimation accuracy results of the proposed algorithm are continuously improved. When the SNR = −10 dB, the HOSVD method still has the ability to estimate the target angle. When the SNR = −5 dB, the algorithm can accurately estimate the angle of five multipath targets.

5.2. Comparison of Algorithm Simulation

According to Table 2, for HOSVD and Music, the Tensor Decomposition algorithm estimated the AoA and ToF with 150 snapshots. The tested AoA and ToF are $(-68.4$°,$4.0\times {10}^{-8}$ ns), $(-33.1$°$,1.8\times {10}^{-8}$ ns), $(10.1$°$,2.8\times {10}^{-8}$ ns), $(40.3$°$,7.2\times {10}^{-8}$ ns), and $(59.8$°$,5.8\times {10}^{-8}$ ns) with an SNR = −5 dB, and the spectral peak plots are as follows.

From Figure 8, we can see that Tensor Decomposition cannot accurately estimate the five multipath targets with the low SNR, and the spectral peaks are not sharp enough to make the accuracy high. HOSVD and MUSIC can still accurately estimate the number of the five multipath targets and have a high accuracy with the low SNR. However, the spectral peaks searched by the MUSIC algorithm are not sharp, which can prove that the HOSVD method is more suitable for indoor complex environments.

We perform 500 Monte Carlo experiments on HOSVD, the Tensor Decomposition algorithm, and the MUSIC algorithm respectively to calculate the Root Mean Squared Error (RMSE) [24,25] for estimating five different target angles with different SNRs. Monte Carlo experiments are performed using the parameters within Table 2. The tested AoA and ToF are $(-68.4$°,$4.0\times {10}^{-8}$ ns), $(-33.1$°$,1.8\times {10}^{-8}$ ns), $(10.1$°$,2.8\times {10}^{-8}$ ns), $(40.3$°$,7.2\times {10}^{-8}$ ns), $(59.8$°$,5.8\times {10}^{-8}$ ns). The SNR is set to −5 dB to 20 dB. The expression of the RMSE can be expressed as

$$RMS{E}_{AoA}=\sqrt{\frac{1}{P}\frac{1}{L}\sum _{p=1}^P\sum _{l=1}^L{({\theta }_l-{\widehat{\theta }}_{l,p})}^2}$$

(36)

where ${\theta }_l$ is the l-th true DOA. ${\widehat{\theta }}_{l,p}$ is the angle on the l-th path measured by the p-th Monte Carlo experiment. For example, the error between the estimated AoA and the correct AoA for each loop was recorded when using an algorithm with SNR = −5 dB. The errors obtained under this loop were squared and summed to obtain a sum of errors. Finally, the sum of errors is divided by the number of loops and the total number of targets to obtain the value of RMSE for this algorithm with an SNR = −5 dB. From Figure 8c, the Tensor Decomposition algorithm [24] cannot estimate the exact target number with an SNR = −5 dB. The method of the spectral peak search used in the text is abandoned, and the method of phase extraction is chosen for estimation.

As seen from Figure 9, HOSVD has a lower RMSE than the MUSIC algorithm and the Tensor Decomposition algorithm for all SNRs. This shows that HOSVD has a theoretically better performance in complex environments with a low SNR.

5.3. Experiments

After verification of the simulation experiments, we now perform field tests. The receiver is a computer with an intel 5300 (Manufacturer: Shenzhen Fenvi Technology Co. The city and country of source of the equipment are: Shenzhen and China, respectively.) NIC, operating system version unbuntu 14.04, and the intel 5300 NIC has three ULA antennas. The transmitter is a Xiaomi router with 2 × 2G antennas and 2 × 5G antennas; we mainly choose the 5G antenna as the transmitting antenna.

We first installed the CSI tool on the desktop computer’s unbuntu system. As seen from Figure 10, the desktop computer is equipped with an intel 5300 NIC, and an external antenna with SMA cable is used as the receiving end to facilitate the collection of CSI data. The Xiaomi 4a router is the transmitter. We set the router in the 161 band of the 5G signal with a center frequency of 5.805 GHZ and a bandwidth of 40 MHZ. We connect the router wirelessly with a computer and open two terminals, one to start the CSI tool and the other to continue data interaction with the router, thus enabling data collection.All data processing in the experiments was performed using MATLAB R2021a.

As Figure 11 shows, the tensor $\mathcal{X}$ stores multi-packet CSI data, but one of the dimensions of the tensor $\mathcal{X}$ is smaller than the multipath number L. The new tensor $\tilde{\mathcal{X}}$ needs to be reconstructed again by Equation (10). Then, each dimension of the tensor $\tilde{\mathcal{X}}$ is larger than L. Decomposition of this large tensor $\tilde{\mathcal{X}}$ with HOSVD takes a long time. In this case, the tensor can be compressed by means of a compression matrix. With this method, it is possible to preserve most of the information of the signal and also to keep the tensor as small as possible. On the other hand, the intel 5300 NIC is not sufficient for ToF resolution and cannot be used to estimate the distance from the transmitter to the receiver. Therefore, only the AoA needs to be estimated. We use the RD-MUSIC method to estimate the required AoA, which can change the time-consuming two-dimensional spectral peak search into a one-dimensional spectral peak search and reduce the computation time significantly.

5.4. Experimental Testing in Different Scenarios

We chose the equipment mentioned in the previous section to conduct localization experiments on the playground and laboratory. First, we conduct the positioning experiment on the playground. The experimental area is an open area of 6 m × 8 m. The experimental equipment is placed as follows in Figure 12.

As shown in Figure 12, three AP points are placed, where the red point is the test point, and data are collected 10 times on each test point, each time for 150 packets. After performing data processing with different algorithms, we chose the HOSVD algorithm, Tensor Decomposition algorithm [24], MoLA algorithm [19], and SpotFi algorithm [6] for comparison. The Cumulative Distribution Function (CDF) plot of the overall error of the AoA is obtained as follows in Figure 13.

According to Figure 13, considering the estimation error at a ratio of 0.667, it can be seen that the HOSVD algorithm, Tensor Decomposition algorithm, MoLA algorithm, and SpotFi algorithm can achieve AoA estimation errors of 3.25°, 4.57°, 6.87°, and 7.45°. The HOSVD algorithm has better estimation performances at test locations with a high SNR. There is no significant multipath signal in the spatial playground; therefore, the collected CSI data will only contain a direct path angle. The calculated angles deviate slightly from the true angles, which proves that the measured data and the algorithms are convincing. The CDF of the distance gap between the calculated position and the true test point position for triangulation using Equation (35) are plotted as follows in Figure 14.

According to Figure 14, localization experiments were carried out on the playground. The HOSVD, the Tensor Decomposition algorithm, the MoLA algorithm, and the SpotFi algorithm could achieve localization errors of 0.39 m, 0.50 m, 0.72 m, and 0.85 m. The second experimental location is the laboratory, which has many tables and pieces of equipment, so the multipath effect is evident here. The distribution scheme of experimental AP points and test points is in Figure 15.

As shown in Figure 15, there are three AP points set up in the laboratory. The red dots are the test points, and there are 23 test points in total. At each test point, 10 sets of data are collected; each set contains 150 packets of data. The heights of the receiver and transmitter are at the same level. The true position and orientation of the ground was confirmed by using a building protractor and the geometric relationship between the transmitter and receiver. In the experiments, the proposed algorithm was compared with the Tensor Decomposition algorithm [24], the MoLA algorithm [19], and the SpotFi algorithm [6] to evaluate the parameter estimation accuracy and localization accuracy. Firstly, the overall AoA accuracy and localization accuracy of the different algorithms are evaluated in the same scenario. Secondly, we compare the program execution times. Everything is the same in each experiment except that the algorithms were different.

In the experiment, 10 sets of CSI data are collected from the same AP point. The true angle of this AP point to the receiver was 33.3°. The AoA results are obtained by calculating these 10 sets of CSI data by the HOSVD, Tensor Decomposition, MoLA, and SpotFi algorithms, respectively. The results are shown in

Table 4. Experimental results comparison.

True Angle	HOSVD	Tensor Decomposition	MoLA	SpotFi
33.3°	38°	41°	43°	51°
33.3°	37°	40°	42°	50°
33.3°	39°	42°	45°	48°
33.3°	38°	41°	42°	53°
33.3°	27°	25°	21°	13°
33.3°	38°	41°	44°	49°
33.3°	39°	40°	44°	49°
33.3°	40°	41°	45°	52°
33.3°	29°	28°	25°	16°
33.3°	26°	25°	22°	11°
average error	5.38°	7.48°	10.38°	17.88°

. According to Table 4, the HOSVD and Tensor Decomposition algorithms still retain certain robustnesses in the complex laboratory environment, and the average angle errors of the 10 measurements are 5.38° and 7.6°, respectively. The MoLA and SpotFi algorithms are single snapshot two-dimensional algorithms with poor robustnesses in complex environments, and the average angle errors of the 10 measurements are 10.38° and 17.88°, respectively. The MoLA algorithm employs source estimation and I-MUSIC to decorrelate the signal as much as possible and improve the accuracy of the estimation. According to Table 4, the HOSVD algorithm exhibits a better accuracy of the estimation, which also indicates that the estimation accuracy in complex environments can be improved by tensor storage of multi-packet data.

According to Figure 16, considering the estimation error when the ratio is 0.667, it can be seen that the AoA estimation accuracies of the HOSVD algorithm, Tensor Decomposition algorithm, MoLA algorithm, and SpotFi algorithm can reach 5.81°, 7.62°, 12.33°, and 17.12°, respectively, which proves that HOSVD has the best angle estimation in the complex laboratory indoor conditions. The CDF diagram of the distance between the true position and the calculated position with triangulation in Equation (35) is shown in Figure 17.

The results of the localization experiments conducted in the laboratory can be seen in Figure 17. The HOSVD algorithm, Tensor Decomposition algorithm, MoLA algorithm, and SpotFi algorithm can achieve localization errors of 0.68 m, 0.89 m, 1.21 m, and 1.33 m, considering an estimation error at a ratio of 0.667. The laboratory environment is more complex. MoLA and SpotFi use single packet estimation, resulting in a larger overall estimation error. Both the Tensor Decomposition algorithm and HOSVD use packet fusion, leading to a smaller overall estimation error. The final results show that HOSVD performs better in terms of its estimation performance with a low SNR.

5.5. Calculation Time

The short calculation time of the angle estimate benefits real-time localization. The HOSVD algorithm, Tensor Decomposition algorithm [24], MoLA algorithm [19], and SpotFi algorithm [6], respectively, perform angle estimations for different numbers of packets. The specific computation times are shown in

Table 5. Calculation time with different packet counts.

Number of Data Packets	HOSVD	Tensor Decomposition	MoLA	SpotFi
50	0.38 s	1.04 s	8.24 s	7.94 s
70	0.40 s	1.05 s	9.34 s	9.01 s
90	0.41 s	1.06 s	12.14 s	11.57 s
110	0.41 s	1.18 s	13.74 s	13.40 s
130	0.41 s	1.25 s	15.64 s	15.37 s
150	0.42 s	1.34 s	17.94 s	17.27 s

.

According to Table 5, on the one hand, the increase in packets does not cause a significant time change to HOSVD due to the tensor compression technology. The time spent by the Tensor Decomposition algorithm is mainly for the 2D spectrum peak search. MoLA and SpotFi need to perform a spectral peak search on each packet, which takes a lot of time. The tensor constructed in this paper is also smaller in dimension than the one constructed by the Tensor Decomposition algorithm. As a result, the proposed algorithm takes the least amount of time when faced with large amounts of data. In addition, the computation time does not change rapidly as the amount of data grows.

6. Conclusions

In this paper, we propose a HOSVD-based algorithm for indoor WiFi array localization. This algorithm utilizes HOSVD to estimate the AoA of the target with higher accuracy while making full use of the multi-dimensional structural properties characteristics of CSI data. This algorithm uses data reconstruction, tensor compression, and RD-MUSIC to improve the efficiency of localization. As can be seen from the resulting figure of the simulation experiment, the HOSVD algorithm not only accurately estimates the precise number of sources but also has a smaller AoA error compared to the other algorithms. In addition, the computation time of the HOSVD algorithm is not only the shortest but also does not increase with the increase in the number of packets. The experimental results in the playground and laboratory show that the HOSVD algorithm has a lower AoA error and localization error compared to other algorithms. In conclusion, the HOSVD-based WiFi array indoor localization algorithm has some advantages in terms of accuracy and computation time compared to other algorithms.