# RIS-Aided Near-Field Localization and Channel Estimation for the Sub-
Terahertz System

Yijin Pan, Cunhua Pan, Shi Jin and Jiangzhou Wang Y. Pan, C. Pan and S. Jin
are with the National Mobile Communications Research Laboratory, Southeast
University, Nanjing 211111, China. (Corresponding author: Cunhua Pan) J. Wang
is with the School of Engineering, University of Kent, Canterbury, CT2 7NT,
UK.

###### Abstract

The low hardware cost makes ultra-large (XL) reconfigurable intelligent
surface (RIS) an attractive solution for the performance enhancement of
localization and communication systems, but it results in near-field
propagation channels, especially for the high-frequency communication systems,
e.g. mmWave and sub-terahertz (subTHz) systems. This makes localization and
channel estimation in the near field of the RIS-aided subTHz system more
challenging. In this paper, we consider the spherical wavefront propagation in
the near field of the subTHz system with the assistance of a RIS. A near-field
channel estimation and localization (NF-JCEL) algorithm is proposed based on
the derived second-order Fresnel approximation of the near-field channel
model. To be specific, we first decouple the user equipment (UE) distances and
angles of arrival (AoAs) through a down-sampled Toeplitz covariance matrix, so
that the vertical and azimuth AoAs in the array steering vectors can be
estimated separately with low complexities. Then, the UE distance can be
estimated by the simple one-dimensional search, and the channel attenuation
coefficients are obtained through the orthogonal matching pursuit (OMP)
method. Simulation results validate the superiority of the proposed NF-JCEL
algorithm to the conventional far-field algorithm, and show that higher
resolution accuracy can be obtained by the proposed algorithm.

###### Index Terms:

Near-field transmission, localization, intelligent reflecting surface (IRS),
reconfigurable intelligent surface (RIS), terahertz (THz) communications.

##  I Introduction

The sixth-generation (6G) mobile networks are expected to cater to various
Internet of Things (IoT) applications such as healthcare, smart homes and
intelligent manufacture[1]. These applications and services entail decision
making on the basis of high-resolution location data from the involved
collections of sensors, machines and other items. For 6G services, sub
teraherze (subTHz) spectrum, i.e. from 90 GHz to 300 GHz, has primarily been a
driving force for high-resolution sensing with centimeter-level positioning
accuracy[2, 3], resulting from its high frequency, sufficiently large
bandwidth and increased spectral efficiency[4, 5].

However, radio transmission at high frequency, for instance subTHz spectrum,
suffer from high path loss, bad propagation conditions and is easy to be
blocked. To cope with this issue, the reconfigurable intelligent surfaces
(RISs) has been regarded as a promising solution to overcome the blockage
issue of the subTHz localization system. By carefully tuning the phase shifts
of the reflecting elements, RIS can help constructively accumulate increased
power at the target receiver, and thus shows a great potential to construct
the intelligent radio environment. Thus, it was leveraged to help enhance
performance of the simultaneous wireless information and power transfer
(SWIPT) system [6, 7] and multicell network [8]. Specially, in localization
system, it is shown that the RIS can help make the localization more
accurate[9, 10]. It was shown in [11] that proper RIS phase shifts
configuration can also help improve the Cram$\rm{\acute{e}}$r-Rao lower bound
(CRLB) of the position/orientation estimation.

Conventionally, the UE localization problems aim to utilize the received
signal strength (RSS), the time-difference-of-arrival (TDoA)[12], and the
angle of arrival (AoA) of multiple sensors to distinguish the multiple signal
directions from different emitters. The main idea is that all the geometric
information for localization is included in these channel state information
(CSI) measurements. Thus, the RIS-assisted localization accuracy is highly
dependent on the necessity to acquire side channel measurements [13]. However,
the channels between the RIS and UEs, as well as the channel between the
access point (AP) and the RIS, are cascaded[14]. In addition, the RIS has
little signal processing capabilities, making channel estimation of RIS-
assisted channels a challenging task[15]. Many researches successfully managed
to obtain the cascaded channel gains, such as [16], which have fully exploited
the correlations between UEs and different paths to reduce the pilot overhead.
In the RIS-assisted mmWave localization system, beam training was investigated
in [17] to estimate AoAs and angles of departure (AoDs) of the line of sight
(LoS) paths. In [18], the RSS was utilized in the RIS aided localization
system, and the RIS phase shifts were optimized to minimize the weighted
probabilities of false localization. The impact of phase shifts and the number
of reflecting elements on the CSI parameter estimation, positioning error
bound (PEB) and orientation error bound (OEB) was studied in [19].

However, the above researches were mainly based on the assumption of far field
channel model, where the planar wavefront was adopted. Thanks to the low
hardware cost and reduced power consumption, a large number of passive and
reflecting elements can be integrated into an RIS panel, so that its size can
reach several meters. In some envisioned industrial IoT application scenarios,
the RIS panels can be installed to cover the entire roof and walls of a
building. However, with the extra-large RIS panel, the spherical wavefront
feature should be considered, which is the so-called the near-field effect.
According to [20], to maintain a maximum phase difference of $\pi/8$ rad, the
observation distance from the UE/AP to the RIS must be no less than the
Fraunhofer distance $\frac{2L^{2}}{\lambda}$, where $L$ and $\lambda$ denote
the maximal aperture of the RIS panel (i.e. the diagonal length of a
rectangular panel) and radio wavelength, respectively. For instance, for an
RIS panel with aperture $50$ cm, the Fraunhofer distance is $83.33$ m for the
electrical magnetism (EM) waves at $50$ GHz ($\lambda=6$ mm), and it is
increased to $166.67$ m for the radio frequency increasing to $100$ GHz
($\lambda=3$ mm) in THz band. Therefore, in this case, the UEs are very likely
in the near field of the RIS panel and the general spherical wavefront should
be considered. With the spherical wavefront radio, the emitted EM wave will
arrive at each reflecting element in the RIS array with different AoAs. This
indicates that all array elements share a common AoA on the same path in the
far-field model is no longer valid [21]. The near-field effect makes the
channel modeling and CSI estimation for localization more challenging, and the
near-field communications will profoundly degrade the localization
performance[22], as also to be proved in this paper.

The study on the RIS behavior in the near-field is just in its infancy. The
RIS reflected power behavior was analyzed and measured in [23] under near/far-
field conditions. The near-field channel modeling for active antenna arrays
and RIS was investigated in [24]. In [25], a generic communication model with
XL-array/surface was investigated with the consideration of the variations of
signal phase, amplitude and aperture across array elements. The channel model
mismatch was addressed in [26] by leveraging the MCRB (misspecified Cramer-Rao
lower bound) to lower bound the localization error using a simplified
mismatched model. A near-field codebook was developed in [27] for the XL-RIS
beam training by dividing the three-dimensional (3D) space into sampled points
in the x-y-z coordinate system. The authors of [28] demonstrated that a RIS
can serve as an anomalous mirror in the antenna array’s near field, by using
antenna theory to calculate the electric field of a finite-size RIS. As for
the radio localization approaches with RIS, there are only a few works with
the consideration of spherical wavefront in the near-field. In [29], the
Fisher information matrix (FIM) was analyzed for an uplink localization system
using a RIS-based lens, and the PEB and OEB were also evaluated by exploiting
the wavefront curvature. The CRLB for RIS-assisted localization was
investigated in [30] and the RIS phase shift was optimized to maximize the
received signal to noise ratio (SNR). Similar researches were also presented
in [31], where FIM was analyzed for the multipath-aided localization with only
RIS-assisted non line of sight (NLoS) links.

Most existing works such as [24, 25, 26] considered the near-field channel
modeling and system performance analysis but did not address the issue of
angle estimation and positioning algorithm in the near field. Although several
RIS-assisted localization schemes have been proposed for the near-field
case[29, 30, 32], they were focused on the estimation error bound analysis,
such as FIM/PEB/CRLB, and there is a lack of the specific positioning
algorithm. In addition, some works such as [33, 34] considered the uniform
linear array (ULA) of RIS. However, the proposed localization solutions in
these works cannot be directly applied to the UPA case, as only either the
azimuth or elevation AoA/AoD was considered. These researches are not
sufficient for the more general 3D localization system, and cannot be directly
applied to the RIS channel with a uniform planar array (UPA). In practice, RIS
is normally deployed as a panel, and thus the angles in both azimuth or
elevation directions should be jointly estimated, however, which is not a
simple extension to the linear case. According to our work, the down-sampled
Toeplitz covariance matrix should be carefully re-investigated, and angles in
different directions need to be further decoupled to avoid two-dimensional
grid searching, and these tasks for the UPA RIS system are not trivial. In
addition, the benefits of the adjustable phase shifts of the RIS panel are not
fully exploited in current work. For instance, the RIS phase shifts were
considered to improve the error bounds in [29, 30, 32, 33], and to maximize
the SNR in [34]. In fact, the RIS can be employed to efficiently increase the
received signal samples from the same pilot data, thereby it has the potential
to effectively improve the accuracy in the process of the localization
algorithm itself. Unfortunately, how this approach can be integrated into the
localization algorithm has not been addressed. Overall, it is still unclear
about how to efficiently obtain near-field UE locations and the involved CSI
parameters, such as the AoA/AoDs, in the subTHz localization system with RIS.

In this paper, we investigate the localization and CSI estimation scheme in
the near field of the subTHz system with the assistance of a RIS panel. Our
contributions are summarized as follows:

  * •

By considering the propagation of the spherical wavefront across the RIS
array, the near-field UE-RIS uplink channel is modeled. To reduce the
complexity introduced by the different spatial paths at each reflecting
element, the channel formulation based the second-order Taylor approximation
is derived, which involves both the UE distance and the AoAs. Then, the
cascaded AP-RIS-UE channel in the near field is presented111The simplified
approximation of near-field channel model is presented in conference
version[35], but the following localization and channel estimation algorithm
is only presented in this work..

  * •

The RIS training phase shifts and pilots are carefully designed to increase
the channel rank so that the channel covariance matrix can be obtained by the
least square (LS) estimation. Then, a down-sampled covariance matrix is
derived to decouple the UE distance and the AoAs separately.

  * •

Based on the down-sampled Toeplitz covariance matrix, the vertical and azimuth
AoAs are estimated separately. Then, the UE distance can be estimated by the
simple one-dimensional search, and the orthogonal matching pursuit (OMP)
method is leveraged to obtain the channel attenuation coefficients.

  * •

Simulation results validate the effectiveness of the proposed near field joint
channel estimation and localization (NF-JCEL) algorithm. Compared with the
far-field algorithm, the proposed algorithm shows superiority in terms of
localization and channel estimation accuracies.

The remainder of this paper is organized as follows: Section II derives the
near/far-field RIS-aided channel model and formulates the channel estimation
problem. Section III develops the detailed algorithm for joint channel and
localization estimation. In Section IV, the simulation results are presented
to show the performance gain and the impact of system parameters, and Section
V concludes the paper.

Notation: For a vector $\bm{x}$, $|\bm{x}|$ and $(\bm{x})^{T}$ denote its
Euclidean norm and its transpose, respectively. For matrix $\bm{A}$,
$\bm{A}^{H}$ $\bm{A}^{-1}$ and $\bm{A}^{{\dagger}}$ represent the conjugate
transpose, the inverse and Moore-Penrose pseudoinverse operator, respectively.
$\mathbb{C}^{M\times N}$ denotes the set of $M\times N$ complex matrix.
$\text{diag}(\bm{x})$ represents the diagonal matrix $\bm{X}$ obtained from
vector $\bm{x}$. $\bm{a}\otimes\bm{b}$ represents the kronecker product of
$\bm{a}$ and $\bm{b}$.

##  II Channel Model

Consider the uplink transmission of a subTHz localization system, as shown in
Fig. 1. In the considered scenario, each UE is equipped with a single antenna,
the AP is located in the XOZ plane, and the number of UEs is denoted as $U$.
It is assumed that the direct links between the AP and the UEs do not exist
due to the blockage or unfavorable propagation environments. Thus, an RIS is
leveraged to construct the alternative AP-RIS-UE links for localization
service.

Figure 1: System model and communication scenario

The layouts of the AP antenna array and the RIS panel are shown in Fig. (2a)
and Fig. (2b), respectively. The receive antenna array of the AP is assumed to
be a ULA, and the number of AP’s antenna elements is denoted as $N_{A}$. The
distance between two antenna elements of the AP is denoted as $\Delta_{A}$. As
shown in Fig. (2b), the RIS is installed on the wall at the YOZ-plane, and the
RIS is equipped with $N_{R}$ reflecting elements, where the number of the
reflecting elements along the Y-axis and Z-axis are denoted as
$(2N^{Y}_{R}+1)$ and $(2N^{Z}_{R}+1)$, respectively. Thus,
$N_{R}=(2N^{Z}_{R}+1)\times(2N^{Y}_{R}+1)$. The distance between two
reflecting elements is denoted as $\Delta_{R}$.

Figure 2: The layout of the AP antenna arrays and RIS arrays.

Due to the short wavelength of the subTHz EM waves, second-bounce or more
reflections and scattering components are severely attenuated and thus they
are negligible. Similar to [36, 37, 30], we only consider the LOS paths in the
UE-RIS link and the AP-RIS link in the subTHz channels. Different from most of
the existing works, we consider the channel model with the near-field effects.
In the following, the formulation of signal propagation that accounts for the
spherical curvature of the wavefront is provided.

###  II-A UE-RIS link

We first consider the LoS path between UE $u$ and the RIS. With the spherical
wavefront, the receiving array response at the RIS’s $m$-th element is given
by

$\bm{a}_{R}[m]=\exp\left(-j\frac{2\pi}{\lambda}(d_{R,u}^{m}-d_{R,u}^{0})\right),$
(1)

where $d_{R,u}^{m}$ represents the distance from the RIS’s $m$-th element to
UE $u$, and $d_{R,u}^{0}$ represents the distance from the central reflecting
element to UE $u$. Then, the complex channel gain of UE $u$-RIS link is
denoted as

$\bm{h}_{R,u}=g_{R,u}\bm{a}_{R},$ (2)

where $g_{R,u}$ represents the complex channel attenuation coefficient, and
$\bm{a}_{R}$ collects the $N_{R}$ elements as
$\bm{a}_{R}=[\bm{a}_{R}[1],\cdots,\bm{a}_{R}[N_{R}]]^{T}$.

Denote the center of the RIS panel as the phase reference point, and its
coordinate is $(x_{R},y_{R},z_{R})$. Suppose that the $m$-th reflecting
element is located at $(x_{R},y_{R}+m_{y}\Delta_{R},z_{R}+m_{z}\Delta_{R})$,
where $m_{y}=-N^{Y}_{R},\cdots,0,N^{Y}_{R}$, and
$m_{z}=-N^{Z}_{R},\cdots,0,N^{Z}_{R}$. The index $m$ of the reflecting element
can be expressed as

$m=m_{z}(2N^{Y}_{R}+1)+m_{y},$ (3)

and the value range of $m$ is given as
$m=-N^{Z}_{R}(2N^{Y}_{R}+1)-N^{Y}_{R},\cdots,N^{Z}_{R}(2N^{Y}_{R}+1)+N^{Y}_{R}$.

Figure 3: Signal Propagation Model.

As shown in Fig. 3, the projection of distance $d_{R,u}^{m}$ on the XOY plane
is denoted as $l_{R,u}^{m}$, given by

$l_{R,u}^{m}=\sqrt{(m_{y}\Delta_{R})^{2}+(l_{R,u}^{0})^{2}-2\sin\theta_{u,0}m_{y}\Delta_{R}l_{R,u}^{0}},$
(4)

where $l_{R,u}^{0}$ denotes the projection of distance $d_{R,u}^{0}$ on the
XOY plane, given by

$l_{R,u}^{0}=d_{R,u}^{0}\cos\phi_{u,0}.$ (5)

Then, the distance $d_{R,u}^{m}$ is calculated as

$d_{R,u}^{m}=\sqrt{(d_{R,u}^{0}\sin\phi_{u,0}+m_{z}\Delta_{R})^{2}+(l_{R,u}^{m})^{2}}.$
(6)

Substituting (4) and (5) into (6), one obtains

$d_{R,u}^{m}=((m_{y}\Delta_{R})^{2}+(m_{z}\Delta_{R})^{2}+(d_{R,u}^{0})^{2}\\\
-2d_{R,u}^{0}m_{y}\Delta_{R}\cos\phi_{u,0}\sin\theta_{u,0}+2m_{z}\Delta_{R}d_{R,u}^{0}\sin\phi_{u,0})^{\frac{1}{2}}.$
(7)

Letting $y=\frac{m_{y}\Delta_{R}}{d_{R,u}^{0}}$, and
$z=\frac{m_{z}\Delta_{R}}{d_{R,u}^{0}}$, the distance between the $m$-th
reflecting element and UE $u$ can be represented as

$d_{R,u}^{m}=d_{R,u}^{0}(y^{2}+z^{2}+1-2y\cos\phi_{u,0}\sin\theta_{u,0}\\\
+2z\sin\phi_{u,0})^{\frac{1}{2}}\triangleq d_{R,u}^{0}F(y,z).$ (8)

In order to compute the distance $d_{R,u}^{m}$ in (8), we consider the
following simplifications. First of all, for ease of exposition, the first and
second order derivatives of $F(y,z)$ with respect to y and/or z are given by

$\frac{\partial F}{\partial
y}=\frac{y-\cos\phi_{u,0}\sin\theta_{u,0}}{F(y,z)},\frac{\partial F}{\partial
z}=\frac{z+\sin\phi_{u,0}}{F(y,z)},$ (9)
$\displaystyle\frac{\partial^{2}F}{\partial
y^{2}}=\frac{1+z^{2}+2z\sin\theta_{u,0}-(\cos\phi_{u,0}\sin\theta_{u,0})^{2}}{F(y,z)^{3}},$
(10) $\displaystyle\frac{\partial^{2}F}{\partial
z^{2}}=\frac{1+y^{2}-2y\cos\phi_{u,0}\sin\theta_{u,0}-(\sin\phi_{u,0})^{2}}{F(y,z)^{3}},$
(11) $\displaystyle\frac{\partial^{2}F}{\partial z\partial
y}=\frac{-(2y-2\cos\phi_{u,0}\sin\theta_{u,0})(z+\sin\phi_{u,0})}{2F(y,z)^{3}}.$
(12)

####  II-A1 Far field approximation

From [20], the far field of the RIS can be defined as the set of observation
distances $R$ that are greater than the Fraunhofer distance $R_{f}$, i.e.,

$R\geq R_{f}=\frac{2L^{2}}{\lambda}.$ (13)

This definition corresponds to the maximum phase difference, $\pi/8$, of radio
waves between any two reflecting elements. In this paper, we select the center
of the RIS panel as the reference point, as shown in (1), thus the aperture of
the RIS UPA should be the distance from panel vertex to the panel center as

$L=\sqrt{(M_{y}\Delta_{R})^{2}+(M_{z}\Delta_{R})^{2}}.$ (14)

As a result, when UE $u$ is located in the far field of the RIS panel, i.e.,

$d_{R,u}^{M}>2\frac{L^{2}}{\lambda}=2\frac{(M_{y}\Delta_{R})^{2}+(M_{z}\Delta_{R})^{2}}{\lambda},$

where $d_{R,u}^{M}$ denotes the distance of UE $u$ to the nearest RIS panel
vertex. Then, (7) can be approximated by the first-order Taylor expansion as

$\displaystyle d_{R,u}^{m}$ $\displaystyle\approx
d_{R,u}^{0}(F(y,z)|_{y=0,z=0}+\frac{\partial F}{\partial
y}|_{y=0}(y)+\frac{\partial F}{\partial z}|_{z=0}(z))$
$\displaystyle=d_{R,u}^{0}(1-\cos\phi_{u,0}\sin\theta_{u,0}y+\sin\phi_{u,0}z)$
$\displaystyle=d_{R,u}^{0}\left(\frac{m_{z}\Delta_{R}}{d_{R,u}^{0}}\sin\phi_{u,0}-\frac{m_{y}\Delta_{R}}{d_{R,u}^{0}}\sin\theta_{u,0}\cos\phi_{u,0}+1\right)$
$\displaystyle=m_{z}\Delta_{R}\sin\phi_{u,0}-m_{y}\Delta_{R}\sin\theta_{u,0}\cos\phi_{u,0}+d_{R,u}^{0}.$
(15)

Therefore, in this case, the array response is characterized only by the
azimuth and vertical AoAs. For ease of exposure, we use
$\omega_{u}=\sin\phi_{u,0}$ and $\varphi_{u}=\sin\theta_{u,0}\cos\phi_{u,0}$
to denote the AoAs. Then, the array response for the far field can be denoted
as

$\bm{a}^{F}_{R}(\omega_{u},\varphi_{u})[m]=\exp\left(-j\frac{2\pi}{\lambda}(m_{z}\Delta_{R}\omega_{u}-m_{y}\Delta_{R}\varphi_{u})\right).$
(16)

where the relationship between index $m$ and $(m_{y},m_{z})$ was given by (3).

According to (13), the Fraunhofer distance would be very large if an extra-
large RIS is utilized. If the RIS is working at high-frequency bands, such as
sub-THz, the Fraunhofer distance also increases due to the smaller wavelength
[38]. Therefore, it is very likely that some UEs are in the near field when
the RIS is deployed for indoor applications such as localization. In this
case, the spherical wavefront of the EM waves should be considered.

####  II-A2 Near field formulation

When UE $u$ is located at the near field of the RIS panel, i.e., the Fresnel
region

$d_{R,u}^{m}\in\left[0.62\frac{L^{3}}{\sqrt{\lambda}},2\frac{L^{2}}{\lambda}\right],$

a good approximation is based on the Fresnel approximation, which corresponds
to the second-order Taylor expansion[39]. So far, the second-order Fresnel
approximation results have been developed for the array with ULA geometry[39,
40, 41], but they cannot be directly applied to the two-dimensional RIS UPA
case. Thus, we derive the two-dimensional Fresnel approximation as

$F(y,z)\approx F(y,z)|_{y=0,z=0}+\frac{\partial F}{\partial y}|_{y=0}\times
y+\frac{\partial F}{\partial z}|_{z=0}\times z\\\
+\frac{1}{2}[y,z]\times\Delta^{2}\bm{F}\times[y,z]^{T},$ (17)

where

$\begin{split}\Delta^{2}\bm{F}&=\begin{bmatrix}\frac{\partial^{2}F}{\partial
y^{2}}|_{y=0,z=0}&\frac{\partial^{2}F}{\partial z\partial y}|_{y=0,z=0}\\\
\frac{\partial^{2}F}{\partial z\partial
y}|_{y=0,z=0}&\frac{\partial^{2}F}{\partial z^{2}}|_{y=0,z=0}\end{bmatrix}\\\
&=\begin{bmatrix}1-(\cos\phi_{u,0}\sin\theta_{u,0})^{2}&\cos\phi_{u,0}\sin\theta_{u,0}\sin\phi_{u,0}\\\
\cos\phi_{u,0}\sin\theta_{u,0}\sin\phi_{u,0}&1-(\sin\phi_{u,0})^{2}\end{bmatrix}.\end{split}$

Then, (17) can be reformulated as

$F(y,z)\approx
1+\omega_{u}z-\varphi_{u}y+\frac{1}{2}(\overline{\omega}_{u}y^{2}+\overline{\varphi}_{u}z^{2}+2\varphi_{u}\omega_{u}yz),$
(18)

where $\overline{\omega}_{u}=1-\omega_{u}^{2}$, and
$\overline{\varphi}_{u}=1-\varphi_{u}^{2}$. Then, the distance of the $m$-th
reflecting element $d_{R,u}^{m}$ can be approximated as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\quad d_{R,u}^{m}$ $\displaystyle\approx
m_{z}\Delta_{R}\omega_{u}-m_{y}\Delta_{R}\varphi_{u}+\frac{1}{2}\overline{\omega}_{u}\frac{(m_{y}\Delta_{R})^{2}}{d_{R,u}^{0}}+\frac{1}{2}\overline{\varphi}_{u}\frac{(m_{z}\Delta_{R})^{2}}{d_{R,u}^{0}}+\varphi_{u}\omega_{u}\frac{m_{y}\Delta_{R}(m_{z}\Delta_{R})}{d_{R,u}^{0}}+d_{R,u}^{0}$
$\displaystyle=m_{z}\Delta_{R}\omega_{u}-m_{y}\Delta_{R}\varphi_{u}-\frac{\left(m_{z}\Delta_{R}\omega_{u}-m_{y}\Delta_{R}\varphi_{u}\right)^{2}}{2d_{R,u}^{0}}+\frac{1}{2d_{R,u}^{0}}\left((m_{z}\Delta_{R})^{2}+(m_{y}\Delta_{R})^{2}\right)+d_{R,u}^{0}$
$\displaystyle\triangleq
J_{m}(\omega_{u},{\varphi_{u}})+Q_{m}(\omega_{u},{\varphi_{u}},d^{0}_{R,u})+d_{R,u}^{0},$
(19)

where
$J_{m}(\omega_{u},{\varphi_{u}})=m_{z}\Delta_{R}\omega_{u}-m_{y}\Delta_{R}\varphi_{u}$,
and

$\begin{split}Q_{m}(\omega_{u},{\varphi_{u}},d^{0}_{R,u})=\frac{1}{2d^{0}_{R,u}}\left((m_{z}\Delta_{R})^{2}+(m_{y}\Delta_{R})^{2}\right)-\frac{1}{2d^{0}_{R,u}}\left(m_{z}\Delta_{R}\omega_{u}-m_{y}\Delta_{R}\varphi_{u}\right)^{2}.\end{split}$

In the near field case, apart from the azimuth and vertical angles, the array
response is also dependent on the distance $d_{R,u}^{0}$, which is the
distance from UE $u$ to the center of RIS. Then, the array response for the
near field UE can be approximated as

$\bm{a}^{N}_{R}(\omega_{u},\varphi_{u},d_{R,u}^{0})[m]=\exp\left(-j\frac{2\pi}{\lambda}(J_{m}(\omega_{u},{\varphi_{u}})+Q_{m}(\omega_{u},{\varphi_{u}},d_{R,u}^{0}))\right).$
(20)

It can be verified that
$\bm{a}_{R}(\omega_{u},\varphi_{u},d_{R,u}^{0})[m]\to\bm{a}_{R}(\omega_{u},\varphi_{u})[m]$
when $d_{R,u}^{0}$ is sufficiently large. In other words, the near field case
reduces to the far-field case. Then, in the near field, the complex channel
gain from UE $u$ to the RIS is denoted as

$\bm{h}_{R,u}=g_{R,u}\bm{a}^{N}_{R}(\omega_{u},\varphi_{u},d_{R,u}^{0}),$ (21)

where the $m$-th element of
$\bm{a}^{N}_{R}(\omega_{u},\varphi_{u},d_{R,u}^{0})$ is given by (20).

(a) Approximation errors of distances.

(b) Cosine of the angles between steering vectors.

Figure 4: Performance of near/far-field approximation.

Overall, the distance $d_{R,u}^{m}$ can be approximated by the near-field
model (19) and the far-field model (II-A1), and the approximation errors are
shown in Fig. (4(a)) with labels “NF” and “FF”, respectively. In Fig. (4(a)),
the Y-axis denotes the summation of approximation errors of $d_{R,u}^{m}$ for
all
$m=-N_{R}^{Z}(N_{R}^{Y}+1)-N_{R}^{Y},\cdots,N_{R}^{Z}(N_{R}^{Y}+1)+N_{R}^{Y}$.
As expected, it is observed that the proposed near-field model is more
accurate than the far field model, as the approximation errors of “NF” are
much less than that of “FF”. Also, when the distance between UE and the RIS
panel increases, the approximation error of the near-field model also
decreases, which indicates that the “NF” model becomes more accurate for
longer transmission distance. Also, the array response given in (1) can be
approximated by the near-field steering vector (20) and the far-field steering
vector (16), respectively, i.e., $\bm{a}_{R}\approx\bm{a}^{N}_{R}$ and
$\bm{a}_{R}\approx\bm{a}^{F}_{R}$. Fig. (4(b)) shows the cosine of the angles
between $\bm{a}_{R}$ and $\bm{a}^{N}_{R}$, and the cosine of the angles
between $\bm{a}^{F}_{R}$ and $\bm{a}_{R}$222The simulation configurations of
Fig. (4(a)) and Fig. (4(b)) as the same as that in Section IV, where the
details of configurations can be found.. It is observed that the near-field
$\bm{a}^{N}_{R}$ is very close to $\bm{a}_{R}$, i.e. the value of cosine is
equal to $1$, while the approximation error of $\bm{a}^{F}_{R}$ is much
larger. In addition, the cosine values of the $\bm{a}^{F}_{R}$ slightly
increases with the distance as the near-field effects decreases.

###  II-B RIS-AP link

It is assumed that the LoS channel between the RIS and the AP is varying
slowly. Due to the spherical nature of the wavefront, the transmission
distances from each element of the RIS to each antenna element of the AP array
determine the phases of the received signals at the AP.

According to [42, 43, 44], a high-rank channel model of the RIS-AP link can be
obtained with a proper deployment. With the spherical wave modeling, the
normalized channel gain of the RIS-AP link is denoted as $\bm{G}_{A,R}$, which
collects the channel elements of the RIS-AP LoS link. The $(m,n)$-th element
of the $\bm{G}_{A,R}$ is given by [42]

$\bm{G}_{A,R}(m,n)=\exp\left(j\frac{2\pi}{\lambda}r_{m,n}\right),$ (22)

where $r_{m,n}$ denotes the path length between the $m$th antenna element of
the AP and the $n$-th reflecting element. Letting $g_{A}$ denote the common
large-scale path loss attenuation, the complex channel gain is then given
by[43]

$\bm{H}_{A,R}=g_{A}\bm{G}_{A,R}.$ (23)

###  II-C Cascaded Channel

Let $\bm{E}_{RIS}=\text{diag}\\{\bm{e}\\}\in\mathbb{C}^{N_{R}\times 1}$ denote
the phase shift matrix of the RIS, where $\bm{e}$ denotes the phase shift
vector of the RIS. Then, the RIS-assisted channel of UE $u$ is expressed as

$\bm{h}_{u}=g_{R,u}\sqrt{P_{u}^{t}}\bm{H}_{A,R}\bm{E}_{RIS}\bm{a}_{R,u}=g_{u}\sqrt{P_{u}^{t}}\bm{G}_{A,R}\bm{E}_{RIS}\bm{a}_{R,u},$
(24)

where $\sqrt{P_{u}^{t}}$ denotes the transmit power of UE $u$, $g_{u}$ is the
cascaded complex channel attenuation defined as $g_{u}=g_{A}g_{R,u}$, and
$\bm{a}_{R,u}=\bm{a}^{N}_{R}(\omega_{u},\varphi_{u},d_{R,u}^{0})$.

In time slot $t$, the transmit signal vector of all the $U$ UEs is denoted as
$\bm{x}_{t}=[x_{1}(t),\cdots,x_{U}(t)]^{T}$ $\in\mathbb{C}^{U\times 1}$. With
the assistance of the RIS, the received signal at the AP is modeled as

$\displaystyle\bm{y}(t)$
$\displaystyle=\sqrt{P_{u}^{t}}\bm{G}_{A,R}\bm{E}_{RIS}(t)\sum_{u=1}^{U}g_{u}\bm{a}_{R,u}x_{u}(t)+\bm{n}(t)$
$\displaystyle=\sqrt{P_{u}^{t}}\bm{G}_{A,R}\bm{E}_{RIS}(t)\bm{A}\bm{x}_{t}+\bm{n}(t),$
(25)

where
$\bm{A}=[g_{1}\bm{a}_{R,1},\cdots,g_{u}\bm{a}_{R,u},\cdots,g_{U}\bm{a}_{R,U}]$,
and $\bm{n}(t)$ denotes the Gaussian noise vector with the distribution of
$\mathcal{CN}(0,\sigma^{2}\bm{I})$.

In the AP-RIS link, as the locations of the RIS and the AP are known as prior,
the distances $\\{r_{m,n}\\}$ between the elements of the AP and that of the
RIS can be readily obtained. First, the unknown channel matrix $\bm{A}$ should
be estimated. Then, to locate UE $u$, the unknown parameters in $\bm{A}$ that
need to be estimated are

  * •

$\omega_{u}$, $\varphi_{u}$ and $d_{R,u}^{0}$, which depends on the location
of UE $u$ in the RIS-UE $u$ link.

  * •

$g_{u}$, which is the complex channel attenuation coefficient.

##  III Joint Channel and Localization Estimation

###  III-A RIS Training Phase Shifts and Pilot Design

First, we need to obtain the unknown channel matrix $\bm{A}$. The simplest
method is to adopt the LS method. However, in the considered scenario, the
number of antenna elements of the AP is less than that of the RIS panel, i.e.,
$N_{A}<N_{R}$, resulting in that $\bm{G}_{A,R}\in\mathbb{C}^{N_{A}\times
N_{R}}$ has the rank of $N_{A}$. Then, the LS estimation cannot be directly
applied to (25). Therefore, to obtain the unique estimation of $\bm{A}$, we
utilize different RIS training phase shift vectors $\bm{e}(t)$ and pilot data
$\bm{x}_{t}$.

Figure 5: Design of RIS phase shift vectors and the pilot data.

Fig. 5 shows the numbers of required RIS training phase shift vectors and
pilot data sets. As shown in Fig. 5, the RIS phase shift vector changes for
$S$ times while the pilot data keeps fixed, and the set of phase shift vectors
is denoted as $\\{\bm{e}_{1},\cdots,\bm{e}_{S}\\}$. The number of RIS training
phase shift vectors is denoted as $S$, and the required number of pilot data
is denoted as $\tau$. Let $\bm{y}_{t}$ denote the composited signal for the
$t$-th pilot data duration, which collects the received signal with $S$
different phase shift vectors as

$\bm{y}_{t}=\bm{G}_{RIS}\bm{A}\bm{x}_{t}+\bm{n}_{t},$ (26)

where the received vector $\bm{y}_{t}$ has the size of $N_{A}S\times 1$,
$\bm{y}_{t}=[\bm{y}(t_{1})^{T},\cdots,\bm{y}(t_{S})^{T}]^{T}$, and
$\bm{n}_{t}=[\bm{n}(t_{1})^{T},\cdots,\bm{n}(t_{S})^{T}]^{T}$. The matrix
$\bm{G}_{RIS}$ with size $N_{A}S\times N_{R}$ collects the composite channels
introduced by $S$ different phase shift vectors, which are
$\bm{G}_{RIS}=[(\bm{G}_{A,R}\text{diag}{(\bm{e}_{1})})^{T},\bm{G}_{A,R}\text{diag}{(\bm{e}_{2})}^{T},\cdots,\\\
\bm{G}_{A,R}\text{diag}{(\bm{e}_{S})}^{T}]^{T}$.

Then, we have the LS estimation for $\bm{A}\bm{x}_{t}$ as

$\widehat{\bm{A}\bm{x}}_{t}=(\bm{G}_{RIS}^{H}\bm{G}_{RIS})^{-1}\bm{G}_{RIS}^{H}\bm{y}_{t}=\bm{H}^{{\dagger}}_{RIS}\bm{y}_{t}.$
(27)

Remark 1: It is observed that the rank of cascaded channel $\bm{G}_{RIS}$ can
be increased by applying different RIS phase shifts, so that the Moore-Penrose
pseudoinverse of $\bm{G}_{RIS}$ denoted as $\bm{H}^{{\dagger}}_{RIS}$ can be
obtained. As demonstrated in [45], the Discrete Fourier Transform (DFT) RIS
phase training scheme can achieve the nearly optimal performance in terms of
the mean squared error (MSE). Thus, the training phase shift matrix
${\bm{E}_{S}=[\bm{e}_{1},\bm{e}_{2},\cdots,\bm{e}_{S}]}$ can be set as the
first $S$ columns of the $N_{R}\times N_{R}$ DFT matrix.

For transmit pilot design, we assume that transmit signal from UE $u$ is
denoted as $x_{u}$ following the distribution of $\mathbb{E}[|x_{u}|^{2}]=1$
and $\mathbb{E}[x_{u}x_{k}]=0$, for all $k\neq u$. Consequently, the array
covariance matrix can be estimated as

$\tilde{\bm{R}}=\mathbb{E}[\bm{H}^{{\dagger}}_{IRS}\bm{y}_{t}\bm{y}_{t}^{H}(\bm{H}^{{\dagger}}_{IRS})^{H}]=\bm{A}\bm{A}^{H}+\sigma^{2}(\bm{G}_{RIS}^{H}\bm{G}_{RIS})^{-1}.$
(28)

By collecting the received signal $\bm{y}_{t}$ with different $\tau$ pilot
durations, the estimated array covariance matrix $\hat{\bm{R}}$ of
$\tilde{\bm{R}}$ is given by

$\hat{\bm{R}}=\frac{1}{\tau}\sum_{t=1}^{\tau}\widehat{\bm{A}\bm{x}}_{t}(\widehat{\bm{A}\bm{x}}_{t})^{H},$
(29)

where the total number of snapshots is $\tau$.

Remark 2: In order to ensure that the matrix $\bm{G}_{RIS}$ is full column
rank, the number of training phase shift vectors should satisfy
$S>\frac{N_{R}}{N_{A}}$. Moreover, the output noise of the LS estimator is
amplified by $(\bm{G}_{RIS}^{H}\bm{G}_{RIS})^{-1}$ as shown in (28).
Theoretically, the Hermitian matrix $\bm{G}_{H}=\bm{G}_{RIS}^{H}\bm{G}_{RIS}$
may have rank deficiency, i.e. the minimum singular values of $\bm{G}_{H}$ is
sufficiently small so that $\rho_{min}(\bm{G}_{H})\to 0$. In this case, the
noise will be significantly amplified. Therefore, the singular values of the
constructed Hermitian matrix should be large enough. Thus, in order to
suppress noise, the number of RIS training phase shift vectors should be
sufficiently large, and its impact on the performance will be give in Section
IV.

###  III-B Estimation of AoAs

Without noise, we have the following explicit expression

$\displaystyle\bm{R}\\!\\!=\\!\\!\bm{A}\bm{A}^{H}\\!\\!\\!=\sum_{u=1}^{U}|g_{R,u}|^{2}\bm{a}_{R,u}\bm{a}^{H}_{R,u}.$
(30)

Then, the $(p,q)$-th element of covariance $\bm{R}$ is given by

$\bm{R}(p,q)=\sum_{u=1}^{U}|g_{R,u}|^{2}\exp\left(-j\frac{2\pi}{\lambda}(J_{p,q}+Q_{p,q})\right),$
(31)

where
$J_{p,q}=J_{p}(\omega_{u},{\varphi_{u}})-J_{q}(\omega_{u},{\varphi_{u}})$,
$Q_{p,q}=Q_{p}(\omega_{u},{\varphi_{u}},d_{R,u}^{0})-Q_{q}(\omega_{u},{\varphi_{u}},d_{R,u}^{0})$.
To be specific,

$\displaystyle Q_{p,q}$
$\displaystyle=\frac{1}{2d_{R,u}^{0}}\left((p_{z}\Delta_{R})^{2}+(p_{y}\Delta_{R})^{2}\right)-\frac{1}{2d_{R,u}^{0}}\left(p_{z}\Delta_{R}\omega_{u}-p_{y}\Delta_{R}\varphi_{u}\right)^{2}$
$\displaystyle-\frac{1}{2d_{R,u}^{0}}\left((q_{z}\Delta_{R})^{2}+(q_{y}\Delta_{R})^{2}\right)+\frac{1}{2d_{R,u}^{0}}\left(q_{z}\Delta_{R}\omega_{u}-q_{y}\Delta_{R}\varphi_{u}\right)^{2}.$
$\displaystyle J_{p,q}$
$\displaystyle=(p_{z}-q_{z})\Delta_{R}\omega_{u}-(p_{y}-q_{y})\Delta_{R}\varphi_{u},$

where indices are given by
$p(p_{z},p_{y})=N^{Z}_{R}(2N^{Y}_{R}+1)+N^{Y}_{R}+p_{z}(2N^{Y}_{R}+1)+p_{y}$,
$q(q_{z},q_{y})=N^{Z}_{R}(2N^{Y}_{R}+1)+N^{Y}_{R}+q_{z}(2N^{Y}_{R}+1)+q_{y}$,
$p_{y},q_{y}=-N^{Y}_{R},\cdots,N^{Y}_{R},p_{z},q_{z}=-N^{Z}_{R},\cdots,N^{Z}_{R}$.
Then, it is found that for all $p_{y}=-q_{y}$, and $p_{z}=-q_{z}$, one obtains

$J_{p,q}=2p_{z}\Delta_{R}\omega_{u}-2p_{y}\Delta_{R}\varphi_{u},\qquad
Q_{p,q}=0.$

In other words, the quadratic term is eliminated for the element of
$\bm{R}(p(p_{z},p_{y}),q(-p_{z},-p_{y}))$ in the covariance matrix $\bm{R}$,
which are

$\bm{R}(p(p_{z},p_{y}),q(-p_{z},-p_{y}))=\sum_{u=1}^{U}\\!|g_{R,u}|^{2}\\!\exp\\!\left(\\!-j\frac{2\pi}{\lambda}(2J_{p}(\omega_{u},{\varphi_{u}}))\\!\right).$

We provide the following approach to separate the distance and the AoAs in the
covariance matrix. First, we define steering vectors as

$\displaystyle\bm{v}_{u}(\omega_{u})[z]=\exp\left(-j\frac{2\pi}{\lambda}(2z\Delta_{R}\omega_{u})\right),z=0,\cdots,N^{Z}_{R},$
$\displaystyle\bm{s}_{u}(\varphi_{u})[y]=\exp\left(-j\frac{2\pi}{\lambda}(-2y\Delta_{R}\varphi_{u})\right),y=0,\cdots,N^{Y}_{R},$
$\displaystyle\bm{b}_{u}(\omega_{u},\varphi_{u})=\bm{v}_{u}(\omega_{u})\otimes\bm{s}_{u}(\varphi_{u})=[\bm{v}_{u}(\omega_{u})[0]\bm{s}^{T}_{u}(\varphi_{u}),\cdots,\bm{v}_{u}(\omega_{u})[N^{Z}_{R}]\bm{s}^{T}_{u}(\varphi_{u})]^{T}.$

Then, we define the correlation matrix $\bm{R}_{b}$ as

$\bm{R}_{b}=\sum_{u=1}^{U}|g_{R,u}|^{2}\bm{b}_{u}(\omega_{u},\varphi_{u})\bm{b}^{H}_{u}(\omega_{u},\varphi_{u}),$
(32)

Defining the following function

$e(z,y)=\sum_{u=1}^{U}|g_{R,u}|^{2}\exp\left(-j\frac{4\pi}{\lambda}(z\Delta_{R}\omega_{u}-y\Delta_{R}\varphi_{u})\right),$
(33)

we can represent the covariance matrix $\bm{R}_{b}$ as (34). Note that

$\bm{R}(p(z,y),q(-z,-y))=e(z,y)=(e(-z,-y))^{*}.$

For instance, the first element of the matrix $\bm{R}_{b}$, i.e., $e(0,0)$ is
the same as $\bm{R}(p(0,0),q(0,0))$, where
$p(0,0)=N^{Z}_{R}(2N^{Y}_{R}+1)+N^{Y}_{R}=\frac{N_{R}-1}{2},q(0,0)=N^{Z}_{R}(2N^{Y}_{R}+1)+N^{Y}_{R}=\frac{N_{R}-1}{2}$.

$\bm{R}_{b}=\begin{bmatrix}e(0,0)&\\!e(0,-1)&\cdots&e(0,-N^{Y}_{R})\\!&\\!\\!\\!\\!\\!\\!e(-1,0)&\\!\\!\cdots\\!\\!&e(-N^{Z}_{R},-N^{Y}_{R})\\\
e(0,1)&\\!e(0,0)&\cdots&e(0,-N^{Y}_{R}+1)\\!&\\!\\!\\!\\!\\!\\!e(-1,1)&\\!\\!\cdots\\!\\!&e(-N^{Z}_{R},-N^{Y}_{R}+1)\\\
\vdots&\\!\vdots&\ddots&\vdots\\!&\\!\\!\\!\\!\\!\\!\vdots\\!&\\!\\!\ddots\\!\\!&\vdots\\\
e(0,N^{Y}_{R})&\\!e(0,N^{Y}_{R}-1)&\cdots\\!&\\!\\!e(0,0)\\!&\\!\\!\\!\\!\\!\\!e(-1,N^{Y}_{R})\\!\\!&\cdots&e(-N^{Z}_{R},0)\\\
\vdots&\vdots&\\!\ddots&\vdots&\\!\\!\\!\\!\\!\\!\vdots\\!&\\!\\!\ddots\\!\\!&\vdots\\\
e(N^{Z}_{R},N^{Y}_{R})&\\!e(N^{Z}_{R},N^{Y}_{R}-1)\\!&\\!\\!\cdots\\!&\\!\\!e(N^{Z}_{R},0)\\!&\\!\\!\\!\\!\\!\\!e(N^{Z}_{R}-1,N^{Y}_{R})&\\!\\!\cdots\\!\\!&e(0,0)\\\
\end{bmatrix}.$ (34)

$\bm{S}_{i}\\!=\\!\\!\\!\begin{bmatrix}\bm{R}(N_{i},\overline{N}_{i})&\bm{R}(N_{i}-1,\overline{N}_{i}+1)&\\!\\!\\!\cdots\\!\\!\\!&\\!\\!\\!\bm{R}(N_{i}-N^{Y}_{R},\overline{N}_{i}+N^{Y}_{R})\\\
\bm{R}(N_{i}+1,\overline{N}_{i}-1)&\bm{R}(N_{i},\overline{N}_{i})&\\!\\!\\!\cdots\\!\\!\\!&\\!\\!\\!\bm{R}(N_{i}-N^{Y}_{R}-1,\overline{N}_{i}+N^{Y}_{R}+1)\\\
\vdots&\\!\\!\\!\vdots\\!\\!\\!&\\!\\!\\!\ddots&\vdots\\\
\bm{R}(N_{i}+N^{Y}_{R},\overline{N}_{i}-N^{Y}_{R})&\bm{R}(N_{i}+N^{Y}_{R}-1,\overline{N}_{i}-N^{Y}_{R}+1)&\\!\\!\\!\cdots\\!\\!\\!&\\!\\!\\!\bm{R}(N_{i},\overline{N}_{i})\\\
\end{bmatrix}.$ (35)

Let $N_{b}$ represent the length of vector $\bm{b}_{u}$, and
$N_{b}=(N^{Z}_{R}+1)\times(N^{Y}_{R}+1)$. The index of the center reflecting
element is denoted as $N_{0}$, and $N_{0}=\frac{N_{R}-1}{2}$. Then, the index
functions can be represented as $p(z,y)=N_{0}+z(2N^{Y}_{R}+1)+y$, and
$q(-z,-y)=N_{0}-z(2N^{Y}_{R}+1)-y$. Letting $N_{i}=N_{0}+i(2N^{Y}_{R}+1)$ and
$\overline{N}_{i}=N_{0}-i(2N^{Y}_{R}+1)$, the element matrix $\bm{S}_{i}$ with
size $(N^{Y}_{R}+1)\times(N^{Y}_{R}+1)$ is given by (35).

Utilizing $\bm{T}_{i}$, we construct the down-sampled Toeplitz matrix with
size $N_{b}\times N_{b}$ as

$\bm{T}=\begin{bmatrix}\bm{S}_{0}&\bm{S}_{-1}&\cdots&\bm{S}_{-N^{Z}_{R}}\\\
\bm{S}_{1}&\bm{S}_{0}&\cdots&\bm{S}_{-N^{Z}_{R}+1}\\\
\vdots&\vdots&\ddots&\vdots\\\
\bm{S}_{N^{Z}_{R}}&\bm{S}_{N^{Z}_{R}-1}&\cdots&\bm{S}_{0}\\\ \end{bmatrix}.$
(36)

We can derive that

$\bm{T}=\sum_{u=1}^{U}|g_{R,u}|^{2}\bm{b}_{u}(\omega_{u},\varphi_{u})\bm{b}^{H}_{u}(\omega_{u},\varphi_{u}).$
(37)

It is observed that the sampled correlation matrix $\bm{T}$ only depends on
the azimuth and vertical angles, and which is the same as the far-field case.

Remark 3: Through this method, the distance and the AoAs can be decoupled by
the sampled correlation matrix $\bm{T}$. Then, the conventional far-field
angle estimation approaches can be leveraged in this case, such as the two-
dimensional MUSIC-like spectrum peak searching[46]. However, note that the
MUSIC algorithm requires eigendecomposition. In the following, we propose a
computationally efficient subspace-based method for the AoA estimation.

First, we define matrix $\bm{B}=[g_{R,1}\bm{b}_{1},\cdots,g_{R,U}\bm{b}_{U}]$,
and divide the $N_{b}\times U$ matrix into two parts as

$\bm{B}\triangleq\left[\begin{array}[]{c}\bm{B}_{1}\in\mathbb{C}^{\small{U}\times\small{U}}\\\
\bm{B}_{2}\in\mathbb{C}^{(\small{N_{R}-U})\times\small{U}}\end{array}\right]$
(38)

where matrix $\bm{B}_{1}$ contains the first $U$ rows of $\bm{B}$ and
$\bm{B}_{2}$ contains the remaining $N_{b}-U$ rows. Note that $\bm{B}_{1}$ is
a Vandermonde matrix. Assuming that the $U$ UEs are from distinct directions,
$\bm{B}_{1}$ has the full rank of $U$, and the rows of $\bm{B}_{2}$ can be
expressed as a linear combination of linearly independent rows of
$\bm{B}_{1}$[47, 48]. Equivalently, there is a $U\times(N_{b}-U)$ linear
operator $\bm{P}_{1}$ between $\bm{B}_{1}$ and $\bm{B}_{2}$ as

$\bm{P}_{1}^{H}\bm{B}_{1}=\bm{B}_{2},$ (39)

and the linear operator $\bm{P}_{1}$ can be calculated as

$\bm{P}_{1}=(\bm{B}_{1}^{H})^{-1}\bm{B}_{2}^{H}.$ (40)

Then, we utilize the correlation matrix to obtained linear operator
$\bm{P}_{1}$. The correlation matrix $\bm{T}$ can be represented as

$\bm{T}=\bm{B}\bm{B}^{H}=\begin{bmatrix}\bm{B}_{1}\bm{B}^{H}\\\
\bm{B}_{2}\bm{B}^{H}\end{bmatrix}=\begin{bmatrix}\bm{T}_{1}\in\mathbb{C}^{U\times
N_{R}}\\\ \bm{T}_{2}\in\mathbb{C}^{(N_{R}-U)\times
N_{R}}\end{bmatrix}\begin{matrix}\end{matrix}$ (41)

where $\bm{T}_{1}$ consists of the first $U$ rows and $\bm{T}_{2}$ consists of
the remaining $(N_{b}-U)$ rows. Then, the relationship between the submatrices
of $\bm{T}_{1}$ and $\bm{T}_{2}$ is

$\bm{P}_{1}^{H}\bm{T}_{1}=\bm{T}_{2}.$ (42)

Then, the linear operator $\bm{P}_{1}$ can be found from $\bm{T}_{1}$ and
$\bm{T}_{2}$ as

$\bm{P}_{1}=(\bm{T}_{1}\bm{T}_{1}^{H})^{-1}\bm{T}_{1}\bm{T}_{2}^{H}.$ (43)

With the obtained linear operator $\bm{P}_{1}$, we have

$\bm{Q}^{H}\bm{B}=\bm{0}_{(N_{b}-U)\times N_{b}},$ (44)

where $\bm{Q}^{H}=[\bm{P}_{1}^{H},-\bm{I}_{N_{b}-U}]$. As a result, the
columns of $\bm{Q}$ are in the null space of $\bm{B}$, i.e.
$\mathcal{N}(\bm{B})$ and the orthogonal projector onto this subspace
$\mathcal{N}(\bm{B})$ is given by

$\bm{\Pi}_{Q}\triangleq\bm{Q}(\bm{Q}^{H}\bm{Q})^{-1}\bm{Q}^{H}.$ (45)

Then, we can infer that $\bm{\Pi}_{Q}\bm{b}_{u}=\bm{0},\forall u=1,\cdots,U$.
Define
$\bm{v}(\omega)[z]=\exp(-j\frac{2\pi}{\lambda}(2z\Delta_{R}\omega)),z=0,\cdots,N^{Z}_{R}$,
and
$\bm{s}(\varphi)[y]=\exp\left(-j\frac{2\pi}{\lambda}(-2y\Delta_{R}\varphi)\right),y=0,\cdots,N^{Y}_{R}$.
Then, the angles can be estimated in a manner similar to the MUSIC method by
minimizing the following cost function as

$f(\omega,\varphi)=\bm{b}^{H}(\omega,\varphi)\bm{\Pi}_{Q}\bm{b}(\omega,\varphi),$
(46)

where $\bm{b}(\omega,\varphi)=\bm{v}(\omega)\otimes\bm{s}(\varphi)$.

Directly searching the two-dimensional angle space requires a significant
computation cost. In the following, we provide a low-complexity algorithm to
solve the angle estimation problem. Then, we formulate the AoA estimation
problem as

$\displaystyle\underset{\omega,\varphi}{\text{min}}$ $\displaystyle\quad
f(\omega,\varphi)=(\bm{v}(\omega)\otimes\bm{s}(\varphi))^{H}\bm{\Pi}_{Q}(\bm{v}(\omega)\otimes\bm{s}(\varphi))$
(47a) s.t. $\displaystyle\quad\bm{e}^{H}\bm{s}(\varphi)=1,$ (47b)

where $\bm{e}=[1,0,\cdots,0]^{T}$, and the constraint (47b) is introduced to
avoid the trivial solution $\bm{s}(\varphi)=\bm{0}$. Then, the objective
function is reformulated as

$\displaystyle f(\omega,\varphi)$
$\displaystyle=\bm{s}(\varphi)^{H}[\bm{v}(\omega)^{H}\otimes\bm{I}_{N^{Y}_{R}+1}]\bm{\Pi}_{Q}[\bm{v}(\omega)\otimes\bm{I}_{N^{Y}_{R}+1}]\bm{s}(\varphi)$
(48) $\displaystyle=\bm{s}(\varphi)^{H}\bm{\Theta}(\omega)\bm{s}(\varphi),$
(49)

where
$\bm{\Theta}(\omega)=[\bm{v}(\omega)^{H}\otimes\bm{I}_{N^{Y}_{R}+1}]\bm{\Pi}_{Q}[\bm{v}(\omega)\otimes\bm{I}_{N^{Y}_{R}+1}]$.

We construct the following Lagrange function

$\mathcal{L}=\bm{s}(\varphi)^{H}\bm{\Theta}(\omega)\bm{s}(\varphi)-\lambda(\bm{e}^{H}\bm{s}(\varphi)-1),$
(50)

where $\lambda$ is the introduced Lagrange multiplier. Then, one obtains

$\frac{\partial\mathcal{L}}{\partial\bm{s}(\varphi)}=2\bm{\Theta}(\omega)\bm{s}(\varphi)-\lambda\bm{e}.$
(51)

We can infer that $\bm{s}(\varphi)=\zeta\bm{\Theta}^{-1}(\omega)\bm{e}$, where
$\zeta$ is a non-zero constant. As $\bm{e}^{H}\bm{s}(\varphi)=1$, we have
$\zeta=(\bm{e}^{H}\bm{\Theta}^{-1}(\omega)\bm{e})^{-1}$. Then, one obtains

$\bm{s}(\varphi)=\frac{\bm{\Theta}^{-1}(\omega)\bm{e}}{\bm{e}^{H}\bm{\Theta}^{-1}(\omega)\bm{e}}.$
(52)

Then, substituting (52) into the objective function (48), Problem (47) is
transformed into

$\underset{\omega}{\text{min}}\quad
f(\omega)=\frac{1}{\bm{e}^{H}\bm{\Theta}^{-1}(\omega)\bm{e}}.$ (53)

Then, the optimal solution to (53) is

$\omega^{*}=\arg\max\bm{e}^{H}\bm{\Theta}^{-1}(\omega)\bm{e}.$ (54)

Remark 4: In practice, two-layer grids can be used for searching the optimal
solution $\omega^{*}$. To be specific, the grid with large step-size is used
for the preliminary coarse search. Then, based on the preliminary coarse
results, the finer grid is then applied to obtain the fine results.
Furthermore, due to the presence of noise, the energy peak will leakage to
adjacent grid points. Therefore, to enhance the position accuracy, it is
suggested to take several largest values as estimation candidates, and then
use the clustering method, such as the K-mean method[49], to obtain a more
accurate estimation.

Recall that $\omega_{u}=\sin\phi_{u}$ and
$\varphi_{u}=\sin\theta_{u}\cos\phi_{u}$, where $\phi_{u}$ and $\theta_{u}$
represent the vertical and azimuth AoAs of UE $u$, respectively. Then, we can
search the interval $\omega_{u}\in[-1,1]$, and obtain $U$ largest peaks of the
$(1,1)$-th element of $\bm{\Theta}^{-1}(\omega)$. These $U$ peaks correspond
to the vertical AoAs of UEs, denoted as
$[\hat{\omega}_{1},\hat{\omega}_{2},\cdots,\hat{\omega}_{U}]$. Then, for the
composite angle $\varphi_{u}$, according to (52), we can obtain the $U$
vectors
$[\bm{\hat{s}}(\varphi_{1}),\bm{\hat{s}}(\varphi_{2}),\cdots,\bm{\hat{s}}(\varphi_{U})]$.

Define

$\bm{\hat{q}}_{u}=\angle\bm{\hat{s}}(\varphi_{u})=[0,\frac{4\pi}{\lambda}\Delta_{R}\varphi_{u},\cdots,\frac{4\pi}{\lambda}\Delta_{R}N_{R}^{Y}\varphi_{u}]^{T},$
(55)

and vector
$\bm{p}=[0,\frac{4\pi}{\lambda}\Delta_{R},\cdots,\frac{4\pi}{\lambda}\Delta_{R}N_{R}^{Y}]^{T}$.
Then, we utilize the LS estimator to obtain the estimation of $\varphi_{u}$ as

$\hat{\varphi}_{u}^{LS}=(\bm{p}^{T}\bm{p})^{-1}\bm{p}^{T}\bm{\hat{q}}_{u}.$
(56)

###  III-C Estimation of Distances

With the estimated AoAs $\\{\hat{\varphi}_{u}\\}$ and
$\\{\hat{\omega}_{u}\\}$, we reformulate the array response as

$\bm{a}_{R}(\hat{\omega}_{u},\hat{\varphi}_{u},d_{R,u}^{0})=\text{diag}\\{\bm{p}_{u}(\hat{\omega}_{u},\hat{\varphi}_{u})\\}\bm{q}_{u}(\hat{\omega}_{u},\hat{\varphi}_{u},d_{R,u}^{0}),$
(57)

where
$\bm{p}_{u}(\hat{\omega}_{u},\hat{\varphi}_{u})=\left[\exp\left(-j\frac{2\pi}{\lambda}J_{-\frac{N_{R}-1}{2}}(\hat{\omega}_{u},\hat{\varphi}_{u})\right),\cdots,\right.$
$\left.\exp\left(-j\frac{2\pi}{\lambda}J_{\frac{N_{R}-1}{2}}(\hat{\omega}_{u},\hat{\varphi}_{u})\right)\right]^{T}$.
Also,

$\bm{q}_{u}(\hat{\omega}_{u},\hat{\varphi}_{u},d_{R,u}^{0})=\left[\exp\left(-j\frac{2\pi}{\lambda
d_{R,u}^{0}}q_{-\frac{N_{R}-1}{2}}\right),\cdots,\exp\left(-j\frac{2\pi}{\lambda
d_{R,u}^{0}}q_{\frac{N_{R}-1}{2}}\right)\right]$

with the
$q_{m}=\left((m_{z}\Delta_{R})^{2}+(m_{y}\Delta_{R})^{2}-\left(m_{z}\Delta_{R}\hat{\omega}_{u}-m_{y}\Delta_{R}\hat{\varphi}_{u}\right)^{2}\right)$,
and index $m$ was given by (3).

Similar to (38) and (41), we divide the matrix $\bm{A}$ and its covariance
matrix $\bm{R}$ into two partitions as
$\bm{A}^{H}=[\bm{A}^{H}_{1},\bm{A}^{H}_{2}]$ and
$\bm{R}^{H}=[\bm{R}^{H}_{1},\bm{R}^{H}_{2}]$ respectively, where
$\bm{A}_{1},\bm{R}_{1}\in\mathbb{C}^{U\times N_{R}}$ collects the first $U$
rows and $\bm{A}_{2},\bm{R}_{2}\in\mathbb{C}^{(N_{R}-U)\times N_{R}}$ contains
the remaining $(N_{R}-U)$ rows. As $\bm{A}_{1}$ and $\bm{R}^{H}_{1}$ are of
full rank $U$, the rows of $\bm{R}_{2}(\bm{A}_{2})$ can be expressed as a
linear combination of linear independent rows of $\bm{R}_{1}(\bm{A}_{1})$ with
the $U\times(N_{R}-U)$ linear operator $\bm{P}_{2}$. Similar to (40)-(42), one
obtains

$\bm{P}_{2}^{H}\bm{A}_{1}=\bm{A}_{2},\bm{P}_{2}^{H}\bm{R}_{1}=\bm{R}_{2}.$
(58)

Then, $\bm{P}_{2}$ is given by

$\bm{P}_{2}=(\bm{R}_{1}\bm{R}_{1}^{H})^{-1}\bm{R}_{1}\bm{R}_{2}^{H}.$ (59)

Defining $\bm{Q}^{H}_{p}=[\bm{P}_{2}^{H},-\bm{I}_{N_{R}-U}]$, the orthogonal
projector of $\bm{A}$ is given by

$\displaystyle\bm{Q}^{H}_{p}\bm{A}=\bm{0},\bm{\Pi}_{p}\triangleq\bm{Q}_{p}(\bm{Q}_{p}^{H}\bm{Q}_{p})^{-1}\bm{Q}_{p}^{H}.$
(60)

Then, similar to (46), we have the following distance estimation problem

$\displaystyle{d^{0}_{R,u}}^{*}\\!\\!\\!$
$\displaystyle=\underset{d_{R,u}^{0}}{\text{min}}f(d_{R,u}^{0})$ (61)
$\displaystyle=\underset{d_{R,u}^{0}}{\text{min}}(\text{diag}\\{\bm{p}_{u}\\}\bm{q}_{u}(d_{R,u}^{0}))^{H}\bm{\Pi}_{p}(\text{diag}\\{\bm{p}_{u}\\}\bm{q}_{u}(d_{R,u}^{0})).$

Then, the distance can be estimated by conducting one-dimensional search in
the interval $[d_{min},d_{max}]$. The maximal $U$ largest peaks of the
searching results are readily obtained as the estimated distances.

###  III-D Estimation of Locations and Channel Gains

With the AoAs and distances, the position of UE $u$ denoted as
$(x_{u},y_{u},z_{u})$ can be obtained according to its geometric relationship
with RIS position $(x_{R},y_{R},z_{R})$, given by

$\displaystyle
x_{u}=\sqrt{(1-\hat{\omega}_{u}^{2})\hat{d}^{2}_{R,u}-\varphi^{2}_{u}\hat{d}^{2}_{R,u}}.$
(62) $\displaystyle y_{u}=\varphi_{u}\hat{d}_{R,u}$ (63) $\displaystyle
z_{u}=z_{R}-\hat{\omega}_{u}\hat{d}_{R,u}$ (64)

It is worth pointing out that , in the far-field model, with one RIS panel,
only the AoAs can be obtained, but the 3D position of the UE cannot be
obtained due to the lack of dimensions. To obtain the channel fading
coefficients $\bm{g}=[g_{1},\cdots,g_{U}]$, we utilize the received signal
$\bm{y}_{t}$ in the $t$-th pilot data duration, which collects the received
signal with these $L$ different phase shift vectors as

$\bm{y}_{t}=\bm{G}_{RIS}\bm{A}\bm{x}_{t}+\bm{n}_{t}\\\
=\sum_{u=1}^{U}\bm{G}_{RIS}\bm{a}_{R}({\omega}_{u},{\varphi}_{u},{d}_{R,u})\bm{x}_{t}(u)+\bm{n}_{t},$
(65)

where the received vector $\bm{y}_{t}$ has the size of $N_{A}L\times 1$.

Without loss of generality, we sort the UEs according to the ascending order
of their distances, i.e., $d_{1}\leq\cdots\leq d_{U}$. The reason is that a
far UE normally experiences a severe channel fading gain. Then, the received
residual signal is initialized as $\bm{y}^{Res}_{0}=\bm{y}_{t}$, and the
residual signal for UE $u$ is calculated as

$\bm{y}^{Res}_{u}=\bm{y}_{t}-\sum_{k=1}^{u-1}\hat{g}_{k,t}\bm{G}_{RIS}\bm{a}_{R}(\hat{\omega}_{k},\hat{\varphi}_{k},\hat{d}_{R,k})\bm{x}_{t}(k).$
(66)

According to the OMP algorithm, the complex channel gain can be formulated as
the solution to

$\arg\underset{g_{u,t}}{\text{min}}||\bm{y}^{Res}_{u}-\hat{g}_{u,t}\bm{G}_{RIS}\bm{a}_{R}(\hat{\omega}_{u},\hat{\varphi}_{u},\hat{d}_{R,u})\bm{x}_{t}(u)||^{2}.$
(67)

Then, the estimated channel fading parameter is

$\hat{g}_{u,t}=\frac{\bm{a}_{R,u}^{H}\bm{y}^{Res}_{u}}{\sqrt{P_{u}^{t}}\bm{a}_{R,u}^{H}\bm{a}_{R,u}},$
(68)

where
$\bm{a}_{R,u}=\bm{G}_{RIS}\bm{a}_{R}(\hat{\omega}_{u},\hat{\varphi}_{u},\hat{d}_{R,u})$.
After we have collected $\tau$ different slot pilot data, the estimated
channel fading gain for each time slot should be averaged as

$\hat{g}_{u}=\frac{1}{\tau}\sum_{t=1}^{\tau}\hat{g}_{u,t}.$ (69)

Overall, the complex channel gain of the RIS-UE $u$ link can be obtained as

$\hat{\bm{h}}_{u}=\hat{g}_{u}\bm{G}_{A,R}\bm{E}_{RIS}\bm{G}_{RIS}\bm{a}_{R}(\hat{\omega}_{u},\hat{\varphi}_{u},\hat{d}_{R,u}).$
(70)

###  III-E Complexity Analysis

Algorithm 1 Near Field Joint Channel Estimation and Localization (NF-JCEL)
Algorithm

0: The received signal $\bm{y}_{t}\in N_{A}L\times 1$ ; pilot data $x_{t}$,
$t=1,\cdots,\tau$.

0: The AoAs of the RIS-UE links $\omega_{u}$, $\varphi_{u}$; The distance
$d_{R,u}^{0}$ of the $u$-th RIS-UE link; The channel fading coefficients
$g_{u}$; The UE locations$(x_{u},y_{u},z_{u})$; Complex channel gains
$\hat{\bm{h}}_{u}$.

1: Estimate array covariance matrix $\hat{\bm{R}}$ according to (29) ;

2: Construct matrix $\bm{T}$ according to (36);

3: Calculate AoAs $\omega_{u}$ and $\varphi_{u}$ by solving Problem (47),
where the optimal $\omega_{u}^{*}$ is obtained according to (54) and the
estimated $\varphi_{u}$ is obtained according to (56);

4: Search optimal distance ${d_{R,u}^{0}}^{*}$ by solving Problem (61);

5: Calculate locations$(x_{u},y_{u},z_{u})$ according to (64) -(62);

6: Obtain the channel fading coefficients $\bm{g}=[g_{1},\cdots,g_{U}]$ by
using the OMP method;

In summary, a NF-JCEL algorithm is proposed to estimate the UE locations and
the complex channel coefficients. The detailed algorithm is presented in
Algorithm 1. The proposed NF-JCEL algorithm solves three sub-problems
sequentially: 1) the estimation of array covariance matrix based on the LS
estimation; 2) the estimation of locations, including the AoAs and distance;
3) the channel fading coefficients estimation. Therefore, the complexity of
the NF-JCEL algorithm can be evaluated from these three parts, respectively.

The complexity of step 1 mainly lies in the matrix calculation. The matrix
$\bm{G}_{RIS}$ has size $N_{A}L\times N_{R}$, and the complexity of product
$\bm{G}_{RIS}^{H}\bm{G}_{RIS}$ is $\mathcal{O}(N_{R}^{2}N_{A}L)$. Its
inversion requires the complexity of $\mathcal{O}(N_{R}^{3})$. Thus, the
calculation of $\bm{H}^{{\dagger}}_{RIS}$ has the total complexity of
$\mathcal{O}(N_{R}^{3}+2N_{R}^{2}N_{A}L)$. Then, the complexity of LS
estimation $\widehat{\bm{A}\bm{x}}_{t}$ is given by
$\mathcal{O}(N_{R}^{3}+2N_{R}^{2}N_{A}L+N_{R}N_{A}L)$. Finally, the complexity
of array covariance matrix estimation is
$\mathcal{O}(\tau(N_{R}^{3}+2N_{R}^{2}N_{A}L+N_{R}N_{A}L+N_{R}))=\mathcal{O}(\tau(N_{R}^{3}+2N_{R}^{2}N_{A}L))$.

The complexity of UE localization mainly lies in the orthogonal projectors and
optimization objectives. Similarly, according to (43), $\bm{P}_{1}$ is
calculated with complexity $\mathcal{O}(U^{3}+2U^{2}N_{b}+(N_{b}-U)UN_{b})$,
and the complexity of calculating $\bm{\Pi}_{Q}$ is
$\mathcal{O}((N_{b}-U)^{3}+2(N_{b}-U)^{2}N_{b})$. Then, the complexity of
calculating $\bm{\Theta}^{-1}(\omega)$ is
$\mathcal{O}((N_{R}^{Y}+1)N_{b}^{2}+(N_{R}^{Y}+1)^{2}N_{b}+(N_{R}^{Y}+1)^{3}+(N_{R}^{Y}+1)^{2}+N_{R}^{Y})=\mathcal{O}((N_{R}^{Y}+1)N_{b}^{2}+(N_{R}^{Y}+1)^{2}N_{b})$.
Consequently, the total complexity of estimating $\omega_{u}$ is
$\mathcal{O}(N_{\omega}((N_{R}^{Y}+1)^{2}N_{b}^{2}+(N_{R}^{Y}+1)^{3}N_{b}))$,
where $N_{\omega}$ denotes the number of searching grids for solving Problem
(54). In addition, the complexity of estimating ${\phi}_{u}$ is
$\mathcal{O}(2(N_{Y}+1))$. Furthermore, the complexity of solving Problem (61)
is given by $\mathcal{O}(N_{d}(2N_{R}^{2}+N_{R}))$, where $N_{d}$ denotes the
number of searching grids. Overall, the complexity of UE localization is
$\mathcal{O}(UN_{\omega}((N_{R}^{Y}+1)^{2}N_{b}^{2}+(N_{R}^{Y}+1)^{3}N_{b})+2(N_{Y}+1)+N_{d}(2N_{R}^{2}+N_{R}))$.
Finally, the complexity of channel fading estimation is $\mathcal{O}(2\tau
N_{A}U)$. Then, the overall complexity of NF-JCEL algorithm is given by
$\mathcal{O}(\tau(N_{R}^{3}+2N_{R}^{2}N_{A}L)+UN_{\omega}((N_{R}^{Y}+1)^{2}N_{b}^{2}+(N_{R}^{Y}+1)^{3}N_{b})+2(N_{Y}+1)+N_{d}(2N_{R}^{2}+N_{R})+2\tau
N_{A}U)$.

##  IV Simulation results

In this section, some representative simulation results are presented for
performance validation. In the considered scenario, UEs are distributed in a
$5$ m $\times$ $5$ m square area. The AP is installed on the wall at
$(1.3,0,2.7)$ with $25$ antenna elements, and the separation between antenna
elements is set as $2$ mm. The RIS is installed on the wall at $(0,1,2.5)$
with $3.33$ mm separations of reflecting elements. The noise power is set to
$-120$ dBm. The transmission carrier frequency is $90$ GHz, $G_{t}=40$ dBi,
$G_{r}=50$ dBi and transmit power is $27$ dBm. Unless specified otherwise, the
RIS panel is equipped with $11\times 101=1111$ reflecting elements with
Fraunhofer distance of $17.58$ m, the number of RIS training phase shift
vectors is $S=150$, the number of pilot data is $\tau=20$. All results are
obtained by averaging over $200$ random generations of UE locations. For
performance comparison, we utilize the conventional far-field model with the
MUSIC method as the performance benchmark, which is labeled as “FF”. Note that
the UEs are supposed to be on the XOY plane, i.e., $z_{u}=0$, as in this case,
the far-field method can only obtain the 2D localization. The proposed NF-JCEL
algorithm is labeled as “Prop”.

Figure 6: RMSE of AoA $\omega$ versus number of RIS elements. Figure 7: RMSE
of AoA $\varphi$ versus number of RIS elements. Figure 8: RMSE of location
versus number of RIS elements. Figure 9: RMSE of channel gain versus number of
RIS elements.

###  IV-A Impact of the Number of RIS Elements

Fig. 6-Fig.9 show the RMSE performances as functions of the numbers of
reflecting elements, where the Fraunhofer distance corresponding to the number
of reflecting elements is from $0.4823$ m to $8.1$ m. The number of UE is $2$,
one is a close UE and the other is a far UE, resulting in three RMSE curves:
the respective RMSEs of these two UEs and the sum of two RMSEs, which are
labeled as “RMSE of near UE”, “RMSE of far UE” and “Sum-RMSE”, respectively.
Unless specified otherwise, the same configuration and labels are adopted in
the following subsections.

Fig. 6 shows the RMSEs of the estimated AoA $\omega_{u}=\sin\phi_{u}$. It can
be seen from Fig. 6 that the sum-RMSE of the proposed algorithm decreases with
increasing the number of reflecting elements, while the sum-RMSE of the far-
field algorithm increases with the number of reflecting elements. This means
that the near-field effect increases with increasing the number of reflecting
elements, i.e. the panel size of RIS. To be more specific, the channel model
approximation error of the far-field algorithm increases with the number of
reflecting elements. Fig. 7 shows the RMSEs of the estimated AoA
$\varphi_{u}=\sin\theta_{u}\cos\phi_{u}$. From Fig. 7, similar trends in the
sum RMSE can be observed to those from Fig. 6.

Fig. 8 shows the localization performance in terms of RMSE. It can be seen
from Fig. 8 that the proposed algorithm outperforms the corresponding far-
field cases, and the performance gap increases with the number of reflecting
elements. This implies that the proposed algorithm can fully leverage the
information provided in the near-field scenario. The RMSE of the near UE is
less than that of the far UE for both far-field and proposed near-field
algorithms, as the localization RMSE is also heavily affected by distance-
dependent fading. Fig. 9 shows the RMSE performance of the channel fading
coefficients. Similar trends to Fig. 8 can be seen from Fig. 9. In Fig. 9, it
is observed that the RMSE of the proposed algorithm is dominated by the far
UEs, as the channel fading is more severe for the longer transmission
distance.

###  IV-B Impact of Number of UEs

Figure 10: RMSE of AoA $\omega$ versus number of UEs. Figure 11: RMSE of AoA
$\varphi$ versus number of UEs.

Fig. 10 shows the RMSE of $\omega_{u}=\sin\phi_{u}$ versus the number of UEs,
where the number of UEs gradually increases by adding UEs moving away from the
RIS panel. The sum of RMSEs for all UEs, the worst RMSE among all UEs, and the
average RMSE per UE are shown, which are labeled as “RMSE Sum”, “RMSE Max” and
“RMSE AvrgUE”, respectively. It can be seen that the proposed algorithm always
performs better than the far-field case. The sum RMSEs of the proposed
algorithm increase with the number of UEs, as more UEs need to be located. In
addition, the sum-RMSE of AoA estimation is dominated by the worst RMSE in the
proposed algorithm, and the performance gap between the UE’s average RMSE and
the sum-RMSE increases with the number of far-away UEs. Meanwhile, the“RMSE
Max” and “RMSE AvrgUE” of the far-field algorithm slightly decreases with
increasing number of UEs. The reason is that more far-away UEs are added in
the simulation, as the increased number of distant UEs do not result in more
severe near-field effects.

Fig. 11 shows the RMSE of $\varphi_{u}=\sin\theta_{u}\cos\phi_{u}$. Similar
trends of the far-field algorithm can be seen from Fig. 11 to those from Fig.
10. Also, the sum RMSEs of the proposed algorithm increase with the number of
UEs, and the main reason is that UEs’ separations in the angular domain
decreases with the number of UEs.

Figure 12: RMSE of location versus number of UEs.

Fig. 12 shows the localization performance. The localization RMSEs obtained by
the proposed algorithm increase with the number of UEs, but the proposed
algorithm shows better performance than the far-field algorithm. Meanwhile,
the“RMSE Max” and “RMSE AvrgUE” of the far-field algorithm decreases slightly.
This is consistent with Fig. 10 and Fig. 11. Furthermore, the performance gap
between the far-field algorithm and the proposed algorithm decreases. This is
because the added UE is gradually far away from RIS panels so that the near
field effect gradually becomes insignificant.

###  IV-C Impact of Transmit Power

Figure 13: RMSE of location Figure 14: RMSE of channel fading cofficients

Fig.13and Fig. 14 show the impact of the transmit power on the RMSE
performance. It can be seen from Fig. 13 that the proposed algorithm
outperforms the far-field algorithm for all cases. The RMSEs of the proposed
algorithm decrease with the transmit power, while the RMSE of the far-field
algorithm keeps flat. The sum RMSE of the proposed algorithm is dominated by
that of the far UE, and the RMSE of the near UE decreases faster than the far
UE. Fig. 14 shows the RMSE of the channel fading coefficients, and similar
trends to those in Fig. 13 are illustrated.

###  IV-D Impact of Pilot Overhead

Fig. 15 and Fig. 16 show the impact of the number of RIS training phase shift
vectors, i.e. $S$, on the RMSE performance. It can be seen from Fig. 15 that
the proposed algorithm outperforms the far-field algorithm. The RMSE of
locations obtained by the proposed algorithm decreases rapidly when the number
of training phase shift vectors increases to $120$. Then, the decent speeds of
RMSE slow down for the proposed algorithm. Meanwhile, the RMSE of the far-
field algorithm keeps flat with the number of training phase shift vectors.
Fig. 16 shows the RMSE of the channel fading coefficients. This is consistent
with Fig. 15. The reason can be explained as follows. With more training phase
shift vectors, i.e., the increased number of received pilot signals, the
proposed algorithm can fully exploit the near-field information, while the
performance of the far-field algorithm cannot be improved.

Figure 15: RMSE of location versus the number of RIS training phase shift
vectors. Figure 16: RMSE of channel fading coefficients versus the number of
RIS training phase shift vectors.

##  V Conclusion

In this paper, we have investigated the localization and CSI estimation scheme
for the near-field subTHz system with RIS. The following conclusions are
drawn:

  * •

The proposed NF-JCEL algorithm shows attractive performance in terms of
localization and CSI estimation RMSE, especially for the near-field
transmission introduced by the XL-RIS panel. Meanwhile, the conventional far-
field algorithm shows severe performance degradation in the near-field
channel.

  * •

The complexity in the near-field CSI estimation is highly dependent on the
array steering vector formulation, where the UE distance and the AoAs are
different by reflecting elements and jointly coupled together. However, this
also leads to higher resolution accuracy. Thus, the near-field channel model
can provide 3D localization with a single RIS panel, which is not applicable
for the far-field model with one RIS panel in the same case.

  * •

The localization accuracy is also highly dependent on the angle separations
between different UEs, which can be improved with a large RIS panel with more
elements. However, without considering the near-field effect, the more
reflecting elements cannot bring more benefits but severe performance
degradation. Thus, it is a necessity to consider the spherical wavefront
feature for the high precision localization system.

## References

  * [1] X. You, C.-X. Wang, J. Huang, X. Gao, Z. Zhang, M. Wang, Y. Huang, C. Zhang, Y. Jiang, J. Wang _et al._ , “Towards 6G wireless communication networks: Vision, enabling technologies, and new paradigm shifts,” _Science China Information Sciences_ , vol. 64, no. 1, pp. 1–74, 2021.
  * [2] O. Kanhere and T. S. Rappaport, “Position location for futuristic cellular communications: 5G and beyond,” _IEEE communications magazine_ , vol. 59, no. 1, pp. 70–75, 2021.
  * [3] Y. Xing and T. S. Rappaport, “Millimeter wave and terahertz urban microcell propagation measurements and models,” _IEEE Communications Letters_ , vol. 25, no. 12, pp. 3755–3759, 2021.
  * [4] B. Ning, Z. Chen, W. Chen, Y. Du, and J. Fang, “Terahertz multi-user massive MIMO with intelligent reflecting surface: Beam training and hybrid beamforming,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 2, pp. 1376–1393, 2021.
  * [5] Y. Pan, K. Wang, C. Pan, H. Zhu, and J. Wang, “Sum-rate maximization for intelligent reflecting surface assisted terahertz communications,” _IEEE Transactions on Vehicular Technology_ , 2022.
  * [6] C. Pan, H. Ren, K. Wang, M. Elkashlan, A. Nallanathan, J. Wang, and L. Hanzo, “Intelligent reflecting surface aided MIMO broadcasting for simultaneous wireless information and power transfer,” _IEEE Journal on Selected Areas in Communications_ , vol. 38, no. 8, pp. 1719–1734, 2020.
  * [7] Y. Pan, K. Wang, C. Pan, H. Zhu, and J. Wang, “Self-sustainable reconfigurable intelligent surface aided simultaneous terahertz information and power transfer (stipt),” _IEEE Transactions on Wireless Communications_ , 2022\.
  * [8] C. Pan, H. Ren, K. Wang, W. Xu, M. Elkashlan, A. Nallanathan, and L. Hanzo, “Multicell MIMO communications relying on intelligent reflecting surfaces,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 8, pp. 5218–5233, 2020.
  * [9] T. Wu, C. Pan, Y. Pan, S. Hong, H. Ren, M. Elkashlan, F. Shu, and J. Wang, “Two-step mmwave positioning scheme with ris-part i: Angle estimation and analysis,” _arXiv preprint arXiv:2208.07602_ , 2022.
  * [10] C. Pan, G. Zhou, K. Zhi, S. Hong, T. Wu, Y. Pan, H. Ren, M. D. Renzo, A. Lee Swindlehurst, R. Zhang, and A. Y. Zhang, “An overview of signal processing techniques for RIS/IRS-aided wireless systems,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 16, no. 5, pp. 883–917, 2022.
  * [11] R. Wang, Z. Xing, and E. Liu, “Joint location and communication study for intelligent reflecting surface aided wireless communication system.” [Online]. Available: http://arxiv.org/pdf/2103.01063v1
  * [12] Y. T. Chan and K. C. Ho, “A simple and efficient estimator for hyperbolic location,” _IEEE Transactions on Signal Processing_ , vol. 42, no. 8, pp. 1905–1915, 1994.
  * [13] H. Wymeersch, J. He, B. Denis, A. Clemente, and M. Juntti, “Radio localization and mapping with reconfigurable intelligent surfaces.” [Online]. Available: https://arxiv.org/pdf/1912.09401.pdf
  * [14] H. Zhang, Di Boya, K. Bian, Z. Han, H. V. Poor, and L. Song, “Towards ubiquitous sensing and localization with reconfigurable intelligent surfaces.” [Online]. Available: http://arxiv.org/pdf/2201.10101v1
  * [15] C. Pan, G. Zhou, K. Zhi, S. Hong, T. Wu, Y. Pan, H. Ren, M. Di Renzo, A. L. Swindlehurst, R. Zhang, and A. Y. Zhang, “An overview of signal processing techniques for RIS/IRS-aided wireless systems.” [Online]. Available: https://arxiv.org/pdf/2112.05989
  * [16] G. Zhou, C. Pan, H. Ren, P. Popovski, and A. L. Swindlehurst, “Channel estimation for RIS-aided multiuser millimeter-wave systems.” [Online]. Available: http://arxiv.org/pdf/2106.14792v2
  * [17] W. Wang and W. Zhang, “Joint beam training and positioning for intelligent reflecting surfaces assisted millimeter wave communications,” _IEEE Transactions on Wireless Communications_ , p. 1, 2021.
  * [18] H. Zhang, H. Zhang, B. Di, K. Bian, Z. Han, and L. Song, “Metalocalization: Reconfigurable intelligent surface aided multi-user wireless indoor localization,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 12, pp. 7743–7757, 2021.
  * [19] J. He, H. Wymeersch, L. Kong, O. Silvén, and M. Juntti, “Large intelligent surface for positioning in millimeter wave MIMO systems,” in _2020 IEEE 91st Vehicular Technology Conference (VTC2020-Spring)_ , 2020, pp. 1–5.
  * [20] K. T. Selvan and R. Janaswamy, “Fraunhofer and fresnel distances: Unified derivation for aperture antennas.” _IEEE Antennas and Propagation Magazine_ , vol. 59, no. 4, pp. 12–15, 2017.
  * [21] J. C. Chen, R. E. Hudson, and K. Yao, “Maximum-likelihood source localization and unknown sensor location estimation for wideband signals in the near-field,” _IEEE Transactions on Signal Processing_ , vol. 50, no. 8, pp. 1843–1854, 2002.
  * [22] B. Friedlander, “Localization of signals in the near-field of an antenna array,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 15, pp. 3885–3893, 2019.
  * [23] W. Tang, M. Z. Chen, X. Chen, J. Y. Dai, Y. Han, M. Di Renzo, Y. Zeng, S. Jin, Q. Cheng, and T. J. Cui, “Wireless communications with reconfigurable intelligent surface: Path loss modeling and experimental measurement,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 1, pp. 421–439, 2021.
  * [24] E. Björnson, Ö. T. Demir, and L. Sanguinetti, “A primer on near-field beamforming for arrays and reconfigurable intelligent surfaces.” [Online]. Available: https://arxiv.org/pdf/2110.06661
  * [25] H. Lu and Y. Zeng, “Communicating with extremely large-scale array/surface: Unified modelling and performance analysis,” _IEEE Transactions on Wireless Communications_ , pp. 1–1, 2021.
  * [26] H. Chen, A. Elzanaty, R. Ghazalian, M. F. Keskin, R. Jäntti, and H. Wymeersch, “Channel model mismatch analysis for XL-MIMO systems from a localization perspective,” _arXiv preprint arXiv:2205.15417_ , 2022.
  * [27] X. Wei, L. Dai, Y. Zhao, G. Yu, and X. Duan, “Codebook design and beam training for extremely large-scale RIS: Far-field or near-field?” [Online]. Available: http://arxiv.org/pdf/2109.10143v1
  * [28] J. C. B. Garcia, A. Sibille, and M. Kamoun, “Reconfigurable intelligent surfaces: Bridging the gap between scattering and reflection,” _IEEE Journal on Selected Areas in Communications_ , vol. 38, no. 11, pp. 2538–2547, 2020.
  * [29] Z. Abu-Shaban, K. Keykhosravi, M. F. Keskin, G. C. Alexandropoulos, G. Seco-Granados, and H. Wymeersch, “Near-field localization with a reconfigurable intelligent surface acting as lens,” in _ICC 2021 - IEEE International Conference on Communications_. IEEE, 2021, pp. 1–6.
  * [30] A. Elzanaty, A. Guerra, F. Guidi, and M. Alouini, “Reconfigurable intelligent surfaces for localization: Position and orientation error bounds,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 5386–5402, 2021.
  * [31] M. Rahal, B. Denis, K. Keykhosravi, B. Uguen, and H. Wymeersch, “RIS-enabled localization continuity under near-field conditions.” [Online]. Available: http://arxiv.org/pdf/2109.11965v1
  * [32] X. Zhang and H. Zhang, “Hybrid reconfigurable intelligent surfaces-assisted near-field localization,” _IEEE Communications Letters_ , pp. 1–1, 2022\.
  * [33] Z. Chen, P. Chen, Z. Guo, and X. Wang, “A RIS-based passive DOA estimation method for integrated sensing and communication system,” _arXiv preprint arXiv:2204.11626_ , 2022.
  * [34] O. Rinchi, A. Elzanaty, and M.-S. Alouini, “Compressive near-field localization for multipath RIS-aided environments,” _IEEE Communications Letters_ , 2022.
  * [35] Y. Pan, C. Pan, S. Jin, and J. Wang, “Localization in the near field of a RIS-assisted mmwave/subTHz system,” _To apprear in IEEE Globecom_ , 2022\.
  * [36] W. Hao, G. Sun, M. Zeng, Z. Zhu, Z. Chu, O. A. Dobre, and P. Xiao, “Robust design for intelligent reflecting surface assisted MIMO-OFDMA terahertz communications.” [Online]. Available: https://arxiv.org/pdf/2009.05893
  * [37] C. Chaccour, M. N. Soorki, W. Saad, M. Bennis, and P. Popovski, “Can terahertz provide high-rate reliable low latency communications for wireless VR?” [Online]. Available: https://arxiv.org/pdf/2005.00536
  * [38] M. Di Renzo, A. Zappone, M. Debbah, M.-S. Alouini, C. Yuen, J. de Rosny, and S. Tretyakov, “Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead,” _IEEE Journal on Selected Areas in Communications_ , vol. 38, no. 11, pp. 2450–2525, 2020.
  * [39] Z. Zheng, M. Fu, W.-Q. Wang, S. Zhang, and Y. Liao, “Localization of mixed near-field and far-field sources using symmetric double-nested arrays,” _IEEE Transactions on Antennas and Propagation_ , vol. 67, no. 11, pp. 7059–7070, 2019.
  * [40] J.-W. Tao, L. Liu, and Z.-Y. Lin, “Joint DOA, range, and polarization estimation in the fresnel region,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 47, no. 4, pp. 2657–2672, 2011.
  * [41] W. Zuo, J. Xin, N. Zheng, and A. Sano, “Subspace-based localization of far-field and near-field signals without eigendecomposition,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 17, pp. 4461–4476, 2018.
  * [42] F. Bohagen, P. Orten, and G. E. Oien, “Design of optimal high-rank line-of-sight MIMO channels,” _IEEE Transactions on Wireless Communications_ , vol. 6, no. 4, pp. 1420–1425, 2007.
  * [43] F. Bohagen, P. Orten, and G. Oien, “On spherical vs. plane wave modeling of line-of-sight MIMO channels,” _IEEE Transactions on Communications_ , vol. 57, no. 3, pp. 841–849, 2009.
  * [44] F. Bøhagen, P. Orten, and G. Øien, “Optimal design of uniform rectangular antenna arrays for strong line-of-sight MIMO channels,” _EURASIP Journal on Wireless Communications and Networking_ , vol. 2007, pp. 1–10, 2007.
  * [45] N. K. Kundu and M. R. McKay, “Channel estimation for reconfigurable intelligent surface aided MISO communications: From LMMSE to deep learning solutions,” _IEEE Open J. Commun. Soc._ , vol. 2, pp. 471–487, Mar. 2021.
  * [46] J. Liang and D. Liu, “Passive localization of mixed near-field and far-field sources using two-stage MUSIC algorithm,” _IEEE Transactions on Signal Processing_ , vol. 58, no. 1, pp. 108–120, 2010.
  * [47] J. Xin and A. Sano, “Computationally efficient subspace-based method for direction-of-arrival estimation without eigendecomposition,” _IEEE Transactions on Signal Processing_ , vol. 52, no. 4, pp. 876–893, 2004.
  * [48] S. Marcos, A. Marsal, and M. Benidir, “The propagator method for source bearing estimation,” _Signal Processing_ , vol. 42, no. 2, pp. 121–138, 1995\. [Online]. Available: https://www.sciencedirect.com/science/article/pii/016516849400122G
  * [49] J. Yadav and M. Sharma, “A review of k-mean algorithm,” _Int. J. Eng. Trends Technol_ , vol. 4, no. 7, pp. 2972–2976, 2013.