RReferences of Rebuttal

# Cross-Field Channel Estimation for Ultra Massive-MIMO THz Systems

Simon Tarboush, Anum Ali, , and Tareq Y. Al-Naffouri This work was supported
by the KAUST Office of Sponsored Research. A preliminary version of this work
was presented at the IEEE International Conference on Acoustics, Speech and
Signal Processing (ICASSP 2023) [1]. S. Tarboush is a researcher from
Damascus, Syria (e-mail: simon.w.tarboush@gmail.com). A. Ali is with the
Standards and Mobility Innovation Laboratory, Samsung Research America, Plano,
TX 75023, USA (e-mail: anum.ali@samsung.com). T. Y. Al-Naffouri is with the
Department of Computer, Electrical and Mathematical Sciences and Engineering
(CEMSE), King Abdullah University of Science and Technology (KAUST), Thuwal,
Makkah Province, Kingdom of Saudi Arabia, 23955-6900 (e-mail:
tareq.alnaffouri@kaust.edu.sa).

###### Abstract

The large bandwidth combined with ultra-massive multiple-input multiple-output
(UM-MIMO) arrays enables terahertz (THz) systems to achieve terabits-per-
second throughput. The THz systems are expected to operate in the near,
intermediate, as well as the far-field. As such, channel estimation strategies
suitable for the near, intermediate, or far-field have been introduced in the
literature. In this work, we propose a cross-field, i.e., able to operate in
near, intermediate, and far-field, compressive channel estimation strategy.
For an array-of-subarrays (AoSA) architecture, the proposed method compares
the received signals across the arrays to determine whether a near,
intermediate, or far-field channel estimation approach will be appropriate.
Subsequently, compressed estimation is performed in which the proximity of
multiple subarrays (SAs) at the transmitter and receiver is exploited to
reduce computational complexity and increase estimation accuracy. Numerical
results show that the proposed method can enhance channel estimation accuracy
and complexity at all distances of interest.

###### Index Terms:

Near-field, far-field, intermediate-field, spherical wave model, planar wave
model, hybrid spherical-planar wave model, array-of-subarrays.

##  I Introduction

The terahertz (THz)-band, i.e., $0.3\\!-\\!10\,\mathrm{THz}$, is a promising
candidate to enable the diverse use cases in the sixth generation (6G) mobile
networks [2, 3, 4]. Specifically, THz communications can support a plethora of
novel services and applications, e.g., ubiquitous connectivity, inter-chip
communications, accurate localization, and sensing [5, 6]. Moreover, the large
bandwidths available in THz can satisfy the extreme throughput requirement (up
to terabits-per-second (Tbps)) for super-high-definition (SHD), extremely
high-definition (EHD) video, and full sensory experience (i.e., hearing,
vision, smell, touch, and taste) [7].

The THz channels have various peculiarities [8, 9, 10, 11]. First, the
propagation is dominated by the line-of-sight (LoS), and non-LoS paths can be
$\sim$5-15 dB weaker than LoS due to the increased reflection and diffraction
loss [8]. In a typical THz indoor scenario, the number of channel paths is
typically less than ten, and this number decreases with high-gain antennas or
massive arrays [10]. Therefore, the THz channels are extremely sparse in the
time and angle domains [9, 11]. Second, the THz channels have a frequency and
distance dependent frequency selectivity (even in LoS) due to molecular
absorption [10]. Third, due to the very high free space path loss, the
communication distance is limited [12]. Thus, ultra-massive multiple-input
multiple-output (UM-MIMO) antenna arrays at the transmitter (Tx) and receiver
(Rx) [12] are necessary to get high beamforming gain and achieve longer
communication distances. The THz system design, however, becomes challenging
due to large bandwidth and antenna arrays. As such, several hybrid beamforming
THz transceiver architectures have been proposed [13, 14]. The sub-connected
array-of-subarrays (AoSA) architecture is one promising candidate [15, 16,
13]. In this architecture, the Tx and Rx UM-MIMO arrays are divided into
subarrays (SAs), where each SA is connected to one radio frequency (RF)-chain
exclusively. This provides several advantages, e.g., reduced complexity and
power consumption, and also makes channel estimation simpler [1, 17]. On the
downside, some spectral efficiency is lost compared to the fully-connected
structures.

The unique THz channel characteristics and UM arrays add a new dimension to
channel modeling and estimation. It is well understood that while the planar
wave model (PWM) is sufficient for far-field, the appropriate channel model
for near-field is the spherical wave model (SWM). Though SWM is accurate for
channel modeling, it makes the channel estimation more complex [18]. An
alternative to SWM is a hybrid spherical-planar wave model (HSPWM) [10, 17].
The HSPWM provides, compared to PWM and SWM, a good trade-off in modeling
accuracy and complexity, and is particularly suitable for hybrid beamforming
AoSA architecture [10, 17].

The prior work [19, 10] has shown that the HSPWM approximates the SWM with an
approximation error of less than $-10\,\mathrm{dB}$, across a wide range of
communication distances, array apertures, and carrier frequencies. In the
context of localization, the angle and distance Cramér-Rao bounds for SWM and
HSPWM, have been studied in [20]. In this work, we refer to the region where
HSPWM accurately approximates SWM - but PWM does not - as the intermediate-
field.

With different channel models applicable to different distances, namely SWM
for near-field, HSPWM for intermediate-field, and PWM for far-field, the
cross-field channel estimation problem arises. Specifically, the cross-field
channel estimation problem entails $\mathrm{(i)}$ determining whether the user
is currently in the near, intermediate, or far-field, and $\mathrm{(ii)}$
using the appropriate channel model and estimation method accordingly.

###  I-A Contributions

In this work, we develop a channel estimation strategy for a sub-connected
AoSA THz system. The main contributions of this work are

  * •

Formulating and solving the cross-field channel estimation problem. This
problem arises as the users transition between near, intermediate, and far-
field scenarios, and hence it is crucial to employ the suitable channel
estimation approach depending on their current location.

  * •

Developing a method to determine if the user is in the near, intermediate, or
far-field. This determination permits the use of an appropriate channel
estimation strategy which in turn leads to a good trade-off in computational
complexity and channel estimation accuracy.

  * •

Proposing a reduced dictionary (RD) approach to enhance the estimation
accuracy and reduce the computational complexity of channel estimation in the
near-, intermediate-, and far-field (i.e., for SWM, HSPWM, and PWM). The
proposed RD method exploits the proximity of different SAs in an AoSA system
to reduce the search space after an initial channel estimate is available.

###  I-B Related Work

For far-field channel estimation - based on PWM - the sparsity of the THz
channels is exploited using compressed sensing (CS). Specifically, a modified
simultaneous orthogonal matching pursuit (SOMP) is proposed in [21] for
wideband THz channel estimation.

For near-field channel estimation - based on SWM - a polar codebook is
introduced, leveraging uniform angle sampling and non-uniform distance
sampling [22]. Alternatively, [23] presents a distance-parameterized angle-
domain codebook. The mixed near-field channel scenario separately models LoS
and non-LoS components [24]. Channel estimation is divided into two sub-
problems in [24]. Firstly, the LoS path component is estimated by searching
through a collection of distance and angle points to acquire coarse on-grid
parameters. This initial estimation is then refined using an off-grid
approach, utilizing an iterative gradient descent optimization method. The
performance of the first part heavily depends on the coarse estimation. In
practice, acquiring these initial distance and angle estimates before the
channel estimation phase begins can be quite challenging. Secondly, the non-
LoS paths are estimated using the OMP algorithm in combination with polar-
domain codebooks.

The intermediate-field channel estimation - based on HSPWM - has been
considered in only a handful of articles. Specifically, a CS-based approach is
taken in [25], and codebooks are designed for SAs, whereas data-driven
strategies based on deep learning are used in [19].

The hybrid-field channel means the coexistence of both far-field and near-
field path components [26]. The hybrid-field channel estimation approaches
leverage the sparsity of both its far- and near-field components in the
angular and polar domains. In [26], a two-stage CS algorithm is introduced for
hybrid-field channel estimation. In the first phase, the algorithm leverages
OMP combined with the angular-domain codebook to estimate the far-field
components. Subsequently, after eliminating the influence of these components
from the received signal, the second phase detects near-field paths using OMP
in conjunction with a polar-domain codebook. An alternative approach in [27]
introduces a support detection mechanism into the far-field estimation while
retaining the procedure of [26] for near-field estimation. These methods
exhibit high computational complexity since they do not take into account
where one estimation approach might outperform the other. A hybrid-field
channel estimation approach, suitable for AoSA architecture, is introduced in
[28] by combining the orthogonal approximate message passing (OAMP) and deep
learning fixed-point neural networks into a novel estimator. This method
requires significant data training and uses complex networks to achieve a good
performance. Moreover, all of the aforementioned hybrid-field solutions do not
consider the intermediate-field which helps to reduce the estimation
complexity, and do not take into account the prior information of previous
estimations to enhance the estimation performance.

To the best of the authors’ knowledge, this is the first article to solve the
cross-field channel estimation problem using mixed data/model-based solution.
Though [17] mentions the term "cross-field channel estimation", it uses HSPWM
regardless of the distance. This is equivalent to assuming intermediate-field
for all distances and makes [17] akin to other work on HSPWM channel
estimation e.g., [19, 25]. We use SWM in near-field, HSPWM in intermediate-
field, and PWM in far-field. We take a data-driven approach to determine
whether the user is in a near, intermediate, or far-field and subsequently use
a CS-based approach suitable for the determined field. In addition, we propose
a RD method to enhance estimation accuracy and reduce computational complexity
in CS-based channel estimation. The idea of RD for HSPWM was given in the
preliminary version of this work [1]. The prior work [25] also outlined a
similar idea for HSPWM and termed it dictionary shrinkage. In this work, we
extend the RD idea of [1] to SWM and PWM.

###  I-C Organization and Notation

TABLE I: Summary of Notations Notation | Description | Notation | Description
---|---|---|---
$B_{\mathrm{sys}}$ | system bandwidth | $K$ | number of subcarriers
$f_{c}$ | center frequency | $f_{k}$ |  ${k}{\text{th}}$ subcarrier frequency
$Q_{\mathrm{T}}/Q_{\mathrm{R}}$ | number of Tx/Rx SAs | $\bar{Q}_{\mathrm{T}}/\bar{Q}_{\mathrm{R}}$ | number of Tx/Rx AEs per SA
$q_{\mathrm{T}}/q_{\mathrm{R}}$ | index over Tx/Rx SAs | $\bar{q}_{\mathrm{T}}/\bar{q}_{\mathrm{R}}$ | index over Tx/Rx AEs per SA
$N_{\mathrm{T}}/N_{\mathrm{R}}$ | number of UM-MIMO Tx/Rx AEs | $\Delta_{\mathrm{T}}/\Delta_{\mathrm{R}}$ |  | distance between centers of two
---
Tx/Rx adjacent SAs
$\delta_{\mathrm{T}}/\delta_{\mathrm{R}}$ | Tx/Rx AE spacing | $\mathbf{s}[k]$ |  pilot/data symbol at ${k}{\text{th}}$ subcarrier
$\mathbf{x}[k]$ |  transmitted signal at ${k}{\text{th}}$ subcarrier | $\mathbf{y}[k]$ |  received signal at ${k}{\text{th}}$ subcarrier
$\mathbf{F}_{\mathrm{RF}}/\mathbf{W}_{\mathrm{RF}}$ | analog RF beamformer/combiner | $\zeta$ | quantized phase-shifter angle
$N_{\mathrm{S}}$ | number of data streams | $P_{\mathrm{T}}$ | pilot/data transmitted power
$Q_{\mathrm{T}}^{\mathrm{quant}}/Q_{\mathrm{R}}^{\mathrm{quant}}$ |  | number of Tx/Rx phase-shifter
---
quantization bits
$\mathcal{A}$ |  | set of all possible quantized
---
phase-shifter angles
$\mathbf{H}[k]$ |  | UM-MIMO channel response
---
at ${k}{\text{th}}$ subcarrier
$\mathbf{F}_{\mathrm{BB}}[k]/\mathbf{W}_{\mathrm{BB}}[k]$ |  | baseband digital precoder/combiner
---
at ${k}{\text{th}}$ subcarrier
$\mathbf{n}[k]$ |  noise vector at ${k}{\text{th}}$ subcarrier | $\sigma_{n}^{2}$ | noise power
$\dot{\alpha},\dot{\beta},\dot{\gamma}$ |  rotation angles around $\mathrm{ZYX}$-axis | $\theta/\phi$ | angle-of-departure/arrival
$L$ | number of paths | $\alpha^{\ell}$ |  complex gain of the ${\ell}{\text{th}}$ path
$\lambda_{k}$ |  wavelength at ${k}{\text{th}}$ subcarrier | $c_{0}$ | speed of light
$d^{\ell}$ |  length of ${\ell}{\text{th}}$ path | $d^{\mathrm{los}}$ | distance of LoS path
$\mathbf{b}(\cdot,\cdot)$ | near-field array response vector | $\mathbf{a}(\cdot)$ | far-field array response vector
$\mathcal{K}(\cdot)$ | molecular absorption coefficient | $\Gamma^{\ell}$ |  reflection coefficient of ${\ell}{\text{th}}$ path
$\kappa_{\mathrm{T}}$ | refractive index | $\sigma_{\mathrm{rough}}$ | roughness coefficient
$\varphi_{\mathrm{inc}}^{\ell}$ |  angle of incident wave of the ${\ell}{\text{th}}$ path | $\varphi_{\mathrm{ref}}^{\ell}$ |  angle of refracted wave of the ${\ell}{\text{th}}$ path
$\bar{\mathbf{P}}$ | polar-domain near-field dictionary | $\bar{\mathbf{A}}$ | angular-domain far-field dictionary
$G$ |  | number of samples in
---
angular-domain dictionary
$G^{\mathrm{pol}}$ |  | number of samples in
---
polar-domain dictionary
$\mathbf{\Sigma}$ | beamspace representation | $\psi$ | spatial angle
$\bar{\mathbf{\Sigma}}$ | quantized beamspace representation | $\bar{\psi}$ | discrete spatial angle
$M_{\mathrm{T}}/M_{\mathrm{R}}$ | number of Tx/Rx pilots (measurements) | $\eta$ | model selection metric
$\mathbf{C}$ | training combining matrix | $\mathbf{Z}$ | training beamforming matrix
$E$ | number of trials | $R$ | number of rotations
$\gamma_{\mathrm{S}\text{-}\mathrm{H}}$ | near-intermediate fields threshold | $\gamma_{\mathrm{H}\text{-}\mathrm{P}}$ | intermediate-far fields threshold
$\mathbf{\Psi}$ | measurement matrix | $\bar{\mathbf{\Theta}}$ | dictionary matrix
$\mathbf{\Upsilon}$ | sensing matrix | $\hat{L}$ | number of paths to be estimated

The remainder of this paper is organized as follows: We first introduce the
system and channel model in Sec. II. In Sec. III, we formulate the cross-field
channel estimation problem and outline the proposed solution. The simulation
results are presented in Sec. IV, and the conclusion and directions for future
work are given in Sec. V.

We use the following notation throughout the paper. Non-bold lower and upper
case letters $a,A$ denote scalars, bold lower case letters $\mathbf{a}$ denote
vector, and bold upper case letters $\mathbf{A}$ denote matrices.
$[\mathbf{A}]_{n,:}$, and $[\mathbf{A}]_{:,m}$ denote ${n}{\text{th}}$ row of
$\mathbf{A}$ and ${m}{\text{th}}$ column of $\mathbf{A}$, respectively.
$\left|{\cdot}\right|$ shows the absolute value of $a$, absolute value of each
entry of $\mathbf{a}$, and determinant of $\mathbf{A}$.
$\left\|{\mathbf{A}}\right\|_{\mathrm{F}}$ is the Frobenius norm of
$\mathbf{A}$ and $\left\|{\mathbf{a}}\right\|$ is the Euclidean norm of
$\mathbf{a}$. $\mathbf{0}_{N,1}$, $\mathbf{0}_{N,M}$, and $\mathbf{I}_{N}$ are
the zero vector of size $N\times 1$, zero matrix of size $N\times M$, and
identity matrix of size $N\times N$, respectively. The superscripts
${(\cdot)}^{\sf T}$, ${(\cdot)}^{\sf*}$, ${(\cdot)}^{\sf H}$, and
${(\cdot)}^{\sf-1}$ represent the transpose, conjugate, Hermitian (conjugate
transpose), and inverse operators, respectively. $\mathbf{A}\otimes\mathbf{B}$
denotes the Kronecker product of $\mathbf{A}$, and $\mathbf{B}$. $\mbox{\bb
E}[\cdot]$ is the expectation operator, and
$\mathrm{vec}\left(\mathbf{A}\right)$ is the vectorized version of
$\mathbf{A}$ obtained by stacking the column of $\mathbf{A}$.
$\mathrm{diag}\left(a_{1},a_{2},\dots,a_{N}\right)$ is an $N\times N$ diagonal
matrix with diagonal entries $\\{a_{1},a_{2},\dots,a_{N}\\}$,
$\mathrm{blkdiag}\left(\mathbf{a}_{1},\mathbf{a}_{2},\dots,\mathbf{a}_{N}\right)$
denotes the block diagonal matrix with diagonal entries of vectors
$\\{\mathbf{a}_{1},\mathbf{a}_{2},\dots,\mathbf{a}_{N}\\}$. $j=\sqrt{-1}$
denotes the imaginary unit. $\mathcal{CN}(\mathbf{a},\mathbf{A})$ denotes a
complex circularly symmetric Gaussian random vector of mean $\mathbf{a}$ and
covariance $\mathbf{A}$. $\mathcal{U}(a_{1},a_{2})$ is a Uniform distribution
over the interval of $a_{1}$ and $a_{2}$. The subscripts
$\\{\cdot\\}_{\mathrm{T}}$ and $\\{\cdot\\}_{\mathrm{R}}$ denote Tx and Rx
parameters, respectively. The subscript
$\\{\cdot\\}_{(n_{\mathrm{R}},n_{\mathrm{T}})}$ represents the indexing of the
$({n_{\mathrm{T}}}{\text{th}},{n_{\mathrm{R}}}{\text{th}})$ antenna element
(AE) pair within the UM-MIMO system’s dimension. The rest of the notations are
presented in Table I.

##  II System and Channel Models

This section describes the AoSA architecture and the system model. We define
the Tx and Rx signals, the RF analog beamformer/combiner, and the related
hardware constraints. Moreover, we discuss three channel models, i.e., SWM,
PWM, and HSPWM, suitable for near, far, and intermediate-field, respectively.

###  II-A AoSA UM-MIMO System Model and Received Signal

Figure 1: Block diagram of a wideband UM-MIMO THz communication system with an
AoSA architecture.

Consider a downlink multi-carrier [29] THz-band system with an ultra-wide
bandwidth $B_{\mathrm{sys}}$ and UM arrays at the Tx and Rx. The system has
$K$ subcarriers, with the center frequency $f_{c}$ and the frequency of the
${k}{\text{th}}$ subcarrier
$f_{k}=f_{c}+\frac{B_{\mathrm{sys}}}{K}(k-\frac{K-1}{2}),k=\\{1,\cdots,K\\}$.
The system has a sub-connected AoSA hybrid beamforming architecture, in which
each SA contains many AEs, as shown in Fig. 1. We assume that the Tx is
equipped with $Q_{\mathrm{T}}$ SAs and every Tx SA is composed of a uniform
linear array (ULA) consisting of $\bar{Q}_{\mathrm{T}}$ tightly-packed AEs
with inter-element spacing $\delta_{\mathrm{T}}$. Thus, the total number of
AEs at the Tx is $N_{\mathrm{T}}=Q_{\mathrm{T}}\times\bar{Q}_{\mathrm{T}}$.
The centers of two adjacent Tx SAs are separated by a distance
$\Delta_{\mathrm{T}}$. Similar AoSA architecture is used at the Rx, and the
total UM-MIMO system dimension is $N_{\mathrm{R}}\times N_{\mathrm{T}}$.

The discrete-time Tx complex baseband signal at the ${k}{\text{th}}$
subcarrier is expressed as

$\mathbf{x}[k]=\mathbf{F}[k]\mathbf{s}[k]=\mathbf{F}_{\mathrm{RF}}\mathbf{F}_{\mathrm{BB}}[k]\mathbf{s}[k],\vspace{-1mm}$
(1)

where the vector of Tx symbols is
$\mathbf{s}[k]={\left[{s_{1},s_{2},\cdots,s_{N_{\mathrm{S}}}}\right]}^{\sf
T}\in\mbox{\bb C}^{N_{\mathrm{S}}\times 1}$ and $N_{\mathrm{S}}$ is the number
of data streams $(N_{\mathrm{S}}\leq\min\\{Q_{\mathrm{T}},Q_{\mathrm{R}}\\})$.
The Tx symbol vector $\mathbf{s}[k]$ is first precoded by a baseband digital
precoder $\mathbf{F}_{\mathrm{BB}}[k]\in\mbox{\bb C}^{Q_{\mathrm{T}}\times
N_{\mathrm{S}}}$, and subsequently by an analog RF beamformer
$\mathbf{F}_{\mathrm{RF}}\in\mbox{\bb C}^{N_{\mathrm{T}}}\times
Q_{\mathrm{T}}$. We assume frequency-flat analog RF beamformers and get the
analog beamforming block diagonal matrix
$\mathbf{F}_{\mathrm{RF}}=\mathrm{blkdiag}\left(\mathbf{f}_{\mathrm{RF}}^{1},\mathbf{f}_{\mathrm{RF}}^{2},\cdots,\mathbf{f}_{\mathrm{RF}}^{Q_{\mathrm{T}}}\right)$
as follows

$\mathbf{F}_{\mathrm{RF}}=\begin{bmatrix}\mathbf{f}_{\mathrm{RF}}^{1}&\mathbf{0}_{\bar{Q}_{\mathrm{T}}\times
1}&\cdots&\mathbf{0}_{\bar{Q}_{\mathrm{T}}\times 1}\\\
\mathbf{0}_{\bar{Q}_{\mathrm{T}}\times
1}&\mathbf{f}_{\mathrm{RF}}^{2}&\dots&\mathbf{0}_{\bar{Q}_{\mathrm{T}}\times
1}\\\ \vdots&\vdots&\ddots&\vdots\\\ \mathbf{0}_{\bar{Q}_{\mathrm{T}}\times
1}&\mathbf{0}_{\bar{Q}_{\mathrm{T}}\times
1}&\dots&\mathbf{f}_{\mathrm{RF}}^{Q_{\mathrm{T}}}\end{bmatrix},$ (2)

where $\mathbf{f}_{\mathrm{RF}}^{q_{\mathrm{T}}}\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{T}}\times 1},~{}\forall
q_{\mathrm{T}}\in\\{1,\cdots,Q_{\mathrm{T}}\\}$. Note that one RF chain drives
one disjoint SA, where the total number of Tx and Rx RF chains is equal to
$Q_{\mathrm{T}}$ and $Q_{\mathrm{R}}$, respectively. Each AE is connected to a
wideband analog finite-resolution phase-shifter. Therefore, the
${\bar{q}_{\mathrm{T}}}{\text{th}}$ element of the beamforming vector
$\mathbf{f}_{\mathrm{RF}}^{q_{\mathrm{T}}}=\left[f_{\mathrm{RF}}^{q_{\mathrm{T}}}[1],\dots,f_{\mathrm{RF}}^{q_{\mathrm{T}}}[\bar{q}_{\mathrm{T}}],\dots,f_{\mathrm{RF}}^{q_{\mathrm{T}}}[\bar{Q}_{\mathrm{T}}]\right]^{\sf
T}$, i.e., $f_{\mathrm{RF}}^{q_{\mathrm{T}}}[\bar{q}_{\mathrm{T}}]$, has a
constant magnitude and a variable phase, defined by

$f_{\mathrm{RF}}^{q_{\mathrm{T}}}[\bar{q}_{\mathrm{T}}]={\frac{1}{\sqrt{\bar{Q}_{\mathrm{T}}}}}e^{j\zeta_{q_{\mathrm{T}},\bar{q}_{\mathrm{T}}}},\bar{q}_{\mathrm{T}}\in\\{1,\dots,\bar{Q}_{\mathrm{T}}\\}.$
(3)

For a $Q_{\mathrm{T}}^{\mathrm{quant}}$-bit phase-shifter,
$\zeta_{q_{\mathrm{T}},\bar{q}_{\mathrm{T}}}$ has
$2^{Q_{\mathrm{T}}^{\mathrm{quant}}}$ possible values given in the set

$\mathcal{A}=\\{0,\frac{2\pi}{2^{Q_{\mathrm{T}}^{\mathrm{quant}}}},\dots,\frac{2\pi(2^{Q_{\mathrm{T}}^{\mathrm{quant}}}-1)}{2^{Q_{\mathrm{T}}^{\mathrm{quant}}}}\\}.$
(4)

The hybrid precoders satisfy the total power constraint
$\sum_{k=1}^{K}\left(\left\|{\mathbf{F}_{\mathrm{RF}}\mathbf{F}_{\mathrm{BB}}[k]}\right\|_{\mathrm{F}}^{2}\right)=KN_{\mathrm{S}}$.
The Rx side RF and baseband combiners have the same structure and constraints
as the Tx.

The Tx signal passes through the frequency-selective time-invariant THz
channel. Denoting the overall complex UM-MIMO channel matrix at the
${k}{\text{th}}$ subcarrier by $\mathbf{H}[k]\in\mbox{\bb
C}^{N_{\mathrm{R}}\times N_{\mathrm{T}}}$, the channel matrix is

$\mathbf{H}[k]=\begin{bmatrix}\mathbf{H}_{1,1}[k]&\cdots&\mathbf{H}_{1,Q_{\mathrm{T}}}[k]\\\
\vdots&\ddots&\vdots\\\
\mathbf{H}_{Q_{\mathrm{R}},1}[k]&\cdots&\mathbf{H}_{Q_{\mathrm{R}},Q_{\mathrm{T}}}[k]\\\
\end{bmatrix},$ (5)

where $\mathbf{H}_{q_{\mathrm{R}},q_{\mathrm{T}}}[k]\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{R}}\times\bar{Q}_{\mathrm{T}}}$ is the frequency-domain
channel sub-matrix between the ${q_{\mathrm{T}}}{\text{th}}$ Tx SA and the
${q_{\mathrm{R}}}{\text{th}}$ Rx SA. Further details on the channel model will
be discussed in Sec. II-B. By assuming perfect time and frequency
synchronization, the Rx signal at the ${k}{\text{th}}$ subcarrier can be
expressed as

$\mathbf{y}[k]=\mathbf{W}^{\sf H}_{\mathrm{BB}}[k]\mathbf{W}^{\sf
H}_{\mathrm{RF}}\mathbf{H}[k]\mathbf{x}[k]+\mathbf{W}^{\sf
H}_{\mathrm{BB}}[k]\mathbf{W}^{\sf H}_{\mathrm{RF}}\mathbf{n}[k],$ (6)

where
$\mathbf{n}[k]\sim\mathcal{CN}\left(\mathbf{0},\sigma_{n}^{2}\mathbf{I}_{N_{\mathrm{R}}}\right)$
denotes the additive white Gaussian noise (AWGN) vector, and $\sigma_{n}^{2}$
is the noise power. In the next subsection, we describe the used UM-MIMO THz
channel models in detail.

###  II-B Channel Models

Figure 2: Geometry for the UM-MIMO Tx and Rx arrays in the near-field channel
with arbitrary rotation at the Rx.

In practice, both the Tx and Rx are positioned arbitrarily with unknown
orientations. We assume that the Tx location and rotation vector remain fixed.
To represent the array’s 3D orientation, we utilize Euler angles
$\mathbf{o}=[\dot{\alpha},\dot{\beta},\dot{\gamma}]^{\top}$
($\dot{\alpha}\in(-\pi,\pi]$, $\dot{\beta}\in[-\pi/2,\pi/2]$, and
$\dot{\gamma}\in(-\pi,\pi]$) [10]. For modeling purposes, we consider both the
Tx and Rx ULAs aligned along the $\mathrm{Z}$-axis, with propagation occurring
along the $\mathrm{X}$-axis. Since our primary focus is on the relative
orientation between Tx and Rx, we only need to account for rotation around the
$\mathrm{Y}$-axis (represented by $\dot{\beta}$) [10]. The centers of the
Tx/Rx AoSA serve as the coordinate origin. Direction vectors for the angle-of-
departure (AoD) and angle-of-arrival (AoA) are determined based on the
position vectors and are denoted as $\theta\in[0,\pi]$ and $\phi\in[0,\pi]$,
respectively.

####  II-B1 SWM

The SWM is the most accurate channel model for any communication distance and
models the channel between all Tx-Rx AE pairs individually. Typically this
model is used in the near-field. We can define the UM-MIMO channel response of
(5) using SWM, by collecting all channels
$h_{(n_{\mathrm{R}},n_{\mathrm{T}})}^{\mathrm{SWM}}[k]$ from the
${n_{\mathrm{T}}}{\text{th}}$ Tx AE to the ${n_{\mathrm{R}}}{\text{th}}$ Rx AE
at subcarrier $k$ in a matrix $\mathbf{H}^{\mathrm{SWM}}[k]\\!\in\mbox{\bb
C}^{N_{\mathrm{R}}\times N_{\mathrm{T}}}$ as

$\displaystyle\mathbf{H}^{\mathrm{SWM}}[k]$
$\displaystyle=\begin{bmatrix}\mathbf{H}^{\mathrm{SWM}}_{1,1}[k]&\cdots&\mathbf{H}^{\mathrm{SWM}}_{1,Q_{\mathrm{T}}}[k]\\\
\vdots&\ddots&\vdots\\\
\mathbf{H}^{\mathrm{SWM}}_{Q_{\mathrm{R}},1}[k]&\cdots&\mathbf{H}^{\mathrm{SWM}}_{Q_{\mathrm{R}},Q_{\mathrm{T}}}[k]\\\
\end{bmatrix}$ (7)
$\displaystyle=\begin{bmatrix}h_{(1,1)}^{\mathrm{SWM}}[k]&\cdots&h_{(1,N_{\mathrm{T}})}^{\mathrm{SWM}}[k]\\\
\vdots&h_{(n_{\mathrm{R}},n_{\mathrm{T}})}^{\mathrm{SWM}}[k]&\vdots\\\
h_{(N_{\mathrm{R}},1)}^{\mathrm{SWM}}[k]&\cdots&h_{(N_{\mathrm{R}},N_{\mathrm{T}})}^{\mathrm{SWM}}[k]\\\
\end{bmatrix},$

where
$\mathbf{H}^{\mathrm{SWM}}_{q_{\mathrm{R}},q_{\mathrm{T}}}[k]\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{R}}\times\bar{Q}_{\mathrm{T}}}$ is the channel response
between the ${q_{\mathrm{T}}}{\text{th}}$ Tx SA and
${q_{\mathrm{R}}}{\text{th}}$ Rx SA based on SWM and
$h_{(n_{\mathrm{R}},n_{\mathrm{T}})}^{\mathrm{SWM}}[k]$ is defined as

$h_{(n_{\mathrm{R}},n_{\mathrm{T}})}^{\mathrm{SWM}}[k]=\sqrt{\frac{1}{L}}\sum_{\ell=1}^{L}\alpha^{\ell}(f_{k},d^{\ell}_{(n_{\mathrm{R}},n_{\mathrm{T}})})e^{-j\frac{2\pi}{\lambda_{k}}d^{\ell}_{(n_{\mathrm{R}},n_{\mathrm{T}})}},$
(8)

where $L$ is the number of paths, $\alpha^{\ell}$ denotes the complex gain of
the ${\ell}{\text{th}}$ path (further details will be discussed in Sec. II-C),
$\lambda_{k}=\frac{c_{0}}{f_{k}}$ is the wavelength, $c_{0}$ is the speed of
light, and $d^{\ell}_{(n_{\mathrm{R}},n_{\mathrm{T}})}$ is the length of the
${\ell}{\text{th}}$ propagation path between the ${n_{\mathrm{T}}}{\text{th}}$
Tx AE and ${n_{\mathrm{R}}}{\text{th}}$ Rx AE. Here, $\ell=1$ denotes the LoS
path, while $\ell>1$ represents non-LoS paths. The distance
$d^{\mathrm{los}}_{(n_{\mathrm{R}},n_{\mathrm{T}})}$ denotes the communication
distance of the LoS path and $d^{\ell}_{(n_{\mathrm{R}},n_{\mathrm{T}})}$ \-
with $\ell>1$ \- includes the summation of two distances between the Tx-to-
scatterer and scatterer-to-Rx, i.e., single-bounce model. The aforementioned
distances are shown in Fig. 2 in addition to the distance between the
reference, i.e., 1st, AEs at Tx and Rx for the LoS component, denoted as
$d^{\mathrm{los}}_{(1,1)}$.

The array response vector (ARV) $\mathbf{b}(\theta^{\ell},d^{\ell})$ of the
SWM is defined as a function of both distance and angle, i.e., specific
location [22]. This near-field ARV, for a Tx equipped with $\bar{Q}$-AE ULA
and an Rx with single AE, can be expressed as [22]

$\mathbf{b}(\theta^{\ell},d^{\ell})=\frac{1}{\sqrt{\bar{Q}}}\left[e^{-j\frac{2\pi}{\lambda_{k}}(d^{\ell}_{(1,1)}-d^{\ell})},\cdots,e^{-j\frac{2\pi}{\lambda_{k}}(d^{\ell}_{(1,\bar{Q})}-d^{\ell})}\right]^{\sf
T},$ (9)

where
$d^{\ell}_{(1,\bar{q})}=\sqrt{(d^{\ell})^{2}+\bar{\delta}_{\bar{q}}^{2}\delta^{2}-2\bar{\delta}_{\bar{q}}\delta
d^{\ell}\cos{\theta^{\ell}}}$ is the distance of the ${\ell}{\text{th}}$
scatterer from the ${\bar{q}}{\text{th}}$ Tx AE,
$\bar{\delta}_{\bar{q}}=\bar{q}-1-\frac{\bar{Q}-1}{2}(\bar{q}\in\\{1,\cdots,\bar{Q}\\})$
is the AE index, and $d^{\ell}$ denotes the distance of the
${\ell}{\text{th}}$ path (Rx for $\ell=1$ and scatterer for $\ell>1$) from the
Tx center. In this work, we use UM arrays at both Tx and Rx. We can then write
the SWM-based channel response between ${q_{\mathrm{T}}}{\text{th}}$ Tx SA and
${q_{\mathrm{R}}}{\text{th}}$ Rx SA, utilizing the near-field ARVs, as [24]

$\mathbf{H}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{SWM}}[k]\\!=\\!\sqrt{\frac{\bar{Q}_{\mathrm{R}}\bar{Q}_{\mathrm{T}}}{L}}\\!\sum_{\ell=1}^{L}\\!\alpha^{\ell}(\\!f_{k},d^{\ell}_{(n_{\mathrm{R}},n_{\mathrm{T}})})\mathbf{b}_{\mathrm{R}}(\phi^{\ell},d^{\ell}_{\mathrm{R}})\mathbf{b}^{\sf
T}_{\mathrm{T}}(\theta^{\ell},d^{\ell}_{\mathrm{T}}).$ (10)

The model in (10) describes the channel as an outer product of Tx and Rx ARVs.
This description is valid for non-LoS paths, but is only approximately valid
for the LoS path. We refer the interested reader to [24] for a detailed
derivation of (10).

####  II-B2 PWM

If the distance between the Tx and Rx is much larger than the array size, PWM
can be applied. For an AoSA architecture, the PWM implies that the propagation
model is the same across all the Tx and Rx SAs and the distance
$d^{\ell}_{(n_{\mathrm{R}},n_{\mathrm{T}})}$ in (8) is approximated by the
far-field assumptions [19]. Consequently, the UM-MIMO channel
$\mathbf{H}^{\mathrm{PWM}}[k]\in\mbox{\bb C}^{N_{\mathrm{R}}\times
N_{\mathrm{T}}}$ is composed of all sub-channels
$\mathbf{H}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{PWM}}[k]\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{R}}\times\bar{Q}_{\mathrm{T}}}$ between Tx SA
$q_{\mathrm{T}}$ and Rx SA $q_{\mathrm{R}}$ at subcarrier $k$, where the
channel between each Tx and Rx SA is based on PWM, and is given as

$\mathbf{H}^{\mathrm{PWM}}[k]=\begin{bmatrix}\mathbf{H}^{\mathrm{PWM}}_{1,1}[k]&\cdots&\mathbf{H}^{\mathrm{PWM}}_{1,Q_{\mathrm{T}}}[k]\\\
\vdots&\ddots&\vdots\\\
\mathbf{H}^{\mathrm{PWM}}_{Q_{\mathrm{R}},1}[k]&\cdots&\mathbf{H}^{\mathrm{PWM}}_{Q_{\mathrm{R}},Q_{\mathrm{T}}}[k]\\\
\end{bmatrix},\vspace{-1mm}$ (11)

with
$\mathbf{H}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{PWM}}[k]=[h_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{PWM}}[k]]_{\bar{q}_{\mathrm{R}}=1,\cdots,\bar{Q}_{\mathrm{R}},\bar{q}_{\mathrm{T}}=1,\cdots,\bar{Q}_{\mathrm{T}}}$
and the channel response
$[h_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{PWM}}[k]]_{\bar{q}_{\mathrm{R}},\bar{q}_{\mathrm{T}}}$
between the ${\bar{q}_{\mathrm{T}}}{\text{th}}$ Tx AE and
${\bar{q}_{\mathrm{R}}}{\text{th}}$ Rx AE within a SA is defined as

$\displaystyle]$
$\displaystyle{}_{\bar{q}_{\mathrm{R}},\bar{q}_{\mathrm{T}}}=\sqrt{\frac{1}{L}}\sum_{\ell=1}^{L}\alpha^{\ell}(f_{k},d^{\ell}_{(1,1)})e^{-j\frac{2\pi}{\lambda_{k}}d_{(1,1)}^{\ell}}\times$
(12) $\displaystyle
e^{j\frac{2\pi}{\lambda_{k}}\delta_{\mathrm{T}}(\bar{q}_{\mathrm{T}}-1)\cos{\theta_{(1,1)}^{\ell}}}e^{j\frac{2\pi}{\lambda_{k}}\delta_{\mathrm{R}}(\bar{q}_{\mathrm{R}}-1)\cos{\phi_{(1,1)}^{\ell}}},$

where for the ${\ell}{\text{th}}$ path, the channel between two AEs in the
far-field only differs in phase. This phase depends on the angle and the
inter-element spacing. The AoD/AoA are calculated, for all SAs, based on the
reference AE inside the reference, i.e., 1st, SA at the Tx and Rx and the
location of the scatterer, while the path gains are the same across the SAs
following far-field approximations and calculated using the reference
distances $d^{\ell}_{(1,1)}$ (see Fig. 2). For a ULA on $\mathrm{Z}$-axis and
by dropping the subscripts/superscripts, the far-field ARV is

$\mathbf{a}(\theta)=\frac{1}{\sqrt{\bar{Q}}}\left[e^{j\frac{2\pi}{\lambda_{k}}\delta\cos(\theta)\mathbf{\bar{q}}}\right],$
(13)

where $\mathbf{\bar{q}}=\left[0,1,\cdots,\bar{Q}-1\right]^{\sf T}$. Thus, we
can reformulate the sub-channels in (11), using the Tx/Rx ARV
$\mathbf{a}_{\mathrm{T}}(\cdot)\in\mbox{\bb C}^{\bar{Q}_{\mathrm{T}}\times
1}/\mathbf{a}_{\mathrm{R}}(\cdot)\in\mbox{\bb C}^{\bar{Q}_{\mathrm{R}}\times
1}$ of (13), to get

$\displaystyle\mathbf{H}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{PWM}}[k]=\sqrt{\frac{\bar{Q}_{\mathrm{R}}\bar{Q}_{\mathrm{T}}}{L}}\sum_{\ell=1}^{L}$
$\displaystyle\alpha^{\ell}(f_{k},d_{(1,1)}^{\ell})e^{-j\frac{2\pi}{\lambda_{k}}d^{\ell}_{(1,1)}}\times$
(14)
$\displaystyle\mathbf{a}_{\mathrm{R}}\left(\phi_{(1,1)}^{\ell}\right)\mathbf{a}_{\mathrm{T}}^{\sf
T}\left(\theta_{(1,1)}^{\ell}\right).$

####  II-B3 HSPWM

In an UM-MIMO system, a SA usually contains substantially fewer AEs than the
whole array. With the HSPWM channel model, the channel response between each
Tx and Rx SA is based on PWM, whereas the channel variation across SAs is
captured by SWM. Consequently, the MIMO channel between Tx SA $q_{\mathrm{T}}$
and Rx SA $q_{\mathrm{R}}$ at subcarrier $k$,
$\mathbf{H}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{HSPWM}}[k]\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{R}}\times\bar{Q}_{\mathrm{T}}}$, is

$\displaystyle\mathbf{H}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{HSPWM}}[k]=\sqrt{\frac{\bar{Q}_{\mathrm{R}}\bar{Q}_{\mathrm{T}}}{L}}\sum_{\ell=1}^{L}$
$\displaystyle\alpha^{\ell}(f_{k},d_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\ell})e^{-j\frac{2\pi}{\lambda_{k}}d^{\ell}_{q_{\mathrm{R}},q_{\mathrm{T}}}}\times$
(15)
$\displaystyle\mathbf{a}_{\mathrm{R}}\left(\phi^{\ell}_{q_{\mathrm{T}},q_{\mathrm{R}}}\right)\mathbf{a}_{\mathrm{T}}^{\sf
T}\left(\theta^{\ell}_{q_{\mathrm{R}},q_{\mathrm{T}}}\right),$

where AoDs, AoAs, and distances are calculated between the centers of
${q_{\mathrm{T}}}{\text{th}}$ Tx SA and ${q_{\mathrm{R}}}{\text{th}}$ Rx SA
based on the locations of the SAs and scatterers. Moreover, the path gains
$\alpha^{\ell}$ vary across SAs even for the LoS path. HSPWM channel model is
different from the PWM-based AoSA model of (14) where the AoDs and AoAs are
the same for all of Tx-Rx SA channels, and the path gains and the delays are
calculated based on the reference distance. HSPWM channel model (15) is
different from the SWM channel model (10), as the near-field ARV
$\mathbf{b}(\theta,d)$ is replaced by the far-field ARV $\mathbf{a}(\theta)$.

###  II-C THz Specific Channel Characteristics

The path loss for the THz channels consists of two main parts the spreading
loss and molecular absorption loss, both of which vary with carrier frequency
and communication distance [10, 8]. For the LoS path, i.e. $\ell=1$, this loss
is expressed as [10]

$\alpha^{\mathrm{los}}(f_{k},d)=\left(\frac{c_{0}}{4\pi
f_{k}d}\right)^{\frac{\upsilon}{2}}e^{-\frac{1}{2}\mathcal{K}(f_{k})d},$ (16)

where $\upsilon$ is the path loss exponent (its best-fit value is around 2 in
many measurement-based sub-THz campaigns), and $\mathcal{K}(f_{k})$ is the
frequency-dependent molecular absorption coefficient (see [10] for detailed
expressions and calculations).

For the non-LoS paths $(\ell>1)$, we only consider single-bounce reflected
rays, since the first- and second-order reflected paths are attenuated by an
average of $5\\!-\\!10\,\mathrm{dB}$ and more than $15\,\mathrm{dB}$ [10, 21,
8], respectively. Moreover, the scattered and diffracted rays are assumed to
add insignificant contributions to the received signal power, due to the
extreme attenuation with distances [10, 21, 8]. To this end, the reflection
coefficient, $\Gamma^{\ell}\ (\ell>1)$, is given by [21, 8]

$\Gamma^{\ell}=\frac{\cos{\varphi_{\mathrm{inc}}^{\ell}}-\kappa_{\mathrm{T}}\cos{\varphi_{\mathrm{ref}}^{\ell}}}{\cos{\varphi_{\mathrm{inc}}^{\ell}}+\kappa_{\mathrm{T}}\cos{\varphi_{\mathrm{ref}}^{\ell}}}e^{-\left(\frac{8\pi^{2}f_{k}^{2}\sigma_{\mathrm{rough}}^{2}\cos^{2}{\varphi_{\mathrm{inc}}^{\ell}}}{c_{0}^{2}}\right)},$
(17)

where $\kappa_{\mathrm{T}}$ is the refractive index, $\sigma_{\mathrm{rough}}$
models the roughness of the reflecting material,
$\varphi_{\mathrm{inc}}^{\ell}$ is the angle of incident wave of the
${\ell}{\text{th}}$ path, and
$\varphi_{\mathrm{ref}}^{\ell}=\arcsin{\left((\sin{\varphi_{\mathrm{inc}}^{\ell}})/\kappa_{\mathrm{T}}\right)}$
is the angle of refracted wave of the ${\ell}{\text{th}}$ path. Therefore, we
can define the path loss for non-LoS paths as

$\alpha^{\ell}(f_{k},d^{\ell})=\left|{\Gamma^{\ell}}\right|\alpha^{\mathrm{los}}(f_{k},d^{\ell})e^{j\vartheta^{\ell}},$
(18)

where $\vartheta^{\ell}$ represents the random phase shift of the
${\ell}{\text{th}}$ path.

###  II-D Beamspace Representation

We now give the beamspace representation of the THz channels of (10), (14),
and (15), to facilitate the formulation of CS channel estimation problem.

####  II-D1 SWM

Inspired by the fact that both angle and distance information is embedded in
the phase, [22] proposed a polar transform to exploit the channel sparsity in
the near-field. By utilizing the polar-domain matrix $\mathbf{P}\in\mbox{\bb
C}^{\bar{Q}\times L}$, we can write (10) as [22, 24]

$\mathbf{H}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{SWM}}[k]=\mathbf{P}_{\mathrm{R}}\mathbf{\Sigma}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{SWM}}[k]\mathbf{P}_{\mathrm{T}}^{\sf
T},$ (19)

where each column of the Tx/Rx polar-domain matrix $\mathbf{P}$ is a near-
field ARV of (9) and
$\mathbf{\Sigma}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{SWM}}[k]=\sqrt{\frac{\bar{Q}_{\mathrm{R}}\bar{Q}_{\mathrm{T}}}{L}}\mathrm{diag}(\alpha^{\mathrm{los}}(f_{k},d^{\mathrm{los}}_{(n_{\mathrm{R}},n_{\mathrm{T}})}),\cdots,\alpha^{L}(f_{k},d^{L}_{(n_{\mathrm{R}},n_{\mathrm{T}})}))\in\mbox{\bb
C}^{L\times L}$ is the near-field beamspace representation. Due to the
extremely sparse nature of THz channels, CS strategies are appropriate for THz
channel estimation. Since the angle or distance of the channel is not known,
we try to approximate the channel based on pre-determined codebooks and apply
on-grid CS techniques to solve the sparse channel estimation problem. By
linearly sampling the angles and non-linearly sampling the distance, the
polar-domain dictionary $\bar{\mathbf{P}}$ is given as [22, 24]

$\displaystyle\bar{\mathbf{P}}=\Bigg{[}\Bigg{.}$
$\displaystyle\bar{\mathbf{b}}\left(\bar{\theta}_{1},\bar{d}_{(1,1)}^{1}\right),\cdots,\bar{\mathbf{b}}\left(\bar{\theta}_{1},\bar{d}_{(1,1)}^{G^{\mathrm{pol}}_{1}}\right),\cdots,$
(20)
$\displaystyle\bar{\mathbf{b}}\left(\bar{\theta}_{\bar{Q}},\bar{d}^{1}_{(1,\bar{Q})}\right),\cdots,\bar{\mathbf{b}}\left(\bar{\theta}_{\bar{Q}},\bar{d}_{(1,\bar{Q})}^{G^{\mathrm{pol}}_{\bar{Q}}}\right)\Bigg{.}\Bigg{]}.$

The number of samples in the dictionary is
$G^{\mathrm{pol}}=\sum_{\bar{q}=1}^{\bar{Q}}G^{\mathrm{pol}}_{\bar{q}}$, where
$G^{\mathrm{pol}}_{\bar{q}}$ is the number of distance samples for the sampled
angle $\bar{\theta}_{\bar{q}}$. Using (20), we can write (10), by neglecting
the grid quantization error, as

$\mathbf{H}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{SWM}}[k]=\bar{\mathbf{P}}_{\mathrm{R}}\bar{\mathbf{\Sigma}}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{SWM}}[k]\bar{\mathbf{P}}_{\mathrm{T}}^{\sf
T},$ (21)

where $\bar{\mathbf{P}}_{\mathrm{R}}\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{R}}\times
G_{\mathrm{R}}^{\mathrm{pol}}},\bar{\mathbf{P}}_{\mathrm{T}}\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{T}}\times G_{\mathrm{T}}^{\mathrm{pol}}}$, and
$\bar{\mathbf{\Sigma}}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{SWM}}[k]\in\mbox{\bb
C}^{G_{\mathrm{R}}^{\mathrm{pol}}\times G_{\mathrm{T}}^{\mathrm{pol}}}$.

####  II-D2 PWM

By concatenating all of the Tx/Rx far-field ARVs and the channel coefficients
in a matrix form, the PWM channel matrix of (14) can be equivalently written
as [30]

$\mathbf{H}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{PWM}}[k]=\mathbf{A}^{\mathrm{PWM}}_{\mathrm{R}}\mathbf{\Sigma}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{PWM}}[k](\mathbf{A}_{\mathrm{T}}^{\mathrm{PWM}})^{\sf
T},$ (22)

where
$\mathbf{\Sigma}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{PWM}}[k]\in\mbox{\bb
C}^{L\times L}$ is the far-field beamspace representation composed of the
complex path gains defined in (14),
$\mathbf{A}^{\mathrm{PWM}}_{\mathrm{R}}=\left[\mathbf{a}_{\mathrm{R}}\left(\phi_{(1,1)}^{1}\right),\cdots,\mathbf{a}_{\mathrm{R}}\left(\phi_{(1,1)}^{L}\right)\right]\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{R}}\times L}$, and
$\mathbf{A}_{\mathrm{T}}^{\mathrm{PWM}}=\left[\mathbf{a}_{\mathrm{T}}\left(\theta_{(1,1)}^{1}\right),\cdots,\mathbf{a}_{\mathrm{T}}\left(\theta_{(1,1)}^{L}\right)\right]\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{T}}\times L}$.

For far-field PWM CS-based channel estimation, we quantize the angular region
to a grid of size $G$. The ARV of (13), can be defined by using the spatial
angle $\psi\triangleq\frac{\delta}{\lambda_{k}}\cos(\theta)$ as

$\mathbf{a}(\psi)=\frac{1}{\sqrt{\bar{Q}}}\left[e^{j2\pi\psi\mathbf{\bar{q}}}\right].$
(23)

Thus, by defining the discrete spatial angles as
$\bar{\psi}_{\bar{q}}=\frac{1}{\bar{Q}}\left(\bar{q}-\frac{\bar{Q}+1}{2}\right),\bar{q}=\\{1,\cdots,\bar{Q}\\}$,
the matrices $\bar{\mathbf{A}}^{\mathrm{PWM}}_{\mathrm{R}}$ and
$\bar{\mathbf{A}}^{\mathrm{PWM}}_{\mathrm{T}}$, can be expressed using the
far-field dictionary matrix

$\mathbf{\bar{A}}=\left[\mathbf{\bar{a}}(\bar{\psi}_{1}),\mathbf{\bar{a}}(\bar{\psi}_{2}),\cdots,\mathbf{\bar{a}}(\bar{\psi}_{\bar{Q}})\right].$
(24)

Consequently, by neglecting the grid quantization error, we can write (22) as

$\mathbf{H}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{PWM}}[k]=\bar{\mathbf{A}}^{\mathrm{PWM}}_{\mathrm{R}}\bar{\mathbf{\Sigma}}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{PWM}}[k](\bar{\mathbf{A}}_{\mathrm{T}}^{\mathrm{PWM}})^{\sf
T},$ (25)

where
$\bar{\mathbf{\Sigma}}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{PWM}}[k]\in\mbox{\bb
C}^{G_{\mathrm{R}}\times G_{\mathrm{T}}}$ and
$\bar{\mathbf{A}}^{\mathrm{PWM}}_{\mathrm{R}}\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{T}}\times G_{\mathrm{R}}}$ and
$\bar{\mathbf{A}}^{\mathrm{PWM}}_{\mathrm{T}}\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{R}}\times G_{\mathrm{T}}}$ are the angular-domain far-
field dictionaries following the discrete spatial angles of (24) with
$G_{\mathrm{R}}$ and $G_{\mathrm{T}}$ levels, respectively.

####  II-D3 HSPWM

The HSPWM channel of (15) can be written using the same matrix form of (22)
with some changes. The angles of $\mathbf{A}_{\mathrm{R}}$ and
$\mathbf{A}_{\mathrm{T}}$ are changed across the channels between the
${q_{\mathrm{T}}}{\text{th}}$ Tx SA and ${q_{\mathrm{R}}}{\text{th}}$ Rx SA as
well as the complex path gains ${\alpha^{\ell}},\ell\in\\{1,\cdots,L\\}$.
Then, the HSPWM channel matrix of (15) can be expressed as

$\mathbf{H}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{HSPWM}}[k]=\mathbf{A}^{\mathrm{HSPWM}}_{\mathrm{R}}\mathbf{\Sigma}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{HSPWM}}[k](\mathbf{A}_{\mathrm{T}}^{\mathrm{HSPWM}})^{\sf
T},$ (26)

where
$\mathbf{\Sigma}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{HSPWM}}[k]\in\mbox{\bb
C}^{L\times L}$ is the beamspace representation defined following (15),
$\mathbf{A}^{\mathrm{HSPWM}}_{\mathrm{R}}=\left[\mathbf{a}_{\mathrm{R}}\left(\phi^{1}_{q_{\mathrm{T}},q_{\mathrm{R}}}\right),\cdots,\mathbf{a}_{\mathrm{R}}\left(\phi^{L}_{q_{\mathrm{T}},q_{\mathrm{R}}}\right)\right]\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{R}}\times L}$, and
$\mathbf{A}_{\mathrm{T}}^{\mathrm{HSPWM}}=\left[\mathbf{a}_{\mathrm{T}}\left(\theta^{1}_{q_{\mathrm{R}},q_{\mathrm{T}}}\right),\cdots,\mathbf{a}_{\mathrm{T}}\left(\theta^{L}_{q_{\mathrm{R}},q_{\mathrm{T}}}\right)\right]\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{T}}\times L}$. Similar to (25), we can express (26) as

$\mathbf{H}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{HSPWM}}[k]=\bar{\mathbf{A}}^{\mathrm{HSPWM}}_{\mathrm{R}}\bar{\mathbf{\Sigma}}_{q_{\mathrm{R}},q_{\mathrm{T}}}^{\mathrm{HSPWM}}[k](\bar{\mathbf{A}}_{\mathrm{T}}^{\mathrm{HSPWM}})^{\sf
T}.$ (27)

It is important to emphasize that we use the same dictionaries for
intermediate- and far-field channel estimation, i.e.,
$\bar{\mathbf{A}}^{\mathrm{HSPWM}}_{\mathrm{R}}=\bar{\mathbf{A}}^{\mathrm{PWM}}_{\mathrm{R}}$
and
$\bar{\mathbf{A}}_{\mathrm{T}}^{\mathrm{HSPWM}}=\bar{\mathbf{A}}_{\mathrm{T}}^{\mathrm{PWM}}$.
The unknown parameters, i.e., AoDs, AoAs, and path gains, however, are
different across SAs for HSPWM, whereas, they are the same for PWM.

##  III Proposed Strategy for Cross-Field Channel Estimation

As the Tx and Rx can be located in the near, intermediate, or far-field, the
cross-field channel estimation problem arises. Specifically, the appropriate
channel model and estimation strategy need to be selected since the estimation
strategies differ for the three fields. Consequently, we seek a rough estimate
of the distance or an indicator for which model (SWM, HSPWM, or PWM) should be
used for the channel estimation. For this purpose, we introduce a metric
(which will be presented in Sec. III-B1) that - when compared with pre-
determined thresholds - provides this information.

To facilitate the presentation of the overall solution based on CS, we first
outline the pilot transmission using random beamformers/combiners.

###  III-A Pilot Transmission

To begin the training procedure, we transmit pilots using random beamforming
from the reference Tx SA while simultaneously receiving pilots through random
combining at all Rx SAs. If necessary, the pilot transmission is repeated
sequentially through all the $Q_{\mathrm{T}}$ Tx SAs. The pre-computed random
beam weights are stored in the beam codebook and are obtained by uniformly
distributed random sampling of the entries of (4) (see (3) and the following
discussion). In the training phase, the beamforming weight vector
$\mathbf{z}_{m_{\mathrm{T}}}$ is applied on the ${q_{\mathrm{T}}}{\text{th}}$
Tx SA, and the combining weight vector $\mathbf{c}_{m_{\mathrm{R}}}$ is
applied to all Rx SAs. Thus, the received signal on the ${k}{\text{th}}$
subcarrier is

$y_{m_{\mathrm{R}},m_{\mathrm{T}}}^{q_{\mathrm{R}},q_{\mathrm{T}}}[k]=\mathbf{c}^{\sf
H}_{m_{\mathrm{R}}}\mathbf{H}^{\mathrm{SWM}}_{q_{\mathrm{R}},q_{\mathrm{T}}}[k]\mathbf{z}_{m_{\mathrm{T}}}s_{m_{\mathrm{T}}}[k]+\mathbf{c}^{\sf
H}_{m_{\mathrm{R}}}\mathbf{n}_{m_{\mathrm{R}},m_{\mathrm{T}}}[k],$ (28)

where $s_{m_{\mathrm{T}}}[k]$ is the training symbol on subcarrier $k$ and we
assume $\mbox{\bb E}[s_{m_{\mathrm{T}}}[k]s_{m_{\mathrm{T}}}^{\sf
H}[k]]=P_{\mathrm{T}}$, with $P_{\mathrm{T}}$ is the total Tx power used per
transmission during the training phase. The training pilots are transmitted by
applying $M_{\mathrm{T}}$ beamforming vectors sequentially, and for each
beamforming vector $\mathbf{z}_{m_{\mathrm{T}}}$
($m_{\mathrm{T}}=1,\cdots,M_{\mathrm{T}}$), $M_{\mathrm{R}}$ distinct
combining vectors are employed sequentially. Therefore, the number of total
training pilots for the ${q_{\mathrm{T}}}{\text{th}}$ Tx SA is
$M_{\mathrm{R}}\times M_{\mathrm{T}}$, and the resulting received matrix,
after taking away the impact of $s_{m_{\mathrm{T}}}[k]$ using a matched
filter, will be [31]

$\mathbf{Y}^{q_{\mathrm{R}},q_{\mathrm{T}}}[k]=\sqrt{P_{\mathrm{T}}}\mathbf{C}^{\sf
H}\mathbf{H}^{\mathrm{SWM}}_{q_{\mathrm{R}},q_{\mathrm{T}}}[k]\mathbf{Z}+\bar{\mathbf{N}}[k],$
(29)

where
$\mathbf{C}={\left[{\mathbf{c}_{1},\mathbf{c}_{2},\cdots,\mathbf{c}_{M_{\mathrm{R}}}}\right]}^{\sf
T}$ is the $\bar{Q}_{\mathrm{R}}\times M_{\mathrm{R}}$ combining matrix used
at the ${q_{\mathrm{R}}}{\text{th}}$ Rx SA,
$\mathbf{Z}={\left[{\mathbf{z}_{1},\mathbf{z}_{2},\cdots,\mathbf{z}_{M_{\mathrm{T}}}}\right]}^{\sf
T}$ is the $\bar{Q}_{\mathrm{T}}\times M_{\mathrm{T}}$ beamforming matrix used
at the ${q_{\mathrm{T}}}{\text{th}}$ Tx SA, and
$\bar{\mathbf{N}}[k]=\mathbf{C}^{\sf H}\mathbf{N}[k]$ is an
$M_{\mathrm{R}}\times M_{\mathrm{T}}$ noise matrix.

###  III-B Appropriate Channel Model Selection

####  III-B1 Proposed model selection metric

Our objective is to select a suitable channel estimation strategy using an
indicator about the region before starting the channel estimation procedure.
To do this, we can analyze the variation of received signals across the SAs
and make an informed decision about the most suitable channel estimation
approach. Specifically, we expect a lot of variation in received signals
across SAs for near-field, some variation in intermediate-field, and only a
slight - if any - variation in the far-field. The proposed model selection
metric captures the variation across SAs by $\mathrm{(i)}$ collecting the
received signals from all Rx SAs, $\mathrm{(ii)}$ calculating the power of the
difference of received signals for all Rx SA pairs, and $\mathrm{(iii)}$
taking the maximum as an indicator of variation.

Fig. 3 provides an overview of the underlying principles guiding the
development of our proposed model selection metric, denoted as $\eta$. During
the training phase, consider a specific beamforming weight vector, such as
$\mathbf{z}_{m_{\mathrm{T}}}$ with $m_{\mathrm{T}}=1$. It’s important to note
that we employ the same random combining weight vector
$\mathbf{c}_{m_{\mathrm{R}}}$ for a given $m_{\mathrm{R}}$ (Sec. III-A). For
instance, let’s take the case of $m_{\mathrm{R}}=1$. In the near-field, the
reference SA has higher power than the last SA, resulting in a large
difference between the received power of the two. In the far-field, however,
the AoA becomes comparable on all SAs, and the array gain in the direction of
AoA is also comparable giving a small value of $\eta$. In the case of
$m_{\mathrm{R}}=M_{\mathrm{R}}$, the reference SA has a small gain, but the
last SA has a large gain, keeping the difference in power of these two SAs
large. In the far-field, however, the difference still remains small. By
considering all possible Tx and Rx measurements, the proposed metric allows to
effectively designate an indicator for the specific region.

Figure 3: Illustration of the rationale behind the proposed model selection
metric $\eta$.

To formalize, let us collect the vectorized version of the received signal in
(29) for all subcarriers in a vector
$\bar{\boldsymbol{\chi}}^{q_{r},1}=[\mathrm{vec}\left(\mathbf{Y}^{q_{\mathrm{R}},1}[1]\right)^{\sf
T},\cdots,\mathrm{vec}\left(\mathbf{Y}^{q_{\mathrm{R}},1}[K]\right)^{\sf
T}]^{\sf T}$. Subsequently, we obtain a normalized vector
$\boldsymbol{\chi}^{q_{r},1}=\frac{\left|{\bar{\boldsymbol{\chi}}^{q_{r},1}}\right|}{\left\|{\bar{\boldsymbol{\chi}}^{q_{r},1}}\right\|}$
where $r=\\{1,\cdots,Q_{\mathrm{R}}\\}$. Note that we collect the received
signals only from the reference Tx SA as it is sufficient for determining the
region. We now define

$\eta_{r,c}=(\boldsymbol{\chi}^{q_{r},1}-\boldsymbol{\chi}^{q_{c},1})^{\sf
H}(\boldsymbol{\chi}^{q_{r},1}-\boldsymbol{\chi}^{q_{c},1}),$ (30)

where $c=\\{r+1,\cdots,Q_{\mathrm{R}}-1\\}$. Subsequently, we define the
proposed metric $\eta$ as

$\eta=\max{\\{[\eta_{1,2},\\!\cdots\\!,\eta_{r,c},\\!\cdots\\!,\eta_{Q_{\mathrm{R}}-1,Q_{\mathrm{R}}}]_{r=1,\cdots,Q_{\mathrm{R}};c=r+1,\cdots,Q_{\mathrm{R}}-1}\\}}.$
(31)

####  III-B2 Determining the Thresholds for Different Regions

The model selection metric $\eta$ introduced earlier needs to be compared to
pre-computed thresholds to determine the appropriate channel model. The
proposed method to determine the thresholds for model selection metric is
described in Algorithm 1. This procedure is conducted over distances ranging
from near-field to far-field that are stored in a vector
$\mathbf{d}_{\mathrm{offline}}$. For each distance, the metric is calculated
$R_{\mathrm{offline}}\times E_{\mathrm{offline}}$ times.
$R_{\mathrm{offline}}$ represents the number of rotations, and
$E_{\mathrm{offline}}$ represents the number of trials for each rotation.
After that, we collect all $R_{\mathrm{offline}}\times E_{\mathrm{offline}}$
values of $\eta$ for every communication distance and plot the cumulative
distribution function (CDF) to get the thresholds. The threshold
$\gamma_{\mathrm{S}\text{-}\mathrm{H}}$ is used to distinguish between near
and intermediate fields, whereas $\gamma_{\mathrm{H}\text{-}\mathrm{P}}$ is
used to distinguish between intermediate and far fields. Specifically, if
$\eta\geq\gamma_{\mathrm{S}\text{-}\mathrm{H}}$, we select SWM-based channel
estimation, if
$\gamma_{\mathrm{H}\text{-}\mathrm{P}}\leq\eta<\gamma_{\mathrm{S}\text{-}\mathrm{H}}$,
we consider HSPWM-based channel estimation, and if
$\eta<\gamma_{\mathrm{H}\text{-}\mathrm{P}}$, we use PWM-based channel
estimation. These pre-computed thresholds are used subsequently during online
operation.

Algorithm 1 Determining the thresholds

1:Input: Offline distances vector $\mathbf{d}_{\mathrm{offline}}$ of length
$N_{\mathrm{offline}}$, $R_{\mathrm{offline}}$, $E_{\mathrm{offline}}$

2:Output: Thresholds $\gamma_{\mathrm{S}\text{-}\mathrm{H}}$ and
$\gamma_{\mathrm{H}\text{-}\mathrm{P}}$

3:for $\iota=1,\cdots,N_{\mathrm{offline}}$ do

4:for $r=1,\cdots,R_{\mathrm{offline}}$ do

5:for $e=1,\cdots,E_{\mathrm{offline}}$ do

6:Generate $\mathbf{H}^{\mathrm{SWM}}$ over $K$ subcarriers using (7)

7:Select $q_{\mathrm{T}}=1$ and generate random $\mathbf{Z}$ and $\mathbf{C}$

8:for $q_{\mathrm{R}}=1,\cdots,Q_{\mathrm{R}}$ do

9:Compute $\mathbf{Y}^{q_{\mathrm{R}},1}_{\iota,r,e}$ using (29) and store

10:end for

11:Compute $\eta_{\iota,r,e}$ of (31) based on
$\\{\mathbf{Y}^{1,1}_{\iota,r,e},\cdots,\mathbf{Y}^{Q_{\mathrm{R}},1}_{\iota,r,e}\\}$

12:end for

13:end for

14:Compute CDF of $\eta_{\iota,:,:}$

15:end for

16:Determine the thresholds $\gamma_{\mathrm{S}\text{-}\mathrm{H}}$ and
$\gamma_{\mathrm{H}\text{-}\mathrm{P}}$

###  III-C Compressed Channel Estimation using the Selected Model

We vectorize the received pilot signal matrix of (29) to get [31]

$\displaystyle\mathbf{y}^{q_{\mathrm{R}},q_{\mathrm{T}}}[k]$
$\displaystyle=\mathrm{vec}\left(\mathbf{Y}^{q_{\mathrm{R}},q_{\mathrm{T}}}[k]\right)$
(32) $\displaystyle=\sqrt{P_{\mathrm{T}}}\left(\mathbf{Z}^{\sf
T}\otimes\mathbf{C}^{\sf
H}\right)\mathrm{vec}\left(\mathbf{H}^{\mathrm{SWM}}_{q_{\mathrm{R}},q_{\mathrm{T}}}[k]\right)+\mathrm{vec}\left(\bar{\mathbf{N}}[k]\right),$

where $\mathbf{\Psi}=\left(\mathbf{Z}^{\sf T}\otimes\mathbf{C}^{\sf
H}\right)\in\mbox{\bb
C}^{M_{\mathrm{T}}M_{\mathrm{R}}\times\bar{Q}_{\mathrm{T}}\bar{Q}_{\mathrm{R}}}$
is the measurement matrix and
$\mathbf{y}^{q_{\mathrm{R}},q_{\mathrm{T}}}[k]\in\mbox{\bb
C}^{M_{\mathrm{T}}M_{\mathrm{R}}\times 1}$. With the received pilot signal
model (32), we now outline the details of the channel estimation procedure for
SWM, PWM, and HSPWM, respectively.

####  III-C1 SWM-based channel estimation with RD

In the near-field, using (10), (21), and (32), we get

$\displaystyle\mathbf{y}^{q_{\mathrm{R}},q_{\mathrm{T}}}[k]$
$\displaystyle=\sqrt{P_{\mathrm{T}}}\mathbf{\Psi}\left(\bar{\mathbf{P}}_{\mathrm{T}}\otimes\bar{\mathbf{P}}_{\mathrm{R}}\right)\mathrm{vec}\left(\bar{\mathbf{\Sigma}}^{\mathrm{SWM}}_{q_{\mathrm{R}},q_{\mathrm{T}}}[k]\right)+\mathrm{vec}\left(\bar{\mathbf{N}}[k]\right)$
(33)
$\displaystyle=\sqrt{P_{\mathrm{T}}}\mathbf{\Psi}\bar{\mathbf{\Theta}}^{\mathrm{SWM}}\bar{\boldsymbol{\alpha}}^{\mathrm{SWM}}_{(n_{\mathrm{R}},n_{\mathrm{T}})}[k]+\mathrm{vec}\left(\bar{\mathbf{N}}[k]\right)$
$\displaystyle=\sqrt{P_{\mathrm{T}}}\mathbf{\Upsilon}^{\mathrm{SWM}}\bar{\boldsymbol{\alpha}}^{\mathrm{SWM}}_{(n_{\mathrm{R}},n_{\mathrm{T}})}[k]+\mathrm{vec}\left(\bar{\mathbf{N}}[k]\right),$

where
$\bar{\boldsymbol{\alpha}}^{\mathrm{SWM}}_{(n_{\mathrm{R}},n_{\mathrm{T}})}[k]=\mathrm{vec}\left(\bar{\mathbf{\Sigma}}^{\mathrm{SWM}}_{q_{\mathrm{R}},q_{\mathrm{T}}}[k]\right)\in\mbox{\bb
C}^{G^{\mathrm{pol}}_{\mathrm{T}}G^{\mathrm{pol}}_{\mathrm{R}}\times 1}$ is a
sparse vector with only $\bar{L}$ non-zero elements
$(\bar{L}\ll\min\\{G^{\mathrm{pol}}_{\mathrm{R}},G^{\mathrm{pol}}_{\mathrm{T}}\\})$,
the dictionary matrix
$\bar{\mathbf{\Theta}}^{\mathrm{SWM}}=\left(\bar{\mathbf{P}}_{\mathrm{T}}\otimes\bar{\mathbf{P}}_{\mathrm{R}}\right)\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{T}}\bar{Q}_{\mathrm{R}}\times
G^{\mathrm{pol}}_{\mathrm{T}}G^{\mathrm{pol}}_{\mathrm{R}}}$, and hence the
SWM-based sensing matrix is
$\mathbf{\Upsilon}^{\mathrm{SWM}}=\mathbf{\Psi}\bar{\mathbf{\Theta}}^{\mathrm{SWM}}\in\mbox{\bb
C}^{M_{\mathrm{T}}M_{\mathrm{R}}\times
G^{\mathrm{pol}}_{\mathrm{T}}G^{\mathrm{pol}}_{\mathrm{R}}}$. Moreover, the CS
problem of (33) needs to be solved for each subcarrier $k$ using some CS tool,
e.g., OMP [32]. Alternatively, by exploiting the common support property [33],
a single CS problem can be solved across all subcarriers, e.g., using SOMP
[33]. The CS problem (33) is solved for all Tx-Rx SA pairs.

We propose a RD method for channel estimation by exploiting the geometric
structure in the problem and using the prior estimation information to reduce
the complexity and enhance the estimation accuracy. The key idea of RD is to
reduce the search space around an estimate of the angles/distances obtained
from an initial search. Specifically, only $G_{\mathrm{RD}}$ columns of the
full-size Tx/Rx dictionaries are kept and searched over. We use the full-size
dictionary for the reference Tx-Rx SA and RD for the remaining Rx SAs.

Fig. 4 provides a visual representation of the underlying concept behind the
proposed SWM-RD method. The complete search space in the polar domain, for the
near-field, is a two-dimensional space, involving non-uniform samples across
the distance domain and uniform samples of the angle within the range of
$[0^{\circ},180^{\circ}]$ in the angular domain. The ground truth polar-domain
representation of each Rx SA is uniquely marked with a star symbol in
different colors. During the first Tx-Rx SA channel estimation, which involves
no prior estimation information, the entire search region will be used which
spans a wide range of angles and distances. Subsequently, as we reduce the
search space, both the angles and distances are trimmed down. The distances
and angles relevant to the ${q_{\mathrm{R}}}{\text{th}}$ Rx SA are depicted in
purple, while those for the $Q_{\mathrm{R}}$ Rx SA are represented in orange.

Figure 4: A visual explanation of the proposed reduced dictionary (RD)
combined with the SWM estimation method to result in SWM-RD in the near-field
case where for the reference Tx-Rx SA the complete search space is used at the
Rx SA. For the ${q_{\mathrm{R}}}{\text{th}}$ Rx SA, a reduced search space is
used, where the reduction is based on the ${(q_{\mathrm{R}}-1)}{\text{th}}$
SA’s estimation. The search space for Rx SA $Q_{\mathrm{R}}$ is also reduced,
but the reduced search space is different since it is based on the channel
estimate from SA number $Q_{\mathrm{R}}-1$.

To do this, we start by solving the CS problem (33) for the reference Tx-Rx
SAs, e.g., by using SOMP algorithm and set the number of paths to be estimated
as $\hat{L}$. The algorithm’s output will contain the estimated support. We
then generate the reduced dictionaries using the procedure summarized in
Algorithm 2. The input dictionary contains $\hat{L}$ columns selected from the
dictionary $\bar{\mathbf{P}}_{\mathrm{T}}$ (or
$\bar{\mathbf{P}}_{\mathrm{R}}$) based on the support of the previously
estimated channel. Recall that $\mathcal{L}$ is the index set associated with
the input dictionary. Hence $[\bar{\mathbf{P}}]_{:,\mathcal{L}}$ are the
desired dictionary columns to limit the search space. We construct an
oversampled - by a factor of $G_{\mathrm{ovs}}$ \- polar-domain dictionary
$\bar{\mathbf{P}}_{\mathrm{ovs}}\in\mbox{\bb C}^{\bar{Q}\times
G_{\mathrm{ovs}}G^{\mathrm{pol}}}$ according to Algorithm 1 in [22]. The
desired RD then consists of the $G^{\mathrm{pol}}_{\mathrm{RD}}$ columns
$(G^{\mathrm{pol}}_{\mathrm{RD}}\leq G^{\mathrm{pol}})$ from the oversampled
dictionary that have the highest correlation with
$[\bar{\mathbf{P}}]_{:,\mathcal{L}}$, and corresponds to specific locations
(angles and distances) centered around the estimated location. Subsequently,
for each remaining Rx SAs ($2,\cdots,Q_{\mathrm{R}}$), we solve (33) by using
$\bar{\mathbf{\Theta}}_{\mathrm{RD}}^{\mathrm{SWM}}=\left(\bar{\mathbf{P}}_{\mathrm{T},\mathrm{RD}}\otimes\bar{\mathbf{P}}_{\mathrm{R},\mathrm{RD}}\right)\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{T}}\bar{Q}_{\mathrm{R}}\times
G^{\mathrm{pol}}_{\mathrm{T},\mathrm{RD}}G^{\mathrm{pol}}_{\mathrm{R},\mathrm{RD}}}$
which is the RD matrix obtained from the support of the previous channel
estimation. The proposed RD method is repeated for all $Q_{\mathrm{T}}$ Tx
SAs.

Algorithm 2 Reduced Dictionary Method

1:Input: $\mathcal{L}$, $\mathcal{G}$, $\bar{\mathbf{A}}$

2:Output: $\bar{\mathbf{A}}_{\mathrm{RD}}$

3:Construct an oversampled (by a factor of $G_{\mathrm{ovs}}$) dictionary
$\bar{\mathbf{A}}_{\mathrm{ovs}}$.

4:Let
$\mathbf{U}=\bar{\mathbf{A}}_{\mathrm{ovs}}^{\sf*}[\bar{\mathbf{A}}]_{:,\mathcal{L}}$.
Populate the index set $\mathcal{G}$ with the indices of $G_{\mathrm{RD}}$
columns of $\mathbf{U}$ that have the largest value.

5:Create the RD matrix
$\bar{\mathbf{A}}_{\mathrm{RD}}=[\bar{\mathbf{A}}_{\mathrm{ovs}}]_{:,\mathcal{G}}$.

####  III-C2 PWM-based channel estimation with RD

Using the beamspace representation of the PWM channel (25), (32) can be
written as

$\displaystyle\mathbf{y}^{q_{\mathrm{R}},q_{\mathrm{T}}}[k]$
$\displaystyle=\sqrt{P_{\mathrm{T}}}\mathbf{\Psi}\bar{\mathbf{\Theta}}^{\mathrm{PWM}}\bar{\boldsymbol{\alpha}}^{\mathrm{PWM}}_{(1,1)}[k]+\mathrm{vec}\left(\bar{\mathbf{N}}[k]\right)$
(34)
$\displaystyle=\sqrt{P_{\mathrm{T}}}\mathbf{\Upsilon}^{\mathrm{PWM}}\bar{\boldsymbol{\alpha}}^{\mathrm{PWM}}_{(1,1)}[k]+\mathrm{vec}\left(\bar{\mathbf{N}}[k]\right),$

where
$\bar{\boldsymbol{\alpha}}^{\mathrm{PWM}}_{(1,1)}[k]=\mathrm{vec}\left(\bar{\mathbf{\Sigma}}^{\mathrm{PWM}}_{q_{\mathrm{R}},q_{\mathrm{T}}}[k]\right)\in\mbox{\bb
C}^{G_{\mathrm{T}}G_{\mathrm{R}}\times 1}$,
$\bar{\mathbf{\Theta}}^{\mathrm{PWM}}=\left(\bar{\mathbf{A}}_{\mathrm{T}}^{\mathrm{PWM}}\otimes\bar{\mathbf{A}}_{\mathrm{R}}^{\mathrm{PWM}}\right)\in\mbox{\bb
C}^{\bar{Q}_{\mathrm{T}}\bar{Q}_{\mathrm{R}}\times
G_{\mathrm{T}}G_{\mathrm{R}}}$, and
$\mathbf{\Upsilon}^{\mathrm{PWM}}\in\mbox{\bb
C}^{M_{\mathrm{T}}M_{\mathrm{R}}\times G_{\mathrm{T}}G_{\mathrm{R}}}$ is the
PWM-based sensing matrix. The greedy search - used in SOMP algorithm - depends
on the number of the sensing matrix columns. Hence, the PWM-based CS
estimation is expected to be less complex than SWM-based since
$G_{\mathrm{T}}G_{\mathrm{R}}\ll
G_{\mathrm{T}}^{\mathrm{pol}}G_{\mathrm{R}}^{\mathrm{pol}}$. The additional
grid points, due to distance sampling, while constructing the polar-domain
codebook of (20) results in this increased complexity.

The conventional PWM estimation involves solving the problem defined in (34)
separately for each Tx-Rx SA pair. However, drawing inspiration from the
analysis of the far-field PWM channel model and the corresponding beamspace
representation, it becomes evident that the channels across different SAs are
the same. Thus, we introduce a more efficient approach by estimating only the
reference Tx-Rx SA channel, which we refer to as PWM-RD. To illustrate this
concept, Fig. 5 presents a visual representation where it can be observed that
all of the Rx SAs’ beamspace representations fall within the same beam in the
far-field angular domain. As a result, we perform the estimation only once for
the reference Tx-Rx SA channel, as demonstrated in Fig. 5. The complete search
space, in this case, is within the range of $[0^{\circ},180^{\circ}]$
(depicted in gold). For the remaining Rx SAs, there is no need for a search,
as we can utilize the beams obtained from the reference Tx-Rx SA estimation.
PWM-RD has low computational complexity as only the channel between reference
SAs is estimated. It also has low pilot overhead since pilots are transmitted
only from the reference Tx SA.

Figure 5: A visual representation of the proposed reduced dictionary (RD)
method in the far-field case, referred to as PWM-RD.

####  III-C3 HSPWM-based channel estimation with RD

Using beamspace representation (27), the Rx signal (32) in the intermediate
field can be written as

$\displaystyle\mathbf{y}^{q_{\mathrm{R}},q_{\mathrm{T}}}[k]$
$\displaystyle=\sqrt{P_{\mathrm{T}}}\mathbf{\Psi}\bar{\mathbf{\Theta}}^{\mathrm{HSPWM}}\bar{\boldsymbol{\alpha}}^{\mathrm{HSPWM}}_{q_{\mathrm{R}},q_{\mathrm{T}}}[k]+\mathrm{vec}\left(\bar{\mathbf{N}}[k]\right)$
(35)
$\displaystyle=\sqrt{P_{\mathrm{T}}}\mathbf{\Upsilon}^{\mathrm{HSPWM}}\bar{\boldsymbol{\alpha}}^{\mathrm{HSPWM}}_{q_{\mathrm{R}},q_{\mathrm{T}}}[k]+\mathrm{vec}\left(\bar{\mathbf{N}}[k]\right).$

In HSPWM, the angles and distances vary across the SAs, which constitutes the
primary distinction in channel estimation between HSPWM and PWM. This
difference is manifested in the fact that the best beam for each channel
corresponds to different angular beams in Fig.6. Furthermore, for HSPWM-RD,
the complete search space aligns with that of PWM, denoted in gold. However,
the reduced search space narrows down to specific angles for the
${q_{\mathrm{R}}}{\text{th}}$ Rx SA, depicted in purple, and the Rx SA
$Q_{\mathrm{R}}$, highlighted in orange.

Figure 6: A visual explanation of the proposed reduced dictionary (RD) method
in the intermediate-field case, referred to as HSPWM-RD.

For HSPWM-RD, we use the same procedure to get the RD as for the SWM case but
there are four differences: $(\mathrm{i})$ we extract the indices of the
dominant beamspace elements, referred to as $L_{\mathrm{dom}}$, from the
support of the estimated channel. These indices correspond to the estimated
channel elements responsible for $95\%$ of the total power. We then modify the
index set $\mathcal{L}$ to exclusively include these $L_{\mathrm{dom}}$
elements instead of all $\hat{L}$ indices, $(\mathrm{ii})$ In scenarios
involving both LoS and non-LoS paths, we anticipate that the index set
$\mathcal{G}$ will comprise many different clusters centered around the beams
corresponding to each path component in the beamspace. Therefore, after
obtaining $U$ in step 2 of Algorithm 2, we detect each set of consecutive
indices and group them into a cluster
$\mathcal{B}^{\mathrm{RD}}_{n_{\mathrm{clu}}}$, where $n_{\mathrm{clu}}$
ranges from 1 to $N_{\mathrm{clu}}$, representing the number of detected
clusters. We calculate the center and the number of indices,
$L_{n_{\mathrm{clu}}}$, for each cluster. The center of each cluster
corresponds to the central point of the corresponding reduced space search,
and the number of columns for this cluster is determined by
$\lfloor\frac{L_{n_{\mathrm{clu}}}}{L_{\mathrm{dom}}}G_{\mathrm{RD}}\rfloor$.
We keep only the unique columns, preventing redundant searches over the same
columns in the case of overlapping clusters. Furthermore, for a LoS-dominant
channel, this approach centers the reduced search space around the first beam
(the best beam-pair) of the unique cluster we obtained, with the reduction
focused on specific spatial angles centered around the best estimated AoD/AoA,
$\mathrm{(iii)}$ the oversampled dictionary $\bar{\mathbf{A}}_{\mathrm{ovs}}$
is obtained according to (24) with $G_{\mathrm{ovs}}\times G$ levels, and
$(\mathrm{iv})$ the resultant reduced dictionaries for HSPWM are not
necessarily equal in size to the reduced dictionaries for SWM. This is
because: $(\mathrm{a})$ the polar dictionaries used for SWM have a large size,
due to sampling in both angles and distances, compared to angular dictionaries
used for HSPWM, and $(\mathrm{b})$ further, the reduction factors for HSPWM
and SWM can also be different.

###  III-D Overall Channel Estimation Algorithm

Algorithm 3 Cross-Field RD-based CS Channel Estimation

1:Input: Number of paths $\hat{L}$, thresholds
$\gamma_{\mathrm{S}\text{-}\mathrm{H}}$ and
$\gamma_{\mathrm{H}\text{-}\mathrm{P}}$, sensing matrices
$\mathbf{\Upsilon}^{\mathrm{SWM}}$, $\mathbf{\Upsilon}^{\mathrm{HSPWM}}$, and
$\mathbf{\Upsilon}^{\mathrm{PWM}}$

2:Output: Estimated UM-MIMO channel $\hat{\mathbf{H}}[k]$ for all subcarriers

3:Pilot transmission from reference Tx SA

4:Compute $\eta_{\mathrm{online}}$ using (31) based on the received signals
$\mathbf{y}^{q_{\mathrm{R}},q_{\mathrm{T}}}[k]$ in (32) for all Rx SAs

5:if $\eta_{\mathrm{online}}\geq\gamma_{\mathrm{S}\text{-}\mathrm{H}}$ then
start SWM estimation

6:if $q_{\mathrm{R}}=1$

7:Perform SWM-based channel estimation for
$\mathbf{y}^{q_{\mathrm{R}},q_{\mathrm{T}}}[k]$ in (32) using SOMP algorithm
(based on $\hat{L}$ and $\mathbf{\Upsilon}^{\mathrm{SWM}}$)

8:else for the remaining Rx SAs
$q_{\mathrm{R}}\in\\{2,\cdots,Q_{\mathrm{R}}\\}$

9:Perform SWM-based channel estimation for
$\mathbf{y}^{q_{\mathrm{R}},q_{\mathrm{T}}}[k]$ in (32) using SOMP-RD
algorithm (based on $\hat{L}$ and
$\mathbf{\Upsilon}_{\mathrm{RD}}^{\mathrm{SWM}}$)

10:end if

11:Send a one-time flag to the Tx to continue the pilot transmission

12:Repeat the steps (4 $\rightarrow$ 8) for the remaining $(Q_{\mathrm{T}}-1)$
Tx SAs.

13:else if
$\gamma_{\mathrm{H}\text{-}\mathrm{P}}\leq\eta_{\mathrm{online}}<\gamma_{\mathrm{S}\text{-}\mathrm{H}}$
then start HSPWM estimation

14:if $q_{\mathrm{R}}=1$

15:Perform HSPWM-based channel estimation for
$\mathbf{y}^{q_{\mathrm{R}},q_{\mathrm{T}}}[k]$ in (32) using SOMP algorithm
(based on $\hat{L}$ and $\mathbf{\Upsilon}^{\mathrm{HSPWM}}$)

16:else for the remaining Rx SAs
$q_{\mathrm{R}}\in\\{2,\cdots,Q_{\mathrm{R}}\\}$

17:Perform HSPWM-based channel estimation for
$\mathbf{y}^{q_{\mathrm{R}},q_{\mathrm{T}}}[k]$ in (32) using SOMP-RD
algorithm (based on $\hat{L}$ and
$\mathbf{\Upsilon}_{\mathrm{RD}}^{\mathrm{HSPWM}}$)

18:end if

19:Send a one-time flag to the Tx to continue the pilot transmission

20:Repeat the steps (12 $\rightarrow$ 16) for the remaining
$(Q_{\mathrm{T}}-1)$ Tx SAs.

21:else start PWM-RD estimation

22:Perform PWM-RD channel estimation for
$\mathbf{y}^{q_{\mathrm{R}},q_{\mathrm{T}}}[k]$ in (32) using SOMP algorithm
(based on $\hat{L}$ and $\mathbf{\Upsilon}^{\mathrm{PWM}}$) for the reference
channel

23:Send a one-time flag to the Tx to stop the pilot transmission

24:end if

The complete proposed channel estimation algorithm is given in Algorithm 3. By
comparing the value of the proposed model selection metric
$\eta_{\mathrm{online}}$ to the thresholds
$\gamma_{\mathrm{S}\text{-}\mathrm{H}}$ and
$\gamma_{\mathrm{H}\text{-}\mathrm{P}}$, we can determine whether SWM, HSPWM,
or PWM based estimation should be used. Feedback to the Tx is required to
continue or stop the pilot transmission after transmission from the reference
Tx SA.

After estimating the channel, detecting the largest absolute entry of (33)
((34) or (35)) determines the best beam-pair [31]. The Rx needs to feed the
estimated index of the transmit beam for each Tx SA back to the Tx. Both Tx
and Rx use the estimated Tx/Rx beam indices to configure the analog
beamformers $\hat{\mathbf{F}}_{\mathrm{RF}}$ and combiners
$\hat{\mathbf{W}}_{\mathrm{RF}}$, respectively. Furthermore,
$\tilde{\mathbf{A}}^{\mathrm{PWM}}$, $\tilde{\mathbf{A}}^{\mathrm{HSPWM}}$,
and $\tilde{\mathbf{P}}$ are used for RF beamforming/combining and are
obtained by replacing the phase of each entry in
$\bar{\mathbf{A}}^{\mathrm{PWM}}$, $\bar{\mathbf{A}}^{\mathrm{HSPWM}}$, and
$\bar{\mathbf{P}}$ by their closest phase in the quantized phase set due to
finite-resolution phase-shifter.

The signals from multiple SAs are combined in the baseband after down
conversion. For the baseband-to-baseband equivalent channel
$\hat{\mathbf{W}}_{\mathrm{RF}}^{\sf
H}\mathbf{H}[k]\hat{\mathbf{F}}_{\mathrm{RF}}$, the singular value
decomposition (SVD)-based precoders $\hat{\mathbf{F}}_{\mathrm{BB}}[k]$ and
combiners $\hat{\mathbf{W}}_{\mathrm{BB}}[k]$ can be obtained [34], using only
the first $N_{\mathrm{S}}$ columns of the decomposition.

###  III-E Complexity Analysis

The computational complexity of Algorithm 3 is dominated by SOMP and SOMP-RD
complexity. Specifically, for SOMP adopting angular-domain dictionaries, for
one Tx-Rx SA channel, the complexity can be described as
$\mathcal{C}_{\mathrm{SOMP}}=\mathcal{O}\left(G_{\mathrm{R}}G_{\mathrm{T}}KM_{\mathrm{T}}M_{\mathrm{R}}\hat{L}\right)$
[33]. In cases where polar-domain dictionaries are utilized, the complexity
takes on the form of
$\mathcal{C}_{\mathrm{SOMP}}^{\mathrm{pol}}=\mathcal{O}\left(G_{\mathrm{R}}^{\mathrm{pol}}G_{\mathrm{T}}^{\mathrm{pol}}KM_{\mathrm{T}}M_{\mathrm{R}}\hat{L}\right)$
[22]. Consequently, for an AoSA architecture, the complexity of SOMP for SWM,
HSPWM, PWM, and their RD counterparts can be generally written as

$\mathcal{C}=\mathcal{O}\left(N_{\mathrm{est}}KM_{\mathrm{T}}M_{\mathrm{R}}\hat{L}\right),$
(36)

where $N_{\mathrm{est}}$ differs between different strategies, and is given
below

  * •

SWM:
$N_{\mathrm{est}}=Q_{\mathrm{R}}Q_{\mathrm{T}}G_{\mathrm{R}}^{\mathrm{pol}}G_{\mathrm{T}}^{\mathrm{pol}}$
equals the number of Tx and Rx SAs multiplied by the Tx and Rx polar
dictionaries. This complexity arises from our use of SOMP and polar
dictionaries for the estimation of all Tx-Rx SA channels.

  * •

SWM-RD:
$N_{\mathrm{est}}=Q_{\mathrm{T}}G_{\mathrm{R}}^{\mathrm{pol}}G_{\mathrm{T}}^{\mathrm{pol}}+(Q_{\mathrm{R}}-1)Q_{\mathrm{T}}G_{\mathrm{R},\mathrm{RD}}^{\mathrm{pol}}G_{\mathrm{T},\mathrm{RD}}^{\mathrm{pol}}$.
The first term is because we use the SOMP with the full-size SWM sensing
matrix for estimating the first Tx-Rx SA channel. In contrast, for the
remaining Rx SAs, we employ the RD approach, which employs a SWM sensing
matrix with fewer columns, as detailed in Sec. III-C1. Since
$G_{\mathrm{RD}}^{\mathrm{pol}}\leq G^{\mathrm{pol}}$, the complexity is
reduced compared to SWM.

  * •

HSPWM:
$N_{\mathrm{est}}=Q_{\mathrm{R}}Q_{\mathrm{T}}G_{\mathrm{R}}G_{\mathrm{T}}$
since we use SOMP combined with angular-domain dictionaries for estimating all
Tx-Rx SA channels.

  * •

HSPWM-RD:
$N_{\mathrm{est}}=Q_{\mathrm{T}}G_{\mathrm{R}}G_{\mathrm{T}}+(Q_{\mathrm{R}}-1)Q_{\mathrm{T}}G_{\mathrm{R},\mathrm{RD}}G_{\mathrm{T},\mathrm{RD}}$.
In this case, the two terms are the same as SWM-RD but based on the angular-
domain codebooks instead of the polar-domain. Since $G_{\mathrm{RD}}\leq G$,
the complexity is reduced compared to HSPWM.

  * •

PWM: the complexity is equal to the HSPWM case.

  * •

PWM-RD: $N_{\mathrm{est}}=G_{\mathrm{R}}G_{\mathrm{T}}$ since the estimation
is performed only once.

The complexity analysis provided in (36) reveals that PWM-RD estimation has
the lowest complexity. SWM estimation complexity is the highest since
additional grid points are used in the polar codebooks. Finally, the
complexity of the proposed cross-field RD channel estimation will be a
weighted sum of three cases, i.e., SWM-RD, HSPWM-RD, and PWM-RD, where the
weights correspond to the fraction of times SWM-RD, HSPWM-RD, or PWM-RD is
utilized.

An additional component is added for the SWM-RD complexity is equal to
$\mathcal{C}_{\mathrm{SWM}}^{\mathrm{RD}}=\mathcal{C}_{\mathrm{SWM}}^{\mathrm{T},\mathrm{RD}}+\mathcal{C}_{\mathrm{SWM}}^{\mathrm{R},\mathrm{RD}}=\mathcal{O}\left(G_{\mathrm{ovs}}G_{\mathrm{T}}^{\mathrm{pol}}\hat{L}\left(Q_{\mathrm{R}}-1\right)Q_{\mathrm{T}}\right)+\mathcal{O}\left(G_{\mathrm{ovs}}G_{\mathrm{R}}^{\mathrm{pol}}\hat{L}\left(Q_{\mathrm{R}}-1\right)Q_{\mathrm{T}}\right)$.
These two terms are incorporated to account for the complexity of obtaining
the Tx and Rx reduced dictionaries. These terms are primarily dominated by
step 2 in Algorithm 2. This same analysis extends to HSPWM-RD, resulting in
$\mathcal{C}_{\mathrm{HSPWM}}^{\mathrm{RD}}=\mathcal{O}\left(G_{\mathrm{ovs}}G_{\mathrm{T}}L_{\mathrm{dom}}\left(Q_{\mathrm{R}}-1\right)Q_{\mathrm{T}}\right)+\mathcal{O}\left(G_{\mathrm{ovs}}G_{\mathrm{R}}L_{\mathrm{dom}}\left(Q_{\mathrm{R}}-1\right)Q_{\mathrm{T}}\right)$.

##  IV Simulation Results and Discussion

In this section, we present simulation results and numerical analysis of the
proposed strategy.

###  IV-A Performance Analysis Metrics

We use the following metrics in our evaluations and performance analysis.

####  IV-A1 Normalized Mean Square Error

The first metric used for performance comparison is the NMSE, defined as

$\mathrm{NMSE}=\frac{1}{ER}\sum_{e=1}^{E}\sum_{r=1}^{R}\frac{\sum_{k=1}^{K}\left\|{\hat{\mathbf{H}}_{e,r}[k]-\mathbf{H}_{e,r}[k]}\right\|_{\mathrm{F}}^{2}}{\sum_{k=1}^{K}\left\|{\mathbf{H}_{e,r}[k]}\right\|_{\mathrm{F}}^{2}},$
(37)

where $\hat{\mathbf{H}}[k]$ is the estimated channel at subcarrier $k$, $E$ is
the number of simulation trials conducted for ensemble averaging, and $R$ is
the number of rotations.

####  IV-A2 Effective Achievable Rate

The second metric is the effective achievable rate (EAR),
$\mathcal{R}_{\mathrm{eff}}$, which is defined as

$\mathcal{R}_{\mathrm{eff}}=\frac{\varrho}{ERK}\sum_{e=1}^{E}\sum_{r=1}^{R}\sum_{k=1}^{K}\log_{2}\left|{\mathbf{I}_{Q_{\mathrm{R}}}+\frac{P_{\mathrm{T}}}{KN_{\mathrm{S}}}\mathbf{R}_{n}^{\sf-1}\mathbf{C}_{e,r}^{\sf
H}[k]\mathbf{C}_{e,r}[k]}\right|,$ (38)

where $\mathbf{C}[k]=\hat{\mathbf{F}}^{\sf H}[k]\mathbf{H}^{\sf
H}[k]\hat{\mathbf{W}}[k]$,
$\hat{\mathbf{F}}[k]=\hat{\mathbf{F}}_{\mathrm{RF}}\hat{\mathbf{F}}_{\mathrm{BB}}[k]$
($\hat{\mathbf{W}}[k]$ is defined in a similar way), and
$\mathbf{R}_{n}=\sigma_{n}^{2}\hat{\mathbf{W}}^{\sf H}[k]\hat{\mathbf{W}}[k]$,
$\varrho\\!\triangleq\\!\max\left(0,1-\frac{M_{\mathrm{T}}^{\mathrm{tr}}M_{\mathrm{R}}^{\mathrm{tr}}}{T_{\mathrm{coh}}}\right)$
where $M_{\mathrm{T}}^{\mathrm{tr}}/M_{\mathrm{R}}^{\mathrm{tr}}$ is the
number of used Tx/Rx pilots during the training (pilot overhead), and
$T_{\mathrm{coh}}$ is the THz channel coherence time in symbols. Furthermore,
$\varrho$ captures the loss in achievable rate due to the training [31]. Note
that by setting $\varrho=1$, we get the AR $\mathcal{R}$.

####  IV-A3 Computational Complexity

We evaluate the computational complexity numerically based on two indicators
the complexity analysis of Sec. III-E and the wall-clock running time.

####  IV-A4 Benchmark Reference Work

We have adopted four benchmark methods for comparison with our proposed work:

  * •

SOMP NF-polar [22]: The SOMP algorithm is combined with polar-domain codebooks
for all Tx-Rx SA channel estimations regardless of the distances. The
complexity of this method is equal to the SWM case in Sec. III-E. We have also
introduced a modified version with reduced complexity, achieved by employing
only $\hat{L}/2$ SOMP iterations. This modified implementation is referred to
as "SOMP NF-polar $\hat{L}/2$" in our simulation results.

  * •

SOMP FF-angular [33]: The polar-domain codebooks used in SOMP NF-Polar [22]
are replaced by angular-domain codebooks. The complexity is similar to HSPWM
in Sec. III-E.

  * •

Hybrid-field FF-NF [26]: This is the same implementation as Algorithm 1 in
[26] for all distances. However, we replaced the OMP with SOMP to be
compatible with our work. The number of far-field and near-field SOMP
iterations is $\hat{L}/2$ for each. Hence, the complexity of this algorithm is
equal to the average of SOMP NF-polar and SOMP FF-angular.

  * •

Hybrid-field NF-FF [26]: For this implementation, we reversed the order of
applying the estimation, i.e., we start with near-field SOMP and polar
codebooks and then use the far-field part. The complexity remains the same as
Hybrid-Field FF-NF.

  * •

Mixed-field NF [24]: The LoS component is estimated using Algorithm 1 in [24]
and the non-LoS paths are estimated using Algorithm 2 in [24] which is the
same as SOMP NF-polar. The complexity is equal to
$\mathcal{O}\left(\left(\left(S_{\mathrm{los}}+IM_{\mathrm{T}}\right)M_{\mathrm{R}}\bar{Q}_{\mathrm{T}}\bar{Q}_{\mathrm{R}}\right)Q_{\mathrm{T}}Q_{\mathrm{R}}\right)$
in addition to SOMP NF-polar complexity. $S_{\mathrm{los}}$ is the size of the
parameters collection for coarse estimation and $I$ is the number of
iterations of the gradient descent method 111We refer the reader to [24] for
further details. We set $I=20$ and $S_{\mathrm{los}}=45$..

The proposed solution is denoted _Proposed cross-field RD_ in the subsequent
results. Moreover, we evaluate the estimation lower bound by computing the
Oracle least-square [32, 33].

###  IV-B Simulation Parameters and Methodology

The default simulation settings are listed in Table II (modifications are
declared explicitly, THz-specific channel parameters are taken form [10, 21].

TABLE II: Simulation parameters Parameters | Values | Parameters | Values
---|---|---|---
Operating frequency $f_{c}$ | $0.3\,\mathrm{THz}$ |  Refractive index $\kappa_{\mathrm{T}}$ | $2.24-j0.025$
System bandwidth $B_{\mathrm{sys}}$ | $10\,\mathrm{GHz}$ |  Roughness factor $\sigma_{\mathrm{rough}}$ | $0.088\times 10^{-3}\,\mathrm{m}$
Number of subcarriers $K$ | $16$ |  Angle of incidence $\varphi_{\mathrm{inc}}^{\ell}(\ell>1)$ | $\mathcal{U}(0,\pi/2)$
Tx/Rx SAs $Q_{\mathrm{T}}/Q_{\mathrm{R}}$ | $4/4$ |  non-LoS AoD/AoA $\theta^{\ell}/\phi^{\ell}$ | $\mathcal{U}(0,\pi)$
Tx/Rx AEs $\bar{Q}_{\mathrm{T}}/\bar{Q}_{\mathrm{R}}$ | $256/16$ |  non-LoS phase shift $\vartheta^{\ell}$ | $\mathcal{U}(0,2\pi)$
AEs spacing $\delta_{\mathrm{T}}=\delta_{\mathrm{R}}$ | $\frac{\lambda_{c}}{2}$ |  Number of rotations (offline-training) $R_{\mathrm{offline}}$ | $31$
SAs spacing $\Delta_{\mathrm{T}}/\Delta_{\mathrm{R}}$ (Scenario 1) | $128\lambda_{c}/8\lambda_{c}$ |  Number of trials (offline-training) $E_{\mathrm{offline}}$ | $45$
SAs spacing $\Delta_{\mathrm{T}}/\Delta_{\mathrm{R}}$ (Scenario 2) | $256\lambda_{c}/80\lambda_{c}$ |  Number of rotations (online-estimation) $R$ | $9$
Number of channel paths $L$ | $3$ |  Number of trials (online-estimation) $E$ | $45$
Number of paths for SOMP $\hat{L}$ | $10$ |  Noise variance $\sigma^{2}_{n}$ | $-73.8\,\mathrm{dBm}$

We simulate a THz channel with both LoS and non-LoS paths with $L=3$. The Rx
is located at an arbitrary distance $d$. For the LoS path, the distances
between each Tx-Rx SAs are computed based on the exact location of each AE.
Moreover, for non-LoS paths $(\ell>1)$, the distance $d^{\ell}$ is modeled as
a uniform distribution
$d^{\ell}\sim\mathcal{U}(d_{\mathrm{min}},d_{\mathrm{max}})$, where
$d_{\mathrm{min}}=d+\frac{d}{1000}$ and $d_{\mathrm{max}}=6d$. Moreover, since
we model the reflections as a single-bounce (see Fig. 2), the length of the
path from the Tx to the reflectors $d_{\mathrm{T}}^{\ell}$ is modeled as
$d^{\ell}\times\mathcal{U}(0,1)$ and
$d_{\mathrm{R}}^{\ell}=d^{\ell}-d_{\mathrm{T}}^{\ell}$ (see (10)).
Consequently, the time delay for LoS and non-LoS paths are computed based on
$d$ and $d^{\ell}$. We use a compression ratio of $0.5$ at the Tx, i.e.
$M_{\mathrm{T}}=\frac{\bar{Q}_{\mathrm{T}}}{2}$, and full measurements at the
Rx side, i.e. $M_{\mathrm{R}}=\bar{Q}_{\mathrm{R}}$. Two THz scenarios are
simulated where the main difference is that the Tx/Rx UM-MIMO arrays have
different SA spacing, i.e., different apertures as listed in Table II. For
every distance value, in the offline-training, the results are averaged over
$R_{\mathrm{offline}}$ rotations around $\mathrm{Y}$-axis using a random
rotation vector
$\mathbf{o}_{\mathrm{R}}=[180^{\circ},\dot{\beta}_{\mathrm{R}}^{\circ},0^{\circ}]^{\sf
T}$ (where $\dot{\beta}_{\mathrm{R}}\in\mathcal{U}[-30,30]$). For every random
rotation, we use $E_{\mathrm{offline}}=45$ trials. In the online estimation,
we simulate the case of
$\dot{\beta}_{\mathrm{R}}\in\\{-30^{\circ},-22^{\circ},-15^{\circ},-7^{\circ},0^{\circ},7^{\circ},15^{\circ},22^{\circ},30^{\circ}\\}$
and use $E=45$.

We emphasize that the SWM channel model is used in all simulations as it is
the most accurate model for all communication distances. Further, we seek to
minimize the impact of the variation in received signal-to-noise ratio (SNR)
with distance, as we evaluate for different distances and with different
channel models suitable for those distances. As such, we vary the transmit
power $P_{\mathrm{T}}$ so that the average received power remains the same
regardless of the distance. Effectively this is similar to channel
normalization as considered in [22]. The SNR is defined as
$\frac{P_{\mathrm{T}}}{K\sigma^{2}_{n}\rho}$ where $\rho$ is the path loss
defined based on $d$ and (16). The number of data streams is
$N_{\mathrm{S}}=4$. Regarding the Tx/Rx intermediate-/far-field dictionaries,
we set $G_{\mathrm{T}}=\bar{Q}_{\mathrm{T}}$ and
$G_{\mathrm{R}}=\bar{Q}_{\mathrm{R}}$, i.e., no oversampling in the angular
domain. For the polar dictionaries, after optimizing the design parameters
(coherence, minimum allowable distance, etc. [22]) to get a good trade-off in
NMSE and complexity, we get $G_{\mathrm{T}}^{\mathrm{pol}}=2077$ and
$G_{\mathrm{R}}^{\mathrm{pol}}=\bar{Q}_{\mathrm{R}}$.

For the proposed RD, we choose $G_{\mathrm{ovs}}=2$ to construct the Tx/Rx
oversampled intermediate-field (angle-domain)/near-field (polar-domain)
dictionaries, while
$G_{\mathrm{T},\mathrm{RD}}=\frac{\bar{Q}_{\mathrm{T}}}{2}$,
$G_{\mathrm{R},\mathrm{RD}}=\bar{Q}_{\mathrm{R}}$,
$G_{\mathrm{T},\mathrm{RD}}^{\mathrm{pol}}=\frac{G_{\mathrm{T}}^{\mathrm{pol}}}{4}$,
and $G_{\mathrm{R},\mathrm{RD}}^{\mathrm{pol}}=G_{\mathrm{R}}^{\mathrm{pol}}$.

###  IV-C Analysis of Proposed Model Selection Metric

Fig. 7 shows the model selection metric $\eta$ as a function of distance. In
this analysis, we investigate how the number of channel paths and SNR values
impact this metric. The simulation parameters for scenarios 1 and 2 are
detailed in Table II, and each result is obtained by averaging over
$R_{\mathrm{offline}}\times E_{\mathrm{offline}}=1395$ trials. Recall that SWM
is always the most accurate model regardless of the distance and is used to
generate these results.

(a) $\eta$ versus distance for different $L$ and SNR values.

(b) CDFs of $\eta$, and thresholds
$\gamma_{\mathrm{S}\text{-}\mathrm{H}}/\gamma_{\mathrm{H}\text{-}\mathrm{P}}$
for SNR $=6$ dB and $L=3$.

(c) $\eta$ versus distance for different $L$ and SNR values.

(d) CDFs of $\eta$, and thresholds
$\gamma_{\mathrm{S}\text{-}\mathrm{H}}/\gamma_{\mathrm{H}\text{-}\mathrm{P}}$
for SNR $=6$ dB and $L=3$.

Figure 7: Analysis of the proposed model selection metric $\eta$ in the
offline training stage: (a) and (b) for Scenario 1, (c) and (d) for Scenario
2.

The results presented in Fig. 7 clearly indicate that, as anticipated, the
proposed metric $\eta$ exhibits a decreasing trend with distance, even in the
presence of non-LoS paths. As shown in Figs. 7a and 7c, we can identify three
regions based on the $\eta$ values, as intuitively expected from Sec. III-B1.
First, in the near-field, we observe a notable degree of variation in $\eta$.
Moving to the intermediate-field, we witness a sharp decrease in $\eta$. This
drop is a result of some variations in the Rx power as we traverse this
distance range. As we transition to the third region, far-field, the metric is
almost constant.

Moreover, these trends play out differently across SNR levels. In regions with
low SNR, $\eta$ exhibits minimal variations, indicating that it becomes
challenging to favor one channel estimation method over another, as all
methods are expected to result in low estimation accuracy under these
conditions. However, in scenarios with moderate and high SNR levels, our
proposed metric, $\eta$, serves as an effective indicator for selecting the
most suitable channel estimation method. We also specify the regions’
thresholds, to be used in the online operation, in Figs. 7b and 7d for $L=3$
and a moderate SNR value ($6\,\mathrm{dB}$) based on the behavior of $\eta$ in
Figs. 7a and 7c. Note that choosing other threshold values will allow the
control of the fraction of times when one method is preferred over the other.

###  IV-D Evaluation of the Proposed Solution

The performance evaluation of our proposed method relies on several metrics,
which collectively provide a comprehensive assessment. These metrics encompass
decision accuracy (by showing the CDFs of the proposed model selection metric,
$\eta$, during online operations and the pre-selected offline thresholds
$\gamma_{\mathrm{S}\text{-}\mathrm{H}}$ and
$\gamma_{\mathrm{H}\text{-}\mathrm{P}}$), NMSE, AR, and the computational
complexity. To offer a more holistic view, we provide the running time and
EAR. With the used Tx/Rx arrays for both scenarios, each CDF curve is obtained
by averaging over $R\times E=405$ trials.

We also incorporate distances corresponding to the offline thresholds
$d_{\gamma_{\mathrm{S}\text{-}\mathrm{H}}}$ and
$d_{\gamma_{\mathrm{H}\text{-}\mathrm{P}}}$, SA-level MIMO Rayleigh distance
$d_{\mathrm{SA-MIMO-rd}}$, and MIMO Rayleigh distance $d_{\mathrm{MIMO-rd}}$
[24] to emphasize the demarcation between regions. It’s essential to note that
$\eta$ remains relatively constant in the far-field region, and therefore,
$d_{\gamma_{\mathrm{H}\text{-}\mathrm{P}}}$ serves as a coarse indicator
primarily designed to control the probability of selecting PWM-RD.

(a) CDFs of the proposed $\eta$, and thresholds
$\gamma_{\mathrm{S}\text{-}\mathrm{H}}/\gamma_{\mathrm{H}\text{-}\mathrm{P}}$.

(b) NMSE versus distance.

(c) AR versus distance.

(d) Computational Complexity versus distance.

Figure 8: Performance evaluation and comparison for the proposed method with
reference methods as a function of distance for Scenario 1 at SNR $=6$ dB.

The performance evaluation of the proposed methods versus the communication
distance for Scenario 1 is shown in Fig. 8. The CDF of the proposed model
selection metric obtained during online estimations, along with the pre-
computed offline thresholds, are depicted in Fig. 8a. The NMSE, AR analysis,
and numerical evaluation of computational complexity for all estimation
methods are shown in Figs. 8b, 8c, and 8d, respectively.

The first region is the near-field (i.e., encompassing distances less than
$8.5\,\mathrm{m}$ $(d_{\gamma_{\mathrm{S}\text{-}\mathrm{H}}})$), and
showcases the superior performance of the proposed cross-field RD method. Our
estimation method outperforms all other approaches, surpassing them by
approximately $1\,\mathrm{dB}$ or even more at certain distances in terms of
estimation accuracy. For instance, at $d=7\,\mathrm{m}$, cross-field RD
estimation yields a NMSE of $-12\,\mathrm{dB}$, surpassing the best NMSE
achieved by any benchmark method by approximately $2\,\mathrm{dB}$. Taking a
closer look at the additional metrics, we can observe that the proposed method
entails AR values either superior or equivalent to those of all other
reference works. Moreover, the cross-field RD approach demonstrates
significantly lower computational complexity compared to all other prior
works, with the exception of SOMP FF-angular. Note that this marginal
difference in computational complexity is counterbalanced by the substantial
advantages in NMSE and AR achieved by the proposed method. For example, at
short distances, the performance gains are quite significant, with
improvements of approximately $5.5\,\mathrm{dB}$ in NMSE and around
$10\,\mathrm{bits/sec/Hz}$ in AR when comparing cross-field RD to SOMP FF-
angular. Within this region, the SWM-RD dominates the performance, and the
accuracy of selecting it as the proper estimation method is higher than
$98\%$. These notable performance gains clearly underscore the robustness and
demonstrate the effectiveness and superiority of the cross-field RD approach
in near-field scenarios. We introduced perfect channel state information
(PCSI) to AR evaluation, providing an upper bound for this metric. It’s worth
noting that for distances greater than $3\,\mathrm{m}$, the difference between
the proposed method and PCSI is less than $2\,\mathrm{bits/sec/Hz}$. We also
assessed the lower bound for NMSE using the Oracle least-square method [32,
33]. This bound provides a constant value of $-34.1\,\mathrm{dB}$ for NMSE,
irrespective of the communication distance. For clarity, we have omitted the
NMSE lower bound from the results.

The second region is the intermediate-field (spanning approximately
$8.5\,\mathrm{m}$ to $200\,\mathrm{m}$). In this region, the proposed cross-
field RD method consistently outperforms all benchmark works across all
metrics within the distance range of
$d_{\gamma_{\mathrm{S}\text{-}\mathrm{H}}}$ to
$d_{\gamma_{\mathrm{H}\text{-}\mathrm{P}}}$. The selection between
intermediate and far fields can be erroneous - as clear from the CDFs. For
distances greater than $d_{\gamma_{\mathrm{H}\text{-}\mathrm{P}}}$, the
proposed method offers almost the same estimation accuracy as the SOMP FF-
angular estimation, which is the best-performing benchmark in terms of NMSE.
However, it achieves this with a reduced computational complexity. Therefore,
in this region, the proposed cross-field RD method emerges as the preferred
choice. In the far-field region (beyond $200\,\mathrm{m}$), many estimation
strategies exhibit similar or comparable accuracy. Here, we observe a slight
reduction in terms of NMSE (approximately $0.3\,\mathrm{dB}$) and AR (around
$0.25\,\mathrm{bits/sec/Hz}$) for the proposed cross-field RD method when
compared to SOMP FF-angular estimation at distances beyond $d_{\mathrm{MIMO-
rd}}$. Nevertheless, the proposed method maintains lower computational
complexity. The polar-domain dictionaries, used in SOMP NF-polar estimation,
show slight performance degradation compared with angular-domain dictionaries
at large distances since the design of these codebooks is specific for the
near-field. We also note that the performance gap of the proposed method in
terms of NMSE, compared to other methods, decreases as the distance increases.
However, other advantages will be maintained such as lower complexity. In
conclusion, the proposed cross-field RD method outperforms the benchmarks, for
all communication distances, by at least one performance metric. This
emphasizes the validity of our proposal to benefit from the geometric
relations among SAs by using the channel estimate from a SA as prior
information in estimating the channel of subsequent SA. Fig. 9 shows the same
evaluation metrics as in Fig. 8, but for Scenario 2. Different threshold
values have been chosen to determine the regions’ boundaries as illustrated in
Fig. 9a. Results confirm that our proposal can maintain the promised
performance gain for different array configurations in terms of NMSE, AR, and
computational complexity.

(a) CDFs of the proposed $\eta$, and thresholds
$\gamma_{\mathrm{S}\text{-}\mathrm{H}}/\gamma_{\mathrm{H}\text{-}\mathrm{P}}$.

(b) NMSE versus distance.

(c) AR versus distance.

(d) Computational Complexity versus distance.

Figure 9: Performance evaluation and comparison for the proposed methods with
reference methods as a function of distance for Scenario 2 at SNR
$=6\mathrm{dB}$.

We have incorporated the wall-clock running time into our analysis to gain
further insights regarding the computational complexity. The results presented
in Table III were obtained by averaging the running time for each distance,
rotation, and simulation trial. Furthermore, since we consider an AoSA, we
further averaged the time taken by each method per Tx-Rx SA channel estimation
222We have not included the running time of the proposed cross-field RD (PWM-
RD) for two main reasons: (i) it is related to the probability of detecting
PWM in the online operation, and (ii) the estimation is performed only once
for the AoSA.. Results illustrated in Table III indicate that the proposed
cross-field RD algorithm can reduce the running time compared to other
methods.

TABLE III: Comparison of the running time for different methods Method | Running Time (s) | Method | Running Time (s)
---|---|---|---
SOMP NF-polar [22] | 3.725 |  SOMP FF-angular [33] | 0.470
SOMP NF-polar $\hat{L}/2$ [22] | 2.461 |  Mixed-field NF [24] | 3.962
Hybrid-field NF-FF [26] | 2.764 | Proposed cross-field RD (SWM-RD) | 1.144
Hybrid-field FF-NF [26] | 2.764 | Proposed cross-field RD (HSPWM-RD) | 0.225

The EAR results in Fig. 10 for both scenarios show that the proposed method -
in the near-field region - outperforms both SOMP FF-angular and SOMP NF-polar.
The performance gain of cross-field is large compared to SOMP FF-angular. For
most of the intermediate-field distances, the proposed cross-field RD also
does better than the conventional methods. In large communication distances,
all of the methods have almost the same performance, however, the proposed
cross-field RD has a lower computational complexity.

(a) Scenario 1

(b) Scenario 2

Figure 10: The EAR of the proposed methods with reference methods as a
function of distance for different numbers of pilots at SNR $6\,\mathrm{dB}$
and
$T_{\mathrm{coh}}=2Q_{\mathrm{R}}Q_{\mathrm{T}}\bar{Q}_{\mathrm{R}}\bar{Q}_{\mathrm{T}}$
symbols.

##  V Conclusion and Future Work

This paper formulates and solves the cross-field channel estimation problem
for the AoSA-based UM-MIMO THz system. We determine whether near,
intermediate, or far-field channel estimation is appropriate using a model
selection metric and pre-computed thresholds for different regions. In
addition, we use reduced dictionaries to exploit the geometric relation
between multiple SAs. The proposed method outperforms conventional channel
estimation strategies and has (i) lower NMSE, (ii) lower computational
complexity, and (iii) higher AR.

In the future, the cross-field channel estimation problem could be considered
with the beam split effect experienced in wideband systems. Even though the
proposed metric $\eta$ makes intuitive sense - as explained in Sec. III-B1 \-
formal theoretical analysis remains for future work. Currently, we use reduced
dictionaries to exploit the geometric proximity of SAs. In addition, the
design of the training codebooks - e.g., structured codebooks [31] \- as well
as the estimation algorithm - e.g., one with non-uniform prior [35] \- can be
used for enhancements. Further, the cross-field channel estimation problem
needs to be investigated for other wideband hybrid beamforming THz
architectures such as dynamic AoSA [13] and dynamic SA with fixed true-time
delay [14].

## References

  * [1] S. Tarboush, A. Ali, and T. Y. Al-Naffouri, “Compressive estimation of near field channels for ultra massive-MIMO wideband THz systems,” in _Proc. IEEE Int. Conf. Acoustics, Speech, and Signal Process. (ICASSP)_ , 2023, pp. 1–5.
  * [2] Z. Zhang  _et al._ , “6G wireless networks: Vision, requirements, architecture, and key technologies,” _IEEE Veh. Technol. Mag._ , vol. 14, no. 3, pp. 28–41, 2019.
  * [3] T. S. Rappaport  _et al._ , “Wireless communications and applications above 100 GHz: Opportunities and challenges for 6G and beyond,” _IEEE Access_ , vol. 7, pp. 78 729–78 757, 2019.
  * [4] I. F. Akyildiz  _et al._ , “Terahertz band communication: An old problem revisited and research directions for the next decade,” _IEEE Trans. Commun._ , vol. 70, no. 6, pp. 4250–4285, 2022.
  * [5] H. Sarieddeen  _et al._ , “Next generation terahertz communications: A rendezvous of sensing, imaging, and localization,” _IEEE Commun. Mag._ , vol. 58, no. 5, pp. 69–75, 2020.
  * [6] H. Chen et al., “A tutorial on terahertz-band localization for 6G communication systems,” _IEEE Commun. Surveys Tuts._ , vol. 24, no. 3, pp. 1780–1815, May 2022.
  * [7] H. Sarieddeen, M.-S. Alouini, and T. Y. Al-Naffouri, “An overview of signal processing techniques for terahertz communications,” _Proc. IEEE_ , vol. 109, no. 10, pp. 1628–1665, August 2021.
  * [8] C. Han, A. O. Bicen, and I. F. Akyildiz, “Multi-ray channel modeling and wideband characterization for wireless communications in the terahertz band,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 5, pp. 2402–2412, 2014.
  * [9] S. Ju  _et al._ , “Millimeter wave and sub-terahertz spatial statistical channel model for an indoor office building,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 6, pp. 1561–1575, 2021.
  * [10] S. Tarboush  _et al._ , “TeraMIMO: A channel simulator for wideband ultra-massive MIMO terahertz communications,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 12, pp. 12 325–12 341, October 2021.
  * [11] C. Han et al., “Terahertz wireless channels: A holistic survey on measurement, modeling, and analysis,” _IEEE Commun. Surveys Tuts._ , vol. 24, no. 3, pp. 1670–1707, June 2022.
  * [12] I. F. Akyildiz, C. Han, and S. Nie, “Combating the distance problem in the millimeter wave and terahertz frequency bands,” _IEEE Commun. Mag._ , vol. 56, no. 6, pp. 102–108, 2018.
  * [13] C. Han, L. Yan, and J. Yuan, “Hybrid beamforming for terahertz wireless communications: Challenges, architectures, and open problems,” _IEEE Wireless Commun._ , vol. 28, no. 4, pp. 198–204, 2021.
  * [14] L. Yan, C. Han, and J. Yuan, “Energy-efficient dynamic-subarray with fixed true-time-delay design for terahertz wideband hybrid beamforming,” _IEEE J. Sel. Areas Commun._ , vol. 40, no. 10, pp. 2840–2854, 2022.
  * [15] C. Lin and G. Y. Li, “Terahertz communications: An array-of-subarrays solution,” _IEEE Commun. Mag._ , vol. 54, no. 12, pp. 124–131, 2016.
  * [16] ——, “Energy-efficient design of indoor mmWave and sub-THz systems with antenna arrays,” _IEEE Trans. Wireless Commun._ , vol. 15, no. 7, pp. 4660–4672, 2016.
  * [17] C. Han  _et al._ , “Cross far-and near-field wireless communications in terahertz ultra-large antenna array systems,” _arXiv preprint arXiv:2301.03035_ , 2023.
  * [18] M. Cui  _et al._ , “Near-field MIMO communications for 6G: Fundamentals, challenges, potentials, and future directions,” _IEEE Commun. Mag._ , vol. 61, no. 1, pp. 40–46, 2023.
  * [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._ , vol. 69, no. 10, pp. 7063–7076, October 2021.
  * [20] S. Yang  _et al._ , “Performance bounds for near-field localization with widely-spaced multi-subarray mmWave/THz MIMO,” _arXiv preprint arXiv:2309.05944_ , 2023.
  * [21] K. Dovelos  _et al._ , “Channel estimation and hybrid combining for wideband terahertz massive MIMO systems,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 6, pp. 1604–1620, 2021.
  * [22] 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, 2022.
  * [23] X. Zhang, H. Zhang, and Y. C. Eldar, “Near-field sparse channel representation and estimation in 6G wireless communications,” _arXiv preprint arXiv:2212.13527_ , 2022.
  * [24] Y. Lu and L. Dai, “Near-field channel estimation in mixed LoS/NLoS environments for extremely large-scale MIMO systems,” _IEEE Trans. Commun._ , 2023, Early access.
  * [25] Y. Chen  _et al._ , “Hybrid spherical-and planar-wave channel modeling and estimation for terahertz integrated UM-MIMO and IRS systems,” _IEEE Trans. Wireless Commun._ , 2023, Early access.
  * [26] X. Wei and L. Dai, “Channel estimation for extremely large-scale massive MIMO: Far-field, near-field, or hybrid-field?” _IEEE Commun. Lett._ , vol. 26, no. 1, pp. 177–181, 2021.
  * [27] Z. Hu  _et al._ , “Hybrid-field channel estimation for extremely large-scale massive MIMO system,” _IEEE Commun. Lett._ , vol. 27, no. 1, pp. 303–307, 2022.
  * [28] W. Yu  _et al._ , “An adaptive and robust deep learning framework for THz ultra-massive MIMO channel estimation,” _IEEE J. Sel. Topics Signal Process._ , 2023.
  * [29] S. Tarboush  _et al._ , “Single- versus multicarrier terahertz-band communications: A comparative study,” _IEEE Open J. of the Commun. Soc._ , vol. 3, pp. 1466–1486, August 2022.
  * [30] A. M. Sayeed, “Deconstructing multiantenna fading channels,” _IEEE Trans. Signal Process._ , vol. 50, no. 10, pp. 2563–2579, 2002.
  * [31] A. Ali, N. González-Prelcic, and R. W. Heath, “Millimeter wave beam-selection using out-of-band spatial information,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 2, pp. 1038–1052, 2017.
  * [32] 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._ , vol. 64, no. 6, pp. 2370–2386, 2016.
  * [33] J. Rodríguez-Fernández  _et al._ , “Frequency-domain compressive channel estimation for frequency-selective hybrid millimeter wave MIMO systems,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 5, pp. 2946–2960, 2018.
  * [34] O. El Ayach  _et al._ , “Spatially sparse precoding in millimeter wave MIMO systems,” _IEEE Trans. Wireless Commun._ , vol. 13, no. 3, pp. 1499–1513, 2014.
  * [35] J. Scarlett, J. S. Evans, and S. Dey, “Compressed sensing with prior information: Information-theoretic limits and practical decoders,” _IEEE Trans. Signal Process._ , vol. 61, no. 2, pp. 427–439, 2012.