# RIS-Position and Orientation Estimation in MIMO-OFDM Systems with Practical
Scatterers

Sheng Hong, Minghui Li, Cunhua Pan, Marco Di Renzo, , Wei Zhang, , Lajos Hanzo
Sheng Hong and Minghui Li are with the School of Information Engineering of
Nanchang University, Nanchang 330031, China. (email<EMAIL_ADDRESS>minghuili@email.ncu.edu.cn). Cunhua Pan is with the National Mobile
Communications Research Laboratory, Southeast University, Nanjing, 210096,
China. (e-mail:cpan@seu.edu.cn). M. Di Renzo is with Université Paris-Saclay,
CNRS, CentraleSupélec, Laboratoire des Signaux et Systèmes, 3 Rue Joliot-
Curie, 91192 Gif-sur-Yvette, France. (marco.di-renzo@universite-paris-
saclay.fr). Wei Zhang is with the school of electrical engineering and
telecommunications, the University of New South Wales, Sydney, Australia
(e-mail: wzhang@ee.unsw.edu.au). Lajos Hanzo is with the University of
Southampton. (email: lh@ecs.soton.ac.uk). (Corresponding author: Cunhua Pan)

###### Abstract

In this paper, we investigate the problem of estimating the position and the
angle of rotation of a mobile station (MS) in a millimeter wave (mmWave)
multiple-input-multiple-output (MIMO) system aided by a reconfigurable
intelligent surface (RIS). The virtual line-of-sight (VLoS) link created by
the RIS and the non-line-of-sight (NLoS) links that originate from scatterers
in the considered environment are utilized to facilitate the estimation. A
two-step positioning scheme is exploited, where the channel parameters are
first acquired, and the position-related parameters are then estimated. The
channel parameters are obtained through a coarser and a subsequent finer
estimation processes. As for the coarse estimation, the distributed compressed
sensing orthogonal simultaneous matching pursuit (DCS-SOMP) algorithm, the
maximum likelihood (ML) algorithm, and the discrete Fourier transform (DFT)
are utilized to separately estimate the channel parameters. The obtained
channel parameters are then jointly refined by using the space-alternating
generalized expectation maximization (SAGE) algorithm, which circumvents the
high-dimensional optimization issue of ML estimation. Departing from the
estimated channel parameters, the positioning-related parameters are
estimated. The performance of estimating the channel-related and position-
related parameters is theoretically quantified by using the Cramer-Rao lower
bound (CRLB). Simulation results demonstrate the superior performance of the
proposed positioning algorithms.

###### Index Terms:

Reconfigurable Intelligent Surface, two-step positioning, channel parameter
estimation, position and orientation estimation, Cramer-Rao lower bound.

##  I Introduction

Reconfigurable intelligent surface (RIS), has recently emerged as a promising
technology to realize reconfigurable wireless propagation environments for
future communication networks [1] [2] [3]. Besides the RIS-assisted
communications [4, 5, 6, 7, 8], the RIS also triggers the exciting new
opportunities for localization in future 6G systems [9]. The RISs can be
integrated into the existing radio localization systems with low cost and high
energy efficiency, and enhance the positioning accuracy and reliability. When
the positioning signal falls into blind area, or the line-of-sight (LoS) link
is obstructed by obstacles, the RIS can restore the high-precision positioning
ability by creating a virtual LoS (VLoS) link. Through cooperating with other
anchor nodes (e.g., BS), the RIS can expand the localization range of agent
nodes (e.g., MS). Furthermore, the large physical aperture of RIS plane
enables the high-precise estimations of angle of arrival (AOA) and angle of
departure (AOD) at the RIS, which equivalently increases the angle information
for positioning.

Due to the potential benefits of RIS in the context of radio localization,
research on RIS-assisted wireless positioning has attracted increasing
attention. To evaluate the potential of RIS on positioning, the positioning
performance measured by position error bound (PEB) and orientation error bound
(OEB) is analyzed in [10, 11, 12] with different kinds of system models.

The RIS-assisted positioning algorithms are also investigated in several
works. In [13], a single-antenna MS was located on a three-dimensional (3D)
space in a multiple-RIS-aided multiple-input-single-output (MISO) system.
Therein, the maximum likelihood estimation (MLE) of the AOA and AOD was
obtained first, and then the location coordinates were estimated. In [13],
however, only the angles were exploited for positioning, and the non-line-of-
sight (NLoS) links that originate from the scatterers in the environment were
regarded as Gaussian noise. In [14], a single-antenna user was located on a
two-dimensional (2D) space with the aid of a single RIS in MISO systems.
However, the considered system model was simple, the multipaths originated
from practical scatters were also not considered, and the formulated MLE
problem for direct positioning had a very low dimensionality. Other research
works [15, 16, 17] considered the RIS-aided single-input-single-output (SISO)
system, restricted the analysis to a few channel parameters to formulate low-
dimensional MLE problems, and retrieved the coordinates through simple
geometric calculations. The work in [18] first considered using the multipath
parameters in RIS-aided multiple-input-multiple-output orthogonal frequency
division multiplexing (MIMO-OFDM) systems to realize user localization,
however, the proposed mechanism relied on the specifically designed twin-RIS
structures, which will increase the power consumption and the cost of
deployment in practical scenarios.

There also exist works considering using the multipaths to aid localization in
MIMO-OFDM systems without RIS [19, 20]. However, these methods cannot be
directly extended into the case that RIS is integrated into MIMO-OFDM systems.
Different from active arrays, it is more challenging to retrieve parameters at
the RIS. Especially for a planar RIS and a scattering scenario, the parameters
at RIS to be estimated will increase greatly, which results a high-dimensional
parameter estimation problem difficult to be solved.

Motivated by these considerations, we study the positioning algorithm for
application to RIS-aided MIMO-OFDM systems with practical scatterers, where an
uplink MIMO-OFDM system assisted by a single planar RIS is considered, and all
available channel parameters, including those on the VLoS link made available
by the RIS and the NLoS links due to the presence of scatterers, are exploited
to locate an MS in a 3D space. The traditional two-step positioning scheme is
utilized, according to which the channel parameters are first estimated, and
the location coordinates and rotation angle of the MS are then estimated. For
estimating the channel parameters, we formulate an MLE problem, which results
in two main challenges: The need of good initial values to execute the
algorithms and the high dimensionality of the problem at hand. To obtain good
initial values, we design a RIS phase shift profile to enable decoupling the
channel parameters, and estimating these parameters sequentially.
Specifically, the AODs at the MS are estimated by using the distributed
compressed sensing simultaneous orthogonal matching pursuit (DCS-SOMP)
algorithm [21] and the maximum likelihood (ML) algorithm. The AOAs at the RIS
are estimated by using the DCS-SOMP algorithm as well, and the times of
arrival (TOAs) are estimated by using the discrete Fourier transform (DFT)
approach. To tackle the high dimensionality of the MLE problem, we utilize the
space-alternating generalized expectation maximization (SAGE) algorithm [22],
which allows us to optimize the parameters associated to different channel
links alternately, while offering good estimation performance approaching that
of the MLE. Once the channel parameters are estimated, the location
coordinates and the rotation angle of the MS are calculated in a closed-form
expression by using geometric relationships. By taking these closed-form
solutions as initial values, the positioning problem is formulated as an
equivalent MLE problem, and is solved by using the Levenberg-Marquardt
algorithm (LMA) [23].

Specifically, the main contributions made by this paper are the following.

  1. 1.

Instead of treating the NLoS links from the deleterious scatterers as noise or
interference [13, 14, 15, 16, 17], we beneficially harness them for 3D
positioning in RIS-aided MIMO-OFDM systems.

  2. 2.

We propose a position and orientation estimation algorithm based on MLE in
RIS-aided MIMO-OFDM systems in the face of practical scatterers. Compared to
[13, 14, 15, 16, 17], the system model considered is more general and
practical, explicitly ULAs are assumed at the BS and the MS, the RIS is
modeled as a uniform planar array (UPA), a 3D positioning problem is
considered, and a practical scattering scenario is assumed. Thus, a very high-
dimensional parameter estimation problem is solved by the MLE, and all the
available channel parameters are beneficially harnessed for positioning.

  3. 3.

The channel parameters are first estimated coarsely, and are then refined by
using the SAGE algorithm. A bespoke RIS phase shift profile is designed to
achieve parameter decoupling in the coarse estimation process. The two-prong
strategy of coarse estimation complemented by refined estimation mitigates the
error propagation and ensures a high positioning accuracy, despite using some
parameter dimensionality reduction processes.

  4. 4.

The CRLB for the estimates of the channel parameters, and the PEB/OEB for the
estimates of the position are derived as performance benchmarks for the
proposed algorithms. Simulation results show that the estimation errors of the
proposed algorithms approach the error bounds, when the transmit power is
sufficiently high.

The remainder of this paper is organized as follows. Section II gives the
system and channel model. In Section III, the estimation problems are
formulated, and the estimation framework is illustrated. The channel
parameters are coarsely estimated in Section IV, and are refined jointly in
Section V. The position coordinates and orientation angle of the MS are
estimated in Section VI. Section VII derives the CRLB, PEB and OEB. In Section
VIII, simulation results are provided. Section IX concludes this work.

#### Notations

The symbols $[\cdot]^{*}$ , $[\cdot]^{T}$ and $[\cdot]^{H}$ denote the
conjugate, transpose and Hermitian transpose, respectively.
$[\mathbf{X}]_{u,v}$ denotes the element in the $u$-th row and $v$-th column
of matrix $\mathbf{X}$. $[\mathrm{\mathbf{x}}]_{u}$ denotes the $u$-th element
of vector $\mathbf{x}$. ${\mathbf{I}}_{M}$ represents an $M\times M$ identity
matrix. $\mathrm{diag}(\mathbf{x})$ represents a diagonal matrix by placing
the elements of vector $\mathbf{x}$ on the diagonal, while
$\mathrm{diag}(\mathrm{\mathbf{X}}$) means extracting the diagonal elements of
matrix $\mathbf{X}$ to construct a vector. The Hadamard, Kronecker and Khatri-
Rao products are represented by $\circ$, $\otimes$ and $\odot$, respectively.
The symbols $\mathrm{Tr}(\cdot)$, $\mathfrak{Re}(\cdot)$,
$\mathfrak{Im}(\cdot)$, $\mathrm{det}(\cdot)$ and $\|\cdot\|$ respectively
denote the trace operation, extracting the real part, extracting the imaginary
part, computing the determinant and norm. $\mathbb{E}[\cdot]$ denotes the
expectation. $\left\lceil\cdot\right\rceil$ represents the ceiling integer.
$\angle$ denotes the phase angle. $\left\langle\cdot\right\rangle$ denotes the
inner product.

##  II System and Channel Model

###  II-A System Model

Figure 1: The RIS-aided positioning system.

As shown in Fig. 1, we consider an uplink RIS-aided mmWave MIMO system
operating at a carrier frequency of $f_{c}$ (corresponding to wavelength of
$\lambda=c/f_{c}$) and bandwidth of $B$, in which an RIS is adopted to locate
the MS. The BS is equipped with a uniform linear array (ULA) of $N_{{\rm b}}$
antennas, and the MS is equipped with a ULA of $N_{{\rm m}}$ antennas. The
reflecting elements on the RIS constitute an UPA of $N_{{\rm
r}}=N_{{\rm{a}}}N_{{\rm{e}}}$ antennas, where $N_{{\rm{a}}}$ passive
reflecting elements are on the horizontal direction of the UPA, and
$N_{{\rm{e}}}$ passive reflecting elements are on the vertical direction of
the UPA.

The ULA of the BS is placed parallel to the $x$ axis with the center located
at ${\bf b}=[b_{x},b_{y},b_{z}]^{T}\in\mathbb{R}^{3}$. The UPA of the RIS is
placed parallel to the ${y-o-z}$ plane with their centers located at ${\bf
r}=[r_{x},r_{y},r_{z}]^{T}$. The MS moves on a plane parallel to the ${x-o-y}$
plane, and the center of the ULA of the MS is located at ${\bf
m}=[m_{x},m_{y},m_{z}]^{T}$. The ULA of the MS rotates by an angle
$\alpha\in[0,\pi)$ relative to the $x$ axis. We assume that the RIS and the BS
are relatively close and thus there is only one path between them. On the
other hand, since the MS is generally deployed on the ground with rich
scatterers, we assume that there are $Q$ scatterers between the MS and the
RIS. The location of the $q$-th scatterer is denoted by ${\bf
s}^{q}=[s_{x}^{q},s_{y}^{q},s_{z}^{q}]^{T},q=1,2,\cdots,Q$, which is unknown.
Generally, once the RIS and the BS are deployed, the coordinates of
$\mathbf{r}$ and ${\bf b}$ are known and invariant. We only need to locate the
MS, which is equivalent to estimating the coordinate ${\bf m}$ and the
rotation angle $\alpha$.

###  II-B Transmitter Model

The OFDM signals are exploited for an uplink position estimation in the RIS-
aided mmWave MIMO system. At the MS, $T$ OFDM symbols with $N$ subcarriers are
transmitted sequentially. The transmitted signal over subcarrier $n$ at the
$t$-th time slot can be expressed as
$\mathbf{x}_{t}[n]=\sqrt{P_{\mathrm{m}}}\mathbf{\bar{x}}_{t}[n]\in\mathbb{C}^{N_{m}\times
1}$, where $n=1,2,\cdots,N$ and $t=1,2,\cdots,T$. $P_{\mathrm{m}}$ is the
transmit power of the MS, and $\mathbf{\bar{x}}_{t}[n]$ is the power-
normalized pilot signal. We assume bandwidth $B\ll f_{c}$ for a narrow-band
model.

###  II-C Wireless Channel Model

The mmWave propagation channels can be described by the array steering
vectors. Generally, the steering vector of an ULA can be represented as

$\mathbf{f}(u_{{\rm sf}},N_{{\rm ant}})=[1,e^{-j2\pi u_{{\rm
sf}}},\cdots,e^{-j2\pi(N_{{\rm ant}}-1)u_{{\rm
sf}}}]^{T}\in\mathbb{C^{\mathit{N_{{\rm ant}}\times\mathrm{1}}}},$ (1)

where $u_{{\rm sf}}=\frac{d}{\lambda}\sin\theta$ is the spatial frequency at
the ULA corresponding to the physical angle $\theta$, $d$ is the distance
between adjacent antennas, and $N_{{\rm ant}}$ is the number of antennas. We
assume that $d\leq\frac{\lambda}{2}$ so that there is a one-to-one
relationship between the spatial frequencies and the physical angles on either
the ULA or UPA.

As shown in Fig. 1, there is only an LoS path from the RIS to the BS. There
are ($Q+1$) paths from the MS to the RIS, in which the ($q=0$)-th path
represents the LoS path, and the others are the NLoS paths. We term the
cascaded path of the MS-RIS link and the RIS-BS link as the VLoS path. Then,
the channel of the RIS-BS link ${\bf H}_{{\rm RB}}[n]\in\mathbb{C}^{N_{{\rm
b}}\times N_{\mathrm{r}}}$ and the channel of the MS-RIS link ${\bf H}_{{\rm
MR}}[n]\in\mathbb{C}^{N_{\mathrm{r}}\times N_{\mathrm{m}}}$ can be
respectively modeled as

$\displaystyle{\bf H}_{{\rm RB}}[n]$ $\displaystyle=\delta_{{\rm
RB},0}e^{-j2\pi\tau_{{\rm RB},0}\frac{(n-1)B}{N}}{\bf a}_{{\rm
B}}(\vartheta_{{\rm r},0}){\bf a}_{{\rm R}}^{H}(\omega_{{\rm out},0}^{{\rm
a}},\omega_{{\rm out},0}^{{\rm e}}),$ (2) $\displaystyle{\bf H}_{{\rm MR}}[n]$
$\displaystyle=\sum\limits_{q=0}^{Q}\delta_{{\rm MR},q}e^{-j2\pi\tau_{{\rm
MR},q}\frac{(n-1)B}{N}}{\bf a}_{{\rm R}}(\omega_{{\rm in},q}^{{\rm
a}},\omega_{{\rm in},q}^{{\rm e}}){\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm
t},q}),$ (3)

where $\delta_{{\rm RB},0}$ and $\tau_{{\rm RB},0}$ are respectively the
channel gain and TOA of the direct RIS-BS path, while $\delta_{{\rm MR},q}$
and $\tau_{{\rm MR},q}$ are respectively the channel gain and TOA of the
$q$-th path from the MS to the RIS. By considering a narrow-band model, the
signal wavelength at the $n$-th subcarrier satisfies
$\lambda_{n}\approx\lambda$, and the steering vectors ${\bf a}_{{\rm
B}}(\vartheta_{{\rm r},0})\in\mathbb{C^{\mathit{N_{{\rm
b}}\times\mathrm{1}}}}$ at the BS, ${\bf a}_{{\rm M}}(\vartheta_{{\rm
t},q})\in\mathbb{C^{\mathit{N_{{\rm m}}\times\mathrm{1}}}}$ at the MS, and
${\bf a}_{{\rm
R}}(\omega_{x,y}^{{\rm{a}}},\omega_{x,y}^{{\rm{e}}})\in\mathbb{C^{\mathit{N_{{\rm
r}}\times\mathrm{1}}}}$ at the RIS are the same for all $N$ subcarriers. These
array steering vectors can be represented by ${\bf a}_{{\rm
B}}(\vartheta_{{\rm r},0})=\mathbf{f}(\vartheta_{{\rm r},0},N_{\mathrm{b}})$,
${\bf a}_{{\rm M}}(\vartheta_{{\rm t},q})=\mathbf{f}(\vartheta_{{\rm
t},q},N_{\mathrm{m}})$ and ${\bf a}_{{\rm
R}}(\omega_{x,y}^{{\rm{a}}},\omega_{x,y}^{{\rm{e}}})=\mathbf{f}(\omega_{x,y}^{{\rm{e}}},N_{{\rm{e}}})\otimes\mathbf{f}(\omega_{x,y}^{{\rm{a}}},N_{{\rm{a}}})$,
where $\vartheta_{{\rm
r},0}=\frac{d_{\textrm{{\rm{BS}}}}}{\lambda}\sin\theta_{{\rm r},0}$,
$\vartheta_{{\rm t},q}=\frac{d_{\textrm{{\rm{MS}}}}}{\lambda}\sin\theta_{{\rm
t},q}$,
$\omega_{x,y}^{{\rm{a}}}=\frac{d_{\textrm{{\rm{RIS,a}}}}}{\lambda}\sin\psi_{x,y}\sin\phi_{x,y}$,
$\omega_{x,y}^{{\rm{e}}}=\frac{d_{\textrm{{\rm{RIS,e}}}}}{\lambda}\cos\phi_{x,y}$,
$x\in\\{{\rm out,in}\\}$, $y\in\\{0,q\\}$ are the spatial frequencies.

The physical angles are defined as follows. The positive elevation AOD
$\phi_{{\rm out,}0}$ at the RIS of the RIS-BS direct link is the
anticlockwise-rotation angle of the wave vector with respect to the positive
$z$ axis. The positive azimuth AOD $\psi_{{\rm out,}0}$ at the RIS of the RIS-
BS direct link is the anticlockwise-rotation angle of the vector on the
$x-o-y$ plane projected from the wave vector with respect to the positive $x$
axis. The positive AOA $\theta_{{\rm r,}0}$ at the BS of the RIS-BS direct
link is the anticlockwise-rotation angle of the wave vector with respect to
the normal direction of the ULA at the BS, which is on the plane formed by the
line collinear with the ULA and the line collinear with the wave vector, and
is denoted by the green dotted line in Fig. 1. Similarly, the positive AOD
$\theta_{\mathrm{t},q}$ at the MS of the MS-Scatterer-$q$-RIS link is the
anticlockwise-rotation angle of the wave vector with respect to the normal
direction of the ULA at the MS. The positive elevation AOA
$\phi_{\mathrm{in},q}$ at the RIS of the MS-Scatterer-$q$-RIS link is the
anticlockwise-rotation angle of the wave vector with respect to the positive
$z$ axis. The positive azimuth AOA $\psi_{\mathrm{in},q}$ at the RIS of the
MS-Scatterer-$q$-RIS link is the anticlockwise-rotation angle of the vector on
the $x-o-y$ plane projected from the wave vector with respect to the positive
$x$ axis. The angles $\phi_{{\rm out,}0}$, $\psi_{{\rm out,}0}$, $\theta_{{\rm
r,}0}$, $\theta_{\mathrm{t},q}$, $\phi_{\mathrm{in},q}$ and
$\psi_{\mathrm{in},q}$ become negative when the angles defined above are
rotated in clockwise.

By using (2) and (3), the complete channel from the MS to the BS via the RIS
is expressed as

$\displaystyle{\bf H}_{t}[n]$ $\displaystyle={\bf H}_{{\rm
RB}}[n]\text{$\bm{{\Omega}}$}_{t}{\bf H}_{{\rm MR}}[n]$
$\displaystyle=\sum\limits_{q=0}^{Q}\delta_{{\rm RB},0}\delta_{{\rm
MR},q}e^{-j2\pi(\tau_{{\rm RB},0}+\tau_{{\rm MR},q})\frac{(n-1)B}{N}}{\bf
a}_{{\rm B}}(\vartheta_{{\rm r},0}){\bf a}_{{\rm R}}^{H}(\omega_{{\rm
out},0}^{{\rm a}},\omega_{{\rm out},0}^{{\rm e}})\text{$\bm{\Omega}$}_{t}{\bf
a}_{{\rm R}}(\omega_{{\rm in},q}^{{\rm a}},\omega_{{\rm in},q}^{{\rm e}}){\bf
a}_{{\rm M}}^{H}(\vartheta_{{\rm t},q})$
$\displaystyle\triangleq\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\bf
a}_{{\rm B}}(\vartheta_{{\rm r},0}){\bf a}_{{\rm R}}^{H}(\omega_{{\rm
out},0}^{{\rm a}},\omega_{{\rm out},0}^{{\rm e}})\text{$\bm{\Omega}$}_{t}{\bf
a}_{{\rm R}}(\omega_{{\rm in},q}^{{\rm a}},\omega_{{\rm in},q}^{{\rm e}}){\bf
a}_{{\rm M}}^{H}(\vartheta_{{\rm t},q}),$ (4)

where $\delta_{q}\triangleq\delta_{{\rm RB},0}\delta_{{\rm MR},q}$ and
$\tau_{q}\triangleq\tau_{{\rm RB},0}+\tau_{{\rm MR},q}$. For later use, we
introduce
$\tilde{\delta}_{q}[n]\triangleq\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}$.
$\text{$\bm{{\Omega}}$}_{t}=\textrm{diag}\\{{\bf g}_{t}\\}$ is a diagonal
matrix with ${\bf g}_{t}$ on the diagonal. ${\bf g}_{t}$ contains the
reflection coefficients of the RIS reflecting elements at the $t$-th time
slot, and can be represented by ${\bf
g}_{t}=[e^{j\varphi_{t}^{(1)}},e^{j\varphi_{t}^{(2)}},\cdots,e^{j\varphi_{t}^{(N_{\textrm{r}})}}]^{T}$,
where $\varphi_{t}^{(n)},n=1,2,\cdots,N_{\textrm{r}}$ are the phase shifts of
reflection coefficients. In (4), there are $(Q+1)$ paths, and the $q$-th path
can be characterized by the channel parameters
$[\tau_{q},\delta_{q},\theta_{{\rm r},0},\phi_{{\rm out},0},\psi_{{\rm
out},0},\phi_{{\rm in},q},\psi_{{\rm in},q},\theta_{{\rm
t},q}]^{T},q=0,1,\cdots,Q$.

###  II-D Sparse Representation of mmWave Channel

Based on (2) and (3), the two channel matrices can be re-expressed in an
angular domain representation as

$\displaystyle{\bf H}_{{\rm RB}}[n]$ $\displaystyle={\bf A}_{{\rm
B}}\bm{{\varLambda}}_{\textrm{RB}}[n]\mathbf{A}_{{\rm R}}^{H},$ (5a)
$\displaystyle{\bf H}_{{\rm MR}}[n]$ $\displaystyle={\bf A}_{{\rm
R}}\bm{{\varLambda}}_{\textrm{MR}}[n]{\bf A}_{{\rm M}}^{H},$ (5b)

where the three matrices of ${\bf A}_{{\rm B}}\in\mathbb{C}^{N_{{\rm b}}\times
G_{\mathrm{b}}}$ and ${\bf A}_{{\rm R}}\in\mathbb{C}^{N_{{\rm r}}\times
G_{\mathrm{r}}}$, ${\bf A}_{{\rm M}}\in\mathbb{C}^{N_{{\rm m}}\times
G_{\mathrm{m}}}$ are over-complete dictionary matrices in the angle domain.
The columns of ${\bf A}_{{\rm B}}$, ${\bf A}_{{\rm R}}$ and ${\bf A}_{{\rm
M}}$ consist of the array steering vectors ${\bf a}_{{\rm B}}(\vartheta)$,
${\bf a}_{{\rm R}}(\omega^{{\rm a}},\omega^{{\rm e}})$ and ${\bf a}_{{\rm
M}}(\vartheta)$ respectively at the specific angles. The $G_{\mathrm{b}}$,
$G_{\mathrm{r}}$ and $G_{\mathrm{m}}$ are the numbers of the per-discretized
grids. Since there is only one path between the RIS and the BS, the
$\bm{{\varLambda}}_{\textrm{RB}}[n]\in\mathbb{C}^{G_{\mathrm{b}}\times
G_{\mathrm{r}}}$ is a sparse matrix containing only one non-zero element given
by $\tilde{\delta}_{{\rm RB,}0}[n]\triangleq\delta_{{\rm
RB,}0}e^{-j2\pi\tau_{{\rm RB,}0}\frac{(n-1)B}{N}}$. The
$\bm{{\varLambda}}_{\textrm{MR}}[n]\in\mathbb{C}^{G_{\mathrm{r}}\times
G_{\mathrm{m}}}$ is a sparse matrix containing ($Q+1$) non-zero elements given
by $\\{\tilde{\delta}_{{\rm MR,}q}[n]\triangleq\delta_{{\rm
MR,}q}e^{-j2\pi\tau_{{\rm MR,}q}\frac{(n-1)B}{N}}\\}_{q=0,1,\cdots,Q}$. By
using (1), ${\bf A}_{{\rm B}}$ and ${\bf A}_{{\rm M}}$ can be respectively
written as

$\displaystyle{\bf A}_{{\rm B}}$
$\displaystyle=[\mathbf{f}(-1\frac{d_{\textrm{{\rm{BS}}}}}{\lambda},N_{\mathrm{b}}),\mathbf{f}((-1+\frac{2}{G_{\mathrm{b}}})\frac{d_{\textrm{{\rm{BS}}}}}{\lambda},N_{\mathrm{b}}),\cdots,\mathbf{f}((1-\frac{2}{G_{\mathrm{b}}})\frac{d_{\textrm{{\rm{BS}}}}}{\lambda},N_{\mathrm{b}})],$
(6a) $\displaystyle{\bf A}_{{\rm M}}$
$\displaystyle=[\mathbf{f}(-1\frac{d_{\textrm{{\rm{MS}}}}}{\lambda},N_{\mathrm{m}}),\mathbf{f}((-1+\frac{2}{G_{\mathrm{m}}})\frac{d_{\textrm{{\rm{MS}}}}}{\lambda},N_{\mathrm{m}}),\cdots,\mathbf{f}((1-\frac{2}{G_{\mathrm{m}}})\frac{d_{\textrm{{\rm{MS}}}}}{\lambda},N_{\mathrm{m}})].$
(6b)

Since the UPA is utilized at the RIS, ${\bf A}_{{\rm R}}$ can be respectively
written as

$\displaystyle{\bf A}_{{\rm R}}$ $\displaystyle={\bf A}_{{\rm R}}^{{\rm
e}}\otimes{\bf A}_{{\rm R}}^{{\rm a}},$ (7)

where ${\bf A}_{{\rm R}}^{{\rm e}}\in\mathbb{C}^{N_{{\rm{e}}}\times
G_{\mathrm{e}}}$ and ${\bf A}_{{\rm R}}^{{\rm
a}}\in\mathbb{C}^{N_{{\rm{a}}}\times G_{\mathrm{a}}}$ can be written as

$\displaystyle{\bf A}_{{\rm R}}^{{\rm a}}$
$\displaystyle=[\mathbf{f}(-1\frac{d_{\textrm{{\rm{RIS,a}}}}}{\lambda},N_{{\rm{a}}}),\mathbf{f}((-1+\frac{2}{G_{\mathrm{a}}})\frac{d_{\textrm{{\rm{RIS,a}}}}}{\lambda},N_{{\rm{a}}}),\cdots,\mathbf{f}((1-\frac{2}{G_{\mathrm{a}}})\frac{d_{\textrm{{\rm{RIS,a}}}}}{\lambda},N_{{\rm{a}}})],$
(8a) $\displaystyle{\bf A}_{{\rm R}}^{{\rm e}}$
$\displaystyle=[\mathbf{f}(-1\frac{d_{\textrm{{\rm{RIS,e}}}}}{\lambda},N_{{\rm{e}}}),\mathbf{f}((-1+\frac{2}{G_{\mathrm{e}}})\frac{d_{\textrm{{\rm{RIS,e}}}}}{\lambda},N_{{\rm{e}}}),\cdots,\mathbf{f}((1-\frac{2}{G_{\mathrm{e}}})\frac{d_{\textrm{{\rm{RIS,e}}}}}{\lambda},N_{{\rm{e}}})],$
(8b)

and $G_{\mathrm{r}}=G_{\mathrm{a}}\times G_{\mathrm{e}}$.

###  II-E Received Signal Model

We assume that the direct path from the MS to the BS is obstructed by
obstacles. Then, the received signal at the $t$-th time slot of the $n$-th
subcarrier can be written as

$\displaystyle\mathbf{y}_{t}[n]$ $\displaystyle={\bf
H}_{t}[n]\mathbf{x}_{t}[n]+\mathbf{z}_{t}[n]$ $\displaystyle={\bf H}_{{\rm
RB}}[n]\text{$\bm{{\Omega}}$}_{t}{\bf H}_{{\rm
MR}}[n]\mathbf{x}_{t}[n]+\mathbf{z}_{t}[n],$ (9)

where $t=1,2,\cdots,T$. $\mathbf{x}_{t}[n]$ is the pilot signal transmitted at
the $t$-th time slot, which is known at the BS. The transmit power of the
pilot signal is given by
$P_{\mathrm{m}}=\left\|\mathbf{x}_{t}[n]\right\|^{2}$.
$\mathbf{z}_{t}[n]\sim\mathcal{CN}(0,\sigma^{2}\mathbf{I}_{N_{{\rm b}}})$ is
the noise vector, where $\sigma^{2}$ denotes the noise power.

##  III Estimation Problem Formulation

###  III-A Parameters to be Estimated

From the received signal $\mathbf{y}_{t}[n]$ in (9), we should first estimate
all the unknown channel parameters, which are collected in the vector
$\bm{{\eta}}\in\mathbb{C}^{6(Q+1)}$ as

$\displaystyle\bm{{\eta}}$
$\displaystyle=[\bm{{\eta}}_{0}^{T},\bm{{\eta}}_{1}^{T},\cdots,\bm{{\eta}}_{q}^{T},\cdots,\bm{{\eta}}_{Q}^{T}]^{T},$
(10)

where $\bm{{\eta}}_{q}=[\tau_{q},\bm{{h}}_{q}^{T},\bm{{\theta}}_{q}^{T}]^{T}$,
$\bm{{h}}_{q}=[\delta_{q,\textrm{R}},\delta_{q,\textrm{I}}]^{T}$,
$\delta_{q,\textrm{R}}=\mathfrak{Re}(\delta_{q})$,
$\delta_{q,\textrm{I}}=\mathfrak{Im}(\delta_{q})$ and
$\bm{{\theta}}_{q}=[\theta_{{\rm t},q},\phi_{{\rm in},q},\psi_{{\rm
in},q}]^{T}$. We term these channel parameters as the intermediate parameters.

Then, from these intermediate parameters, we estimate the location coordinates
of the MS, the rotation angle of the MS, the location coordinates of the
scatterer, etc., which are termed as the final parameters. We collect the
final parameters in a vector $\tilde{\bm{{\eta}}}\in\mathbb{C}^{2(Q+1)+3Q+4}$
, which can be represented as

$\displaystyle\tilde{\bm{{\eta}}}$
$\displaystyle=[\bm{{h}}_{0}^{T},\cdots,\bm{{h}}_{q}^{T},\cdots,\bm{{h}}_{Q}^{T},\tilde{\bf{m}}^{T},({\bf
s}^{1})^{T},\cdots,({\bf s}^{q})^{T},\cdots,({\bf s}^{Q})^{T}]^{T},$ (11)

where $\tilde{\bf{m}}=[{\bf m}^{T},\alpha]^{T}$. The final parameters in
$\tilde{\bm{{\eta}}}$ are related with the intermediate parameters
$\bm{{\bf\eta}}$ according to the geometric relationship as follows.

###  III-B The Geometric Relationship

Obviously, the channel parameters $\\{\theta_{{\rm r},0},\phi_{{\rm
out},0},\psi_{{\rm out},0}\\}$ in the RIS-BS link can be calculated directly
from the coordinates of $\mathbf{b}$ and $\mathbf{r}$. Other channel
parameters $\\{\tau_{q},\delta_{q},\phi_{{\rm in},q},\psi_{{\rm
in},q},\theta_{{\rm t},q}\\}_{q=0,1,\cdots,Q}$ are determined by the location
coordinates ${\bf m}$ of the MS, the rotation angle $\alpha$ of the MS, and
the location coordinates ${\bf s}^{q}$ of the scatterer according to the
geometric relationship described in Fig. 1, which is illustrated by equations
as follows.

####  III-B1 The Angles in the MS-RIS Link

The angles of the LoS path from the MS to the RIS are $\theta_{{\rm t},0}$,
$\psi_{{\rm in},0}$ and $\phi_{{\rm in},0}$. These angles are related to
$\tilde{\bf{m}}$ as

$\displaystyle\theta_{{\rm t},0}$
$\displaystyle=\arcsin\frac{(r_{x}-m_{x})\cos\alpha-(r_{y}-m_{y})\sin\alpha}{\sqrt{(r_{x}-m_{x})^{2}+(r_{y}-m_{y})^{2}+(r_{z}-m_{z})^{2}}},$
(12a) $\displaystyle\psi_{{\rm in},0}$
$\displaystyle=\pi-\arcsin\frac{r_{y}-m_{y}}{\sqrt{(r_{x}-m_{x})^{2}+(r_{y}-m_{y})^{2}}},$
(12b) $\displaystyle\phi_{{\rm in},0}$
$\displaystyle=\arccos\frac{r_{z}-m_{z}}{\sqrt{(r_{x}-m_{x})^{2}+(r_{y}-m_{y})^{2}+(r_{z}-m_{z})^{2}}}.$
(12c)

The angles of the NLoS path from the MS to the RIS via scatterers are
$\theta_{{\rm t},q}$, $\psi_{{\rm in},q}$ and $\phi_{{\rm in},q}$. These
angles are related to $\tilde{\bf{m}}$ and ${\bf s}^{q}$ as

$\displaystyle\theta_{{\rm t},q}$
$\displaystyle=\arcsin\frac{(s_{x}^{q}-m_{x})\cos\alpha-(s_{y}^{q}-m_{y})\sin\alpha}{\sqrt{(s_{x}^{q}-m_{x})^{2}+(s_{y}^{q}-m_{y})^{2}+(s_{z}^{q}-m_{z})^{2}}},$
(13a) $\displaystyle\psi_{{\rm in},q}$
$\displaystyle=\pi-\arcsin\frac{r_{y}-s_{y}^{q}}{\sqrt{(r_{x}-s_{x}^{q})^{2}+(r_{y}-s_{y}^{q})^{2}}},$
(13b) $\displaystyle\phi_{{\rm in},q}$
$\displaystyle=\arccos\frac{r_{z}-s_{z}^{q}}{\sqrt{(r_{x}-s_{x}^{q})^{2}+(r_{y}-s_{y}^{q})^{2}+(r_{z}-s_{z}^{q})^{2}}},$
(13c)

where $q=1,2,\cdots,Q$.

####  III-B2 The Angles in the RIS-BS Link

The angles of the LoS path from the RIS to the BS are $\theta_{{\rm r},0}$,
$\psi_{{\rm out},0}$ and $\phi_{{\rm out},0}$. These angles are related to
$\mathbf{b}$ and $\mathbf{r}$ as

$\displaystyle\theta_{{\rm r},0}$
$\displaystyle=\arcsin\frac{b_{x}-r_{x}}{\sqrt{(b_{x}-r_{x})^{2}+(b_{y}-r_{y})^{2}+(b_{z}-r_{z})^{2}}},$
(14a) $\displaystyle\psi_{{\rm out},0}$
$\displaystyle=\arcsin\frac{b_{y}-r_{y}}{\sqrt{(b_{x}-r_{x})^{2}+(b_{y}-r_{y})^{2}}},$
(14b) $\displaystyle\phi_{{\rm out},0}$
$\displaystyle=\arccos\frac{b_{z}-r_{z}}{\sqrt{(b_{x}-r_{x})^{2}+(b_{y}-r_{y})^{2}+(b_{z}-r_{z})^{2}}}.$
(14c)

It is noted that $\theta_{{\rm r},0}$, $\psi_{{\rm out},0}$ and $\phi_{{\rm
out},0}$ can be calculated directly from $\mathbf{b}$ and $\mathbf{r}$, and
can be taken as known values.

####  III-B3 The TOAs in All Paths

The TOA is equivalent to the distance. Thus, the TOAs are related to the
coordinates of all nodes as

$\displaystyle\tau_{0}$ $\displaystyle=\left\|{\bf r}-{\bf
b}\right\|_{2}/c+\left\|{\bf m}-{\bf r}\right\|_{2}/c,$ (15a)
$\displaystyle\tau_{q}$ $\displaystyle=\left\|{\bf r}-{\bf
b}\right\|_{2}/c+\left\|{\bf s}^{q}-{\bf r}\right\|_{2}/c+\left\|{\bf m}-{\bf
s}^{q}\right\|_{2}/c,{\ }q=1,\cdots,Q.$ (15b)

###### Remark 1

According to the definition of physical angles above, the positive azimuth
angle at the RIS is defined as the anticlockwise rotation angle of the
projection vector with respect to the positive $x$ axis, while the clockwise
rotation angle of the projection vector with respect to the positive $x$ axis
is the negative azimuth angle. Thus, from the system model described in Fig.
1, the ranges of azimuth angles are $\psi_{{\rm in},0}\in[\pi,3\pi/2]$,
$\psi_{{\rm in},q}\in[\pi/2,\pi]$ and $\psi_{{\rm out},0}\in[-\pi/2,0]$,
respectively. Since the function of $\arcsin(\cdotp)$ can only generate angles
within $[-\pi/2,\pi/2]$, the expressions in (12b) and (13b) are utilized to
ensure the correct angle range based on the periodicity of $\sin(\cdotp)$
function.

###  III-C Estimation Framework

The classic two-step positioning strategy is utilized in this paper. Our task
is to first estimate the intermediate parameters in $\bm{{\bf\eta}}$ from the
received signal, then estimate the final parameters in $\tilde{\bm{{\eta}}}$
from the estimated ${\bm{{\bf\eta}}}$. Thus, two parameter estimation problems
should be formulated and solved.

Figure 2: The configuration for the RIS phase shifts. Figure 3: The flow chart
of the proposed positioning algorithm.

The first parameter estimation problem is to estimate the channel parameters
in $\bm{{\bf\eta}}$ from the received signal, which is solved by a coarse
estimation process and a refining process sequentially. In the coarse
estimation process, to circumvent the high-dimensional parameter estimation
issue, we propose to estimate the channel parameters separately by compression
sensing algorithms, thus the phase shifts of the RIS are specifically
configured as shown in Fig. 2. Assume that the total number of time slots
utilized for localization is $T$. In the first $T_{1}$ time slots, the phase
shifts of the RIS are time-invariant. The remaining $T_{2}=T-T_{1}$ time slots
are divided into $\Upsilon$ time blocks, each of which has $V$ time slots. The
phase shifts of the RIS are time-invariant in each time block, and are time-
variant over different time blocks. The coarsely-estimated channel parameters
are provided as the initial values of the refining process, where the SAGE
algorithm is utilized. The SAGE algorithm belongs to the ML algorithms, and
approaches the performance of the ML algorithms with greatly reduced
complexity.

The second parameter estimation problem is to estimate the position-related
parameters in $\tilde{\bm{{\eta}}}$ from the estimated channel parameters
$\hat{\bm{{\bf\eta}}}$. Due to the unavoidable estimation error in
$\hat{\bm{{\bf\eta}}}$, the geometric relationship described in Section III-B
can only be utilized to calculate the initial values for
$\tilde{\bm{{\eta}}}$, and then an equivalent ML algorithm is exploited to
refine the estimation for $\tilde{\bm{{\eta}}}$.

The flow chart of the proposed positioning algorithm is given in Fig. 3.

##  IV Coarse Estimation of Channel Parameters

The main purpose of coarse estimation of channel parameters in this section is
to provide good initial values for the SAGE algorithm in the refining process.
Based on the designed configuration of RIS phase shifts, the channel
parameters of AODs, AOAs and TOAs can be estimated separately.

The channel parameters of ${\bf\bm{\eta}}$ are coarsely estimated in three
sub-steps sequentially. In the first sub-step, we estimate the AOD
$\theta_{{\rm t},q},q=0,1,2,\cdots,Q$ at the MS in the first $T_{1}$ time
slots, where the phase shifts of the RIS are time-invariant. Then, by using
the estimated AODs $\theta_{{\rm t},q},q=0,1,2,\cdots,Q$, we estimate
$\omega_{q}\triangleq\cos\phi_{{\rm in},q}-\cos\phi_{{\rm out},0}$,
$\gamma_{q}\triangleq\sin\psi_{{\rm in},q}\sin\phi_{{\rm in},q}-\sin\psi_{{\rm
out},0}\sin\phi_{{\rm out},0}$ and $\tilde{\delta}_{q}[n]$ in the overall $T$
time slots in the second sub-step, where $q=0,1,2,\cdots,Q$. Since $\phi_{{\rm
out},0}$ and $\psi_{{\rm out},0}$ can be calculated directly according to
(14b) and (14c), $\phi_{{\rm in},q}$ and $\psi_{{\rm in},q}$ can be estimated
from $\omega_{q}$ and $\gamma_{q}$ in the second sub-step. In the third sub-
step, by using the estimated values of $\tilde{\delta}_{q}[n],n=1,2,\cdots,N$
in the second sub-step, the TOAs $\tau_{q}$ and channel gains $\delta_{q}$,
$q=0,1,2,\cdots,Q$ are estimated by the DFT approach and a time rotation
operation.

###  IV-A Estimation of the AODs at the MS

####  IV-A1 Estimate AODs by the DCS-SOMP Algorithm

By substituting (5b) into (9), the pilot signal received during the first
$T_{1}$ time slots can be written as

$\displaystyle\mathbf{y}_{t}[n]$ $\displaystyle={\bf H}_{{\rm
RB}}[n]\text{$\bm{{\Omega}}$}_{0}{\bf A}_{{\rm
R}}\bm{{\varLambda}}_{\textrm{MR}}[n]{\bf A}_{{\rm
M}}^{H}\mathbf{x}_{t}[n]+\mathbf{z}_{t}[n],\forall n=1,\cdots,N,$ (16)

where $t=1,2,\cdots,T_{1}$. The conjugate transpose of $\mathbf{y}_{t}[n]$ can
be expressed as

$\displaystyle\mathbf{y}_{t}^{H}[n]=\mathbf{x}_{t}^{H}[n]{\bf A}_{{\rm
M}}\text{$\bm{{\Gamma}}$}_{\textrm{MR}}[n]+\mathbf{z}_{t}^{H}[n],$ (17)

where
$\text{$\bm{{\Gamma}}$}_{\textrm{MR}}[n]\triangleq\bm{{\varLambda}}_{\textrm{MR}}^{H}[n]{\bf
A}_{{\rm R}}^{H}\text{$\bm{{\Omega}}$}_{0}^{H}{\bf H}_{{\rm RB}}^{H}[n]$.
Since $\bm{{\varLambda}}_{\textrm{MR}}[n]$ is a sparse matrix with $(Q+1)$
nonzero elements,
$\text{$\bm{{\Gamma}}$}_{\textrm{MR}}[n]\in\mathbb{C}^{G_{\textrm{r}}\times
N_{\textrm{b}}}$ is a row sparse matrix with $(Q+1)$ nonzero rows. By
collecting the received $\mathbf{y}_{t}^{H}[n]$ in the $T_{1}$ time slots, we
have

$\displaystyle\mathbf{Y}_{1}^{H}[n]$ $\displaystyle=\mathbf{X}_{1}^{H}[n]{\bf
A}_{{\rm M}}\text{$\bm{{\Gamma}}$}_{\textrm{MR}}[n]+\mathbf{Z}_{1}^{H}[n]$
$\displaystyle=\mathbf{\bm{{\Theta}}}_{{\rm{M}}}[n]\text{$\bm{{\Gamma}}$}_{\textrm{MR}}[n]+\mathbf{Z}_{1}^{H}[n],$
(18)

where
$\mathbf{Y}_{1}[n]=[\mathbf{y}_{1}[n],\mathbf{y}_{2}[n],\cdots,\mathbf{y}_{T_{1}}[n]]\in\mathbb{C}^{N_{\textrm{b}}\times
T_{1}}$,
$\mathbf{X}_{1}[n]=[\mathbf{x}_{1}[n],\mathbf{x}_{2}[n],\cdots,\mathbf{x}_{T_{1}}[n]]\in\mathbb{C}^{N_{\textrm{m}}\times
T_{1}}$,
$\mathbf{Z}_{1}[n]=[\mathbf{z}_{1}[n],\mathbf{z}_{2}[n],\cdots,\mathbf{z}_{T_{1}}[n]]\in\mathbb{C}^{N_{\textrm{b}}\times
T_{1}}$, and
$\mathbf{\bm{{\Theta}}}_{{\rm{M}}}[n]\triangleq\mathbf{X}_{1}^{H}[n]{\bf
A}_{{\rm M}}$. Note that the $k$-th row of
$\text{$\bm{{\Gamma}}$}_{\textrm{MR}}[n]$ is nonzero only when the nonzero
element of $\bm{{\varLambda}}_{\textrm{MR}}[n]$ exists on the $k$-th row of
$\bm{{\varLambda}}_{\textrm{MR}}^{H}[n]$. The $k$-th row of
$\bm{{\varLambda}}_{\textrm{MR}}^{H}[n]$ is associated with the $k$-th column
of the dictionary matrix ${\bf A}_{{\rm M}}$, which corresponds to one of the
transmit array steering vectors $\\{{\bf a}_{{\rm M}}^{H}(\theta_{{\rm
t},q})\\}_{q=0,1,\cdots,Q}$. Thus, the task of the first sub-step is to
determine the $(Q+1)$ columns of ${\bf A}_{{\rm M}}$ corresponding to non-zero
rows of $\text{$\bm{{\Gamma}}$}_{\textrm{MR}}[n]$ based on the received signal
matrix $\mathbf{Y}_{1}^{H}[n]$ in (18), which is a sparse recovery problem
that can be solved by the DCS-SOMP algorithm [21] given in Algorithm 1. For
simplicity, we assume that the pilot signal $\mathbf{X}_{1}[n]$ is the same
over all subcarriers, which can be denoted by $\mathbf{X}_{1}$. To achieve the
best performance of the DCS-SOMP algorithm, the columns in the equivalent
dictionary matrix $\mathbf{\bm{{\Theta}}}_{{\rm{M}}}$ should be mutually-
orthogonal, thus the pilot signal $\mathbf{X}_{1}$ should be properly
designed. The analytic designing method for $\mathbf{X}_{1}$ can be found in
[24], and a simpler designing method is to randomly generate the elements in
$\mathbf{X}_{1}$ from $\\{-1,1\\}$ with equal probability [25].

Algorithm 1 DCS-SOMP algorithm for AOD estimation at the MS

1: Parameter Setting: Initialize ${\bar{\bm{{\Theta}}}}_{{\rm{M}}}$ and
$\bar{{\bf A}}_{{\rm M}}$ as empty matrices. $\mathcal{K}_{0}$ is chosen to be
an empty set. Set $\mathbf{R}_{\textrm{res}}[n]=\mathbf{Y}_{1}^{H}[n]$,
$\forall n$ and ${\bm{{\Theta}}}_{{\rm{M}}}=\mathbf{X}_{1}^{H}{\bf A}_{{\rm
M}}$.

2: for $q\leq Q$ do

3:
$\bm{{\Psi}}[n]=\mathbf{\bm{{\Theta}}}_{{\rm{M}}}^{H}\mathbf{R}_{\textrm{res}}[n]$,
$\forall n$;

4:
$k=\arg\underset{l}{\max}{\sum}_{n=1}^{N}\\{\bm{{\Psi}}[n]\bm{{\Psi}}^{H}[n]\\}_{l,l}$;

5: Update the set of AOD indices $\mathcal{K}_{0}=\mathcal{K}_{0}\cup\\{k\\}$;

6:
${\bar{\bm{{\Theta}}}}_{{\rm{M}}}=[{\bar{\bm{{\Theta}}}}_{{\rm{M}}},{\bm{{\Theta}}}_{{\rm{M}}}(:,k)]$;

7: $\bar{{\bf A}}_{{\rm M}}=[\bar{{\bf A}}_{{\rm M}},{\bf A}_{{\rm M}}(:,k)]$;

8:
$\text{$\bm{{\Gamma}}$}_{\textrm{MR}}[n]=({\bar{\bm{{\Theta}}}}_{{\rm{M}}}^{H}{\bar{\bm{{\Theta}}}}_{{\rm{M}}})^{-1}{\bar{\bm{{\Theta}}}}_{{\rm{M}}}^{H}\mathbf{R}_{\textrm{res}}[n]$,
$\forall n$;

9:
$\mathbf{R}_{\textrm{res}}[n]=\mathbf{Y}_{1}^{H}[n]-{\bar{\bm{{\Theta}}}}_{{\rm{M}}}\text{$\bm{{\Gamma}}$}_{\textrm{MR}}[n]$,
$\forall n$;

10: end for

11: Output: $\mathcal{K}_{0},\bar{{\bf A}}_{{\rm
M}}(:,q)=\mathbf{f}((-1+\frac{2(k-1)}{G_{\mathrm{m}}})\frac{d_{\textrm{{\rm{MS}}}}}{\lambda},N_{\mathrm{m}})$
and $\hat{\theta}_{{\rm t},q}=\arcsin[-1+\frac{2(k-1)}{G_{\mathrm{m}}}]$.

###### Remark 2

According to [26], to find an $m$-sparse complex signal (vector) of dimension
$n$, the number of measurements $T$ must satisfy $T>8m-2$. The dimension of
the equivalent sensing matrix ${\bm{{\Theta}}}_{{\rm{M}}}$ is $T_{1}\times
G_{\mathrm{m}}$, and the corresponding sparsity level is $(Q+1)$, thus the
pilot overhead required for the AOD estimation should satisfy $T_{1}\geq
8(Q+1)-2$.

####  IV-A2 Refining the AODs by the MLE

Due to the fact that the resolution of the DCS-SOMP algorithm is limited by
the number of the pre-discretized grids in the dictionary matrix, there exists
a mismatch between the discrete estimated angle and the real continuous angle.
Moreover, the estimation errors in AODs will greatly impact the accuracy of
estimations in the next two sub-steps. To reduce the error propagation in the
next two sub-steps, we adopt the maximum likelihood estimation (MLE) to
further refine the AOD estimations, where the first $T_{1}$ time slots are
utilized, and the AODs estimated by the DCS-SOMP algorithm are taken as the
initial values.

According to (4) and (9), the received signal over the first $T_{1}$ time
slots can be rewritten as

$\displaystyle\mathbf{y}_{t}[n]$
$\displaystyle=\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\bf
a}_{{\rm B}}(\vartheta_{{\rm r},0}){\bf a}_{{\rm R}}^{H}(\omega_{{\rm
out},0}^{{\rm a}},\omega_{{\rm out},0}^{{\rm
e}})\text{$\bm{{\Omega}}$}_{0}{\bf a}_{{\rm R}}(\omega_{{\rm in},q}^{{\rm
a}},\omega_{{\rm in},q}^{{\rm e}}){\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm
t},q}){\bf x}_{t}[n]+{\bf z}_{t}[n]$
$\displaystyle=\sum\limits_{q=0}^{Q}\check{\delta}_{q}[n]\tilde{\bf{H}}_{q}{\bf
x}_{t}[n]+{\bf z}_{t}[n]$
$\displaystyle=\tilde{\bf{X}}_{t}[n]\check{\bm{{\bf\delta}}}[n]+{\bf
z}_{t}[n],t=1,2,\cdots,T_{1},$ (19)

where
$\check{\delta}_{q}[n]\triangleq\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\bf
a}_{{\rm R}}^{H}(\omega_{{\rm out},0}^{{\rm a}},\omega_{{\rm out},0}^{{\rm
e}})\text{$\bm{{\Omega}}$}_{0}{\bf a}_{{\rm R}}(\omega_{{\rm in},q}^{{\rm
a}},\omega_{{\rm in},q}^{{\rm e}})$, $\tilde{\bf{H}}_{q}\triangleq{\bf
a}_{{\rm B}}(\vartheta_{{\rm r},0}){\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm
t},q})$,
$\check{\bm{{\bf\delta}}}[n]\triangleq[\check{\delta}_{0}[n],\check{\delta}_{1}[n],\cdots,\check{\delta}_{Q}[n]]^{T}$
and $\tilde{\bf{X}}_{t}[n]\triangleq[\tilde{\bf{H}}_{0}{\bf
x}_{t}[n],\tilde{\bf{H}}_{1}{\bf x}_{t}[n],\cdots,\tilde{\bf{H}}_{Q}{\bf
x}_{t}[n]]$. The parameters to be estimated from (19) include
$\\{\check{\bm{{\bf\delta}}}[n],\theta_{{\rm t},q}\\}$,
${q=0,1,\cdots,Q,n=1,\cdots,N}$. The joint probability density function of the
received signal in (19) can be expressed as

$\displaystyle f(\\{{\bf y}_{t}[n]\\}_{n=1,\cdots,N}^{t=1,\cdots,T_{1}})$
$\displaystyle=\prod\limits_{n=1}^{N}\prod\limits_{t=1}^{T_{1}}\frac{1}{\pi^{N_{{\rm
b}}}{\rm det}(\sigma^{2}{\bf I}_{N_{{\rm
b}}})}\exp\\{-\frac{1}{\sigma^{2}}({\bf
y}_{t}[n]-\tilde{\bf{X}}_{t}[n]\check{\bm{{\bf\delta}}}[n])^{H}({\bf
y}_{t}[n]-\tilde{\bf{X}}_{t}[n]\check{\bm{{\bf\delta}}}[n])\\}.$ (20)

By taking the natural logarithm of (20) and omitting the constant term, the
logarithmic likelihood function becomes

$\displaystyle\ln f$
$\displaystyle=-\sum\limits_{n=1}^{N}\sum\limits_{t=1}^{T_{1}}\left\|{\bf
y}_{t}[n]-\tilde{\bf{X}}_{t}[n]\check{\bm{{\bf\delta}}}[n]\right\|_{2}^{2}.$
(21)

By fixing $\theta_{{\rm t},q},q=0,1,\cdots,Q$, we take the derivative of $\ln
f$ with respect to the conjugation of
$\check{\bm{{\bf\delta}}}[n],n=1,\cdots,N$ as [27]

$\displaystyle\frac{\partial\ln f}{\partial\check{\bm{{\bf\delta}}}^{*}[n]}$
$\displaystyle=-\frac{\partial\sum\limits_{t=1}^{T_{1}}[({\bf
y}_{t}[n]-\tilde{\bf{X}}_{t}[n]\check{\bm{{\bf\delta}}}[n])^{H}({\bf
y}_{t}[n]-\tilde{\bf{X}}_{t}[n]\check{\bm{{\bf\delta}}}[n])]}{\partial\check{\bm{{\bf\delta}}}^{*}[n]}$
$\displaystyle=-\frac{\partial\sum\limits_{t=1}^{T_{1}}({\bf y}_{t}^{H}[n]{\bf
y}_{t}[n]-{\bf
y}_{t}^{H}[n]\tilde{\bf{X}}_{t}[n]\check{\bm{{\bf\delta}}}[n]-\check{\bm{{\bf\delta}}}^{H}[n]\tilde{\bf{X}}_{t}^{H}[n]{\bf
y}_{t}[n]+\check{\bm{{\bf\delta}}}^{H}[n]{\tilde{\bf
X}}_{t}^{H}[n]\tilde{\bf{X}}_{t}[n]\check{\bm{{\bf\delta}}}[n])}{\partial\check{\bm{{\bf\delta}}}^{*}[n]}$
$\displaystyle=\sum\limits_{t=1}^{T_{1}}(\tilde{\bf{X}}_{t}^{H}[n]{\bf
y}_{t}[n]-\tilde{\bf{X}}_{t}^{H}[n]\tilde{\bf{X}}_{t}[n]\check{\bm{{\bf\delta}}}[n]).$
(22)

Then, by setting the derivative in (22) to $\mathbf{0}$, we have

$\displaystyle\sum\limits_{t=1}^{T_{1}}({\tilde{\bf X}}_{t}^{H}[n]{\bf
y}_{t}[n]-{\tilde{\bf X}}_{t}^{H}[n]{\tilde{\bf
X}}_{t}[n]\check{\bm{{\bf\delta}}}[n])$ $\displaystyle=\mathbf{0}.$ (23)

By solving (23), we have

$\displaystyle\check{\bm{{\bf\delta}}}[n]$
$\displaystyle=[\sum\limits_{t=1}^{T_{1}}({\tilde{\bf
X}}_{t}^{H}[n]{\tilde{\bf
X}}_{t}[n])]^{-1}\sum\limits_{t=1}^{T_{1}}({\tilde{\bf X}}_{t}^{H}[n]{\bf
y}_{t}[n]).$ (24)

By substituting (24) into (21), the logarithmic likelihood function becomes

$\displaystyle\ln f$
$\displaystyle=-\sum\limits_{n=1}^{N}\sum\limits_{t=1}^{T_{1}}\left\|{\bf
y}_{t}[n]-{\tilde{\bf X}}_{t}[n][\sum\limits_{t=1}^{T_{1}}({\tilde{\bf
X}}_{t}^{H}[n]{\tilde{\bf
X}}_{t}[n])]^{-1}\sum\limits_{t=1}^{T_{1}}({\tilde{\bf X}}_{t}^{H}[n]{\bf
y}_{t}[n])\right\|_{2}^{2}.$ (25)

Then, the maximum likelihood estimate of
${{\bm{\theta}}_{\rm{t}}}\triangleq{[{\theta_{{\rm{t}},0}},{\theta_{{\rm{t}},1}},\cdots,{\theta_{{\rm{t}},Q}}]^{T}}$
becomes

$\displaystyle\hat{\bm{\theta}}_{\rm{t}}$
$\displaystyle=\arg\mathop{\max}\limits_{{\bm{\theta}}_{\rm{t}}}\\{-\sum\limits_{n=1}^{N}\sum\limits_{t=1}^{T_{1}}\left\|{\bf
y}_{t}[n]-{\tilde{\bf X}}_{t}[n][\sum\limits_{t=1}^{T_{1}}({\tilde{\bf
X}}_{t}^{H}[n]{\tilde{\bf
X}}_{t}[n])]^{-1}\sum\limits_{t=1}^{T_{1}}({\tilde{\bf X}}_{t}^{H}[n]{\bf
y}_{t}[n])\right\|_{2}^{2}\\}$
$\displaystyle=\arg\mathop{\min}\limits_{{\bm{\theta}}_{\rm{t}}}\sum\limits_{n=1}^{N}\sum\limits_{t=1}^{T_{1}}\left\|{\bf
y}_{t}[n]-{\tilde{\bf X}}_{t}[n][\sum\limits_{t=1}^{T_{1}}({\tilde{\bf
X}}_{t}^{H}[n]{\tilde{\bf
X}}_{t}[n])]^{-1}\sum\limits_{t=1}^{T_{1}}({\tilde{\bf X}}_{t}^{H}[n]{\bf
y}_{t}[n])\right\|_{2}^{2}.$ (26)

The maximum likelihood function in (26) can be further simplified to reduce
the computational complexity, which can be found in Appendix A. The MLE in
(26) becomes a $(Q+1)$-dimensional search problem, and we adopt the AOD
estimations obtained by the DCS-SOMP algorithm to serve as the initial values
for the MLE.

###  IV-B Estimation of the AOAs at the RIS

According to the RIS phase profile designed in Fig. 2, besides the first
$T_{1}$ time slots, the remaining $T_{2}$ time slots are divided into
$\Upsilon$ time blocks, and each time block has $V$ time slots. By
substituting the sparse representation of channels in (5a) and (5b) into (9),
the received signal at the $t$-th time slot of the $i$-th time block can be
written in the angle domain as

$\displaystyle\mathbf{y}_{t,i}[n]$ $\displaystyle={\bf H}_{{\rm
RB}}[n]\text{$\bm{{\Omega}}$}_{i}{\bf H}_{{\rm
MR}}[n]\mathbf{x}_{t}[n]+\mathbf{z}_{t,i}[n]$ $\displaystyle={\bf A}_{{\rm
B}}\bm{{\varLambda}}_{\textrm{RB}}[n]\mathbf{A}_{{\rm
R}}^{H}\text{$\bm{{\Omega}}$}_{i}{\bf A}_{{\rm
R}}\bm{{\varLambda}}_{\textrm{MR}}[n]{\bf A}_{{\rm
M}}^{H}\mathbf{x}_{t}[n]+\mathbf{z}_{t,i}[n],$ (27)

where $i=0,1,\cdots,\Upsilon$, $t=T_{1}+(i-1)V+1,\cdots,T_{1}+iV$ if $i\neq
0$, and $t=1,2,\cdots,T_{1}$ if $i=0$.

By using the AODs estimated in Section IV-A, we can obtain a
$N_{\textrm{b}}\times(Q+1)$ matrix $\bar{{\bf A}}_{{\rm M}}\triangleq$
$[{\bf a}_{{\rm M}}(\hat{\vartheta}_{{\rm t},0}),{\bf a}_{{\rm
M}}(\hat{\vartheta}_{{\rm t},1}),\cdots,{\bf a}_{{\rm
M}}(\hat{\vartheta}_{{\rm t},Q})]$ whose columns correspond to the non-zero
columns of $\bm{{\varLambda}}_{\textrm{MR}}[n]$. We define a
$N_{\textrm{r}}\times(Q+1)$ matrix $\bar{{\bf A}}_{{\rm R,in}}\triangleq[{\bf
a}_{{\rm R}}(\omega_{{\rm in},0}^{{\rm a}},\omega_{{\rm in},0}^{{\rm e}}),{\bf
a}_{{\rm R}}(\omega_{{\rm in},1}^{{\rm a}},\omega_{{\rm in},1}^{{\rm
e}}),\cdots,{\bf a}_{{\rm R}}(\omega_{{\rm in},Q}^{{\rm a}},\omega_{{\rm
in},Q}^{{\rm e}})]$, which contains the columns in the dictionary matrix ${\bf
A}_{{\rm R}}$ corresponding to the non-zero rows of
$\bm{{\varLambda}}_{\textrm{MR}}[n]$. We collect the non-zero elements of
$\bm{{\varLambda}}_{\textrm{MR}}[n]$ in a diagonal matrix as
${\bar{\bm{\varLambda}}}_{{\rm{MR}}}[n]\triangleq\textrm{diag}\\{\tilde{\delta}_{{\rm
MR,}0}[n],\tilde{\delta}_{{\rm MR,}1}[n],\cdots,\tilde{\delta}_{{\rm
MR,}Q}[n]\\}$, where $\tilde{\delta}_{{\rm MR,}q}[n]\triangleq\delta_{{\rm
MR,}q}e^{-j2\pi\tau_{{\rm MR,}q}\frac{(n-1)B}{N}}$, $q=0,1,\cdots,Q$. The
matrix
${\bar{\bm{\varLambda}}}_{\textrm{MR}}[n]\in\mathbb{C}^{(Q+1)\times(Q+1)}$ can
be also obtained by eliminating all zero rows and columns of
$\bm{{\varLambda}}_{\textrm{MR}}[n]\in\mathbb{C}^{G_{\mathrm{r}}\times
G_{\mathrm{m}}}$. The matrix $\bar{{\bf A}}_{{\rm M}}$ is estimated in Section
IV-A, while the matrices $\bar{{\bf A}}_{{\rm R,in}}$ and
${\bar{\bm{\varLambda}}}_{{\rm{MR}}}[n]$ are unknown. Then, the subchannel
matrix ${\bf H}_{{\rm MR}}[n]$ can be compactly expressed as

$\displaystyle{\bf H}_{{\rm MR}}[n]$ $\displaystyle=\bar{{\bf A}}_{{\rm
R,in}}{\bar{\bm{\varLambda}}}_{{\rm{MR}}}[n]{\bar{\bf A}}_{{\rm M}}^{H}.$ (28)

The vector ${\bar{\bf a}}_{{\rm B}}\triangleq{\bf a}_{{\rm B}}(\vartheta_{{\rm
r},0})$ is a column in the dictionary matrix ${\bf A}_{{\rm B}}$ corresponding
to the only non-zero row of $\bm{{\varLambda}}_{\textrm{RB}}[n]$, and the
vector $\bar{{\bf a}}_{{\rm R,out}}\triangleq{\bf a}_{{\rm R}}(\omega_{{\rm
out},0}^{{\rm a}},\omega_{{\rm out},0}^{{\rm e}})$ is a column in the
dictionary matrix $\mathbf{A}_{{\rm R}}$ corresponding to the only non-zero
column of $\bm{{\varLambda}}_{\textrm{RB}}[n]$. Since the AOA $\theta_{{\rm
r},0}$ at the BS and the AOD $\phi_{{\rm out},0},\psi_{{\rm out},0}$ at the
RIS can be calculated in advance, the steering vectors ${\bar{\bf a}}_{{\rm
B}}$ and $\bar{{\bf a}}_{{\rm R,out}}$ are known. The only one non-zero
element in $\bm{{\varLambda}}_{\textrm{RB}}[n]$ is $\tilde{\delta}_{{\rm
RB,}0}[n]\triangleq\delta_{{\rm RB,}0}e^{-j2\pi\tau_{{\rm
RB,}0}\frac{(n-1)B}{N}}$, which is unknown. Then, the subchannel matrix ${\bf
H}_{{\rm RB}}[n]$ can be compactly expressed as

$\displaystyle{\bf H}_{{\rm RB}}[n]$ $\displaystyle={\bar{\bf a}}_{{\rm
B}}\tilde{\delta}_{{\rm RB,}0}[n]\bar{{\bf a}}_{{\rm R,out}}^{H}.$ (29)

By substituting (28) and (29) into (27), the received signal can be rewritten
as

$\displaystyle\mathbf{y}_{t,i}[n]$ $\displaystyle={\bar{\bf a}}_{{\rm
B}}\tilde{\delta}_{{\rm RB,}0}[n]\bar{{\bf a}}_{{\rm
R,out}}^{H}\textrm{diag}\\{\mathbf{g}_{i}\\}\bar{{\bf A}}_{{\rm
R,in}}{\bar{\bm{\varLambda}}}_{{\rm{MR}}}[n]{\bar{\bf A}}_{{\rm
M}}^{H}\mathbf{x}_{t}[n]+\mathbf{z}_{t,i}[n]$ $\displaystyle={\bar{\bf
a}}_{{\rm B}}\tilde{\delta}_{{\rm RB,}0}[n]\mathbf{g}_{i}^{T}\bar{{\bf
A}}_{{\rm R,in-out}}{\bar{\bm{\varLambda}}}_{{\rm{MR}}}[n]{\bar{\bf A}}_{{\rm
M}}^{H}\mathbf{x}_{t}[n]+\mathbf{z}_{t,i}[n],$ (30)

where $\bar{{\bf A}}_{{\rm R,in-out}}\triangleq\textrm{diag}\\{\bar{{\bf
a}}_{{\rm R,out}}^{H}\\}\bar{{\bf A}}_{{\rm R,in}}$. The $q$-th column of
$\bar{{\bf A}}_{{\rm R,in-out}}$ can be written as

$\displaystyle\bar{{\bf A}}_{{\rm R,in-out}}(:,q)=$ $\displaystyle{\bf
a}_{{\rm R}}(\omega_{{\rm in},q}^{{\rm a}},\omega_{{\rm in},q}^{{\rm
e}})\circ{\bf a}_{{\rm R}}^{*}(\omega_{{\rm out},0}^{{\rm a}},\omega_{{\rm
out},0}^{{\rm e}})$ $\displaystyle=$ $\displaystyle{\bf a}_{{\rm
R}}(\triangle\omega_{q}^{{\rm a}},\triangle\omega_{q}^{{\rm e}}),$ (31)

where $\triangle\omega_{q}^{{\rm a}}\triangleq\omega_{{\rm in},q}^{{\rm
a}}-\omega_{{\rm out},0}^{{\rm a}}$ and $\triangle\omega_{q}^{{\rm
e}}\triangleq\omega_{{\rm in},q}^{{\rm e}}-\omega_{{\rm out},0}^{{\rm e}}$.

By multiplying $\bar{\bf{a}}_{{\rm B}}^{H}/N_{\textrm{b}}$ on the left hand
side of $\mathbf{y}_{t,i}[n]$, we have

$\displaystyle\check{y}_{t,i}[n]$ $\displaystyle={\bar{\bf a}}_{{\rm
B}}^{H}\mathbf{y}_{t,i}[n]/N_{\textrm{b}}$
$\displaystyle=\mathbf{g}_{i}^{T}\bar{{\bf A}}_{{\rm R,in-
out}}{\bar{\bm{\varLambda}}}_{{\rm{MRB}}}[n]{\bar{\bf A}}_{{\rm
M}}^{H}\mathbf{x}_{t}[n]+\check{z}_{t,i}[n],$ (32)

where $\check{z}_{t,i}[n]={\bar{\bf a}}_{{\rm
B}}^{H}\mathbf{z}_{t,i}[n]/N_{\textrm{b}}$, and
${\bar{\bm{\varLambda}}}_{{\rm{MRB}}}[n]\triangleq\tilde{\delta}_{{\rm
RB,}0}[n]{\bar{\bm{\varLambda}}}_{{\rm{MR}}}[n]$ is a $(Q+1)\times(Q+1)$
diagonal matrix with the $q$-th element given by
$\tilde{\delta}_{q}[n]=\tilde{\delta}_{{\rm RB,}0}[n]\tilde{\delta}_{{\rm
MR,}q}[n]$.

For $i\neq 0$, by collecting the processed signal ${\bar{\bf a}}_{{\rm
B}}^{H}\mathbf{y}_{t,i}[n]/N_{\textrm{b}}$, $t=T_{1}+(i-1)V+1,\cdots,T_{1}+iV$
over the $V$ time slots of the $i$-th time block in a row vector, we have

$\displaystyle({\check{\mathbf{y}}}_{i}[n])^{T}$
$\displaystyle=\mathbf{g}_{i}^{T}\bar{{\bf A}}_{{\rm R,in-
out}}\bar{\bm{{\varLambda}}}_{{\rm{MRB}}}[n][\bar{{\bf A}}_{{\rm
M}}^{H}\mathbf{x}_{T_{1}+(i-1)V+1}[n],\cdots,\bar{{\bf A}}_{{\rm
M}}^{H}\mathbf{x}_{T_{1}+iV}[n]]+\check{\mathbf{z}}_{i}^{T}[n]$
$\displaystyle\triangleq\mathbf{g}_{i}^{T}\bar{{\bf A}}_{{\rm R,in-
out}}\bar{\bm{{\varLambda}}}_{{\rm{MRB}}}[n]\bar{{\bf
B}}_{i}+\check{\mathbf{z}}_{i}^{T}[n],i\neq 0$ (33)

where
$\check{\mathbf{y}}_{i}[n]=[\check{y}_{T_{1}+(i-1)V+1,i}[n],\cdots,\check{y}_{T_{1}+iV,i}[n]]^{T}$,
$\check{\mathbf{z}}_{i}[n]=[\check{z}_{T_{1}+(i-1)V+1,i}[n],\cdots,\check{z}_{T_{1}+iV,i}[n]]^{T}$
and $\bar{{\bf B}}_{i}\triangleq\bar{{\bf A}}_{{\rm
M}}^{H}[\mathbf{x}_{T_{1}+(i-1)V+1}[n],\cdots,\mathbf{x}_{T_{1}+iV}[n]]$.

For $i=0$, by collecting the processed signal ${\bar{\bf a}}_{{\rm
B}}^{H}\mathbf{y}_{t,0}[n]/N_{\textrm{b}}$, $t=1,2,\cdots,T_{1}$ over the
$T_{1}$ time slots of the $0$-th time block in a row vector, we have

$\displaystyle({\check{\mathbf{y}}}_{0}[n])^{T}$
$\displaystyle=\mathbf{g}_{0}^{T}\bar{{\bf A}}_{{\rm R,in-
out}}{\bar{\bm{\varLambda}}}_{{\rm{MRB}}}[n][\bar{{\bf A}}_{{\rm
M}}^{H}\mathbf{x}_{1}[n],\bar{{\bf A}}_{{\rm
M}}^{H}\mathbf{x}_{2}[n],\cdots,\bar{{\bf A}}_{{\rm
M}}^{H}\mathbf{x}_{T_{1}}[n]]+\check{\mathbf{z}}_{0}^{T}[n]$
$\displaystyle\triangleq\mathbf{g}_{0}^{T}\bar{{\bf A}}_{{\rm R,in-
out}}{\bar{\bm{\varLambda}}}_{{\rm{MRB}}}[n]\bar{{\bf
B}}_{0}+\check{\mathbf{z}}_{0}^{T}[n],$ (34)

where
${\check{\mathbf{y}}}_{0}[n]=[\check{y}_{1,0}[n],\check{y}_{2,0}[n],\cdots,\check{y}_{T_{1},0}[n]]^{T}$,
$\check{\mathbf{z}}_{0}[n]=[\check{z}_{1,0}[n],\check{z}_{2,0}[n],\cdots,\check{z}_{T_{1},0}[n]]^{T}$,
and $\bar{{\bf B}}_{0}\triangleq\bar{{\bf A}}_{{\rm
M}}^{H}[\mathbf{x}_{1}[n],\mathbf{x}_{2}[n],\cdots,\mathbf{x}_{T_{1}}[n]]$. We
further find the inverse matrix of $\bar{{\bf B}}_{i}\in\mathbb{C}^{Q\times
V^{(i)}}$, where $V^{(i)}=V$ if $i\neq 0$ and $V^{(0)}=T_{1}$. The number of
time slots within each block $V^{(i)}$ should satisfy $V^{(i)}\geq Q$,
otherwise the right inverse matrix $\bar{{\bf B}}_{i}^{\dagger}=\bar{{\bf
B}}_{i}^{H}(\bar{{\bf B}}_{i}\bar{{\bf B}}_{i}^{H})^{-1}$ does not exist.

The signal $({\check{\mathbf{y}}}_{i}[n])^{T}\in\mathbb{C}^{1\times V^{(i)}}$
can be further processed as

$\displaystyle({\check{\mathbf{y}}}_{i}[n])^{T}\bar{{\bf B}}_{i}^{\dagger}$
$\displaystyle=\mathbf{g}_{i}^{T}\bar{{\bf A}}_{{\rm R,in-
out}}{\bar{\bm{\varLambda}}}_{{\rm{MRB}}}[n]+\check{\mathbf{z}}_{i}^{T}[n]\bar{{\bf
B}}_{i}^{\dagger}.$ (35)

By collecting the processed signal $({\check{\mathbf{y}}}_{i}[n])^{T}\bar{{\bf
B}}_{i}^{\dagger}$ of all the $(\Upsilon+1)$ time blocks, we have

$\displaystyle{\check{\mathbf{Y}}}[n]$ $\displaystyle=[(\bar{{\bf
B}}_{0}^{\dagger})^{T}{\check{\mathbf{y}}}_{0}[n],(\bar{{\bf
B}}_{1}^{\dagger})^{T}{\check{\mathbf{y}}}_{1}[n],\cdots,(\bar{{\bf
B}}_{\Upsilon}^{\dagger})^{T}{\check{\mathbf{y}}}_{\Upsilon}[n]]^{T}$
$\displaystyle=\mathbf{G}^{T}\bar{{\bf A}}_{{\rm R,in-
out}}{\bar{\bm{\varLambda}}}_{{\rm{MRB}}}[n]+{\check{\mathbf{Z}}}[n],$ (36)

where
$\mathbf{G}^{T}=[\mathbf{g}_{0},\mathbf{g}_{1},\cdots,\mathbf{g}_{\Upsilon}]^{T}$
and ${\check{\mathbf{Z}}}[n]=[(\bar{{\bf
B}}_{0}^{\dagger})^{T}\check{\mathbf{z}}_{0},(\bar{{\bf
B}}_{1}^{\dagger})^{T}\check{\mathbf{z}}_{1},\cdots,(\bar{{\bf
B}}_{\Upsilon}^{\dagger})^{T}\check{\mathbf{z}}_{\Upsilon}]^{T}$.

Then, we further rewrite
${\check{\mathbf{Y}}}[n]=[\mathbf{r}_{0}[n],\mathbf{r}_{1}[n],\cdots,\mathbf{r}_{Q}[n]]$,
${\bar{\bm{\varLambda}}}_{{\rm{MRB}}}[n]=[\bm{{\lambda}}_{0}[n],\bm{{\lambda}}_{1}[n],\cdots,\bm{{\lambda}}_{Q}[n]]$,
and
${\check{\mathbf{Z}}}[n]=[\mathbf{n}_{0}[n],\mathbf{n}_{1}[n],\cdots,\mathbf{n}_{Q}[n]]$,
and have

$\displaystyle\mathbf{r}_{q}[n]$ $\displaystyle=\mathbf{G}^{T}\bar{{\bf
A}}_{{\rm R,in-out}}\bm{{\lambda}}_{q}[n]+\mathbf{n}_{q}[n],$ (37a)
$\displaystyle=\mathbf{G}^{T}{\bf a}_{{\rm R}}(\triangle\omega_{q}^{{\rm
a}},\triangle\omega_{q}^{{\rm e}})\tilde{\delta}_{q}[n]+\mathbf{n}_{q}[n].$
(37b)

Equation (37a) can be cast as a sparse signal recovery problem, which can be
solved by the DCS-SOMP given in Algorithm 1. The dimension of the equivalent
dictionary matrix $\mathbf{G}^{T}\mathbf{A}_{{\rm R}}$ is $(\Upsilon+1)\times
G_{\mathrm{r}}$, and the corresponding sparsity level is 1, thus the pilot
overhead required should satisfy $(\Upsilon+1)\geq 8\times 1-2$. Since $Q\leq
V^{(i)}$ and $T=\sum_{i=0}^{\Upsilon}V^{(i)}$, we have $T\geq
8(Q+1)-2+\Upsilon Q\geq 13Q+6$. The phase shift matrix $\mathbf{G}$ of the RIS
should be properly designed, similar to how the pilot signal $\mathbf{X}_{1}$
was designed in Section IV-A, to produce an orthogonal dictionary matrix.

Equation (37a) gives the sparse signal recovery problem for the $q$-th path.
To obtain the cascaded cosine $\triangle\omega_{q}^{{\rm e}}$, the cascaded
sine $\triangle\omega_{q}^{{\rm a}}$ and the hybrid channel gain
$\tilde{\delta}_{q}[n]$ for $q=0,1,2,\cdots,Q$, the DCS-SOMP algorithm can be
performed by $(Q+1)$ times to recover $(Q+1)$ independent sparse variables
based on (37a). If the $k^{q}$-th column of the sensing matrix
$\mathbf{G}^{T}\mathbf{A}_{{\rm R}}$ is selected for the $q$-th path, the
column ${\bf A}_{{\rm R}}(:,k^{q})$ corresponds to the RIS array response
vectors ${\bf a}_{{\rm R}}(\triangle\omega_{q}^{{\rm
a}},\triangle\omega_{q}^{{\rm e}})$ according to (31), and the projection
coefficient is the hybrid complex channel gain $\tilde{\delta}_{q}[n]$. Based
on (7), the column ${\bf A}_{{\rm R}}(:,k^{q})$ can be written as

$\displaystyle{\bf A}_{{\rm R}}(:,k^{q})$
$\displaystyle=\mathbf{f}([-1+\frac{2(k_{\textrm{e}}^{q}-1)}{G_{\mathrm{e}}}]\frac{d_{\textrm{{\rm{RIS,e}}}}}{\lambda},N_{\textrm{e}})\otimes\mathbf{f}([-1+\frac{2(k_{\textrm{a}}^{q}-1)}{G_{\mathrm{a}}}]\frac{d_{\textrm{{\rm{RIS,a}}}}}{\lambda},N_{\textrm{a}}),$
(38)

where $k^{q}=(k_{\textrm{e}}^{q}-1)G_{\mathrm{a}}+k_{\textrm{a}}^{q}$.
Equivalently, we have $k_{\textrm{e}}^{q}=\left\lceil
k^{q}/G_{\mathrm{a}}\right\rceil$,
$k_{\textrm{a}}^{q}=k^{q}-(k_{\textrm{e}}^{q}-1)G_{\mathrm{a}}$. Since (31) is
equivalent to (38), we arrive at

$\displaystyle\cos\phi_{{\rm in},q}-\cos\phi_{{\rm out},0}$
$\displaystyle=[-1+\frac{2(k_{\textrm{e}}^{q}-1)}{G_{\mathrm{e}}}],$ (39a)
$\displaystyle\sin\psi_{{\rm in},q}\sin\phi_{{\rm in},q}-\sin\psi_{{\rm
out},0}\sin\phi_{{\rm out},0}$
$\displaystyle=[-1+\frac{2(k_{\textrm{a}}^{q}-1)}{G_{\mathrm{a}}}].$ (39b)

Since $\phi_{{\rm out},0}$ and $\psi_{{\rm out},0}$ are known, the $\phi_{{\rm
in},q}$ and $\psi_{{\rm in},q}$ can be calculated from (39).

###  IV-C Estimation of the TOAs and Channel Gains

By using the DCS-SOMP algorithm, the nonzero element at the $q$-th row and
$q$-th column of ${\bar{\bm{\varLambda}}}_{{\rm{MRB}}}[n]$ (i.e.,
$({\bar{\bm{\varLambda}}}_{{\rm{MRB}}}[n])_{q,q}$) is the estimated value for
the complex channel gain $\tilde{\delta}_{q}[n]$ of the $q$-th cascaded path.
Since $\tilde{\delta}_{q}[n]=\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}$, by
stacking the $\tilde{\delta}_{q}[n]$ of all the $N$ subcarriers in a vector,
we have
$\tilde{\bm{{\delta}}}_{q}=[\tilde{\delta}_{q}[1],\tilde{\delta}_{q}[2],\cdots,\tilde{\delta}_{q}[N]]^{T}$
and
$\tilde{\bm{{\delta}}}_{q}=\delta_{q}\mathbf{t}(\upsilon_{q}),q=0,1,\cdots,Q$,
where
$\mathbf{t}(\upsilon_{q})=[1,e^{-j2\pi\upsilon_{q}},\cdots,e^{-j2\pi(N-1)\upsilon_{q}}]^{T}$
and $\upsilon_{q}\triangleq\tau_{q}\frac{B}{N}$. We assume that $\upsilon_{q}$
satisfies $0<\upsilon_{q}<1$.

By defining a normalized DFT matrix $\mathbf{U}_{N}$ with
$[\mathbf{U}_{N}]_{m^{\prime},m}=\frac{1}{\sqrt{N}}e^{-j\frac{2\pi}{N}(m^{\prime}-1)(m-1)}$,
we have
$\mathbf{z}_{q}\triangleq\mathbf{U}_{N}^{H}\tilde{\bm{{\delta}}}_{q}=\delta_{q}\mathbf{U}_{N}^{H}\mathbf{t}(\upsilon_{q})$,
and the $m$-th element of $\mathbf{z}_{q}$ can be written as

$\displaystyle[\mathbf{z}_{q}]_{m}$
$\displaystyle=\frac{\delta_{q}}{\sqrt{N}}\sum_{m^{\prime}=1}^{N}e^{j\frac{2\pi}{N}(m^{\prime}-1)(m-1)}e^{-j2\pi(m^{\prime}-1)\upsilon_{q}}$
$\displaystyle=\frac{\delta_{q}}{\sqrt{N}}e^{-j2\pi\\{[\upsilon_{q}-\frac{(m-1)}{N}]\frac{N-1}{2}\\}}\frac{\sin\\{2\pi[\upsilon_{q}-\frac{(m-1)}{N}]\frac{N}{2}\\}}{\sin\\{2\pi[\upsilon_{q}-\frac{(m-1)}{N}]\frac{1}{2}\\}}.$
(40)

From (40), it is noted that if $N\rightarrow\infty$, there always exists some
integers $m_{q}=N\upsilon_{q}+1=\tau_{q}B+1$ such that
$[\mathbf{z}_{q}]_{m}=\frac{\delta_{q}}{\sqrt{N}}$ while the other elements of
$\mathbf{z}_{q}$ are all zeros. Thus, we can find $m_{q}$ by

$\displaystyle m_{q}$
$\displaystyle=\arg\underset{m}{\max}\left|\mathbf{e}_{m}^{T}\mathbf{U}_{N}^{H}\hat{\tilde{\bm{{\delta}}}}_{q}\right|,$
(41)

where $\mathbf{e}_{m}$ is a vector comprising all zeros, except a 1 in the
$m$-th entry, and $\hat{\tilde{\bm{{\delta}}}}_{q}$ is the estimated value for
$\tilde{\bm{{\delta}}}_{q}$. Then we obtain the coarsely-estimated TOA as
$\hat{\tau}_{q}=(m_{q}-1)/B$.

The accuracy of TOA estimations obtained above are limited by the resolution
of the DFT, thus we further refine the TOA estimation by utilizing the
rotation property. We define a rotation matrix as
$\bm{{\Phi}}_{N}(\triangle\tau_{q})=\textrm{diag}\\{1,e^{j2\pi\triangle\tau_{q}\frac{B}{N}},\cdots,e^{j2\pi(N-1)\triangle\tau_{q}\frac{B}{N}}\\}$,
then we have

$\displaystyle[\mathbf{U}_{N}^{H}\bm{{\Phi}}_{N}^{H}(\triangle\tau_{q})\hat{\tilde{\bm{{\delta}}}}_{q}]_{m}$
$\displaystyle=\frac{\delta_{q}}{\sqrt{N}}e^{-j2\pi[\tau_{q}\frac{B}{N}+\triangle\tau_{q}\frac{B}{N}-\frac{(m-1)}{N}]\frac{N-1}{2}\\}}\frac{\sin\\{2\pi[\tau_{q}\frac{B}{N}+\triangle\tau_{q}\frac{B}{N}-\frac{(m-1)}{N}]\frac{N}{2}\\}}{\sin\\{2\pi[\tau_{q}\frac{B}{N}+\triangle\tau_{q}\frac{B}{N}-\frac{(m-1)}{N}]\frac{1}{2}\\}}.$
(42)

It is noted that the power will concentrate on the $m_{q}$th row without power
leakage when the time rotation parameters satisfy
$\triangle\tau_{q}=(m_{q}-1)/{B}-\tau_{q}$. The rotation parameter
$\triangle\tau_{q}$ can be found by a one-dimensional search by solving the
following problem

$\displaystyle\triangle\tau_{q}$
$\displaystyle=\arg\underset{\triangle\tau_{q}\in[-\frac{1}{2B},\frac{1}{2B}]}{\max}\left|\mathbf{U}_{m_{q},:}^{H}\bm{{\Phi}}_{N}^{H}(\triangle\tau_{q})\hat{\tilde{\bm{{\delta}}}}_{q}\right|.$
(43)

When $\triangle\tau_{q}$ is obtained from (43), the TOA $\tau_{q}$ can be
estimated by

$\displaystyle\hat{\tau}_{q}$
$\displaystyle=\frac{(m_{q}-1)}{B}-\triangle\tau_{q}.$ (44)

Then, according to the least squares (LS) estimation, the complex channel gain
can be estimated by

$\displaystyle\hat{\delta}_{q}$
$\displaystyle=\frac{\mathbf{t}^{H}(\hat{\upsilon}_{q})\hat{\tilde{\bm{{\delta}}}}_{q}}{N},$
(45)

where $\hat{\upsilon}_{q}=\hat{\tau}_{q}\frac{B}{N}$.

##  V Refining Estimation of Channel Parameters

Although the AOD estimations are refined by the MLE after using the DCS-SOMP
algorithm, only the first $T_{1}$ time slots are utilized for estimation. The
estimations for AOD can be further refined by exploiting all the $T$ time
slots. Although the TOA estimations are refined based on the time rotation
operation after the DFT approach, these estimations rely on the hybrid complex
channel gains in $\tilde{\bm{{\delta}}}_{q}$ which are obtained from the AOA
estimation at the RIS by the DCS-SOMP algorithm. Thus, the TOA estimations
should be further refined due to the lack of high accuracy of the coarsely-
estimated hybrid channel gains. Moreover, since the accuracy of azimuth AOA
$\psi_{\mathrm{in},q}$ and elevation AOA $\phi_{{\rm in},q}$ at the RIS is
limited by the DCS-SOMP algorithm, the AOAs at the RIS should also be further
refined. Thus, we consider to exploit the MLE to jointly refine all these
channel parameters by utilizing the received signal over all the $T$ time
slots.

###  V-A MLE of All Channel Parameters

Similar to the manipulations in (30), the channel matrix in (4) can be
rewritten as

$\displaystyle{\bf H}_{t}[n]$
$\displaystyle=\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\bf
a}_{{\rm B}}(\vartheta_{{\rm r},0})\mathbf{g}_{t}^{T}\textrm{diag}\\{{\bf
a}_{{\rm R}}^{H}(\omega_{{\rm out},0}^{{\rm a}},\omega_{{\rm out},0}^{{\rm
e}})\\}{\bf a}_{{\rm R}}(\omega_{{\rm in},q}^{{\rm a}},\omega_{{\rm
in},q}^{{\rm e}}){\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm t},q})$
$\displaystyle=\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}\mathbf{g}_{t}^{T}{\bf
a}_{{\rm R}}(\triangle\omega_{q}^{{\rm a}},\triangle\omega_{q}^{{\rm e}}){\bf
a}_{{\rm B}}(\vartheta_{{\rm r},0}){\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm
t},q}).$ (46)

By substituting (46) into (9), the received signal can be written as

$\displaystyle{\bf y}_{t}[n]$
$\displaystyle=\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}\mathbf{g}_{t}^{T}{\bf
a}_{{\rm R}}(\triangle\omega_{q}^{{\rm a}},\triangle\omega_{q}^{{\rm
e}}){\tilde{{\bf H}}}_{q}{\bf x}_{t}[n]+{\bf z}_{t}[n],$ (47)

where ${\tilde{{\bf H}}}_{q}\triangleq{\bf a}_{{\rm B}}(\vartheta_{{\rm
r},0}){\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm t},q})$.

By stacking the received signal ${\bf y}_{t}[n]$ at different time slots
(i.e., $t=1,\cdots,T$) in a columnwise order, we have

${\bf
Y}[n]=\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q}+{\bf Z}[n],$ (48)

where ${\bf Y}[n]\triangleq[{\bf y}_{1}[n],\cdots,{\bf y}_{T}[n]]$, ${\bf
X}[n]\triangleq[{\bf x}_{1}[n],\cdots,{\bf x}_{T}[n]]$, ${\bf
Z}[n]\triangleq[{\bf z}_{1}[n],\cdots,{\bf z}_{T}[n]]$ and
$\bm{{\Sigma}}_{q}\triangleq{\rm diag}\\{\mathbf{g}_{1}^{T}{\bf a}_{{\rm
R}}(\triangle\omega_{q}^{{\rm a}},\triangle\omega_{q}^{{\rm
e}}),\cdots,\mathbf{g}_{T}^{T}{\bf a}_{{\rm R}}(\triangle\omega_{q}^{{\rm
a}},\triangle\omega_{q}^{{\rm e}})\\}$. Then we vectorize the ${\bf Y}[n]$ in
(48), and have

$\displaystyle{\bf y}[n]$
$\displaystyle=\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}(\bm{{\Sigma}}_{q}^{T}\otimes{\tilde{{\bf
H}}}_{q}){\bf x}[n]+{\bf z}[n],$ (49)

where ${\bf y}[n]\triangleq{\rm vec}({\bf Y}[n])$, ${\bf x}[n]\triangleq{\rm
vec}({\bf X}[n])$ and ${\bf z}[n]\triangleq{\rm vec}({\bf Z}[n])$. The (49) is
obtained from (48) by using the identity ${{\rm vec}({\bf ABC})=({\bf
C}^{T}\otimes{\bf A}){\rm vec}({\bf B})}$. Then, based on the received signal
in (49), the likelihood function of the random vector $\\{{\bf
y}[n]\\}_{n=1,2,\cdots,N}$ conditioned on $\bm{{\bf\eta}}$ can be written as

$\displaystyle f(\\{{\bf y}[n]\\}_{n=1,2,\cdots,N}|\bm{{\bf\eta}})$
$\displaystyle=\prod\limits_{n=1}^{N}\frac{1}{\pi^{TN_{{\rm b}}}{\rm
det}(\sigma^{2}{\bf I}_{TN_{{\rm b}}})}\exp\left[-({\bf
y}[n]-\bm{\mu}[n])^{H}{\bf\bm{{\Xi}}}^{-1}({\bf y}[n]-\bm{\mu}[n])\right].$
(50)

where
$\bm{\mu}[n]\triangleq\sum_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}(\bm{{\Sigma}}_{q}^{T}\otimes{\tilde{{\bf
H}}}_{q}){\bf x}[n]$ and ${\bf\bm{{\Xi}}}=\sigma^{2}{\bf I}_{TN_{{\rm b}}}$.

By taking the natural logarithm of (50) and omitting the constant term, we
have the log-likelihood functon $\Lambda(\bm{{\bf\eta}};\\{{\bf
y}[n]\\}_{n=1,2,\cdots,N})$ of $\bm{{\bf\eta}}$ given the observed signal
${\bf y}[n],n=1,2,\cdots,N$. After some mathmatical simplifications, the
$\Lambda(\bm{{\bf\eta}};\\{{\bf y}[n]\\}_{n=1,2,\cdots,N})$ can be finally
expressed as

$\displaystyle\Lambda(\bm{{\bf\eta}};\\{{\bf
Y}[n]\\}_{n=1,2,\cdots,N})=2{\mathfrak{Re}}\left[\sum\limits_{q=0}^{Q}{\bf
a}_{{\rm M}}^{H}(\vartheta_{{\rm
t},q})\left(\sum\limits_{n=1}^{N}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\bf
X}[n]\bm{{\Sigma}}_{q}{\bf Y}^{H}[n]\right){\bf a}_{{\rm B}}(\vartheta_{{\rm
r},0})\right]$ $\displaystyle\quad-
N_{\mathrm{b}}\sum\limits_{q_{1}=0}^{Q}\sum\limits_{q_{2}=0}^{Q}\left[\delta_{q_{1}}^{*}\delta_{q_{2}}{\bf
a}_{{\rm M}}^{H}(\vartheta_{{\rm
t},q_{2}})\left(\sum\limits_{n=1}^{N}e^{j2\pi(\tau_{q_{1}}-\tau_{q_{2}})\frac{(n-1)B}{N}}{\bf
X}[n]\bm{{\Sigma}}_{q_{2}}\bm{{\Sigma}}_{q_{1}}^{H}{\bf X}^{H}[n]\right){\bf
a}_{{\rm M}}(\vartheta_{{\rm t},q_{1}})\right].$ (51)

The detailed derivation process is given in (115) of Appendix B-A. Then the
MLE for the channel parameters in $\bm{{\bf\eta}}$ can be formulated by
$\hat{\bm{{\bf\eta}}}=\arg\underset{\bm{{\bf\eta}}}{\max}{\>}\Lambda(\bm{{\bf\eta}};\\{{\bf
Y}[n]\\}_{n=1,2,\cdots,N})$. Obviously, the expression of
$\Lambda(\bm{{\bf\eta}};\\{{\bf Y}[n]\\}_{n=1,2,\cdots,N})$ is complex, and
the high dimensional nonlinear optimization for
$\bm{{\bf\eta}}\in\mathbb{C}^{6(Q+1)}$ is computationally prohibitive. Since
the MLE is intractable, we propose to exploit the space alternating
generalized expectation (SAGE) algorithm[28] to update the parameters
sequentially. The SAGE algorithm has been widely utilized in the estimation
problem for channel parameters [22], where the high-dimensional optimization
process is circumvented by several sequentially-performed low-dimensional
maximization procedures.

###  V-B Solving the MLE by SAGE algorithm

The derivation of the SAGE algorithm relies on two key notions of the complete
(unobservable) and incomplete (observable) data. The received signal ${\bf
y}_{t}[n]$ in (47) is identified as the incomplete data. The individual
received signals ${\bf y}_{q,t}[n]$ from the $q$-th path, $q=0,1,\cdots,Q$,
constitute a set of complete data, which can be written as

$\displaystyle{\bf y}_{q,t}[n]$
$\displaystyle=\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}\mathbf{g}_{t}^{T}{\bf
a}_{{\rm R}}(\triangle\omega_{q}^{{\rm a}},\triangle\omega_{q}^{{\rm
e}}){\tilde{{\bf H}}}_{q}{\bf x}_{t}[n]+{\bf z}_{q,t}[n],q=0,1,\cdots,Q,$ (52)

where ${\bf z}_{q,t}[n]\in\mathbb{C}^{N_{\mathrm{b}}}$ represents the noise of
the $q$-th path satisfying ${\bf
z}_{q,t}[n]\sim\mathcal{CN}(0,\sigma_{q}^{2}\mathbf{I}_{N_{{\rm b}}})$, and
$\sigma_{q}^{2}$ denotes the noise power. The mutually-independent ${\bf
z}_{0,t}[n],{\bf z}_{1,t}[n],\cdots,{\bf z}_{Q,t}[n]$ form a decomposition of
${\bf z}_{t}[n]$ in (47), thus we have
$\sum_{q=0}^{Q}\sigma_{q}^{2}=\sigma^{2}$. The incomplete data ${\bf
y}_{t}[n]$ are related to the complete data $\\{{\bf
y}_{q,t}[n]\\}_{q=0,1,\cdots,Q}$ by ${\bf y}_{t}[n]=\sum_{q=0}^{Q}{\bf
y}_{q,t}[n]$. Since ${\bf y}_{0,t}[n],{\bf y}_{1,t}[n],\cdots,{\bf
y}_{Q,t}[n]$ are independent, the signals ${\bf y}_{q^{\prime},t}[n]$,
$q^{\prime}\neq q$, are irrelevant for the estimation of $\bm{{\bf\eta}}_{q}$.
If the complete data can be observed, the MLE of the parameters of the $q$-th
path in $\bm{{\bf\eta}}_{q}$ can be obtained from the observed ${\bf
y}_{q,t}[n]$ through $(\hat{\bm{{\bf\eta}}}_{q})_{\textrm{ML}}({\bf
y}_{q,t}[n])=\arg\max L(\bm{{\bf\eta}}_{q};{\bf y}_{q,t}[n])$, where
$L(\bm{{\bf\eta}}_{q};{\bf y}_{q,t}[n])$ denotes the likelihood function of
$\bm{{\bf\eta}}_{q}$ given an observation ${\bf y}_{q,t}[n]$, and
$(\hat{\bm{{\bf\eta}}}_{q})_{\textrm{ML}}({\bf y}_{q,t}[n])$ means the MLE of
$\bm{{\bf\eta}}_{q}$ is a function of the observation ${\bf y}_{q,t}[n]$.
Since the ${\bf y}_{q,t}[n]$ cannot be observed, it is usually estimated from
the incomplete data ${\bf y}_{t}[n]$ and a previous estimation of
$\bm{{\bf\eta}}_{q}$. Thus, in the SAGE algorithms, we first formulate the MLE
problem of a single path, then estimate the complete data ${\bf y}_{q,t}[n]$,
and finally solve the MLE problem by utilizing the estimated ${\bf
y}_{q,t}[n]$. The three key steps are illustrated as follows.

####  V-B1 MLE of a Single Path

We first derive the log-likelihood function of the channel parameters
$\bm{{\bf\eta}}_{q}$ assuming that ${\bf y}_{q,t}[n]$ is observable. By
stacking ${\bf y}_{q,t}[n]$ over $T$ time slots in a columnwise order, we can
obtain ${\bf Y}_{q}[n]\triangleq[{\bf y}_{q,1}[n],\cdots,{\bf y}_{q,T}[n]]$.
By using (52) and similar to (48), ${\bf Y}_{q}[n]$ can be written as

${\bf Y}_{q}[n]=\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q}+{\bf Z}_{q}[n],$ (53)

where ${\bf Z}_{q}[n]\triangleq[{\bf z}_{q,1}[n],\cdots,{\bf z}_{q,T}[n]]$. We
vectorize the equality in (53) to obtain

$\displaystyle{\bf y}_{q}[n]$
$\displaystyle=\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}(\bm{{\Sigma}}_{q}^{T}\otimes{\tilde{{\bf
H}}}_{q}){\bf x}[n]+{\bf z}_{q}[n],$ (54)

where ${\bf y}_{q}[n]\triangleq{\rm vec}({\bf Y}_{q}[n])$ and ${\bf
z}_{q}[n]\triangleq\mathrm{vec}({\bf Z}_{q}[n])$.

According to (54), the joint probability density function of ${\bf y}_{q}[n]$
over $N$ subcarriers is similar to that in (50). And the corresponding log-
likelihood function of $\bm{{\bf\eta}}_{q}$ for the observations $\\{{\bf
y}_{q}[n]\\}_{n=1,\cdots,N}$ can be written as

$\displaystyle L(\bm{{\bf\eta}}_{q};\\{{\bf y}_{q}[n]\\}_{n=1,\cdots,N})$
$\displaystyle\propto-\sum\limits_{n=1}^{N}\left\|{\bf
y}_{q}[n]-\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}(\bm{{\Sigma}}_{q}^{T}\otimes\tilde{{\bf{H}}}_{q}){\bf
x}[n]\right\|_{2}^{2}$ $\displaystyle=-\sum\limits_{n=1}^{N}\left\|{\bf
Y}_{q}[n]-\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}\tilde{{\bf{H}}}_{q}{\bf
X}[n]\bm{{\Sigma}}_{q}\right\|_{F}^{2}.$ (55)

where $\propto$ denotes equality up to irrelevant constants. By further
omitting the constant terms in (55), the $L(\bm{{\bf\eta}}_{q};\\{{\bf
y}_{q}[n]\\}_{n=1,\cdots,N})$ can be simplified into

$\displaystyle L(\bm{{\bf\eta}}_{q};\\{{\bf Y}_{q}[n]\\}_{n=1,\cdots,N})$
$\displaystyle=2{\mathfrak{Re}}\left[\delta_{q}{\bf a}_{{\rm
M}}^{H}(\vartheta_{{\rm
t},q})\left(\sum\limits_{n=1}^{N}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\bf
X}[n]\bm{{\Sigma}}_{q}{\bf Y}_{q}^{H}[n]\right){\bf a}_{{\rm
B}}(\vartheta_{{\rm r},0})\right]$ $\displaystyle\quad-N_{{\rm
b}}\delta_{q}^{*}\delta_{q}{\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm
t},q})\left(\sum\limits_{n=1}^{N}{\bf
X}[n]\bm{{\Sigma}}_{q}\bm{{\Sigma}}_{q}^{H}{\bf X}^{H}[n]\right){\bf a}_{{\rm
M}}(\vartheta_{{\rm t},q}).$ (56)

The detailed derivation process for (56) is similar as that in Appendix B-A.
Or equivalently, the (56) is a special case of (51), thus can be obtained from
(51) directly.

Then the MLE of $\bm{{\bf\eta}}_{q}$ can be formulated as

$\displaystyle(\hat{\bm{\eta}}_{q})_{\textrm{ML}}(\\{{\bf
Y}_{q}[n]\\}_{n=1,\cdots,N})$
$\displaystyle=\arg\underset{\bm{{\bf\eta}}_{q}}{\max}{\>}L(\bm{{\bf\eta}}_{q};\\{{\bf
Y}_{q}[n]\\}_{n=1,\cdots,N}).$ (57)

####  V-B2 Estimation of the Complete Data

From (57), we find that to obtain the MLE of $\bm{{\bf\eta}}_{q}$, the
received signal ${\bf Y}_{q}[n]$ (equivalently the ${\bf y}_{q,t}[n]$) should
be given. Since ${\bf y}_{q,t}[n]$ is unobservable, it can be estimated based
on the observation ${\bf y}_{t}[n]$ of the incomplete data and the previous
estimation of $\bm{{\bf\eta}}$. Specifically, in the $k$-th iteration, the
received signal of the $q$-th path can be estimated by its conditional
expectation given the incomplete data ${\bf y}_{t}[n]$ and the last estimation
of ${\hat{\bm{\eta}}}^{(k-1)}$ as

$\displaystyle{\hat{{\bf y}}}_{q,t}^{(k)}[n]$ $\displaystyle=\mathbb{E}({\bf
y}_{q,t}[n]\left|{\bf y}_{t}[n],{\hat{\bm{\eta}}}^{(k-1)}\right.).$ (58)

In the SAGE algorithm, (58) can be explicitly calculated by[29]

$\displaystyle{\hat{{\bf y}}}_{q,t}^{(k)}[n]=$ $\displaystyle{\bf
y}_{t}[n]-\sum\limits_{q^{\prime}=0}^{q-1}\delta_{q^{\prime}}^{(k)}e^{-j2\pi\tau_{q^{\prime}}^{(k)}\frac{(n-1)B}{N}}{\bf
g}^{T}_{t}{\bf a}_{{\rm R}}((\triangle\omega_{q^{\prime}}^{{\rm
a}})^{(k)},(\triangle\omega_{q^{\prime}}^{{\rm e}})^{(k)}){\tilde{{\bf
H}}}_{q^{\prime}}^{(k)}{\bf x}_{t}[n]$
$\displaystyle-\sum\limits_{q^{\prime\prime}=q+1}^{Q}\delta_{q^{\prime\prime}}^{(k-1)}e^{-j2\pi\tau_{q^{\prime\prime}}^{(k-1)}\frac{(n-1)B}{N}}{\bf
g}^{T}_{t}{\bf a}_{{\rm R}}((\triangle\omega_{q^{\prime\prime}}^{{\rm
a}})^{(k-1)},(\triangle\omega_{q^{\prime\prime}}^{{\rm
e}})^{(k-1)}){\tilde{{\bf H}}}_{q^{\prime\prime}}^{(k-1)}{\bf x}_{t}[n],$ (59)

where ${{\tilde{{\bf H}}}_{q^{\prime}}^{(k)}={\bf a}_{{\rm B}}(\vartheta_{{\rm
r},0}^{(k)}){\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm t},q^{\prime}}^{(k)})}$ and
${\tilde{{\bf H}}}_{q^{\prime\prime}}^{(k-1)}={\bf a}_{{\rm
B}}(\vartheta_{{\rm r},0}^{(k-1)}){\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm
t},q^{\prime\prime}}^{(k-1)})$.

####  V-B3 Solving the MLE by Coordinate Ascent Method

With the estimated ${\hat{{\bf y}}}_{q,t}^{(k)}[n]$ above, we have $\hat{{\bf
Y}}_{q}^{(k)}[n]=[{\hat{{\bf y}}}_{q,1}^{(k)}[n],{\hat{{\bf
y}}}_{q,2}^{(k)}[n],\cdots,{\hat{{\bf y}}}_{q,T}^{(k)}[n]]$. By substituting
the estimated ${\hat{{\bf Y}}}_{q}^{(k)}[n]$ into (57), the MLE of
$\bm{{\bf\eta}}_{q}$ can be rewritten as

$\displaystyle(\hat{\bm{\eta}}_{q})_{\textrm{ML}}(\\{{\hat{{\bf
Y}}}_{q}^{(k)}[n]\\}_{n=1,\cdots,N})$
$\displaystyle=\arg\underset{\bm{{\bf\eta}}_{q}}{\max}{\>}L(\bm{{\bf\eta}}_{q};\\{{\hat{{\bf
Y}}}_{q}^{(k)}[n]\\}_{n=1,\cdots,N}).$ (60)

When solving the MLE problem in (60), we can first find the closed-form
solution for the channel gain $\delta_{q}$ in
${\bm{\eta}}_{q}=[\tau_{q},\delta_{q,\textrm{R}},\delta_{q,\textrm{I}},\theta_{{\rm
t},q},\phi_{{\rm in},q},\psi_{{\rm in},q}]^{T}$. Specifically, by fixing
$\\{\tau_{q},\theta_{{\rm t},q},\phi_{{\rm in},q},\psi_{{\rm in},q}\\}$, we
take the derivative of $L(\bm{{\bf\eta}}_{q};\\{{\hat{{\bf
Y}}}_{q}^{(k)}[n]\\}_{n=1,\cdots,N})$ with respect to the conjugate of
$\delta_{q}$ as

$\displaystyle\frac{\partial L(\bm{{\bf\eta}}_{q};\\{{\hat{{\bf
Y}}}_{q}^{(k)}[n]\\}_{n=1,\cdots,N})}{\partial\delta_{q}^{*}}=-\frac{\partial\sum\limits_{n=1}^{N}\left\|{\bf\hat{Y}}_{q}^{(k)}[n]-\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q}\right\|_{F}^{2}}{\partial\delta_{q}^{*}},$
(61a)
$\displaystyle\overset{(a)}{=}\sum\limits_{n=1}^{N}\left\\{e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}{\rm
tr}[({\tilde{{\bf H}}}_{\mathit{q}}{\bf X}[n]\bm{{\Sigma}}_{q})^{H}{\hat{{\bf
Y}}}_{q}^{(k)}[n]]\right\\}-\delta_{q}\sum\limits_{n=1}^{N}{\rm
tr}[({\tilde{{\bf H}}}_{\mathit{q}}{\bf
X}[n]\bm{{\Sigma}}_{q})^{H}({\tilde{{\bf H}}}_{\mathit{q}}{\bf
X}[n]\bm{{\Sigma}}_{q})],$ (61b)

where the (61a) is due to (56) $\propto$ (55), and the detailed derivation of
step $(a)$ can be found in (116) of Appendix B-B. By setting the derivative in
(61b) to zero, we can build an equation [27], and the solution for
$\delta_{q}$ can be readily obtained as

$\displaystyle\delta_{q}$
$\displaystyle=\frac{\sum\limits_{n=1}^{N}\left\\{e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}{\rm
tr}[({\tilde{{\bf H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q})^{H}{\hat{{\bf
Y}}}_{q}^{(k)}[n]]\right\\}}{\sum\limits_{n=1}^{N}{\rm tr}[({\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q})^{H}({\tilde{{\bf H}}}_{q}{\bf
X}[n]\bm{{\Sigma}}_{q})]}$ $\displaystyle\overset{(b)}{=}\frac{{\bf a}_{{\rm
B}}^{H}(\vartheta_{{\rm
r,0}})\left\\{\sum\limits_{n=1}^{N}(e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}{\hat{{\bf
Y}}}_{q}^{(k)}[n]\bm{{\Sigma}}_{q}^{H}{\bf X}^{H}[n])\right\\}{\bf a}_{{\rm
M}}(\vartheta_{{\rm t,\mathit{q}}})}{N_{{\rm b}}{\bf a}_{{\rm
M}}^{H}(\vartheta_{{\rm t,\mathit{q}}})\left\\{\sum\limits_{n=1}^{N}({\bf
X}[n]\bm{{\Sigma}}_{q}\bm{{\Sigma}}_{q}^{H}{\bf X}^{H}[n])\right\\}{\bf
a}_{{\rm M}}(\vartheta_{{\rm t,\mathit{q}}})},$ (62)

where the detailed simplification of step $(b)$ in (62) can be found in (117)
of Appendix B-C. By substituting (62) into $L(\bm{{\bf\eta}}_{q};\\{{\hat{{\bf
Y}}}_{q}^{(k)}[n]\\}_{n=1,\cdots,N})$, we obtain the likelihood function of
$\bar{\bm{\eta}}_{q}=[\tau_{q},\theta_{{\rm t},q},\phi_{{\rm in},q},\psi_{{\rm
in},q}]^{T}$ as follows:

$\displaystyle F(\bar{\bm{\eta}}_{q};\\{{\hat{{\bf
Y}}}_{q}^{(k)}[n]\\}_{n=1,\cdots,N})$
$\displaystyle\overset{(c)}{=}\frac{\left|{\bf a}_{{\rm
B}}^{H}(\vartheta_{{\rm
r},0})\left\\{\sum\limits_{n=1}^{N}(e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}{\hat{{\bf
Y}}}_{q}^{(k)}[n]\bm{{\Sigma}}_{q}^{H}{\bf X}^{H}[n])\right\\}{\bf a}_{{\rm
M}}(\vartheta_{{\rm t},q})\right|^{2}}{N_{{\rm b}}{\bf a}_{{\rm
M}}^{H}(\vartheta_{{\rm t},q})\left\\{\sum\limits_{n=1}^{N}({\bf
X}[n]\bm{{\Sigma}}_{q}\bm{{\Sigma}}_{q}^{H}{\bf X}^{H}[n])\right\\}{\bf
a}_{{\rm M}}(\vartheta_{{\rm t},q})},$ (63)

where the detailed simplification of step $(c)$ can be found in (119) of
Appendix B-D. Thus the maximum likelihood estimation of $\bar{\bm{\eta}}_{q}$
becomes

$\hat{\bar{\bm{\eta}}}_{q}^{(k)}=\arg\mathop{\max}\limits_{\bar{\bm{\eta}}_{q}}F(\bar{\bm{\eta}}_{q};\\{{\hat{{\bf
Y}}}_{q}^{(k)}[n]\\}_{n=1,\cdots,N}).$ (64)

To further reduce the complexity in (64), we aim at replacing the multiple-
dimension optimization to compute the MLE of the parameters in
$\bar{\bm{\eta}}_{q}$ by multiple separate one-dimensional optimizations,
where each parameter in $\bar{\bm{\eta}}_{q}$ is estimated individually. By
denoting the objective function of (64) by $F(\bar{\bm{\eta}}_{q})$ for
notation convenience, the coordinate-wise updating procedure of the parameter
estimates of the $q$-th path can be performed sequentially as follows:

$\displaystyle\hat{\tau}_{q}^{(k)}$
$\displaystyle=\arg\mathop{\max}\limits_{\tau_{q}}F(\tau_{q},\hat{\theta}_{{\rm
t},q}^{(k-1)},\hat{\phi}_{{\rm in},q}^{(k-1)},\hat{\psi}_{{\rm
in},q}^{(k-1)}),$ (65a) $\displaystyle\hat{\theta}_{{\rm t},q}^{(k)}$
$\displaystyle=\arg\mathop{\max}\limits_{\theta_{\mathrm{t},q}}F(\hat{\tau}_{q}^{(k)},\theta_{\mathrm{t},q},\hat{\phi}_{{\rm
in},q}^{(k-1)},\hat{\psi}_{{\rm in},q}^{(k-1)}),$ (65b)
$\displaystyle\hat{\phi}_{{\rm in},q}^{(k)}$
$\displaystyle=\arg\mathop{\max}\limits_{\phi_{{\rm
in},q}}F(\hat{\tau}_{q}^{(k)},\hat{\theta}_{{\rm t},q}^{(k)},\phi_{{\rm
in},q},\hat{\psi}_{{\rm in},q}^{(k-1)}),$ (65c) $\displaystyle\hat{\psi}_{{\rm
in},q}^{(k)}$ $\displaystyle=\arg\mathop{\max}\limits_{\psi_{{\rm
in},q}}F(\hat{\tau}_{q}^{(k)},\hat{\theta}_{{\rm t},q}^{(k)},\hat{\phi}_{{\rm
in},q}^{(k)},\psi_{{\rm in},q}),$ (65d) $\displaystyle\hat{\delta}_{q}^{(k)}$
$\displaystyle=\frac{{\bf a}_{{\rm
B}}^{H}(\vartheta_{\mathrm{r},0})\left\\{\sum\limits_{n=1}^{N}(e^{j2\pi\hat{\tau}_{q}^{(k)}\frac{(n-1)B}{N}}\hat{{\bf{Y}}}_{q}^{(k)}[n]({\bf\Sigma}_{q}^{(k)})^{H}{\bf
X}^{H}[n])\right\\}{\bf a}_{{\rm M}}(\hat{\vartheta}_{{\rm
t},q}^{(k)})}{N_{{\rm b}}{\bf a}_{{\rm M}}^{H}(\hat{\vartheta}_{{\rm
t},q}^{(k)})\left\\{\sum\limits_{n=1}^{N}({\bf
X}[n]({\bf\Sigma}_{q}^{(k)})({\bf\Sigma}_{q}^{(k)})^{H}{\bf
X}^{H}[n])\right\\}{\bf a}_{{\rm M}}(\hat{\vartheta}_{{\rm t},q}^{(k)})},$
(65e)

where the $\hat{\delta}_{q}^{(k)}$ in (65e) is obtained by sustituting the
estimated values in (65a), (65b), (65c) and (65d) into (62).

The SAGE algorithm first allows the splitting of the joint optimization for
superimposed $(Q+1)$ paths into separate optimization for a single path, then
for each single path, not all parameters but only a single parameter related
to the path is estimated while keeping the estimates of the other parameters
fixed. The iteration process in the SAGE algorithm is terminated when the
distance between the estimated $\bm{{\bf\eta}}$ returned by the SAGE algorithm
at two consecutive iteration steps is below a predefined threshold or when the
value sequence of ML $\Lambda(\bm{{\bf\eta}};\\{{\bf
Y}[n]\\}_{n=1,2,\cdots,N})$ has stabilized. The overall SAGE algorithm to
estimate the channel parameters in $\bm{{\bf\eta}}$ is shown in Algorithm 2,
where the channel parameters in $\bm{\eta}$ estimated in Section IV are taken
as the initial values.

Algorithm 2 The SAGE algorithm to estimate parameters $\bm{{\bf\eta}}$

0: (1) The received signal $\\{{\bf
y}_{t}[n]\\}_{t=1,\cdots,T}^{n=1,\cdots,N}$; (2) Initial values of all channel
parameters in
${\hat{\bm{\eta}}}^{(0)}=[({\hat{\bm{\eta}}}_{0}^{(0)})^{T},({\hat{\bm{\eta}}}_{1}^{(0)})^{T},\cdots,({\hat{\bm{\eta}}}_{q}^{(0)})^{T},\cdots,({\hat{\bm{\eta}}}_{Q}^{(0)})^{T}]^{T}$
, where
${\hat{\bm{\eta}}}_{q}^{(0)}=[\tau_{q}^{(0)},\delta_{q,\textrm{R}}^{(0)},\delta_{q,\textrm{I}}^{(0)},\theta_{{\rm
t},q}^{(0)},\phi_{{\rm in},q}^{(0)},\psi_{{\rm in},q}^{(0)}]^{T}$; (3) The
error thresholds $\bm{{\bf\varepsilon}}_{1}\in\mathbb{C}^{6(Q+1)\times 1}$ and
$\varepsilon_{2}$; (4) $k=0$;

1: while $1$ do

2: $q=k{\ }{\rm mod}(Q+1)$;

3: Estimate $\hat{\bf{y}}_{q,t}^{(k+1)}[n],t=1,\cdots,T,n=1,\cdots,N$, from
$\\{{\bf y}_{t}[n]\\}_{t=1,\cdots,T}^{n=1,\cdots,N}$ with
$\hat{\bm{\eta}}^{(k)}$ according to (59);

4: Estimate the channel parameters in
${\hat{\bm{\eta}}}_{q}^{(k+1)}=[\tau_{q}^{(k+1)},\delta_{q,\textrm{R}}^{(k+1)},\delta_{q,\textrm{I}}^{(k+1)},\theta_{{\rm
t},q}^{(k+1)},\phi_{{\rm in},q}^{(k+1)},\psi_{{\rm in},q}^{(k+1)}]^{T}$ based
on the estimated $\hat{\bf{y}}_{q,t}^{(k)}[n]$ sequentially according to (65);

5: Update
$\hat{\bm{\eta}}^{(k+1)}\\!\\!=\\![(\hat{\bm{\eta}}_{0}^{(k)}\\!)^{T}\\!,\\!\cdots\\!,\\!(\hat{\bm{\eta}}_{q-1}^{(k)}\\!)^{T}\\!,\\!(\hat{\bm{\eta}}_{q}^{(k+1)}\\!)^{T}\\!,\\!(\hat{\bm{\eta}}_{q+1}^{(k)}\\!)^{T},\cdots,(\hat{\bm{\eta}}_{Q}^{(k)})^{T}]^{T}$;

6: if
$\left|{\hat{\bm{\eta}}}^{(k+1)}\\!-\\!{\hat{\bm{\eta}}}^{(k)}\right|\\!\preceq\\!\bm{{\bf\varepsilon}}_{1}$
or $\left|\Lambda({\hat{\bm{\eta}}}^{(k+1)};\\{{\bf
Y}[n]\\}_{n=1,2,\cdots,N})\\!-\\!\Lambda({\hat{\bm{\eta}}}^{(k)};\\{{\bf
Y}[n]\\}_{n=1,2,\cdots,N})\right|\\!\leqslant\\!\varepsilon_{2}$ then

7: break;

8: else

9: ${k}\leftarrow{k}+1$;

10: end if

11: end while

11: Estimated channel parameters
$\hat{\bm{\eta}}=[\hat{\bm{\eta}}_{0}^{T},\hat{\bm{\eta}}_{1}^{T},\cdots,\hat{\bm{\eta}}_{q}^{T},\cdots,\hat{\bm{\eta}}_{Q}^{T}]^{T}$.

##  VI Estimation of Position-Related Parameters

Based on the refined estimates of channel parameters in $\hat{\bm{\eta}}$, the
final parameters in $\tilde{\bm{{\eta}}}$ are estimated, from which the
coordinates ${\bf m}$ of the MS and rotation angle $\alpha$ can be recovered.
The estimation for $\tilde{\bm{{\eta}}}$ also consists of a coarse estimation
process and a refining process, where the former provides the initial values
for the latter.

###  VI-A Coarse Estimation

Regardless of the estimation errors, the final parameters in
$\tilde{\bm{{\eta}}}$ can be calculated directly from the intermediate
parameters $\bm{{\bf\eta}}$ based on their relations. According to the
geometric relationship described in Fig. 1, we can find $\hat{{\bf m}}$ and
$\hat{\alpha}$ in closed-form expressions by using the estimated channel
parameters in $\hat{\bm{{\eta}}}$ as

$\displaystyle{\hat{{\bf m}}}$ $\displaystyle={\bf
r}+(\hat{\tau}_{0}c-\left\|{\bf r}-{\bf b}\right\|_{2})[-\sin\hat{\phi}_{{\rm
in},0}\cos\hat{\psi}_{{\rm in},0},-\sin\hat{\phi}_{{\rm
in},0}\sin\hat{\psi}_{{\rm in},0},-\cos\hat{\phi}_{{\rm in},0}]^{T},$ (66)
$\displaystyle{\hat{\alpha}}$ $\displaystyle=(2\pi-\hat{\psi}_{{\rm
in},0})-\arccos[(\left\|{\bf m}-{\bf r}\right\|_{2}\sin\hat{\theta}_{{\rm
t},0})/(\left\|{\bf m}-{\bf r}\right\|_{2}\sin\hat{\phi}_{{\rm in},0})]$
$\displaystyle=(2\pi-\hat{\psi}_{{\rm in},0})-\arccos(\sin\hat{\theta}_{{\rm
t},0}/\sin\hat{\phi}_{{\rm in},0}).$ (67)

To find the closed-form expression for ${\bf s}^{q}$, we can build the
following set of equations:

$\displaystyle{\hat{{\bf s}}}^{q}$ $\displaystyle={\bf
r}+{\hat{d}}_{sr}[-\sin\hat{\phi}_{{\rm in},q}\cos\hat{\psi}_{{\rm
in},q},-\sin\hat{\phi}_{{\rm in},q}\sin\hat{\psi}_{{\rm
in},q},-\cos\hat{\phi}_{{\rm in},q}]^{T},$ (68a)
$\displaystyle(\hat{\tau}_{q}c-\left\|{\bf r}-{\bf b}\right\|_{2})$
$\displaystyle={\hat{d}}_{sr}+{\hat{d}}_{ms},$ (68b)
$\displaystyle\left\|{\hat{{\bf s}}}^{q}-{\hat{{\bf
m}}}\right\|_{2}\sin{\hat{\theta}}_{{\rm t},q}$
$\displaystyle=({\hat{s}}_{x}^{q}-{\hat{m}}_{x})\cos{\hat{\alpha}}-({\hat{s}}_{y}^{q}-{\hat{m}}_{y})\sin{\hat{\alpha}}.$
(68c)

The set of equations in (68) can be readily solved. After some mathematical
derivations, the closed-form solution for the scatterer coordinates in ${\bf
s}^{q}$ can be represented as

$\displaystyle\hat{s}_{x}^{q}$
$\displaystyle=\frac{(AD+r_{x})\sin\hat{\theta}_{{\rm
t},q}+\hat{m}_{x}A\cos\hat{\alpha}+(r_{y}A-r_{x}B-\hat{m}_{y}A)\sin\hat{\alpha}}{\sin\hat{\theta}_{{\rm
t},q}+A\cos\hat{\alpha}-B\sin\hat{\alpha}},$ (69a)
$\displaystyle\hat{s}_{y}^{q}$
$\displaystyle=r_{y}+(\hat{s}_{x}^{q}-r_{x})\frac{B}{A},$ (69b)
$\displaystyle\hat{s}_{z}^{q}$
$\displaystyle=r_{z}+(\hat{s}_{x}^{q}-r_{x})\frac{C}{A},$ (69c)

where $A=-\sin{\hat{\phi}_{{\rm{in}},q}}\cos{\hat{\psi}_{{\rm{in}},q}}$,
$B=-\sin{\hat{\phi}_{{\rm{in}},q}}\sin{\hat{\psi}_{{\rm{in}},q}}$,
$C=-\cos{\hat{\phi}_{{\rm{in}},q}}$ and
$D={\hat{\tau}_{q}}c-{\left\|{{\bf{r}}-{\bf{b}}}\right\|_{2}}$.

Although the closed-form solutions for the coordinates of the MS, rotation
angle of the MS, coordinates of the scatterer can be obtained from (66), (67)
and (69), these estimations are not accurate due to the estimation error of
the channel parameters of $\bm{{\bf\eta}}$. Moreover, the closed-form
estimations for the final parameters only exploit some of the channel
parameters, and do not extract the full information in all the channel
parameters. Thus, refining estimation for the final parameters in
$\tilde{\bm{{\eta}}}$ should be further performed by taking the closed-form
solutions as the initial values.

###  VI-B Refining Estimation

The equations in Section III-B describe the relationship between $\bm{{\eta}}$
and $\tilde{\bm{{\eta}}}$. Since the coordinates of the RIS and the BS are
assumed to be known, the angles of the RIS-BS link are taken as known values,
which means that the equations in (14) are not estimated but calculated in
advance. The expressions in (12), (13) and (15) describe a mapping
$\bm{{\eta}}=\mathcal{G}(\tilde{\bm{{\eta}}})$. Thus, the estimated
$\hat{\tilde{\bm{{\eta}}}}$ can be obtained by

$\displaystyle\hat{\tilde{\bm{{\eta}}}}$
$\displaystyle=\arg\underset{\tilde{\bm{{\eta}}}}{\min}{\
}[\hat{\bm{{\eta}}}-\mathcal{G}(\tilde{\bm{{\eta}}})]^{T}\mathbf{J}_{\hat{\bm{{\eta}}}}[\hat{\bm{{\eta}}}-\mathcal{G}(\tilde{\bm{{\eta}}})],$
(70)

where $\mathbf{J}_{\hat{\bm{{\eta}}}}$ is the fisher information matrix (FIM)
for the channel parameters $\bm{{\bf\eta}}$ calculated by substituting the
estimated value $\hat{\bm{{\eta}}}$. The FIM for $\bm{{\bf\eta}}$ is given in
(74) of Section VII-A. The estimator in (70) is asymptotically (w.r.t.
$T\times N$) equivalent to the ML estimation of the transformed parameter
$\tilde{\bm{{\eta}}}$[30][31]. The problem in (70) is a nonlinear least
squares problem, which can be solved by the Levenberg-Marquardt Method [32,
23], which can be initialized in Section VI-A.

##  VII Position and Orientation Estimation Fundamental Bounds

To provide a benchmark for the accuracy of the proposed estimation algorithms,
we derive the error bounds for the formulated estimation problems.

###  VII-A CRB for Channel Parameters

We first derive the FIM ${\bf J}_{\bm{\eta}}$ of the channel parameters in
$\bm{{\bf\eta}}$. The joint probability density function of the vectorized
signal ${\bf y}[n]$ is given in (50). By taking the natural logarithm of (50)
and omitting the constant terms,the log-likelihood functon
$\Lambda(\bm{{\bf\eta}};\\{{\bf y}[n]\\}_{n=1,\cdots,N})$ of $\bm{{\bf\eta}}$
given the observed signal $\\{{\bf y}[n]\\}_{n=1,\cdots,N}$ can be written as

$\displaystyle\Lambda(\bm{{\bf\eta}};\\{{\bf y}[n]\\}_{n=1,\cdots,N})$
$\displaystyle\propto-\sum\limits_{n=1}^{N}\left\\{({\bf
y}[n]-\bm{\mu}[n])^{H}{\bf\bm{{\Xi}}}^{-1}({\bf y}[n]-\bm{\mu}[n])\right\\},$
(71)

where $\bm{\mu}[n]=\tilde{\bf{H}}[n]{\bf x}[n]$ and
$\tilde{\bf{H}}[n]\triangleq\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}(\bm{{\Sigma}}_{q}^{T}\otimes\tilde{\bf{H}}_{q})$.

Defining $\hat{\bm{\eta}}$ as the unbiased estimator of $\bm{\eta}$, the mean
squared error (MSE) of $\hat{\bm{\eta}}$ is bounded as

$\displaystyle\mathbb{E}\left[(\hat{\bm{\eta}}-\bm{\eta})(\hat{\bm{\eta}}-\bm{\eta})^{T}\right]$
$\displaystyle\succeq{\bf J}_{\bm{\eta}}^{-1},$ (72)

where ${\bf J}_{\bm{\eta}}\in\mathbb{C}^{6(Q+1)\times 6(Q+1)}$ is the FIM of
$\bm{\eta}$. The element at the $u$-th row and the $v$-th column of ${\bf
J}_{\bm{\eta}}$ can be represented as [33]

$\displaystyle[{\bf J}_{\bm{\eta}}]_{u,v}$
$\displaystyle=-\mathbb{E}\left\\{\frac{\partial^{2}\Lambda(\bm{{\bf\eta}};\\{{\bf
y}[n]\\}_{n=1,\cdots,N})}{\partial[\bm{\eta}]_{u}\partial[\bm{\eta}]_{v}}\right\\}$
$\displaystyle=\mathbb{E}\left\\{\frac{\partial\Lambda(\bm{{\bf\eta}};\\{{\bf
y}[n]\\}_{n=1,\cdots,N})}{\partial[\bm{\eta}]_{u}}\frac{\partial\Lambda(\bm{{\bf\eta}};\\{{\bf
y}[n]\\}_{n=1,\cdots,N})}{\partial[\bm{\eta}]_{v}}\right\\}.$ (73)

It is noted that the $(u,u)$-th element of ${\bf J}_{\bm{\eta}}^{-1}$ (i.e.,
$[{\bf J}_{\bm{\eta}}^{-1}]_{u,u}$) is the CRB for the $u$-th parameter of
${\bf\bm{\eta}}$.

By substituting (71) into (73), (73) can be rewritten as

$\displaystyle[{\bf J}_{\bm{\eta}}]_{u,v}$
$\displaystyle\overset{(d)}{=}\sum\limits_{n=1}^{N}\left\\{2\mathfrak{Re}[\frac{\partial\bm{\mu}^{H}[n]}{\partial[\bm{\eta}]_{u}}\bm{{\Xi}}^{-1}\frac{\partial\bm{\mu}[n]}{\partial[\bm{\eta}]_{v}}]+{\rm
Tr}[\bm{{\Xi}}^{-1}\frac{\partial{\bf\bm{{\Xi}}}}{\partial[\bm{\eta}]_{u}}\bm{{\Xi}}^{-1}\frac{\partial{\bf\bm{{\Xi}}}}{\partial[\bm{\eta}]_{v}}]\right\\}$
$\displaystyle\overset{(e)}{=}\frac{2}{\sigma^{2}}\sum\limits_{n=1}^{N}\mathfrak{Re}\left\\{{\bf
x}^{H}[n]\frac{\partial\tilde{\bf{H}}^{H}[n]}{\partial[\bm{\eta}]_{u}}\frac{\partial\tilde{\bf{H}}[n]}{\partial[\bm{\eta}]_{v}}{\bf
x}[n]\right\\},$ (74)

where the step $(d)$ can refer to Appendix 15C of [33]. The step $(e)$ is
obtained because ${\bf\bm{{\Xi}}}$ is a constant matrix. From (74), we find
that the key to obtain the FIM ${\bf J}_{\bm{\eta}}$ is to calculate
$\frac{\partial\tilde{\bf{H}}[n]}{\partial[\bm{\eta}]_{u}},u=1,2,\cdots,6(Q+1)$,
which are derived in Appendix C-A.

###  VII-B Position and Orientation Error Bounds

We further derive the FIM ${\bf J}_{\tilde{\bm{{\eta}}}}$ of the final
position parameters in $\tilde{\bm{{\eta}}}$. Through a transformation of
variables from $\bm{{\bf\eta}}$ to $\tilde{\bm{{\eta}}}$, the ${\bf
J}_{\tilde{\bm{{\eta}}}}\in\mathbb{C}^{[5(Q+1)+1]\times[5(Q+1)+1]}$ can be
determined from ${\bf J}_{\bm{\eta}}\in\mathbb{C}^{[6(Q+1)]\times[6(Q+1)]}$ as

$\displaystyle{\bf J}_{\tilde{\bm{{\eta}}}}$ $\displaystyle=\mathbf{T}{\bf
J}_{\bm{\eta}}\mathbf{T}^{T},$ (75)

where
$\mathbf{T}\triangleq\frac{\partial\bm{\eta}^{T}}{\partial\tilde{\bm{{\eta}}}}$
is the transformation matrix. To obtain
$\mathbf{T}\in\mathbb{C}^{[5(Q+1)+1]\times[6(Q+1)]}$, the relationships
between $\bm{{\bf\eta}}$ and $\tilde{\bm{{\eta}}}$ described in Section III-B
are utilized. The specific derivations for the elements of $\mathbf{T}$ can be
found in Appendix C-B.

The PEB is obtained by inverting ${\bf J}_{\tilde{\bm{{\eta}}}}$, adding the
diagonal entries of the $3\times 3$ sub-matrix, and taking the root square as

$\displaystyle\mathrm{PEB}$ $\displaystyle=\sqrt{\mathrm{tr}\left\\{[{\bf
J}_{\tilde{\bm{{\eta}}}}^{-1}]_{[2(Q+1)+1]:[2(Q+1)+3],[2(Q+1)+1]:[2(Q+1)+3]}\right\\}},$
(76)

and the OEB is obtained as

$\displaystyle\mathrm{OEB}$ $\displaystyle=\sqrt{{[{\bf
J}_{\tilde{\bm{{\eta}}}}^{-1}]_{[2(Q+1)+4],[2(Q+1)+4]}}}.$ (77)

##  VIII Simulation Results

In this section, we present simulation results to evaluate the performance of
the proposed estimation algorithms.

###  VIII-A Simulation Setup

We set $f_{c}=4.9$ GHz, $B=20$ MHz, $c=3\times 10^{8}$ m/s and $N=20$. The
noise power density is $-174$ dBm/Hz. The numbers of transmit and receive
antennas are set to $N_{\mathrm{m}}=16$ and $N_{\mathrm{b}}=40$, respectively.
Their array element spacing is
$d_{\textrm{{\rm{BS}}}}=d_{\textrm{{\rm{MS}}}}=\lambda/2$. The number of
reflecting elements of the RIS is
$N_{\mathrm{r}}=N_{{\rm{a}}}N_{{\rm{e}}}=100$, where
$N_{{\rm{a}}}=N_{{\rm{e}}}=10$. The reflection element spacing of the RIS is
$d_{\textrm{{\rm{RIS,a}}}}=d_{\textrm{{\rm{RIS,e}}}}=\lambda/3$. The number of
total time slots utilized for positioning is $T=37$, in which the number of
time slots in the first stage is $T_{1}=16$ and the remaining number of time
slots is $T_{2}=21$. In the remaining $T_{2}$ time slots, there are
$\Upsilon=7$ time blocks, and each time block consists of $V=3$ time slots.
The BS is located at $\mathbf{b}=[0,0,28]^{T}$ m, the RIS is located at
$\mathbf{r}=[-6,8,20]^{T}$ m, and the MS is located at
$\mathbf{m}=[22,35,1.5]^{T}$ m with the rotation angle of $\alpha=75^{\circ}$.
We assume that there is a scatterer in the link from the MS to the RIS, and
the location of the scatterer is $\mathbf{s}=[6,5,3]^{T}$ m. The relative
positions of the BS, RIS, MS and scatterer are described in Fig. 4.

Figure 4: Relative position of different nodes.

The complex channel gain $\delta_{0}=\delta_{{\rm RB},0}\delta_{{\rm MR},0}$
of the VLoS path follows the distribution of
$\mathcal{CN}(0,10^{\frac{-{\rm{PL}}_{0}}{10}})$, where ${\rm{PL}}_{0}$ is the
path loss of the MS-RIS-BS link. Considering the urban UMa scenario[34], the
large-scale fading of the MS-RIS-BS reflection link can be modeled as
${\rm{PL}}_{0}=28+40\mathrm{log10}(f_{c})+10\alpha_{\mathrm{r}}\mathrm{log10}(d_{\mathrm{MR}}d_{\mathrm{RB}})+\xi$[35,
36], where $\alpha_{\mathrm{r}}$ is the path loss exponent of the reflection
link, $d_{\mathrm{MR}}$ is the three dimensional distance (in meters) of the
MS-RIS link, $d_{\mathrm{RB}}$ is the three dimensional distance (in meters)
of the RIS-BS link, and $\xi$ is the lognormal shadowing. $\xi$ follows the
distribution of $\mathcal{N}(0,\sigma_{\mathrm{SF}}^{2})$, where
$\sigma_{\mathrm{SF}}^{2}$ is the shadow fading std (standard deviation) in dB
and $\sigma_{\mathrm{SF}}=4$. The $\alpha_{\mathrm{r}}$ for the VLoS is set to
$\alpha_{\mathrm{r}}=2.2$. The complex channel gain $\delta_{q}=\delta_{{\rm
RB},0}\delta_{{\rm MR},q}$ of the NLoS path follows the distribution of
$\mathcal{CN}(0,10^{\frac{-{\rm{PL}}_{q}}{10}})$, where ${\rm{PL}}_{q}$ is the
path loss of the MS-Scatterer-$q$-RIS-BS link. Here, we set
${\rm{PL}}_{q}={\rm{PL}}_{0}+3\textrm{ dB}$. The numbers of the per-
discretized grids are set as $G_{\mathrm{m}}=128$, $G_{\mathrm{a}}=10$ and
$G_{\mathrm{e}}=10$. Unless otherwise specified, the RMSE of the estimation
algorithm is calculated over 1000 Monte Carlo simulations, where the channel
gain and noise are random.

###  VIII-B Comparison between the Coarse and Refined Estimation of the
Channel Parameters

(a) Real part of $\delta$

(b) Imaginary part of $\delta$

(c) TOA $\tau$

(d) AOD $\theta_{\rm t}$

(e) Elevation AOA $\phi_{\rm in}$ at the RIS

(f) Azimuth AOA $\psi_{\rm in}$ at the RIS

Figure 5: The RMSE of channel parameters estimated by different algorithms
versus transmit power.

The channel parameters in $\bm{\eta}$ include the complex channel gain, the
AOD at the MS, the azimuth AOA at the RIS, the elevation AOA at the RIS, and
the TOA. There parameters are first coarsely estimated, then refined by the
SAGE algorithm. To evaluate the estimation performance for channel parameters,
the RMSEs of different estimation algorithms versus transmit power are shown
in Fig. 5, where CRLB of different channel parameters are given as benchmarks.
Fig. 5(a) and Fig. 5(b) show the estimation performance for the real part and
imaginary part of the complex channel gain, respectively. Fig. 5(c) evaluates
the estimation performance for the TOA, and Fig. 5(d) evaluates the estimation
performance for the AOD. The estimation performance of the azimuth AOA and
elevation AOD at the RIS is respectively shown in Fig. 5(e) and 5(f).

From these figures, we can have the following observations and conclusions.
Firstly, it is observed that in all sub-figures of Fig. 5, the RMSE of the
channel parameters for both the coarse estimation and refining estimation
decreases with the transmit power, which means that the estimation accuracy
can be improved with transmit power. Secondly, in Fig. 5(d), by comparing the
black line marked by $\star$ (or $\rhd$) with the blue line marked by $\circ$
(or $\times$), we find that the RMSEs of the refined AODs by MLE are lower
than those of the coarsely-estimated AODs by the DCS-SOMP algorithm, and the
performance gap becomes significant at high transmit power. This verifies the
effectiveness of refining AODs by MLE, which mitigates the error propagation
into AOA/TOA estimation, and generate better initial values for the SAGE
algorithm. Thirdly, by comparing the blue line marked by $\circ$ (or $\times$)
with the green line marked by $\triangledown$ (or $\lozenge$) in all the
subfigures of Fig. 5, it can be seen that the RMSEs of channel parameters by
coarse estimation are greatly reduced after the refining estimation by the
SAGE algorithm, and the performance gap increases with the transmit power.
This verifies the effectiveness of the refining estimation by the SAGE
algorithm. Fourthly, by comparing the green line marked by $\triangledown$ (or
$\lozenge$) with the red line marked by $\square$ (or $\triangle$) of all
subfigures in Fig. 5, it is observed that after the refining estimation
concatenated by the coarse estimation, the RMSEs of all channel parameters by
the SAGE algorithm gradually approach their corresponding CRLBs as the
transmit power increases, which validates the good estimation performance of
our proposed estimation algorithms.

(a) Position error (b) Orientation error

Figure 6: The RMSE of position and rotation angle estimation versus transmit
power.

###  VIII-C Positioning Performance Evaluation

From the estimated channel parameters, we can further estimate the coordinates
and angle rotation of the MS. Thus, the RMSE of the position and rotation
estimation can be evaluated. Fig. 6 shows the RMSEs of coordinates and the
rotation angle of the MS versus the transmit power. It is observed that as the
transmit power increases, RMSE($\hat{\mathbf{m}}$) and RMSE($\hat{\alpha}$)
gradually decrease, and finally approach the PEB and OEB respectively. The
RMSEs of $\hat{\mathbf{m}}$ and $\hat{\alpha}$ become close to the PEB and OEB
at about $P_{{\rm m}}=10$ dBm, which is consistent with the trend of the RMSE
curves for channel parameters. The curves in Fig. 6 demonstrate that the RMSE
of the proposed positioning algorithm can approach the lower bound of the
positioning performance, thus verifies the superiority of our proposed
algorithms.

##  IX Summary and Conclusion

We proposed algorithms for estimating the position coordinates and orientation
angle of an MS in an RIS-aided mmWave MIMO system over dispersive channels.
The two-step positioning scheme is utilized, where the channel parameters are
accurately estimated first, then the position-related parameters are
estimated. For the channel parameter estimation, the refining estimation is
cascaded with the coarse estimation. A series of algorithms are proposed for
coarse estimation based on the sparsity of the mm-wave channel. The refining
estimation is realized by the SAGE algorithm. Finally, the position and
orientation parameters are estimated from the estimated channel parameters.
The PEB and OEB are derived to provide performance benchmark. Simulation
results demonstrate the superior performance of our proposed algorithms for
both the channel parameter estimation and the position estimation.

##  Appendix A Simplification of the ML Function for Refining AODs

Since we assume that the pilot signal is the same over all subcarriers,
${\tilde{{\bf X}}}_{t}[n]$ can be written as ${\tilde{{\bf
X}}}_{t}=[{\tilde{{\bf H}}}_{0}{\bf x}_{t},{\tilde{{\bf H}}}_{1}{\bf
x}_{t},\cdots,{\tilde{{\bf H}}}_{Q}{\bf x}_{t}]$. Then (26) could be rewritten
as

$\displaystyle\hat{\bm{\theta}}_{{\rm t}}$
$\displaystyle=\arg\mathop{\min}\limits_{{\bm{\theta}}_{{\rm
t}}}\sum\limits_{n=1}^{N}\sum\limits_{t=1}^{T_{1}}\left\|{\bf
y}_{t}[n]-{\tilde{{\bf X}}}_{t}[\sum\limits_{t=1}^{T_{1}}({\tilde{{\bf
X}}}_{t}^{H}{\tilde{{\bf X}}}_{t})]^{-1}\sum\limits_{t=1}^{T_{1}}({\tilde{{\bf
X}}}_{t}^{H}{\bf y}_{t}[n])\right\|_{2}^{2}.$ (78)

The term $\sum\limits_{t=1}^{T_{1}}({\tilde{{\bf X}}}_{t}^{H}{\bf y}_{t}[n])$
in (78) can be simplified as follows:

$\displaystyle\sum\limits_{t=1}^{T_{1}}({\tilde{{\bf X}}}_{t}^{H}{\bf
y}_{t}[n])$
$\displaystyle=\sum\limits_{t=1}^{T_{1}}\left[\begin{array}[]{c}{\bf
x}_{t}^{H}{\tilde{{\bf H}}}_{0}^{H}{\bf y}_{t}[n]\\\ {\bf
x}_{t}^{H}{\tilde{{\bf H}}}_{1}^{H}{\bf y}_{t}[n]\\\ \vdots\\\ {\bf
x}_{t}^{H}{\tilde{{\bf H}}}_{Q}^{H}{\bf y}_{t}[n]\end{array}\right]$ (83)
$\displaystyle=\sum\limits_{t=1}^{T_{1}}[{\bf y}_{t}^{H}[n]{\tilde{{\bf
H}}}_{0}{\bf x}_{t},{\bf y}_{t}^{H}[n]{\tilde{{\bf H}}}_{1}{\bf
x}_{t},\cdots,{\bf y}_{t}^{H}[n]{\tilde{{\bf H}}}_{Q}{\bf x}_{t}]^{H}$
$\displaystyle=\sum\limits_{t=1}^{T_{1}}[{\bf y}_{t}^{H}[n]{\bf a}_{{\rm
B,0}}{\bf a}_{{\rm M,0}}^{H}{\bf x}_{t},{\bf y}_{t}^{H}[n]{\bf a}_{{\rm
B,0}}{\bf a}_{{\rm M,1}}^{H}{\bf x}_{t},\cdots,{\bf y}_{t}^{H}[n]{\bf a}_{{\rm
B,0}}{\bf a}_{{\rm M,}Q}^{H}{\bf x}_{t}]^{H}$ $\displaystyle=[{\bf a}_{{\rm
M,0}}^{H}(\sum\limits_{t=1}^{T_{1}}{\bf x}_{t}{\bf y}_{t}^{H}[n]){\bf a}_{{\rm
B,0}},{\bf a}_{{\rm M,1}}^{H}(\sum\limits_{t=1}^{T_{1}}{\bf x}_{t}{\bf
y}_{t}^{H}[n]){\bf a}_{{\rm B,0}},\cdots,{\bf a}_{{\rm
M,}Q}^{H}(\sum\limits_{t=1}^{T_{1}}{\bf x}_{t}{\bf y}_{t}^{H}[n]){\bf a}_{{\rm
B,0}}]^{H}$ $\displaystyle=\left[\begin{array}[]{c}{\bf a}_{{\rm
B,0}}^{H}(\sum\limits_{t=1}^{T_{1}}{\bf y}_{t}[n]{\bf x}_{t}^{H}){\bf a}_{{\rm
M,0}}\\\ {\bf a}_{{\rm B,0}}^{H}(\sum\limits_{t=1}^{T_{1}}{\bf y}_{t}[n]{\bf
x}_{t}^{H}){\bf a}_{{\rm M,1}}\\\ \vdots\\\ {\bf a}_{{\rm
B,0}}^{H}(\sum\limits_{t=1}^{T_{1}}{\bf y}_{t}[n]{\bf x}_{t}^{H}){\bf a}_{{\rm
M,}Q}\end{array}\right]$ (88) $\displaystyle=({\bf a}_{{\rm
B,0}}^{H}(\sum\limits_{t=1}^{T_{1}}{\bf y}_{t}[n]{\bf x}_{t}^{H}){\bf A}_{{\rm
M,}Q})^{T}$ $\displaystyle={\bf A}_{{\rm
M,}Q}^{T}(\sum\limits_{t=1}^{T_{1}}{\bf x}_{t}^{*}{\bf y}_{t}^{T}[n]){\bf
a}_{{\rm B,0}}^{*},$ (89)

where ${\tilde{{\bf H}}}_{q}{\\!}\triangleq{\\!}{\bf a}_{{\rm
B}}(\vartheta_{{\rm r},0}){\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm t},q})$, ${\bf
a}_{{\rm B,0}}{\\!}\triangleq{\\!}{\bf a}_{{\rm B}}(\vartheta_{{\rm r},0})$,
${\bf a}_{{\rm M},q}{\\!}\triangleq{\\!}{\bf a}_{{\rm M}}(\vartheta_{{\rm
t},q})$ and ${\bf A}_{{\rm M,}Q}{\\!}\triangleq{\\!}[{\bf a}_{{\rm M,0}},{\bf
a}_{{\rm M,1}},\cdots,{\bf a}_{{\rm M,}Q}]$.

The term ${\tilde{{\bf X}}}_{t}[\sum\limits_{t=1}^{T_{1}}({\tilde{{\bf
X}}}_{t}^{H}{\tilde{{\bf X}}}_{t})]^{-1}$ in (78) can be simplified as

$\displaystyle{\tilde{{\bf X}}}_{t}^{H}{\tilde{{\bf X}}}_{t}$
$\displaystyle=\left[\begin{array}[]{c}{\bf x}_{t}^{H}{\bf a}_{{\rm M,0}}\\\
{\bf x}_{t}^{H}{\bf a}_{{\rm M,1}}\\\ \vdots\\\ {\bf x}_{t}^{H}{\bf a}_{{\rm
M,}Q}\end{array}\right]{\bf a}_{{\rm B,0}}^{H}{\bf a}_{{\rm B,0}}[{\bf
a}_{{\rm M,0}}^{H}{\bf x}_{t},{\bf a}_{{\rm M,1}}^{H}{\bf x}_{t},\cdots,{\bf
a}_{{\rm M,}Q}^{H}{\bf x}_{t}]$ (94)
$\displaystyle{\rm=}N_{\mathrm{b}}\left[\begin{array}[]{cccccc}{\bf
x}_{t}^{H}{\bf a}_{{\rm M,0}}{\bf a}_{{\rm M,0}}^{H}{\bf x}_{t}&{\bf
x}_{t}^{H}{\bf a}_{{\rm M,0}}{\bf a}_{{\rm M,1}}^{H}{\bf x}_{t}&\cdots&{\bf
x}_{t}^{H}{\bf a}_{{\rm M,0}}{\bf a}_{{\rm M,}j}^{H}{\bf x}_{t}&\cdots&{\bf
x}_{t}^{H}{\bf a}_{{\rm M,0}}{\bf a}_{{\rm M,}Q}^{H}{\bf x}_{t}\\\ {\bf
x}_{t}^{H}{\bf a}_{{\rm M,1}}{\bf a}_{{\rm M,0}}^{H}{\bf x}_{t}&{\bf
x}_{t}^{H}{\bf a}_{{\rm M,1}}{\bf a}_{{\rm M,1}}^{H}{\bf x}_{t}&\cdots&{\bf
x}_{t}^{H}{\bf a}_{{\rm M,1}}{\bf a}_{{\rm M,}j}^{H}{\bf x}_{t}&\cdots&{\bf
x}_{t}^{H}{\bf a}_{{\rm M,1}}{\bf a}_{{\rm M,}Q}^{H}{\bf x}_{t}\\\
\vdots&\vdots&\ddots&\vdots&\ddots&\vdots\\\ {\bf x}_{t}^{H}{\bf a}_{{\rm
M,}i}{\bf a}_{{\rm M,0}}^{H}{\bf x}_{t}&{\bf x}_{t}^{H}{\bf a}_{{\rm M,}i}{\bf
a}_{{\rm M,1}}^{H}{\bf x}_{t}&\cdots&{\bf x}_{t}^{H}{\bf a}_{{\rm M,}i}{\bf
a}_{{\rm M,}j}^{H}{\bf x}_{t}&\cdots&{\bf x}_{t}^{H}{\bf a}_{{\rm M,}i}{\bf
a}_{{\rm M,}Q}^{H}{\bf x}_{t}\\\ \vdots&\vdots&\ddots&\vdots&\ddots&\vdots\\\
{\bf x}_{t}^{H}{\bf a}_{{\rm M,}Q}{\bf a}_{{\rm M,0}}^{H}{\bf x}_{t}&{\bf
x}_{t}^{H}{\bf a}_{{\rm M,}Q}{\bf a}_{{\rm M,1}}^{H}{\bf x}_{t}&\cdots&{\bf
x}_{t}^{H}{\bf a}_{{\rm M,}Q}{\bf a}_{{\rm M,}j}^{H}{\bf x}_{t}&\cdots&{\bf
x}_{t}^{H}{\bf a}_{{\rm M,}Q}{\bf a}_{{\rm M,}Q}^{H}{\bf
x}_{t}\end{array}\right]$ (101)
$\displaystyle\mathop{{\rm=}}\limits^{(g)}N_{\mathrm{b}}\left[\begin{array}[]{cccccc}{\bf
a}_{{\rm M,0}}^{H}{\bf x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,0}}&{\bf a}_{{\rm
M,1}}^{H}{\bf x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,0}}&\cdots&{\bf a}_{{\rm
M,}j}^{H}{\bf x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,0}}&\cdots&{\bf a}_{{\rm
M,}Q}^{H}{\bf x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,0}}\\\ {\bf a}_{{\rm
M,0}}^{H}{\bf x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,1}}&{\bf a}_{{\rm
M,1}}^{H}{\bf x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,1}}&\cdots&{\bf a}_{{\rm
M,}j}^{H}{\bf x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,1}}&\cdots&{\bf a}_{{\rm
M,}Q}^{H}{\bf x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,1}}\\\
\vdots&\vdots&\ddots&\vdots&\ddots&\vdots\\\ {\bf a}_{{\rm M,0}}^{H}{\bf
x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,}i}&{\bf a}_{{\rm M,1}}^{H}{\bf
x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,}i}&\cdots&{\bf a}_{{\rm M,}j}^{H}{\bf
x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,}i}&\cdots&{\bf a}_{{\rm M,}Q}^{H}{\bf
x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,}i}\\\
\vdots&\vdots&\ddots&\vdots&\ddots&\vdots\\\ {\bf a}_{{\rm M,0}}^{H}{\bf
x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,}Q}&{\bf a}_{{\rm M,1}}^{H}{\bf
x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,}Q}&\cdots&{\bf a}_{{\rm M,}j}^{H}{\bf
x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,}Q}&\cdots&{\bf a}_{{\rm M,}Q}^{H}{\bf
x}_{t}{\bf x}_{t}^{H}{\bf a}_{{\rm M,}Q}\end{array}\right]$ (108)
$\displaystyle=N_{\mathrm{b}}({\bf A}_{{\rm M,}Q}^{H}{\bf x}_{t}{\bf
x}_{t}^{H}{\bf A}_{{\rm M,}Q})^{T}$ $\displaystyle=N_{\mathrm{b}}{\bf A}_{{\rm
M,}Q}^{T}{\bf x}_{t}^{*}{\bf x}_{t}^{T}{\bf A}_{{\rm M,}Q}^{*},$ (109)

where the step $(g)$ is obtained by
${\bf{x}}_{t}^{H}{\bf{a}}_{{\rm{M,}}i}{\bf{a}}_{{\rm{M,}}j}^{H}{{\bf{x}}_{t}}={\rm{Tr}}({\bf{a}}_{{\rm{M,}}j}^{H}{{\bf{x}}_{t}}{\bf{x}}_{t}^{H}{\bf{a}}_{{\rm{M,}}i})={\bf{a}}_{{\rm{M,}}j}^{H}{{\bf{x}}_{t}}{\bf{x}}_{t}^{H}{\bf{a}}_{{\rm{M,}}i}$.
Then ${\tilde{{\bf X}}}_{t}[\sum\limits_{t=1}^{T_{1}}({\tilde{{\bf
X}}}_{t}^{H}{\tilde{{\bf X}}}_{t})]^{-1}$ could be rewritten as

$\displaystyle{\tilde{{\bf X}}}_{t}[\sum\limits_{t=1}^{T_{1}}({\tilde{{\bf
X}}}_{t}^{H}{\tilde{{\bf X}}}_{t})]^{-1}$
$\displaystyle=\frac{1}{N_{\mathrm{b}}}{\bf a}_{{\rm B,0}}[{\bf a}_{{\rm
M,0}}^{H}{\bf x}_{t},{\bf a}_{{\rm M,1}}^{H}{\bf x}_{t},\cdots,{\bf a}_{{\rm
M,}Q}^{H}{\bf x}_{t}][{\bf A}_{{\rm M,}Q}^{T}\sum\limits_{t=1}^{T_{1}}({\bf
x}_{t}^{*}{\bf x}_{t}^{T}){\bf A}_{{\rm M,}Q}^{*}]^{-1}$
$\displaystyle=\frac{1}{N_{\mathrm{b}}}{\bf a}_{{\rm B,0}}{\bf x}_{t}^{T}{\bf
A}_{{\rm M,}Q}^{*}[{\bf A}_{{\rm M,}Q}^{T}\sum\limits_{t=1}^{T_{1}}({\bf
x}_{t}^{*}{\bf x}_{t}^{T}){\bf A}_{{\rm M,}Q}^{*}]^{-1},$ (110)

where $[{\bf a}_{{\rm M,1}}^{H}{\bf x}_{t},{\bf a}_{{\rm M,2}}^{H}{\bf
x}_{t},\cdots,{\bf a}_{{\rm M,}Q}^{H}{\bf x}_{t}]=({\bf A}_{{\rm M,}Q}^{H}{\bf
x}_{t})^{T}={\bf x}_{t}^{T}{\bf A}_{{\rm{M,}}Q}^{*}$.

By using (89) and (110), the term ${\tilde{{\bf
X}}}_{t}[\sum\limits_{t=1}^{T_{1}}({\tilde{{\bf X}}}_{t}^{H}{\tilde{{\bf
X}}}_{t})]^{-1}\sum\limits_{t=1}^{T_{1}}({\tilde{{\bf X}}}_{t}^{H}{\bf
y}_{t}[n])$ in (78) could be simplified as

$\displaystyle{\tilde{{\bf X}}}_{t}[\sum\limits_{t=1}^{T_{1}}({\tilde{{\bf
X}}}_{t}^{H}{\tilde{{\bf X}}}_{t})]^{-1}\sum\limits_{t=1}^{T_{1}}({\tilde{{\bf
X}}}_{t}^{H}{\bf y}_{t}[n])$ $\displaystyle=$
$\displaystyle\frac{1}{N_{\mathrm{b}}}{\bf a}_{{\rm B,0}}{\bf x}_{t}^{T}{\bf
A}_{{\rm M,}Q}^{*}[{\bf A}_{{\rm M,}Q}^{T}\sum\limits_{t=1}^{T_{1}}({\bf
x}_{t}^{*}{\bf x}_{t}^{T}){\bf A}_{{\rm M,}Q}^{*}]^{-1}[{\bf A}_{{\rm
M,}Q}^{T}\sum\limits_{t=1}^{T_{1}}({\bf x}_{t}^{*}{\bf y}_{t}^{T}[n]){\bf
a}_{{\rm B,0}}^{*}]$ $\displaystyle=$
$\displaystyle\frac{1}{N_{\mathrm{b}}}{\bf a}_{{\rm B,0}}{\bf a}_{{\rm
B,0}}^{H}(\sum\limits_{t=1}^{T_{1}}{\bf y}_{t}[n]{\bf x}_{t}^{H}){\bf A}_{{\rm
M,}Q}[{\bf A}_{{\rm M,}Q}^{H}\sum\limits_{t=1}^{T_{1}}({\bf x}_{t}{\bf
x}_{t}^{H}){\bf A}_{{\rm M,}Q}]^{-1}{\bf A}_{{\rm M,}Q}^{H}{\bf x}_{t}$
$\displaystyle=$ $\displaystyle\frac{1}{N_{\mathrm{b}}}{\bf a}_{{\rm B,0}}{\bf
a}_{{\rm B,0}}^{H}{\bf B}[n]{\bf A}_{{\rm M,}Q}({\bf A}_{{\rm M,}Q}^{H}{\bf
C}{\bf A}_{{\rm M,}Q})^{-1}{\bf A}_{{\rm M,}Q}^{H}{\bf x}_{t}$
$\displaystyle=$ $\displaystyle\frac{1}{N_{\mathrm{b}}}{\bf a}_{{\rm B,0}}{\bf
a}_{{\rm B,0}}^{H}{\bf B}[n]\mathbf{D}(\bm{\theta}_{\mathrm{t}}){\bf x}_{t},$
(111)

where

$\displaystyle{\bf B}[n]$
$\displaystyle\triangleq\sum\limits_{t=1}^{T_{1}}{\bf y}_{t}[n]{\bf
x}_{t}^{H}=[{\bf y}_{1}[n],{\bf y}_{2}[n],\cdots,{\bf
y}_{T_{1}}[n]]\left[\begin{array}[]{c}{\bf x}_{1}^{H}\\\ {\bf x}_{2}^{H}\\\
\vdots\\\ {\bf x}_{T_{1}}^{H}\end{array}\right]={\bf Y}_{1}[n]{\bf
X}_{1}^{H},$ (112e) $\displaystyle{\bf C}$
$\displaystyle\triangleq\sum\limits_{t=1}^{T_{1}}{\bf x}_{t}{\bf
x}_{t}^{H}={\bf X}_{1}{\bf X}_{1}^{H},$ (112f) $\displaystyle{\bf
D}(\bm{\theta}_{{\rm t}})$ $\displaystyle\triangleq{\bf A}_{{\rm M,}Q}({\bf
A}_{{\rm M,}Q}^{H}{\bf CA}_{{\rm M,}Q})^{-1}{\bf A}_{{\rm M,}Q}^{H}.$ (112g)

By substituting (111) into (78), and defining
$\mathbf{E}\triangleq\frac{1}{N_{\mathrm{b}}}{\bf a}_{{\rm B,0}}{\bf a}_{{\rm
B,0}}^{H}$, we have

$\displaystyle\sum\limits_{n=1}^{N}\sum\limits_{t=1}^{T_{1}}\left\|{\bf
y}_{t}[n]-{\tilde{{\bf X}}}_{t}[\sum\limits_{t=1}^{T_{1}}({\tilde{{\bf
X}}}_{t}^{H}{\tilde{{\bf X}}}_{t})]^{-1}\sum\limits_{t=1}^{T_{1}}({\tilde{{\bf
X}}}_{t}^{H}{\bf y}_{t}[n])\right\|_{2}^{2}$ $\displaystyle=$
$\displaystyle\sum\limits_{n=1}^{N}\sum\limits_{t=1}^{T_{1}}\left\|{\bf
y}_{t}[n]-\frac{1}{N_{\mathrm{b}}}{\bf a}_{{\rm B,0}}{\bf a}_{{\rm
B,0}}^{H}{\bf B}[n]{\bf D}(\bm{\theta}_{{\rm t}}){\bf x}_{t}\right\|_{2}^{2}$
$\displaystyle=$
$\displaystyle\sum\limits_{n=1}^{N}\sum\limits_{t=1}^{T_{1}}\left\|{\bf
y}_{t}[n]-{\bf EB}[n]{\bf D}(\bm{\theta}_{{\rm t}}){\bf
x}_{t}\right\|_{2}^{2}$ $\displaystyle=$
$\displaystyle\sum\limits_{n=1}^{N}\sum\limits_{t=1}^{T_{1}}({\bf
y}_{t}[n]-{\bf EB}[n]{\bf D}(\bm{\theta}_{{\rm t}}){\bf x}_{t})^{H}({\bf
y}_{t}[n]-{\bf EB}[n]{\bf D}(\bm{\theta}_{{\rm t}}){\bf x}_{t})$
$\displaystyle=$
$\displaystyle\sum\limits_{n=1}^{N}\sum\limits_{t=1}^{T_{1}}\\{{\bf
y}_{t}^{H}[n]{\bf y}_{t}[n]-2\mathfrak{Re}({\bf y}_{t}^{H}[n]{\bf EB}[n]{\bf
D}(\bm{\theta}_{{\rm t}}){\bf x}_{t})+{\bf x}_{t}^{H}{\bf
D}^{H}(\bm{\theta}_{{\rm t}})({\bf EB}[n])^{H}{\bf EB}[n]{\bf
D}(\bm{\theta}_{{\rm t}}){\bf x}_{t}\\}$ $\displaystyle\propto$
$\displaystyle\sum\limits_{n=1}^{N}\sum\limits_{t=1}^{T_{1}}\\{-2\mathfrak{Re}({\bf
y}_{t}^{H}[n]{\bf EB}[n]{\bf D}(\bm{\theta}_{{\rm t}}){\bf x}_{t})+{\bf
x}_{t}^{H}{\bf D}^{H}(\bm{\theta}_{{\rm t}})({\bf EB}[n])^{H}{\bf EB}[n]{\bf
D}(\bm{\theta}_{{\rm t}}){\bf x}_{t}\\}$ $\displaystyle=$
$\displaystyle-2\mathfrak{Re}\\{\sum\limits_{n=1}^{N}\sum\limits_{t=1}^{T_{1}}{\rm
Tr}({\bf D}(\bm{\theta}_{{\rm t}}){\bf x}_{t}{\bf y}_{t}^{H}[n]{\bf
EB}[n])\\}+\sum\limits_{n=1}^{N}\sum\limits_{t=1}^{T_{1}}{\rm Tr}[({\bf
EB}[n])^{H}{\bf EB}[n]{\bf D}(\bm{\theta}_{{\rm t}}){\bf x}_{t}{\bf
x}_{t}^{H}{\bf D}^{H}(\bm{\theta}_{{\rm t}})]$ $\displaystyle=$
$\displaystyle{\\!}-{\\!}2\mathfrak{Re}\\{{\rm Tr}({\bf D}(\bm{\theta}_{{\rm
t}})\sum\limits_{n=1}^{N}\sum\limits_{t=1}^{T_{1}}{\bf x}_{t}{\bf
y}_{t}^{H}[n]{\bf EB}[n])\\}{\\!}+{\\!}{\rm Tr}\\{[\sum\limits_{n=1}^{N}({\bf
EB}[n])^{H}{\bf EB}[n]]{\bf D}(\bm{\theta}_{{\rm
t}})[\sum\limits_{t=1}^{T_{1}}{\bf x}_{t}{\bf x}_{t}^{H}]{\bf
D}^{H}(\bm{\theta}_{{\rm t}})\\}$ $\displaystyle\mathop{{\rm=}}\limits^{(h)}$
$\displaystyle-2\mathfrak{Re}\\{{\rm Tr}({\bf D}(\bm{\theta}_{{\rm
t}})\sum\limits_{n=1}^{N}({\bf B}[n])^{H}{\bf EB}[n])\\}+{\rm
Tr}\\{(\sum\limits_{n=1}^{N}({\bf B}[n])^{H}{\bf EB}[n]){\bf
D}(\bm{\theta}_{{\rm t}})\mathbf{C}{\bf D}^{H}(\bm{\theta}_{{\rm t}})\\},$
(113)

where $\mathbf{E}^{H}\mathbf{E}=\mathbf{E}$ is utilized in step $(h)$.

##  Appendix B Simplifications and Derivations in the SAGE Algorithm

###  B-A Simplification of the ML Function for All Channel Parameters in (51)

By taking the natural logarithm of (50), and omitting the constant term, we
get

$\displaystyle{\rm ln}f({\bf y}[1],{\bf y}[2],\cdots,{\bf
y}[N]|\bm{{\bf\eta}})$ $\displaystyle\propto$
$\displaystyle-\sum\limits_{n=1}^{N}({\bf
y}[n]-\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}(\bm{{\Sigma}}_{q}^{T}\otimes{\tilde{{\bf
H}}}_{q}){\bf x}[n])^{H}({\bf
y}[n]-\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}(\bm{{\Sigma}}_{q}^{T}\otimes{\tilde{{\bf
H}}}_{q}){\bf x}[n])$ $\displaystyle=$
$\displaystyle-\sum\limits_{n=1}^{N}\left\|{\bf
y}[n]-\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}(\bm{{\Sigma}}_{q}^{T}\otimes{\tilde{{\bf
H}}}_{q}){\bf x}[n]\right\|_{2}^{2}$ $\displaystyle=$
$\displaystyle-\sum\limits_{n=1}^{N}\left\|{\bf
Y}[n]-\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q}\right\|_{F}^{2}$ $\displaystyle=$
$\displaystyle-\sum\limits_{n=1}^{N}{\rm tr}\left[({\bf
Y}[n]-\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q})^{H}({\bf
Y}[n]-\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q})\right]$ $\displaystyle=$
$\displaystyle-\sum\limits_{n=1}^{N}{\rm tr}({\bf Y}^{H}[n]{\bf
Y}[n])+\sum\limits_{n=1}^{N}{\rm
tr}\left[\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\bf
Y}^{H}[n]{\tilde{{\bf H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q}\right]$
$\displaystyle+\sum\limits_{n=1}^{N}{\rm
tr}\left[\sum\limits_{q=0}^{Q}\delta_{q}^{*}e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}\bm{{\Sigma}}_{q}^{H}{\bf
X}^{H}[n]{\tilde{{\bf H}}}_{q}^{H}{\bf Y}[n]\right]$
$\displaystyle-\sum\limits_{n=1}^{N}{\rm
tr}\left[(\sum\limits_{q=0}^{Q}\delta_{q}^{*}e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}\bm{{\Sigma}}_{q}^{H}{\bf
X}^{H}[n]{\tilde{{\bf
H}}}_{q}^{H})(\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q}))\right].$ (114)

By further omitting the constant terms in (114), the (114) can be simplified
into

$\displaystyle\Lambda(\bm{{\bf\eta}};\\{{\bf Y}[n]\\}_{n=1,2,\cdots,N})$
$\displaystyle=$ $\displaystyle\sum\limits_{n=1}^{N}{\rm
tr}\left[\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\bf
Y}^{H}[n]{\tilde{{\bf H}}}_{q}{\bf
X}[n]\bm{{\Sigma}}_{q}\right]+\sum\limits_{n=1}^{N}{\rm
tr}\left[\sum\limits_{q=0}^{Q}\delta_{q}^{*}e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}\bm{{\Sigma}}_{q}^{H}{\bf
X}^{H}[n]{\tilde{{\bf H}}}_{q}^{H}{\bf Y}[n]\right]$
$\displaystyle-\sum\limits_{n=1}^{N}{\rm
tr}\left[(\sum\limits_{q=0}^{Q}\delta_{q}^{*}e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}\bm{{\Sigma}}_{q}^{H}{\bf
X}^{H}[n]{\tilde{{\bf
H}}}_{q}^{H})(\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q}))\right]$
$\displaystyle\mathop{{\rm=}}\limits^{(i)}$ $\displaystyle
2{\mathfrak{Re}}\left\\{\sum\limits_{n=1}^{N}{\rm
tr}\left[\sum\limits_{q=0}^{Q}(\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\bf
Y}^{H}[n]{\tilde{{\bf H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q})\right]\right\\}$
$\displaystyle-\sum\limits_{n=1}^{N}{\rm
tr}(\sum\limits_{q_{1}=0}^{Q}\sum\limits_{q_{2}=0}^{Q}\delta_{q_{1}}^{*}\delta_{q_{2}}e^{j2\pi(\tau_{q_{1}}-\tau_{q_{2}})\frac{(n-1)B}{N}}\bm{{\Sigma}}_{q_{1}}^{H}{\bf
X}^{H}[n]{\tilde{{\bf H}}}_{q_{1}}^{H}{\tilde{{\bf H}}}_{q_{2}}{\bf
X}[n]\bm{{\Sigma}}_{q_{2}})$ $\displaystyle=$ $\displaystyle
2{\mathfrak{Re}}\left\\{\sum\limits_{n=1}^{N}\sum\limits_{q=0}^{Q}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\rm
tr(}{\bf Y}^{H}[n]{\tilde{{\bf H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q})\right\\}$
$\displaystyle-\sum\limits_{n=1}^{N}\sum\limits_{q_{1}=0}^{Q}\sum\limits_{q_{2}=0}^{Q}\delta_{q_{1}}^{*}\delta_{q_{2}}e^{j2\pi(\tau_{q_{1}}-\tau_{q_{2}})\frac{(n-1)B}{N}}{\rm
tr}(\bm{{\Sigma}}_{q_{1}}^{H}{\bf X}^{H}[n]{\tilde{{\bf
H}}}_{q_{1}}^{H}{\tilde{{\bf H}}}_{q_{2}}{\bf X}[n]\bm{{\Sigma}}_{q_{2}})$
$\displaystyle=$ $\displaystyle
2{\mathfrak{Re}}\left\\{\sum\limits_{n=1}^{N}\sum\limits_{q=0}^{Q}(\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\bf
a}_{{\rm M}}^{H}(\vartheta_{{\rm t},q}){\bf X}[n]\bm{{\Sigma}}_{q}{\bf
Y}^{H}[n]{\bf a}_{{\rm B}}(\vartheta_{{\rm r},0}))\right\\}$
$\displaystyle-\sum\limits_{n=1}^{N}\sum\limits_{q_{1}=0}^{Q}\sum\limits_{q_{2}=0}^{Q}(\delta_{q_{1}}^{*}\delta_{q_{2}}e^{j2\pi(\tau_{q_{1}}-\tau_{q_{2}})\frac{(n-1)B}{N}}{\bf
a}_{{\rm B}}^{H}(\vartheta_{{\rm r},0}){\bf a}_{{\rm B}}(\vartheta_{{\rm
r},0}){\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm t},q_{2}}){\bf
X}[n]\bm{{\Sigma}}_{q_{2}}\bm{{\Sigma}}_{q_{1}}^{H}{\bf X}^{H}[n]{\bf a}_{{\rm
M}}(\vartheta_{{\rm t},q_{1}}))$ $\displaystyle=$ $\displaystyle
2{\mathfrak{Re}}\left\\{\sum\limits_{q=0}^{Q}{\bf a}_{{\rm
M}}^{H}(\vartheta_{{\rm
t},q})\left(\sum\limits_{n=1}^{N}\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\bf
X}[n]\bm{{\Sigma}}_{q}{\bf Y}^{H}[n]\right){\bf a}_{{\rm B}}(\vartheta_{{\rm
r},0})\right\\}$ $\displaystyle-
N_{\mathrm{b}}\sum\limits_{q_{1}=0}^{Q}\sum\limits_{q_{2}=0}^{Q}\left\\{\delta_{q_{1}}^{*}\delta_{q_{2}}{\bf
a}_{{\rm M}}^{H}(\vartheta_{{\rm
t},q_{2}})\left(\sum\limits_{n=1}^{N}e^{j2\pi(\tau_{q_{1}}-\tau_{q_{2}})\frac{(n-1)B}{N}}{\bf
X}[n]\bm{{\Sigma}}_{q_{2}}\bm{{\Sigma}}_{q_{1}}^{H}{\bf X}^{H}[n]\right){\bf
a}_{{\rm M}}(\vartheta_{{\rm t},q_{1}})\right\\},$ (115)

where the step $(i)$ refers to the identity ${\rm tr}({\bf A}^{H})=[{\rm
tr}({\bf A})]^{*}$.

###  B-B Derivative of ML function of a Single Path with Respect to Channel
Gain in (61b)

The specific derivation process to obtain (61b) is given as follows.

$\displaystyle\frac{\partial L(\bm{{\bf\eta}}_{q};\\{{\hat{{\bf
Y}}}_{q}^{(k)}[n]\\}_{n=1,\cdots,N})}{\partial\delta_{q}^{*}}$
$\displaystyle\mathop{{\rm=}}\limits$
$\displaystyle-\frac{\partial\sum\limits_{n=1}^{N}\left\|{\hat{{\bf
Y}}}_{q}^{(k)}[n]-\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q}\right\|_{F}^{2}}{\partial\delta_{q}^{*}}$
$\displaystyle=$ $\displaystyle-\frac{\partial\sum\limits_{n=1}^{N}{\rm
tr}\left\\{({\hat{{\bf
Y}}}_{q}^{(k)}[n]-\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q})^{H}({\hat{{\bf
Y}}}_{q}^{(k)}[n]-\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q})\right\\}}{\partial\delta_{q}^{*}}$
$\displaystyle=$ $\displaystyle-\frac{\partial\sum\limits_{n=1}^{N}{\rm
tr}\left\\{({\hat{{\bf Y}}}_{q}^{(k)}[n])^{H}{\hat{{\bf
Y}}}_{q}^{(k)}[n]\right\\}}{\partial\delta_{q}^{*}}+\frac{\partial\sum\limits_{n=1}^{N}{\rm
tr}\left\\{\delta_{q}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}({\hat{{\bf
Y}}}_{q}^{(k)}[n])^{H}{\tilde{{\bf H}}}_{q}{\bf
X}[n]\bm{{\Sigma}}_{q}\right\\}}{\partial\delta_{q}^{*}}$
$\displaystyle+\frac{\partial\sum\limits_{n=1}^{N}{\rm
tr}\left\\{\delta_{q}^{*}e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}({\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q})^{H}{\hat{{\bf
Y}}}_{q}^{(k)}[n]\right\\}}{\partial\delta_{q}^{*}}$
$\displaystyle-\frac{\partial\sum\limits_{n=1}^{N}{\rm
tr}\left\\{\delta_{q}^{*}\delta_{q}e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}({\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q})^{H}{\tilde{{\bf H}}}_{q}{\bf
X}[n]\bm{{\Sigma}}_{q}\right\\}}{\partial\delta_{q}^{*}}$ $\displaystyle=$
$\displaystyle-\frac{\partial\sum\limits_{n=1}^{N}{\rm tr}\left\\{({\hat{{\bf
Y}}}_{q}^{(k)}[n])^{H}{\hat{{\bf
Y}}}_{q}^{(k)}[n]\right\\}}{\partial\delta_{q}^{*}}+\frac{\partial\delta_{q}\sum\limits_{n=1}^{N}{\rm
tr}\left\\{e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}({\hat{{\bf
Y}}}_{q}^{(k)}[n])^{H}{\tilde{{\bf H}}}_{q}{\bf
X}[n]\bm{{\Sigma}}_{q}\right\\}}{\partial\delta_{q}^{*}}$
$\displaystyle+\frac{\partial\delta_{q}^{*}\sum\limits_{n=1}^{N}{\rm
tr}\left\\{e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}({\tilde{{\bf H}}}_{q}{\bf
X}[n]\bm{{\Sigma}}_{q})^{H}{\hat{{\bf
Y}}}_{q}^{(k)}[n]\right\\}}{\partial\delta_{q}^{*}}-\frac{\partial\delta_{q}^{*}\delta_{q}\sum\limits_{n=1}^{N}{\rm
tr}\left\\{({\tilde{{\bf H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q})^{H}{\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q}\right\\}}{\partial\delta_{q}^{*}}$
$\displaystyle\mathop{{\rm=}}\limits^{(j)}$
$\displaystyle\sum\limits_{n=1}^{N}\left\\{e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}{\rm
tr}[({\tilde{{\bf H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q})^{H}{\hat{{\bf
Y}}}_{q}^{(k)}[n]]\right\\}-\delta_{q}\sum\limits_{n=1}^{N}{\rm
tr}[({\tilde{{\bf H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q})^{H}({\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q})],$ (116)

where the step $(j)$ is obtained by using
$\frac{\partial\delta}{\partial\delta^{*}}=0$,
$\frac{\partial\delta^{*}}{\partial\delta^{*}}=1$ and
$\frac{\partial\delta^{*}\delta}{\partial\delta^{*}}=\frac{\partial\delta^{*}}{\partial\delta^{*}}\delta+\delta^{*}\frac{\partial\delta}{\partial\delta^{*}}=\delta$,
which can be referred to the formula (3.4.7) on page 173 in [27].

###  B-C Simplification for the Expression of Channel Gain in (62)

The specific simplification process of ${\delta_{q}}$ in (62) is as follows.

$\displaystyle\delta_{q}$
$\displaystyle=\frac{\sum\limits_{n=1}^{N}\left\\{e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}{\rm
tr}(\bm{{\Sigma}}_{q}^{H}{\bf X}^{H}[n]{\tilde{{\bf H}}}_{q}^{H}{\hat{{\bf
Y}}}_{q}^{(k)}[n])\right\\}}{\sum\limits_{n=1}^{N}{\rm
tr}(\bm{{\Sigma}}_{q}^{H}{\bf X}^{H}[n]{\tilde{{\bf H}}}_{q}^{H}{\tilde{{\bf
H}}}_{q}{\bf X}[n]\bm{{\Sigma}}_{q}^{H})}$
$\displaystyle=\frac{\sum\limits_{n=1}^{N}\left\\{e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}{\rm
tr}({\bf a}_{{\rm M}}(\vartheta_{{\rm t},q}){\bf a}_{{\rm
B}}^{H}(\vartheta_{{\rm r},0}){\hat{{\bf
Y}}}_{q}^{(k)}[n]\bm{{\Sigma}}_{q}^{H}{\bf
X}^{H}[n])\right\\}}{\sum\limits_{n=1}^{N}{\rm tr}({\bf a}_{{\rm
M}}(\vartheta_{{\rm t},q}){\bf a}_{{\rm B}}^{H}(\vartheta_{{\rm r},0}){\bf
a}_{{\rm B}}(\vartheta_{{\rm r},0}){\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm
t},q}){\bf X}[n]\bm{{\Sigma}}_{q}\bm{{\Sigma}}_{q}^{H}{\bf X}^{H}[n])}$
$\displaystyle=\frac{\sum\limits_{n=1}^{N}\left\\{e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}{\bf
a}_{{\rm B}}^{H}(\vartheta_{{\rm r},0}){\hat{{\bf
Y}}}_{q}^{(k)}[n]\bm{{\Sigma}}_{q}^{H}{\bf X}^{H}[n]{\bf a}_{{\rm
M}}(\vartheta_{{\rm t},q})\right\\}}{\sum\limits_{n=1}^{N}[N_{{\rm b}}{\bf
a}_{{\rm M}}^{H}(\vartheta_{{\rm t},q}){\bf
X}[n]\bm{{\Sigma}}_{q}\bm{{\Sigma}}_{q}^{H}{\bf X}^{H}[n]{\bf a}_{{\rm
M}}(\vartheta_{{\rm t},q})]}$ $\displaystyle=\frac{{\bf a}_{{\rm
B}}^{H}(\vartheta_{{\rm
r},0})\left\\{\sum\limits_{n=1}^{N}(e^{j2\pi\tau_{q}\frac{(n-1)B}{N}}{\hat{{\bf
Y}}}_{q}^{(k)}[n]\bm{{\Sigma}}_{q}^{H}{\bf X}^{H}[n])\right\\}{\bf a}_{{\rm
M}}(\vartheta_{{\rm t},q})}{N_{{\rm b}}{\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm
t},q})\left\\{\sum\limits_{n=1}^{N}({\bf
X}[n]\bm{{\Sigma}}_{q}\bm{{\Sigma}}_{q}^{H}{\bf X}^{H}[n])\right\\}{\bf
a}_{{\rm M}}(\vartheta_{{\rm t},q})}.$ (117)

###  B-D Simplification of the ML Function in (63)

By substituting ${\hat{{\bf Y}}}_{q}^{(k)}[n]$ into (56), we can obtain

$\displaystyle L(\bm{{\bf\eta}}_{q};\\{{\hat{{\bf
Y}}}_{q}^{(k)}[n]\\}_{n=1,\cdots,N})=$ $\displaystyle
2{\mathfrak{Re}}\left[\delta_{q}{\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm
t},q})\left(\sum\limits_{n=1}^{N}e^{-j2\pi\tau_{q}\frac{(n-1)B}{N}}{\bf
X}[n]\bm{{\Sigma}}_{q}(\hat{\bf{Y}}_{q}^{(k)}[n])^{H}\right){\bf a}_{{\rm
B}}(\vartheta_{{\rm r},0})\right]$ $\displaystyle-N_{{\rm
b}}\delta_{q}^{*}\delta_{q}{\bf a}_{{\rm M}}^{H}(\vartheta_{{\rm
t},q})\left(\sum\limits_{n=1}^{N}{\bf
X}[n]\bm{{\Sigma}}_{q}\bm{{\Sigma}}_{q}^{H}{\bf X}^{H}[n]\right){\bf a}_{{\rm
M}}(\vartheta_{{\rm t},q}).$ (118)

Then, by substituting $\delta_{q}$ in (62) (i.e., (117)) into
$L(\bm{{\bf\eta}}_{q};\\{{\hat{{\bf Y}}}_{q}^{(k)}[n]\\}_{n=1,\cdots,N})$ in
(118), we can have the simplified expression of
$F(\bar{\bm{\eta}}_{q};\\{{\hat{{\bf Y}}}_{q}^{(k)}[n]\\}_{n=1,\cdots,N})$ as