# Near-Field Channel Estimation
in Mixed LoS/NLoS Environments
for Extremely Large-Scale MIMO Systems

Yu Lu, , and Linglong Dai All authors are with the Department of Electronic
Engineering, Tsinghua University, Beijing 100084, China, and also with the
Beijing National Research Center for Information Science and Technology
(BNRist), Beijing 100084, China. (e-mails<EMAIL_ADDRESS>daill@tsinghua.edu.cn).This work was supported in part by the National Key
Research and Development Program of China (Grant No. 2020YFB1807201) and in
part by the National Natural Science Foundation of China (Grant No. 62031019),
and in part by the European Commission through the H2020-MSCA-ITN META
WIRELESS Research Project under Grant 956256. (Corresponding author: Linglong
Dai.)

###### Abstract

Accurate channel model and channel estimation are essential to empower
extremely large-scale MIMO (XL-MIMO) in 6G networks with ultra-high spectral
efficiency. With the sharp increase in the antenna array aperture of the XL-
MIMO scenario, the electromagnetic propagation field will change from far-
field to near-field. Unfortunately, due to the near-field effect, most of the
existing XL-MIMO channel models fail to describe mixed line-of-sight (LoS) and
non-line-of-sight (NLoS) path components simultaneously. In this paper, a
mixed LoS/NLoS near-field XL-MIMO channel model is proposed to match the
practical near-field XL-MIMO scenario, where the LoS path component is modeled
by the geometric free space propagation assumption while NLoS path components
are modeled by the near-field array response vectors. Then, to define the
range of near-field for XL-MIMO, the MIMO Rayleigh distance (MIMO-RD) and MIMO
advanced RD (MIMO-ARD) is derived. Next, a two stage channel estimation
algorithm is proposed, where the LoS path component and NLoS path components
are estimated separately. Moreover, the Cramér-Rao lower bound (CRLB) of the
proposed algorithm is derived in this paper. Numerical simulation results
demonstrate that, the proposed two stage scheme is able to outperform the
existing methods in both the theoretical channel model and the QuaDRiGa
channel emulation platform.

###### Index Terms:

6G, extremely large-scale MIMO, channel estimation, near-field.

##  I Introduction

Due to the emergence of new applications such as digital twin, holographic
imaging, and extended reality, the spectral efficiency of 6G is expected to
grow tenfold [1, 2, 3]. To achieve higher spectral efficiency, the extremely
large-scale multiple-input multiple-output (XL-MIMO) is regarded as one of the
most important technologies for 6G [4, 5]. Compared with the massive MIMO
technology for 5G, the sharp increase of antenna aperture in XL-MIMO for 6G
induces the fundamental change of the electromagnetic field property. The
electromagnetic field can be divided into the far-field region and the near-
field region. The boundary between these two regions is the Rayleigh distance
(RD) $Z=2D^{2}/\lambda$ [6], which is defined in a multiple-input single-
output (MISO) scenario. The RD is proportional to the square of the array
aperture $D$ and the inverse of wavelength $\frac{1}{\lambda}$. Since the
antenna number in massive MIMO systems for 5G is not large enough (e.g., 64
antennas in [7]), the near-field range of up to a few tens of meters is
negligible. For example, with the carrier frequency at 28 GHz, the RD for an
array with array aperture $D$ as 0.1 $\rm m$ is only 1.9 $\rm m$ [8]. Thus,
users are likely located in the far-field region, where the far-field channel
can be modeled with planar wavefronts assumption. On the contrary, as the
antenna number dramatically increases in XL-MIMO systems (e.g., 1024 antennas
in [9]), the near-field range will expand by orders of magnitude, which can be
up to several hundreds of meters. For instance, with the array aperture as 1
$\rm m$ at 28 GHz, the RD is 187 $\rm m$, which nearly covers a region of a
cellular cell. In this case, the XL-MIMO channel should be modeled under near-
field assumption with spherical wavefronts.

In XL-MIMO communications, the transmitter equipped with an extremely large-
scale antenna array (ELAA) serves multiple users by exploiting the high
spatial multiplexing gain. Except for antenna array, reconfigurable
intelligent surface (RIS) is one of the main implementation methods of XL-MIMO
[3], which is an energy efficient alternative of ELAA. In order to obtain
spatial multiplexing gain, XL-MIMO should generate a directional beam with
high array gain by beamforming. Before realizing beamforming, the channel
state information (CSI) should be acquired in advance by channel estimation
[10, 11]. Since the electromagnetic field properties change from far-field to
near-field in the XL-MIMO, the existing far-field channel model mismatches the
practical near-field XL-MIMO channel feature. In this case, the existing far-
field channel estimation methods suffer from serious performance loss in the
near-field XL-MIMO channel model. Thus, it is important to delicately model
the XL-MIMO channel and design a near-field channel estimation scheme by means
of analyzing the propagation property in the near-field XL-MIMO scenario [12].

###  I-A Prior Contributions

The existing near-field channel estimation methods for ELAA based
communication systems can be divided into two categories, i.e., MISO channel
estimation [13, 14] and MIMO channel estimation. For the near-field MISO
scenario, the receiver with a single antenna is in the near-field region of
the transmitter with ELAA. Considering the near-field effects, the XL-MISO
channel should be modeled under the condition of the spherical wave assumption
instead of the planar wave assumption as far-field to ensure accuracy in the
near-field range [15, 14]. For example, [13] considers the scenario where the
transmitter employs an ELAA serving multiple single-antenna receivers. In this
scenario, the receivers are located in the transmitter’s near-field region.
Thus, the near-field MISO channel between the transmitter and receiver is
modeled based on a near-field array response vectors accurately, which relates
not only to the angle but also to the distance due to the spherical wave
assumption. To estimate the near-field channel, the whole two-dimensional
distance-angle plane is uniformly partitioned into multiple grids, and then
the corresponding near-field transform matrix can be constructed by multiple
near-field array response vectors associated with different grids. By means of
the constructed transform matrix, the near-field channel shows sparsity in the
transform domain, which can be estimated by compressive sensing (CS)-based
methods with low pilot overhead. The authors of [14] also consider the same
scenario as that in [13], interestingly, it is proved that the distance should
be non-uniformly divided to reduce the correlation among the near-field array
response vectors of the transform matrix in [14]. Based on the improved
transform matrix, a new sparse representation in polar-domain is proposed in
[14] for XL-MIMO. By utilizing this polar-domain sparsity, a CS-based
algorithm has been proposed to increase the estimation accuracy.

Figure 1: The XL-MIMO communication system with ELAA at the transmitter, where
the receiver is a multi-antenna user or a vehicle equipped with an ELAA.

The second category is based on the near-field XL-MIMO scenario, where the
transmitter employs an ELAA that serves receivers with multiple antennas or
even an ELAA as shown in Fig. 1. For instance, an ELAA is installed on the top
of a train [16, 17], and the transmitter also deploys an ELAA. Additionally,
ELAA can also be integrated into large infrastructures as revisers, such as
the roof of airports and the walls of stadiums. Additionally, XL-MIMO-based
wireless backhaul for the ultra-dense network (UDN) is also a promising use
case in the future 6G [18]. Specifically, in UDN, macrocell base stations
(BSs) and many ultra-dense small-cell BSs equipped with ELAAs cooperate to
provide a high data rate for users [19]. A premise to deploy UDN is the
reliable and large bandwidth backhaul connecting ultra-dense small-cell BSs
and macrocell BSs. In this case, the signal propagation between these two
kinds BSs should be modeled under near-field assumption instead of far-field
assumption. Therefore, the backhaul connecting in UDN is also a prospective
application for near-field XL-MIMO communications. There are two typical kinds
of methods for modeling near-field XL-MIMO channel, i.e., utilizing
approximation by the far-field channel model [20] and the product of the near-
field array response vectors similar to far-field MIMO model [21, 22, 23]. For
the first kind, in [20], the ELAA at transmitter and receiver are separated
into subarrays with a small aperture. The small aperture of the subarray leads
to a small near-field region. Thus, the subchannels between these subarrays
can be regarded as far-field channels due to the small apertures. In this
case, the whole channel between the transmitter and receiver can be
approximated by these far-field subchannels. However, the accuracy of this
method highly relies on the number of subarrays and the distance between the
transmitter and receiver. As for the second kind of method based on the near-
field array response vectors, the authors of [21] draw on the experience of
the far-field MIMO channel model and construct the near-field MIMO channel as
the product of the transmitter and receiver near-field array response vectors.
By utilizing the sparsity in polar-domain, the CS-based method can be applied
to the channel estimation problem under this channel model.

However, the existing near-field XL-MIMO channel model can only accurately
describe the non-line-of-sight (NLoS) path component of the near-field MIMO
while leading to an inaccurate description of the line-of-sight (LoS) path
component. The reasons are explained as below. 1) For the NLoS path component,
the signal radiates from the transmitter to the scatterers and then to the
receiver. Since the scatterers can be regarded as single-antenna transceivers,
the signal actually propagates through two MISO channels. Therefore, the NLoS
path component can be presented as the product of transmitter and receiver
array response vectors. This modeling method can apply to both far-field and
near-field NLoS components. 2) As for the LoS path component, in the far-field
scenario, the LoS path component is modeled in the same way as the NLoS path
component. This is because the transmitter with a small-aperture antenna array
can be viewed as a point from the receiver and vice versa, where the
transmitting and receiving processes are equivalent to two MISO propagation
processes. Thus, in the far-field scenario, the LoS path component can also be
represented by the product of transmitter and receiver array response vectors.
However, in a near-field XL-MIMO system, due to the large aperture of ELAA,
the receiver cannot be viewed as a point from the transmitter and vice versa,
and the transmitting and receiving processes cannot be approximated by the
previous theoretical model of two MISO propagation processes. In this case,
there exists a huge gap, which is defined as the accuracy loss, between the
propagation distance of the practical LoS path component and that of the
theoretical model based on two MISO processes. Intuitively, the accuracy loss
will increase as the array aperture increases, especially when the transmitter
and receiver employ the ELAAs simultaneously.

In this case, the existing near-field XL-MIMO channel model based on near-
field array response vectors mismatches the XL-MIMO scenario in practice.
Therefore, the existing near-field channel estimation schemes cannot estimate
the XL-MIMO channel accurately. Unfortunately, to the best of our knowledge,
there is no study of this problem for XL-MIMO in the current literature.

###  I-B Our Contributions

In order to solve this problem, an accurate mixed LoS/NLoS near-field XL-MIMO
channel model is first proposed. Then, a two stage XL-MIMO channel estimation
scheme is proposed for the mixed LoS/NLoS near-field XL-MIMO channel
model111Simulation codes are provided to reproduce the results in this paper:
http://oa.ee.tsinghua.edu.cn/dailinglong/publications/publications.html.. Our
specific contributions are listed below.

  1. 1.

In order to accurately describe the near-field XL-MIMO channel, a mixed
LoS/NLoS near-field XL-MIMO channel model is proposed. Specifically, the LoS
path components and NLoS path components are modeled separately. The NLoS path
components are modeled based on the near-field array response vectors, which
is similar to the far-field array response vectors based on the far-field
channel model. By contrast, the LoS path component is modeled based on the
precise geometric free space assumption.

  2. 2.

Based on the proposed mixed LoS/NLoS near-field XL-MIMO channel model, we
derive the MIMO Rayleigh distance (MIMO-RD) and MIMO advanced Rayleigh
distance (MIMO-ARD) to provide the near-field region for the XL-MIMO scenario.
In specific, similar to the derivation of MISO-RD, by considering the
condition that the maximum phase discrepancy between the proposed channel
model and the far-field array response vectors based XL-MIMO channel model in
the free space cannot exceed $\pi/8$, we first define the MIMO-RD, which is
proportional to the square of the sum of the antenna array apertures of the
transmitter and receiver. Then, by considering the condition that the maximum
phase discrepancy between the proposed channel model and the existing near-
field channel model in the free space is no more than $\pi/8$, the MIMO-ARD is
calculated, which is proportional to the product of the transmitter and
receiver array apertures. The derived MIMO-RD and MIMO-ARD show the bound of
the proposed channel model between the far-field array response vectors based
channel model and near-field array response vectors based channel model,
separately.

  3. 3.

We propose a two stage channel estimation algorithm for the proposed channel
model to estimate the XL-MIMO channel accurately. In our proposed two stage
channel estimation algorithm, the LoS and NLoS path components are estimated
separately. The LoS path component estimation is first realized by searching
collection with coarse on-grid parameters, and then refined by iteration
optimization, while the NLoS path components are estimated by orthogonal
matching pursuit (OMP)-based estimation with their polar-domain sparsity.
Moreover, the Cramér-Rao lower bound (CRLB) and the complexity of the proposed
channel estimation algorithm are derived in this paper. Finally, we provide
numerical simulation results to illustrate the effectiveness of our scheme in
both the proposed channel model and the QuaDRiGa channel emulation platform.

###  I-C Organization and Notation

Organization: The rest of the paper is organized as follows. In Section II, we
first introduce the signal model, and then review existing near-field XL-MISO
and XL-MIMO channel models. In Section III, we propose the mixed LoS/NLoS
near-field XL-MIMO channel model. Then, the MIMO-RD and MIMO-ARD are derived
to determine the range of the XL-MIMO near-field region in Section IV. In
Section V, we propose the corresponding channel estimation scheme based on the
proposed channel model, and analyze CRLB and its complexity. Simulation
results is provided in Section VI, and finally, conclusions are provided in
Section VII.

Notation: Lower-case and upper-case boldface letters ${\bf{a}}$ and ${\bf{A}}$
denote a vector and a matrix, respectively; $||{\bf{A}}||_{F}$ denotes the
Frobenius norm; ${{{\bf{a}}^{H}}}$ and ${{{\bf{A}}^{H}}}$ denote the conjugate
transpose of vector $\bf{a}$ and matrix $\bf{A}$, respectively. The circularly
symmetric complex Gaussian distribution is denotes by ${\cal
CN}\left(\mu,\sigma^{2}\right)$, with $\mu$ as mean set and $\sigma^{2}$ as
variance, and ${{\cal U}(-a,a)}$ denotes the uniform distribution on $(-a,a)$.
$\otimes$ denotes the Kronecker product. $\bm{I}$ denotes the identity matrix.
$\rm floor$ denotes round down operation.

##  II System Model

In this section, the signal model of the near-field XL-MIMO system used in
this paper will be introduced first. Then, we will review the existing near-
field XL-MISO and XL-MIMO channel models.

###  II-A Signal Model

In this work, we consider that the transmitter and receiver are equipped with
$N_{1}$-element and $N_{2}$-element antenna arrays, respectively. Since the
antenna arrays are usually implemented in a digital-analog hybrid manner with
a few RF (radio frequency) chains, we assume that RF chain numbers of the
transmitter and receiver are ${N_{\rm t}^{\rm RF}}$ and ${N_{\rm r}^{\rm
RF}}$. Let ${\bf{H}}\in\mathbb{C}^{{N_{2}}\times{N_{1}}}$ denotes the channel
from transmitter to receiver. The corresponding signal model can be presented
as

${\bf{y}}_{m}={\bf{W}}{\bf{H}}{\bf{Q}}{\bf{s}}_{m}+{\bf{n}}_{m},$ (1)

where ${\bf{y}}_{m}\in\mathbb{C}^{{N_{\rm r}^{\rm RF}}\times 1}$,
${\bf{W}}\in\mathbb{C}^{{N_{\rm r}^{\rm RF}}\times{N_{2}}}$,
${{\bf{Q}}}\in\mathbb{C}^{{N_{1}}\times{N_{\rm t}^{\rm RF}}}$ and
${\bf{s}}_{m}\in\mathbb{C}^{{N_{\rm t}^{\rm RF}}\times 1}$ denote the received
pilots signal, combining matrix, the hybrid precoding matrix, and transmitted
signal in $m$-th times slots, and ${{\bf{n}}_{m}}\sim{\cal C}{\cal
N}\left({{\bf{0}},\sigma^{2}{\bf{I}}_{N_{\rm r}^{\rm RF}}}\right)$ denotes the
${{N_{\rm r}^{\rm RF}}\times 1}$ received noise with ${\sigma^{2}}$
representing the noise power after combining in the $m$-th times slots.

Denote ${\bf p}_{m}={\bf{Q}}{\bf{s}}_{m}\in\mathbb{C}^{{N_{1}}\times 1}$,
where the $i$-th element of the ${\bf p}_{m}$ is the signal transmitted by the
$i$-th antenna at transmitter in $m$-th time slots. By collecting the received
pilots in $M$ time slots, we have

${\bf{Y}}={\bf{W}}{\bf{H}}{\bf{P}}+{\bf{N}},$ (2)

where ${\bf Y}=\left[{\bf{y}}_{1},{\bf{y}}_{2},\ldots,{\bf{y}}_{M}\right]$,
and ${\bf P}=\left[{\bf{p}}_{1},{\bf{p}}_{2},\ldots,{\bf{p}}_{M}\right]$ and
${\bf N}=\left[{\bf{n}}_{1},{\bf{n}}_{2},\ldots,{\bf{n}}_{M}\right]$.
${{\bf{Y}}}\in\mathbb{C}^{{N_{\rm r}^{\rm RF}}\times M}$,
${\bf{W}}\in\mathbb{C}^{{N_{\rm r}^{\rm RF}}\times{N_{2}}}$ and
${{\bf{P}}}\in\mathbb{C}^{{N_{1}}\times M}$ denote the received pilots signal,
combining matrix and the transmitted pilots signal in $M$ times slots in a
coherence interval, and ${{\bf{N}}}\sim{\cal C}{\cal
N}\left({{\bf{0}},\sigma^{2}{\bf{I}}_{{N_{\rm r}^{\rm
RF}}}\otimes{\bf{I}}_{M}}\right)$ denotes the ${{N_{2}}\times M}$ received
noise with ${\sigma^{2}}$ representing the noise power in $M$ times slots. In
channel estimation problem, we need to estimate $\bf{H}$ with given
${\bf{P}}$, ${\bf{W}}$ and ${\bf{Y}}$. The number of transmitter $N_{1}$ is
usually large in the XL-MIMO system. Thus, to reduce the pilot overhead in a
practical communication system, the channel estimation scheme with low
overhead should be utilized $\left(M<N_{1}\right)$. Since the channel model is
significant for designing a channel estimation scheme, we will review the
current XL-MIMO channel models next.

###  II-B Existing Near-Field XL-MISO Channel Model

Most of the existing near-field channel estimation works consider that only
the transmitter employs the ELAA, while the receivers are usually equipped
with a single antenna, i.e., near-field XL-MISO scenario. The near-field XL-
MISO channel between the receiver and the transmitter
${\bf{h}}_{\rm{n\mbox{-}f}}$ can be represented as

${\bf{h}}_{\rm{n\mbox{-}f}}=\sqrt{\frac{N_{1}}{L}}\sum\limits_{l=1}^{L}{\alpha_{l}}{\bf{b}}\left({\theta}_{l},r_{l}\right).$
(3)

The near-field array response vector ${\bf{b}}\left({\theta}_{l},r_{l}\right)$
in (3) is derived on the base of spherical wave assumption instead of planar
wave assumption, which can presented as [14]

${\bf{b}}(\theta_{l},r_{l})=\frac{1}{\sqrt{N_{1}}}[e^{-j{\frac{2\pi}{\lambda}}(r_{l}^{(1)}-r_{l})},\cdots,e^{-j{\frac{2\pi}{\lambda}}(r_{l}^{(N_{1})}-r_{l})}]^{H},$
(4)

where $r_{l}$ denotes the distance of the $l$-th scatterer from the center of
the transmitter antenna array,
$r_{l}^{({n_{1}})}=\sqrt{r_{l}^{2}+\delta_{n_{1}}^{2}d^{2}-2r_{l}\delta_{n_{1}}d\sin\theta_{l}}$
represents the distance of the $l$-th scatterer from the $n_{1}$-th
transmitter antenna, $\delta_{n_{1}}=\frac{2{n_{1}}-{N_{1}}-1}{2}$ with
${n_{1}}=1,2,\cdots,{N_{1}}$, $\theta_{l}\in(-\pi/2,\pi/2)$ are the practical
physical angles, $d$ is the antenna spacing of half wavelength to avoid
coupling between the antennas. It is worth pointing out that a compact array
with below half a wavelength spacing can occur in XL-MIMO, which will cause
coupling between the antennas. In this case, the coupling effect should be
measured first and then modeled as a coupling coefficient matrix with tunable
loads [24] in the channel model.

By utilizing the Taylor expansion $\sqrt{(1+x)}\approx
1+\frac{1}{2}x-\frac{1}{8}x^{2}+\mathcal{O}(x^{2})$, the distance difference
$r_{l}^{({n_{1}})}-r_{l}$ in (4) can be approximated by:

$r_{l}^{({n_{1}})}-r_{l}\approx=-\delta_{n_{1}}d{\sin\theta_{l}}+\delta_{n_{1}}^{2}d^{2}\frac{{\cos^{2}\theta_{l}}}{2{r_{l}}}.$
(5)

To exploit the sparsity of near-field channel, [14] proposed a new transform
matrix $\bf{D}$ to change the channel ${{\bf{h}}_{{\rm{n\mbox{-}f}}}}$ in (3)
into polar-domain, which can be presented as

$\displaystyle{\bf{D}}=[{\bf{b}}(\theta_{1},r_{1}^{1}),\cdots,{\bf{b}}(\theta_{1},r_{1}^{S_{1}}),\cdots,$
(6)
$\displaystyle{\bf{b}}(\theta_{N_{1}},r_{N_{1}}^{1}),\cdots,{\bf{b}}(\theta_{N_{1}},r_{N_{1}}^{S_{N_{1}}})],$

where each column of polar-domain transform matrix $\bf{D}$ is a near-field
array response vector sampled on the grid ($\theta_{n_{1}}$,
$r_{{n_{1}}}^{s_{n_{1}}}$), with $s_{n_{1}}=1,2,\cdots,{S_{n_{1}}}$,
$S_{n_{1}}$ denotes the number of sampled distances at $\theta_{{n_{1}}}$.
Therefore, the number of total sampled grids of the whole propagation
environment is $S=\sum\limits_{{n_{1}}=1}^{N_{1}}S_{n_{1}}$. With help of
$\bf{D}$, ${\bf{h}}_{\rm{n\mbox{-}f}}$ in (3) can be presented as

${\bf{h}}_{\rm{n\mbox{-}f}}={\bf{D}}{{\bf{h}}_{\rm{n\mbox{-}f}}^{P}},$ (7)

where ${{\bf{h}}_{\rm{n\mbox{-}f}}^{P}}$ is the channel represented in polar-
domain of size $S\times 1$. This near-field ${{\bf{h}}_{\rm{n\mbox{-}f}}^{P}}$
in (7) shows certain sparsity in the polar-domain. The authors in [14]
proposed the associated polar-domain based CS algorithm to solve the near-
field channel estimation problem.

###  II-C Existing Near-Field XL-MIMO Channel Models

The existing near-field XL-MIMO channel model refers to the far-field MIMO
channel model. Specifically, the far-field channel model is based on the
production of the far-field array response vectors. The near-field array
response vectors are utilized to replace the far-field array response vectors
to model near-field XL-MIMO channel [21], which is presented as

${\bf{H}}_{\rm{n\mbox{-}f}}=\sum\limits_{l=1}^{L}{{g_{l}}{\bf{b}}(\theta_{r}^{l},d_{r}^{l})}{{\bf{b}}^{H}}(\theta_{t}^{l},d_{t}^{l}),$
(8)

where $L$ denotes the number of path components, ${g_{l}}$ represents the
complex gain. ${\bf{b}}(\theta_{t}^{l},d_{t}^{l})$ and
${\bf{b}}(\theta_{r}^{l},d_{r}^{l})$ denote the near-field array response
vectors at transmitter and receiver on the base of spherical wave assumption,
which is denoted by

$\displaystyle{\bf{b}}(\theta_{t}^{l},d_{t}^{l})=\frac{1}{\sqrt{N_{1}}}[e^{-j{\frac{2\pi}{\lambda}}(d_{t}^{l}(1)-d_{t}^{l})},\cdots,e^{-j{\frac{2\pi}{\lambda}}(d_{t}^{l}(N_{1})-d_{t}^{l})}]^{H},$
(9)
$\displaystyle{\bf{b}}(\theta_{r}^{l},d_{r}^{l})=\frac{1}{\sqrt{N_{2}}}[e^{-j{\frac{2\pi}{\lambda}}(d_{r}^{l}(1)-d_{r}^{l})},\cdots,e^{-j{\frac{2\pi}{\lambda}}(d_{r}^{l}(N_{2})-d_{r}^{l})}]^{H},$

where $\theta_{t}^{l}\ (\theta_{r}^{l})$ represent the angle for the $l$-th
path at transmitter (receiver), and $d_{t}^{l}\ (d_{r}^{l})$ represent the
distance of the $l$-th scatterer from the center of the antenna array of
transmitter (receiver) for the $l$-th path,
$d_{t}^{l}({n_{1}})=\sqrt{{d_{t}^{l}}^{2}+\delta_{n_{1}}^{2}d^{2}-2{d_{t}^{l}}\delta_{n_{1}}d{\sin\theta_{t}^{l}}}$
represents the distance of the $l$-th scatterer from the $n_{1}$-th element on
transmitter antenna array, and $\delta_{n_{1}}=\frac{2{n_{1}}-N_{1}-1}{2}$
with $n=1,2,\cdots,N_{1}$, and
$d_{r}^{l}(n_{2})=\sqrt{{d_{r}^{l}}^{2}+\delta_{n_{2}}^{2}d^{2}-2{d_{r}^{l}}\delta_{n_{2}}d{\sin\theta_{r}^{l}}}$
represents the distance of the $l$-th scatterer from the $n_{2}$-th
transmitter antenna array, and $\delta_{n_{2}}=\frac{2{n_{2}}-N_{2}-1}{2}$
with $n=1,2,\cdots,N_{2}$.

The polar-domain transform matrix mentioned in (6) can be presented as

$\displaystyle{\bf{D}}_{t}=[{\bf{b}}(\theta_{1},d_{1}^{1}),\cdots,{\bf{b}}(\theta_{1},d_{1}^{S_{1}}),\cdots,$
(10)
$\displaystyle{\bf{b}}(\theta_{N_{1}},d_{N_{1}}^{1}),\cdots,{\bf{b}}(\theta_{N_{1}},r_{N_{1}}^{S_{N_{1}}})],$
$\displaystyle{\bf{D}}_{r}=[{\bf{b}}(\theta_{1},r_{1}^{1}),\cdots,{\bf{b}}(\theta_{1},d_{1}^{S_{1}}),\cdots,$
$\displaystyle{\bf{b}}(\theta_{{N_{2}}},r_{{N_{2}}}^{1}),\cdots,{\bf{b}}(\theta_{{N_{2}}},r_{{N_{2}}}^{S_{N_{2}}})],$

where each column of the matrix ${\bf{D}}_{t}$ (${\bf{D}}_{r}$) is a near-
field array response vector sampled at angle $\theta_{n_{1}}(\theta_{n_{2}})$
and distance $d_{n_{1}}^{s_{n_{1}}}(d_{n_{2}}^{s_{n_{2}}})$, with
$s_{n}=1,2,\cdots,{S_{n_{1}}}({S_{n_{2}}})$. ${S_{n_{1}}}({S_{n_{2}}})$
denotes the number of sampled distances at the sampled angle
$\theta_{n_{1}}(\theta_{n_{2}})$. Therefore, we can calculate the total number
of all sampled grids, i.e., the number of ${\bf{D}}_{t}$ (${\bf{D}}_{r}$)
columns, which can be presented as
${S_{1}}=\sum\limits_{{n_{1}}=1}^{N}S_{n_{1}}({S_{2}}=\sum\limits_{{n_{2}}=1}^{N}S_{n_{2}})$.

On the base of this polar-domain transform matrix ${\bf{D}}_{t}$ and
${\bf{D}}_{r}$, the channel ${{\bf{H}}}_{\rm{n\mbox{-}f}}$ can be represented
by

${{\bf{H}}}_{\rm{n\mbox{-}f}}={\bf{D}}_{r}{{\bf{H}}}^{P}_{\rm{n\mbox{-}f}}{\bf{D}}_{t}^{H},$
(11)

where ${{\bf{H}}}^{P}_{\rm{n\mbox{-}f}}$ is the $S_{2}\times S_{1}$ polar-
domain XL-MIMO channel, which also shows sparsity in the polar-domain.

In the existing near-field XL-MIMO channel model above, the LoS path component
is modeled as the same method for NLoS path components, where the MIMO channel
in (8) can be modeled based on two MISO channels. This method is not suitable
for the LoS path component in a near-field XL-MIMO scenario. Therefore, the
existing near-field channel model will mismatch the LoS path component of the
near-field XL-MIMO. Thus, the existing near-field channel estimation
algorithms cannot be utilized to address the XL-MIMO channel estimation
problem accurately. In order to design the channel estimation algorithm for
the XL-MIMO, we should first provide an accurate description of the LoS path
component, which will be described in Section III.

##  III The Proposed Mixed LoS/NLoS Near-Field XL-MIMO Channel Model

In this section, we will first utilize the free space propagation assumption
to accurately model the LoS path component for each transmitter-receiver
antenna pair. Then, the mixed LoS/NLoS near-field XL-MIMO channel model is
provided to capture the different features of LoS and NLoS path components,
which are modeled separately.

Figure 2: The proposed mixed LoS/NLoS near-field channel model for XL-MIMO.

###  III-A LoS Path Component

For the near-field XL-MIMO LoS path component, each transmitter-receiver
antenna pair will experience different propagation paths as shown in Fig. 2.
Thus, in this case, the channel model built with the near-field array response
vectors mismatches the practical feature of XL-MIMO near-field LoS path
component. As a result, instead of utilizing near-field array response vectors
like NLoS path components, we model the LoS path component under the geometric
free space propagation assumption [25] for each transmitter-receiver antenna
pair. Specifically, ${{\bf{H}}}({n_{2}},{n_{1}})$ denotes the LoS path
component of channel between the transmitter’s $n_{1}$-th antenna and the
receiver’s $n_{2}$-th antenna, which can be represented as

$\displaystyle{{\bf{H}}}({n_{2}},{n_{1}})$
$\displaystyle={{\frac{1}{{r_{{n_{2}},{n_{1}}}}}}e^{-j\frac{2\pi}{\lambda}{r_{{n_{2}},{n_{1}}}}}},$
(12)

where ${{r_{n_{2},n_{1}}}}$ denotes the distance of the $n_{1}$-th antenna at
receiver from the $n_{2}$-th antenna at transmitter. It is worth to point out
that the free space path loss of each transmitter-receiver antenna pair has
been normalized as ${\frac{1}{{r_{{n_{2}},{n_{1}}}}}}$ in this paper. The
${{r_{n_{2},n_{1}}}}$ can be represented as

$\displaystyle{r_{n_{2},n_{1}}}$ $\displaystyle=\sqrt{(r\cos\theta-
d_{2}\sin\phi)^{2}+(r\sin\theta+d_{2}\cos\phi-d_{1})^{2}}$ (13)
$\displaystyle=\sqrt{{r^{2}}\\!+\\!d_{1}^{2}\\!+\\!d_{2}^{2}\\!+\\!2\\!\left(r{d_{2}}\sin(\varphi\\!+\\!\theta)\\!\\!-\\!\\!r{d_{1}}\sin\theta\\!\\!-\\!\\!{d_{1}}{d_{2}}\cos\varphi\right)}$

where the $r$ is the distance of the 1-st antenna at receiver from the 1-st
antenna at transmitter, $\varphi$ denotes relative angle between receiver and
transmitter, and $\theta$ denotes the angle of departure (AoD) of the signal.
Moreover, $d_{1}$ and $d_{2}$ are denoted as $d_{1}=n_{1}d$ and
$d_{2}=n_{2}d$, where $d$ is antenna spacing. Thus, by utilizing geometry
relation in free space, the channel can be presented as

${{\bf{H}}_{{\rm{LoS}}}}=\mathbf{H}_{\mathrm{LoS}}(r,\theta,\varphi)={\left[{{\frac{1}{{r_{{n_{2}},{n_{1}}}}}}{e^{-j2\pi{r_{{n_{2}},{n_{1}}}}/\lambda}}}\right]_{{N_{2}}\times{N_{1}}}},$
(14)

Unlike the NLoS path components, the LoS path component cannot be decoupled by
near-field array response vectors. Thus, the LoS path component cannot be
presented by polar-domain channel with transform matrices.

###  III-B Proposed Mixed LoS/NLoS Near-Field XL-MIMO Channel

In order to capture both LoS and NLoS path components features, we propose a
mixed LoS/NLoS near-field XL-MIMO channel model based on (11) and (14). The
proposed XL-MIMO channel model can be presented:

$\displaystyle{\bf{H}}$
$\displaystyle={\bf{H}}_{{\rm{LoS}}}+{\bf{H}}_{{\rm{NLoS}}}$ (15)
$\displaystyle=\mathbf{H}_{\mathrm{LoS}}(r,\theta,\varphi)+{\bf{D}}_{r}\left(\sum_{l=1}^{L}g_{l}\mathbf{b}\left(\theta_{r}^{l},d_{r}^{l}\right)\mathbf{b}^{H}\left(\theta_{t}^{l},d_{t}^{l}\right)\right){\bf{D}}_{t}^{H}.$

It is worth noting that the LoS path component will degenerate into the same
representation as NLoS path components as the distance of the transmitter from
the receiver increases to a certain extent. Furthermore, if the distance
continues to increase, the proposed mixed LoS/NLoS near-field XL-MIMO channel
model will degenerate into a far-field MIMO channel model. To prove this
statement, in the next Section IV, we will derive the boundary between the
proposed mixed LoS/NLoS near-field XL-MIMO channel model and far-field MIMO
channel model and define the Rayleigh distance for the XL-MIMO scenario. Then,
the boundary between the proposed mixed LoS/NLoS near-field XL-MIMO and the
existing near-field XL-MIMO is provided and the advanced Rayleigh distance for
the XL-MIMO scenario is defined.

##  IV The Definitions of Rayleigh Distance for Near-Field XL-MIMO

In this section, we will first define the MIMO Rayleigh distance (MIMO-RD) to
determine the boundary between the proposed mixed LoS/NLoS near-field XL-MIMO
channel model and the far-field MIMO channel model. Then, MIMO advanced
Rayleigh distance (MIMO-ARD) is derived to determine the boundary between the
proposed mixed LoS/NLoS near-field XL-MIMO and the existing near-field XL-MIMO
channel model.

As described in [6], the RD for a MISO scenario, i.e., MISO-RD, is defined as
$Z=2D^{2}/\lambda$, where $D$ is the aperture of the antenna array and
$\lambda $ is the wavelength. MISO-RD is calculated by the condition that the
largest phase discrepancy between the far-field planar wavefront and the near-
field spherical wavefront in the free space is no more than $\pi/8$. However,
the current MISO-RD is identified based on the scenario where only the
transmitter employs the ELAA while the receiver is equipped with a single
antenna. In this case, the MISO-RD is only related to the aperture of the
transmitter ELAA. Thus, MISO-RD is not suitable for the XL-MIMO scenario.  In
this paper, by considering the XL-MIMO scenario, MIMO Rayleigh distance (MIMO-
RD) is defined by the condition that the largest phase discrepancy between the
far-field planar wavefronts and the near-field spherical wavefronts is no more
than $\pi/8$. Similar to MISO-RD, the proposed MIMO-RD describes the boundary
between the far-field and near-field regions for the XL-MIMO scenario.
However, MIMO-RD is unable to capture the feature of the largest phase
discrepancy between the proposed channel model and the existing near-field
array response vectors based channel model. Thus, we further define the MIMO-
ARD by the condition that the largest phase discrepancy between the proposed
mixed LoS/NLoS near-field XL-MIMO and the existing near-field MIMO scenario in
the free space is no more than $\pi/8$. The MIMO-RD and MIMO-RD are derived as
below.

Figure 3: The near-field MIMO scenario: transmitter and receiver are both
equipped with ELAAs with antenna array elements (orange squares) in the x-y
coordinate system

As shown in Fig. 3, we consider the scenario where the transmitter and
receiver are both equipped with antenna arrays of aperture $D_{1}$ and
$D_{2}$, respectively. These two antenna arrays are set in parallel since the
largest phase discrepancy occurs when the wave impinges perpendicularly [26].
The center of transmitter antenna arrays is regarded as the $x-y$ coordinate
origin. The coordinate of the $n_{1}$-th antenna of the transmitter antenna
array, the center of receiver antenna array, the $n_{2}$-th antenna of the
receiver antenna array are $(0,d_{1})$, $(x_{2},y_{2})$, and
$(x_{2},y_{2}+d_{2})$, respectively, where $-\frac{D_{1}}{2}\leq
d_{1}\leq\frac{D_{1}}{2}$ and $-\frac{D_{2}}{2}\leq d_{2}\leq\frac{D_{2}}{2}$.
The polar coordinate of the center of UE’s antenna array can be written as
$(r,\theta)=\left(\sqrt{x_{2}^{2}+y_{2}^{2}},\arctan(\frac{y_{2}}{x_{2}})\right)$,
where $r$, $\theta$ are the distance and angle between the center of
transmitter’s antenna array and the center of receiver’s antenna array,
respectively.

Figure 4: Three categories of XL-MIMO divided by the distance of the
transmitter from the receiver.

To derive the close form of MIMO-RD and MIMO-ARD, we only consider the phase
change caused by LoS path component of the channel. As mentioned in (12),
${{\bf{H}}}({n_{2}},{n_{1}})={{\frac{1}{{r_{{n_{2}},{n_{1}}}}}}{e^{-j2\pi{r_{{n_{2}},{n_{1}}}}/\lambda}}}$
represents the LoS path component of channel between the transmitter’s
$n_{1}$-th antenna and the receiver’s $n_{2}$-th antenna.
${r_{{n_{2}},{n_{1}}}}$ is the distance of the transmitter’s $n_{1}$-th
antenna from the receiver’s $n_{2}$-th antenna. Here, the true phase is
$\phi=\frac{2\pi}{\lambda}r_{{n_{2}},{n_{1}}}$, where ${r_{{n_{2}},{n_{1}}}}$
is presented as

$\displaystyle{r_{{n_{2}},{n_{1}}}}$
$\displaystyle=\sqrt{x_{2}^{2}+(y_{2}+d_{2}-d_{1})^{2}}$ (16)
$\displaystyle=r\sqrt{1+\frac{({d_{2}}-{d_{1}})^{2}}{r^{2}}+\frac{2({d_{2}}-{d_{1}})\sin\theta}{r}},$

where $\theta$ is the practical physical angle of departure. Specifically, by
utilizing the second-order Taylor expansion $\sqrt{1+x}\approx
1+\frac{1}{2}x-\frac{1}{8}x^{2}+\mathcal{O}(x^{2})$, we have

$\displaystyle{r_{{n_{2}},{n_{1}}}}$
$\displaystyle=r(1+\frac{({d_{2}}-{d_{1}})^{2}}{2r^{2}}+\frac{({d_{2}}-{d_{1}})\sin\theta}{r}$
(17)
$\displaystyle-\frac{1}{8}(\frac{2({d_{2}}-{d_{1}})\sin\theta}{r})^{2}+\mathcal{O}(\frac{1}{r^{2}}))$
$\displaystyle\approx
r+{(d_{2}-d_{1})}\sin\theta+\frac{(d_{2}-d_{1})^{2}\cos^{2}\theta}{2r}.$

###  IV-A MIMO Rayleigh Distance (MIMO-RD)

Based on the far-field assumption, ${r_{{n_{2}},{n_{1}}}}$ can be approximated
by its first-order Taylor expansion. Since $\sqrt{1+x}\approx 1+\frac{1}{2}x$,
we have

$\displaystyle r^{\text{far}}_{n_{2},n_{1}}$ $\displaystyle\approx
r(1+\frac{(d_{2}-d_{1})\sin\theta}{r})$
$\displaystyle=r+(d_{2}-d_{1})\sin\theta$ (18)

Thus, the far-field phase becomes
$\phi^{\text{far}}=\frac{2\pi}{\lambda}r^{\text{far}}_{n_{2},n_{1}}$.
Accordingly, the phase discrepancy between the far-field planar and near-field
spherical wavefronts can be presented as

$\Delta=\left|\phi_{n_{2},n_{1}}-\phi_{n_{2},n_{1}}^{\text{far}}\right|=\frac{2\pi}{\lambda}\left|r_{n_{2},n_{1}}-r_{n_{2},n_{1}}^{\text{far}}\right|$
(19)

Notice that $r_{n_{2},n_{1}}^{\text{far}}$ in (18) is the first-order Taylor
expansion of $r_{n_{2},n_{1}}$, so
$\left|r_{n_{2},n_{1}}-r_{n_{2},n_{1}}^{\text{far}}\right|$ in (24) is mainly
determined by _the second-order Taylor expansion term of $r_{n_{2},n_{1}}$_.
Therefore, the phase discrepancy between planar wavefronts and spherical
wavefronts is

$\displaystyle\Delta\approx\frac{\pi(d_{2}-d_{1})^{2}\cos^{2}\theta}{\lambda
r}.$ (20)

Since $\theta\in[-\frac{\pi}{2},\frac{\pi}{2}]$, $-\frac{D_{1}}{2}\leq
d_{1}\leq\frac{D_{1}}{2}$, and $-\frac{D_{2}}{2}\leq
d_{2}\leq\frac{D_{2}}{2}$, it can be easily observed that when $\theta=0$,
$d_{1}=\dfrac{D_{1}}{2}(-\dfrac{D_{1}}{2})$ and
$d_{2}=-\dfrac{D_{2}}{2}(\dfrac{D_{2}}{2})$, phase discrepancy achieves
maximum, where

$\displaystyle{\rm{max}}\frac{\pi(d_{2}-d_{1})^{2}\cos^{2}\theta}{\lambda
r}=\frac{\pi({D_{1}}+{D_{2}})^{2}}{4r\lambda}.$ (21)

Eventually, since the largest phase discrepancy is larger than $\pi/8$ in the
near-field region, the distance $r$ should satisfy
$\frac{\pi({D_{1}}+{D_{2}})^{2}}{4r\lambda}>\frac{\pi}{8}$, where

$r\leq\frac{2(D_{1}+D_{2})^{2}}{\lambda}.$ (22)

Thus, the MIMO-RD systems can be defined as
$\frac{2(D_{1}+D_{2})^{2}}{\lambda}$.

###  IV-B MIMO Advanced Rayleigh Distance (MIMO-ARD)

In the existing near-field XL-MIMO channel model, the LoS and NLoS path
components are both modeled based on the product of two near-field response
vectors. In this case, the phase of the signal from the $n_{1}$-th antenna of
the transmitter to the $n_{2}$-th antenna of the receiver can be acquired by
the sum of the phase of the $n_{1}$-th phase of the array response vector at
the transmitter and the $n_{2}$-th phase of array response vector at the
receiver. As delivered in (5), the phase of the near-field array response
vector can be approximated by its second-order Taylor expansion. We utilize
this approximation to calculate the sum of the phase of the $n_{1}$-th phase
of the transmitter array response vector and the $n_{2}$-th receiver phase of
the array response vector. Thus, the phase of the signal from the $n_{1}$-th
antenna of the transmitter to the $n_{2}$-th antenna of the receiver can be
represented as

$\displaystyle{\phi_{{n_{2}},{n_{1}}}^{\text{near}}}$ $\displaystyle\approx$
$\displaystyle=-{d_{1}}\sin\theta+\frac{\cos^{2}\theta{d_{1}}^{2}}{2r}+{d_{2}}\sin\theta+\frac{\cos^{2}\theta{d_{2}}^{2}}{2r}.$
(23)

Accordingly,based on (17) and (23), the phase discrepancy between the existing
near-field channel model and true near-field channel model can be presented as

$\displaystyle\Delta$
$\displaystyle=\left|\phi_{n_{2},n_{1}}-\phi_{n_{2},n_{1}}^{\text{near}}\right|$
$\displaystyle\approx\frac{2\pi}{\lambda}\left|\frac{\cos^{2}\theta{d_{1}}{d_{2}}}{r}\right|$
(24)

Since $\theta\in[-\frac{\pi}{2},\frac{\pi}{2}]$, $-\frac{D_{1}}{2}\leq
d_{1}\leq\frac{D_{1}}{2}$, and $-\frac{D_{2}}{2}\leq
d_{2}\leq\frac{D_{2}}{2}$, it can be easily observed that when $\theta=0$,
$d_{1}=\dfrac{D_{1}}{2}(-\dfrac{D_{1}}{2})$ and
$d_{2}=\dfrac{D_{2}}{2}(-\dfrac{D_{2}}{2})$, phase discrepancy achieves
maximum, where

$\displaystyle{\rm{max}}\frac{2\pi{\cos^{2}\theta{d_{1}}{d_{2}}}}{\lambda
r}=\frac{\pi{D_{1}}{D_{2}}}{2r\lambda}.$ (25)

Eventually, since the largest phase discrepancy is larger than $\pi/8$ in the
near-field region, the distance $r$ should satisfy
$\frac{\pi{D_{1}}{D_{2}}}{2r\lambda}>\frac{\pi}{8}$, where

$r<\frac{4D_{1}D_{2}}{\lambda}.$ (26)

Thus, the MIMO-ARD can be defined as $\frac{4D_{1}D_{2}}{\lambda}$.

From the derivation above, the XL-MIMO can be divided into three categories as
shown in Fig. 4, where the LoS path component of the XL-MIMO channel has
different features due to the different distances of the transmitter from the
receiver. Specifically, when the distance is bigger than the MIMO-RD, the
phase discrepancy between the existing near-field scenario and the far-field
scenario can be ignored. In this case, the LoS path component can be modeled
by the product of far-field array response vectors at the transmitter and
receiver. Additionally, when the distance of the transmitter from the receiver
is bigger than the MIMO-ARD and smaller than the MIMO-RD, the phase
discrepancy between the LoS path components of the proposed mixed LoS/NLoS
near-field channel and the existing near-field channel can be ignored. In this
case, the LoS path component can be modeled by the multiplication of near-
field response vectors at transmitter and receiver. Furthermore, when the
distance of the transmitter from the receiver is smaller than the MIMO-ARD,
the phase discrepancy between the LoS path components of the proposed mixed
LoS/NLoS near-field channel and the existing near-field channel cannot be
ignored.

In this paper, we mainly focus on the scenario where the distance of the
transmitter from the receiver is smaller than the MIMO-ARD. In the following
Section V, the specific channel estimation scheme will be described.

##  V Proposed Two Stage Channel Estimation Algorithm

In this section, based on the proposed mixed LoS/NLoS near-field XL-MIMO
channel model, we propose the two stage near-field XL-MIMO channel estimation
algorithm, where a parameter estimation algorithm is utilized for the LoS path
component estimation and an OMP-based algorithm is utilized for the NLoS path
components estimation. At last, we analyze the CRLB and the computational
complexity of the proposed two stage channel estimation algorithm.

###  V-A Stage 1: LoS Path Component Estimation

Since the energy of the LoS path component is usually dominant, we will first
conduct the LoS path component estimation. From (14), we can observe that the
LoS path component of near-field MIMO channel is determined by three
parameters, i.e., the distance of the 1-st antenna at receiver from the 1-st
antenna at transmitter $r$, relative angle between receiver and transmitter
$\varphi$, and the AoD $\theta$. Therefore, the LoS path component estimation
of (2) can be recognized as a parameter estimation problem, which can be
presented as

$\displaystyle\mathop{\min}\limits_{r,\theta,\varphi}G(r,\theta,\varphi)\buildrel\Delta\over{=}$
$\displaystyle{\left\|{{\bf{Y}}-\mathbf{W}\mathbf{H}_{\mathrm{LoS}}(r,\theta,\varphi)\bf{P}}\right\|^{2}_{F}}.$
(27)

Algorithm 1 shows the specific procedure to solve the problem (27).

Algorithm 1 LoS path component estimation

Inputs: Received signal ${\bf{Y}}$, pilot ${\bf{P}}$, $R_{max}$, $R_{min}$,
$\theta_{max}$, $\theta_{min}$, $\varphi_{max}$, $\varphi_{min}$, $r_{s}$,
$\theta_{s}$, $\varphi_{s}$, $I$, $\epsilon$.
Initialization: $\Delta r_{1}=\frac{r_{max}-r_{min}}{r_{s}}$,
$\Delta\theta_{1}=\frac{\theta_{max}-\theta_{min}}{\theta_{s}}$,
$\Delta\varphi_{1}=\frac{\varphi_{max}-\varphi_{min}}{\varphi_{s}}$, calculate
$\Xi$ based on (28)
// Estimate coarse on-grid parameters
1\. for $\left({r,\theta,\varphi}\right)\in\Xi$ do
2\. calculate ${\bf{H}}_{{\rm{LoS}}}^{*}$ based on (14)
3\.
${r^{{\rm{int}}},\theta^{{\rm{int}}},\varphi^{{\rm{int}}}}={\mathop{\rm{argmin}}\limits_{r,\theta,\varphi}}\|{\bf{Y}}-{\bf{H}}_{{\rm{LoS}}}{\bf{P}}\|^{2}_{2}$
4\. end for
// Refine off-grid parameters
5\. for $i\in I$ do
6\. update $\hat{{r}}^{(i+1)}$ based on (31)
7\. update $\hat{{\theta}}^{(i+1)}$ based on (32)
8\. update $\hat{{\theta}}^{(i+1)}$ based on (33)
9\. if
$\left|\frac{\hat{{r}}^{(i+1)}-\hat{{r}}^{(i)}}{\hat{{r}}^{(i+1)}}\right|$
$\leq$ $\epsilon$ and
$\left|\frac{\hat{{\theta}}^{(i+1)}-\hat{{\theta}}^{(i)}}{\hat{{\theta}}^{(i+1)}}\right|$
$\leq$ $\epsilon$ and
$\left|\frac{\hat{{\theta}}^{(i+1)}-\hat{{r}}^{(i)}}{\hat{{\theta}}^{(i+1)}}\right|$
$\leq$ $\epsilon$ then
10\. $({r^{\rm{opt}},{\theta}^{\rm{opt}},{\varphi}^{\rm{opt}}})$ =
$({\hat{{r}}^{(i+1)},{\hat{\theta}}^{(i+1)},{\hat{\varphi}}^{(i+1)}})$
11\. jump to Step 15
12\. end if
13.end for
14.$({r^{\rm{opt}},{\theta}^{\rm{opt}},{\varphi}^{\rm{opt}}})$ =
$({\hat{{r}}^{(I)},{\hat{\theta}}^{(I)},{\hat{\varphi}}^{(I)}})$
15.calculate ${\hat{{\bf{H}}}}_{{\rm{LoS}}}$ based on (14) by
$({r^{\rm{opt}},{\theta}^{\rm{opt}},{\varphi}^{\rm{opt}}})$
Output: Estimated LoS path component ${\hat{{\bf{H}}}}_{{\rm{LoS}}}$.

First, in Steps1-4 we obtain on-grid coarse parameters estimates by searching
the collection $\Xi$, which can be presented as

$\displaystyle\Xi=\\{\left({r,\theta,\varphi}\right)|r$
$\displaystyle=r_{\rm{min}},r_{\rm{min}}+\Delta r,\cdots,r_{\rm{max}};$ (28)
$\displaystyle\theta$
$\displaystyle=\theta_{\rm{min}},\theta_{\rm{min}}+\Delta\theta,\cdots,\theta_{\rm{max}};$
$\displaystyle\varphi$
$\displaystyle=\varphi_{\rm{min}},\varphi_{\rm{min}}+\Delta\varphi,\cdots,\varphi_{\rm{max}}\\},$

where $r_{\rm{min}}$, $r_{\rm{max}}$, $\theta_{\rm{min}}$,
$\theta_{\rm{max}}$, $\varphi_{\rm{min}}$, and $\varphi_{\rm{max}}$ represent
the lower and upper boundaries of the distance of the 1-st antenna at receiver
antenna array from the 1-st antenna at transmitter antenna array $r$, relative
angle between receiver and transmitter $\varphi$, and the AoD of the signal
$\theta$, respectively. $\Delta r$, $\Delta\theta$, and $\Delta\varphi$ are
the step sizes of $r$, $\theta$, and $\varphi$.

After searching the collection $\Xi$, we can get
${r^{\rm{int}},{\theta}^{\rm{int}},{\varphi}^{\rm{int}}}$ as the initial on-
grid estimated parameters. Then, to obtain the accurate estimated parameters,
we refine three parameters
${r^{\rm{int}},{\theta}^{\rm{int}},{\varphi}^{\rm{int}}}$ by iterative
optimization method. Specifically, we define the objective function
$G(r,\theta,\varphi)$ as

$\displaystyle\mathop{\min}\limits_{r,\theta,\varphi}G(r,\theta,\varphi)\buildrel\Delta\over{=}$
$\displaystyle{\left\|{{\bf{Y}}-\mathbf{W}\mathbf{H}_{\mathrm{LoS}}(r,\theta,\varphi)\bf{P}}\right\|^{2}_{F}},$
(29) $\displaystyle=$ $\displaystyle\sum_{m=1}^{M}\left({\bf{y}}_{m}-{\bf
W}{\bf{H}}{\bf{p}}_{m}\right)^{H}\left({\bf{y}}_{m}-{\bf
W}{\bf{H}}{\bf{p}}_{m}\right)$ $\displaystyle=$
$\displaystyle\sum_{m=1}^{M}{\bf{p}}_{m}^{H}{\bf{H}}^{H}{\bf W}^{H}{\bf
W}{\bf{H}}{\bf{p}}_{m}\\!\\!+\\!\\!\sum_{m=1}^{M}{\bf{y}}_{m}^{H}{\bf{y}}_{m}$
$\displaystyle-\\!\\!\sum_{m=1}^{M}\left({\bf{p}}_{m}^{H}{\bf{H}}^{H}{\bf
W}^{H}{\bf{y}}_{m}\\!\\!+\\!\\!\left({\bf{p}}_{m}^{H}{\bf{H}}^{H}{\bf
W}^{H}{\bf{y}}_{m}\right)^{H}\right)\\!\\!.\ $

Since the ${\bf{y}}_{m}$ is fixed received signal, the (29) above is equal to
the problem as

$\displaystyle\mathop{\min}\limits_{r,\theta,\varphi}G(r,\theta,\varphi)=\sum_{m=1}^{M}{\bf{p}}_{m}^{H}{\bf{H}}^{H}{\bf
W}^{H}{\bf W}{\bf{H}}{\bf{p}}_{m}$ (30)
$\displaystyle-\sum_{m=1}^{M}\left({\bf{p}}_{m}^{H}{\bf{H}}^{H}{\bf
W}^{H}{\bf{y}}_{m}+\left({\bf{p}}_{m}^{H}{\bf{H}}^{H}{\bf
W}^{H}{\bf{y}}_{m}\right)^{H}\right)$

The objective function $G(r,\theta,\varphi)$ can be optimized with an
iterative gradient descent approach methods. In the $i$-th iteration, we need
to calculate new estimated the three parameters, i.e., $\hat{r}^{(i+1)}$,
$\hat{\theta}^{(i+1)}$, and $\hat{\varphi}^{(i+1)}$, which can be updated as

$\displaystyle\hat{{r}}^{(i+1)}$
$\displaystyle=\hat{{r}}^{(i)}-\eta_{r}\cdot\nabla_{{r}}G_{\mathrm{opt}}^{(i)}\left(\hat{{r}}^{(i)},\hat{{\theta}}^{(i)},\hat{{\varphi}}^{(i)}\right),$
(31) $\displaystyle\hat{{\theta}}^{(i+1)}$
$\displaystyle=\hat{{\theta}}^{(i)}-\eta_{\theta}\cdot\nabla_{{\theta}}G_{\mathrm{opt}}^{(i)}\left(\hat{{r}}^{(i)},\hat{{\theta}}^{(i)},\hat{{\varphi}}^{(i)}\right),$
(32) $\displaystyle\hat{{\varphi}}^{(i+1)}$
$\displaystyle=\hat{{\varphi}}^{(i)}-\eta_{\varphi}\cdot\nabla_{{\varphi}}G_{\mathrm{opt}}^{(i)}\left(\hat{{r}}^{(i)},\hat{{\theta}}^{(i)},\hat{{\varphi}}^{(i)}\right),$
(33)

where the gradients can be calculated according to Appendix A. $\eta_{r}$,
$\eta_{\theta}$, $\eta_{\varphi}$ denote the lengths of step for the distance
and angles to guarantee the
$G_{\mathrm{opt}}^{(i)}\left(\hat{{r}}^{(i+1)},\hat{{\theta}}^{(i+1)},\hat{{\varphi}}^{(i+1)}\right)\leq
G_{\mathrm{opt}}^{(i)}\left(\hat{{r}}^{(i)},\hat{{\theta}}^{(i)},\hat{{\varphi}}^{(i)}\right)$.
As the number of iterations increases, the parameters estimates becomes more
accurate. Since only one variable in each iteration is optimized while the
other variables remain unchanged, the gradient descent method utilized in this
paper can be regarded as block coordinate descent. The iteration process will
suspend until the biggest iteration number achieves or the normalized
difference between the latest estimated parameters
$\left(\hat{{r}}^{(i+1)},\hat{{\theta}}^{(i+1)},\hat{{\varphi}}^{(i+1)}\right)$
and the last estimated parameters
$\left(\hat{{r}}^{(i)},\hat{{\theta}}^{(i)},\hat{{\varphi}}^{(i)}\right)$ are
smaller than $\epsilon$.

Based on the
$\left(\hat{{r}}^{(I)},\hat{{\theta}}^{(I)},\hat{{\varphi}}^{(I)}\right)$, we
can obtain the estimated ${\hat{{\bf{H}}}}_{{\rm{LoS}}}$ by (14). Then, we can
eliminate the influence of ${\hat{{\bf{H}}}}_{{\rm{LoS}}}$ on the received
pilots $\bf{Y}$ and then estimate the ${\hat{{\bf{H}}}}_{{\rm{NLoS}}}$, which
is shown as follows.

Algorithm 2 NLoS path components estimation

Inputs: ${\bf{Y}}_{{\rm{NLoS}}}$, ${\bf{P}}$, ${\bf{W}}$, ${\bf{D}}_{r}$,
${\bf{D}}_{t}$, $L$.
Initialization: ${\Omega_{1}}=\emptyset$, ${\Omega_{2}}=\emptyset$,
${\Omega}=\emptyset$, ${\bf A}={\bf{0}}_{S_{1}S_{2}\times 1}$,
${\bf{R}}={\bf{Y}}_{{\rm{NLoS}}}$, ${\bf A}_{t}={\bf D}_{t}^{H}{\bf P}$, ${\bf
A}_{r}={\bf W}{\bf D}_{r}$,
1\. for $l=1,2,\cdots,L$ do
2\. ${n^{*}}={\mathop{\rm{argmax}}}\|{\rm vec}({{\bf A}_{r}}^{H}{{\bf{R}}{{\bf
A}_{t}}^{H})\|}^{2}_{2}$
3\. $n_{1}={\rm floor}(({n}^{*}-1)/N_{2})+1$, $n_{2}=\mod({n}^{*}-1,N_{2})+1$
4\. ${\Omega}={\Omega}\bigcup{n}^{*}$, ${\Omega_{1}}={\Omega}\bigcup n_{1}$,
${\Omega_{2}}={\Omega}\bigcup n_{2}$
5\. ${\bf A}=\left[{\bf A}\ {\bf
A}_{t}\left(\Omega_{1},:\right)^{H}\otimes{\bf
A}_{r}\left(:,\Omega_{2}\right)\right]$
6\. ${\hat{{\bf{h}}}}_{A}({\Omega})=({\bf A}^{H}{\bf A})^{-1}{\bf A}^{H}{\rm
vec}({\bf Y}_{\rm NLoS})$
7\. reshape ${\hat{{\bf{h}}}}_{A}$ into ${\hat{{\bf{H}}}}_{A}$ of size
$S_{2}\times S_{1}$
8\.
${\bf{R}}={\bf{Y}}_{{\rm{NLoS}}}-{\bf{A}}_{r}{\hat{{\bf{H}}}}_{A}{\bf{A}}_{t}$
9\. end for
10\.
${\hat{{\bf{H}}}}_{{\rm{NLoS}}}={\bf{D}}_{r}{\hat{{\bf{H}}}}_{A}{\bf{D}}_{t}^{H}$
Output: Estimated NLoS path components ${\hat{{\bf{H}}}}_{{\rm{NLoS}}}$.

###  V-B Stage 2: NLoS Path Components Estimation

As we already acquire ${\hat{{\bf{H}}}}_{{\rm{LoS}}}$, the received pilots
without the effect of ${\hat{{\bf{H}}}}_{{\rm{LoS}}}$ can be presented as

${\bf{Y}}_{{\rm{NLoS}}}={\bf{Y}}-{\bf
W}{\hat{{\bf{H}}}}_{{\rm{LoS}}}{\bf{P}}.$ (34)

On the base of the polar-domain representation (11), the
${\bf{Y}}_{{\rm{NLoS}}}$ can be presented as

${\bf{Y}}_{{\rm{NLoS}}}={\bf W}{{\bf{H}}}_{{\rm{NLoS}}}{\bf{P}}+{\bf{N}}={\bf
W}{\bf{D}}_{r}{{\bf{H}}}_{{\rm{NLoS}}}^{P}{\bf{D}}_{t}^{H}{\bf{P}}+{\bf{N}}.$
(35)

As mentioned above, ${\bf{H}}_{{\rm{NLoS}}}^{P}$ is sparse in polar-domain,
thus the NLoS path components estimation is reformulated as a sparse recovery
problem. In this sparse recovery problem, the NLoS path components sensing
matrix at transmitter and receiver sides can be denoted as
${\bf{A}}_{t}={\bf{D}}_{t}^{H}{\bf{P}}$ and ${\bf{A}}_{r}={\bf
W}{\bf{D}}_{r}$. The low-complexity matrix based OMP algorithm to solve this
problem can be summarized in Algorithm 2.

Specifically, since there are $L$ components in the polar-domain, we will
conduct $L$ iterations to find $L$ supports in transmitter sensing matrix
${\bf{A}}_{t}$ and $L$ supports in receiver sensing matrix ${\bf{A}}_{r}$. In
$l$-th iteration, we will calculate the correlation between the transmitter
and receiver sensing matrices ${\bf{A}}_{t}$, ${\bf{A}}_{r}$ and the residual
matrix $\bf{R}$. In Step 4, we obtain the updated support $\Omega$,
$\Omega_{1}$, $\Omega_{2}$, where $\Omega_{1}$, $\Omega_{2}$ denote the
support of transmitter and receiver sides. Then, in Step 6, the currently
estimated near-field NLoS path component ${\hat{{\bf{h}}}}_{A}$ is calculated
by the least square (LS) algorithm. In Step 7, the ${\hat{{\bf{h}}}}_{A}$ need
to be reshaped into ${\hat{{\bf{H}}}}_{A}$ in the polar-domain of the size
$S_{2}\times S_{1}$. Finally, after $L$ iterations are performed, we obtain
the estimated NLoS path components ${\hat{{\bf{H}}}}_{{\rm{NLoS}}}$. Since
there exists an estimation error in the first stage, i.e., LoS path component
estimation, the error of NLoS path components estimation will contain the
error of LoS path component estimation.

After estimating the near-field LoS path component and the NLoS paths
components, the ${\hat{{\bf{H}}}}$ can be written as

${\hat{{\bf{H}}}}={\hat{{\bf{H}}}}_{{\rm{LoS}}}+{\hat{{\bf{H}}}}_{{\rm{NLoS}}}.$
(36)

###  V-C Cramér-Rao Lower Bound

The CRLB bound is a theoretical bound of MSE to evaluate channel estimation
algorithms [27]. For the channel estimation problem (2), we utilize
${\rm{vec}}({\bf{ABC}})=({\bf{C}}^{T}\otimes{\bf{A}}){\rm{vec}}({\bf{B}})$ to
reformulate (2) as

${\bf{y}}=({\bf{P}}^{T}\otimes{\bf{W}}){\bf{h}}+{\bf{n}}={\bf{Q}}{\bf{h}}+{\bf{n}}$
(37)

where ${\bf{y}}={\rm{vec}}({\bf{Y}})\in\mathbb{C}^{{N_{\rm r}^{\rm RF}}M\times
1}$, ${\bf{h}}={\rm{vec}}({{\bf{H}}})\in\mathbb{C}^{N_{1}N_{2}\times 1}$,
${\bf{n}}={\rm{vec}}({\bf{N}})\in\mathbb{C}^{{N_{\rm r}^{\rm RF}}M\times 1}$,
and ${\bf{Q}}=({\bf{P}}^{T}\otimes{\bf{W}})\in\mathbb{R}^{{N_{\rm r}^{\rm
RF}}M\times N_{1}N_{2}}$. Since $\bf Q$ is a real matrix and $\bf y,h,n$ are
complex vectors, we can split problem (37) into real part and imaginary part
as

$\displaystyle{\bf{y}}_{u}={\bf{Q}}{\bf{h}}_{u}+{\bf{n}}_{u}$ (38)
$\displaystyle{\bf{y}}_{v}={\bf{Q}}{\bf{h}}_{v}+{\bf{n}}_{u}$

where ${\bf{y}}_{u}={\rm Re}({\bf{y}})$, ${\bf{y}}_{v}={\rm Im}({\bf{y}})$,
${\bf{h}}_{u}={\rm Re}({\bf{h}})$, ${\bf{h}}_{v}={\rm Im}({\bf{h}})$,
${\bf{n}}_{u}={\rm Re}({\bf{n}})$, ${\bf{n}}_{v}={\rm Im}({\bf{n}})$. We
denote $\hat{{\bf h}}=\hat{{\bf h}}_{u}+\hat{{\bf h}}_{v}$ as the estimated
channel. The CRLB of the unbiased estimator ${\hat{\bf{h}}}$ can also be
divided into two parts, which are shown as

$\displaystyle{\rm CRLB}$ $\displaystyle={\rm CRLB}_{u}+{\rm CRLB}_{v}$ (39)
$\displaystyle=\mathrm{E}\left\\{\left\|\widehat{\mathbf{h}}-\mathbf{h}\right\|^{2}\right\\}$
$\displaystyle=\mathrm{E}\left\\{\left\|\widehat{\mathbf{h}}_{u}-\mathbf{h}_{u}\right\|^{2}\right\\}+\mathrm{E}\left\\{\left\|\widehat{\mathbf{h}}_{v}-\mathbf{h}_{v}\right\|^{2}\right\\}$

In this case, we first consider real part of (38). Since the ${\bf n}_{u}$
follows the distribution of Gaussian distribution with $\bf{0}$ mean and
${\bf{\sigma}}^{2}$ variance, the conditional probability density function of
${\bf y}_{u}$ with the given ${\bf h}_{u}$ is

$\begin{aligned}
p_{\mathbf{y}_{u}\mid\mathbf{h}_{u}}\left(\mathbf{y}_{u};\mathbf{h}_{u}\right)=\frac{1}{\left(2\pi\sigma^{2}\right)^{{N_{\rm
r}^{\rm
RF}}M/2}}\exp\left\\{-\frac{1}{2\sigma^{2}}\left\|{\bf{y}}_{u}-{\bf{Q}}{\bf{h}}_{u}\right\|^{2}\right\\}\end{aligned}.$
(40)

The Fisher information matrix of real part of (38) can then be derived as

$\displaystyle{[\mathbf{J}]_{m,n}}$
$\displaystyle\triangleq-\mathrm{E}\left\\{\frac{p_{\mathbf{y}_{u}\mid\mathbf{h}_{u}}\left(\mathbf{y}_{u};\mathbf{h}_{u}\right)}{\partial
h_{u,m}\partial h_{u,n}}\right\\}$
$\displaystyle=\frac{1}{\sigma^{2}}\left[\mathbf{Q}^{H}\mathbf{Q}\right]_{m,n},$
(41)

where $h_{u,m},h_{u,n}$ denote the $m$-th and $n$-th entry of
$\mathbf{h}_{u}$. Then, the real part ${\rm CRLB}_{u}$ is

$\displaystyle{\rm CRLB}_{u}$
$\displaystyle=\mathrm{E}\left\\{\left\|\widehat{\mathbf{h}}_{u}-\mathbf{h}_{u}\right\|^{2}\right\\}\geq{\rm
Tr}\left\\{{\bf J}_{u}^{-1}\right\\}=\sigma^{2}{\rm Tr}\left\\{({\bf
Q}^{H}{\bf Q})^{-1}\right\\}.$ (42)

Since ${\bf{Q}}={\bf{P}}^{T}\otimes{\bf{W}}$, $({\bf Q}^{H}{\bf Q})^{-1}$ can
be presented as

$\displaystyle({\bf Q}^{H}{\bf Q})^{-1}$
$\displaystyle=\left(({\bf{P}}^{T}\otimes{\bf{W}})^{H}({\bf{P}}^{T}\otimes{\bf{W}})\right)^{-1}$
(43) $\displaystyle=\left(\left({\bf P}{\bf P}^{H}\right)^{T}\otimes\left({\bf
W}^{H}{\bf W}\right)\right)^{-1}$ $\displaystyle=\left(\left({\bf P}{\bf
P}^{H}\right)^{-1}\right)^{T}\otimes\left({\bf W}^{H}{\bf W}\right)^{-1}.$

Thus, ${\rm Tr}\left\\{({\bf Q}^{H}{\bf Q})^{-1}\right\\}$ can be calculated
as

$\displaystyle{\rm Tr}\left(({\bf Q}^{H}{\bf Q})^{-1}\right)$
$\displaystyle={\rm Tr}\left(\left(\left({\bf P}{\bf
P}^{H}\right)^{-1}\right)^{T}\right){\rm Tr}\left(\left({\bf W}^{H}{\bf
W}\right)^{-1}\right)$ (44) $\displaystyle={\rm Tr}\left(\left({\bf P}{\bf
P}^{H}\right)^{-1}\right){\rm Tr}\left(\left({\bf W}^{H}{\bf
W}\right)^{-1}\right)$
$\displaystyle=\left(\sum_{i=1}^{N_{1}}\lambda_{i}^{-1}\right)\left(\sum_{j=1}^{N_{2}}\eta_{j}^{-1}\right)$
$\displaystyle\geq{N_{1}}\left({N_{1}}/\sum_{i=1}^{N_{1}}\lambda_{i}\right){N_{2}}\left({N_{2}}/\sum_{j=1}^{N_{2}}\eta_{j}\right)$
$\displaystyle=\frac{{N_{1}}^{2}}{\operatorname{Tr}\left\\{{\bf P}{\bf
P}^{H}\right\\}}\frac{{N_{2}}^{2}}{\operatorname{Tr}\left\\{{\bf W}^{H}\bf
W\right\\}}.$

${\left\\{\lambda_{i}\right\\}}_{i=1}^{N_{1}}$ and
${\left\\{\eta_{j}\right\\}}_{j=1}^{N_{2}}$ are denoted as the $N_{1}$ and
$N_{2}$ eigenvalues of the matrix of ${\bf P}{\bf P}^{H}$ and ${\bf W}^{H}\bf
W$. The equality in (44) holds when these eigenvalues satisfy
$\lambda_{1}=\lambda_{2}=\cdots\lambda_{N_{1}}$, and
$\eta_{1}=\eta_{2}=\cdots\eta_{N_{2}}$, i.e., the columns of ${\bf P}^{H}$
($\bf W$) and are orthogonal. Thus, the matrix ${\bf P}{\bf P}^{H}$ of size
$N_{1}\times N_{1}$ has identical diagonals equal to $M$, and ${\bf W}^{H}{\bf
W}$ of size $N_{2}\times N_{2}$ has identical diagonals equal to $N^{\rm
RF}_{\rm r}$. In this case, ${\rm Tr}\left\\{({\bf P}{\bf P})^{H}\right\\}{\rm
Tr}\left\\{({\bf W}^{H}\bf W)\right\\}=N_{1}MN_{2}N^{\rm RF}_{\rm r}$.
Finally, the CRLB of the real part of our estimation problem becomes

$\displaystyle{\rm CRLB}_{u}$
$\displaystyle=\mathrm{E}\left\\{\left\|\widehat{\mathbf{h}}_{u}-\mathbf{h}_{u}\right\|^{2}\right\\}=\sigma^{2}\dfrac{N_{1}N_{2}}{MN^{\rm
RF}_{\rm r}}.$ (45)

From the (38), we can observe that the real part and imaginary part have the
same form, ${\rm CRLB}_{v}={\rm CRLB}_{u}=\sigma^{2}\dfrac{N_{1}N_{2}}{MN^{\rm
RF}_{\rm r}}$. At last, thus the CRLB of the (37) is

${\rm CRLB}={\rm CRLB}_{u}+{\rm
CRLB}_{v}=2\sigma^{2}\dfrac{N_{1}N_{2}}{MN^{\rm RF}_{\rm r}}.$ (46)

It is noteworthy that the columns of pilot matrix ${\bf P}^{H}$ and combing
matrix ${\bf W}$ are not orthogonal in practical. Thus, the MSE of the
practical channel estimation algorithm cannot achieve the CRLB bound in (46),
which is verified in the simulation results in Section VI.

###  V-D Computational Complexity Analysis

For the proposed two stage channel estimation algorithm, we analyze its
computational complexity as follows. In the stage of the LoS path component
estimation, we can observe that the complexity comes from two parts, i.e.,
coarse on-grid estimation and refining processes.  In Steps 1-4, for coarse
on-grid estimation, we need to compute the ${\bf{H}}^{*}_{{\rm{LoS}}}$
according to (14) in parameters collection and find the best parameters.
Therefore, the complexity of this part is
$\mathcal{O}\left(S_{{\rm{LoS}}}\left(N_{\rm r}^{\rm RF}N_{1}N_{2}+N_{\rm
r}^{\rm RF}N_{1}M\right)\right)$, where $S_{{\rm{LoS}}}$ is the size of the
parameters collection. Since we only need a coarse on-grid estimation,
$S_{{\rm{LoS}}}$ usually is small. Then, for the refining process in Steps
5-9, the complexity is introduced by the gradient calculation. The complexity
to calculate gradients $\mathcal{O}(MI(N_{\rm r}^{\rm RF}N_{1}N_{2}+N_{\rm
r}^{\rm RF}N_{1}+N_{\rm r}^{\rm RF}N_{2}+N_{2}N_{1}+N_{\rm r}^{\rm RF}))$.
Since the $N_{1}$, $N_{2}$ is usually much larger than $N_{\rm r}^{\rm RF}$
and $M$, the complexity of the gradient calculation can be presented as
$\mathcal{O}\left(\left(S_{{\rm{LoS}}}+MI\right)N_{\rm r}^{\rm
RF}N_{1}N_{2}\right)$. For the NLoS path components estimation, the
computational complexity can be obtained as
$\mathcal{O}(N_{1}N_{2}(S_{1}+S_{2})L)$ by referring to the OMP algorithm
[28].

##  VI Simulation Result

In this section, we conduct the simulations to verify the performance of the
proposed two stage channel estimation algorithm for the proposed mixed
LoS/NLoS near-field XL-MIMO channel model. The system parameters are as
follows: the number of antenna of transmitter is $N_{1}=256$, the number of
antenna of receiver is $N_{2}=128$. The carrier frequency is $f=50$ GHz,
corresponding to $\lambda=0.006\,{\rm m}$. By utilizing (22), the MIMO-RD can
be calculated as
$\frac{2(D_{1}+D_{2})^{2}}{\lambda}=\frac{2({N_{1}\lambda/2}+N_{2}\lambda/2)^{2}}{\lambda}=442.7\,{\rm
m}$ in this scenario. By utilizing (26), the MIMO-ARD can be calculated as
$\frac{4D_{1}D_{2}}{\lambda}=\frac{4({N_{1}\lambda/2})(N_{2}\lambda/2)}{\lambda}=196.6\,{\rm
m}$ in this scenario. The near-field channel in (8) contains $L=3$ NLoS path
components. Meanwhile, the sampled angles of arrival follow the uniform
distribution $\mathcal{U}\left(-\frac{\pi}{3},\frac{\pi}{3}\right)$. The
distances are generated in range of [50, 500] meters. The pilot matrix
$\bf{P}$ and combining matrix $\bf W$ randomly chooses their elements from
$\left\\{-\frac{1}{\sqrt{M}},+\frac{1}{\sqrt{M}}\right\\}$ and
$\left\\{-\frac{1}{\sqrt{N_{2}}},+\frac{1}{\sqrt{N_{2}}}\right\\}$ [29].

TABLE I: Complexity Comparison  Method | Complexity
---|---
| Far-field codebook
---
based OMP
$\mathcal{O}(N_{1}N_{2}(N_{1}+N_{2})L)$
| Near-field codebook
---
based OMP
$\mathcal{O}(N_{1}N_{2}(S_{1}+S_{2})L)$
| Proposed two stage
---
channel estimation
$\mathcal{O}\\!\left(\left(S_{{\rm{LoS}}}\\!+\\!MI\right)N_{\rm r}^{\rm
RF}N_{1}N_{2}\\!+\\!N_{1}N_{2}(S_{1}\\!+\\!S_{2})L\right)$

We compare the proposed two stage channel estimation algorithm and the
existing near-field codebook based OMP method [14] and far-field codebook
based OMP method [28]. The complexities of the three methods are presented in
Table I. The normalized mean square error (NMSE) performance, which is defined
as
$\mathrm{NMSE}=\frac{\mathbb{E}(\|\mathbf{H}-\hat{\mathbf{H}}\|_{2}^{2})}{\mathbb{E}(\|\mathbf{H}\|_{2}^{2})}$,
is used to evaluate the accuracy of different methods. It worth to point out
that we use the $\frac{\rm CRLB}{{\mathbb{E}(\|\mathbf{H}\|_{2}^{2})}}$ as the
bound of NMSE performance.

Figure 5: NMSE performance comparison with respect to the distance of the
transmitter from the receiver.

Fig. 5 depicts the NMSE performance comparison with respect to the distance of
the transmitter from the receiver. The range of distance is from 50 ${\rm m}$
to 500 ${\rm m}$. The SNR is 5dB and the size of pilot matrix is $256\times
128$, the compressive ratio is $\dfrac{M}{N_{1}}=0.5$ [30]. The proposed two
stage method can achieve better NMSE performance than the existing far-field
codebook based OMP method and the near-field codebook based OMP method.
Specifically, the NMSE of the proposed algorithm is robust and remains the
lowest value of all the schemes in the whole range of distance. The NMSE
performance of the far-field codebook based OMP method and the near-field
codebook based OMP method degrades gradually with the decrease of the
distance. The reason is that there are still phase discrepancies between the
proposed channel and the existing far-field array response vectors based
channel as well as the existing near-field channel. In other words, the
existing two channel models mismatch the feature of the piratical near-field
XL-MIMO scenario. Furthermore, we can observe that the MIMO-RD can capture the
turning point of performance loss between the proposed two stage scheme and
the far-field codebook based OMP scheme. Specifically, when the distance is
larger than $442.7\,{\rm m}$, the performance of the proposed two stage scheme
and the far-field codebook based OMP channel estimation show the same
performance. The reason is that when the distance is larger than MIMO-RD, the
proposed channel model degenerates into the far-field channel model.
Similarly, MIMO-ARD can capture the turning point of performance loss between
the proposed two stage scheme and the near-field codebook based OMP scheme.
When the distance is larger than MIMO-ARD, i.e., $196.6\,{\rm m}$, the
performance of the proposed two stage scheme and the near-field codebook based
OMP channel estimation show the same performance. This is because when the
distance is larger than MIMO-ARD, the phase discrepancy between the proposed
channel model and the existing model vanishes.

Figure 6:  NMSE performance with respect to the size of the pilot. (a)
$r=60\,{\rm m}$; (b) $r=100\,{\rm m}$.

Fig. 6 depicts the NMSE performance comparison with respect to the size of the
pilot $\bf{P}$ under the condition of 5 dB SNR. The proposed scheme achieves
the best performance out of all three considered schemes. When the distance
becomes large as shown in Fig. 6 (b), the NMSE achieved by the two stage
channel estimation scheme and other schemes are similar. Since the fewer
pilots, the greater the error in the first stage of the LoS path component,
the gap between near-field codebook-based OMP and the proposed algorithm at
low pilot overhead is smaller than that of low pilot overhead. When the
distance becomes smaller as shown in Fig. 6 (a), the proposed scheme
outperforms all the other schemes.

Figure 7:  NMSE performance comparison with respect to the SNR under different
distances. (a) $r=60\,{\rm m}$; (b) $r=100\,{\rm m}$.

Fig. 7 shows the NMSE performance comparison with respect to the SNR under
different distances, where the size of pilot is $256\times 64$. In Fig. 7 (a)
and (b), the distance is $60\,{\rm m}$ and $100\,{\rm m}$, respectively. We
can see that when the distance is smaller than MIMO-ARD, the proposed two
stage channel estimation scheme outperforms all the other schemes. In
particular, when SNR is 5 dB and $r=60\,{\rm m}$, the proposed scheme can
achieve about 4 dB improvement compared with the near-field codebook based OMP
scheme. The reason is that the existing far-field and near-field codebook
based OMP schemes cannot deal with the LoS path component of near-field XL-
MIMO channel model when the distance is smaller than MIMO-ARD. Furthermore, we
can see that, in the Fig. 7 (b), the performance gap of NMSE of the three
schemes become smaller in the range of all the SNRs.

In order to evaluate the proposed two stage algorithm in the realistic
datasets, we utilize QuaDRiGa channel emulation platform [31, 32] as the
channel generator to test our proposed method. The standard channel cluster
model of 3GPP TR 38.901 [33] is used to generate the channel dataset. The main
parameters of the channel are set as shown in the Table. II.

TABLE II: Parameters of 3GPP TR 38.901 channel  Parameters | Value
---|---
Simulation scenario | UMa
Frequency | 3.5GHz
Number of antennas at transmitter | 64
Number of antennas at receiver | 256
Receiver location |  | Outdoor and Indoor
---
LoS and NLoS
Figure 8: NMSE performance comparison with respect to the SNR under the
QuaDRiGa channel dataset.

The specific simulation result with the GuaDRiGa channel dataset is shown in
Fig. 8. It can be observed that the proposed method still outperforms the
existing far-field codebook-based OMP algorithm and the near-field codebook
based OMP algorithm. Thus, the proposed mixed LoS/NLoS near-field MIMO channel
model is more accurate compared with the existing model. Furthermore, the
proposed two stage channel estimation scheme is an efficient scheme for the
practical near-field XL-MIMO.

##  VII Conclusions

In this paper, the channel estimation of the near-field XL-MIMO scenario was
investigated. We proposed the mixed LoS/NLoS near-field XL-MIMO channel model,
where the LoS and NLoS path components were characterized by geometric free
space assumption and the near-field response vectors, respectively. Then, we
derived the range of the near-field region of XL-MIMO, i.e., MIMO-RD and MIMO-
ARD. It is worth pointing out that the mixed LoS/NLoS near-field XL-MIMO
channel model proposed in this paper can be recognized as the generalization
of both the far-field channel model as well as the existing near-field channel
model for MIMO systems. Simulation results showed that, compared with the far-
field codebook based and near-field codebook based channel estimation schemes,
the proposed two stage channel estimation scheme achieved better NMSE
performance in both theoretical and practical channel models. For future work,
the uniform planar array (UPA) based 3D near-field XL-MIMO channel model will
attract further research interest.

##  Appendix A Optimization of $G^{(i)}$ in with regard to $r^{(i)}$,
$\theta^{(i)}$, $\varphi^{(i)}$

In this appendix, the derivations of the $G^{(i)}$ in with regard to
$r^{(i)}$, $\theta^{(i)}$, $\varphi^{(i)}$ are provided. For the
simplification of expression, we will ignore the superscript $(i)$.

For the distance $r$, the gradient of $G$ is given by

$\displaystyle\frac{{\partial{G}}}{{\partial{r}}}=$
$\displaystyle\sum_{m=1}^{M}{\bf{p}}_{m}^{H}\frac{\partial{\bf{H}}^{H}}{{\partial{r}}}{\bf
W}^{H}{\bf
W}{\bf{H}}{\bf{p}}_{m}+\sum_{m=1}^{M}{\bf{p}}_{m}^{H}{\bf{H}}^{H}{\bf
W}^{H}{\bf W}\frac{\partial{\bf{H}}}{{\partial{r}}}{\bf{p}}_{m}$ (47)
$\displaystyle-\sum_{m=1}^{M}{\bf{p}}_{m}^{H}\frac{\partial{\bf{H}}^{H}}{{\partial{r}}}{\bf
W}^{H}{\bf{y}}_{m}-\sum_{m=1}^{M}{\bf{y}}_{m}^{H}{\bf
W}\frac{\partial{\bf{H}}}{{\partial{r}}}{\bf{p}}_{m}.$

Since $r$ is a real variable,
$\frac{\partial{\bf{H}}^{H}}{{\partial{r}}}=\left(\frac{\partial{\bf{H}}}{{\partial{r}}}\right)^{H}$.
Thus, the (47) can be rewritten as

$\displaystyle\frac{{\partial{G}}}{{\partial{r}}}=\sum_{m=1}^{M}{\bf{p}}_{m}^{H}\left(\frac{\partial{\bf{H}}}{{\partial{r}}}\right)^{H}{\bf
W}^{H}\left({\bf W}{\bf{H}}{\bf{p}}_{m}-{\bf{y}}_{m}\right)$ (48)
$\displaystyle+\sum_{m=1}^{M}\left({\bf{p}}_{m}^{H}{\bf{H}}^{H}{\bf
W}^{H}-{\bf y}_{m}^{H}\right){\bf
W}\frac{\partial{\bf{H}}}{{\partial{r}}}{\bf{p}}_{m}.$

Similar to the derivation of $\frac{{\partial{G}}}{{\partial{r}}}$ above, the
$\frac{{\partial{G}}}{{\partial{\theta}}}$ and
$\frac{{\partial{G}}}{{\partial{\varphi}}}$ can be presented as

$\displaystyle\frac{{\partial{G}}}{{\partial{\theta}}}=\sum_{m=1}^{M}{\bf{p}}_{m}^{H}\left(\frac{\partial{\bf{H}}}{{\partial{\theta}}}\right)^{H}{\bf
W}^{H}\left({\bf W}{\bf{H}}{\bf{p}}_{m}-{\bf{y}}_{m}\right)$ (49)
$\displaystyle+\sum_{m=1}^{M}\left({\bf{p}}_{m}^{H}{\bf{H}}^{H}{\bf
W}^{H}-{\bf y}_{m}^{H}\right){\bf
W}\frac{\partial{\bf{H}}}{{\partial{\theta}}}{\bf{p}}_{m},$

$\displaystyle\frac{{\partial{G}}}{{\partial{\varphi}}}=\sum_{m=1}^{M}{\bf{p}}_{m}^{H}\left(\frac{\partial{\bf{H}}}{{\partial{\varphi}}}\right)^{H}{\bf
W}^{H}\left({\bf W}{\bf{H}}{\bf{p}}_{m}-{\bf{y}}_{m}\right)$ (50)
$\displaystyle+\sum_{m=1}^{M}\left({\bf{p}}_{m}^{H}{\bf{H}}^{H}{\bf
W}^{H}-{\bf y}_{m}^{H}\right){\bf
W}\frac{\partial{\bf{H}}}{{\partial{\varphi}}}{\bf{p}}_{m}.$

Based on the $\left(n_{2},n_{1}\right)$-th element in the $\bf H$, i.e., ${\bf
H}_{n_{2},n_{1}}$ $=$
${{\frac{1}{{r_{{n_{2}},{n_{1}}}}}}{e^{-j2\pi{r_{{n_{2}},{n_{1}}}}/\lambda}}}$,
where ${r_{n_{2},n_{1}}}$ is presented in (13), we can obtain the
$\left(n_{2},\\!n_{1}\right)$-th element in the $\frac{\partial{\bf
H}}{{\partial{r}}}$, $\frac{\partial{\bf H}}{{\partial{\theta}}}$, and
$\frac{\partial{\bf H}}{{\partial{\varphi}}}$ as

$\displaystyle\frac{{\partial{\bf H}_{n_{2},n_{1}}}}{{\partial{r}}}$
$\displaystyle=\left(d_{1}\sin\theta-
d_{2}\sin(\theta+\varphi)-r\right)\left(1\\!+\\!j{r_{{n_{2}},{n_{1}}}}\right)\Gamma,$
(51)

$\displaystyle\frac{{\partial{\bf H}_{n_{2},n_{1}}}}{{\partial{\theta}}}$
$\displaystyle=r\left(d_{1}\cos\theta-
d_{2}\cos(\theta+\varphi)\right)\left(1\\!+\\!j{r_{{n_{2}},{n_{1}}}}\right)\Gamma,$
(52)

$\displaystyle\frac{{\partial{\bf H}_{n_{2},n_{1}}}}{{\partial{\varphi}}}$
$\displaystyle=jd_{2}\left(d_{1}\sin\varphi+r\cos(\theta+\varphi)\right)\left({r_{{n_{2}},{n_{1}}}}\\!-\\!j\right)\Gamma.$
(53)

where $\Gamma=\frac{e^{-j{r_{{n_{2}},{n_{1}}}}}}{{r_{{n_{2}},{n_{1}}}^{3}}}$.

## References

  * [1] P. P. Ray, N. Kumar, and M. Guizani, “A vision on 6G-enabled NIB: Requirements, technologies, deployments, and prospects,” _IEEE Wireless Commun._ , vol. 28, no. 4, pp. 120–127, May 2021.
  * [2] M. Giordani, M. Polese, M. Mezzavilla, S. Rangan, and M. Zorzi, “Toward 6G networks: Use cases and technologies,” _IEEE Commun. Mag._ , vol. 58, no. 3, pp. 55–61, Mar. 2020.
  * [3] M. Cui, Z. Wu, Y. Lu, X. Wei, and L. Dai, “Near-field MIMO communications for 6G: Fundamentals, challenges, potentials, and future directions,” _IEEE Commun. Mag._ , Jan. 2023.
  * [4] E. D. Carvalho, A. Ali, A. Amiri, M. Angjelichinoski, and R. W. Heath, “Non-stationarities in extra-large-scale massive MIMO,” _IEEE Wireless Commun._ , vol. 27, no. 4, pp. 74–80, Aug. 2020.
  * [5] V. Croisfelt, A. Amiri, T. Abrao, E. de Carvalho, and P. Popovski, “Accelerated randomized methods for receiver design in extra-large scale MIMO arrays,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 7, pp. 6788–6799, May 2021.
  * [6] K. T. Selvan and R. Janaswamy, “Fraunhofer and Fresnel distances: Unified derivation for aperture antennas,” _IEEE Antennas Propag. Mag._ , vol. 59, no. 4, pp. 12–15, Aug. 2017.
  * [7] Z. Chen, F. Sohrabi, and W. Yu, “Multi-cell sparse activity detection for massive random access: Massive MIMO versus cooperative MIMO,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 8, pp. 4060–4074, Aug. 2019.
  * [8] A. Pizzo, L. Sanguinetti, and T. L. Marzetta, “Fourier plane-wave series expansion for holographic MIMO communications,” _IEEE Trans. Wireless Commun._ , Sept. 2022.
  * [9] Y. Zhu, H. Guo, and V. K. N. Lau, “Bayesian channel estimation in multi-user massive MIMO with extremely large antenna array,” _IEEE Trans. Signal Process._ , Sept. 2021.
  * [10] X. Ma, Z. Gao, F. Gao, and M. Di Renzo, “Model-driven deep learning based channel estimation and feedback for millimeter-wave massive hybrid MIMO systems,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 8, pp. 2388–2406, Jun. 2021.
  * [11] X. Li, X. Zhang, Y. Zhou, and L. Hanzo, “Optimal massive-MIMO-aided clustered base-station coordination,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 3, pp. 2699–2712, Feb. 2021.
  * [12] S. Sun, T. S. Rappaport, R. W. Heath, A. Nix, and S. Rangan, “MIMO for millimeter-wave wireless communications: beamforming, spatial multiplexing, or both?” _IEEE Commun. Mag._ , vol. 52, no. 12, pp. 110–121, Dec. 2014\.
  * [13] Y. Han, S. Jin, C. Wen, and X. Ma, “Channel estimation for extremely large-scale massive MIMO systems,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 5, pp. 633–637, May 2020.
  * [14] M. Cui and L. Dai, “Channel estimation for extremely large-scale MIMO: Far-field or near-field,” _IEEE Trans. Commun._ , vol. 70, no. 4, pp. 2663–2677, Apr. 2022.
  * [15] X. Yin, S. Wang, N. Zhang, and B. Ai, “Scatterer localization using large-scale antenna arrays based on a spherical wave-front parametric model,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 10, pp. 6543–6556, Jul. 2017.
  * [16] Z. Gong, C. Li, F. Jiang, and M. Z. Win, “Data-aided doppler compensation for high-speed railway communications over mmwave bands,” _IEEE Trans. on Wireless Commun._ , vol. 20, no. 1, pp. 520–534, Jan. 2021.
  * [17] B. Chen, Z. Zhong, B. Ai, and D. G. Michelson, “A geometry-based stochastic channel model for high-speed railway cutting scenarios,” _IEEE Antennas Wireless Propag. Lett._ , vol. 14, pp. 851–854, Dec. 2014.
  * [18] Q. Zhang, W. Ma, Z. Feng, and Z. Han, “Backhaul-capacity-aware interference mitigation framework in 6g cellular internet of things,” _IEEE Internet of Things Journal_ , Jun. 2021.
  * [19] Z. Gao, L. Dai, D. Mi, Z. Wang, M. A. Imran, and M. Z. Shakir, “Mmwave massive-MIMO-based wireless backhaul for the 5G ultra-dense network,” _IEEE Wireless Commun._ , vol. 22, no. 5, pp. 13–21, Oct. 2015.
  * [20] Y. Chen, L. Yan, and C. Han, “Hybrid spherical- and planar-wave modeling and DCNN-powered estimation of terahertz ultra-massive MIMO channels,” _IEEE Trans. Commun._ , vol. 69, no. 10, pp. 7063–7076, Oct. 2021.
  * [21] E. Zhou, H. Jiang, and H. Qi, “4-D parameter estimation in bistatic MIMO radar for near-field target localization,” in _Proc. 2015 IEEE International Wireless Symposium (IWS’2015)_ , 2015, pp. 1–4.
  * [22] P. R. Singh, Y. Wang, and P. Charge, “Near field targets localization using bistatic MIMO system with spherical wavefront based model,” in _Proc. 25th European Signal Processing Conference (EUSIPCO’2017)_ , 2017, pp. 2408–2412.
  * [23] L. Khamidullina, I. Podkurkov, and M. Haardt, “Conditional and unconditional cram r-rao bounds for near-field localization in bistatic MIMO radar systems,” _IEEE Trans. Signal Process._ , vol. 69, pp. 3220–3234, 2021.
  * [24] B. Clerckx, C. Craeye, D. Vanhoenacker-Janvier, and C. Oestges, “Impact of antenna coupling on 2 $\times$ 2 MIMO communications,” _IEEE Trans. Veh. Technol._ , May 2007.
  * [25] J. Sherman, “Properties of focused apertures in the Fresnel region,” _IEEE Trans. Antennas Propag._ , vol. 10, no. 4, pp. 399–408, Jul. 1962.
  * [26] E. Björnson, . T. Demir, and L. Sanguinetti, “A primer on near-field beamforming for arrays and reconfigurable intelligent surfaces,” in _Proc. 2021 55th Asilomar Conference on Signals, Systems, and Computers_ , 2021, pp. 105–112.
  * [27] S. M. Kay, _Fundamentals of Statistical Signal Processing, Volumn I: Estimation Theory_. New Jersey, USA: Prentice-Hall, 1993.
  * [28] J. Lee, G. Gil, and Y. H. Lee, “Channel estimation via orthogonal matching pursuit for hybrid MIMO systems in millimeter wave communications,” _IEEE Trans. Commun._ , vol. 64, no. 6, pp. 2370–2386, Jun. 2016.
  * [29] X. Rao and V. K. N. Lau, “Distributed compressive CSIT estimation and feedback for FDD multi-user massive MIMO systems,” _IEEE Trans. Signal Process._ , vol. 62, no. 12, pp. 3261–3271, Jun. 2014.
  * [30] Z. Qin, J. Fan, Y. Liu, Y. Gao, and G. Y. Li, “Sparse representation for wireless communications: A compressive sensing approach,” _IEEE Signal Process. Mag._ , vol. 35, no. 3, pp. 40–58, May 2018.
  * [31] S. Jaeckel, L. Raschkowski, K. Borner, and L. Thiele, “QuaDRiGa: A 3-D multi-cell channel model with time evolution for enabling virtual field trials,” _IEEE Trans. Antennas Propag._ , Sept. 2014.
  * [32] S. Jaeckel, L. Raschkowski, K. Borner, L. T. F. Burkhardt, and E. Eberlein, “OuaDRiGa-quasi deterministic radio channel generator, user manual and documentation,” _Fraunhofer Heinrich Hertz Institute, Tech. Rep., v2.6.1_ , 2021.
  * [33] 3GPP TR 38.901, “Study on channel model for frequencies from 0.5 to 100 GHz,” _v 16.1.0_ , Dec. 2019.