RESEARCH PAPER 2022

Near-Field Wideband Channel Estimation for Extremely Large-Scale MIMO

<EMAIL_ADDRESS>
Mingyao Cui and Linglong Dai

# Near-Field Wideband Channel Estimation for Extremely Large-Scale MIMO

Mingyao Cui  Linglong Dai  Department of Electronic Engineering, Tsinghua
University,
Beijing National Research Center for Information Science and Technology
(BNRist), Beijing 100084, China

###### Abstract

Extremely large-scale multiple-input-multiple-output (XL-MIMO) at millimeter-
wave (mmWave) and terahertz (THz) bands plays an important role in supporting
extreme high beamforming gain as well as ultra-wideband spectrum resources.
Unfortunately, accurate wideband XL-MIMO channel estimation suffers from the
new challenge called as the near-field beam split effect. Prior works either
neglect the accurate near-field channel model or fail to exploit the beam
split effect, resulting in poor channel estimation accuracy for wideband XL-
MIMO. To tackle this problem, this paper proposes a bilinear pattern detection
(BPD) based approach to accurately recover the wideband XL-MIMO channel.
Specifically, by analyzing the characteristics of near-field wideband
channels, we first reveal the bilinear pattern of the near-field beam split
effect, which implies that the sparse support set of near-field channels in
both the angle and the distance domains can be regarded as a linear function
against frequency. Then, inspired by the classical simultaneously orthogonal
matching pursuit technique, we use the bilinear pattern to estimate the angle-
of-arrival (AoA) and distance parameters of each near-field path component at
all frequencies. In this way, the entire wideband XL-MIMO channel can be
recovered by compressed sensing algorithms. Moreover, we provide the
computational complexity of the proposed algorithm compared with existing
algorithms. Finally, simulation results demonstrate that our scheme can
achieve the accurate estimation of the near-field wideband XL-MIMO channel in
the presence of near-field beam split effect.

###### keywords:

XL-MIMO, wideband, near-field, beam split, channel estimation

##  1 Introduction

Compared with the existing massive multiple-input-multiple-output (MIMO) with
dozens or hundreds of antennas for 5G, extremely large-scale MIMO (XL-MIMO)
embraces a 10-fold increase in the scale of antenna array, which enables the
considerable improvement of the spatial multiplexing gain [1, 2, 3, 4].
Moreover, benefiting from the high beamforming gain achieved by extremely
narrow pencil-like beams, XL-MIMO is essential to alleviate the serious path
loss of millimeter wave (mmWave) and terahertz (THz) signals with abundant
spectrum resources [5, 6, 7]. Consequently, the natural integration of high-
frequency wideband communications and XL-MIMO is regarded as a promising
technology in future 6G networks.

To unleash the performance superiority of wideband XL-MIMO, it is essential to
acquire accurate channel state information (CSI). Unfortunately, due to the
enormous number of antennas, XL-MIMO channel estimation requires an
unacceptable pilot overhead, particularly when the base station (BS) is
deploying the hybrid analog and digital precoding architecture with only a
small number of radio-frequency (RF) chains [8, 9]. Thereby, how to exploit
the channel structure to reduce the pilot overhead for channel estimation has
always been a hot topic from 5G massive MIMO to 6G XL-MIMO.

###  1.1 Prior Works

There are extensive works to explore the sparse structure of mmWave and THz
channels to design channel estimation algorithms with low pilot overhead.
These works can be generally categorized into three types: far-field
narrowband [10, 11, 12, 13, 14], far-field wideband [15, 16, 17], and near-
field narrowband [18, 19, 20] channel estimations.

The first type of _far-field narrowband_ channel estimation schemes mainly
works for 5G massive MIMO, and they usually assume that the bandwidth is
limited, e.g., around several hundreds of MHz. The array aperture of massive
MIMO is not large enough, resulting in the dominance of far-field transmission
environments, where a channel path can be modeled by the planar wave model. In
this case, the multipath far-field channel is built up by the superposition of
several planar waves [9]. As each planar wavefront is dependent on a certain
angle-of-arrival (AoA), and the number of channel paths is usually small, this
far-field channel exhibits sparsity in the angle domain. Thus it can be
recovered by compressed sensing (CS) based algorithms with low pilot overhead.
For instance, [10] utilized orthogonal matching pursuit (OMP) scheme to
estimate the angle-domain sparse channel at a single frequency. Moreover,
several works have considered to improve the OMP algorithm by using the
channel sparse structure at multiple frequencies. Specifically, due to the
limited bandwidth, the common sparse support assumption has been widely used
to transform the multi-frequency channel estimation into a multiple
measurement vector (MMV) problem, where the sparse support sets are the same
at different frequencies. Then, various MMV-specific algorithms can be
utilized to estimate the far-field channel, such as the simultaneously OMP
(SOMP) [12, 13] and the group-sparse Bayesian CS approach [14].

Compared with the _far-field narrowband_ schemes above, the _far-field
wideband_ channel estimation algorithms assume a pretty wide bandwidth. In
this condition, although communications are still operating in far-field
environments, the widely-used common sparse support assumption [12, 13] does
not hold anymore due to the beam split effect [15, 16]. The beam split effect
indicates that owing to the wide bandwidth, the spatial channel directions are
separated from each other at different frequencies, which results in a
frequency-dependent angle-domain sparse structure and severely degrades the
performance of classical MMV-specific solutions. Recently, a few advanced
signal processing techniques have been proposed to cope with this beam split
effect in far-field scenarios. Specifically, [16] proved the beam split
pattern in the angle domain, i.e., the sparse support set of a physical angle
grows linearly with the subcarrier frequency. Then, a beam split pattern
detection (BSPD) based algorithm was proposed in [16] to construct the match
between a certain physical angle and the beam split pattern for recovering the
far-field wideband channel. Besides, in [17], the authors proposed to
construct a series of dictionary matrices for each subcarrier frequency to
match the beam split pattern, based on which a generalized simultaneous OMP
(GSOMP) was applied to overcome the beam split effect in the far field.

For the third category of the _near-field narrowband_ algorithms, they place
heavy emphasis on the fundamental change from 5G massive MIMO to 6G XL-MIMO by
considering the near-field propagation. The deployment of XL-MIMO makes the
near-field propagation highly possible to happen, particularly for high-
frequency communications [21]. To be specific, the radius of near-field areas,
which is also called as Fraunhofer distance, is proportional to the square of
array aperture and frequency [22]. The Frauhofer distance may reach several
dozens or hundreds of meters in high-frequency XL-MIMO systems. As opposed to
the far-field channel model based on the planar wavefront assumption, the
multi-path near-field channel model is built up by the superposition of
multiple _spherical waves_ [23]. Although the number of channel paths is still
limited, the angle-domain near-field channel is not sparse any more, since
each spherical wave is dependent on a range of physical angles [20, 18]. Thus,
existing angle-domain sparsity based techniques cannot work well in near-field
environments. To address this problem, [20] proposed to transform the near-
field channel into the joint angle-distance domain (polar domain). As an
angle-distance sample contains enough information to describe a spherical
wave, the near-field channel can be then sparsely represented in the polar
domain. Accordingly, the polar-domain SOMP algorithms were used in [20] to
recover the near-field channel accurately.

Despite the extensive studies on XL-MIMO channel estimation, the above three
categories of techniques neglect the crucial topic of _near-field wideband_
channel estimation. In other words, they overlook the important integration of
the near-field and beam split effect. For wideband XL-MIMO communications,
both the array aperture and bandwidth are pretty large. Such a wideband XL-
MIMO system poses a _near-field beam split_ effect [24], where the sparse
channel support set in the polar domain differs from each other at different
frequencies. Furthermore, as all prior techniques (far-field narrowband, far-
field wideband, and near-field narrowband approaches) mismatch this polar-
domain frequency-dependent sparse structure, they will suffer from a serious
performance loss in wideband XL-MIMO systems. To the best of our knowledge,
the essential near-field wideband channel estimation has not been studied in
the literature.

###  1.2 Our Contributions

To fill in this gap, we propose a bilinear pattern detection (BPD) based
algorithm to realize near-field wideband channel estimation by exploring the
unique polar-domain sparse structure with near-field beam split. Our
contributions are summarized as follows.

  * •

First, we comprehensively analyze the polar-domain channel sparse structure
over different subcarrier frequencies when the near-field beam split effect is
considered. By analyzing the mutual coherence of near-field array response
vectors at different frequencies, we reveal the phenomenon that both the
angle-domain and distance-domain sparse support sets grow linearly with the
subcarrier frequency, which is defined as the bilinear pattern (BP) of near-
field beam split effect. This phenomenon indicates that the physical location
of a scatter (or user) in a channel path has a one-to-one matching
relationship with a certain pattern of the frequency-dependent sparse support
sets in the polar domain. This relationship can be used to accurately estimate
the physical location of each path.

  * •

Then, a BPD based channel estimation algorithm is proposed according to the
revealed phenomenon above. We first construct the polar-domain bilinear
patterns for all angle-distance samples. Then, for each path component, we use
these bilinear patterns to accumulate the power of each angle-distance sample
from the entire bandwidth and capture the sample with the largest power as the
estimated physical location. We carry out the above procedure several times to
successively detect all channel paths. As the proposed algorithm can jointly
exploit the polar-domain sparsity and the bilinear pattern, it is promising to
achieve accurate channel estimation.

  * •

Finally, we analyze the computational complexity of the BPD based algorithm
and compare it with existing algorithms. Extensive simulation results are
provided to demonstrate the superior performance of the proposed scheme. It is
illustrated that the BPD based channel estimation technique outperforms
existing schemes. More importantly, it is verified that the BPD based
algorithm is capable of realizing accurate channel estimation in all far-
field/near-field/narrowband/wideband conditions111Simulation codes of this
paper are available at
http://oa.ee.tsinghua.edu.cn/dailinglong/publications/publications.html..

_Organization_ : The remainder of this paper is organized as follows. Section
2 presents the system model, including the channel model and problem
formulation. In Section 3, the bilinear pattern is proved, and the BPD based
algorithm is provided. Simulations are carried out in Section 4, and finally
conclusions are drawn in Section 5.

##  2 System Model

###  2.1 Channel Model

In this paper, we investigate the uplink channel estimation for XL-MIMO system
with orthogonal frequency division multiplexing (OFDM), where one base station
(BS) uses an $N$-antenna uniform linear array to serve $K$ single-antenna
users. We adopt the Saleh-Valenzuela multipath channel to model the near-field
channel [25, 20]. As presented in Fig. 1, the frequency-domain channel
$\mathbf{h}_{m}\in\mathbb{C}^{N\times 1}$ at the $m$-th subcarrier
($m\in\\{0,1,\cdots,M-1\\}$) from a certain user is denoted as

$\displaystyle\mathbf{h}_{m}=\sqrt{\frac{N}{L}}\sum_{l=0}^{L-1}g_{l,m}e^{-j\frac{2\pi}{\lambda_{m}}r_{l}}\mathbf{a}(\vartheta_{l},r_{l},f_{m}).$
(1)

Here, $g_{l,m}$ denotes the complex path gain at $m$-th subcarrier of the
$l$-th path, and $L$ represents the number of channel paths. Besides,
$f_{m}=f_{c}+\frac{2m-M}{2M}B$ is the $m$-th subcarrier frequency with $B$ and
$f_{c}$ corresponding to the bandwidth and the carrier frequency, and
$\lambda_{m}=\frac{c}{f_{m}}$ is the wavelength at $f_{m}$ with $c$ being the
light speed. Moreover, as shown in Fig. 1, $\vartheta_{l}$ represents the AOA
of the $l$-th path, $r_{l}$ is the distance from the last-hop scatter to the
BS array’s center, and $\mathbf{a}(\vartheta_{l},r_{l},f_{m})$ is the array
response vector of $(\vartheta_{l},r_{l})$ describing the spherical wavefront
between the last-hop scatter and BS. According to the geometrical structure of
BS array, $\mathbf{a}(\vartheta_{l},r_{l},f_{m})$ is modeled as

$\displaystyle\mathbf{a}(\vartheta_{l},r_{l},f_{m})=\frac{1}{\sqrt{N}}[e^{j\phi_{l,m}^{(0)}},e^{j\phi_{l,m}^{(1)}},\cdots,e^{j\phi_{l,m}^{(N-1)}}]^{T},$
(2)

where $\phi_{l,m}^{(n)}=-\frac{2\pi}{\lambda_{m}}(r_{l}^{(n)}-r_{l})$ and
$r_{l}^{(n)}$ is the distance from the last-hop scatter to the $n$-th BS
antenna. As indicated in [20], $r_{l}^{(n)}$ can be written as
$r_{l}^{(n)}=\sqrt{r_{l}^{2}+\delta_{n}^{2}d^{2}-2r_{l}\delta_{n}d\sin\vartheta_{l}}$,
where $\delta_{n}=n-\frac{N-1}{2}$ with $n\in\\{0,1,\cdots,N-1\\}$. In
addition, $d=\frac{\lambda_{c}}{2}$ is the antena spacing and
$\lambda_{c}=\frac{c}{f_{c}}$ is the carrier wavelengh.

Figure 1:  The near-field channel model for an arbitrary user.

From (2), two crucial aspects are revealed: the near-field propagation and the
beam split effect. First, the far-field scenario is always adopted in
conventional massive MIMO to simplify $\phi_{l,m}^{(n)}$ as
$\frac{2\pi}{\lambda_{m}}\delta_{n}d\sin\vartheta_{l}$, which is independent
of distance $r_{l}$. On the contrary, in XL-MIMO systems, the Fraunhofer
distance $\frac{2D^{2}}{\lambda_{c}}$ is comparable to $r_{l}$, making it
highly possible for signals to propagate in near-field environments [22].
Here, $D=Nd$ is the array aperture. For example, we consider a 0.3-meter-array
working at 100 GHz, the Fraunhofer distance of which reaches 60 meters. In
this case, the spherical wavefront model has to be investigated and the impact
of $r_{l}$ cannot be ignored [21]. Second, in narrowband communications when
$f_{m}\approx f_{c}$, the array response vector is nearly frequency-
independent, leading to the common sparse structure. This structure allows
jointly estimating channels at different subcarriers by using MMV-specific CS
algorithms, such as SOMP and the group-sparse Bayesian approach [12, 13, 14,
20]. However, in wideband communications when $f_{m}\neq f_{c}$, the array
response vector considerably varies over frequencies. More seriously, the
support set of different subcarriers can be quite different from each other in
wideband XL-MIMO systems, undermining the common sparse support structure of
the MMV model. This phenomenon is defined as the beam split effect [16, 24],
which severely degrades the performance of conventional channel estimation
schemes.

To the best of our knowledge, existing channel estimation works fail to
jointly exploit these two underlying characteristics of the wideband XL-MIMO
channel, giving rise to poor estimation performance. To fill in this gap, this
paper simultaneously considers the near-field propagation and beam split
effect (or the near-field beam split effect) to realize high-accuracy channel
estimation..

###  2.2 Problem Formulation

As indicated in Fig. 2, we adopt an uplink time-division-duplex (TDD) channel
estimation scenario, where $K$ users harness orthogonal time or frequency
resources to transmit orthogonal pilots to BS. Thus we can consider channel
estimation for an arbitrary user without loss of generality. We suppose the
length of pilot sequence is $P$. Denote $s_{m,p}$ as the pilot signal at the
$m$-th subcarrier in the $p$-th time slot. Then, the received signals
$\mathbf{y}_{m,p}\in\mathbb{C}^{N_{\text{RF}}\times 1}$ can be presented as

$\displaystyle\mathbf{y}_{m,p}=\mathbf{A}_{p}\mathbf{h}_{m}s_{m,p}+\mathbf{A}_{p}\mathbf{n}_{m,p},$
(3)

where $\mathbf{A}_{p}\in\mathbb{C}^{N_{\text{RF}}\times N}$ denotes the
combining matrix and $\mathbf{n}_{m,p}\in\mathbb{C}^{N\times 1}$ denotes the
complex Gaussian noise following the distribution
$\mathcal{CN}(0,\sigma^{2}\mathbf{I}_{N})$. Assume $s_{m,p}=1$ for
$p=1,2,\cdots,P$. Then, considering the total $P$ time slots for pilot
trasmission, we arrive at

$\displaystyle\mathbf{y}_{m}=\mathbf{A}\mathbf{h}_{m}+\mathbf{n}_{m},$ (4)

where
$\mathbf{y}_{m}=[\mathbf{y}_{m,1}^{T},\mathbf{y}_{m,2}^{T},\cdots,\mathbf{y}_{m,P}^{T}]^{T}\in\mathbb{C}^{PN_{\text{RF}}\times
1}$ and
$\mathbf{n}_{m}=[\mathbf{n}_{m,1}^{T}\mathbf{A}_{1}^{T},\mathbf{n}_{m,2}^{T}\mathbf{A}_{2}^{T},\cdots,\mathbf{n}_{m,P}^{T}\mathbf{A}_{P}^{T}]^{T}\in\mathbb{C}^{PN_{\text{RF}}\times
1}$. Besides,
$\mathbf{A}=[\mathbf{A}_{1}^{T},\mathbf{A}_{2}^{T},\cdots,\mathbf{A}_{P}^{T}]^{T}\in\mathbb{C}^{PN_{\text{RF}}\times
N}$ is the overall observation matrix, whose elements can be independently
generated from the binomial distribution $\frac{1}{\sqrt{N}}\\{-1,1\\}$.

Figure 2:  An illustration of the multi-user uplink transmission system.

To recover the high-dimensional near-field channel $\mathbf{h}_{m}$ from the
low-dimensional observation $\mathbf{y}_{m}$, where $PN_{\text{RF}}<N$, we can
transform the antenna-domain channel to its angle-distance domain (polar
domain) by a polar-domain representation matrix
$\mathbf{W}\in\mathbb{C}^{N\times N_{a}N_{d}}$ [20]. As discussed in [20],
$\mathbf{W}$ is composed of $N_{d}$ submatrices, where

$\displaystyle\mathbf{W}=[\mathbf{W}_{1},\mathbf{W}_{2},\cdots,\mathbf{W}_{N_{d}}].$
(5)

Each submatrix $\mathbf{W}_{n_{d}}\in\mathbb{C}^{N\times N_{a}}$ includes
$N_{a}$ array response vectors at the carrier frequency, where

$\displaystyle\mathbf{W}_{n_{d}}=[\mathbf{a}(\bar{\vartheta}_{1},\bar{r}_{n_{d},1},f_{c}),\mathbf{a}(\bar{\vartheta}_{2},\bar{r}_{n_{d},2},f_{c}),\cdots,\mathbf{a}(\bar{\vartheta}_{n_{a}},\bar{r}_{n_{d},N_{a}},f_{c})].$
(6)

As discussed in [20], the sampled angles and distances, i.e.,
($\bar{\vartheta}_{n_{a}},\bar{r}_{n_{d},n_{a}}$) with
$n_{a}\in[1,2,\cdots,N_{a}]$ and $n_{d}\in[1,2,\cdots,N_{d}]$, should satisfy
$\sin\bar{\vartheta}_{n_{a}}=\frac{2(n_{a}-1)-N_{a}}{N_{a}}$ and
$\bar{r}_{n_{d},n_{a}}=\frac{D^{2}\cos^{2}\bar{\vartheta}_{n_{a}}}{2\beta^{2}\lambda_{c}n_{d}}$,
where $\beta$ is a predefined parameter. Then, the polar-domain channel
$\bar{\mathbf{h}}_{m}$ can be obtained by solving the underdetermined equation
$\mathbf{h}_{m}=\mathbf{W}\bar{\mathbf{h}}_{m}$. As investigated in [20],
matrix $\mathbf{W}$ is capable of extracting the angle-distance information
embedded in each array response vector of (1). In addition, the number of
paths $L$ is usually small, especially for mmWave and THz bands (e.g. $L=5$).
As a consequence, the polar-domain channel $\bar{\mathbf{h}}_{m}$ is sparse
and the near-field channel estimation problem is equivalent to a sparse signal
recovery problem:

$\displaystyle\mathbf{y}_{m}=\mathbf{A}\mathbf{W}\bar{\mathbf{h}}_{m}+\mathbf{n}_{m}.$
(7)

Finally, we attempt to jointly estimate channels for all subcarriers. Hence,
we rewrite (7) as

$\displaystyle\mathbf{Y}=\mathbf{A}\mathbf{W}\bar{\mathbf{H}}+\mathbf{N},$ (8)

where
$\mathbf{Y}=[\mathbf{y}_{1},\mathbf{y}_{2},\cdots,\mathbf{y}_{M}]\in\mathbb{C}^{PN_{\text{RF}}\times
M}$,
$\bar{\mathbf{H}}=[\bar{\mathbf{h}}_{1},\bar{\mathbf{h}}_{2},\cdots,\bar{\mathbf{h}}_{M}]$,
and $\mathbf{N}=[\mathbf{n}_{1},\mathbf{n}_{2},\cdots,\mathbf{n}_{M}]$. Notice
that as the array response vector (2) is frequency-dependent, the sparse
support set of $\bar{\mathbf{h}}_{m}$ for different subcarriers $f_{m}$ is
quite different from each other. Therefore, existing wideband channel
estimation algorithms tailored for MMV problems cannot perfectly match the
model (8). Consequently, an accurate near-field channel estimation scheme for
wideband XL-MIMO is essential.

##  3 Bilinear Pattern Detection Based Channel Estimation

In this section, we first reveal the bilinear pattern of near-field beam split
effect. This pattern reveals that with proper preprocessing, the support set
of near-field channels varies linearly over frequencies in both the angle and
the distance domain. Then, we propose a bilinear pattern detection (BSD) based
near-field wideband channel estimation scheme. Finally, we provide complexity
analysis of our scheme.

###  3.1 Bilinear Pattern of Near-Field Beam Split

To jointly recover the wideband channel for all subcarriers in (8), it is
crucial to investigate the map between the support index and frequency for an
arbitrary path, which is defined as the beam split pattern (BSP). This pattern
has been extensively studied under _far-field_ environments in the literature
[15, 16, 17], while we are trying to discover the pattern of _near-field beam
split_.

Specifically, we investigate the $l$-th scatter located at
$(\vartheta_{l},r_{l})$, which corresponds to a series of array response
vectors $\mathbf{a}(\vartheta_{l},r_{l},f_{m})$ with different frequencies.
According to the CS theory, the support index
$(n_{a,l,m}^{\star},n_{d,l,m}^{\star})$ of
$\mathbf{a}(\vartheta_{l},r_{l},f_{m})$ refers to the index of the strongest
element of $\mathbf{W}^{H}\mathbf{a}(\vartheta_{l},r_{l},f_{m})$, i.e.,
$(n_{a,l,m}^{\star},n_{d,l,m}^{\star})=\arg\max_{n_{a},n_{d}}\|\mathbf{a}^{H}(\bar{\vartheta}_{n_{a}},\bar{r}_{n_{d},n_{a}},f_{c})\mathbf{a}(\vartheta_{l},r_{l},f_{m})\|$.
Unfortunately, as the phase
$\phi_{l,m}^{(n)}=-\frac{2\pi}{\lambda_{m}}(r_{l}^{(n)}-r_{l})$ of the $n$-th
element of $\mathbf{a}(\vartheta_{l},r_{l},f_{m})$ is a complicated radical
function, it is intractable to obtain the close form of
$(n_{a,l,m}^{\star},n_{d,l,m}^{\star})$.

To address this problem, we use the second-order Taylor expansion to
approximate $\phi_{l,m}^{(n)}$ as

$\displaystyle\phi_{l,m}^{(n)}\approx\bar{\phi}_{l,m}^{(n)}=\frac{2\pi}{\lambda_{m}}(\delta_{n}d\sin\vartheta_{l}-\delta_{n}^{2}d^{2}\frac{\cos^{2}\vartheta_{l}}{2r_{l}}),$
(9)

which is much more accurate than the conventional far-field approximation with
first-order Taylor expansion. For expression simplicity, we use parameters
$(\theta_{l},\alpha_{l})=(\sin\vartheta_{l},\frac{\cos^{2}\vartheta_{l}}{2r_{l}})$
to replace parameters $(\vartheta_{l},r_{l})$, so we have
$\bar{\phi}_{l,m}^{(n)}=\frac{2\pi}{\lambda_{m}}(\delta_{n}d\theta_{l}-\delta_{n}^{2}d^{2}\alpha_{l})$.
As a result, we can use
$\bar{\mathbf{a}}(\theta_{l},\alpha_{l},f_{m})={\frac{1}{\sqrt{N}}}[e^{j\bar{\phi}_{l,m}^{(0)}},\cdots,e^{j\bar{\phi}_{l,m}^{(N-1)}}]^{T}$
to simplify $\mathbf{a}(\vartheta_{l},r_{l},f_{m})$. Similarly, vector
$\mathbf{a}(\bar{\vartheta}_{n_{a}},\bar{r}_{n_{d},n_{a}},f_{c})$ can be
approximated as
$\bar{\mathbf{a}}(\bar{\theta}_{n_{a}},\bar{\alpha}_{n_{d}},f_{c})$. Notice
that the sampled angles $\bar{\vartheta}_{n_{a}}$ and distances
$\bar{r}_{n_{d},n_{a}}$ are transformed to
$\bar{\theta}_{n_{a}}=\sin\bar{\vartheta}_{n_{a}}=\frac{2(n_{a}-1)-N_{a}}{N_{a}}$
and
$\alpha_{n_{b}}=\frac{\cos^{2}\bar{\vartheta}_{n_{a}}}{2\bar{r}_{n_{d},n_{a}}}=\frac{\cos^{2}\bar{\vartheta}_{n_{a}}}{2\frac{D^{2}\cos^{2}\bar{\vartheta}_{n_{a}}}{2\beta^{2}\lambda_{c}n_{d}}}=\frac{\beta^{2}\lambda_{c}n_{d}}{D^{2}}$,
with $n_{a}\in\\{1,2,\cdots,N_{a}\\}$ and $n_{d}\in\\{1,2,\cdots,N_{d}\\}$.

Based on the above discussion, the near-field beam split pattern can be
acquired by solving the problem
$(n_{a,l,m}^{\star},n_{d,l,m}^{\star})=\arg\max_{n_{a},n_{d}}\|\bar{\mathbf{a}}^{H}(\bar{\theta}_{n_{a}},\bar{\alpha}_{n_{d}},f_{c})\bar{\mathbf{a}}(\theta_{l},\alpha_{l},f_{m})\|$.
The following Lemma 1 provides a more concise form of
$\|\bar{\mathbf{a}}^{H}(\bar{\theta}_{n_{a}},\bar{\alpha}_{n_{d}},f_{c})\bar{\mathbf{a}}(\theta_{l},\alpha_{l},f_{m})\|$.
We consider two frequencies $f_{1}$ and $f_{2}$, and two arbitrary locations
$(\vartheta_{1},r_{1})$ and $(\vartheta_{2},r_{2})$ with
$(\theta_{1},\alpha_{1})=(\sin\vartheta_{1},\frac{\cos^{2}\vartheta_{1}}{2r_{1}})$
and
$(\theta_{2},\alpha_{2})=(\sin\vartheta_{2},\frac{\cos^{2}\vartheta_{2}}{2r_{2}})$.
When the number of antennas $N\rightarrow+\infty$, then
$\|\bar{\mathbf{a}}^{H}(\theta_{1},\alpha_{1},f_{1})\bar{\mathbf{a}}(\theta_{2},\alpha_{2},f_{2})\|$
can be represented by

$\displaystyle\|\bar{\mathbf{a}}^{H}(\theta_{1},\alpha_{1},f_{1})\bar{\mathbf{a}}(\theta_{2},\alpha_{2},f_{2})\|\overset{N\rightarrow+\infty}{=}\left|\int_{-\frac{1}{2}}^{\frac{1}{2}}e^{j2\pi
x\gamma-j2\pi x^{2}\zeta}{\rm d}x\right|=\Xi(\gamma,\zeta),$ (10)

where
$\gamma=\frac{D}{\lambda_{1}}(\theta_{1}-\frac{f_{2}}{f_{1}}\theta_{2})$,
$\zeta=\frac{D^{2}}{\lambda_{1}}(\alpha_{1}-\frac{f_{2}}{f_{1}}\alpha_{2})$,
and $\lambda_{1}=\frac{c}{f_{1}}$.

Figure 3:  The numerical results of function $\Xi(\gamma,\zeta)$.

_Proof_ : Accoding to the definition of
$\bar{\mathbf{a}}(\theta_{1},\alpha_{1},f_{1})$ and
$\bar{\mathbf{a}}(\theta_{2},\alpha_{2},f_{2})$, we have

$\displaystyle\|\bar{\mathbf{a}}^{H}(\theta_{1},\alpha_{1},f_{1})\bar{\mathbf{a}}(\theta_{2},\alpha_{2},f_{2})\|$
$\displaystyle=\left|\frac{1}{N}\sum_{n=-(N-1)/2}^{(N-1)/2}e^{j\frac{2\pi
nd}{\lambda_{1}}(\theta_{1}-\frac{f_{2}}{f_{1}}\theta_{2})-j\frac{2\pi
n^{2}d^{2}}{\lambda_{1}}(\alpha_{1}-\frac{f_{2}}{f_{1}}\alpha_{2})}\right|$
$\displaystyle\overset{(a)}{=}\left|\sum_{x=-\frac{1}{2}+\frac{1}{2N}}^{\frac{1}{2}-\frac{1}{2N}}e^{j\frac{2\pi
Dx}{\lambda_{1}}(\theta_{1}-\frac{f_{2}}{f_{1}}\theta_{2})-j\frac{2\pi
D^{2}x^{2}}{\lambda_{1}}(\alpha_{1}-\frac{f_{2}}{f_{1}}\alpha_{2})}\frac{1}{N}\right|,$
(11)

where (a) is achieved by letting $x=\frac{n}{N}$ and $D=Nd$. In addition,
since XL-MIMO has a very large number of antennas $N$, the summation (11) can
be rewriten as the following Riemann integral

$\displaystyle\|\bar{\mathbf{a}}^{H}(\theta_{1},\alpha_{1},f_{1})\bar{\mathbf{a}}(\theta_{2},\alpha_{2},f_{2})\|\overset{N\rightarrow+\infty}{=}\left|\int_{-\frac{1}{2}}^{\frac{1}{2}}e^{j\frac{2\pi
Dx}{\lambda_{1}}(\theta_{1}-\frac{f_{2}}{f_{1}}\theta_{2})-j\frac{2\pi
D^{2}x^{2}}{\lambda_{1}}(\alpha_{1}-\frac{f_{2}}{f_{1}}\alpha_{2})}\text{d}x\right|.$
(12)

Eventually, plug
$\gamma=\frac{D}{\lambda_{1}}(\theta_{1}-\frac{f_{2}}{f_{1}}\theta_{2})$ and
$\zeta=\frac{D^{2}}{\lambda_{c}}(\alpha_{1}-\frac{f_{2}}{f_{1}}\alpha_{2})$
into (12), the right-hand side in (10) can be acquired, and the proof is
completed. $\hfill\blacksquare$

Lemma 1 indicates that for XL-MIMO communications with a very large number of
antennas $N$, finding the support index
$(n_{a,l,m}^{\star},n_{d,l,m}^{\star})$ is equivalent to solving the following
problem:

$\displaystyle(n_{a,l,m}^{\star},n_{d,l,m}^{\star})=\arg\max_{n_{a},n_{d}}\Xi\left(\frac{D}{\lambda_{c}}(\bar{\theta}_{n_{a}}-\frac{f_{m}}{f_{c}}\theta_{l}),\frac{D^{2}}{\lambda_{c}}(\bar{\alpha}_{n_{d}}-\frac{f_{m}}{f_{c}}\alpha_{l})\right).$
(13)

To address problem (13), it is of great importance to grasp the property of
function $\Xi(\gamma,\zeta)$. First of all, it is easy to prove that
$\Xi(\gamma,\zeta)$ is an even function with respect to (w.r.t) $\gamma$ and
$\zeta$, i.e.,

$\displaystyle\Xi(\gamma,\zeta)=\Xi(-\gamma,\zeta),\quad\Xi(\gamma,\zeta)=\Xi(\gamma,-\zeta).$
(14)

Therefore, we can rewrite $\Xi(\gamma,\zeta)$ as
$\Xi(\gamma,\zeta)=\Xi(|\gamma|,|\zeta|)$, or we only need to discuss the
domain $\gamma>0$ and $\zeta>0$.

In addition, $\Xi(\gamma,\zeta)$ is a non-parameter function, we can use
numerical integration to acquire its numerical result, which is plotted in
Fig. 3. This numerical result indicates that within the main lobe of
$\Xi(\gamma,\zeta)$ (area surrounded by the black line in Fig. 3 (b)), the
value of $\Xi(\gamma,\zeta)$ declines with the increase of $|\gamma|$ and
$|\zeta|$. That is to say, if $|\gamma_{1}|\leq|\gamma_{2}|$ and
$|\zeta_{1}|\leq|\zeta_{2}|$, then we have
$\Xi(\gamma_{1},\zeta_{1})\geq\Xi(\gamma_{2},\zeta_{2})$. Making use of these
properties, Lemma 2 gives out the close form solution of (13).

If
$\frac{D}{\lambda_{c}}(\theta_{n_{a,l,m}^{\star}}-\frac{f_{m}}{f_{c}}\theta_{l})$
and
$\frac{D^{2}}{\lambda_{c}}(\alpha_{n_{a,l,m}^{\star}}-\frac{f_{m}}{f_{c}}\alpha_{l})$
are within the main lobe of $\Xi(\gamma,\zeta)$, then the optimal solution of
problem (13) is

$\displaystyle
n_{a,l,m}^{\star}=\arg\min_{{n_{a}}}|\bar{\theta}_{n_{a}}-\frac{f_{m}}{f_{c}}\theta_{l}|,$
(15) $\displaystyle
n_{d,l,m}^{\star}=\arg\min_{{n_{d}}}|\bar{\alpha}_{n_{d}}-\frac{f_{m}}{f_{c}}\alpha_{l}|.$
(16)

_Proof_ : In the main lobe, since $\Xi(|\gamma|,|\zeta|)$ is a decreasing
function w.r.t $|\gamma|$ and $|\zeta|$, and
$\Xi(\gamma,\zeta)=\Xi(|\gamma|,|\zeta|)$, hence maximizing
$\Xi\left(\frac{D}{\lambda_{c}}(\bar{\theta}_{n_{a}}-\frac{f_{m}}{f_{c}}\theta_{l}),\frac{D^{2}}{\lambda_{c}}(\bar{\alpha}_{n_{d}}-\frac{f_{m}}{f_{c}}\alpha_{l})\right)$
is equivalent to minimizing
$|\frac{D}{\lambda_{c}}(\bar{\theta}_{n_{a}}-\frac{f_{m}}{f_{c}}\theta_{l})|$
and
$|\frac{D^{2}}{\lambda_{c}}(\bar{\alpha}_{n_{d}}-\frac{f_{m}}{f_{c}}\alpha_{l})|$
separately, which gives rise to the solutions (15) and (16).
$\hfill\blacksquare$

Lemma 2 explicitly displays the bilinear pattern over frequencies of the near-
field beam split effect. As illustrated by the dotted lines in Fig. 4, both
$\frac{f_{m}}{f_{c}}\alpha$ and $\frac{f_{m}}{f_{c}}\theta$ grow linearly to
frequency $f_{m}$. The support indices $n_{a,l,m}^{\star}$ and
$n_{d,l,m}^{\star}$ correspond to the sampled $\theta_{n_{a}}$ and
$\alpha_{n_{d}}$ that are nearest to $\frac{f_{m}}{f_{c}}\theta$ and
$\frac{f_{m}}{f_{c}}\alpha$. Thus, this double linear growth pattern is
defined as the bilinear pattern of near-field beam split. More importantly,
the bilinear pattern can be utilized to capture the support set of the $l$-th
path for all frequencies, so as to accurately estimate the location
$(\vartheta_{l},r_{l})$. Its superiority comes from two aspects. At first,
Lemma 2 implies that the indices $n_{a,l,m}^{\star}$ and $n_{d,l,m}^{\star}$
can gather the largest power of the $l$-th path from the entire bandwidth.
Therefore, compared to conventional common support set assumption in MMV
model, the proposed bilinear pattern enjoys higher recovery accuracy for
wideband channels. In addition, different from existing works that only
explain the beam split pattern (15) in the angle domain [15, 16, 17], the
bilinear pattern reveals the linear-changing structure in the joint angle-
distance domain. Thus, the proposed pattern is capable of estimating both far-
field and near-field channels.

###  3.2 Bilinear Pattern Detection Based Channel Estimation

Figure 4:  The bilinear pattern of the near-field beam split effect.

In this subsection, based on the bilinear pattern discussed before, a bilinear
pattern detection based channel estimation scheme is proposed to solve problem
$(\ref{eq:MMV})$.  This algorithm is inspired by the classical CS-based polar-
domain SOMP method [20]. As illustrated in [20], the polar-domain SOMP method
assumes that the polar-domain channel $\bar{{\mathbf{h}}}_{m}$ is sparse and
the sparse support sets of $\bar{{\mathbf{h}}}_{m},\forall m$ are the same.
Based on this assumption, the polar-domain SOMP method computes the total
power of each row of $\bar{{\mathbf{H}}}$ so as to accumulate the information
from the entire bandwidth. In this way, it can recover the location
$(\vartheta_{l},r_{l})$. However, we have proved previously that the actual
sparse support sets for different frequencies are not the same in wideband
conditions. Indeed, they should satisfy the bilinear pattern presented in
(15)-(16). Therefore, instead of computing the total power of each row of
$\bar{{\mathbf{H}}}$, our BPD based algorithm takes advantage of the bilinear
pattern to collect powers from all frequencies, so as to improve the accuracy
of the estimation for $(\vartheta_{l},r_{l})$.

Algorithm 1 Bilinear pattern detection based channel estimation algorithm

0: Received signal matrix $\mathbf{Y}$; Observation matrix $\mathbf{A}$;
Polar-domain representation matrix $\mathbf{W}$; Number of paths to be
detected $\hat{L}$;

0: Estimated wideband channel $\hat{\mathbf{H}}$;Pre-whitening stage

1:
$\mathbf{C}={\rm{diag}}\\{\mathbf{A}_{1}\mathbf{A}_{1}^{H},\mathbf{A}_{2}\mathbf{A}_{2}^{H},\cdots,\mathbf{A}_{P}\mathbf{A}_{P}^{H}\\}$;

2: $\mathbf{C}=\sigma^{2}{\mathbf{S}}\bm{\Sigma}{\mathbf{S}}^{H}$,
${\mathbf{D}}={\mathbf{S}}\bm{\Sigma}^{\frac{1}{2}}$;

3: Use matrix $\mathbf{D}$ to pre-white the received signal as
$\bar{\mathbf{Y}}=\mathbf{D}^{-1}\mathbf{Y}=\bm{\Psi}\bar{\mathbf{H}}+\bar{\mathbf{N}}=[\bar{{\mathbf{y}}}_{1},\cdots,\bar{{\mathbf{y}}}_{M}]$;Path
detection stage

4: Initialize the residue matrix as
$\mathbf{R}=[\mathbf{r}_{1},\mathbf{r}_{2},\cdots,\mathbf{r}_{M}]=\bar{\mathbf{Y}}$
and the support set as $\Upsilon=\\{\emptyset\\}$;

5: Angle-domain linear pattern:
$\Gamma({n_{a}},m)=\arg$$\min_{n}|\bar{\theta}_{n}-\frac{f_{m}}{f_{c}}\bar{\theta}_{n_{a}}|$
for $n_{a}\in\\{1,\cdots,N_{a}\\}$ and $m\in\\{1,\cdots,M\\}$;

6: Distance-domain linear pattern:
$\Lambda({n_{d}},m)=\arg$$\min_{n}|\bar{\alpha}_{n}-\frac{f_{m}}{f_{c}}\bar{\alpha}_{n_{d}}|$
for $n_{d}\in\\{1,\cdots,N_{d}\\}$ and $m\in\\{1,\cdots,M\\}$;

7: for $l\in\\{1,2,\cdots,\hat{L}\\}$ do

8:
$\mathbf{U}=\bm{\Psi}^{H}\mathbf{R}=[\mathbf{u}_{1},\mathbf{u}_{2},\cdots,\mathbf{u}_{M}]$;

9:
$(n_{a,l}^{\star},n_{d,l}^{\star})=\arg\max_{(n_{a},n_{d})}\sum_{m=1}^{M}\|u_{m}\left((\Lambda(n_{d},m)-1)N_{a}+\Gamma({n_{a}},m)\right)\|^{2}$;

10: $\Upsilon=\Upsilon\cup\\{(n_{a,l}^{\star},n_{d,l}^{\star})\\}$;

11: for $m\in\\{1,\cdots,M\\}$ do

12:
$\Upsilon_{m}=\\{\left(\Lambda(n_{d},m)-1\right)N_{a}+\Gamma({n_{a}},m)|(n_{a},n_{d})\in\Upsilon\\}$;

13: $\hat{\bar{{\mathbf{h}}}}_{m}={\mathbf{0}}_{N_{a}N_{d}\times 1}$,
$\hat{\bar{{\mathbf{h}}}}_{m}(\Upsilon_{m})=\bm{\Psi}^{\dagger}(:,\Upsilon_{m})\bar{{\mathbf{y}}}_{m}$;

14: end for

15: Update the residue
${\mathbf{R}}=\bar{{\mathbf{Y}}}-\bm{\Psi}[\hat{\bar{{\mathbf{h}}}}_{1},\hat{\bar{{\mathbf{h}}}}_{2},\cdots,\hat{\bar{{\mathbf{h}}}}_{M}]$;

16: end for

17:
$\hat{{\mathbf{H}}}={\mathbf{W}}[\hat{\bar{{\mathbf{h}}}}_{1},\hat{\bar{{\mathbf{h}}}}_{2},\cdots,\hat{\bar{{\mathbf{h}}}}_{M}]$.

Specifically, the pseudo-code of the proposed BPD algorithm is provided in
Algorithm 1. The overall algorithm is composed of two stages. In the first
pre-whitening stage, as the additional noise ${\mathbf{n}}_{m}$ is colored
noise with a covariance matrix of
${\mathbf{C}}=\sigma^{2}\text{diag}\\{{\mathbf{A}}_{1}{\mathbf{A}}_{1}^{H},{\mathbf{A}}_{2}{\mathbf{A}}_{2},\cdots,{\mathbf{A}}_{P}{\mathbf{A}}_{P}^{H}\\}$,
we use a matrix ${\mathbf{D}}\in\mathbb{C}^{PN_{\text{RF}}\times
PN_{\text{RF}}}$ to pre-white the received signals. Here, ${\mathbf{D}}$ is
acquired by calculating the eigenvalue decomposition
${\mathbf{C}}=\sigma^{2}{\mathbf{S}}\bm{\Sigma}{\mathbf{S}}^{H}$, where
${\mathbf{D}}={\mathbf{S}}\bm{\Sigma}^{\frac{1}{2}}$. Left multiplying the
received signals ${\mathbf{Y}}$ by ${\mathbf{D}}^{-1}$, we have

$\displaystyle\bar{{\mathbf{Y}}}={\mathbf{D}}^{-1}{\mathbf{Y}}={\mathbf{D}}^{-1}{\mathbf{A}}{\mathbf{W}}\bar{{\mathbf{H}}}+{\mathbf{D}}^{-1}{\mathbf{N}}=\bm{\Psi}\bar{{\mathbf{H}}}+\bar{{\mathbf{N}}},$
(17)

where $\bm{\Psi}={\mathbf{D}}^{-1}{\mathbf{A}}{\mathbf{W}}$ and
$\bar{{\mathbf{N}}}={\mathbf{D}}^{-1}{\mathbf{N}}$. Notice that here, the
covariance of each column of $\bar{{\mathbf{N}}}$ becomes
${\mathbf{D}}^{-1}{\mathbf{C}}{\mathbf{D}}^{-H}=\sigma^{2}\bm{\Sigma}^{-\frac{1}{2}}{\mathbf{S}}^{-1}{\mathbf{S}}\bm{\Sigma}{\mathbf{S}}^{H}{\mathbf{S}}^{-H}\bm{\Sigma}^{-\frac{1}{2}}=\sigma^{2}{\mathbf{I}}_{PN_{\text{RF}}}$,
which is white noise.

In the second stage of path detection, we first initialize the residue matrix
${\mathbf{R}}=\mathbb{C}^{PN_{\text{RF}}\times M}$ as
${\mathbf{R}}=\bar{{\mathbf{Y}}}$. The support set $\Upsilon$ is initialized
as $\Upsilon=\\{\emptyset\\}$, which will be used to store the detected
support indices at carrier frequency $f_{c}$.

In steps 5-6, according to (15) and (16), we generate the angle-domain and
distance-domain linear patterns as

$\displaystyle\Gamma(n_{a},m)=\arg\min_{n}|\bar{\theta}_{n}-\frac{f_{m}}{f_{c}}\bar{\theta}_{n_{a}}|,$
(18)
$\displaystyle\Lambda(n_{d},m)=\arg\min_{n}|\bar{\alpha}_{n}-\frac{f_{m}}{f_{c}}\bar{\alpha}_{n_{d}}|,$
(19)

where $n_{a}\in\\{1,\cdots,N_{a}\\}$, $n_{d}\in\\{1,\cdots,N_{d}\\}$, and
$m\in\\{1,\cdots,M\\}$. Here, all parameters $\bar{\theta}_{n_{a}}$,
$\bar{\theta}_{n}$, $\bar{\alpha}_{n_{d}}$, and $\bar{\alpha}_{n}$ belong to
the sampled locations in ${\mathbf{W}}$. $(\Gamma(n_{a},m),\Lambda(n_{d},m))$
denotes the support index of $(\bar{\theta}_{n_{a}},\bar{\alpha}_{n_{d}})$ at
frequency $f_{m}$. To elaborate on steps 5-6, based on the CS theory, our
target is to find out a polar-domain sampled location
$(\bar{\vartheta}_{n_{a,l}^{\star}},\bar{r}_{n_{d,l}^{\star},{n_{a,l}^{\star}}})$
that is nearest to the physical location $(\vartheta_{l},r_{l})$. As indicated
in Lemma 2, each sampled location
$(\bar{\vartheta}_{n_{a}},\bar{r}_{n_{d},n_{a}})$ corresponds to a certain
bilinear pattern over frequencies. Hence, we use $\Gamma({n_{a}},m)$ and
$\Lambda({n_{d}},m)$ to represent these patterns. Then, $\Gamma({n_{a}},m)$
and $\Lambda({n_{d}},m)$ allow us to accumulate the coherence between
$(\bar{\theta}_{n_{a}},\bar{\alpha}_{n_{d}})$ and a physical location
$(\vartheta_{l},r_{l})$ from the entire bandwidth. The nearest sample
$(\bar{\vartheta}_{n_{a,l}^{\star}},\bar{r}_{n_{d,l}^{\star},{n_{a,l}^{\star}}})$
or $(\bar{\theta}_{n_{a,l}^{\star}},\bar{\alpha}_{n_{d,l}^{\star}})$ can be
selected by capturing the largest coherence.

To be more specific, for the $l$-th path component, the physical location
$(\vartheta_{l},r_{l})$ is estimated relying on the above idea. Inspired by
the OMP and SOMP methods, we use the product of $\bm{\Psi}^{H}$ and
${\mathbf{R}}$ to calculate the correlation matrix
${\mathbf{U}}=\bm{\Psi}^{H}{\mathbf{R}}=[\mathbf{u}_{1},\mathbf{u}_{2},\cdots,\mathbf{u}_{M}]\in\mathbb{C}^{N_{a}N_{d}\times
M}$ in step 8. After that, in step 9, the bilinear pattern is utilized to
capture the power of the correlation matrix for all frequencies. For each
sampled $(\bar{\vartheta}_{n_{a}},\bar{r}_{n_{d},n_{a}})$ with a support index
$(n_{a},n_{d})$ at carrier frequency $f_{c}$, its support index becomes
$(\Gamma({n_{a}},m),\Lambda({n_{d}},m))$ at frequency $f_{m}$. Besides, as the
dimension of each column ${\mathbf{u}}_{m}$ of ${\mathbf{U}}$ is
$N_{a}N_{d}\times 1$, the support index
$(\Gamma({n_{a}},m),\Lambda({n_{d}},m))$ corresponds to the
$\left((\Lambda(n_{d},m)-1)N_{a}+\Gamma({n_{a}},m)\right)$-th element of
${\mathbf{u}}_{m}$. As a consequence, the power of the correlation matrix for
location $(\bar{\vartheta}_{n_{a}},\bar{r}_{n_{d},n_{a}})$ can be written as
$\sum_{m=1}^{M}\|u_{m}\left((\Lambda(n_{d},m)-1)N_{a}+\Gamma({n_{a}},m)\right)\|^{2}$,
and thus the optimal support index of the $l$-th path component is determined
as

$\displaystyle(n_{a,l}^{\star},n_{d,l}^{\star})=\arg\max_{(n_{a},n_{d})}\sum_{m=1}^{M}\|u_{m}\left((\Lambda(n_{d},m)-1)N_{a}+\Gamma({n_{a}},m)\right)\|^{2}.$
(20)

After detecting $(n_{a,l}^{\star},n_{d,l}^{\star})$, we add it to the support
set $\Upsilon$ in step 10. Basically, $\Upsilon$ represents the detected
sparse channel supports at the carrier frequency. Then in step 12, $\Upsilon$
is transformed to the sparse channel support set $\Upsilon_{m}$ at frequency
$f_{m}$. We define each element in $\Upsilon_{m}$ as the column index of a
sparse support index of $\bar{{\mathbf{h}}}_{m}$, which can be written as

$\displaystyle\Upsilon_{m}=\\{\left(\Lambda(n_{d},m)-1\right)N_{a}+\Gamma({n_{a}},m)|(n_{a},n_{d})\in\Upsilon\\}.$
(21)

Next, in step 13, we denote $\hat{\bar{{\mathbf{h}}}}_{m}$ as the estimated
polar-domain channel at frequency $f_{m}$. Its non-zero elements can be
calculated by using $\Upsilon_{m}$ through the LS algorithm as

$\displaystyle\hat{\bar{{\mathbf{h}}}}_{m}={\mathbf{0}}_{N_{a}N_{d}\times
1},\quad\hat{\bar{{\mathbf{h}}}}_{m}(\Upsilon_{m})=\bm{\Psi}^{\dagger}(:,\Upsilon_{m})\bar{{\mathbf{y}}}_{m}.$
(22)

Here, matrix $\bm{\Psi}(:,\Upsilon_{m})$ is composed of the column vectors in
$\bm{\Psi}$ indexed by the set $\Upsilon_{m}$. After that, we remove the
impact of all detected $l$ paths in the current iteration from the received
signals to update the residue matrix as

$\displaystyle{\mathbf{R}}=\bar{{\mathbf{Y}}}-\bm{\Psi}[\hat{\bar{{\mathbf{h}}}}_{1},\hat{\bar{{\mathbf{h}}}}_{2},\cdots,\hat{\bar{{\mathbf{h}}}}_{M}].$
(23)

Steps 8-15 discussed above are repeated $\hat{L}$ times until all path
components are detected. Eventually, in step 17, the antenna-domain wideband
channel $\hat{{\mathbf{H}}}$ is recovered based on these sparse polar-domain
channels as

$\displaystyle\hat{{\mathbf{H}}}={\mathbf{W}}[\hat{\bar{{\mathbf{h}}}}_{1},\hat{\bar{{\mathbf{h}}}}_{2},\cdots,\hat{\bar{{\mathbf{h}}}}_{M}].$
(24)

It is worth mentioning that although the proposed algorithm is inspired by the
conventional CS-based methods, it is able to exploit both the near-field and
beam split effects embedded in the channel model while existing methods
cannot. To be specific, existing far-field wideband channel estimation
algorithms rely on the angle-domain sparsity and only exploit the angle-domain
linear pattern $\Gamma(n_{a},m)$, thus these methods will undergo a serious
performance loss in the near-field environments. Moreover, existing near-field
narrowband channel estimation algorithms certainly take advantage of the
polar-domain sparsity, but they assume the channel sparse support sets are the
same for all frequencies, which will result in poor estimation accuracy when
the bandwidth is large. On the contrary, the proposed algorithm takes good
advantage of the frequency-dependent sparse structure of near-field wideband
channels. Leveraging such a bilinear pattern, our scheme is expected to
achieve higher estimation accuracy in wideband XL-MIMO systems, which will be
demonstrated in the simulation section222 Notice that the proposed scheme can
achieve higher accuracy on the estimation of the location parameter
$(\vartheta_{l},r_{l})$ at all frequencies than existing channel estimation
schemes. Thus, our scheme can also be used to improve the time of arrival
(ToA) estimation accuracy..

###  3.3 Complexity Analysis

Table 1:  Computational Complexity Algorithm |  Computiational Complexity $\mathcal{O}(\cdot)$
---|---
Angle-domain SOMP [12] |  $\mathcal{O}(\hat{L}N_{a}PN_{\text{RF}}M)+\mathcal{O}(\hat{L}^{3}MPN_{RF}+\hat{L}^{4})+\mathcal{O}(\hat{L}^{2}PN_{\text{RF}}M)$
BSPD algorithm [16] |  $\mathcal{O}(\hat{L}N_{a}PN_{\text{RF}}M)+\mathcal{O}(\hat{L}^{3}MPN_{RF}+\hat{L}^{4})+\mathcal{O}(\hat{L}^{2}PN_{\text{RF}}M)$
Polar-domain SOMP [20] |  $\mathcal{O}(\hat{L}N_{a}N_{d}PN_{\text{RF}}M)+\mathcal{O}(\hat{L}^{3}MPN_{RF}+\hat{L}^{4})+\mathcal{O}(\hat{L}^{2}PN_{\text{RF}}M)$
Proposed BPD algorithm |  $\mathcal{O}(\hat{L}N_{a}N_{d}PN_{\text{RF}}M)+\mathcal{O}(\hat{L}^{3}MPN_{RF}+\hat{L}^{4})+\mathcal{O}(\hat{L}^{2}PN_{\text{RF}}M)$

In this subsection, the computational complexity of the BPD based channel
estimation algorithm is analyzed, where we mainly count the number of complex
multiplications. As indicated in Algorithm 1, the complexity mainly comes from
the iteration procedure of steps 8, 9, 13, and 15.

In step 8, the product of matrices
$\bm{\Psi}^{H}\in\mathbb{C}^{N_{a}N_{d}\times PN_{\text{RF}}}$ and
${\mathbf{R}}\in\mathbb{C}^{PN_{\text{RF}}\times M}$ has a complexity in the
order of $\mathcal{O}(N_{a}N_{d}PN_{\text{RF}}M)$. Notice that the number of
sampled angles $N_{a}$ and distances $N_{d}$ is usually proportional to the
number of antennas $N$ [20].

In step 9, we need to calculate the power
$\sum_{m=1}^{M}\|u_{m}\left((\Lambda(n_{d},m)-1)N_{a}+\Gamma({n_{a}},m)\right)\|^{2}$
for $N_{a}N_{d}$ times, whose complexity is in the order of
$\mathcal{O}(N_{a}N_{d}M)$.

The complexity of step 13 is dominated by the calculation of matrix
${\bm{\Psi}}^{\dagger}(:,\Upsilon_{m})\in\mathbb{C}^{l\times PN_{\text{RF}}}$,
which can be regarded as
${\bm{\Psi}}^{\dagger}(:,\Upsilon_{m})=({\bm{\Psi}}^{H}(:,\Upsilon_{m}){\bm{\Psi}}(:,\Upsilon_{m}))^{-1}{\bm{\Psi}}^{H}(:,\Upsilon_{m})$.
Hence, the complexity of calculating ${\bm{\Psi}}^{\dagger}(:,\Upsilon_{m})$
consists of three computation steps: the product of
${\bm{\Psi}}^{H}(:,\Upsilon_{m})$ and ${\bm{\Psi}}(:,\Upsilon_{m})$, the
inverse of ${\bm{\Psi}}^{H}(:,\Upsilon_{m}){\bm{\Psi}}(:,\Upsilon_{m})$, and
the product of
$({\bm{\Psi}}^{H}(:,\Upsilon_{m}){\bm{\Psi}}(:,\Upsilon_{m}))^{-1}$ and
${\bm{\Psi}}^{H}(:,\Upsilon_{m})$. Therefore, step 13 has a complexity in the
order of
$\mathcal{O}(l^{2}PN_{RF}+l^{3}+l^{2}PN_{RF})=\mathcal{O}(l^{2}PN_{RF}+l^{3})$.
As step 13 is carried out $M$ times, its overall complexity is in the order of
$\mathcal{O}(M(l^{2}PN_{RF}+l^{3}))$.

Next, in step 15, each column vector $\hat{\bar{{\mathbf{h}}}}_{m}$ only has
$l$ non-zero elements. As a result, the left multiplying
$[\hat{\bar{{\mathbf{h}}}}_{1},\hat{\bar{{\mathbf{h}}}}_{2},\cdots,\hat{\bar{{\mathbf{h}}}}_{M}]$
by matrix $\bm{\Psi}$ has a complexity of $\mathcal{O}(PN_{\text{RF}}lM)$.

Eventually, steps 8, 9, 13, 15 are executed $\hat{L}$ times with
$l=1,2,\cdots,\hat{L}$. Therefore, the total complexity can be summarized as

$\displaystyle\>\mathcal{O}(\hat{L}N_{a}N_{d}PN_{\text{RF}}M)+\mathcal{O}(\hat{L}N_{a}N_{d}M)+\mathcal{O}(\hat{L}^{3}MPN_{RF}+\hat{L}^{4})+\mathcal{O}(\hat{L}^{2}PN_{\text{RF}}M)$
$\displaystyle=$
$\displaystyle\>\mathcal{O}(\hat{L}N_{a}N_{d}PN_{\text{RF}}M)+\mathcal{O}(\hat{L}^{3}MPN_{RF}+\hat{L}^{4})+\mathcal{O}(\hat{L}^{2}PN_{\text{RF}}M).$
(25)

In Table 1, we have provided the computational complexity comparison for
different channel estimation algorithms, including the angle-domain SOMP
algorithm [12] for far-field narrowband channel estimation, the BSPD algorithm
[16] for far-field wideband channel estimation, and the polar-domain SOMP
algorithm [20] for near-field narrowband channel estimation. It is clear from
Table 1 that the proposed BPD algorithm has the same complexity as the polar-
domain SOMP scheme in [20]. This is because the main difference between these
two algorithms lies in step 9. Specifically, the polar-domain SOMP technique
uses the common sparse support set assumption to accumulate the power while
the proposed algorithm exploits the bilinear pattern to accumulate the power,
both of which have the same complexity $\mathcal{O}(N_{a}N_{d}M)$.

Moreover, the complexity of near-field algorithms is higher than that of far-
field algorithms, since the far-field methods only need to recover the AoA
parameters while the near-field schemes have to detect the AoA and distance
parameters simultaneously. To illustrate, the complexity difference comes from
the first item in (3.3), i.e., $\mathcal{O}(\hat{L}N_{a}N_{d}PN_{\text{RF}}M)$
for near-field methods and $\mathcal{O}(\hat{L}N_{d}PN_{\text{RF}}M)$ for far-
field methods. As the far-field channel estimation solely cares about the AoA
information, the number of sampled distances $N_{d}$ in [12, 16] can be
regarded as 1, giving rise to the complexity of
$\mathcal{O}(\hat{L}N_{d}PN_{\text{RF}}M)$. Despite the lower computational
complexity, the far-field algorithms neglect the essential distance
information, so they can hardly achieve accurate XL-MIMO channel estimation,
which is demonstrated in the next section.

##  4 Simulation Results

Table 2: Simulation Configurations Parameter | Value | Parameter | Value
---|---|---|---
Number of BS antennas $N$ | 256 |  Number of angle-domain samples $N_{a}$ | 256
Number of RF chains $N_{\text{RF}}$ | 4 |  Number of distance-domain samples $N_{d}$ | 14
Number of users $K$ | 4 |  Parameter $\beta$ | 0.8
Number of subcarriers $M$ | 256 |  Number of channel paths $L$ | 6
Carrier frequency $f_{c}$ | 100 GHz |  Number of paths to be detected $\hat{L}$ | 12
The distribution of path gain $g_{l,m}$ | $\mathcal{CN}(0,1)$ |  The distribution of angle $\vartheta_{l}$ | $\mathcal{U}\left(-\frac{\pi}{2},\frac{\pi}{2}\right)$

In this section, we present simulation results to demonstrate the performance
of the proposed algorithm. A wideband XL-MIMO system is considered, and some
of the simulation configurations are shown in Table 2.  According to these
configurations, the number of BS’s antennas is $N=256$ and the carrier
wavelength is $\lambda_{c}=\frac{c}{f_{c}}=3$ mm. Therefore, the antenna
spacing is $d=\frac{\lambda_{c}}{2}=1.5$ mm, and the BS array aperture is
$D=(N-1)d\approx Nd=\frac{N\lambda_{c}}{2}=0.384$ m. Then, it can be derived
that the Fraunhofer distance is $\frac{2D^{2}}{\lambda_{c}}=98.3$ m.  In
addition, the distance parameters $r_{l}$ in (1) are randomly generated from
the uniform distribution $\mathcal{U}(R_{\min},R_{\max})$, where the smallest
distance $R_{\min}$ and the largest distance $R_{\min}$ have different values
for different simulations. The signal-to-noise radio is defined as
$\text{SNR}=\mathbb{E}\left(\|{\mathbf{H}}\|_{F}^{2}/\|\bar{{\mathbf{N}}}\|_{F}^{2}\right)$.
The compared benchmark channel estimation algorithms are as follows: 1) the LS
algorithm; 2) the far-field narrowband channel estimation schemes, including
the angle-domain OMP [10] and SOMP [12] algorithms; 3) the BSPD based far-
field wideband channel estimation [16]; 4) the near-field narrowband channel
estimation schemes, including the polar-domain OMP and SOMP algorithms [20].
We use the performance of normalized mean square error (NMSE) to evaluate
different algorithms, which is defined as
$\text{NMSE}=\mathbb{E}\left(\frac{\|{\mathbf{H}}-\hat{{\mathbf{H}}}\|_{F}^{2}}{\|{\mathbf{H}}\|_{F}^{2}}\right)$.
Besides, 300 Monte Carlo experiments are carried out to plot each figure.

Figure 5:  NMSE performance against link distance.

First of all, the NMSE performance against distance is shown in Fig. 5 to
explain the influence of near-field effect. In this figure, we assume the
distance $r_{l}$ is growing from $5$ meters to 100 meters, where we set
$R_{\min}=R_{\max}=r_{l}$. The other parameter settings are as follows: 1) the
SNR is fixed as $5$ dB, 2) the bandwidth is $B=10$ GHz, 3) and the pilot
overhead is $P=32$, with an observation dimension of $PN_{\text{RF}}=128$. One
can observe from Fig. 5 that the achieved NMSE of all far-field channel
estimation methods (angle-domain OMP, SOMP, and BSPD) degrades with the
decrease in distance. This is because these algorithms tailored for far-field
communications neglect the impact of spherical wavefront. Besides, as the LS
algorithm is able to work on all kinds of channels, its performance is robust
to different distances. Similarly, all near-field channel estimation
algorithms also have a stable NMSE in both far-field and near-field
environments by exploiting the polar-domain sparsity. Despite this distance-
robust advantage, the achieved NMSE of the LS, polar-domain OMP, and SOMP
methods are not satisfactory, because they fail to utilize the underlying beam
split structure in wideband systems. Fortunately, the proposed BPD algorithm
well captures the polar-domain frequency-dependent sparse support structure,
and thus it outperforms all benchmark algorithms.

Figure 6:  NMSE performance against system bandwidth.

Then, Fig. 6 illustrates the NMSE performance against system bandwidth, which
is growing from $100$ MHz to 10 GHz. The other parameters settings are as
follows: $\text{SNR}=5$ dB, the pilot overhead is $P=32$, $R_{\min}=10$
meters, and $R_{\max}=30$ meters. As the common sparse support set
characteristic is not valid in wideband systems, the NMSE performance of
angle-domain and polar-domain SOMP algorithms becomes worse and worse with the
increase of bandwidth. In addition, the LS, angle-domain OMP, and polar-domain
OMP schemes make no assumption on the wideband channel structure, where the
channel of each subcarrier is independently estimated. As a consequence, these
algorithms achieve a stable but non-satisfactory channel estimation accuracy
for different bandwidth conditions. Moreover, the BSPD method accurately
describes the angle-domain beam split pattern but ignores the distance-domain
beam split pattern, so it is able to slightly improve the estimation accuracy
for different bandwidth cases. On the other hand, we can observe from Fig. 6
that the proposed BPD algorithm realizes both stable and the most accurate
estimation performance in all bandwidth cases. This is because the BPD scheme
fully exploits the frequency-dependent sparse structure resulting from the
near-field beam split effect.

Figure 7:  NMSE performance against SNR.

In Fig. 7, we evaluate the achieved NMSE performance in different SNR
conditions. Here, the SNR grows from -5 dB to 15 dB. The other parameter
settings are set as below: $\text{B}=10$ GHz, $R_{\min}=10$ meters,
$R_{\max}=30$ meters, and $P=32$. The accuracy achieved by all considered
algorithms improves with the increase of SNR. It is clear that the proposed
BPD method significantly outperforms all compared benchmarks, especially in
high SNR cases. For example, when the SNR is 7 dB, around 5 dB improvement of
NMSE is realized compared to the BSPD method and the polar-domain OMP method.

Figure 8:  NMSE performance against pilot overhead.

Eventually, in Fig. 8, the NMSE performance with respect to the pilot overhead
$P$ is plotted. As illustrated in Fig. 8, the pilot overhead $P$ is increasing
from 4 to 32, which corresponds to a compression ratio
$\frac{N}{PN_{\text{RF}}}$ declining from $16$ to $2$. Then, the other
parameters are set as: $\text{SNR}=5$ dB, $R_{\min}=10$ meters, $R_{\max}=30$
meters, and $B=10$ GHz. It is clear that the proposed BPD based method
outperforms other existing channel estimation algorithms for all considered
pilot lengths. This fact implies that the proposed method can be used to
reduce the pilot overhead. Specifically, take the NMSE of -9 dB as a baseline,
around 36 pilot overhead is required for the BSPD to achieve this NMSE
baseline. On the contrary, only 12 pilot length is enough for the proposed BPD
method to reach this baseline. In this case, our scheme has the capability of
reducing the pilot overhead by $66\%$.

##  5 Conclusions

In this paper, we have investigated the channel estimation for wideband XL-
MIMO communications in the presence of near-field beam split effect.
Specifically, we first revealed the bilinear pattern of the near-field beam
split effect, which indicated that the polar-domain sparse support set for
each near-field channel path shows a linear structure over frequencies. Then,
we have proposed a BPD based channel estimation algorithm to recover each
near-field channel path, by using the bilinear pattern to accumulate the
largest polar-domain power from the entire bandwidth. Simulation results
demonstrated our scheme is capable of achieving high channel estimation
accuracy in all far-field/near-field/narrowband/wideband conditions. For
future works, the bilinear pattern discovered in this paper can potentially be
extended to tackle the relevant channel estimation issues in various near-
field wideband communication scenarios, such as reconfigurable intelligent
surface (RIS) communications [26] and cell-free massive MIMO communications
[27].

This work was supported in part by the National Key Research and Development
Program of China (Grant No. 2020YFB1805005), 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.

## References

  * [1] T. S. Rappaport, Y. Xing, O. Kanhere, S. Ju, A. Madanayake, S. Mandal, A. Alkhateeb, and G. C. Trichopoulos. Wireless communications and applications above 100 GHz: Opportunities and challenges for 6G and beyond. IEEE Access, 2019, 7: 78729-78757.
  * [2] Z. Zhang, Y. Xiao, Z. Ma, M. Xiao, Z. Ding, X. Lei, G. K. Karagiannidis, and P. Fan. 6G wireless networks: Vision, requirements, architecture, and key technologies. IEEE Veh. Technol, 2019, 14: 28-41.
  * [3] B. Ning, Z. Tian, Z. Chen, C. Han, J. Yuan, and S. Li. Prospective beamforming technologies for ultra-massive MIMO in terahertz communications: A tutorial. arXiv preprint arXiv:2107.03032, 2021.
  * [4] C. Lin and G. Y. L. Li. Terahertz communications: An array-of-subarrays solution. IEEE Commun. Mag., 2016, 54: 124-131.
  * [5] M. R. Akdeniz, Y. Liu, M. K. Samimi, S. Sun, S. Rangan, T. S. Rappaport, and E. Erkip. Millimeter wave channel modeling and cellular capacity evaluation. IEEE J. Sel. Areas Commun., 2014, 32: 1164-1179.
  * [6] I. F. Akyildiz, H. Chong, and S. Nie. Combating the distance problem in the millimeter wave and terahertz frequency bands. IEEE Commun. Mag., 2018, 56: 102-108.
  * [7] I. F. Akyildiz, J. M. Jornet, and H. Chong. Terahertz band: Next frontier for wireless communications. Phys. Commun., 2014, 12: 16-32.
  * [8] F. Sohrabi and Y. Wei. Hybrid digital and analog beamforming design for large-scale antenna arrays. IEEE J. Sel. Topics Signal Process., 2016, 10: 501-513.
  * [9] R. W. Heath, N. González-Prelcic, S. Rangan, W. Roh, and A. M. Sayeed. An overview of signal processing techniques for millimeter wave MIMO systems. IEEE J. Sel. Topics Signal Process., 2016, 10: 436-453.
  * [10] J. Lee, G.-T. Gil, and Y. H. Lee. Channel estimation via orthogonal matching pursuit for hybrid MIMO systems in millimeter wave communications. IEEE Trans. Commun., 2016, 64: 2370-2386.
  * [11] C. Huang, L. Liu, C. Yuen, and S. Sun. Iterative channel estimation using LSE and sparse message passing for mmWave MIMO systems. IEEE Trans. Signal Process., 2019, 67: 245-259.
  * [12] J. Rodríguez-Fernández, N. González-Prelcic, K. Venugopal, and R. W. Heath. Frequency-domain compressive channel estimation for frequency-selective hybrid millimeter wave MIMO systems. IEEE Trans. Wireless Commun., 2018, 17: 2946-2960.
  * [13] Z. Gao, C. Hu, L. Dai, and Z. Wang. Channel estimation for millimeter-wave massive MIMO with hybrid precoding over frequency-selective fading channels. IEEE Commun. Lett., 2016, 20: 1259-1262.
  * [14] Y. C. Eldar. Sampling Theory: Beyond Bandlimited Systems. Cambridge University Press, 2015.
  * [15] A. Liao, Z. Gao, D. Wang, H. Wang, H. Yin, D. W. K. Ng, and M.-S. Alouini. Terahertz ultra-massive MIMO-based aeronautical communications in space-air-ground integrated networks. IEEE J. Sel. Areas Commun., 2021, 39: 1741-1767.
  * [16] J. Tan and L. Dai. Wideband channel estimation for THz massive MIMO. China Commun., 2021, 18: 66-80.
  * [17] K. Dovelos, M. Matthaiou, H. Q. Ngo, and B. Bellalta. Channel estimation and hybrid combining for wideband terahertz massive MIMO systems. IEEE J. Sel. Areas Commun., 2021, 39: 1604-1620.
  * [18] J. Yang, Y. Zeng, S. Jin, C.-K. Wen, and P. Xu. Communication and localization with extremely large lens antenna array. IEEE Trans. Wireless Commun., 2021, 20: 3031-3048.
  * [19] 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., 2021, 69: 7063-7076.
  * [20] M. Cui and L. Dai. Channel estimation for extremely large-scale MIMO: Far-field or near-field?. IEEE Trans. Commun., 2022, 70: 2663-2677.
  * [21] M. Cui, Z. Wu, Y. Lu, X. Wei, and L. Dai. Near-field communications for 6G: Fundamentals, challenges, potentials, and future directions. IEEE Commun. Mag., 2022.
  * [22] K. T. Selvan and R. Janaswamy. Fraunhofer and fresnel distances: Unified derivation for aperture antennas. IEEE Antennas Propag. Mag., 2017, 59: 12-15.
  * [23] H. Zhang, N. Shlezinger, F. Guidi, D. Dardari, M. F. Imani, and Y. C. Eldar. Beam focusing for near-field multi-user MIMO communications. IEEE Trans. Wireless Commun., 2022, 21:7476-7490.
  * [24] M. Cui, L. Dai, Z. Wang, S. Zhou, and N. Ge. Near-field rainbow: Wideband beam training for XL-MIMO. IEEE Trans. Wireless Commun., 2021.
  * [25] Z. Zhou, X. Gao, J. Fang, and Z. Chen. Spherical wave channel and analysis for large linear array in LoS conditions. In: Proceedings IEEE Globecom Workshops, San Diego, 2015, 1–6.
  * [26] L. Wei, C. Huang, G. C. Alexandropoulos, C. Yuen, Z. Zhang, and M. Debbah. Channel estimation for RIS-empowered multi-user MISO wireless communications. IEEE Trans. Commun., 2021, 69: 4144-4157.
  * [27] Z. Zhang and L. Dai. A joint precoding framework for wideband reconfigurable intelligent surface-aided cell-free network. IEEE Trans. Signal Process., 2021, 69: 4085-4101.