# Toward Wireless Localization Using Multiple Reconfigurable Intelligent
Surfaces
††thanks: This work is supported by the National Natural Science Foundation of
China under Grant No.12141107, and the Interdisciplinary Research Program of
HUST (2023JCYJ012). (Corresponding author: Tiebin Mi)

Fuhai Wang, Tiebin Mi, Chun Wang, Rujing Xiong, Zhengyu Wang, and Robert
Caiming Qiu Fuhai Wang is with the School of Electronic Information and
Communications, Huazhong University of Science and Technology, Wuhan 430074,
China, and also with the Institute of Artificial Intelligence, Huazhong
University of Science and Technology, Wuhan 430074, China (e-mail:
wangfuhai@hust.edu.cn).Tiebin Mi, Chun Wang, Rujing Xiong, Zhengyu Wang and
Robert Caiming Qiu are with the School of Electronic Information and
Communications, Huazhong University of Science and Technology, Wuhan 430074,
China (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>caiming@hust.edu.cn).

###### Abstract

This paper investigates the capabilities and effectiveness of backward sensing
centered on reconfigurable intelligent surfaces (RISs). We demonstrate that
the direction of arrival (DoA) estimation of incident waves in the far-field
regime can be accomplished using a single RIS by leveraging configurational
diversity. Furthermore, we identify that the spatial diversity achieved
through deploying multiple RISs enables accurate localization of multiple
power sources. Physically accurate and mathematically concise models are
introduced to characterize forward signal aggregations via RISs. By employing
linearized approximations inherent in the far-field region, the measurement
process for various configurations can be expressed as a system of linear
equations. The mathematical essence of backward sensing lies in solving this
system. A theoretical framework for determining key performance indicators is
established through condition number analysis of the sensing operators. In the
context of localization using multiple RISs, we examine relationships among
the rank of sensing operators, the size of the region of interest (RoI), and
the number of elements and measurements. For DoA estimations, we provide an
upper bound for the relative error of the least squares reconstruction
algorithm. These quantitative analyses offer essential insights for system
design and optimization. Numerical experiments validate our findings. To
demonstrate the practicality of our proposed RIS-centric sensing approach, we
develop a proof-of-concept prototype using universal software radio
peripherals (USRP) and employ a magnitude-only reconstruction algorithm
tailored for this system. To our knowledge, this represents the first trial of
its kind.

###### Index Terms:

Reconfigurable Intelligent Surface (RIS), forward signal aggregations,
backward sensing, localization, direction of arrival (DoA).

##  I Introduction

Reconfigurable intelligent surfaces (RISs) have received considerable
attention due to their potential to enhance both the communication and sensing
capabilities of wireless networks [1, 2, 3, 4, 5, 6]. RISs comprise cost-
effective, well-designed electromagnetic (EM) units, each capable of
independently modifying the characteristics of incident EM waves. These units
enable the control of the scattering and reflection properties of EM waves
through software-based manipulation. Although current research extensively
explores their applications in forward coverage enhancement under line-of-
sight/non line-of-sight (LoS/NLoS) conditions [7, 8, 9, 10, 11, 12, 3, 13],
backward sensing through RISs is also attracting increasing attention [14, 15,
16, 17]. A key advantage of RISs in sensing applications is their
reconfigurability of the electromagnetic environment, which distinguishes them
from other competitors [18, 19, 20] in the design of sensing signals with
temporal diversity or frequency diversity.

The diversity schemes introduced by RISs are twofold: configurational
diversity [21, 22, 23], enabled by their reconfigurability, and spatial
diversity, achieved through deploying multiple RISs. Current research
primarily focuses on exploring the configurational diversity of an individual
RIS. In these studies, the RIS is often regarded as a promising auxiliary
technology for enhancing wireless systems by producing a sharp point spread
function [24, 25, 26, 27, 28], which facilitates the generation of an
unblurred image or clear map (e.g., [29, 15]). Notably, RISs themselves are
well-suited for sensing tasks due to their nature as large and reconfigurable
antenna arrays. In general, RISs can generate a series of high-dimensional
observations for regions of interest (RoIs) [30], which greatly benefit
sensing applications [31]. Nevertheless, the complete functionality of RISs in
this context remains unclear. It is crucial to assess the feasibility and
capabilities of the RIS-centric backward sensing approach.

A fundamental function of sensing is localization, which involves extracting
location-related information such as the direction of arrival (DoA), also
known as the angle of arrival (AoA) [32, 33]. However, there are only a
limited number of articles addressing DoA estimation through RISs [34, 15,
35]. To the best of our knowledge, most research is based on numerical
simulations, and there are no prototype demonstrations of this capability. One
reason is that accurate physics-based models of RISs in backward sensing are
not yet well-developed. In this paper, we first demonstrate that DoA
estimation of incident waves in the far-field can be performed using a single
RIS, leveraging configurational diversity.

Even with successful DoA estimation using RISs, accurate localization remains
challenging [36]. One approach involves deploying multiple RISs [37, 29] to
enhance the spatial diversity. A DoA parameter defines a straight line passing
through the power source and the RIS, meaning that two DoA parameters are
sufficient to locate a single source via triangulation, as shown in Fig. 1.
However, when multiple sources radiate signals simultaneously, ambiguity
arises, as illustrated in Fig. 1. To resolve this uncertainty, additional RISs
positioned at different angles can be used. This approach is analogous to
computed tomography, which reconstructs a field based on dense angular
samplings [38, 39]. However, the cost of such samplings via RISs is
impractical. Fully utilizing the spatial diversity provided by deploying
multiple RISs remains an open question. Addressing these challenges in RIS-
centric backward sensing requires a paradigm shift.

Figure 1: Illustration of the RIS-centric sensing system via triangulation:
(a) A single source can be localized using two DoA parameters. (b) Ambiguity
arises when multiple sources radiate signals simultaneously.

###  I-A Contributions

The main contributions are summarized as follows:

  * •

Physics-based forward signal aggregations and mathematical formulations for
backward sensing: Physically accurate and mathematically concise models are
presented to characterize the reflection properties of RISs. By employing
linearized approximations inherent in the far-field region, we provide
practical models to describe the forward signal aggregation process via RISs.
Additionally, we rewrite the measurements of the RoI for various
configurations as a system of linear equations. This formulation offers
valuable insights into the capabilities of the proposed RIS-centric sensing
approach. Mathematically, the problem of backward sensing can be formulated as
finding a reasonable solution to this system.

  * •

Key indicators dominating the performance of backward sensing through RISs: A
theoretical framework to identify key performance indicators that dominate the
performance of backward sensing is established through condition number
analysis of the sensing operators. In the context of localization using
multiple RISs, we investigate the relationships among the rank of sensing
operators, the size of the RoI, and the number of elements and measurements. A
necessary condition for high-fidelity recovery is that the dimension of the
RoI is less than the number of elements and measurements. For DoA estimations,
we provide an upper bound for the relative error of the least squares
reconstruction algorithm. These quantitative analyses offer crucial insights
for system design and optimization.

  * •

Both numerical simulations and proof-of-concept prototypes demonstrate the
feasibility of the proposed approach: We validate the capabilities and
effectiveness of the proposed RIS-centric backward sensing approach through
extensive experiments. Numerical trials illustrate the impact of several
factors, including the number of elements, measurements, and RISs, as well as
their positions. To demonstrate the practicality of our approach, we develop a
proof-of-concept prototype using universal software radio peripherals (USRP)
and employ a magnitude-only reconstruction algorithm tailored to this system.
To the best of our knowledge, no previous prototypes have been proposed with
such configurations.

###  I-B Outline

The remainder of the paper is organized as follows. In Section II, we present
forward signal aggregation models using the physical optics method. Section
III addresses the RIS-centric backward sensing problem. Section IV discusses
key indicators dominating the performance of backward sensing. We validate the
findings through extensive numerical experiments in Section V. To demonstrate
the practicality, we develop a proof-of-concept prototype using universal
software radio peripherals (USRP) in Section VI. Finally, the paper is
concluded in Section VII.

###  I-C Reproducible Research

Once the paper is accepted, we will release the code.

###  I-D Notations

Unless explicitly specified, bold capital letters and bold small letters
denote matrices and vectors, respectively. The conjugate transpose, transpose,
and pseudo-inverse of $\mathbf{A}$ are denoted by $\mathbf{A}^{*}$,
$\mathbf{A}^{\top}$ and $\mathbf{A}^{{\dagger}}$, respectively. We denote by
$\|\cdot\|$ the Euclidean norm of a vector and the Frobenius norm of a matrix.
The matrix $\text{diag}(\mathbf{a})$ represents a diagonal matrix with
diagonal elements given by $\mathbf{a}$. Additionally, rank{$\mathbf{A}$}
denotes the rank of the matrix $\mathbf{A}$.

##  II Forward Signal Aggregations via A Single RIS

In the realm of RIS-centric inverse or backward sensing, a deep understanding
of the forward measurement process is essential. Mathematically, this process
is represented by the relation $S_{t}=H_{t}(\mathbf{E})$, where the scalar
$S_{t}$ denotes the observation, the vector $\mathbf{E}$ represents the scene
under investigation, and $H_{t}$ is the sensing operator mapping the scene to
the observation. Here, we employ $t$ to index measurements corresponding to
various configurations. For simplicity, we set the configuration to random
sampling mode. Generally, RISs with random phase configurations are referred
to as statistically isotropic RISs. While this configuration may not fully
utilize the available degrees of freedom, it serves as a widely accepted
benchmark for evaluating and comparing the performance of other RIS
configurations [40, 24, 41].

###  II-A Linear Aggregations of Multiple Incident Signals

We now examine the linear aggregation properties of RISs. Our analysis begins
by concentrating on a single RIS situated in the $xoy$-plane, composed of
multiple units located at $\mathbf{p}_{n}=[x_{n},\ y_{n},\ z_{n}]^{\top}$,
$n=1,\ldots,N$. For simplicity, we assume that the scattering pattern of an
isolated unit is denoted by
$\tau(\theta^{\text{s}},\phi^{\text{s}};\theta^{\text{i}},\phi^{\text{i}})$.
This notation signifies the dependence on both the incident angle
$(\theta^{\text{i}},\phi^{\text{i}})$ and the scattered angle
$(\theta^{\text{s}},\phi^{\text{s}})$. We will not delve into the specific
form of $\tau(\cdot)$, as our goal is to illustrate the linearity of RISs.
Using the physical optics method, we can calculate the bistatic scattered
field of a rectangular metallic patch with near-zero thickness, as detailed in
[40].

Let us consider an example illustrated in Fig. 2. The RIS is illuminated by a
single incident EM wave originating from
$(r^{\text{i}},\theta^{\text{i}},\phi^{\text{i}})$. Suppose the incident
electric field arriving at $\mathbf{p}_{n}$ along
$(\theta^{\text{i}}(n),\phi^{\text{i}}(n))$ is denoted by $E^{\text{i}}(n)$.
Consequently, the scattered filed at
$(\theta^{\text{s}}(n),\phi^{\text{s}}(n))$ is represented as
$E^{\text{s}}(n)=\tau(\theta^{\text{s}}(n),\phi^{\text{s}}(n);\theta^{\text{i}}(n),\phi^{\text{i}}(n))E^{\text{i}}(n)$.
In the case of a point source, radiation propagates radially. As the EM wave
propagates toward the unit at $\mathbf{p}_{n}$, the attenuation behavior is
characterized by the factor $r^{\text{i}}(n)e^{j2\pi
r^{\text{i}}(n)/\lambda}$. Therefore,
$E^{\text{i}}(n)=E(r^{\text{i}},\theta^{\text{i}},\phi^{\text{i}})e^{-j2\pi
r^{\text{i}}(n)/\lambda}/r^{\text{i}}(n)$. Similarly, if $E^{\text{s}}(n)$
denotes the scattered electric field of the unit at $\mathbf{p}_{n}$, then the
electric field at the observation point is $S_{n}=E^{\text{s}}(n)e^{-j2\pi
r^{\text{s}}(n)/\lambda}/r^{\text{s}}(n)$.

Figure 2: An RIS illuminated by a single incident EM wave from
$(r^{\text{i}},\theta^{\text{i}},\phi^{\text{i}})$. The scattered waves from
various elements are aggregated at
$(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}})$.

Combining the three components along the propagation trace through the $n$-th
unit at $\mathbf{p}_{n}$, the electric field at the observation point is
expressed as

$S_{n}(r^{\text{s}}_{n},\theta^{\text{s}},\phi^{\text{s}})=\tau_{n}(\theta^{\text{s}}(n),\phi^{\text{s}}(n);\theta^{\text{i}}(n),\phi^{\text{i}}(n))e^{j\Omega_{n}}\\\
\frac{e^{-j2\pi r^{\text{i}}(n)/\lambda}e^{-j2\pi
r^{\text{s}}(n)/\lambda}}{r^{\text{i}}(n)r^{\text{s}}(n)}E(r^{\text{i}},\theta^{\text{i}},\phi^{\text{i}}).$
(1)

Here, $\Omega_{n}$ represents the phase configuration of the $n$-th unit
within the RIS. The aggregated electric field at the observation point is
determined by the superposition of individual fields scattered by $N$ units,
given by

$S(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}})=E(r^{\text{i}},\theta^{\text{i}},\phi^{\text{i}})\sum_{n=1}^{N}\tau_{n}(\theta^{\text{s}}(n),\phi^{\text{s}}(n);\theta^{\text{i}}(n),\phi^{\text{i}}(n))\\\
e^{j\Omega_{n}}\frac{e^{-j2\pi r^{\text{i}}(n)/\lambda}e^{-j2\pi
r^{\text{s}}(n)/\lambda}}{r^{\text{i}}(n)r^{\text{s}}(n)}.$ (2)

When the RIS is exposed to multiple incident waves, the aggregated electric
field is the superposition of the individual fields

$\displaystyle S(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}})=$
$\displaystyle\sum_{m=1}^{M}E(r^{\text{i}}_{m},\theta^{\text{i}}_{m},\phi^{\text{i}}_{m})$
(3)
$\displaystyle\sum_{n=1}^{N}\tau(\theta^{\text{s}}(n),\phi^{\text{s}}(n);\theta^{\text{i}}_{m}(n),\phi^{\text{i}}_{m}(n))$
$\displaystyle e^{j\Omega_{n}}\frac{e^{-j2\pi
r^{\text{i}}_{m}(n)/\lambda}e^{-j2\pi
r^{\text{s}}(n)/\lambda}}{r^{\text{i}}_{m}(n)r^{\text{s}}(n)}.$

$\displaystyle S(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}})\approx$
$\displaystyle\tau\frac{e^{-j2\pi
r^{\text{s}}/\lambda}}{r^{\text{s}}}\begin{bmatrix}e^{j2\pi\mathbf{p}_{1}^{\top}\mathbf{u}(\theta^{\text{s}},\phi^{\text{s}})/\lambda}&\cdots&e^{j2\pi\mathbf{p}_{N}^{\top}\mathbf{u}(\theta^{\text{s}},\phi^{\text{s}})/\lambda}\end{bmatrix}\begin{bmatrix}e^{j\Omega_{1}}&&0\\\
&\ddots&\\\ 0&&e^{j\Omega_{N}}\end{bmatrix}$ (4)
$\displaystyle\begin{bmatrix}e^{j2\pi\mathbf{p}_{1}^{\top}\mathbf{u}(\theta_{1}^{\text{i}},\phi_{1}^{\text{i}})/\lambda}&\cdots&e^{j2\pi\mathbf{p}_{1}^{\top}\mathbf{u}(\theta_{M}^{\text{i}},\phi_{M}^{\text{i}})/\lambda}\\\
\vdots&\ddots&\vdots\\\
e^{j2\pi\mathbf{p}_{N}^{\top}\mathbf{u}(\theta_{1}^{\text{i}},\phi_{1}^{\text{i}})/\lambda}&\cdots&e^{j2\pi\mathbf{p}_{N}^{\top}\mathbf{u}(\theta_{M}^{\text{i}},\phi_{M}^{\text{i}})/\lambda}\\\
\end{bmatrix}\begin{bmatrix}\frac{e^{-j2\pi
r^{\text{i}}_{1}/\lambda}}{r^{\text{i}}_{1}}&&0\\\ &\ddots&\\\
0&&\frac{e^{-j2\pi
r^{\text{i}}_{M}/\lambda}}{r^{\text{i}}_{M}}\end{bmatrix}\begin{bmatrix}E(r_{1}^{\text{i}},\theta_{1}^{\text{i}},\phi_{1}^{\text{i}})\\\
\vdots\\\
E(r_{M}^{\text{i}},\theta_{M}^{\text{i}},\phi_{M}^{\text{i}})\end{bmatrix}.$

The input-output behavior expression can be simplified under specific
scenarios, particularly when the RIS comprises isotropic units. In isotropic
scattering, incident EM waves uniformly scatter in all directions over the
hemisphere of reflection, as depicted in Fig. 3, irrespective of the angle of
incidence or observation. The hypothetical isotropic element is sometimes
referred to as Lambertian scatterers [42]. Consequently, the unit’s scattering
pattern simplifies to
$\tau_{n}(\theta^{\text{s}},\phi^{\text{s}};\theta^{\text{i}},\phi^{\text{i}})=\tau$.
The aggregated electric field is given as

$\displaystyle S(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}})=$
$\displaystyle\sum_{m=1}^{M}E(r^{\text{i}}_{m},\theta^{\text{i}}_{m},\phi^{\text{i}}_{m})$
(5) $\displaystyle\sum_{n=1}^{N}\tau e^{j\Omega_{n}}\frac{e^{-j2\pi
r^{\text{i}}_{m}(n)/\lambda}e^{-j2\pi
r^{\text{s}}(n)/\lambda}}{r^{\text{i}}_{m}(n)r^{\text{s}}(n)}.$

Figure 3: Lambertian scattering [42]. Incident EM waves uniformly scatter in
all directions over the hemisphere of reflection.

######  Remark 1.

In (3) and (5), we illustrate the linearization approximations inherent in the
functionality of the RIS. Notably, we do not impose any specific topology of
RISs or far-field regime restrictions during the derivation, which allows for
a more general description of their behavior. Understanding the linear
aggregation of forward incident signals offers valuable insights into the
capabilities of RISs, laying a foundation for deploying multiple RISs in
sensing applications.

In the far-field regime, the distances $r^{\text{s}}(n)$ and $r^{\text{i}}(n)$
exhibit linear relationships with $r^{\text{s}}$ and $r^{\text{i}}$,
respectively. This enables concise representations of the linear aggregations
in (3) and (5). Following the approach outlined in [40], we can easily compute
the path length differences $r^{\text{s}}(n)-r^{\text{s}}$ and
$r^{\text{i}}(n)-r^{\text{i}}$. When analyzing uniform planar arrays, which
represent the fundamental topology of RISs, we can employ the canonical linear
representation (4) to describe the input/output behaviors. In this context, we
define $\mathbf{u}(\theta,\phi)=[\sin\theta\cos\phi,\ \sin\theta\sin\phi,\
\cos\theta]^{\top}$.

To shorten notations, we let

$\displaystyle l(r^{\text{s}})$ $\displaystyle=\frac{e^{-j2\pi
r^{\text{s}}/\lambda}}{r^{\text{s}}},$
$\displaystyle\mathbf{E}(\mathbf{r}^{\text{i}},\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})$
$\displaystyle=\left[E(r_{1}^{\text{i}},\theta_{1}^{\text{i}},\phi_{1}^{\text{i}}),\cdots,E(r_{M}^{\text{i}},\theta_{M}^{\text{i}},\phi_{M}^{\text{i}})\right]^{\top},$
$\displaystyle l(\mathbf{r}^{\text{i}})$ $\displaystyle=\left[\frac{e^{-j2\pi
r_{1}^{\text{i}}/\lambda}}{r_{1}^{\text{i}}},\cdots,\frac{e^{-j2\pi
r_{M}^{\text{i}}/\lambda}}{r_{M}^{\text{i}}}\right],$ $\displaystyle
e^{j\mathbf{\Omega}}$
$\displaystyle=\left[e^{j\Omega_{1}},\cdots,e^{j\Omega_{N}}\right],$
$\displaystyle\mathbf{v}(\theta^{\text{s}},\phi^{\text{s}})$
$\displaystyle=\left[e^{j2\pi\mathbf{p}_{1}^{\top}\mathbf{u}(\theta^{\text{s}},\phi^{\text{s}})/\lambda},\cdots,e^{j2\pi\mathbf{p}_{N}^{\top}\mathbf{u}(\theta^{\text{s}},\phi^{\text{s}})/\lambda}\right],$
$\displaystyle\mathbf{V}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})$
$\displaystyle=\begin{bmatrix}e^{j2\pi\mathbf{p}_{1}^{\top}\mathbf{u}(\theta_{1}^{\text{i}},\phi_{1}^{\text{i}})/\lambda}&\cdots&e^{j2\pi\mathbf{p}_{1}^{\top}\mathbf{u}(\theta_{M}^{\text{i}},\phi_{M}^{\text{i}})/\lambda}\\\
\vdots&\ddots&\vdots\\\
e^{j2\pi\mathbf{p}_{N}^{\top}\mathbf{u}(\theta_{1}^{\text{i}},\phi_{1}^{\text{i}})/\lambda}&\cdots&e^{j2\pi\mathbf{p}_{N}^{\top}\mathbf{u}(\theta_{M}^{\text{i}},\phi_{M}^{\text{i}})/\lambda}\\\
\end{bmatrix}.$

Note that $\mathbf{v}(\theta^{\text{s}},\phi^{\text{s}})\in\mathbb{C}^{1\times
N}$ and
$\mathbf{V}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})\in\mathbb{C}^{N\times
M}$. Then, we represent the aggregated field at the observation point as

$S(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}})\approx\overbrace{\tau
l(r^{\text{s}})\mathbf{v}(\theta^{\text{s}},\phi^{\text{s}})\text{diag}(e^{j\mathbf{\Omega}})\mathbf{V}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})}^{\mathbf{h}^{\top}_{\mathbf{\Omega}}}\\\
\underbrace{\text{diag}(l(\mathbf{r}^{\text{i}}))\mathbf{E}(\mathbf{r}^{\text{i}},\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})}_{\mathbf{E}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})}.$
(6)

To simplify notation, we denote the measured field as

$S(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}})=\mathbf{h}^{\top}_{\mathbf{\Omega}}\mathbf{E}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})+\nu.$
(7)

Here, $\nu$ represents the noise, which includes modeling inaccuracies
resulting from approximations in the formulation as well as other inherent
uncertainties in the measurements. Without loss of generality, we assume the
noise follows a distribution with a zero mean and variance $\sigma^{2}$.

###  II-B Configurational Diversity

One significant advantage of the RIS-centric backward sensing approach is its
configurational diversity, which enables the generation of multiple
measurement vectors (several snapshots). We can rewrite the measurement of the
aggregated field corresponding to various configurations as

$S_{t}(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}})=\tau
l(r^{\text{s}})e^{j\mathbf{\Omega}_{t}}\text{diag}(\mathbf{v}(\theta^{\text{s}},\phi^{\text{s}}))\\\
\mathbf{V}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})\mathbf{E}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})+\nu_{t},\
t=1,\ldots,T.$

The $T$ snapshots due to various configurations can be stacked into a column
vector. Thus we have

$\mathbf{S}(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}})=\overbrace{\tau
l(r^{\text{s}})\mathbf{e}^{j\mathbf{\Omega}}\text{diag}\left(\mathbf{v}(\theta^{\text{s}},\phi^{\text{s}})\right)\mathbf{V}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})}^{\mathbf{H}(\mathbf{\Omega})}\\\
\text{diag}(l(\mathbf{r}^{\text{i}}))\mathbf{E}(\mathbf{r}^{\text{i}},\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})+\mathbf{n}.$
(8)

Here
$\mathbf{S}(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}})=\left[S_{1}(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}}),\cdots,S_{T}(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}})\right]^{\top}$,
$\mathbf{e}^{j\mathbf{\Omega}}=\left[e^{j\mathbf{\Omega}_{1}};\cdots;e^{j\mathbf{\Omega}_{T}}\right]$
and $\mathbf{n}=\left[\nu_{1},\cdots,\nu_{T}\right]^{\top}$.

Denoting the measurement operator configured by $\mathbf{\Omega}$ as

$\mathbf{H}(\mathbf{\Omega})=\tau
l(r^{\text{s}})\mathbf{e}^{j\mathbf{\Omega}}\text{diag}\left(\mathbf{v}(\theta^{\text{s}},\phi^{\text{s}})\right)\mathbf{V}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}}),$
(9)

we obtain the canonical representation for the input/output behaviors
corresponding to a single RIS

$\mathbf{S}(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}})=\mathbf{H}(\mathbf{\Omega})\mathbf{E}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})+\mathbf{n}.$
(10)

######  Remark 2.

One advantage of the canonical representation (10) is its use of a simple
system of linear equations to describe the forward signal aggregation process
via RISs under reasonable assumptions. This system is both physically accurate
and mathematically tractable, serving as a fundamental module in the RIS-
centric backward sensing approach. This canonical system facilitates the
design of backward sensing algorithms and enables comprehensive performance
analyses.

######  Remark 3.

The matrix $\mathbf{V}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})$
represents the multiplexing of waves from the sources to each element of the
RIS. A characteristic of this matrix is that if two sources have coherent
positions, such as being aligned along the same line when observed via a
uniform linear RIS,
$\mathbf{V}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})$ can become
rank-deficient, leading to the sensing matrix $\mathbf{H}(\mathbf{\Omega})$
becoming ill-conditioned. This deficiency results in degraded performance in
the backward sensing approach. When deploying RIS, we should carefully
consider the geometry and topology to minimize the risk of ill-conditioning in
$\mathbf{V}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})$.

##  III Backward Sensing Using Multiple RISs

It is important to note that while incident distances are accounted for in the
canonical representations (6) and (8),
$\mathbf{E}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})=\text{diag}(l(\mathbf{r}^{\text{i}}))\mathbf{E}(\mathbf{r}^{\text{i}},\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})$
actually represents the incident fields arriving at the RIS along
$(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})$. In other words,
$l(\mathbf{r}^{\text{i}})$ serves as a distance normalization term. The role
of $l(\mathbf{r}^{\text{i}})$ becomes essential in localizations involving
multiple RISs, where large-scale distance attenuation and phase delay factors
can be detected by different RISs. Before exploring the capabilities provided
by deploying multiple RISs, we first clarify the DoA estimation function using
a single RIS.

###  III-A DoA Estimation via a Single RIS

For the DoA estimation task using a single RIS, the vector
$\mathbf{E}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})$ in the
canonical representation (10) is unknown and needs to be estimated. A
fundamental assumption of this scheme is that the sensing operator
$\mathbf{H}(\mathbf{\Omega})$ is known a priori, which necessitates knowledge
of
$l(r^{\text{s}})\text{diag}\left(\mathbf{v}(\theta^{\text{s}},\phi^{\text{s}})\right)$,
$\mathbf{e}^{j\mathbf{\Omega}}$, and
$\mathbf{V}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})$. The term
$l(r^{\text{s}})\text{diag}\left(\mathbf{v}(\theta^{\text{s}},\phi^{\text{s}})\right)$
is relevant to the receiver’s location, while $\mathbf{e}^{j\mathbf{\Omega}}$
is determined by the configurations. Both terms are available within the
backward sensing scenario.

The situation becomes more complicated for
$\mathbf{V}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})$. This matrix
is governed by the directions of the incident waves, the central task of DoA
estimation. One potential approach involves exploring densely sampled grids
across the RoI. While this manipulation significantly expands the solution
space, the extensive reconfiguration capabilities ensure the problem remains
well-defined under modest assumptions (provided the sampling points are not
too closely spaced). This approach is illustrated in Fig. 4 with a linear RIS
for better understanding. The authors have previously investigated the forward
signal aggregations using linear RISs in [40, 8].

Figure 4: Uniform discretization across the azimuth angle domain (RoI). A
total of $T$ snapshots are captured corresponding to various configurations.

When the sensing operator $\mathbf{H}(\mathbf{\Omega})$ is given, ensuring
that a set of measurements
$\mathbf{S}(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}})$ corresponds
uniquely to $\mathbf{E}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})$
requires that the rank of $\mathbf{H}(\mathbf{\Omega})$ be at least equal to
the dimension of RoI. In a typical scenario where an RIS is configured to
produce a well-conditioned full column rank matrix
$\mathbf{H}(\mathbf{\Omega})$, according to (10), a reliable estimation for
the incident field can be given by

$\widehat{\mathbf{E}}(\mathbf{\Theta}^{\text{i}},\mathbf{\Phi}^{\text{i}})=\bigl{(}\mathbf{H}(\mathbf{\Omega})^{*}\mathbf{H}(\mathbf{\Omega})\bigr{)}^{-1}\mathbf{H}(\mathbf{\Omega})^{*}\mathbf{S}(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}}).$
(11)

######  Remark 4.

Both (10) and (11) indicate that a single RIS is limited to determining the
directions of incident waves, making tasks such as location or imaging
infeasible. For accurate localization, multiple strategically positioned RISs
are required.

###  III-B Localization Using Multiple RISs

Similar to the scenario of DoA estimation, where incident angles are
discretized, the localization using multiple RISs with spatial diversity also
requires volume discretization. A common method involves direct discretization
in Cartesian coordinates, as illustrated in Fig. 5. This technique typically
employs uniform voxels, allowing for straightforward implementation. The
advantage of this strategy lies in its versatility and generalizability.
However, such discretization may lead to ill-conditioning in the sensing
operator, similar to the situation where two sources align along the same line
during sensing with uniform linear RISs. We will explore this topic further in
Section V.

Unlike scenarios involving a single RIS, the operational mode of receivers is
crucial in sensing with multiple RISs. The most straightforward mode is each
RIS having its dedicated receiver, as depicted in Fig. 5. Alternatively, there
is a mode where the sensing system uses a single receiver. In this setup, the
reflected signals from all RISs are aggregated to the receiver, as illustrated
in Fig. 5. Next, we will examine the linear representation of measurement
operators in localization using multiple RISs. The form of this operator
varies based on the operational mode of the receivers.

Figure 5:  Two modes in the RIS-centric sensing system with direct
discretization in Cartesian coordinates, where each voxel represents a defined
volume: (a) Each RIS is equipped with a dedicated receiver. (b) The system
employs a single receiver for aggregating reflected signals from all RISs.

In the context of backward sensing, the geometric and topological
relationships among the RoI, receivers, and RISs are predetermined. For
convenience, we denote by $\mathbf{E}(\mathbf{x},\mathbf{y},\mathbf{z})$ the
incident field in Cartesian coordinates and use the subscript $k$ to index
RISs. Utilizing the canonical aggregation (10), we represent the input/output
behavior corresponding to each individual RIS as

$\mathbf{S}(r^{\text{s}}_{k},\theta^{\text{s}}_{k},\phi^{\text{s}}_{k})=\mathbf{H}(\mathbf{\Omega}_{k})\text{diag}(l(\mathbf{r}^{\text{i}}_{k}))\mathbf{E}(\mathbf{x},\mathbf{y},\mathbf{z})+\mathbf{n}_{k}.$

For the mode where each RIS has its dedicated receiver, the measurement
vectors corresponding to each RIS can be stacked into a column vector. Thus,
we can represent the measurement process associated to the whole system as

$\begin{bmatrix}\mathbf{S}(r^{\text{s}}_{1},\theta^{\text{s}}_{1},\phi^{\text{s}}_{1})\\\
\vdots\\\
\mathbf{S}(r^{\text{s}}_{K},\theta^{\text{s}}_{K},\phi^{\text{s}}_{K})\end{bmatrix}=\begin{bmatrix}\mathbf{H}(\mathbf{\Omega}_{1})\text{diag}(l(\mathbf{r}^{\text{i}}_{1}))\\\
\vdots\\\
\mathbf{H}(\mathbf{\Omega}_{K})\text{diag}(l(\mathbf{r}^{\text{i}}_{K}))\end{bmatrix}\mathbf{E}(\mathbf{x},\mathbf{y},\mathbf{z})+\mathbf{n}.$
(12)

On the other hand, when the sensing system comprises only one receiver that
serves all RISs, we can represent the measurement process

$\displaystyle\mathbf{S}(r^{\text{s}},\theta^{\text{s}},\phi^{\text{s}})=$
$\displaystyle\mathbf{S}(r^{\text{s}}_{1},\theta^{\text{s}}_{1},\phi^{\text{s}}_{1})+\cdots+\mathbf{S}(r^{\text{s}}_{K},\theta^{\text{s}}_{K},\phi^{\text{s}}_{K})$
(13) $\displaystyle=$
$\displaystyle\bigl{[}\mathbf{H}(\mathbf{\Omega}_{1})\text{diag}(l(\mathbf{r}^{\text{i}}_{1}))+\cdots$
$\displaystyle+\mathbf{H}(\mathbf{\Omega}_{K})\text{diag}(l(\mathbf{r}^{\text{i}}_{K}))\bigr{]}\mathbf{E}(\mathbf{x},\mathbf{y},\mathbf{z})+\mathbf{n}.$

######  Remark 5.

Both (12) and (13) are systems of linear equations describing forward signal
aggregation. The distinction between these systems arises from whether each
RIS has its dedicated receiver. In scenarios with specific geometries and
topologies, these systems may exhibit distinct structures.

For simplicity, we represent the measurement process associated with multiple
RISs as

$\mathbf{S}=\mathbf{H}\mathbf{E}+\mathbf{n}.$ (14)

The localization problem can be formulated as follows: given the observation
$\mathbf{S}$ and the sensing operator $\mathbf{H}$, find a suitable RoI
$\mathbf{E}$ that satisfies (14).

Similar to DoA estimation via a single RIS, if the length of the observation
vector $\mathbf{S}$ is greater than the degree of RoI $\mathbf{E}$, and the
sensing operator $\mathbf{H}$ is of full column rank and well-conditioned, the
most straightforward estimation of the RoI is the least squares solution to
(14), which is expressed as

$\hat{\mathbf{E}}=\mathbf{H}^{\dagger}\mathbf{S}=\bigl{(}\mathbf{H}^{*}\mathbf{H}\bigr{)}^{-1}\mathbf{H}^{*}\mathbf{S}.$
(15)

######  Remark 6.

In the proposed RIS-centric backward sensing approach, we model the
measurement process in the far-field regime using a system of linear
equations. Various reconstruction algorithms can provide high-fidelity
solutions, but exploring these algorithms exceedes the scope of this paper.
The advantages of the least squares approach lie in its simplicity, broad
applicability, and strong theoretical foundation. There exist well-established
tools for analyzing the properties of least squares solutions.

##  IV Key Performance Indicators

This section is to identify the key indicators that dominate the performance
of the proposed RIS-centric sensing approaches. It’s important to note that
the indicators affecting sensing performance are algorithm-dependent. This
paper focuses on utilizing least squares solutions, as illustrated in (11) and
(15). If the measurement/sensing operator is a well-conditioned full column
rank matrix, the least squares solution provides a reliable estimation for the
incident field. Therefore, we examine sensing performance from both the rank
and conditioning perspectives.

###  IV-A Rank of the Sensing Operators

Suppose the sensing system comprises $K$ RISs, with each RIS composed of
$N_{k}$ elements. These RISs are collaboratively configured to sense an RoI
with dimension $M$. In the mode where each RIS has its dedicated receiver, the
number of snapshots for each RIS is denoted by $T_{k}$. On the other hand, in
the setup where the sensing system features only one receiver, the number of
snapshots is denoted by $T$. The following theorem clarifies the rank of the
sensing operators.

######  Theorem 7.

In the mode where each RIS has its dedicated receiver, the rank of the sensing
matrix is bounded above by the dimension of the RoI, the total number of
measurements, and the total number of elements. Specifically, we have

$\text{rank}(\mathbf{H})\leq\min\left\\{M,\sum_{k=1}^{K}\min\\{T_{k},N_{k}\\}\right\\}.$

Furthermore, in the scenario where the sensing system utilizes only one
receiver, the rank of the sensing matrix is bounded by

$\text{rank}(\mathbf{H})\leq\min\left\\{M,T,\sum_{k=1}^{K}N_{k}\right\\}.$

###### Proof.

We first establish the proof for the mode where each RIS has its dedicated
receiver. Since $\mathbf{H}$ consists of $M$ columns,
$\text{rank}(\mathbf{H})\leq M$. For each $k$, we note that

$\text{rank}\left(\mathbf{H}(\mathbf{\Omega}_{k})\text{diag}(l(\mathbf{r}^{\text{i}}_{k}))\right)=\text{rank}\left(\mathbf{H}(\mathbf{\Omega}_{k})\right)\leq\min\\{T_{k},N_{k}\\}.$

The last inequality arises from the fact that
$\mathbf{H}(\mathbf{\Omega}_{k})$ is a product of three matrices, as
illustrated in (8). The rank of the first matrix
$\mathbf{e}^{j\mathbf{\Omega}_{k}}$ is bounded by $T_{k}$, and the rank of the
second matrix is $N_{k}$. Because $\mathbf{H}$ is the stacking of $K$ such
matrices, we have

$\text{rank}(\mathbf{H})\leq\sum_{k=1}^{K}\min\\{T_{k},N_{k}\\}.$

The proof is similar for the mode where the sensing system utilizes only one
receiver. ∎

######  Remark 8.

In an RIS-centric backward sensing scheme, the rank of the sensing matrix
indicates how many degrees of freedom of the scene can be measured or coded by
this operator. Theorem 7 specifies that, to achieve high-fidelity recovery of
the RoI, a necessary condition is that the dimension of the RoI is less than
both the total number of elements and the number of measurements. This finding
has significant implication for the understanding that when the size of RISs
is predetermined, the dimension of RoI is bounded by this size, even if the
number of measurements can be greatly increased.

###  IV-B The Spatial Resolution Analysis

We now investigate the resolution associated with least squares solutions. For
clarity, we focus solely on exploring the angular resolution corresponding to
the uniform linear RIS. In fact, the relative error of a least squares
solution is given by

$\displaystyle\frac{\lVert\hat{\mathbf{E}}-\mathbf{E}\rVert}{\lVert\mathbf{E}\rVert}=$
$\displaystyle\frac{\lVert{\mathbf{H}}^{{\dagger}}\mathbf{H}\mathbf{E}-\mathbf{E}+{\mathbf{H}}^{{\dagger}}\mathbf{n}\rVert}{\lVert\mathbf{E}\rVert}$
(16) $\displaystyle\leq$
$\displaystyle\frac{\lVert{\mathbf{H}}^{{\dagger}}\mathbf{H}\mathbf{E}-\mathbf{E}\rVert+\lVert{\mathbf{H}}^{{\dagger}}\mathbf{n}\rVert}{\lVert\mathbf{E}\rVert}$
$\displaystyle\leq$
$\displaystyle\lVert{\mathbf{H}}^{{\dagger}}\mathbf{H}-\mathbf{I}\rVert+\lVert{\mathbf{H}}^{{\dagger}}\rVert\lVert\mathbf{n}\rVert/\lVert\mathbf{E}\rVert$
$\displaystyle=$
$\displaystyle\lVert{\mathbf{H}}^{{\dagger}}\mathbf{H}-\mathbf{I}\rVert+\sqrt{T}\lVert{\mathbf{H}}^{{\dagger}}\rVert/\sqrt{\mathrm{SNR}}.$

For simplicity, we define
$\mathrm{SNR}=\lVert\mathbf{E}\rVert^{2}/\sigma^{2}$.

The relative error is bounded above by the summation of two terms. When the
measurement operator $\mathbf{H}$ is a full column rank matrix, achieving
$\lVert{\mathbf{H}}^{{\dagger}}\mathbf{H}-\mathbf{I}\rVert\approx 0$ is
feasible. Furthermore, if $\lVert{\mathbf{H}}^{{\dagger}}\rVert$ remains
bounded by a reasonable value, then recovery using the least squares solution
ensures high fidelity. It should be noted that
$\lVert{\mathbf{H}}^{{\dagger}}\rVert=1/\sigma_{\text{min}}(\mathbf{H})$. We
will then examine the minimum singular value of $\mathbf{H}$.

In a sensing system, the minimum resolvable distance is defined as the spatial
resolution. The underlying motivation of this definition is that when two
point sources are in close proximity, the system cannot distinguish them as
two distinct points. Consider a scenario where only two incident waves, one
along $\theta^{\text{i}}$ and another along $\theta^{\text{i}}+\Delta$, are
impinging on the RIS. According to (9), the measurement operator is reduced to

$\mathbf{H}(\theta^{\text{i}},\theta^{\text{i}}+\Delta)=\tau
l(r^{\text{s}})\mathbf{e}^{j\mathbf{\Omega}}\text{diag}\left(\mathbf{v}(\theta^{\text{s}})\right)\mathbf{V}(\theta^{\text{i}},\theta^{\text{i}}+\Delta).$

When $\Delta\neq 0$, the matrix
$\mathbf{V}(\theta^{\text{i}},\theta^{\text{i}}+\Delta)$ is of full column
rank. And the matrix $\text{diag}\left(\mathbf{v}(\theta^{\text{s}})\right)$
is also of full rank. Therefore, we have

$\displaystyle\sigma_{\text{min}}(\mathbf{H}(\theta^{\text{i}},\theta^{\text{i}}+\Delta))$
(17) $\displaystyle\geq$ $\displaystyle\lvert\tau
l(r^{\text{s}})\rvert\sigma_{\text{min}}(\mathbf{e}^{j\mathbf{\Omega}})\sigma_{\text{min}}\left(\text{diag}(\mathbf{v}(\theta^{\text{s}}))\right)\sigma_{\text{min}}\left(\mathbf{V}(\theta^{\text{i}},\theta^{\text{i}}+\Delta)\right)$
$\displaystyle=$ $\displaystyle\lvert\tau
l(r^{\text{s}})\rvert\sigma_{\text{min}}(\mathbf{e}^{j\mathbf{\Omega}})\sigma_{\text{min}}\left(\mathbf{V}(\theta^{\text{i}},\theta^{\text{i}}+\Delta)\right).$

This inequality is a consequence of the lower bound for the smallest singular
value of a product of two matrices, as illustrated in Lemma 9.

######  Lemma 9.

If $\mathbf{B}$ has full column rank, then
$\sigma_{\text{min}}(\mathbf{A}\mathbf{B})\geq\sigma_{\text{min}}(\mathbf{A})\sigma_{\text{min}}(\mathbf{B})$.

###### Proof.

Given that $\mathbf{B}$ has full column rank, for any $x\neq 0$, we have
$\mathbf{B}x\neq 0$. Then,

$\sigma_{\text{min}}(\mathbf{A}\mathbf{B})=\min_{x\neq
0}\frac{\lVert\mathbf{A}\mathbf{B}x\rVert}{\lVert x\rVert}=\min_{x\neq
0}\frac{\lVert\mathbf{A}\mathbf{B}x\rVert}{\lVert\mathbf{B}x\rVert}\frac{\lVert\mathbf{B}x\rVert}{\lVert
x\rVert}\\\ \geq\min_{y\neq 0}\frac{\lVert\mathbf{A}y\rVert}{\lVert
y\rVert}\min_{x\neq 0}\frac{\lVert\mathbf{B}x\rVert}{\lVert
x\rVert}=\sigma_{\text{min}}(\mathbf{A})\sigma_{\text{min}}(\mathbf{B}).$

∎

######  Lemma 10.

Consider a Vandermonde matrix given by

$\mathbf{V}(\theta^{\mathrm{i}},\theta^{\mathrm{i}}+\Delta)=\begin{bmatrix}1&1\\\
e^{j2\pi d\sin\theta^{\mathrm{i}}/\lambda}&e^{j2\pi
d\sin(\theta^{\mathrm{i}}+\Delta)/\lambda}\\\ \vdots&\vdots\\\
e^{j2\pi(N-1)d\sin\theta^{\mathrm{i}}/\lambda}&e^{j2\pi(N-1)d\sin(\theta^{\mathrm{i}}+\Delta)/\lambda}\end{bmatrix},$
(18)

where $\Delta\neq 0$. The singular values of
$\mathbf{V}(\theta^{\mathrm{i}},\theta^{\mathrm{i}}+\Delta)$ are given by

$\sigma_{\text{max}}(\mathbf{V}(\theta^{\mathrm{i}},\theta^{\mathrm{i}}+\Delta))\\\
=\sqrt{N+\left|\frac{\sin(\pi
Nd(\sin\theta^{\mathrm{i}}-\sin(\theta^{\mathrm{i}}+\Delta))/\lambda)}{\sin(\pi
d(\sin\theta^{\mathrm{i}}-\sin(\theta^{\mathrm{i}}+\Delta))/\lambda)}\right|},$

and

$\sigma_{\text{min}}(\mathbf{V}(\theta^{\mathrm{i}},\theta^{\mathrm{i}}+\Delta))\\\
=\sqrt{N-\left|\frac{\sin(\pi
Nd(\sin\theta^{\mathrm{i}}-\sin(\theta^{\mathrm{i}}+\Delta))/\lambda)}{\sin(\pi
d(\sin\theta^{\mathrm{i}}-\sin(\theta^{\mathrm{i}}+\Delta))/\lambda)}\right|}.$
(19)

###### Proof.

See the Appendix A ∎

######  Lemma 11.

For $x$ sufficiently small, we have

$\frac{\sin(Nx)}{\sin(x)}\approx N-\frac{N(N^{2}-1)x^{2}}{6}.$ (20)

###### Proof.

See the Appendix B ∎

Combining (19) and (20) yields an approximation for the smallest singular
value of $\mathbf{V}(\theta^{\text{i}},\theta^{\text{i}}+\Delta)$

$\displaystyle\sigma_{\text{min}}(\mathbf{V}(\theta^{\text{i}},\theta^{\text{i}}+\Delta))$
(21) $\displaystyle\approx$
$\displaystyle\frac{\pi}{\sqrt{6}}\frac{d}{\lambda}\sqrt{N(N^{2}-1)}\lvert\sin\theta^{\text{i}}-\sin(\theta^{\text{i}}+\Delta)\rvert$
$\displaystyle\approx$
$\displaystyle\frac{\pi}{\sqrt{6}}\frac{d}{\lambda}\sqrt{N(N^{2}-1)}\lvert\Delta\rvert\cos\theta^{\text{i}}.$

In the last step, we utilize the first-order Taylor expansion of
$\sin(\theta^{\text{i}}+\Delta)$ at $\theta^{\text{i}}$.

Next, we will examine $\sigma_{\text{min}}(\mathbf{e}^{j\mathbf{\Omega}})$
using the Marchenko-Pastur law [43], which characterizes the asymptotic
behavior of singular values of large rectangular random matrices. Let’s
consider a random matrix $\mathbf{X}$ of size $T\times N$ $(T>N)$, whose
entries are independent identically distributed random variables with mean 0
and variance 1. We define the empirical spectral distribution of
$\mathbf{Y}_{N}=\mathbf{X}^{*}\mathbf{X}/T$ as

$\mu_{\mathbf{Y}_{N}}\equiv\frac{1}{N}\sum_{i=1}^{N}\delta_{\lambda_{i}(\mathbf{Y}_{N})},$

where $\lambda_{1}(\mathbf{Y}_{N})\leq\cdots\leq\lambda_{N}(\mathbf{Y}_{N})$
are the eigenvalues (including multiplicity) and
$\delta_{\lambda_{i}(\mathbf{Y}_{N})}(x)$ is the indicator function
$\mathbf{1}_{\lambda_{i}(\mathbf{Y}_{N})\leq x}$.

A celebrated theorem in random matrix theory, known as the Marchenko-Pastur
distribution, describes the limiting spectral distribution of $\mathbf{Y}_{N}$
[43, 44]. Under the Kolmogorov condition, i.e., $T,N\to\infty$ and $N/T\to
c\in(0,1)$, $\mu_{\mathbf{Y}_{N}}$ converges weakly and almost surely to
$\mu_{\text{M-P}}$ with a density function

$d\mu_{\text{M-P}}(x)=\frac{1}{2\pi cx}\sqrt{(x-a)^{+}(b-x)^{+}}dx.$

Here $a=(1-\sqrt{c})^{2}$, $b=(1+\sqrt{c})^{2}$, and $a^{+}=\max(0,a)$. The
Marchenko-Pastur distribution implies, as $T,N\to\infty$, the probability of
the eigenvalues of $\mathbf{Y}_{N}$ lying outside the interval $[a,b]$
approaches 0.

We now turn our attention to the matrix $\mathbf{e}^{j\mathbf{\Omega}}$, whose
entries are independent identically distributed random variables with mean 0
and variance 1. By leveraging the Marchenko-Pastur law, we can obtain a good
approximation for its smallest singular value in the regime of large $T$ and
$N$,

$\sigma_{\text{min}}(\mathbf{e}^{j\mathbf{\Omega}})\approx\sqrt{T}\left(1-\sqrt{\frac{N}{T}}\right)=\sqrt{T}-\sqrt{N}.$
(22)

Substituting (21) and (22) into (17) yields

$\displaystyle\sigma_{\text{min}}(\mathbf{H}(\theta^{\text{i}},\theta^{\text{i}}+\Delta))$
$\displaystyle\gtrsim$
$\displaystyle\frac{\pi}{\sqrt{6}}\tau(r^{\text{s}})^{-1}\frac{d}{\lambda}(\sqrt{T}-\sqrt{N})\sqrt{N(N^{2}-1)}\lvert\Delta\rvert\cos\theta^{\text{i}}$
$\displaystyle\approx$
$\displaystyle\frac{\pi}{\sqrt{6}}\tau(r^{\text{s}})^{-1}\frac{d}{\lambda}\left(1-\sqrt{\frac{N}{T}}\right)T^{\frac{1}{2}}N^{\frac{3}{2}}\lvert\Delta\rvert\cos\theta^{\text{i}}.$

It then follows that

$\left\lVert{\mathbf{H}(\theta^{\text{i}},\theta^{\text{i}}+\Delta)}^{+}\right\rVert\lesssim\frac{\sqrt{6}}{\pi}r^{\text{s}}\tau^{-1}\left(\frac{d}{\lambda}\right)^{-1}\left(1-\sqrt{\frac{N}{T}}\right)^{-1}\\\
T^{-\frac{1}{2}}N^{-\frac{3}{2}}\lvert\Delta\rvert^{-1}(\cos\theta^{\text{i}})^{-1}.$

We finally obtain the main theorem, as detailed below.

######  Theorem 12.

If only two closely incident waves along $\theta^{\mathrm{i}}$ and
$\theta^{\mathrm{i}}+\Delta$ are impinging on a uniform linear RIS, the
relative error for the least squares solution is bounded approximately

$\frac{\lVert\hat{\mathbf{E}}-\mathbf{E}\rVert}{\lVert\mathbf{E}\rVert}\lesssim\frac{\sqrt{6}}{\pi}r^{\mathrm{s}}\tau^{-1}\left(\frac{d}{\lambda}\right)^{-1}\left(1-\sqrt{\frac{N}{T}}\right)^{-1}\\\
N^{-\frac{3}{2}}\lvert\Delta\rvert^{-1}(\cos\theta^{\mathrm{i}})^{-1}\mathrm{SNR}^{-1/2}.$
(23)

######  Remark 13.

Theorem 12 clarifies the factors that affect the sensing performance of a
single uniform linear RIS. As seen from (23), in addition to the SNR, three
distinct categories of factors play significant roles in determining the
relative error. The first category relates to the topology and geometry of the
RIS, including parameters such as $\tau$ (the scattering pattern of the
element), $N$ (the number of elements) and $d$ (the spacing between elements).
By rewriting $(d/\lambda)^{-1}N^{-\frac{3}{2}}$ as
$(Nd/\lambda)^{-1}N^{-\frac{1}{2}}$, we observe that the bound is inversely
proportional to $Nd$, representing the aperture of the RIS. Unlike the
situation in forward beamforming, where undesirable grating lobes arise when
$d/\lambda>1/2$, there is no such constraint on the spacing $d$ in backward
sensing. Indeed, increasing $d$ can enhance the aperture, leading to greater
gains. However, it’s essential to note that increasing $N$ (the number of
elements) can potentially provide more gains. This observation can be
attributed to the fact that the power reflected by RISs grows quadratically
with the total collecting area of the elements, as clarified in Section V.

######  Remark 14.

The second category that affects sensing performance is related to incident
angles and the receiver’s position. As the receiver moves closer to the RIS,
it captures more power, resulting in lower recovery errors. This relationship
is detailed in Eq. (23), where the recovery error is directly proportional to
$r^{\text{s}}$. Moreover, the relative error is inversely proportional to
$\lvert\Delta\rvert$, providing insight into the spatial resolution
performance. Mathematically, for small values of $\lvert\Delta\rvert$, the
matrix $\mathbf{V}(\theta^{\text{i}},\theta^{\text{i}}+\Delta)$ becomes ill-
conditioned, which in turn leads to ill-conditioning of the measurement matrix
$\mathbf{H}(\theta^{\text{i}},\theta^{\text{i}}+\Delta)$. Consequently, even a
small noise in the measurement procedure can result in much larger errors in
the recovery. Additionally, the recovery performance worsens as the scanned
range moves away from the normal direction of the RIS.

######  Remark 15.

The third category that affects sensing performance is related to the number
of measurements. The relationship between the bound and the number of
measurements $T$ is described by the term $(1-\sqrt{N/T})^{-1}$. This
indicates that the bound is inversely proportional to $1-\sqrt{N/T}$. When
$T\gg N$, this term asymptotically approaches 1, its minimum value. As a rule
of thumb, increasing the number of measurements well beyond the number of
elements can potentially enhance the performance of backward sensing.

##  V Numerical Validation

Theorems 7 and 12 provide theoretical insights into the key parameters that
affect the performance of the proposed RIS-centric sensing approaches. In this
section, numerical validations are presented. Specifically, our goal is to
sense a region where multiple power sources of various shapes have been
manually placed. We comprehensively analyze the impact of several factors,
including the total number of samples, the total number of elements, the
number of RISs and their positions.

To evaluate the quality of backward sensing, we use two metrics: relative
error and structural similarity index measure (SSIM). The SSIM is a perceptual
metric to quantify the degradation of the recovered field $\hat{\mathbf{E}}$
compared to the ground truth $\mathbf{E}$, defined as

$\text{SSIM}=\frac{\left(2\mu_{\mathbf{E}}\mu_{\hat{\mathbf{E}}}+c_{1}\right)\left(2\sigma_{\mathbf{E}\hat{\mathbf{E}}}+c_{2}\right)}{\left(\mu_{\mathbf{E}}^{2}+\mu_{\hat{\mathbf{E}}}^{2}+c_{1}\right)\left(\sigma_{\mathbf{E}}^{2}+\sigma_{\hat{\mathbf{E}}}^{2}+c_{2}\right)}.$

Here $\mu_{\hat{\mathbf{E}}}$ and $\mu_{\mathbf{E}}$ are the means of
$\hat{\mathbf{E}}$ and $\mathbf{E}$, $\sigma^{2}_{\hat{\mathbf{E}}}$ and
$\sigma^{2}_{\mathbf{E}}$ are the variances, and
$\sigma_{\mathbf{E}\hat{\mathbf{E}}}$ is the covariance between $\mathbf{E}$
and $\hat{\mathbf{E}}$. The constants $c_{1}$ and $c_{2}$ are used to
stabilize the division with weak denominator. The SSIM value lies within the
range of $[0,1]$, where a higher value indicates that the recovered field
exhibits structures more similar to the ground truth.

Theorem 7 asserts that the rank of the sensing matrix is upper-bounded by the
total number of measurements and the total number of elements. To assess the
influence of these factors, we first examine the performance of backward
sensing in the noiseless case.

Figure 6: Scenario configuration used in the simulation. A 10-meter by 10-meter square region is discretized into a 10 $\times$ 10 grid with 1-meter intervals. Four uniform linear RISs are strategically positioned 25 meters around the RoI, with each adjacent RIS separated by an angle of 45∘. TABLE I: Default setup parameters used in the simulation experiments. Parameter | Symbol | Value
---|---|---
RoI | $x\times y$ | 10 m $\times$ 10 m
Spatial resolution (Pixel size) | $\bigtriangleup x\times\bigtriangleup y$ | 1 m $\times$ 1 m
Pixel dimension | $M_{x}\times M_{y}$ | 10 $\times$ 10
Number of RISs | K | 4
Operating frequency | $f$ |  $30$ GHz
RIS element distance | $d$ |  $0.01$ m
Distance from RoI to RIS | $L$ | 25 m

Let’s imagine a scenario where the task is to sense a square region measuring
$10$ meters by $10$ meters. Four uniform linear RISs are strategically
positioned around the RoI. Each RIS is equipped with a dedicated receiver. The
topology and geometry are demonstrated in Fig. 6. The measurement vectors
corresponding to each RIS are stacked into a column vector, as depicted in
(12). The system operates at a frequency of $f=30$ GHz. The RoI is discretized
into $10\times 10$ pixels. The parameters are detailed in Table I.

For a comprehensive illustration, we systematically record the recovered
results as the number of measurements $T_{k}$ for each RIS varies at 10, 20,
50, and 100, with the number of elements fixed at $N_{k}=50$ for each $k$. The
total number of elements sums up to $N=\sum_{k=1}^{4}N_{k}=200$. The
reconstructed scenes are illustrated in Fig.7. We plot the relative error and
SSIM values in Fig. 8. Additionally, we increase the number of elements of
each RIS from 20 to 50, while keeping the number of measurements fixed at
$T_{k}=50$ (yielding a total of $T=\sum_{k=1}^{4}T_{k}=200$ measurements). The
results are presented in Figs. 9 and 10.

Figure 7: Comparison of the recovered scenes as the number of measurements
$T_{k}$ for each RIS varied at 10, 20, 50, and 100, with the number of
elements fixed at $N_{k}=50$ for each $k$.

Figure 8: Relative error and SSIM as the number of measurements $T_{k}$ for
each RIS varied at 10, 20, 50, and 100, with the number of elements fixed at
$N_{k}=50$ for each $k$.

Figure 9: Comparison of the recovered scenes as the number of elements $N_{k}$
for each RIS varied at 20, 30, 40, and 50, with the number of measurements
fixed at $T_{k}=50$ for each $k$.

Figure 10: Relative error and SSIM as the number of elements $N_{k}$ for each
RIS varied at 20, 30, 40, and 50, with the number of measurements fixed at
$T_{k}=50$ for each $k$.

These plots demonstrate that increasing both the number of measurements and
the number of elements enhances the performance of backward sensing. Thess
results are consistent with intuitive expectations. When $N$ or $T$ is not
sufficiently large, a fundamental limitation arises: the sensing matrix
encodes the high-dimensional scene into a low-dimensional observation.
Clearly, this dimensionality reduction makes recovering the original scene
using linear operators impractical. For instance, with $T_{k}=10$, accurately
reconstructing an original vector in $\mathbb{C}^{100}$ from an observation
lying in $\mathbb{C}^{40}$ becomes challenging.

It’s worth emphasizing that even if the sensing matrix encodes the scene
vector into a high-dimensional observation, achieving high-precision recovery
might still be challenging. The presence of noise, particularly in cases of
ill-conditioned sensing matrices, makes high-fidelity reconstruction from
observations unfeasible. For better understanding, we illustrate the singular
values corresponding to all the aforementioned trials in Fig. 11.

Figure 11: Sorted singular values corresponding to all trials. (a) The number
of measurements $T_{k}$ is varied at 10, 20, 50, and 100, with the number of
elements fixed at $N_{k}=50$. (b) $K$ The number of elements $N_{k}$ is varied
at 20, 30, 40, and 50, with the number of measurements fixed at $T_{k}=50$.

In the context of backward sensing using multiple RISs, the conditioning of
the sensing matrix depends not only on the number of elements and measurements
but also on factors such as the system’s geometry and topology. Optimizing the
topology of RISs has the potential to enhance the conditioning of the sensing
matrix, thereby improving the robustness of scene reconstruction.

To demonstrate this improvement, we conduct a series of experiments by
strategically selecting eight landmark positions around the RoI, each of which
serves as a potential location for an RIS. These RISs are positioned
tangentially to a circle encompassing the RoI. The location of these landmarks
are shown in Fig 12. The specific positions are detailed in Table II. To
ensure a fair comparison, we maintain the total number of elements at
$N=\sum_{k}N_{k}=200$ and the total number of measurements at
$T=\sum_{k}T_{k}=200$. To assess the stability of the sensing system, additive
Gaussian noise is introduced in these experiments, with the SNR set to 30 dB.

TABLE II: Deployment strategies for combinations of different candidate positions. Category | Configuration | Description
---|---|---
Strategy I: | {A} |  -Use only a single landmark position
Strategy II: | {A, E} |  -Use two landmark positions perpendi- cular to the RoI
Strategy III: | {A, C, E, G} |  -Use four landmark positions with each adjacent interval of them at an angle of 45∘
Strategy IV: |  {A, B, C, D, E, F, G, H} |  -Use all landmark positions with high landmark density

The results of the backward sensing experiments are presented in Fig.13. In
addition to the reconstructed scenes, we also plot the sorted singular values
and the condition numbers of the sensing matrices in Fig.14. For better
illustration, the y-axis uses a logarithmic scale. It is evident that
recovering the scene with acceptable fidelity is impossible when using only 1
or 2 RISs in this specific scenario. Notably, when using 2 RISs, the condition
number of the sensing matrix is greater than $3.1\times 10^{5}$. The backward
sensing procedure is highly sensitive to perturbations in the observations.

Figure 12: Scenario configuration and landmark positions. Eight landmark
positions are set 25 meters around the RoI, labeled from A to H. These
landmarks are spaced at 22.5∘ intervals relative to the RoI, serving as
candidate locations for RIS deployment.

Figure 13: Comparison of the recovered scenes using various geometries and
topologies of RISs. Achieving acceptable fidelity in scene recovery is
impossible when employing only 1 or 2 RISs in this specific scenario.

Figure 14: (a) Sorted singular values for each deployment strategy, with the
y-axis presented on a logarithmic scale for clarity. (b) Condition numbers of
the measurement matrix for each deployment strategy.

Figure 15: Photographs of the front (a) and back (b) of the fabricated RIS
prototype with 10$\times$ 16 elements. The RIS, operating at a frequency of
5.8 GHz, consist of 10 $\times$ 16 1-bit elements with an element spacing of
0.025 m.

Figure 16: Experimental setup of two RISs sensing system in a microwave
anechoic chamber.

##  VI Proof-of-concept Prototype

To validate the proposed sensing scheme, we implement a proof-of-concept
prototype using a universal software radio peripheral (USRP) in a microwave
anechoic chamber. Two RISs are strategically positioned at one end of the
chamber to detect a power source located at the other end, as illustrated in
Fig. 16. The schematic is depicted as Fig. 17. The RISs, operating at a
frequency of 5.8 GHz, consist of 10 $\times$ 16 1-bit elements with an element
spacing of 0.025 m, as shown in Fig. 15. To evaluate the performance of
uniform linear RISs within a planar RoI, each column is configured
identically. The two RISs, separated by 2.81 m, are employed to detect the
source positioned 6 m away. For each RIS, $T=500$ measurements are acquired to
construct the sensing matrix $\mathbf{H}$.

Figure 17: The geometry and topology of the experiment conducted inside the
anechoic chamber. The RoI is a square with dimensions of 12 m $\times$ 12 m,
discretized into a grid of 6 $\times$ 6, resulting in each pixel being 2 m
$\times$ 2 m.

Figure 18: Results of DoA derived from both simulated and measured data by the
magnitude-only reconstruction algorithm. The first row displays the results of
DoA from the left RIS. The second row displays the results of DoA from the
right RIS. Results demonstrate high accuracy with errors within $2^{\circ}$.

Figure 19: Recovered results using the magnitude-only reconstruction algorithm
in the chamber. The left figure displays the reconstructed scene from
simulated data. The right figure shows the reconstructed scene from measured
data.

It’s important to note that the validation using USRP has a limitation: while
the strength of the carrier wave is easy to measure, accurately capturing the
phase is more challenging. Therefore, the backward sensing process must rely
on magnitude-only reconstruction algorithms. Mathematically, the magnitude-
only backward sensing problem is formulated as follows: given $\mathbf{H}$ and
$\lvert\mathbf{S}\rvert$, find $\mathbf{E}$ such that

$\lvert\mathbf{S}\rvert\approx\lvert\mathbf{H}\mathbf{E}\rvert.$

In contrast to phased measurements described by (14), the missing phase makes
the problem more complex. The phaseless reconstruction problem is highly
nonlinear and significantly more ill-posed than the phased reconstruction
problem [45]. Various approaches have been developed to address the phaseless
reconstruction problem, including gradient descent methods, iterative
projection methods, and convex optimization-based methods [46]. We focuse on
validating our backward sensing approach rather than delving into the
specifics of these algorithms. For our purposes, we employ the reweighted
Wirtinger flow (RWF) algorithm as the solver [47].

We initiate the validation of the RWF algorithm by estimating the DoA using
data independently obtained from each of the two RISs. Fig. 18 shows the DoA
estimations from both RISs. Sepcifically, a single source is localized at
$-8.5^{\circ}$ for the left RIS and $-13.25^{\circ}$ for the right RIS. In
numerical experiments, a distinct peak is observed with DoA estimations of
$-8.5^{\circ}$ and $-13^{\circ}$ for the left RIS and the right RIS,
respectively. In the results estimated using measured data, the DoA
estimations are $-7.5^{\circ}$ for the left RIS and $-14.5^{\circ}$ for the
right RIS, with errors within two degrees. These figures illustrate the
remarkably accurate estimation of incident angles, thereby demonstrating the
accuracy of our modeling and the effectiveness of the RWF algorithm within
this experimental setup.

To further validate the collaborative sensing capabilities of the two RISs, we
define the RoI as a square with dimensions of 12 m $\times$ 12 m, discretized
into a grid of 6 $\times$ 6, resulting in each pixel being 2 m $\times$ 2 m.
This choice of resolution takes into account the limited separation between
the two RISs and the relatively small size of each RIS, which comprises only
16 elements. The final results, illustrated in Fig. 19, provide a
comprehensive overview. The true location of the source is indicated by red
pentagrams. The plots clearly demonstrate that the power source can be
precisely located using the two RISs. These findings validate the feasibility
of the proposed backward sensing approach with multiple RISs.

$\displaystyle{\mathbf{V}}^{*}\mathbf{V}$
$\displaystyle=\begin{bmatrix}1&e^{-j2\pi
d\sin\theta^{\text{i}}/\lambda}&\cdots&e^{-j2\pi(N-1)d\sin\theta^{\text{i}}/\lambda}\\\
1&e^{-j2\pi
d\sin(\theta^{\text{i}}+\Delta)/\lambda}&\cdots&e^{-j2\pi(N-1)d\sin(\theta^{\text{i}}+\Delta)/\lambda}\\\
\end{bmatrix}\begin{bmatrix}1&1\\\ e^{j2\pi
d\sin\theta^{\text{i}}/\lambda}&e^{j2\pi
d\sin(\theta^{\text{i}}+\Delta)/\lambda}\\\ \vdots&\vdots\\\
e^{j2\pi(N-1)d\sin\theta^{\text{i}}/\lambda}&e^{j2\pi(N-1)d\sin(\theta^{\text{i}}+\Delta)/\lambda}\end{bmatrix}$
(24)
$\displaystyle=\begin{bmatrix}N&1+\cdots+e^{j2\pi(N-1)d(\sin(\theta^{\text{i}}+\Delta)-\sin\theta^{\text{i}})/\lambda}\\\
1+\cdots+e^{j2\pi(N-1)d(\sin\theta^{\text{i}}-\sin(\theta^{\text{i}}+\Delta))/\lambda}&N\end{bmatrix}.$

$\displaystyle\lambda_{1},\lambda_{2}$
$\displaystyle=N\pm\sqrt{(1+\dots+e^{j2\pi(N-1)d(\sin\theta^{\text{i}}-\sin(\theta^{\text{i}}+\Delta))/\lambda})(1+\dots+e^{j2\pi(N-1)d(\sin(\theta^{\text{i}}+\Delta)-\sin\theta^{\text{i}})/\lambda})}$
(25) $\displaystyle=N\pm\left|\frac{1-e^{j2\pi
Nd(\sin\theta^{\text{i}}-\sin(\theta^{\text{i}}+\Delta))\lambda}}{1-e^{j2\pi
d(\sin\theta^{\text{i}}-\sin(\theta^{\text{i}}+\Delta))/\lambda}}\right|,$
$\displaystyle=N\pm\left|\frac{\sin(\pi
Nd(\sin\theta^{\text{i}}-\sin(\theta^{\text{i}}+\Delta))/\lambda)}{\sin(\pi
d(\sin\theta^{\text{i}}-\sin(\theta^{\text{i}}+\Delta))/\lambda)}\right|.$

##  VII Conclusion

This paper investigates the capability and practicality of the proposed RIS-
centric sensing approach. We present physically accurate and mathematically
concise models to characterize the reflection properties of RISs. Using
linearized approximations inherent in the far-field region, practical models
are derived to describe the forward signal aggregation process via RISs. We
demonstrate that DoA estimation of incident waves can be achieved with a
single RIS. Additionally, we show that the spatial diversity provided by
deploying multiple RISs allows for accurate localization of multiple power
sources. A theoretical framework for determining key performance indicators is
established through condition number analysis of the sensing operators.
Numerical experiments validate our findings. To showcase the practicality of
our proposed RIS-centric sensing approach, we develop a proof-of-concept using
USRP. To the best of our knowledge, no previous prototypes have been proposed
with such configurations.

##  Appendix A Proof of Lemma 10

###### Proof.

For the matrix $\mathbf{V}$ defined by (18), we first calculate
${\mathbf{V}}^{*}\mathbf{V}$, as shown in (24). It is straightforward to
calculate the eigenvalues of ${\mathbf{V}}^{*}\mathbf{V}$, which are given in
(25). Finally, we establish the singular values of $\mathbf{V}$. ∎

##  Appendix B Proof of Lemma 11

###### Proof.

The main idea of the proof is to find the series expansion

$\frac{\sin Nx}{\sin x}=\sum_{k=0}^{\infty}a_{k}x^{2k}.$ (26)

It is easy to verify that

$\displaystyle\frac{\sin Nx}{\sin x}$
$\displaystyle=\frac{e^{iNx}-e^{-iNx}}{e^{ix}-e^{-ix}}=\sum_{l=0}^{N-1}e^{i(N-1-2l)x}$
(27)
$\displaystyle=\left\\{\begin{array}[]{ll}1+2\sum_{l=1}^{\frac{N-1}{2}}\cos
2lx&N\text{ odd},\\\ 2\sum_{l=1}^{\frac{N}{2}}\cos(2l-1)x&N\text{
even}.\end{array}\right.$

Using the Taylor expansion $\cos
x=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(2k)!}x^{2k}$, we find

$\cos 2lx=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(2k)!}(2l)^{k}x^{2k},$ (28)

and

$\cos(2l-1)x=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(2k)!}(2l-1)^{k}x^{2k}.$ (29)

Substituting (28) and (29) into (27) yields

$a_{k}=\frac{(-1)^{k}}{(2k)!}\sum_{l=0}^{N-1}(N-1-2l)^{2k}.$ (30)

It follows that $a_{0}=N$ and $a_{1}=-\frac{N(N^{2}-1)}{6}$. Finally, for
sufficiently small $x$, we have

$\frac{\sin(Nx)}{\sin x}\approx N-\frac{N(N^{2}-1)x^{2}}{6}.$ (31)

∎

## References

  * [1] J. Du, Y. Cheng, L. Jin, S. Li, and F. Gao, “Nested tensor-based integrated sensing and communication in RIS-assisted THz MIMO systems,” _IEEE Trans. Signal Process._ , vol. 72, pp. 1141–1157, 2024.
  * [2] W. Tang et al., “MIMO transmission through reconfigurable intelligent surface: System design, analysis, and implementation,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2683–2699, Nov. 2020.
  * [3] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 11, pp. 5394–5409, Nov. 2019.
  * [4] E. Basar et al., “Wireless communications through reconfigurable intelligent surfaces,” _IEEE access_ , vol. 7, pp. 116 753–116 773, 2019.
  * [5] X. Song, J. Xu, F. Liu, T. X. Han, and Y. C. Eldar, “Intelligent reflecting surface enabled sensing: Cramér-Rao bound optimization,” _IEEE Trans. Signal Process._ , vol. 71, pp. 2011–2026, 2023.
  * [6] K. Zhong, J. Hu, C. Pan, M. Deng, and J. Fang, “Joint waveform and beamforming design for RIS-aided ISAC systems,” _IEEE Signal Process. Lett._ , vol. 30, pp. 165–169, 2023.
  * [7] X. Liu, Y. Liu, Y. Chen, and H. V. Poor, “RIS enhanced massive non-orthogonal multiple access networks: Deployment and passive beamforming design,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 4, pp. 1057–1071, Apr. 2021.
  * [8] R. Xiong, K. Yin, T. Mi, J. Lu, K. Wan, and R. C. Qiu, “Fair beam allocations through reconfigurable intelligent surfaces,” _IEEE J. Sel. Areas Commun._ , pp. 1–1, 2024.
  * [9] X. Ma, Y. Fang, H. Zhang, S. Guo, and D. Yuan, “Cooperative beamforming design for multiple RIS-assisted communication systems,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 12, pp. 10 949–10 963, Dec. 2022.
  * [10] Z. Chu, Z. Zhu, F. Zhou, M. Zhang, and N. Al-Dhahir, “Intelligent reflecting surface assisted wireless powered sensor networks for Internet of Things,” _IEEE Trans. Commun._ , vol. 69, no. 7, pp. 4877–4889, Jul. 2021.
  * [11] M. Ouyang, F. Gao, Y. Wang, S. Zhang, P. Li, and J. Ren, “Computer vision-aided reconfigurable intelligent surface-based beam tracking: Prototyping and experimental results,” _IEEE Trans. Wireless Commun._ , vol. 22, no. 12, pp. 8681–8693, Dec. 2023.
  * [12] Q. Cheng et al., “Reconfigurable intelligent surfaces: Simplified-architecture transmitters—from theory to implementations,” _Proc. IEEE_ , vol. 110, no. 9, pp. 1266–1289, Sept. 2022.
  * [13] Z. Ma et al., “Modeling and analysis of MIMO multipath channels with aerial intelligent reflecting surface,” _IEEE J. Sel. Areas Commun._ , vol. 40, no. 10, pp. 3027–3040, Oct. 2022.
  * [14] H. Zhang, H. Zhang, B. Di, K. Bian, Z. Han, and L. Song, “Metalocalization: Reconfigurable intelligent surface aided multi-user wireless indoor localization,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 12, pp. 7743–7757, Dec. 2021.
  * [15] M. Lin et al., “Single sensor to estimate DOA with programmable metasurface,” _IEEE Internet Things J._ , vol. 8, no. 12, pp. 10 187–10 197, Jun. 2021.
  * [16] Y. He, D. Zhang, and Y. Chen, “High-resolution WiFi imaging with reconfigurable intelligent surfaces,” _IEEE Internet Things J._ , vol. 10, no. 2, pp. 1775–1786, Jan. 2022.
  * [17] X. Shao, C. You, W. Ma, X. Chen, and R. Zhang, “Target sensing with intelligent reflecting surface: Architecture and performance,” _IEEE J. Sel. Areas Commun._ , vol. 40, no. 7, pp. 2070–2084, Jul. 2022.
  * [18] S. S. Moghadasian and S. Gazor, “Sparsely localized time-frequency energy distributions for multi-component LFM signals,” _IEEE Signal Process. Lett._ , vol. 27, pp. 6–10, 2020.
  * [19] S. Ahmed, D. Wang, J. Park, and S. H. Cho, “UWB-gestures, a public dataset of dynamic hand gestures acquired using impulse radar sensors,” _Sci. Data_ , vol. 8, no. 1, p. 102, Apr. 2021.
  * [20] B. Zhao, K. Hu, F. Wen, S. Cui, and Y. Shen, “TDLoc: Passive localization for MIMO-OFDM system via tensor decomposition,” _IEEE Internet Things J._ , vol. 10, no. 23, pp. 20 819–20 833, Dec. 2023.
  * [21] G. Lan, M. F. Imani, P. Del Hougne, W. Hu, D. R. Smith, and M. Gorlatova, “Wireless sensing using dynamic metasurface antennas: Challenges and opportunities,” _IEEE Commun. Mag._ , vol. 58, no. 6, pp. 66–71, Jun. 2020\.
  * [22] G. C. Alexandropoulos et al., “RIS-enabled smart wireless environments: Deployment scenarios, network architecture, bandwidth and area of influence,” _EURASIP J. Wirel. Commun. Netw._ , vol. 2023, no. 1, p. 103, Oct. 2023.
  * [23] H. Zhao et al., “Intelligent indoor metasurface robotics,” _Natl. Sci. Rev._ , vol. 10, no. 8, p. nwac266, Nov. 2023.
  * [24] Y. Jiang, F. Gao, Y. Liu, S. Jin, and T. Cui, “Near field computational imaging with RIS generated virtual masks,” _IEEE Trans. Antennas Propag._ , vol. 72, no. 5, pp. 4383–4398, May 2024.
  * [25] S. Buzzi, E. Grossi, M. Lops, and L. Venturino, “Radar target detection aided by reconfigurable intelligent surfaces,” _IEEE Signal Process. Lett._ , vol. 28, pp. 1315–1319, 2021.
  * [26] H. Wymeersch, J. He, B. Denis, A. Clemente, and M. Juntti, “Radio localization and mapping with reconfigurable intelligent surfaces: Challenges, opportunities, and research directions,” _IEEE Veh. Technol. Mag._ , vol. 15, no. 4, pp. 52–61, Dec. 2020.
  * [27] Z. Li, T. Jin, D. Guan, and H. Xu, “Metaphys: Contactless physiological sensing of multiple subjects using RIS-based 4-D radar,” _IEEE Internet Things J._ , vol. 10, no. 14, pp. 12 616–12 626, Jul. 2023.
  * [28] Y. Jiang, F. Gao, M. Jian, S. Zhang, and W. Zhang, “Reconfigurable intelligent surface for near field communications: Beamforming and sensing,” _IEEE Trans. Wireless Commun._ , vol. 22, no. 5, pp. 3447–3459, May 2022\.
  * [29] K. Keykhosravi et al., “Leveraging RIS-enabled smart signal propagation for solving infeasible localization problems: Scenarios, key research directions, and open challenges,” _IEEE Veh. Technol. Mag._ , vol. 18, no. 2, pp. 20–28, Jun. 2023.
  * [30] H. Wymeersch and G. Seco-Granados, “Radio localization and sensing—part II: State-of-the-art and challenges,” _IEEE Commun. Lett._ , vol. 26, no. 12, pp. 2821–2825, Dec. 2022.
  * [31] N. Mehrotra and A. Sabharwal, “When does multipath improve imaging resolution?” _IEEE J. Sel. Areas Inf. Theory_ , vol. 3, no. 1, pp. 135–146, Mar. 2022.
  * [32] S. Gezici, “A survey on wireless position estimation,” _Wirel. Pers. Commun._ , vol. 44, no. 3, pp. 263–282, Oct. 2008.
  * [33] I. Guvenc and C.-C. Chong, “A survey on TOA based wireless localization and NLOS mitigation techniques,” _IEEE Commun. Surv. Tutor._ , vol. 11, no. 3, pp. 107–124, 2009.
  * [34] O. Rinchi, A. Elzanaty, and M.-S. Alouini, “Compressive near-field localization for multipath RIS-aided environments,” _IEEE Commun. Lett._ , vol. 26, no. 6, pp. 1268–1272, Jun. 2022.
  * [35] Q. Y. Zhou et al., “Two-dimensional direction-of-arrival estimation based on time-domain-coding digital metasurface,” _Appl. Phys. Lett._ , vol. 121, no. 18, Oct. 2022.
  * [36] A. H. Sayed, A. Tarighat, and N. Khajehnouri, “Network-based wireless location: Challenges faced in developing techniques for accurate wireless location information,” _IEEE Signal Process. Mag._ , vol. 22, no. 4, pp. 24–40, Jul. 2005.
  * [37] G. C. Alexandropoulos, I. Vinieratou, and H. Wymeersch, “Localization via multiple reconfigurable intelligent surfaces equipped with single receive RF chains,” _IEEE Wireless Commun. Lett._ , vol. 11, no. 5, pp. 1072–1076, May 2022.
  * [38] J. Wilson and N. Patwari, “Radio tomographic imaging with wireless networks,” _IEEE Trans. Mob. Comput._ , vol. 9, no. 5, pp. 621–632, May 2010.
  * [39] W. C. Karl, J. E. Fowler, C. A. Bouman, M. Çetin, B. Wohlberg, and J. C. Ye, “The foundations of computational imaging: A signal processing perspective,” _IEEE Signal Process Mag._ , vol. 40, no. 5, pp. 40–53, Jul. 2023.
  * [40] T. Mi, J. Zhang, R. Xiong, Z. Wang, P. Zhang, and R. C. Qiu, “Towards analytical electromagnetic models for reconfigurable intelligent surfaces,” _IEEE Trans. Wireless Commun._ , vol. 23, no. 5, pp. 4170–4185, May 2024\.
  * [41] T. Sleasman, M. F. Imani, J. N. Gollub, and D. R. Smith, “Microwave imaging using a disordered cavity with a dynamically tunable impedance surface,” _Phys. Rev. Appl._ , vol. 6, no. 5, p. 054019, Nov. 2016.
  * [42] H. Maître, _Processing of synthetic aperture radar (SAR) images_. John Wiley & Sons, 2013.
  * [43] V. Marchenko and L. Pastur, “The distribution of eigenvalues in certain sets of random matrices math,” _USSR-st_ , vol. 1, no. 4, pp. 457–483, 1967.
  * [44] G. W. Anderson, A. Guionnet, and O. Zeitouni, _An introduction to random matrices_. Cambridge university press, 2010, no. 118.
  * [45] H. Ammari, Y. T. Chow, and J. Zou, “Phased and phaseless domain reconstructions in the inverse scattering problem via scattering coefficients,” _SIAM J. Appl. Math._ , vol. 76, no. 3, pp. 1000–1030, 2016\.
  * [46] A. Fannjiang and T. Strohmer, “The numerics of phase retrieval,” _Acta Numer._ , vol. 29, pp. 125–228, 2020.
  * [47] Z. Yuan and H. Wang, “Phase retrieval via reweighted Wirtinger flow,” _Appl. Opt._ , vol. 56, no. 9, pp. 2418–2427, 2017.