# Angularly Sparse Channel Estimation in Dual-Wideband Tera-Hertz (THz) Hybrid
MIMO Systems Relying on Bayesian Learning

Abhisha Garg Suraj Srivastava, Nimish Yadav, Aditya K. Jagannatham, and Lajos
Hanzo L. Hanzo would like to acknowledge the financial support of the
Engineering and Physical Sciences Research Council projects EP/W016605/1,
EP/X01228X/1 and EP/Y026721/1 as well as of the European Research Council’s
Advanced Fellow Grant QuantCom (Grant No. 789028).
The work of Aditya K. Jagannatham was supported in part by the Qualcomm
Innovation Fellowship; in part by the Qualcomm 6G UR Gift; and in part by the
Arun Kumar Chair Professorship.
The work is partially supported by the IEEE SPS Scholarship grant for $2023$.
Abhisha Garg, Nimish Yadav and Aditya K. Jagannatham are with the Department
of Electrical Engineering, Indian Institute of Technology Kanpur,
Kanpur-$208016$, India (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>adityaj@iitk.ac.in). Suraj Srivastava is with the Department of Electrical
Engineering, Indian Institute of Technology Jodhpur, Jodhpur, Rajasthan
$342030$, India (email: surajsri@iitj.ac.in). L. Hanzo is with the School of
Electronics and Computer Science, University of Southampton, Southampton SO17
1BJ, U.K. (email<EMAIL_ADDRESS>
###### Abstract

Bayesian learning aided massive antenna array based THz MIMO systems are
designed for spatial-wideband and frequency-wideband scenarios, collectively
termed as the dual-wideband channels. Essentially, numerous antenna modules of
the THz system result in a significant delay in the transmission/ reception of
signals in the time-domain across the antennas, which leads to spatial-
selectivity. As a further phenomenon, the wide bandwidth of THz communication
results in substantial variation of the effective angle of arrival/ departure
(AoA/ AoD) with respect to the subcarrier frequency. This is termed as the
beam squint effect, which renders the channel state information (CSI)
estimation challenging in such systems. To address this problem, initially, a
pilot-aided (PA) Bayesian learning (PA-BL) framework is derived for the
estimation of the Terahertz (THz) MIMO channel that relies exclusively on the
pilot beams transmitted. Since the framework designed can successfully operate
in an ill-posed model, it can verifiably lead to reduced pilot transmissions
in comparison to conventional methodologies. The above paradigm is
subsequently extended to additionally incorporate data symbols to derive a
Data-Aided (DA) BL approach that performs joint data detection and CSI
estimation. We will demonstrate that it is capable of improving the dual-
wideband channel’s estimate, despite further reducing the training overhead.
The Bayesian Cramér-Rao bounds (BCRLBs) are also obtained for explicitly
characterizing the lower bounds on the mean squared error (MSE) of the PA-BL
and DA-BL frameworks. Our simulation results show the improved normalized MSE
(NMSE) and bit-error rate (BER) performance of the proposed estimation schemes
and confirm that they approach their respective BCRLB benchmarks.

###### Index Terms:

Terahertz, MIMO, dual-wideband, pilot-aided, data-aided, beam squint, Bayesian
learning

##  I Introduction

Terahertz (THz) communication, relying on the (0.1-10) THz band, has the
promise of achieving the ultra-high data rates targeted by next generation
networks [1]. When compared to the mmWave band that ranges from (30-100) GHz,
the THz band can provide spectral windows with bandwidths in the range of 10
GHz or even higher [2]. This can in turn lead to the development of rapid data
transfer capabilities links with rates in the range of several hundreds of
Gbps[3]. However, it must be noted that the successful achievement of THz
communication is enormously challenging, as the THz signal suffers from high
atmospheric losses ranging from 5dB to 100dB [4] and path losses ranging from
70dB to 140 dB [2], which can be attributed to the absorption by gases and the
high carrier frequencies, respectively. In addition, THz signals are prone to
blockage by physical barriers such as walls, human bodies etc [5]. Multiple-
input multiple-output (MIMO) technology, which is capable of generating narrow
directional beams having high array gains, is extremely well-suited for
mitigating the losses occurring in the THz band. Fortuitously, numerous
antenna elements can be embedded at both the transceiver ends in a THz system
within a small form-factor [6]. At this juncture, it should be noted that the
conventional transceiver architecture that necessitates a separate radio
frequency (RF) chain per antenna can lead to a an excessive increase in the
hardware complexity required to carry out the signal processing operations
pertaining to modulation and demodulation [7]. A popular solution recommended
in the latest research literature [8, 9] to overcome this impediment is the
novel hybrid signal processing architecture, wherein the signal processing
operations are effectively split between the RF and baseband stages. This is
in stark contrast to an orthodox transceiver design wherein the entire signal
processing burden is borne by the baseband processor. The hybrid architecture
coupled with large MIMO arrays can enormously increase the signal gain via
optimal beamforming, while also limiting the hardware complexity through joint
RF and baseband signal processing.

Moreover, a typical wideband THz channel is frequency-selective in nature, due
to the multipath induced delay spread, which is a characteristic feature of
wireless signal propagation [10]. This phenomenon is termed as frequency-
wideband effect. A pair of crucial factors contributing to frequency
selectivity in the THz band, even for LoS scenarios, are high atmospheric
absorption and free-space path losses [11]. The heightened atmospheric
absorption primarily arises from water vapor molecules. Furthermore, due to
frequency-selective nature of the absorption coefficient, the THz band is
characterized by ultra-wideband transmission windows (TWs) and absorption
coefficient peak regions (ACPRs) [2]. Additionally, the presence of ACPRs in
the THz band leads to pulse broadening in the time domain, presenting
practical communication challenges within this frequency range. Owing to the
short wavelength of the THz signal, a large antenna array leads to a
progressively increasing phase shift in the received/ transmitted signals
across the constituent antennas, corresponding to commensurate delays in the
time domain [12]. This is termed as the spatial-wideband effect. This leads to
the frequency dependent nature of the array response [13]. A general THz
system that is exposed to both frequency-wideband and spatial-wideband effects
is termed as a dual-wideband system, which is considered in our work.
Furthermore, a practical wideband THz system suffers from the beam squint
problem that can be attributed to the frequency dependent nature of the array
response [13]. The availability of accurate channel state information (CSI),
which plays a key role in beamforming, is of paramount importance in such
systems. However, conventional techniques such as the least-squares (LS) and
linear minimum mean squared error (LMMSE) estimators, require an excessive
pilot overhead, owing to the large number of antennas, which reduces the
spectral efficiency (SE). Thus, efficient CSI estimation techniques, which
exploit the properties of the dual-wideband THz MIMO channel, hold the key for
the successful realization of reliable high speed communication in these
systems which also considers the beam squint effect. A concise overview of the
relevant prior research is provided next.

###  I-A Review of the existing literature

Jornet and Akyildiz [14] developed a novel THz channel model accounting also
for the molecular absorption and path loss, which was successfully leveraged
to determine the capacity of THz systems. In [15], the principles of molecular
absorption were combined with radioactive transfer theory to evaluate the path
loss of THz signals for propagation distances in the range of a centimeter to
a few meters, followed by the formulation of the THz channel model. Lin et
al., [16], developed a general array-of-subarrays architecture for a hybrid
THz MIMO system, complemented by an innovative algorithm for adaptive
distance-aware beamforming, which requires knowledge of the CSI. Therefore, it
is important to note at this juncture that the end-to-end performance of a THz
band hybrid MIMO system is crucially contingent on the availability of
accurate CSI. This motivates the design of superior channel estimation
techniques for THz MIMO applications, which can potentially maximize the
performance. The well-known traditional channel estimation approaches, such as
the LS and MMSE, perform poorly in the THz regime owing to their exclusive
pilot overhead, which consequently leads to a low spectral efficiency [17].
Nevertheless, a significant point to note is that the THz frequency range
exhibits a strong directional propagation characteristic owing to using a
large antenna array, hence resulting in a sparsely populated multipath channel
in the angular domain. Naturally, sparse signal recovery techniques, which
have gained popularity in recent research, can lead to verifiably improved CSI
estimates in THz MIMO transceivers. In this context, the authors of [18]
design a superimposed pilot modulation technique, followed by sparse Bayesian
learning (SBL)-based approaches for the estimation of the frequency-selective
wideband channel. More specifically, they proposed a multiple-measurement
vector (MMV)-based SBL algorithm for joint CSI estimation and data-detection.
Venugopal et al., in their cutting-edge work [19], proposed techniques for the
estimation of a frequency-selective mmWave wideband channel in both single-
carrier frequency domain equalization (SC-FDE) and orthogonal frequency-
division multiplexing (OFDM)-based mmWave systems. In their work, the CSI
estimation model was formulated as a sparse recovery problem, which was
subsequently solved via a time-frequency estimation algorithm. The authors of
[20] developed a block least-mean squares (BLMS) sparse adaptive framework for
estimating the channel in a mmWave MIMO SC-FDE system. As a further advance,
Gao et al., in their pioneering treatise [21], proposed a sparse beamspace CSI
estimation technique for the estimation of a wideband mmWave MIMO system
having lens antenna arrays. The novel scheme proposed in their work is based
on successive support detection (SSD), which is inspired by the successive
interference cancellation (SIC) principle of multi-user systems. Wang et al.
[22] conceived new CSI estimation techniques for both frequency division
duplex (FDD) as well as time division duplex (TDD) mmWave massive MIMO
systems. Their framework also considers the spatial-wideband and frequency-
wideband effects inherent in such systems and exploits the sparsity both in
the angular and delay domains, while comprehensively covering both the uplink
as well as the downlink. An attractive sensing-based high-resolution technique
was designed in [23] that can successfully extract the physical parameters of
the uplink CSI, viz., AoAs, time delays and complex channel gains, while
accounting for the beam squint effect in mmWave hybrid MIMO systems. This is
achieved by adaptively updating the dictionary matrix. A few of the latest
treatises have also considered the impact of the spatial-wideband and
frequency-wideband effects in massive MIMO THz channel estimation, which is of
considerable practical importance. Dovelos et al. [24] model an OFDM-based
wideband THz system using a uniform planar array (UPA) under the spatial-
frequency wideband effect. A sparse beamspace representation of the wideband
THz channel is formulated in their work, followed by the application of the
orthogonal matching pursuit (OMP) and the generalized simultaneous OMP (GSOMP)
algorithm popular for the acquisition of the sparse channel. Chou et al., [25]
addressed the spatial-wideband and frequency-wideband effect with respect to
the estimation of a time-varying MIMO-OFDM channel in sub-Terahertz massive
MIMO systems. The authors in their cutting-edge work developed a multiple
measurement vector least squares compressed sensing (MMV-LS-CS) based CSI
estimation technique, which utilizes the MMV-LS to recover the channel
employing the previous channel support and the MMV-CS to detect any time-
varying components in the residual signal.

Table I: Boldly contrasting our contribution to the literature Features | [6] | [8] | [9] | [13] | [17] | [22] | [24] | [25] | [26] | This paper
---|---|---|---|---|---|---|---|---|---|---
THz hybrid MIMO | ✓ |  |  |  |  |  | ✓ | ✓ |  | ✓
Dual-Wideband effect |  |  | ✓ |  |  | ✓ |  | ✓ |  | ✓
Reflection losses | ✓ |  |  |  |  |  | ✓ | ✓ |  | ✓
Molecular absorption losses | ✓ |  |  |  |  |  | ✓ | ✓ |  | ✓
Diffused-ray modeling | ✓ | ✓ |  |  |  |  |  |  |  | ✓
BCRLB | ✓ | ✓ |  | ✓ |  |  | ✓ |  | ✓ | ✓
SC-FDE system |  |  |  |  |  |  |  |  |  | ✓
Data-Aided CSI acquisition |  |  |  |  |  |  |  |  |  | ✓

Although their pioneering work incorporates time-selectivity in sub-THz
systems, no joint channel estimation and data detection is considered, which
can lead to a remarkable enhancement in performance. To fill this knowledge
gap, we conceive a novel DA Bayesian learning framework for reliable sparse
CSI estimation and data-detection in a THz MIMO system, while considering also
the dual-wideband effect. In Table I we boldly contrast our contributions to
the literature. The next section briefly describes the main contributions of
our work in more details.

###  I-B Contributions of this work

  1. 1.

A hybrid frequency-dependent THz MIMO channel model is developed by
considering the dual-wideband effects, which additionally incorporates the
absorption, reflection as well as free space losses for a uniform linear array
(ULA) at the base station (BS) serving a single user. Notably, our THz channel
model incorporates the spatial wideband effect to construct a frequency-
dependent model of the array response, while also considering the frequency-
wideband effect, which accounts for multipath propagation. Existing literature
on THz channel estimation, such as [24] and [25], have solely focused on the
spatial-wideband effect and did not include the frequency-wideband effect by
considering multipath delays and the sampled version of transmit pulse shaping
filter response. Therefore, the novel contribution of our work lies in
considering the dual wideband effects, which are missing in the existing
literature such as [27], [28]. Moreover, an efficient zero-padded (ZP) block
structure is proposed to transform the time domain (TD) signal processing
problem into an equivalent frequency domain (FD) one, which remarkably reduces
the complexity of THz MIMO CSI estimation as well as signal detection.
Additionally, this novel transmission structure requires very few pilot beams
to excite the numerous angular modes of the channel.

  2. 2.

A sparse signal recovery problem is formulated for the THz MIMO channel. A
pair novel pilot-aided Bayesian learning (PA-BL) and data-aided Bayesian
learning (DA-BL)-based methodologies are proposed for estimating the sparse
beamspace dual-wideband channel. Since it can operate in ill-posed scenerios,
it requires significantly fewer pilot transmissions in comparison to the
traditional LS and MMSE schemes.

  3. 3.

The Bayesian Cramér-Rao lower bounds (BCRLB) are derived for the mean-square
error (MSE) of both the PA and DA estimators proposed in the THz regime. These
serve as excellent benchmarks to gauge the performance of the proposed
algorithm.

  4. 4.

Simulation results are presented for quantifying the CSI estimation and data-
detection performance achieved in the dual-wideband THz MIMO system
considered. It is shown that the normalized MSEs (NMSEs) of the PA-BL and DA-
BL algorithms approach their respective BCRLBs, which demonstrates their
efficiency. This is especially significant considering the fact that they do
not require prior knowledge of either the statistics or of the support of the
sparse CSI.

###  I-C Organization of the paper

The paper is structured as follows. Section II derives the system model of a
hybrid massive MIMO system in the THz regime utilizing the SC-FDE framework.
Section III motivates and models the dual-wideband effect in the above system
and also characterizes the loss induced by the beam squint effect. Section IV
determines the model for sparse CSI estimation in the wideband THz system.
Section V describes the PA-BL approach proposed for CSI estimation, with the
requisite mathematical foundation. Next, Section VI chronicles the DA-BL
approach, which yields an improved performance over the previous paradigm.
Section VII outlines our simulation results, while Section VIII concludes the
paper.

###  I-D Notation

The following notation is used throughout this paper. Letters in upper-case
boldface such as $\mathbf{A}$ represent matrices, while letters in lower-case
boldface denote vectors. Superscripts $(.)^{\mathit{T}}$, $(.)^{*}$,
$(.)^{\textrm{-1}}$, $(.)^{\mathit{H}}$, and $(.)^{\dagger}$ represent the
transpose, conjugate, inverse, Hermitian and pseudo-inverse operators,
respectively. The quantity
$\textrm{blkdiag}\left(\mathbf{A}_{\textrm{1}},\mathbf{A}_{\textrm{2}},\cdots,\mathbf{A}_{\mathit{N}}\right)$
represents the block-diagonal matrix with matrices
$\mathbf{A}_{\textrm{1}},\mathbf{A}_{\textrm{2}},\cdots,\mathbf{A}_{\mathit{N}}$
on the principal diagonal, which are not necessarily square. The vectorization
operation that stacks the columns of any matrix $\mathbf{A}$ is represented by
$\textrm{vec}(\mathbf{A})$, while the inverse of the same operation is
represented by $\textrm{vec}^{-1}(\mathbf{a})$. The operators $*$ and
$\otimes$ represent the linear convolution and Kronecker product,
respectively, while $\textrm{Tr}(.)$ represents the matrix trace operator. The
following standard property of the $\textrm{vec}(.)$ operation,
$\textrm{vec}(\mathbf{A}\mathbf{B}\mathbf{C})=(\mathbf{C}^{T}\otimes\mathbf{A})\textrm{vec}(\mathbf{B})$,
is also used in the paper. The operators $\parallel.\parallel_{\textrm{2}}$
and $\parallel.\parallel_{\mathit{F}}$ represent the $l_{\textrm{2}}$-norm and
Frobenius norm, respectively, of vectors and matrices. The expectation
operator is denoted by $\mathbb{E}\left\\{.\right\\}$, while the absolute
value is denoted by $|\ .\ |$; the symmetric complex Gaussian distribution
with mean $\boldsymbol{\mu}$ and variance $\boldsymbol{\Sigma}$ is represented
by $\mathcal{CN}\left(\boldsymbol{\mu},\boldsymbol{\Sigma}\right)$.

##  II SC-FDE Based THz Hybrid MIMO System Model

Figure 1: Schematic representation of SC-FDE-based wideband THz hybrid MIMO
system

Consider a single carrier wideband THz hybrid MIMO system that transmits
$N_{s}$ data streams from a transmitter equipped with $N_{T}$ transmit
antennas (TAs), to a receiver having $N_{R}$ receive antennas (RAs). There are
$N_{RF}$ RF chains at the transmitter and the receiver, where $N_{s}\leq
N_{RF}\ll\mathrm{min}(N_{T},N_{R})$. The THz hybrid MIMO transmitter is
equipped with a baseband precoder in cascade with an RF precoder, represented
by the matrices $\mathbf{F_{\mathit{BB}}}$ $\in$ $\mathbb{C}^{N_{RF}\times
N_{s}}$ and $\mathbf{F_{\mathit{RF}}}$ $\in$ $\mathbb{C}^{N_{T}\times
N_{RF}}$, respectively. In a similar manner, at the RAs, the hybrid processing
module is comprised of an RF beamformer $\mathbf{W_{\mathit{RF}}}$ $\in$
$\mathbb{C}^{N_{R}\times N_{RF}}$ in cascade with a baseband combiner
$\mathbf{W_{\mathit{BB}}}$ $\in$ $\mathbb{C}^{N_{RF}\times N_{s}}$. The hybrid
signal processing architecture of the THz MIMO system is shown in Fig. 1.
Since the RF transmit precoder (TPC) matrix $\mathbf{F_{\mathit{RF}}}$ and RF
receiver combiner (RC) matrix $\mathbf{W_{\mathit{RF}}}$ are constrained to
exclusively comprise digitally controlled phase-shifters, having elements
obeying $|\mathbf{F_{\mathit{RF}}}(i,j)|=\frac{1}{\sqrt{N_{T}}}$,
$|\mathbf{W_{\mathit{RF}}}(i,j)|=\frac{1}{\sqrt{N_{R}}}$ [19]. Furthermore,
$\mathbf{F}_{RF}$ and $\mathbf{W}_{RF}$ can be represented as
$\mathbf{F}_{\mathit{RF}}=\left[\mathbf{f}_{1},\mathbf{f}_{2},\cdots,\mathbf{f}_{N_{RF}}\right]$
and
$\mathbf{W}_{\mathit{RF}}=\left[\mathbf{w}_{1},\mathbf{w}_{2},\cdots,\mathbf{w}_{N_{RF}}\right]$,
respectively. Considering a frequency-selective wideband channel in the THz
band comprised of $\mathit{L}$ delay taps, the $\mathit{l}$th delay tap can be
expressed by the complex channel matrix $\mathbf{H}_{\mathit{l}}$,
$0\leqslant\mathit{l}\leqslant\mathit{L}-1$, of size $N_{R}\times N_{T}$. The
symbol vector that is transmitted at time instant $\mathit{t}$ is denoted by
${\mathbf{s}}(\mathit{t})$ $\in$ $\mathbb{C}^{N_{s}\times 1}$. After baseband
combining, the output signal vector at the $\mathit{t}$th time instant is
given by

$\displaystyle\check{\mathbf{y}}({\mathit{t}})=\mathbf{W}\mathit{{}_{BB}^{H}}\mathbf{W}\mathit{{}_{RF}^{H}}[\mathbf{H}_{\mathit{t}}\ast\mathbf{x}(\mathit{t})]+\mathbf{W_{\mathit{BB}}^{\mathit{H}}}\mathbf{W_{\mathit{RF}}^{\mathit{H}}}\mathbf{v({\mathit{t}})},$
(1)

where
$\mathbf{{x}}(\mathit{t})=\mathbf{F_{\mathit{RF}}}\mathbf{F_{\mathit{BB}}}{\mathbf{s}}(\mathit{t})\in$
$\mathbb{C}^{N_{\mathit{T}}\times 1}$ represents the precoded symbol vector
and $\mathbf{v({\mathit{t}})}\in\mathbb{C}^{N_{\mathit{R}}\times 1}$
represents the additive noise vector that follows the distribution
$\mathcal{CN}(\mathbf{0}_{N_{R}\times
1},\sigma_{\mathit{v}}^{2}\mathbf{I}_{N_{R}})$. The operator $\ast$ represents
linear convolution. The equivalent FD model of this system is developed next.

###  II-A SC-FDE channel estimation model

Again Fig. 1 depicts the architecture of a hybrid THz MIMO system relying on
FD signal processing. Each transmit frame duration in the system is divided
into training and data phases, as shown in Fig. 2. The training phase is
comprised of $\mathit{M}$ blocks with each individual block consisting of
$N_{p}$ pilot vectors. Let
$\mathbf{u}_{\mathit{m}}^{(\mathit{p})}\in\mathbb{C}^{\mathit{N_{RF}}\times
1}$ denote the $\mathit{p}$th complex pilot vector, which follows the
relationship $0\leqslant\mathit{p}\leqslant\mathit{N_{p}-\mathrm{1}}$ for
every $\mathit{m}$th block [20]. Prior to transmitting the pilot sequence on
an individual RF chain, $\mathit{L}-1$ zeros are appended in each
$\mathit{m}$th block to generate a zero-padded (ZP) block of length
$\mathit{K=N_{p}+L}-1$. Similarly, zero-matrices of size
$\mathit{N}_{R}\times\mathit{N}_{T}$ are appended to the MIMO channel taps
$\mathbf{H}_{l}$ to obtain a block of length $\mathit{K}$. Thus, the zero-
padded pilot vectors and THz MIMO channel taps can be described in a compact
fashion as
$\left\\{{\mathbf{u}_{\mathit{m}}^{(\mathit{q})}}\right\\}_{\mathit{q}=0}^{\mathit{K}-1}=\left\\{{{\mathbf{u}_{\mathit{m}}^{(0)}},{\mathbf{u}_{\mathit{m}}^{(1)}},\cdots,{\mathbf{u}_{\mathit{m}}^{(\mathit{\mathit{N}_{p}}-1)}},\underbrace{\mathbf{0},\cdots,\mathbf{0}}_{\mathit{L}-1}}\right\\}$
and
$\left\\{\mathbf{H}_{\mathit{q}}\right\\}_{\mathit{q}=0}^{\mathit{K}-1}=\left\\{\mathbf{H}_{0},\mathbf{H}_{1},\cdots,\mathbf{H}_{\mathit{L}-1},\underbrace{\mathbf{0},\cdots,\mathbf{0}}_{\mathit{N}_{p}-1}\right\\}$,
respectively. For each block $\mathit{m}$, the training RF TPC and RC are
represented as
$\mathbf{F}_{\mathit{RF,m}}\in\mathbb{C}^{\mathit{N_{T}}\times\mathit{N_{RF}}}$
and
$\mathbf{W}_{\mathit{RF,m}}\in\mathbb{C}^{\mathit{N_{R}}\times\mathit{N_{RF}}}$,
respectively. Thus, the output signal vectors
$\widetilde{\mathbf{y}}_{\mathit{m}}(\mathit{q})\in\mathbb{C}^{\mathit{N_{RF}}\times
1}$ corresponding to the $\mathit{m}$th block follow the model given by

$\displaystyle\widetilde{\mathbf{y}}_{\mathit{m}}(\mathit{q})=\mathbf{W}_{\mathit{RF,m}}^{H}(\mathbf{H}_{q}\circledast_{\mathit{K}}\widetilde{\mathbf{u}}_{\mathit{m}}^{\mathit{(q)}})+\mathbf{W}_{\mathit{RF,m}}^{H}\mathbf{v}_{\mathit{m}}(\mathit{q}).$
(2)

There
$\widetilde{\mathbf{u}}_{\mathit{m}}^{\mathit{(q)}}=\mathbf{F}_{\mathit{RF,m}}{\mathbf{u}_{\mathit{m}}^{(\mathit{q})}}$,
$\mathbf{v}_{\mathit{m}}(\mathit{q})$ denotes the complex white AWGN noises
which follows the distribution $\mathcal{CN}(\mathbf{0}_{\mathit{N}_{r}\times
1},\sigma^{2}_{v}\mathbf{I}_{\mathit{N}_{r}})$, and $\circledast_{\mathit{K}}$
is the circular convolution operation between sequences of length $K$.
Consider the $\mathit{m}$th TD output block defined as
$\left\\{\widetilde{\mathbf{y}}_{\mathit{m}}(\mathit{q})\right\\}_{q=0}^{\mathit{K}-1}=\left\\{\widetilde{\mathbf{y}}_{\mathit{m}}(0),\widetilde{\mathbf{y}}_{\mathit{m}}(1),\cdots,\widetilde{\mathbf{y}}_{\mathit{m}}(\mathit{K}-1)\right\\}$.
At the output, employing the $\mathit{K}$-point fast fourier transform (FFT)
operation, the FD representation that corresponds to the hybrid MIMO system in
the THz band may be obtained as

$\displaystyle\left\\{\mathbf{y}_{\mathit{m}}[\mathit{k}]\right\\}_{\mathit{k}=0}^{\mathit{K}-1}$
$\displaystyle=\textrm{FFT}(\left\\{\widetilde{\mathbf{y}}_{\mathit{m}}(\mathit{q})\right\\}_{q=0}^{\mathit{K}-1})$
$\displaystyle=\left\\{{{\mathbf{y}}_{\mathit{m}}[0],{\mathbf{y}}_{\mathit{m}}[1],\cdots,{\mathbf{y}}_{\mathit{m}}[\mathit{K}-1]}\right\\}.$
(3)

Thus, the received FD output vector
${\mathbf{y}}_{\mathit{m}}[\mathit{k}]\in\mathbb{C}^{\mathit{N}_{\mathit{RF}}\times
1}$ at the $\mathit{k}$th subcarrier can be expressed as:

$\displaystyle{\mathbf{y}}_{\mathit{m}}[\mathit{k}]=\mathbf{W}_{\mathit{RF,m}}^{\mathit{H}}\mathbf{H}[\mathit{k}]\mathbf{F}_{\mathit{RF,m}}\mathbf{u}_{\mathit{m}}[\mathit{k}]+\mathbf{W}_{\mathit{RF,m}}^{\mathit{H}}\check{\mathbf{v}}_{\mathit{m}}[\mathit{k}],$
(4)

where $\mathbf{H}[k]$ depicts the THz MIMO channel as frequency response
(CFR), whose model is detailed in the next section. The vector
$\mathbf{u}_{\mathit{m}}[\mathit{k}]\in\mathbb{C}^{\mathit{N}_{RF}\times 1}$
denotes the $\mathit{k}$th FFT point of the pilot sequence
$\mathbf{u}_{\mathit{m}}^{(\mathit{q})}$, and
$\check{\mathbf{v}}_{\mathit{m}}[\mathit{k}]\in\mathbb{C}^{\mathit{N}_{R}\times
1}$ represents the $\mathit{k}$th FFT point of the noise sequence
$\mathbf{v}_{\mathit{m}}^{(\mathit{q})}$. Upon exploiting the property of
vec(.) operator as described in Section-I-C, (4) can be recast as

$\displaystyle\mathbf{{y}}_{\mathit{m}}[\mathit{k}]=\underbrace{(\mathbf{u}_{\mathit{m}}^{\mathit{T}}[\mathit{k}]\mathbf{F}_{\mathit{RF,m}}^{\mathit{T}}\otimes\mathbf{W}_{\mathit{RF,m}}^{\mathit{H}})}_{\mathbf{\Phi}_{\mathit{m}}[\mathit{k}]\in\mathbb{C}^{\mathit{N}_{\mathit{RF}}\times\mathit{N}_{\mathit{R}}\mathit{N}_{\mathit{T}}}}\mathbf{h}[\mathit{k}]+\mathbf{{v}}_{\mathit{m}}[\mathit{k}],$
(5)

where
$\mathbf{h}[\mathit{k}]=\textrm{vec}(\mathbf{H}[\mathit{k}])\in\mathbb{C}^{\mathit{N}_{T}\mathit{N}_{R}\times
1}$ represents the vectorized form of the THz MIMO channel model at the $k$th
subcarrier, $\mathbf{\Phi}_{\mathit{m}}[\mathit{k}]$ represents the sensing
matrix for the $\mathit{m}$th block of the $\mathit{k}$th subcarrier, while
$\mathbf{{v}}_{\mathit{m}}[\mathit{k}]=\mathbf{W}_{\mathit{RF,m}}^{\mathit{H}}\check{\mathbf{v}}_{\mathit{m}}[\mathit{k}]\in\mathbb{C}^{\mathit{N}_{RF}\times
1}$ is the noise vector at the output of the RF combiner of the $\mathit{k}$th
subcarrier. Note also that,

$\displaystyle\mathbb{E}\left\\{\mathbf{{v}}_{\mathit{m}}[\mathit{k}]\mathbf{{v}}_{\mathit{m}}^{H}[\mathit{k}]\right\\}$
$\displaystyle=\mathbb{E}\left\\{\mathbf{W}_{\mathit{RF,m}}^{\mathit{H}}\check{\mathbf{v}}_{\mathit{m}}[\mathit{k}]\left(\mathbf{W}_{\mathit{RF,m}}^{\mathit{H}}\check{\mathbf{v}}_{\mathit{m}}[\mathit{k}]\right)^{H}\right\\}$
$\displaystyle=\sigma^{2}_{v}K\mathbf{W}_{\mathit{RF,m}}^{\mathit{H}}\mathbf{W}_{\mathit{RF,m}}=\mathbf{R}_{m}.$
(6)

The next section describes the THz channel incorporating the dual-wideband
effect.

Figure 2: Frame structure employed in the wideband THz hybrid MIMO system
using SC-FDE.

##  III THz Channel Model with dual-wideband effect

We now develop an explicit model for the THz MIMO channel incorporating dual
wideband effect, comprising of the spatial and frequency wideband effects. We
commence by describing the spatial wideband effect.

###  III-A Spatial-wideband effect

Let us assume that a signal arrives at a ULA having $N$ identical antennas,
indexed as $\left\\{1,2,\cdots,\mathit{N}\right\\}$, at the azimuth angle
$\theta$. The frequency-dependent array steering vector $\mathbf{a}(\theta,f)$
of this system is given by

$\displaystyle\mathbf{a}(\theta,\mathit{f})=\frac{1}{\sqrt{\mathit{N}}}\left[1,e^{-j\frac{2\pi}{\lambda}\mathit{d}_{c}\cos\theta},\cdots,e^{-j(\mathit{N}-1)\frac{2\pi}{\lambda}\mathit{d}_{c}\cos\theta}\right]^{T}.$
(7)

The quantity $\lambda=\frac{c}{f}$ represents the wavelength as a function of
frequency and $d_{c}$ is the fixed inter-antenna spacing, that is typically
set as $\mathit{d}_{c}=\frac{c}{2\times f_{c}}$. In general, the transceiver
is designed based on the central carrier frequency $f_{c}$, which is an
approximation of the frequency over the dynamic frequency range. Therefore,
the approximation $\mathit{d}_{c}\approx\frac{\lambda}{2}$ only holds for a
narrowband channel, i.e., when the bandwidth B obeys
$\mathit{B}\ll\mathit{f}_{\mathit{c}}$. However, this condition does not hold
for a wideband signal spanning a wide frequency range. Let
$\triangle\varphi_{\mathit{n},k}$ denote the spatial phase offset for the
$\mathit{n}$th antenna with the AoA/AoD equal to $\theta$ at subcarrier $k$.
Let $\mathit{f}_{\mathit{k}}$ represent the frequency of the $\mathit{k}$-th
subcarrier that is defined as
$\mathit{f}_{\mathit{k}}\triangleq\mathit{f}_{\mathit{c}}+\left(\mathit{k}-\frac{\mathit{K}+1}{2}\right)\frac{\mathit{B}}{\mathit{K}}.$
From (7), the expression for $\triangle\varphi_{\mathit{n,k}}$ is given by

$\displaystyle\triangle\varphi_{\mathit{n,k}}$
$\displaystyle=\frac{2\pi}{\lambda_{k}}(\mathit{n}-1)\mathit{d}_{c}\cos\theta=\frac{2\pi}{\frac{c}{f_{k}}}(n-1)\frac{c}{f_{c}}\times\frac{1}{2}\cos\theta$
$\displaystyle=\pi(\mathit{n}-1)\delta_{k}\cos\theta,$ (8)

where $\delta_{k}=\frac{\mathit{f}_{k}}{\mathit{f}_{\mathit{c}}}$ represents
the relative frequency of the subcarrier. Thus, the effective spatial AoA
$\hat{\theta}(f_{k})$ at subcarrier frequency $f_{k}$ is given by

$\displaystyle\hat{\theta}(\mathit{f_{k}})=\textrm{arccos}(\delta_{k}\cos\theta).$
(9)

Note that the spatial phase offset $\triangle\varphi_{n,k}$ varies by the
subcarrier index as seen in Equation (8), which implies that the effective
spatial AoA/ AoD $\hat{\theta}(f_{k})$ as defined in Equation (9) also varies
with respect to the subcarrier index $k$. Thus, the spatial phase offset
depends on the subcarrier frequency. Due to the large antenna aperture in the
THz band, the propagation delay between the first and last antenna elements is
non-negligible, hence the transmit symbol reaches different antennas with
different delays. This is termed the spatial-wideband effect. This spatial
wideband effect is aggravated in the THz band due to having numerous antenna
elements, coupled with the large bandwidth. By the same token, the variation
in the frequency across the subcarriers leads to variations in the effective
AoA/ AoD across the frequency band, which is known as the beam squint effect.
Thus, the beam squint effect is the manifestation of the spatial wideband
effect in the frequency domain and therefore, the two phenomena are analogous
[22], [29]. Moreover, the magnitude of the beam squint is directly
proportional to the extent of spatial wideband effect occurring in such
systems, as discussed in Theorem-$1$ of [23]. Note that a majority of the
literature ignores the above effect. Additionally, the THz channel experiences
the frequency wideband effect as a consequence of the multipath delay spread.
Therefore, a combination of spatial and frequency wideband effects in THz
hybrid MIMO systems is termed the dual-wideband effect. From (9), one can
further see that the effective spatial AoA/ AoD $\hat{\theta}(f_{k})$ varies
with respect to the subcarrier frequency, which depends on the relative
frequency $\delta_{k}$. Hence, the normalized array gain $G(\theta,f_{k})$
[24] for a wideband THz system for the $k$th subcarrier and angular direction
$\theta$ is
$G(\theta,f_{k})=|\tilde{\mathbf{a}}^{H}(\theta,f_{c})\tilde{\mathbf{a}}(\theta,f_{k})|^{2}.$
The next section develops a model for the wideband THz channel suffering from
the dual-wideband effect, which arises due to a combination of the spatial-
wideband and frequency-wideband effects.

###  III-B THz channel model with dual wideband effect

The array steering vector of the wideband THz channel exhibiting beam squint,
characterized by frequency dependent angles, is given by

$\displaystyle\begin{split}\tilde{\mathbf{a}}(\theta,\mathit{f_{k}})=\frac{1}{\sqrt{N}}\left[1,e^{-j\pi\delta_{k}\cos(\theta)},\cdots,e^{{-j(\mathit{N}-1)\pi\delta_{k}\cos(\theta)}}\right]^{T}.\end{split}$
(10)

Thus, the dual wideband channel at the $\mathit{k}$th subcarrier, where
$\mathit{k}=1,\cdots,\mathit{K}$ can be modeled as

$\displaystyle\mathbf{H}[{\mathit{k}}]=\mathbf{H}^{\textrm{LoS}}[{\mathit{k}}]+\mathbf{H}^{\textrm{NLoS}}[{\mathit{k}}],$
(11)

where the LoS and NLoS components $\mathbf{H}^{\mathrm{LoS}}[k]$,
$\mathbf{H}^{\mathrm{NLoS}}[k]$ are defined as

$\displaystyle\mathbf{H}^{\text{LoS}}[k]=\sqrt{N_{T}N_{R}}\alpha_{\text{L}}(f_{k},d)\times\beta_{\tau_{\text{L}}}\times$
$\displaystyle A_{T}A_{R}\tilde{\mathbf{a}}_{R}(\phi_{\text{L}}^{R},f_{k})$
$\displaystyle\tilde{\mathbf{a}}_{T}^{H}(\phi_{\text{L}}^{T},f_{k}),$ (12)
$\displaystyle\mathbf{H}^{\textrm{NLoS}}[{\mathit{k}}]=$
$\displaystyle{\sqrt{\frac{\mathit{N}_{\mathit{T}}\mathit{N}_{\mathit{R}}}{\mathit{N}_{\textrm{NLoS}}\mathit{N}_{\mathit{ray}}}}}\sum_{\mathit{z}=1}^{\mathit{N}_{\textrm{NLoS}}}\sum_{\mathit{j}=1}^{\mathit{N}_{\mathit{ray}}}\alpha_{\mathit{z},\mathit{j}}(\mathit{f_{\mathit{k}}},\mathit{d})\times\beta_{\tau_{\mathit{z},\mathit{j}}}\times$
$\displaystyle\mathit{A}_{\mathit{T}}\mathit{A}_{\mathit{R}}\tilde{\mathbf{a}}_{\mathit{R}}(\phi_{\mathit{z},\mathit{j}}^{\mathit{R}},\mathit{f_{\mathit{k}}})\tilde{\mathbf{a}}_{\mathit{T}}^{\mathit{H}}(\phi_{\mathit{z},\mathit{j}}^{\mathit{T}},\mathit{f_{\mathit{k}}}).$
(13)

The quantity $\beta_{\tau_{z,j}}$ is defined as
$\beta_{\tau_{\mathit{z},\mathit{j}}}=\sum_{\mathit{l}=0}^{\mathit{K}-1}\mathit{p}(\mathit{l}\mathit{T}_{\mathit{s}}-\tau_{\mathit{z},\mathit{j}})e^{-j\frac{2\pi\mathit{k}\mathit{l}}{\mathit{K}}},\forall\mathit{k},\mathit{l}$.
The quantities $N_{R}$ and $N_{T}$ represent the number of RAs and TAs,
$N_{ray}$ is the number of diffused rays from a scatterer, while
$N_{\text{NLoS}}$ represents the number of clusters/ dominant scatterers
present in the propagation environment. Note that, here a cluster represents a
group of rays that arrive at the receiver around the same time and from the
same direction [30], which arise due to reflection and scattering from a
dominant scatterer. The pulse shaping filter is represented by $p(.)$ and the
sampling time by $T_{s}$, whereas $A_{R}$ and $A_{T}$ denote the RA and TA
gains. The quantities $\tau_{(.)}$, $\phi^{R}_{(.)}$, $\phi^{T}_{(.)}$,
$\alpha(.)$ represent the delay, AoA, AoD and complex path gain for the $l$th
delay tap, the models for which are shown in the next section.

####  III-B1 Calculation of the complex-path gains

The typical THz wireless channel incurs a high loss owing to various factors
such as molecular absorption losses and molecular noise due to water vapour,
among others. These aspects can be incorporated in the modeling of the
complex-path gains $\alpha_{L}(f_{k},d)$, $\alpha_{z,j}(f_{k},d)$, which are a
function of both the subcarrier frequency $f_{k}$ and transmission distance
$d$. The complex path gain $\alpha_{(.)}(f_{k},d)$ can be characterized by its
magnitude and angle components as
$\alpha_{(.)}(f_{k},d)=|\alpha_{(.)}(f_{k},d)|e^{j\omega}$. According to [14],
the path-gain of the LoS component can be modeled as

$\displaystyle|\alpha_{\textrm{L}}(\mathit{f}_{\mathit{k}},\mathit{d})|^{2}=\mathit{L}_{\textrm{spread}}(\mathit{f}_{\mathit{k}},\mathit{d})\mathit{L}_{\textrm{abs}}(\mathit{f}_{\mathit{k}},\mathit{d}),$
(14)

where $\mathit{L}_{\textrm{spread}}(\mathit{f}_{\mathit{k}},\mathit{d})$
represents the free-space loss, while
$\mathit{L}_{\textrm{abs}}(\mathit{f}_{\mathit{k}},\mathit{d})$ represents the
molecular absorption loss [24], at frequency $f_{k}$, which are respectively
given by

$\displaystyle\mathit{L}_{\textrm{spread}}(\mathit{f}_{\mathit{k}},\mathit{d})=\left(\frac{\mathit{c}}{4\pi\mathit{f}_{\mathit{k}}\mathit{d}}\right)^{2},\mathit{L}_{\textrm{abs}}(\mathit{f}_{\mathit{k}},\mathit{d})=e^{-\mathit{k}_{\textrm{abs}}(\mathit{f}_{\mathit{k}})\mathit{d}}.$
(15)

The quantity $k_{\mathrm{abs}}(\mathit{f}_{\mathit{k}})$ represents the
molecular absorption coefficient that is given
$\mathit{k}_{\textrm{abs}}(\mathit{f}_{\mathit{k}})=\sum_{\mathit{a},\mathit{g}}\mathit{k}_{\textrm{abs}}^{\mathit{a},\mathit{g}}(\mathit{f}_{\mathit{k}})$,
where
$\mathit{k}_{\textrm{abs}}^{\mathit{a},\mathit{g}}(\mathit{f}_{\mathit{k}})$
represents the measure of radiation absorption by the $a$th isotopologue of
gas $g$. The parameter
$\mathit{k}_{\textrm{abs}}^{\mathit{a},\mathit{g}}(\mathit{f}_{\mathit{k}})$
can be calculated from the HITRAN database given in [14]. Similarly, for the
complex path gain of the $\mathit{z}$th NLoS component, the magnitude of the
$\mathit{j}$th diffuse-ray is mathematically given by

$\displaystyle|\alpha_{\mathit{z},\mathit{j}}(\mathit{f}_{\mathit{k}},\mathit{d})|^{2}=\Gamma^{2}_{\mathit{z,\mathit{j}}}(\mathit{f}_{\mathit{k}})\mathit{L}_{\textrm{spread}}(\mathit{f}_{\mathit{k}},\mathit{d})\mathit{L}_{\textrm{abs}}(\mathit{f}_{\mathit{k}},\mathit{d}),$
(16)

where $\Gamma_{\mathit{z,\mathit{j}}}(\mathit{f}_{\mathit{k}})$ represents the
first-order reflection coefficient for the $\mathit{j}$th diffuse ray of the
$z$th NLoS cluster [6], which is defined as the product of the Fresnel
Reflection Coefficient $\gamma_{z,j}(\mathit{f}_{\mathit{k}})$ and of the
Rayleigh roughness factor $\varrho_{z,j}(\mathit{f}_{\mathit{k}})$ [31], that
are respectively given by

$\Gamma_{z,j}(\mathit{f}_{\mathit{k}})=\gamma_{z,j}(\mathit{f}_{\mathit{k}})\varrho_{z,j}(\mathit{f}_{\mathit{k}}),$
$\gamma_{z,j}(\mathit{f}_{\mathit{k}})=\frac{Z(\mathit{f}_{\mathit{k}})\cos(\theta_{\text{in}_{z,j}})-Z_{0}\cos(\theta_{\text{ref}_{z,j}})}{Z(\mathit{f}_{\mathit{k}})\cos(\theta_{\text{in}_{z,j}})+Z_{0}\cos(\theta_{\text{ref}_{z,j}})},$
$\varrho_{z,j}(\mathit{f}_{\mathit{k}})=e^{-\frac{1}{2}\left(\frac{4\pi\mathit{f}\sigma\cos\left(\theta_{\text{in}_{z,j}}\right)}{\mathit{c}}\right)^{2}},$
(17)

where $\theta_{\text{in}_{z,j}}$ denotes the angle of incidence,
$\theta_{\text{ref}_{z,j}}=\sin^{-1}\left(\sin(\theta_{\text{in}_{z,j}})\frac{\mathit{Z}(\mathit{f}_{\mathit{k}})}{\mathit{Z}_{0}}\right)$
represents the angle of refraction, $\mathit{Z}_{0}=377\Omega$ refers to the
intrinsic impedance of free space and $Z(\mathit{f}_{\mathit{k}})$ represents
the characteristic impedance of the reflecting medium given [31] by
$Z(\mathit{f}_{\mathit{k}})=\sqrt{\frac{\mu_{0}}{\varepsilon_{0}\left(\mathit{n}^{2}-(\frac{\varsigma\mathit{c}}{4\pi\mathit{f}_{\mathit{k}}})^{2}-\mathit{j}\frac{2\mathit{n}\varsigma\mathit{c}}{4\pi\mathit{f}_{\mathit{k}}}\right)}}.$
The constants $\mu_{0}$ and $\varepsilon_{0}$ represent the free-space
permeability and permittivity respectively, $\mathit{n}$ denotes the index of
refraction and $\varsigma$ represents the absorption coefficient of the
reflecting medium. Note that, the traditional CSI estimation techniques such
as LS/ MMSE cannot exploit the spatial-sparsity of the THz MIMO channel that
occurs due to the higher link directionality. The virtual channel model for
the THz system, that enables improved estimation, is described in the
subsequent section.

##  IV Sparse Channel Estimation Model for an SC-FDE THz System

As discussed in Section II-A, the output vector at the RC for the
$\mathit{m}$th block at the $\mathit{k}$th subcarrier is given by

$\displaystyle\mathbf{{y}}_{\mathit{m}}[\mathit{k}]=\mathbf{\Phi}_{\mathit{m}}[\mathit{k}]\mathbf{h}[\mathit{k}]+\mathbf{{v}}_{\mathit{m}}[\mathit{k}].$
(18)

To derive an equivalent model [32] for the $\mathit{k}$th subcarrier, we
vertically stack the outputs $\mathbf{y}_{\mathit{m}}[\mathit{k}]$ for all the
blocks, i.e, $1\leq\mathit{m}\leq\textit{M}$

$\displaystyle\underbrace{\begin{bmatrix}{\mathbf{y}}_{\text{1}}[k]\\\
\vdots\\\
{\mathbf{y}}_{\mathit{M}}[k]\end{bmatrix}}_{{\mathbf{y}_{\mathit{p}}}[k]}=\underbrace{\begin{bmatrix}\mathbf{\Phi}_{\text{1}}[\mathit{k}]\\\
\vdots\\\
\mathbf{\Phi}_{\mathit{M}}[\mathit{k}]\end{bmatrix}}_{\mathbf{\Phi_{\mathit{p}}[\mathit{k}]}}\mathbf{h}[\mathit{k}]+\underbrace{\begin{bmatrix}\mathbf{v}_{\text{1}}[k]\\\
\vdots\\\
\mathbf{v}_{\mathit{M}}[k]\end{bmatrix}}_{\mathbf{v}_{\mathit{p}}[k]},$ (19)

where,
$\mathbf{y}_{\mathit{p}}[\mathit{k}]\in\mathbb{C}^{\mathit{M}\mathit{N}_{\mathit{RF}}\times
1}$ represents the stacked pilot output vector,
$\mathbf{\Phi}_{\mathit{p}}[\mathit{k}]\in\mathbb{C}^{\mathit{M}\mathit{N}_{\mathit{RF}}\times\mathit{N}_{\mathit{R}}\mathit{N}_{\mathit{T}}}$
denotes the stacked pilot sensing matrix and
$\mathbf{v}_{\mathit{p}}[\mathit{k}]\in\mathbb{C}^{\mathit{M}\mathit{N}_{\mathit{RF}}\times
1}$ represents the stacked noise corresponding to all the $\mathit{M}$ blocks
for subcarrier $k$. Thus, the stacked system model at the combiner output is
given by

$\displaystyle\mathbf{y}_{\mathit{p}}[\mathit{k}]=\mathbf{\Phi}_{\mathit{p}}[\mathit{k}]\mathbf{h}[\mathit{k}]+\mathbf{v}_{\mathit{p}}[\mathit{k}].$
(20)

At this point, one can formulate the traditional ZF and MMSE estimators for
the wideband THz MIMO channel at each subcarrier as

$\displaystyle\hat{\mathbf{h}}_{\textrm{LS}}[\mathit{k}]=(\mathbf{\Phi}_{\mathit{p}}^{\mathit{H}}[\mathit{k}]\mathbf{\Phi}_{\mathit{p}}[\mathit{k}])^{-1}\mathbf{\Phi}_{\mathit{p}}^{\mathit{H}}[\mathit{k}]\mathbf{y}_{\mathit{p}}[\mathit{k}]$
(21)
$\displaystyle\widehat{\mathbf{h}}_{\textrm{MMSE}}[\mathit{k}]=(\mathbf{R}_{\mathit{h}}^{-\textrm{1}}[\mathit{k}]+\mathbf{\Phi}_{\mathit{p}}^{\mathit{H}}[\mathit{k}]\mathbf{R}_{\mathit{v}_{\mathit{p}}}^{-\textrm{1}}\mathbf{\Phi}_{\mathit{p}}[\mathit{k}])^{-1}\mathbf{\Phi}_{\mathit{p}}^{\mathit{H}}[\mathit{k}]\mathbf{y}_{\mathit{p}}[\mathit{k}],$
(22)

where, the noise covariance matrix
$\mathbf{R}_{\mathit{v}_{\mathit{p}}}=\sigma^{2}_{v}K\>\text{blkdiag}\left(\left\\{\mathbf{W}^{H}_{RF,m}\mathbf{W}_{RF,m}\right\\}_{m=1}^{M}\right)$
and
$\mathbf{R}_{\mathit{v}_{\mathit{p}}}\in\mathbb{C}^{\mathit{M}\mathit{N}_{RF}\times\mathit{M}\mathit{N}_{RF}}$.
Furthermore,
$\mathbf{R}_{\mathit{h}}[\mathit{k}]=\mathbb{E}[\mathbf{h}[\mathit{k}]\mathbf{h}^{\mathit{H}}[\mathit{k}]]\in\mathbb{C}^{\mathit{N}_{R}\mathit{N}_{T}\times\mathit{N}_{R}\mathit{N}_{T}}$
denotes the channel covariance matrix. Note that the conventional LS and MMSE-
based estimators require excessive pilot overheads in THz MIMO systems since
they require an over-determined model, i.e,
$\mathit{M}\mathit{N}_{\mathit{RF}}\geq\mathit{N}_{\mathit{T}}\mathit{N}_{\mathit{R}}$.
This substantially degrades the spectral efficiency. Moreover, these orthodox
techniques listed above do not leverage the specific properties of the THz
channel. An important observation in this context is the fact that the
directional nature of propagation in the THz band leads to an angularly sparse
multipath channel [6]. In such scenarios, the recently developed sparse signal
processing paradigm can yield excellent results for channel estimation in THz
band hybrid MIMO systems, even in ill-posed scenarios, wherein we have
$\mathit{M}\mathit{N}_{\mathit{RF}}\ll\mathit{N}_{\mathit{T}}\mathit{N}_{\mathit{R}}$.
This feature may lead to substantial enhancement in the accuracy of channel
estimation as well as bandwidth-efficiency. Therefore, the following section
presents a sparse channel estimation model in THz band communication.

###  IV-A Sparse channel estimation model for a wideband THz hybrid MIMO
system

Let $\mathit{G}_{\mathit{T}}$ and $\mathit{G}_{\mathit{R}}$ represent the
angular grid-sizes, which obey the relationship
$(\mathit{G}_{\mathit{T}},\mathit{G}_{\mathit{R}})\geq\textrm{max}(\mathit{N}_{\mathit{T}},\mathit{N}_{\mathit{R}})$.
The transmit and receive angular grids $\Theta_{T}$ and $\Theta_{R}$,
respectively, are created by considering the directional-cosines to be
uniformly spaced between $-1$ to $1$, mathematically described as

$\displaystyle\Theta_{\mathit{T}}=\left\\{\theta_{\mathit{t}}:\cos(\theta_{\mathit{t}})=\frac{\textrm{2}}{\mathit{G}_{\mathit{T}}}(\mathit{t}-\textrm{1})-1,1\leq\mathit{t}\leq\mathit{G}_{\mathit{T}}\right\\},$
(23)
$\displaystyle\Theta_{\mathit{R}}=\left\\{\theta_{\mathit{r}}:\cos(\theta_{\mathit{r}})=\frac{\textrm{2}}{\mathit{G}_{\mathit{R}}}(\mathit{r}-\textrm{1})-1,1\leq\mathit{r}\leq\mathit{G}_{\mathit{R}}\right\\}.$
(24)

Let
$\mathbf{A}_{\mathit{T}}(\Theta_{\mathit{T}},\mathit{f}_{k})\in\mathbb{C}^{\mathit{N}_{\mathit{T}}\times\mathit{G}_{\mathit{T}}}$
and
$\mathbf{A}_{\mathit{R}}(\Theta_{\mathit{R}},\mathit{f}_{k})\in\mathbb{C}^{\mathit{N}_{\mathit{R}}\times\mathit{G}_{\mathit{R}}}$
represent the dictionary matrices for the transmit and receive array manifold
vectors using the angular grids $\Theta_{\mathit{T}}$ and
$\Theta_{\mathit{R}}$ [22], which are modeled as

$\displaystyle\mathbf{A}_{\mathit{T}}(\Theta_{\mathit{T}},\mathit{f}_{k})=[\tilde{\mathbf{a}}_{\mathit{T}}(\theta_{\textrm{1}},\mathit{f}_{k}),\tilde{\mathbf{a}}_{\mathit{T}}(\theta_{\textrm{2}},\mathit{f}_{k}),\cdots,\tilde{\mathbf{a}}_{\mathit{T}}(\theta_{\mathit{G}_{\mathit{T}}},\mathit{f}_{k})],$
(25)
$\displaystyle\mathbf{A}_{\mathit{R}}(\Theta_{\mathit{R}},\mathit{f}_{k})=[\tilde{\mathbf{a}}_{\mathit{R}}(\theta_{\textrm{1}},\mathit{f}_{k}),\tilde{\mathbf{a}}_{\mathit{R}}(\theta_{\textrm{2}},\mathit{f}_{k}),\cdots,\tilde{\mathbf{a}}_{\mathit{R}}(\theta_{\mathit{G}_{\mathit{R}}},\mathit{f}_{k})].$
(26)

Because of the particular selection of the AoD and AoA grids, the array
response dictionary matrices [18] satify the conditions

$\displaystyle\mathbf{A}_{\mathit{T}}(\Theta_{\mathit{T}},\mathit{f}_{k})\mathbf{A}^{\mathit{H}}_{\mathit{T}}(\Theta_{\mathit{T}},\mathit{f}_{k})=\frac{\mathit{G_{\mathit{T}}}}{\mathit{N}_{\mathit{T}}}\mathbf{I}_{\mathit{N}_{\mathit{T}}},$
$\displaystyle\mathbf{A}_{\mathit{R}}(\Theta_{\mathit{R}},\mathit{f}_{k})\mathbf{A}^{\mathit{H}}_{\mathit{R}}(\Theta_{\mathit{R}},\mathit{f}_{k})=\frac{\mathit{G_{\mathit{R}}}}{\mathit{N}_{\mathit{R}}}\mathbf{I}_{\mathit{N}_{\mathit{R}}}.$
(27)

It is noteworthy that the number of multipath components is substantially
lower in a THz system due to the reflections, molecular absorption and free-
space losses, as described earlier. This fact, coupled with the presence of
only a few highly-directional beams, renders the THz MIMO channel
$\mathbf{H}_{\mathit{b}}[k]$ angularly sparse in nature. Subsequently, the THz
MIMO channel $\mathbf{H}[\mathit{k}]$ can be represented using the virtual
channel model [32] as

$\displaystyle\mathbf{H}[\mathit{k}]=\mathbf{A}_{\mathit{R}}(\Theta_{\mathit{R}},\mathit{f}_{\mathit{k}}){\mathbf{H}_{\mathit{b}}[\mathit{k}]}\mathbf{A}_{\mathit{T}}^{\mathit{H}}(\Theta_{\mathit{T}},\mathit{f}_{\mathit{k}}),$
(28)

where
$\mathbf{H}_{\mathit{b}}[\mathit{k}]\in\mathbb{C}^{\mathit{G}_{\mathit{R}}\times\mathit{G}_{\mathit{T}}}$
is the beamspace CFR matrix that pertains to $\mathbf{H}[\mathit{k}]$.
Exploiting the property of the vec(.) operator described in Section-I-C, the
vectorized THz MIMO channel can be represented as

$\displaystyle\mathbf{h}[\mathit{k}]=\textrm{vec}(\mathbf{H}[\mathit{k}])=\left[\mathbf{A}^{*}_{\mathit{T}}(\Theta_{\mathit{T}},\mathit{f}_{\mathit{k}})\otimes\mathbf{A}_{\mathit{R}}(\Theta_{\mathit{R}},\mathit{f}_{\mathit{k}})\right]\mathbf{h}_{\mathit{b}}[\mathit{k}],$
(29)

where
$\mathbf{h}_{\mathit{b}}[\mathit{k}]=\textrm{vec}(\mathbf{H}_{\mathit{b}}[\mathit{k}])\in\mathbb{C}^{\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}\times
1}$. Let
$\mathbf{\Psi}[\mathit{k}]=\left[\mathbf{A}^{*}_{\mathit{T}}(\Theta_{\mathit{T}},\mathit{f}_{\mathit{k}})\otimes\mathbf{A}_{\mathit{R}}(\Theta_{\mathit{R}},\mathit{f}_{\mathit{k}})\right]$
denote the sparsifying-dictionary, where
$\mathbf{\Psi}[\mathit{k}]\in\mathbb{C}^{\mathit{N}_{\mathit{R}}\mathit{N}_{\mathit{T}}\times\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}}$.
Thus, upon substituting (29) in to (20), the expression for the output vector
can be recast as

$\displaystyle\mathbf{y}_{\mathit{p}}[\mathit{k}]=\underbrace{\mathbf{\Phi}_{\mathit{p}}[\mathit{k}]\mathbf{\Psi}[\mathit{k}]}_{\tilde{\mathbf{\Psi}}_{\mathit{p}}[\mathit{k}]}\mathbf{h}_{\mathit{b}}[k]+\mathbf{v}_{\mathit{p}}[\mathit{k}],$
(30)

where
$\tilde{\mathbf{\Psi}}_{\mathit{p}}[\mathit{k}]\in\mathbb{C}^{\mathit{M}\mathit{N}_{\mathit{RF}}\times\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}}$
represents an equivalent sensing matrix. It is noteworthy that other sparse
estimation techniques, such as Basis-Pursuit [33], FOCUSS [34], and OMP [35],
suffer from structural and convergence-related issues. For instance, in the
case of Basis-Pursuit, there is no guarantee that the sparsest solution aligns
with the global minimum of its cost function, leading to structural errors.
Similarly, the FOCUSS algorithm frequently converges to suboptimal local
minima, resulting in convergence errors. Additionally, the OMP algorithm
presents challenges; a low cessation criterion leads to a larger number of
iterations and more non-zero components in the sparse vector, causing
structural errors. Conversely, a high cessation criterion results in fewer
iterations, failing to capture all the dominant modes of the channel and
causing convergence errors. By contrast, the proposed Bayesian learning
methods assure global convergence with a high probability due to the inherent
attributes of their cost function and the use of the EM algorithm.
Furthermore, in a noiseless scenario, the global minimum of the Bayesian
learning cost function corresponds to the sparsest solution [36]. Even in a
noisy scenario, all the local minima are sparse and have improved
representation over non-sparse solutions. The next section describes the
framework that utilizes Bayesian learning to estimate $\mathbf{h}_{b}[k]$,
that intends to exploit the sparsity of the THz MIMO channel as described
above.

##  V PA-BL for Sparse Channel Estimation in THz MIMO Systems

This section introduces a Pilot Aided-Bayesian Learning (PA-BL) technique for
CSI estimation in dual-wideband channel. This can also be applied for the case
of an under-determined system, i.e,
$\mathit{M}\mathit{N}_{\mathit{RF}}\ll\mathit{G}_{\mathit{T}}\mathit{G}_{\mathit{R}}$
[18]. At first, the parameterized Gaussian prior
$\mathit{f}(\mathbf{h}_{\mathit{b}}[\mathit{k}];\mathbf{\Gamma}_{\mathit{k}})$
is allocated to the $\mathit{k}$th pilot subcarrier for the sparse vector
$\mathbf{h}_{\mathit{b}}[\mathit{k}]$ [36], which is given by

$\displaystyle\mathit{f}(\mathbf{h}_{\mathit{b}}[\mathit{k}];\mathbf{\Gamma}_{\mathit{k}})=\prod_{\mathit{i}=1}^{\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}}(\pi\gamma_{\mathit{k},\mathit{i}})^{-1}\textrm{exp}\left(-\frac{|\mathbf{h}_{\mathit{b}}[\mathit{k}](\mathit{i})|^{2}}{\gamma_{\mathit{k},\mathit{i}}}\right),$
(31)

where $\gamma_{\mathit{k},\mathit{i}}$ represents $i$th hyper-parameter for
the $\mathit{k}$th subcarrier;
$1\leq\mathit{i}\leq\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}$, and the
matrix
$\mathbf{\Gamma}_{\mathit{k}}\in\mathbb{R}^{\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}\times\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}}$
is the diagonal matrix of hyper-parameters is given by
$\mathbf{\Gamma}_{\mathit{k}}=\textrm{diag}(\gamma_{\mathit{k},\textrm{1}},\gamma_{\mathit{k},\textrm{2}},\cdots,\gamma_{\mathit{k},\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}}).$
The MMSE estimate of $\hat{\mathbf{h}}_{\mathit{b}}[\mathit{k}]$ is given by

$\displaystyle\hat{\mathbf{h}}_{\mathit{b}}[\mathit{k}]=\left(\tilde{\mathbf{\Psi}}_{\mathit{p}}^{\mathit{H}}[\mathit{k}]\mathbf{R}_{\mathit{v}_{\mathit{p}}}^{-1}\tilde{\mathbf{\Psi}}_{\mathit{p}}[\mathit{k}]+\mathbf{\Gamma}_{\mathit{k}}^{-1}\right)^{-1}\tilde{\mathbf{\Psi}}_{\mathit{p}}^{\mathit{H}}[\mathit{k}]\mathbf{R}_{\mathit{v}_{\mathit{p}}}^{-1}\mathbf{y}_{\mathit{p}}[\mathit{k}].$
(32)

It can be noted that MMSE relies on the hyper-parameter matrix
$\mathbf{\Gamma}_{\mathit{k}}$. Therefore, estimating
$\mathbf{\Gamma}_{\mathit{k}}$ holds the key toward determining the sparse
channel $\hat{\mathbf{h}}_{\mathit{b}}[\mathit{k}]$. The log-likelihood
$\textrm{log}\left[\mathit{f}\left(\mathbf{y}_{\mathit{p}}[\mathit{k}];\mathbf{\Gamma}_{\mathit{k}}\right)\right]$
of the hyper-parameter matrix $\mathbf{\Gamma}_{\mathit{k}}$ can be formulated
as

$\displaystyle\textrm{log}\left[\mathit{f}\left(\mathbf{y}_{\mathit{p}}[\mathit{k}];\mathbf{\Gamma}_{\mathit{k}}\right)\right]=\mathit{c}_{\textrm{1}}-$
$\displaystyle\textrm{log}\left[\textrm{det}\left(\mathbf{R}_{\mathbf{y}_{\mathit{p}}}[\mathit{k}]\right)\right]-$
$\displaystyle\mathbf{y}_{\mathit{p}}^{\mathit{H}}[\mathit{k}]\mathbf{R}_{\mathbf{y}_{\mathit{p}}}^{\textrm{-1}}[\mathit{k}]\mathbf{y}_{\mathit{p}}[\mathit{k}],$
(33)

where
$\mathit{c}_{\textrm{1}}=-\mathit{M}\mathit{N}_{\mathit{RF}}\textrm{log}(\pi)$
and the pilot covariance matrix
$\mathbf{R}_{\mathbf{y}_{\mathit{p}}}[\mathit{k}]\in\mathbb{C}^{\mathit{M}\mathit{N}_{\mathit{RF}}\times\mathit{M}\mathit{N}_{\mathit{RF}}}$
is given by

$\displaystyle\mathbf{R}_{\mathbf{y}_{\mathit{p}}}[\mathit{k}]=\mathbb{E}\left\\{\mathbf{y}_{\mathit{p}}[\mathit{k}]\mathbf{y}_{\mathit{p}}^{\mathit{H}}[\mathit{k}]\right\\}.$
(34)

Substituting equation (30) in to (34) yields the pilot covariance matrix in
the form of

$\displaystyle\mathbf{R}_{\mathbf{y}_{\mathit{p}}}[\mathit{k}]=\mathbf{R}_{\mathit{v}_{\mathit{p}}}+\tilde{\boldsymbol{\Psi}}_{\mathit{p}}[\mathit{k}]\mathbf{\Gamma}_{\mathit{k}}\tilde{\boldsymbol{\Psi}}_{\mathit{p}}^{\mathit{H}}[\mathit{k}],$
(35)

where
$\mathbf{v}_{\mathit{p}}[\mathit{k}]\sim\mathcal{CN}(\mathbf{0},\mathbf{R}_{\mathit{v}_{\mathit{p}}})$.
Maximization of the log-likelihood with respect to
$\mathbf{\Gamma}_{\mathit{k}}$ is mathematically intractable. Thus, the
Expectation Maximization (EM) algorithm is a suitable choice for iterative
maximization, as it ensures convergence to a local optimum [37]. The EM-based
procedure of determining the sparse channel estimate is described below. Let
the estimate of the $\mathit{i}$th hyperparameter at the $\mathit{k}$th
subcarrier obtained from the $(\mathit{j}-1)$ st EM iteration be represented
as $\hat{\gamma}_{\mathit{k},\mathit{i}}^{\left(\mathit{j}-1\right)}$ and
$\hat{\mathbf{\Gamma}}_{\mathit{k}}^{(\mathit{j}-1)}$ denotes the
hyperparameter matrix at the $\mathit{k}$th subcarrier, which is further
defined as
$\hat{\mathbf{\Gamma}}_{\mathit{k}}^{(\mathit{j}-1)}=\textrm{diag}\left(\hat{\gamma}_{\mathit{k},1}^{(\mathit{j}-1)},\hat{\gamma}_{\mathit{k},2}^{(\mathit{j}-1)},\cdots,\hat{\gamma}_{\mathit{k},\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}}^{(\mathit{j}-1)}\right).$
The steps involved in updating the estimate
$\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}\right)}$ in the
$\mathit{j}$th iteration of the EM algorithm can be outlined as follows. Let
the $\mathit{i}$th hyper-parameter estimate obtained in iteration $(j-1)$st
for the $\mathit{k}$th subcarrier be denoted by
$\hat{\gamma}_{\mathit{k},\mathit{i}}^{(\mathit{j}-1)}$. The hyper-parameter
update $\hat{\gamma}_{\mathit{k},\mathit{i}}^{(\mathit{j})}$ in the
$\mathit{j}$th EM-iteration is obtained by maximizing the log-likelihood
function
$\mathcal{L}\left(\mathbf{\Gamma}_{\mathit{k}}|\hat{\mathbf{\Gamma}}_{\mathit{k}}^{(\mathit{j}-1)}\right)$,
defined as

$\displaystyle\mathcal{L}\left(\mathbf{\Gamma}_{\mathit{k}}|\hat{\mathbf{\Gamma}}_{\mathit{k}}^{(\mathit{j}-1)}\right)=\mathbb{E}_{\mathbf{h}_{\mathit{b}}[\mathit{k}]\mathbf{|}\mathbf{y}_{\mathit{p}}[\mathit{k}];\hat{\mathbf{\Gamma}}_{\mathit{k}}^{(\mathit{j}-1)}}\left\\{{\mathrm{log}\
\mathit{f}\left(\mathbf{y}_{\mathit{p}}[\mathit{k}],\mathbf{h}_{\mathit{b}}[\mathit{k}];\mathbf{\Gamma_{\mathit{k}}}\right)}\right\\}.$
(36)

The conditional expectation of the log-likelihood function
$\mathcal{L}\left(\mathbf{\Gamma}_{\mathit{k}}|\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}-1\right)}\right)$
in the E-step is given

$\displaystyle\mathcal{L}\left(\mathbf{\Gamma}_{\mathit{k}}|\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}-1\right)}\right)=$
$\displaystyle\mathbb{E}_{\mathbf{h}_{\mathit{b}}[\mathit{k}]\mathbf{|}\mathbf{y}_{\mathit{p}}[\mathit{k}];\hat{\mathbf{\Gamma}}_{\mathit{k}}^{(\mathit{j}-1)}}\left\\{{\textrm{log}\
\mathit{f}\left(\mathbf{y}_{\mathit{p}}[\mathit{k}],\mathbf{h}_{\mathit{b}}[\mathit{k}];\mathbf{\Gamma_{\mathit{k}}}\right)}\right\\}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\mathbf{h}_{\mathit{b}}[\mathit{k}]\mathbf{|}\mathbf{y}_{\mathit{p}}[\mathit{k}];\hat{\mathbf{\Gamma}}_{\mathit{k}}^{(\mathit{j}-1)}}\left\\{\textrm{log}\left[\mathit{f}(\mathbf{y}_{\mathit{p}}[\mathit{k}]|\mathbf{h}_{\mathit{b}}[\mathit{k}])\right]\right\\}+$
$\displaystyle\mathbb{E}_{\mathbf{h}_{\mathit{b}}[\mathit{k}]\mathbf{|}\mathbf{y}_{\mathit{p}}[\mathit{k}];\hat{\mathbf{\Gamma}}_{\mathit{k}}^{(\mathit{j}-1)}}\left\\{\textrm{log}\left[\mathit{f}(\mathbf{h}_{\mathit{b}}[\mathit{k}];\mathbf{\Gamma}_{\mathit{k}})\right]\right\\}.$
(37)

The expression enclosed within the first expectation operator in the
aforementioned equation can be simplified as

$\displaystyle\textrm{log}\
\mathit{f}(\mathbf{y}_{\mathit{p}}[\mathit{k}]|\mathbf{h}_{\mathit{b}}[\mathit{k}])=-\mathit{M}\mathit{N}_{\textrm{RF}}\
\textrm{log}(\pi)-\textrm{log}\left[\textrm{det}(\mathbf{R}_{\mathit{v}_{\mathit{p}}})\right]-$
(38)
$\displaystyle\left(\mathbf{y}_{\mathit{p}}[\mathit{k}]-\tilde{\boldsymbol{\Psi}}_{\mathit{p}}[\mathit{k}]\mathbf{h}_{\mathit{b}}[\mathit{k}]\right)^{\mathit{H}}\mathbf{R}_{\mathit{v}_{\mathit{p}}}^{-1}\left(\mathbf{y}_{\mathit{p}}[\mathit{k}]-\tilde{\boldsymbol{\Psi}}_{\mathit{p}}[\mathit{k}]\mathbf{h}_{\mathit{b}}[\mathit{k}]\right).$

It can be seen that the above expression is independent of the hyper-parameter
matrix $\mathbf{\Gamma}_{\mathit{k}}$, which implies that it can be overlooked
in the ensuing M-step of the CSI estimation algorithm. Thus, the equivalent
optimization problem of determining the beamspace CSI estimate can be
formulated

$\displaystyle\hat{\boldsymbol{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}\right)}=\underset{\mathbf{\Gamma}_{\mathit{k}}}{\textrm{arg
max}}\
\mathbb{E}_{\mathbf{h}_{\mathit{b}}[\mathit{k}]|\mathbf{y}_{\mathit{p}}[\mathit{k}];\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}-1\right)}}\left\\{\textrm{log}\left[\mathit{f}\left(\mathbf{h}_{\mathit{b}}[\mathit{k}];\boldsymbol{\Gamma}_{\mathit{k}}\right)\right]\right\\}.$
(39)

On substituting
$\mathit{f}(\mathbf{h}_{\mathit{b}}[\mathit{k}];\mathbf{\Gamma}_{\mathit{k}})$
from equation (31) into the optimization objective, one can observe that the
resultant maximization problem can be separated in reference to the hyper-
parameters, and the individual optimization problems are given as

$\displaystyle\hat{\gamma}_{\mathit{k},\mathit{i}}^{\left(\mathit{j}\right)}=\underset{\gamma_{\mathit{k},\mathit{i}}}{\textrm{arg
max}}\left[-\textrm{log}(\gamma_{\mathit{k},\mathit{i}})-\frac{\mathbb{E}_{\mathbf{h}_{\mathit{b}}[\mathit{k}]|\mathbf{y}_{\mathit{p}}[\mathit{k}];\hat{\boldsymbol{\Gamma}}_{\mathit{k}}^{(\mathit{j}-1)}}\left\\{|\mathbf{h}_{\mathit{b}}[\mathit{k}](\mathit{i})|^{2}\right\\}}{\gamma_{\mathit{k},\mathit{i}}}\right].$
(40)

On solving the aforementioned equation, the estimates
$\hat{\gamma}_{\mathit{k},\mathit{i}}^{\left(\mathit{j}\right)}$ are
determined as

$\displaystyle\hat{\gamma}_{\mathit{k},\mathit{i}}^{\left(\mathit{j}\right)}=\mathbb{E}_{\mathbf{h}_{\mathit{b}}[\mathit{k}]|\mathbf{y}_{\mathit{p}}[\mathit{k}];\hat{\boldsymbol{\Gamma}}_{\mathit{k}}^{(\mathit{j}-1)}}\left\\{|\mathbf{h}_{\mathit{b}}[\mathit{k}](\mathit{i})|^{2}\right\\},$
(41)

where, the a posteriori distribution of $\mathbf{h}_{\mathit{b}}[\mathit{k}]$,
required to simplify the conditional expectation
$\mathbb{E}_{\mathbf{h}_{b}[k]|\mathbf{y}_{p}[k];\hat{\mathbf{\Gamma}}^{(j-1)}_{k}}\left\\{.\right\\}$
above, can be determined as

$\displaystyle
f\left(\mathbf{h}_{\mathit{b}}[\mathit{k}]|\mathbf{y}_{\mathit{p}}[\mathit{k}];\hat{\boldsymbol{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}-1\right)}\right)=\mathcal{CN}\left(\boldsymbol{\mu}_{\mathit{b}}^{(\mathit{j})}[\mathit{k}],\mathbf{R}_{\mathit{b}}^{(\mathit{j})}[\mathit{k}]\right).$
(42)

The quantity
$\boldsymbol{\mu}_{\mathit{b}}^{(\mathit{j})}[\mathit{k}]\in\mathbb{C}^{\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}\times
1}$ represents the a posteriori mean, while the quantity
$\mathbf{R}_{\mathit{b}}^{(\mathit{j})}[\mathit{k}]\in\mathbb{C}^{\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}\times\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}}$
represents the covariance matrix. These expressions are given below

$\boldsymbol{\mu}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}]=\mathbf{R}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}]\tilde{\mathbf{\Psi}}_{\mathit{p}}^{\mathit{H}}[\mathit{k}]\mathbf{R}_{\mathit{v}_{\mathit{p}}}^{-\textrm{1}}\mathbf{y}_{\mathit{p}}[\mathit{k}],$
$\mathbf{R}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}]=\left[\tilde{\mathbf{\Psi}}_{\mathit{p}}^{\mathit{H}}[\mathit{k}]\mathbf{R}_{\mathit{v}_{\mathit{p}}}^{-\textrm{1}}\tilde{\mathbf{\Psi}}_{\mathit{p}}[\mathit{k}]+\left(\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}-1\right)}\right)^{-\textrm{1}}\right]^{-\textrm{1}}.$
(43)

Using the a posteriori PDF from (42), the expression in (41) for the hyper-
parameter estimate reduces to

$\displaystyle\hat{\gamma}_{\mathit{k},\mathit{i}}^{\left(\mathit{j}\right)}=\mathbf{R}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}](\mathit{i},\mathit{i})+|\boldsymbol{\mu}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}](\mathit{i})|^{2}.$
(44)

Algorithm 1 PA-BL for sparse channel estimation in THz MIMO systems

Input: Pilot output $\mathbf{y}_{p}[\mathit{k}]$, equivalent sensing matrix
$\tilde{\mathbf{\Psi}}_{p}[\mathit{k}]$, noise covariance
$\mathbf{R}_{\mathit{v_{p}}}$, dictionary matrices
$\mathbf{A}_{\mathit{R}}(\Theta_{\mathit{R}},\mathit{f}_{\mathit{k}})$ and
$\mathbf{A}_{\mathit{T}}(\Theta_{\mathit{T}},\mathit{f}_{\mathit{k}})$,
stopping parameters $\epsilon$ and $\mathit{K}_{\textrm{max}}$
Initialization:
$\gamma_{\mathit{k},\mathit{i}}^{\left(0\right)}=\textrm{1},\forall\
\textrm{1}\leq\mathit{i}\leq\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}$
where $\textrm{1}\leq\mathit{k}\leq\mathit{K}$
$\Rightarrow\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\textrm{0}\right)}=\mathbf{I}_{\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}}$,
$\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(-1\right)}=\mathbf{0}$ and counter
$\mathit{j}=\textrm{0}$
while
$\left(\parallel\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}\right)}-\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}-1\right)}\parallel_{\mathit{F}}\
>\epsilon\ \textrm{and}\ \mathit{j}<\mathit{K}_{\textrm{max}}\right)$ do

  1. 1.

$j\leftarrow j+1$

  2. 2.

E-step:  Calculate the a posteriori mean and covariance

$\mathbf{R}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}]=\left[\tilde{\mathbf{\Psi}}_{\mathit{p}}^{\mathit{H}}[\mathit{k}]\mathbf{R}_{\mathit{v}_{\mathit{p}}}^{-\textrm{1}}\tilde{\mathbf{\Psi}}_{\mathit{p}}[\mathit{k}]+\left(\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}-1\right)}\right)^{-\textrm{1}}\right]^{-\textrm{1}}$
$\boldsymbol{\mu}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}]=\mathbf{R}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}]\tilde{\mathbf{\Psi}}_{p}^{\mathit{H}}[\mathit{k}]\mathbf{R}_{\mathit{v_{p}}}^{-\textrm{1}}\mathbf{y}_{p}[\mathit{k}]$

  3. 3.

M-step: Update the hyperparameters
beta for $\mathit{i}=1,\cdots,\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}$
do

$\hat{\gamma}_{\mathit{k},\mathit{i}}^{\left(\mathit{j}\right)}=\mathbf{R}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}](\mathit{i},\mathit{i})+|\boldsymbol{\mu}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}](\mathit{i})|^{2}$

beta end for

end while
$\hat{\mathbf{h}}_{\mathit{b}}[\mathit{k}]=\boldsymbol{\mu}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}]$
Output: $\widehat{\mathbf{H}}_{\textrm{PA-
BL}}[\mathit{k}]=\mathbf{A}_{\mathit{R}}(\Theta_{\mathit{R}},\mathit{f}_{\mathit{k}})\textrm{vec}^{-1}(\widehat{\mathbf{h}}_{\mathit{b}}[\mathit{k}])\mathbf{A}_{\mathit{T}}^{\mathit{H}}(\Theta_{\mathit{T}},\mathit{f}_{\mathit{k}})$

Algorithm 1 provides an overview of the sparse channel estimation technique
for THz hybrid MIMO systems, employing the PA-BL framework [6]. The procedure
is repeated until the hyperparameter estimate converges, i.e,
$\parallel\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}\right)}-\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}-1\right)}\parallel_{\mathit{F}}^{\textrm{2}}\
<\ \epsilon_{\textrm{0}}$ or the number of iterations reached its maximum
affordable value $\mathit{K}_{\textrm{max}}$. The PA-BL based channel estimate
after convergence is given by

$\displaystyle\hat{\mathbf{h}}_{\mathit{b}}[\mathit{k}]=\boldsymbol{\mu}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}].$
(45)

###  V-A Bayesian Cramér-Rao lower bound for pilot-aided learning

Consider the parameterized Gaussian prior assigned to the $\mathit{i}$th
hyper-parameter of the $\mathit{k}$th pilot subcarrier corresponding to
$\mathbf{h}_{\mathit{b}}[\mathit{k}]$, as mentioned in (31). The corresponding
PDF is given by

$\displaystyle\mathit{f}(\mathbf{h}_{\mathit{b}}[\mathit{k}];\mathbf{\Gamma}_{\mathit{k}})=\frac{1}{(\pi)^{\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}}\
\textrm{det}(\mathbf{\Gamma}_{\mathit{k}})}\textrm{exp}\left(-\mathbf{h}_{\mathit{b}}^{\mathit{H}}[\mathit{k}]\mathbf{\Gamma}_{\mathit{k}}^{-1}\mathbf{h}_{\mathit{b}}[\mathit{k}]\right).$
(46)

The log-likelihood of the beamspace channel
$\mathbf{h}_{\mathit{b}}[\mathit{k}]$ is obtained via applying the log
function to both sides of the above equation, yielding

$\displaystyle\mathcal{L}\left(\mathit{f}\left[\mathbf{h}_{\mathit{b}}[\mathit{k}];\mathbf{\Gamma}_{\mathit{k}}\right]\right)=$
$\displaystyle-\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}\textrm{log}(\pi)-\textrm{log}\left[\textrm{det}(\mathbf{\Gamma}_{\mathit{k}})\right]-$
$\displaystyle\mathbf{h}_{\mathit{b}}^{\mathit{H}}[\mathit{k}]\mathbf{\Gamma}_{\mathit{k}}^{-1}\mathbf{h}_{\mathit{b}}[\mathit{k}].$
(47)

The conditional PDF of the pilot output vector
$\mathbf{y}_{\mathit{p}}[\mathit{k}]$, as defined in equation (30), is given
by the complex normal distribution
$\mathcal{CN}\left(\tilde{\mathbf{\Psi}}_{\mathit{p}}[\mathit{k}]\mathbf{h}_{\mathit{b}}[\mathit{k}],\mathbf{R}_{\mathit{v}_{\mathit{p}}}\right)$.
Hence, the expression of the conditional PDF is

$\displaystyle\mathit{f}\left(\mathbf{y}_{\mathit{p}}[\mathit{k}]|\mathbf{h}_{\mathit{b}}[\mathit{k}]\right)$
$\displaystyle=\vartheta[\mathit{k}]\textrm{exp}\big{(}-\left(\mathbf{y}_{\mathit{p}}[\mathit{k}]-\tilde{\boldsymbol{\Psi}}_{\mathit{p}}[\mathit{k}]\mathbf{h}_{\mathit{b}}[\mathit{k}]\right)^{\mathit{H}}$
$\displaystyle\mathbf{R}_{\mathit{v}_{\mathit{p}}}^{-1}\big{(}\mathbf{y}_{\mathit{p}}[\mathit{k}]-\tilde{\boldsymbol{\Psi}}_{\mathit{p}}[\mathit{k}]\mathbf{h}_{\mathit{b}}[\mathit{k}]\big{)}\big{)},$
(48)

where
$\vartheta[\mathit{k}]=\frac{1}{\left(\pi\right)^{\mathit{M}\mathit{N}_{\mathit{RF}}}|\mathbf{R}_{\mathit{v}_{\mathit{p}}}|}$.
Taking the logarithm of both sides of the equation above gives us

$\displaystyle\mathcal{L}\left(\mathbf{y}_{\mathit{p}}[\mathit{k}]|\mathbf{h}_{\mathit{b}}[\mathit{k}]\right)=-\mathit{M}\mathit{N}_{RF}\
\textrm{log}(\pi)-\textrm{log}\left[\textrm{det}(\mathbf{R}_{\mathit{v}_{\mathit{p}}})\right]-$
(49)
$\displaystyle\left(\mathbf{y}_{\mathit{p}}[\mathit{k}]-\tilde{\mathbf{\Psi}}_{\mathit{p}}[\mathit{k}]\mathbf{h}_{\mathit{b}}[\mathit{k}]\right)^{\mathit{H}}\mathbf{R}_{\mathit{v}_{\mathit{p}}}^{-1}\left(\mathbf{y}_{\mathit{p}}[\mathit{k}]-\tilde{\mathbf{\Psi}}_{\mathit{p}}[\mathit{k}]\mathbf{h}_{\mathit{b}}[\mathit{k}]\right).$

Let
$\mathbf{J}[\mathit{k}]\in\mathbb{C}^{\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}\times\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}}$
denote the Bayesian Fisher Information Matrix (FIM), formulated as

$\displaystyle\mathbf{J}[\mathit{k}]=\mathbf{J}_{\mathit{p}}[\mathit{k}]+\mathbf{J}_{\mathit{b}}[\mathit{k}],$
(50)

where
$\mathbf{J}_{\mathit{p}}[\mathit{k}]\in\mathbb{C}^{\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}\times\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}}$
denotes the FIM associated with the pilot output $\mathbf{y}_{p}[k]$, while
$\mathbf{J}_{\mathit{b}}[\mathit{k}]\in\mathbb{C}^{\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}\times\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}}$
represents the FIM corresponding to the a prior information of the beamspace
CSI $\mathbf{h}_{b}[k]$. The above quantities can be formulated as

$\displaystyle\mathbf{J}_{\mathit{p}}[\mathit{k}]=-\mathbb{E}_{\mathbf{y}_{\mathit{p}}[\mathit{k}],\mathbf{h}_{\mathit{b}}[\mathit{k}]}\left\\{\frac{\partial^{2}\mathcal{L}(\mathbf{y}_{\mathit{p}}[\mathit{k}]|\mathbf{h}_{\mathit{b}}[\mathit{k}])}{\partial\mathbf{h}_{\mathit{b}}[\mathit{k}]\partial\mathbf{h}_{\mathit{b}}^{\mathit{H}}[\mathit{k}]}\right\\},$
$\displaystyle\mathbf{J}_{\mathit{b}}[\mathit{k}]=-\mathbb{E}_{\mathbf{h}_{\mathit{b}}[\mathit{k}]}\left\\{\frac{\partial^{2}\mathcal{L}(\mathbf{h}_{\mathit{b}}[\mathit{k}];\mathbf{\Gamma}_{\mathit{k}})}{\partial\mathbf{h}_{\mathit{b}}[\mathit{k}]\partial\mathbf{h}_{\mathit{b}}^{\mathit{H}}[\mathit{k}]}\right\\}.$
(51)

Upon using the conditional PDF in (49), the FIM
$\mathbf{J}_{\mathit{p}}[\mathit{k}]$ can be derived as

$\displaystyle\mathbf{J}_{\mathit{p}}[\mathit{k}]=\tilde{\mathbf{\Psi}}_{\mathit{p}}^{\mathit{H}}[\mathit{k}]\mathbf{R}_{\mathit{v}_{\mathit{p}}}^{-1}\tilde{\mathbf{\Psi}}_{\mathit{p}}[\mathit{k}].$
(52)

Similarly, the FIM associated with $\mathbf{h}_{\mathit{b}}[\mathit{k}]$ is
obtained by using the a prior PDF from equation (46) in equation (51). The
corresponding expression for $\mathbf{J}_{\mathit{b}}[\mathit{k}]$ can be
formulated as

$\displaystyle\mathbf{J}_{\mathit{b}}[\mathit{k}]=\hat{\mathbf{\Gamma}}_{\mathit{k}}^{-1}.$
(53)

Using the expressions in (52) and (53), the Bayesian FIM is given by

$\displaystyle\mathbf{J}[\mathit{k}]=\tilde{\mathbf{\Psi}}_{\mathit{p}}^{\mathit{H}}[\mathit{k}]\mathbf{R}_{\mathit{v}_{\mathit{p}}}^{-1}\tilde{\mathbf{\Psi}}_{\mathit{p}}[\mathit{k}]+\hat{\mathbf{\Gamma}}_{\mathit{k}}^{-1}.$
(54)

Therefore, the BCRLB that pertains to the MSE of the estimate
$\hat{\mathbf{h}}_{\mathit{b}}[\mathit{k}]$ can be expressed as

$\displaystyle\begin{split}\textrm{MSE}(\hat{\mathbf{h}}_{\mathit{b}}[\mathit{k}])&=\mathbb{E}\left\\{\parallel\hat{\mathbf{h}}_{\mathit{b}}[\mathit{k}]-\mathbf{h}_{\mathit{b}}[\mathit{k}]\parallel^{2}\right\\}\geq\textrm{Tr}\left\\{\mathbf{J}^{-1}[\mathit{k}]\right\\}\\\
&=\textrm{Tr}\left\\{\left(\tilde{\mathbf{\Psi}}_{\mathit{p}}^{\mathit{H}}[\mathit{k}]\mathbf{R}_{\mathit{v}_{\mathit{p}}}^{-1}\tilde{\mathbf{\Psi}}_{\mathit{p}}[\mathit{k}]+\hat{\mathbf{\Gamma}}_{\mathit{k}}^{-1}\right)^{-1}\right\\}.\end{split}$
(55)

Upon exploiting the relationship between the beamspace representation of the
THz MIMO channel and its vectorized version, the BCRLB for the estimated CSI
is finally given by

$\displaystyle\textrm{MSE}(\hat{\mathbf{H}}[\mathit{k}])\geq\textrm{Tr}\left\\{\mathbf{\Psi}[\mathit{k}]\mathbf{J}^{-1}[\mathit{k}]\mathbf{\Psi}^{\mathit{H}}[\mathit{k}]\right\\}.$
(56)

A drawback of THz systems is the fact that the received signal is severely
affected by noise in the training phase. Therefore, it is difficult to obtain
an accurate CSI estimate using exclusively pilot symbols [38]. In order to
overcome this limitation, the DA-BL framework is proposed next which leverages
the data along with the pilot symbols for enhancing the estimation accuracy.

##  VI DA-BL For Sparse Channel Estimation in THz Hybrid MIMO systems

The data-model is developed next for our THz MIMO system, which facilitates
joint channel estimation and data-detection. Let
$\mathbf{x}_{\mathit{d},m}[\mathit{k}]$ denote the data symbol vector whose
elements are derived from an appropriate constellation. The data model
corresponding to the $\mathit{k}$-th subcarrier is given by

$\displaystyle\mathbf{y}_{\mathit{d},m}[\mathit{k}]=\mathbf{W}_{\mathit{RF},m}^{\mathit{H}}\mathbf{H}[\mathit{k}]\mathbf{F}_{\mathit{RF},m}\mathbf{x}_{\mathit{d},m}[\mathit{k}]+\mathbf{W}_{\mathit{RF},m}^{\mathit{H}}\mathbf{v}_{\mathit{d},m}[\mathit{k}],$
(57)

where we have $\mathbf{y}_{d,m}[k]\in\mathbb{C}^{N_{RF}\times
1}$,$\mathbf{W}_{\mathit{RF},m}\in\mathbb{C}^{\mathit{N}_{\mathit{R}}\times\mathit{N}_{\mathit{RF}}}$
and
$\mathbf{F}_{\mathit{RF},m}\in\mathbb{C}^{\mathit{N}_{\mathit{T}}\times\mathit{N}_{\mathit{RF}}}$.
The construction of $\mathbf{W}_{\mathit{RF},m}$ and
$\mathbf{F}_{\mathit{RF},m}$ is explained in Section VII-A. In (57)
$\mathbf{x}_{\mathit{d},m}[\mathit{k}]\in\mathbb{C}^{\mathit{N}_{\mathit{RF}}\times
1}$ and $\mathbf{v}_{\mathit{d},m}[\mathit{k}]\in\mathbb{C}^{N_{R}\times 1}$
represents the $\mathit{K}$-point FFT of $\mathbf{v}_{\mathit{d},m}$, where
$\mathbf{v}_{\mathit{d},m}\sim\mathcal{CN}(\mathbf{0}_{N_{R}\times
1},\sigma^{2}_{d}\mathbf{I}_{N_{R}})$. Let
$\tilde{\mathbf{v}}_{\mathit{d},m}[\mathit{k}]=\mathbf{W}_{\mathit{RF},m}^{\mathit{H}}\mathbf{v}_{\mathit{d},m}[\mathit{k}]$
represent the equivalent noise after receiver combining and
$\mathbf{H}_{\textrm{eq},m}[\mathit{k}]=\mathbf{W}_{\mathit{RF},m}^{\mathit{H}}\mathbf{H}[\mathit{k}]\mathbf{F}_{\mathit{RF},m}\in\mathbb{C}^{\mathit{N}_{\mathit{RF}}\times\mathit{N}_{\mathit{RF}}}$
represent the equivalent channel matrix. Using equation (57), one can readily
derive the expressions for the conventional ZF and MMSE based detectors, which
are respectively given by

$\displaystyle\hat{\mathbf{x}}_{\mathit{d},m}^{\textrm{ZF}}[\mathit{k}]=\left(\mathbf{H}_{\textrm{eq},m}^{\mathit{H}}[\mathit{k}]\mathbf{H}_{\textrm{eq},m}[\mathit{k}]\right)^{-1}\mathbf{H}_{\textrm{eq},m}^{\mathit{H}}[\mathit{k}]\mathbf{y}_{\mathit{d},m}[\mathit{k}],$
(58)
$\displaystyle\hat{\mathbf{x}}_{\mathit{d},m}^{\textrm{MMSE}}[\mathit{k}]=$
$\displaystyle\left(\mathbf{H}_{\textrm{eq},m}^{\mathit{H}}[\mathit{k}]\mathbf{H}_{\textrm{eq},m}[\mathit{k}]+\mathbf{R}_{\mathit{d},m}\right)^{-1}\mathbf{H}_{\textrm{eq},m}^{\mathit{H}}[\mathit{k}]\mathbf{y}_{\mathit{d},m}[\mathit{k}],$
(59)

where $\mathbf{R}_{\mathit{d},m}$ represents the noise covariance matrix for
$\mathit{k}$-th subcarrier, which is defined as

$\displaystyle\mathbf{R}_{d,m}=\mathbb{E}\left\\{\tilde{\mathbf{v}}_{\mathit{d},m}[\mathit{k}]\tilde{\mathbf{v}}_{\mathit{d},m}^{\mathit{H}}[\mathit{k}]\right\\}=\sigma^{2}_{d}\mathit{K}\mathbf{W}_{\mathit{RF},m}^{\mathit{H}}\mathbf{W}_{\mathit{RF},m}.$
(60)

The CSI estimate obtained in the previous section using the PA-BL algorithm is
used for obtaining the initial data symbol estimate of the MMSE estimator.
Note that during the estimation process, $\mathbf{W}^{H}_{RF,m}$ is determined
as per Section VII-A, while $\mathbf{W}^{H}_{BB,m}[k]$ is set as the identity
matrix. During the following data detection phase, $\mathbf{W}^{H}_{BB,m}[k]$
is designed as the LMMSE equalizer
$\left(\mathbf{H}_{\textrm{eq},m}^{\mathit{H}}[\mathit{k}]\mathbf{H}_{\textrm{eq},m}[\mathit{k}]+\mathbf{R}_{\mathit{d},m}\right)^{-1}\mathbf{H}_{\textrm{eq},m}^{\mathit{H}}[\mathit{k}]$.
By expliting the characterstics of the $\textrm{vec}(.)$ operator as discussed
in Section-I-C, one can recast the system model in equation (57) as

$\displaystyle\mathbf{y}_{\mathit{d},m}[\mathit{k}]=\mathbf{\Phi}_{\mathit{d},m}[\mathit{k}]\mathbf{h}[\mathit{k}]+\tilde{\mathbf{v}}_{\mathit{d},m}[\mathit{k}],$
(61)

where the equivalent sensing matrix
$\mathbf{\Phi}_{\mathit{d},m}[\mathit{k}]\in\mathbb{C}^{\mathit{N}_{\mathit{RF}}\times\mathit{N}_{\mathit{T}}\mathit{N}_{\mathit{R}}}$
is given by

$\displaystyle\mathbf{\Phi}_{\mathit{d},m}[\mathit{k}]=(\mathbf{x}_{\mathit{d},m}^{\mathit{T}}[\mathit{k}]\mathbf{F}_{\mathit{RF},m}^{\mathit{T}}\otimes\mathbf{W}_{\mathit{RF},m}^{\mathit{H}}).$
(62)

Finally, upon stacking the outputs $\mathbf{y}_{\mathit{d},m}[\mathit{k}]$
across all the data blocks, one obtains

$\displaystyle\underbrace{\begin{bmatrix}{\mathbf{y}}_{d,\text{1}}[k]\\\
\vdots\\\
{\mathbf{y}}_{d,N_{d}}[k]\end{bmatrix}}_{{\mathbf{y}_{\mathit{d}}}[k]}=\underbrace{\begin{bmatrix}\mathbf{\Phi}_{d,\text{1}}[\mathit{k}]\\\
\vdots\\\
\mathbf{\Phi}_{d,N_{d}}[\mathit{k}]\end{bmatrix}}_{\mathbf{\Phi_{\mathit{d}}[\mathit{k}]}}\mathbf{h}[\mathit{k}]+\underbrace{\begin{bmatrix}\tilde{\mathbf{v}}_{d,\text{1}}[k]\\\
\vdots\\\
\tilde{\mathbf{v}}_{d,N_{d}}[k]\end{bmatrix}}_{\mathbf{v}_{\mathit{d}}[k]},$
(63)

where $\mathbf{y}_{d}[k]\in\mathbb{C}^{N_{d}N_{RF}\times 1}$,
$\mathbf{\Phi}_{d}[k]\in\mathbb{C}^{N_{d}N_{RF}\times N_{T}N_{R}}$ and
$\mathbf{v}_{d}[k]\in\mathbb{C}^{N_{d}N_{RF}\times 1}$. Thus, the noise
covariance matrix obeys $\mathbf{R}_{v_{d}}\in\mathbb{C}^{N_{d}N_{RF}\times
N_{d}N_{RF}}$ where
$\mathbf{R}_{v_{d}}=\sigma^{2}_{d}K\>\text{blkdiag}\left(\left\\{\mathbf{W}_{RF,m}^{H}\mathbf{W}_{RF,m}\right\\}_{m=1}^{M}\right)$.
Eventually, to develop the joint CSI estimation model and data-detection, one
can further concatenate the outputs corresponding to the data and pilot models
i.e $\mathbf{y}_{d}[k]$ and $\mathbf{y}_{p}[k]$ respectively [18], from
equation (63) and (19) as

$\displaystyle\underbrace{\begin{bmatrix}{\mathbf{y}}_{\mathit{d}}[\mathit{k}]\\\
{\mathbf{y}}_{\mathit{p}}[\mathit{k}]\end{bmatrix}}_{{\mathbf{y}}[\mathit{k}]\
\in\
\mathbb{C}^{\left(\mathit{N}_{\mathit{d}}+\mathit{M}\right)\mathit{N}_{\mathit{RF}}\times
1}}=$
$\displaystyle\underbrace{\begin{bmatrix}\mathbf{\Phi}_{\mathit{d}}[\mathit{k}]\\\
\mathbf{\Phi}_{\mathit{p}}[\mathit{k}]\end{bmatrix}}_{\mathbf{\Phi[\mathit{k}]\
\in\
\mathbb{C}^{(\mathit{N}_{\mathit{d}}+\mathit{M})\mathit{N}_{\mathit{RF}}\times\mathit{N}_{\mathit{T}}\mathit{N}_{\mathit{R}}}}}\mathbf{h}[\mathit{k}]+$
$\displaystyle\underbrace{\begin{bmatrix}\mathbf{v}_{\mathit{d}}[\mathit{k}]\\\
\mathbf{v}_{\mathit{p}}[\mathit{k}]\end{bmatrix}}_{\mathbf{v}[k]\ \in\
\mathbb{C}^{(\mathit{N}_{\mathit{d}}+\mathit{M})\mathit{N}_{\mathit{RF}}\times
1}}.$ (64)

Furthermore, the equivalent sensing matrix for the data-aided CSI estimation
model is given by
$\tilde{\mathbf{\Psi}}[\mathit{k}]=\mathbf{\Phi}[\mathit{k}]\mathbf{\Psi}[\mathit{k}]\in\mathbb{C}^{\left(\mathit{N}_{\mathit{d}}+\mathit{M}\right)\mathit{N}_{\mathit{RF}}\times\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}}$.
The covariance matrix
$\mathbf{R}_{v}\in\mathbb{C}^{\left(\mathit{N}_{\mathit{d}}+\mathit{M}\right)\mathit{N}_{\mathit{RF}}\times\left(\mathit{N}_{\mathit{d}}+\mathit{M}\right)\mathit{N}_{\mathit{RF}}}$
of the noise vector $\mathbf{v}[\mathit{k}]$ is

$\displaystyle\mathbf{R}_{v}=\textrm{blkdiag}\left(\mathbf{R}_{{v}_{\mathit{d}}},\mathbf{R}_{{v}_{\mathit{p}}}\right).$
(65)

The proposed DA-BL scheme estimates the hyper-parameter matrix
$\mathbf{\Gamma}_{\mathit{k}}$ using the EM procedure described below. Let the
entire information set for this be defined as
$\boldsymbol{\eta}_{k}=\left\\{\mathbf{X}_{\mathit{d}}[\mathit{k}],\mathbf{\Gamma}_{\mathit{k}}\right\\}$,
where $\mathbf{X}_{d}[k]\in\mathbb{C}^{N_{RF}\times N_{d}}$ obtained by
horizontally concatenating the transmit data vectors across all blocks and the
parameter estimate at iteration $(j-1)$ be given as

$\displaystyle\hat{\boldsymbol{\eta}}_{k}^{\left(\mathit{j}-1\right)}=\left\\{\hat{\mathbf{X}}_{\mathit{d}}^{\left(\mathit{j}-1\right)}[\mathit{k}],\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}-1\right)}\right\\},$
(66)

where $\hat{\mathbf{X}}_{\mathit{d}}^{\left(\mathit{j}-1\right)}[\mathit{k}]$
denotes the estimate of the data in the $(\mathit{j}-1)$st iteration, whereas
$\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}-1\right)}$ denotes the
hyper-parameter matrix obtained in the $(\mathit{j}-1)$st iteration. The
conditional expectation of the log-likelihood function
$\mathcal{L}(\boldsymbol{\eta}_{k}|\hat{\boldsymbol{\eta}_{k}}^{\left(\mathit{j}-1\right)})$
calculated in the E-step for the $j$th iteration is given by

$\displaystyle\mathcal{L}(\boldsymbol{\eta}_{k}|\hat{\boldsymbol{\eta}_{k}}^{\left(\mathit{j}-1\right)})$
$\displaystyle=\mathbb{E}_{\mathbf{h}_{\mathit{b}}[\mathit{k}]|\mathbf{y}[\mathit{k}];\hat{\boldsymbol{\eta}_{k}}^{\left(\mathit{j}-1\right)}}\left\\{\textrm{log}\left[\mathit{p}\left(\mathbf{y}[\mathit{k}],\mathbf{h}_{\mathit{b}}[\mathit{k}];\boldsymbol{\eta}_{k}\right)\right]\right\\}$
$\displaystyle=\mathbb{E}\left\\{\textrm{log}\left[\mathit{p}\left(\mathbf{y}[\mathit{k}]|\mathbf{h}_{\mathit{b}}[\mathit{k}];\mathbf{X}_{\mathit{d}}[\mathit{k}]\right)\right]\right\\}+$
$\displaystyle\quad\quad\quad\quad\mathbb{E}\left\\{\textrm{log}\left[\mathit{p}\left(\mathbf{h}_{\mathit{b}}[\mathit{k}];\mathbf{\Gamma}_{\mathit{k}}\right)\right]\right\\}.$
(67)

By simplifying the first term in the above expectation, we obtain

$\displaystyle\mathcal{L}(\mathit{p}\left(\mathbf{y}[\mathit{k}]|\mathbf{h}_{\mathit{b}}[\mathit{k}];\mathbf{X}_{\mathit{d}}[\mathit{k}]\right))=-(\mathit{N}_{\mathit{d}}+\mathit{M})\mathit{N}_{\mathit{RF}}\textrm{log}(\pi)-$
$\displaystyle\textrm{log}\left[\textrm{det}(\mathbf{R}_{v})\right]-\left(\mathbf{y}[\mathit{k}]-\tilde{\mathbf{\Psi}}[\mathit{k}]\mathbf{h}_{\mathit{b}}[\mathit{k}]\right)^{\mathit{H}}\mathbf{R}_{v}^{-1}\times$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\left(\mathbf{y}[\mathit{k}]-\tilde{\mathbf{\Psi}}[\mathit{k}]\mathbf{h}_{\mathit{b}}[\mathit{k}]\right).$
(68)

Algorithm 2 BL-based channel estimation and robust data detection for THz MIMO
systems

Input: Observation vector $\mathbf{y}[\mathit{k}]$, Data output
$\mathbf{Y}_{\mathit{d}}[\mathit{k}]$, Pilot sensing matrix
$\mathbf{\Phi}_{\mathit{p}}[\mathit{k}]$, Equivalent sensing matrix
$\tilde{\mathbf{\Psi}}[\mathit{k}]$, Noise Covariance
$\mathbf{R}_{\mathit{v}}$, dictionary matrices
$\mathbf{A}_{\mathit{R}}(\Theta_{\mathit{R}},\mathit{f}_{\mathit{k}})$ and
$\mathbf{A}_{\mathit{T}}(\Theta_{\mathit{T}},\mathit{f}_{\mathit{k}})$,
stopping parameters $\epsilon$ and $\mathit{J}_{\textrm{max}}$
Initialization:
$\hat{\mathbf{X}}_{\mathit{d}}^{(\textrm{0})}[\mathit{k}]=\hat{\mathbf{X}}_{\mathit{d},\mathit{p}}^{\textrm{MMSE}}[\mathit{k}]$
estimated using the CSI obtained from PA-BL algorithm, counter $\mathit{j}=0$
while$(\parallel\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}\right)}-\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}-1\right)}\parallel_{\mathit{F}}>\epsilon)$
and $(\mathit{j}<\mathit{J}_{\textrm{max}})$ do

  * $\mathit{j}\leftarrow\mathit{j}+1$

  * E-step: Compute the aposteriori covariance and mean

$\boldsymbol{\Sigma}^{\left(\mathit{j}\right)}[\mathit{k}]=\left[(\tilde{\mathbf{\Psi}}^{(j-1)}[\mathit{k}])^{H}\mathbf{R}_{\mathit{v}}^{-\textrm{1}}\tilde{\mathbf{\Psi}}^{(j-1)}[\mathit{k}]+\left(\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}-1\right)}\right)^{-\textrm{1}}\right]^{-\textrm{1}}$
$\hat{\mathbf{h}}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}]=\boldsymbol{\Sigma}^{\left(\mathit{j}\right)}[\mathit{k}](\tilde{\mathbf{\Psi}}^{\left(\mathit{j}-1\right)}[\mathit{k}])^{\mathit{H}}\mathbf{R}_{\mathit{v}}^{-\textrm{1}}\mathbf{y}[\mathit{k}]$

  * M-step1: Update the hyperparameters

$\begin{split}\hat{\gamma}_{\mathit{k},\mathit{i}}^{\left(\mathit{j}\right)}=\boldsymbol{\Sigma}^{\left(\mathit{j}\right)}[\mathit{k}](\mathit{i},\mathit{i})+|\hat{\mathbf{h}}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}](\mathit{i})|^{2};\quad
1\leq\mathit{i}\leq\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}\end{split}$

  * M-step2:

    1. 1.

Evaluate $\overline{\mathbf{Y}}_{d}[k]$ and
$\hat{\boldsymbol{\mathcal{H}}}_{eq}^{(j)}[k]$ using equation (76)

    2. 2.

Update the data estimate $\hat{\mathbf{X}}_{d}^{\left(j\right)}[k]$ using (82)

    3. 3.

Demodulate the estimated data and use it to update
$\tilde{\mathbf{\Psi}}^{\left(\mathit{j}\right)}[\mathit{k}]$

end while
$\hat{\mathbf{h}}_{\mathit{b}}[\mathit{k}]=\hat{\mathbf{h}}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}]$
Output: $\hat{\mathbf{H}}_{\textrm{DA-
BL}}[\mathit{k}]=\mathbf{A}_{\mathit{R}}(\Theta_{\mathit{R}},\mathit{f}_{\mathit{k}})\textrm{vec}^{-1}(\hat{\mathbf{h}}_{\mathit{b}}[\mathit{k}])\mathbf{A}_{\mathit{T}}^{\mathit{H}}(\Theta_{\mathit{T}},\mathit{f}_{\mathit{k}})$,
$\hat{\mathbf{X}}_{\mathit{d}}[\mathit{k}]=\hat{\mathbf{X}}_{\mathit{d}}^{\left(\mathit{j}\right)}[\mathit{k}]$

Observe that the initial term relies only on the data matrix
$\mathbf{X}_{d}[k]$ through the matrix $\tilde{\mathbf{\Psi}}[\mathit{k}]$. In
a similar fashion, one can see that the second term is solely determined by
the hyper-parameter matrix $\boldsymbol{\Gamma}_{\mathit{k}}$

$\displaystyle\mathcal{L}(p\left[\mathbf{h}_{b}[k];\mathbf{\Gamma}_{k}\right])=-G_{R}G_{T}\mathrm{log}(\pi)-$
$\displaystyle\mathrm{log}[\mathrm{det}(\mathbf{\Gamma}_{k})]-$
$\displaystyle\mathbf{h}^{H}_{b}[k]\mathbf{\Gamma}^{-1}_{k}\mathbf{h}_{b}[k].$
(69)

Thus, the joint maximization of the log-likelihood
$\mathcal{L}(\boldsymbol{\eta}_{k}|\hat{\boldsymbol{\eta}_{k}}^{\left(\mathit{j}-1\right)})$
with respect to $\boldsymbol{\eta}_{k}$ in the M-step breaks down to two
independent maximization steps, one with respect to $\mathbf{X}_{d}[k]$ and
the other with $\boldsymbol{\Gamma}_{\mathit{k}}$. The optimization problem in
the M-step corresponding to the hyper-parameter update is given by

$\displaystyle\hat{\gamma}_{\mathit{k},\mathit{i}}^{\left(\mathit{j}\right)}=\underset{\gamma_{\mathit{k,i}}\geq
0}{\textrm{arg max}}\ $
$\displaystyle\mathbb{E}_{\mathbf{h}_{\mathit{b}}[\mathit{k}]|\mathbf{y}[\mathit{k}];\hat{\boldsymbol{\eta}_{k}}^{\left(\mathit{j}-1\right)}}\left\\{\textrm{log}\left[\mathit{p}\left(\mathbf{h}_{\mathit{b}}[\mathit{k}];\boldsymbol{\Gamma}_{\mathit{k}}\right)\right]\right\\}$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad\quad\forall\
1\leq\mathit{i}\leq\mathit{G}_{\mathit{R}}\mathit{G}_{\mathit{T}}.$ (70)

Along similar lines to the PA-BL technique, the estimate of
$\gamma_{\mathit{k},\mathit{i}}$ in iteration $j$ is given by

$\displaystyle\hat{\gamma}_{\mathit{k},\mathit{i}}^{\left(\mathit{j}\right)}$
$\displaystyle=\mathbb{E}_{\mathbf{h}_{\mathit{b}}[\mathit{k}]|\mathbf{y}[\mathit{k}];\hat{\boldsymbol{\Gamma}}_{k}^{\left(\mathit{j-\textrm{1}}\right)}}\left\\{|\mathbf{h}_{\mathit{b}}[\mathit{k}](\mathit{i})|^{2}\right\\}$
$\displaystyle=\boldsymbol{\Sigma}^{\left(\mathit{j}\right)}[\mathit{k}](\mathit{i},\mathit{i})+|\hat{\mathbf{h}}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}](\mathit{i})|^{2},$
(71)

where the a posteriori mean vector
$\hat{\mathbf{h}}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}]\in\mathbb{C}^{G_{R}G_{T}\times
1}$ is formulated as

$\displaystyle\hat{\mathbf{h}}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}]=\boldsymbol{\Sigma}^{\left(\mathit{j}\right)}[\mathit{k}]\left(\tilde{\boldsymbol{\Psi}}^{\left(\mathit{j-\textrm{1}}\right)}[\mathit{k}]\right)^{\mathit{H}}\mathbf{R}_{v}^{-1}\mathbf{y}[\mathit{k}],$
(72)

and its covariance matrix
$\boldsymbol{\Sigma}^{\left(\mathit{j}\right)}[\mathit{k}]\in\mathbb{C}^{G_{R}G_{T}\times
G_{R}G_{T}}$ is given by the expression

$\displaystyle\boldsymbol{\Sigma}^{\left(\mathit{j}\right)}[\mathit{k}]=\left[(\tilde{\mathbf{\Psi}}^{(j-1)}[\mathit{k}])^{H}\mathbf{R}_{\mathit{v}}^{-\textrm{1}}\tilde{\mathbf{\Psi}}^{(j-1)}[\mathit{k}]+\left(\hat{\mathbf{\Gamma}}_{\mathit{k}}^{\left(\mathit{j}-1\right)}\right)^{-\textrm{1}}\right]^{-\textrm{1}}.$
(73)

The matrix $\tilde{\mathbf{\Psi}}^{(j-\textrm{1})}[\mathit{k}]$ that appears
in the above expression depends on the data symbol matrix
$\hat{\mathbf{X}}_{d}^{(j-1)}[k]$ in the $\left(j-\textrm{1}\right)$st
iteration. Therefore, the M-step of the optimization procedure toward
determining the data matrix update $\hat{\mathbf{X}}_{d}^{\left(j\right)}[k]$
can be formulated as

$\displaystyle\hat{\mathbf{X}}_{d}^{(j)}[k]=\arg\max_{\mathbf{X}_{d}[k]}\mathbb{E}\bigg{\\{}\log\left[p\left(\mathbf{y}_{d}[k]\mid\mathbf{h}_{b}[k];\mathbf{X}_{d}[k]\right)\right]\bigg{\\}}.$
(74)

Using equation (57), the above maximization objective can be reformulated as
seen in equation (77). This can be further simplified as seen in (78)-(80).
After further simplification, (74) can be reframed as

$\displaystyle\hat{\mathbf{X}}_{d}^{\left(j\right)}[k]=$
$\displaystyle\arg\min_{\begin{subarray}{c}\mathbf{X}_{d}[k]\end{subarray}}\left\\{\left\|\overline{\mathbf{Y}}_{d}[k]-\widehat{\boldsymbol{\mathcal{H}}}_{eq}^{(j)}[k]\mathbf{X}_{d}[k]\right\|_{F}^{2}\right\\},$
(75)

where, we have

$\displaystyle\overline{\mathbf{Y}}_{d}[k]=\begin{bmatrix}\mathbf{Y}_{d}[k]\\\
\mathbf{0}\\\ \end{bmatrix},\ \
\hat{\boldsymbol{\mathcal{H}}}_{eq}^{(j)}[k]=\begin{bmatrix}\hat{\mathbf{H}}_{\textrm{eq}}^{(j)}[k]\\\
\left(\boldsymbol{\Xi}^{(j)}[k]\right)^{\frac{1}{2}}\\\ \end{bmatrix},$ (76)

and $\mathbf{Y}_{d}[k]\in\mathbb{C}^{N_{RF}\times N_{d}}$ is obtained by
horizontally concatenating the received data vectors across all blocks. The
quantity

$\displaystyle\hat{\mathbf{X}}_{d}^{\left(j\right)}[k]=$
$\displaystyle\arg\max_{\begin{subarray}{c}\mathbf{X}_{d}[k]\end{subarray}}\mathbb{E}_{\mathbf{H}_{\textrm{eq}}[k]\mid\mathbf{Y}_{d}[k];\hat{\mathbf{X}}_{d}^{\left(j-1\right)}[k]}\Bigg{\\{}\log\left[p\left(\mathbf{Y}_{d}[k]\mid\mathbf{H}_{\textrm{eq}}[k];\mathbf{X}_{d}[k]\right)\right]\Bigg{\\}}$
(77) $\displaystyle=$
$\displaystyle\arg\min_{\begin{subarray}{c}\mathbf{X}_{d}[k]\end{subarray}}\mathbb{E}_{\mathbf{H}_{\textrm{eq}}[k]\mid\mathbf{Y}_{d}[k];\hat{\mathbf{X}}_{d}^{\left(j-1\right)}[k]}\Bigg{\\{}\left\|\mathbf{Y}_{d}[k]-\mathbf{H}_{\textrm{eq}}[k]\mathbf{X}_{d}[k]\right\|_{F}^{2}\Bigg{\\}}$
(78) $\displaystyle=$
$\displaystyle\arg\min_{\begin{subarray}{c}\mathbf{X}_{d}[k]\end{subarray}}\left\\{\left\|\mathbf{Y}_{d}[k]-\hat{\mathbf{H}}_{\textrm{eq}}^{(j)}[k]\mathbf{X}_{d}[k]\right\|_{F}^{2}+\left\|\left(\boldsymbol{\Xi}^{(j)}[k]\right)^{\frac{1}{2}}\mathbf{X}_{d}[k]\right\|_{F}^{2}\right\\}$
(79) $\displaystyle=$
$\displaystyle\arg\min_{\begin{subarray}{c}\mathbf{X}_{d}[k]\end{subarray}}\left\\{\left\|\begin{bmatrix}\mathbf{Y}_{d}[k]\\\
\mathbf{0}_{N_{RF}\times N_{d}}\\\
\end{bmatrix}-\begin{bmatrix}\hat{\mathbf{H}}_{\textrm{eq}}^{(j)}[k]\\\
\left(\boldsymbol{\Xi}^{(j)}[k]\right)^{\frac{1}{2}}\\\
\end{bmatrix}\mathbf{X}_{d}[k]\right\|_{F}^{2}\right\\}$ (80)

$\hat{\mathbf{H}}_{\textrm{eq}}^{\left(\mathit{j}\right)}[\mathit{k}]\in\mathbb{C}^{\mathit{N}_{\mathit{RF}}\times\mathit{N}_{\mathit{RF}}}$
above is defined as

$\displaystyle\hat{\mathbf{H}}_{\textrm{eq}}^{\left(\mathit{j}\right)}[\mathit{k}]$
$\displaystyle=\mathbb{E}\left\\{\mathbf{W}_{RF}^{H}\mathbf{H}[k]\mathbf{F}_{RF}\right\\}$
$\displaystyle=\mathbf{W}_{\mathit{RF}}^{\mathit{H}}\mathbf{A}_{\mathit{R}}(\Theta_{\mathit{R}},\mathit{f}_{\mathit{k}})\textrm{vec}^{-1}\left(\hat{\mathbf{h}}_{\mathit{b}}^{\left(\mathit{j}\right)}[\mathit{k}]\right)\mathbf{A}_{\mathit{T}}^{\mathit{H}}(\Theta_{\mathit{T}},\mathit{f}_{\mathit{k}})\mathbf{F}_{\mathit{RF}},$
(81)

and the procedure of dividing $\boldsymbol{\Xi}^{(j)}[k]$ is detailed in
Appendix-A. The ZF based detector used for updating the data symbol matrix is
obtained as,

$\displaystyle\hat{\mathbf{X}}_{d}^{\mathrm{ZF}(j)}[k]=\left\\{\left(\hat{\boldsymbol{\mathcal{H}}}_{\mathrm{eq}}^{(j)}[k]\right)^{\mathit{H}}\hat{\boldsymbol{\mathcal{H}}}_{\mathrm{eq}}^{(j)}[k]\right\\}^{-1}\left(\hat{\boldsymbol{\mathcal{H}}}_{\mathrm{eq}}^{(j)}[k]\right)^{\mathit{H}}\overline{\mathbf{Y}}_{d}[k].$
(82)

The nearest neighbour detection rule is subsequently used to map each element
in the matrix above to a symbol in the transmit constellation. The E-step and
the M-step described above are executed iteratively until the algorithm
converges. The various steps of the DA-BL framework are summarized in
Algorithm $2$. The BCRLB for DA-BL framework is derived in Appendix B. In
general, for a sparse signal recovery problem, when considering an over-
complete dictionary, the performance of the proposed Bayesian learning schemes
such as PA-BL and DA-BL seem to remain consistent even with a reduced number
of observations. By contrast, the performance of non-Bayesian learning
techniques such as Basis-Pursuit, FOCUSS, OMP degrade significantly, as
observed in our simulation results.

###  VI-A Computational Complexity and Memory Occupancy

This subsection derives the computational complexity of the proposed PA-BL and
DA-BL techniques for THz hybrid MIMO systems. Due to space constraints, the
comprehensive calculations of the computational complexities for both the
proposed and existing THz hybrid MIMO channel estimation methods are provided
in our technical report [39]. The key findings from Tables I to III of the
technical report are discussed next. It is evident from Table I that the OMP
technique exhibits a worst-case complexity order of
$\mathcal{O}(M^{3}N_{RF}^{3})$, attributed to the necessity of an intermediate
least squares (LS) estimate at each iteration. As shown in Table II, the
computational complexity of the PA-BL technique is of the order
$\mathcal{O}(MN_{RF}^{3}+G_{R}^{3}G_{T}^{3})$, which arises due to the matrix
inversion of size $[G_{R}G_{T}\times G_{R}G_{T}]$ for computing the a
posteriori covariance $\mathbf{R}_{b}^{(j)}[k]$. On the other hand, from Table
III, one can see that the computational complexity of the proposed DA-BL
technique is of the order
$\mathcal{O}\big{(}(N_{d}+M)^{3}N_{RF}^{3}+G_{R}^{3}G_{T}^{3}\big{)}$, which
is higher than that of the PA-BL and OMP-based techniques. However, the
performance, as shown in the simulation results of Section VII, indicates
superior performance for DA-BL. The OMP technique has the lowest complexity,
but its performance is also poor. Furthermore, due to space constraints, we
have similarly derived the memory occupancy for different channel estimation
schemes, including OMP, PA-BL, and DA-BL in Table IV of our technical report
[39]. It is evident from the memory and computational complexity analysis,
that the proposed DA-BL technique has the highest memory requirement and
computational complexity, but it also exhibits superior performance. Hence,
one can see a clear trade-off performance vs. computational cost/ memory
between the proposed and existing CSI estimation schemes.

##  VII Simulation Results

This section provides the performance evaluation of the proposed PA-BL and DA-
BL techniques conceived for THz MIMO systems in the face of the dual-wideband
effect. The amplitudes of the path coefficients corresponding to the LoS and
NLoS elements are evaluated using equations (14) and (16), respectively, while
the corresponding phase shifts are produced using a probability distribution
that is uniform over the range (-$\pi$,$\pi$]. The coefficient of molecular
absorption is calculated using equation (15), and its parameters have been
taken from the HITRAN database [40]. The THz MIMO channel incorporating the
dual-wideband effect is generated using $N_{cl}=4$ clusters, namely one LoS
and $\mathit{N}_{\textrm{NLoS}}=3$ clusters. The carrier frequency utilized is
$f_{c}=1$ THz and the transmission distance considered is $d$ = 10 metres. We
consider an indoor office scenario, where the molecular composition of air is
given as water vapour $1\%$, oxygen $20.9\%$ and nitrogen $78.1\%$. The system
temperature and pressure, respectively, are considered to be 1 atm and 296K
for calculating the molecular absorption loss using equation (15). The AoAs/
AoDs associated with each multipath component have been generated from a
Laplacian distribution, where the multipath component’s angle exhibits a
standard deviation of $\frac{1}{10}$ radian from the mean, which in turn is
picked with uniform probability from an angular grid. Three reflecting media
have been considered with the standard deviations of their roughness set as
$\sigma\in{\left\\{0.05,0.13,0.15\right\\}}$ mm [31]. The transmit and receive
antenna gains are set to $A_{T}=31$ dB, $A_{R}=31$ dB, respectively. Two THz
MIMO systems are considered for simulations and their parameters are given in
Table II. The grid sizes are chosen so that $G_{T},G_{R}\geq$
max$(N_{T},N_{R})$. Both the systems have the number of pilot blocks $M$ set
as $M\in\left\\{30,50,70,90\right\\}$. The spacings at the transmit and
receive antenna ends are taken to be $d_{t}=d_{r}=\lambda/2$. Moreover, we
consider the utilization of a raised cosine pulse shaping filter with a roll-
off factor of 0.80 and an up-sampling factor of 20. The constituent elements
of the transmitted pilot symbol vector $\mathbf{u}_{m,n}^{(p)}$ are chosen at
unity average power from a 8-PSK constellation.

###  VII-A Wideband THz hybrid MIMO channel estimation

The RF TPC and RC components, which are frequency-independent as shown in
equation (2), are generated as [20]

$\displaystyle\mathbf{F}_{RF,m}\left(i,j\right)=\frac{1}{\sqrt{N_{t}}}\textrm{exp}\left(j\varphi_{i,j}\right),$
$\displaystyle\mathbf{W}_{RF,m}\left(i,j\right)=\frac{1}{\sqrt{N_{r}}}\textrm{exp}\left(j\psi_{i,j}\right).$
(83)

The phases $\varphi_{i,j}$ ,$\psi_{i,j}$ for the set $(i,j)$ are selected
randomly with a uniform probability distribution across the set
$\mathcal{A}=\left\\{\textrm{0},\frac{2\pi}{2^{N_{Q}}},\cdots,\frac{\left(2^{N_{Q}}-1\right)2\pi}{2^{N_{Q}}}\right\\}$
with $N_{Q}$ = 4 denoting the angle quantization parameter. Fig. 3(a) and (b)
compares the normalized MSE (NMSE) of the channel estimates for the proposed
PA-BL technique and the OMP technique for a dual-wideband THz MIMO channel,
with parameters as per for Systems I and II. The tuning parameter for the OMP
[35] algorithm which impacts the convergence of the algorithm is set as
$\epsilon_{o}=1$ for both the systems, while the stopping criterion and the
upper limit on the number of PA-BL iterations are defined as $\epsilon=1$ and
$K_{max}=20$, respectively. The NMSE metric is defined as
$\textrm{NMSE}=\frac{\sum_{\mathit{k}=0}^{\mathit{K}-1}\parallel\hat{\mathbf{H}}[k]-\mathbf{H}[k]\parallel_{F}^{2}}{\sum_{\mathit{k}=0}^{\mathit{K}-1}\parallel\mathbf{H}[k]\parallel_{F}^{\textrm{2}}}.$
Observe from Fig. 3(a) and (b), that the NMSE performance improves
significantly as one increases the number of training blocks $M$. The
explanation for this lies in the use of greater number of pilots which leads
to more reliable CSI estimates. It is evident from the plot that the proposed
PA-BL technique yields a considerably improved performance over the OMP
algorithm for the entire SNR range under consideration. The existing OMP
approach performs poorly owing to its susceptibility to changes in the
stopping criterion and dictionary matrix. As a consequence, even slight
modifications to the cessation criterion or deviations from the dictionary
matrix constraints can result in severe structural and convergence problems.

Table II: Simulation parameters for the dual-wideband THz system

Parameter | System-I | System-II
---|---|---
$\\#$ of TAs ($N_{T}$) | 32 | 64
$\\#$ of RAs ($N_{R}$) | 32 | 64
$\\#$ of RF chains ($N_{RF}$) | 6 | 12
$\\#$ of diffused ray ($N_{ray}$) | 1 | 1
$\\#$ of subcarriers ($K$) | 16 | 32
$\\#$ of data vectors ($N_{d}$) | 100 | 200
Transmit angular grid size | 64 | 128
Receive angular grid size | 64 | 128

(a)  (b)  (c)  (d)

Figure 3: Performance of the proposed PA-BL technique with the existing OMP
technique for $\left(a\right)$System-I $\left(b\right)$System-II.
$\left(c\right)$ NMSE performance of the proposed PA-BL algorithm versus SNR,
for various values of the bandwidth $B\in\left\\{20,50,100\right\\}$ GHz,
$\left(d\right)$ Normalized array gain versus number of antenna elements with
bandwidth $B=20$GHz at $f_{c}=1$ THz.

The OMP approach also suffers from error propagation, because errors in the
index selection cannot be corrected in subsequent iterations. It can be seen
that the NMSE of the proposed PA-BL scheme improves by approximately $\approx$
30 dB, when the SNR is varied from -10 dB to 20 dB. Upon increasing the number
of pilot blocks $M$, the performance of the PA-BL scheme improves
considerably. Interestingly, for a higher number of training frames such as
$M=90$, the proposed PA-BL algorithm performs close to the corresponding
BCRLB. This is of crucial importance, as the BCRLB is calculated under the
assumption of perfect knowledge of the AoAs/AoDs, which is an ideal case. In
stark contrast, the PA-BL technique does not require any such prior
information. Further, as demonstrated in [36], upon convergence of the PA-BL
algorithm, only the hyperparameters corresponding to the locations of the
active AoA/ AoD pairs retain significant values, while the rest of the
hyperparameters are driven to zero. By contrast, the performance of OMP and
many other sparse recovery algorithms depends highly on the cessation
criteria, which may be different for different propagation scenarios, thus
leading to performance degradation, if not set carefully. Hence the proposed
PA-BL algorithm is capable of learning the propagation environment in
different network deployments. Fig. 3(c), shows the performance of the PA-BL
technique for increasing signal bandwidth $B$. As discussed in Section III-C,
when the user bandwidth increases, the impact of the beam squint progressively
increases, which in turn degrades the performance of CSI estimation. It is
apparent from Fig. 3(c) that the increased NMSE of the channel estimate at
ultra-high bandwidths in the range $B=100$ GHz in a dual-wideband THz system
can be offset by augmenting the number of pilot blocks $M$. Moreover, the
proposed PA-BL technique is seen to yield robust performance over a wide
bandwidth range from $20$ GHz to $100$ GHz, even without any compensation for
the dual-wideband effect, which is convenient for practical implementation.
The proposed PA-BL framework is able to estimate the
$N_{R}\in\left\\{32,64\right\\}$ $\times$ $N_{T}\in\left\\{32,64\right\\}$
dual-wideband THz MIMO channel using $M\in\left\\{30,50,70,90\right\\}$ pilot
blocks and $N_{RF}\in\left\\{6,16\right\\}$, with corresponding $MN_{RF}$
values of $\left\\{180,300,420,540\right\\}$ for $N=6$ and
$\left\\{480,800,1120,1440\right\\}$ for $N=16$. Note that all these values
satisfy $MN_{RF}\ll N_{T}N_{R},$ where $N_{T}N_{R}=1024$ for System-I and
$N_{T}N_{R}=4096$ for System-II. This leads to an under-determined system in
Equation (30), wherein the conventional LS/MMSE estimators cannot be applied
for CSI estimation. In fact, the LS and MMSE require $M\approx 130$ pilot
blocks for System-I and $M\approx 250$ pilot blocks for System-II, which
increases the pilot overhead and reduces the spectral efficiency. This clearly
illustrates the unique advantage of the proposed novel scheme for THz CSI
estimation over the orthodox estimation paradigms. Fig. 3(d) depicts the
variation of the normalized array gain with respect to the number of antenna
elements $N_{R}$. As $N_{R}$ increases, the spatial wideband effect in the THz
system exhibits a significant impact due to the notable delay between the
first and last antenna elements, leading to a decrease in the array gain.
Consequently, it can be inferred that as $N_{R}\rightarrow\infty$,
$G(\theta,f_{k})\rightarrow 0$. Furthermore, Fig. 4(c) illustrates the
variation of the normalized array gain with respect to the bandwidth $B$. A
large bandwidth induces a significant shift in the relative subcarrier
frequency $\delta_{k}=\frac{f_{k}}{f_{c}}$, leading to a significant
degradation in normalized array gain. The results for the DA-BL framework that
can lead to a further improvement in the performance of the CSI estimation are
presented next.

(a)  (b)  (c)

Figure 4:  $\left(a\right)$ NMSE versus SNR comparison of the proposed PA-BL
and DA-BL techniques with the Genie-aided estimator, LS, OMP, FOCUSS and the
BCRLB, $\left(b\right)$ BER versus SNR comparison of the proposed PA-BL and
DA-BL techniques with the Genie-aided estimator for System-I. $\left(c\right)$
Normalized array gain versus bandwidth with $N_{T}=128$ and $N_{R}=128$ at
$f_{c}=1$THz.

###  VII-B DA-BL wideband THz hybrid MIMO channel estimation

Fig. 4 and 5 demonstrate the performance improvement of the proposed PA-BL and
DA-BL based techniques over traditional LS and other sparse estimation
techniques such as OMP [35] and FOCUSS [34] along with a hypothetical
scenario, where the true AoAs and AoDs are known, i.e, a Genie-aided
transceiver for THz hybrid MIMO systems with parameter settings corresponding
to System-I and System-II, as shown in Table II of the paper. Using the
concatenated model for all the $M$ blocks as discussed in the context of
equation (64), the hybrid TPC $\mathbf{F}_{RF}$ and hybrid RC
$\mathbf{W}_{RF}$ are chosen to be block-independent, and are set as described
in equation (83). As discussed, the performance of the traditional LS-based
scheme is inferior, because the technique does not leverage the sparsity
occurring in the THz MIMO channel. The performance of the OMP algorithm is
adversely affected by its sensitivity to the selection of the stopping
criteria and the choice of the dictionary matrix. Furthermore, FOCUSS suffers
from convergence-related problems and it is highly sensitive to the
regularization parameter. On the other hand, the PA-BL framework leverages
sparsity through EM-based hyperparameter estimation. Moreover, it eliminates
the need for tuning or regularization parameters and it is robust to the
choice of dictionary matrix. However, the exclusive reliance on pilot symbols
makes it challenging to attain a precise CSI estimate. Therefore, the proposed
DA-BL demonstrates a superior performance compared to other sparse estimation
techniques due to its ability to utilize the data symbols coupled with a
minimal pilot overhead. Furthermore the DA-BL framework beneficially leverages
the data symbols in addition to the pilot subcarriers to improve the channel
estimate by iteratively maximizing the joint likelihood, as discussed in
Section VI. The number of data vectors $N_{d}$ that contain the information
symbols for the DA-BL framework of System I and II is set as
$N_{d}=\left\\{100,200\right\\},$ respectively. The number of Monte Carlo
iterations for the DA-BL framework is considered to be $DL_{max}=20$ with the
threshold of $\epsilon=1$ as discussed in Algorithm 2. Additionally, it is
clear from the figure that the performance of the proposed DA-BL framework
approaches that of the Genie-aided system in the high SNR range, which
evidences the efficiency of the proposed scheme. To benchmark the performance,
the BCRLB is also plotted, as determined in Section V-A. The accuracy of data
detection at the receiver using the MMSE receiver for the proposed PA-BL, DA-
BL techniques and Genie-aided transceiver has been demonstrated by plotting
the BER in Fig. 4(b) and similarly, for System-II, in Fig. 5 (b). From 4(b)
one can further observe that there is performance improvement of approximately
$4$ dB is attained for the DA-BL framework over the PA-BL technique across the
SNR range of $-10$ dB to $15$ dB.

(a)  (b)  (c)

Figure 5:  $\left(a\right)$ NMSE versus SNR comparison of the proposed PA-BL
and DA-BL techniques with the Genie-aided estimator, LS, OMP, FOCUSS and the
BCRLB, $\left(b\right)$ BER versus SNR comparison of the proposed PA-BL and
DA-BL techniques with the Genie-aided estimator for System-II.
$\left(c\right)$ NMSE of the DA-BL scheme with increasing number of iterations
of DA-BL framework.

For both the settings, the DA-BL-based receiver can be seen to yield a
markedly improved performance in comparison to its PA-BL counterpart.
Furthermore, the performance of the former, when System-II is considered, is
also extremely close to that of a Genie-aided solution, and almost overlaps in
the high-SNR regime, as it is clearly seen in Fig. 5(a). Fig. 5(c) illustrates
the impact of the number of data-aided iterations on NMSE performance for SNR
set as $10$ dB. Importantly, the performance of the PA-BL framework and Genie-
aided estimator are insensitive to the number of data decoding iterations and
are thus constant in this figure. One can also observe that the NMSE of the
proposed DA-BL scheme achieves that of the Genie-aided receiver for as few as
$8$ DA-BL iterations. This is of considerable importance, as the DA-BL
framework achieves this without the need for prior knowledge of the true AoA/
AoD. This underscores its suitability for practical implementation in THz
systems, where prior information is typically unavailable. Another noteworthy
aspect is that, similar to the PA-BL receiver considered previously, for the
setting with $M=30$ pilot blocks and $N_{d}=100$ data vectors, the relation
$(M+N_{d})N_{RF}=780\ll N_{T}N_{R}=1024$ once again holds true for the DA-BL
receiver in System-I, which results in an under-determined estimation model as
per Equation (30). Thus, this clearly shows that even in such a scenario,
wherein the conventional LS and MMSE estimators are ineffective, the dual-
wideband channel can be estimated precisely, lower complexity and a notably
smaller training overhead, employing the proposed Data-Aided framework.
Moreover, the Bayesian learning algorithm has demonstrated superior
performance compared to other sparse recovery techniques in scenarios where
the spatial wideband effect is negligible [6], [18]. However, when the spatial
wideband effect is significant, both the CSI estimation accuracy and the
resultant BER performance degrade. This deterioration is attributed to the
subcarrier-dependent nature of the transmit and receive dictionary matrices
i.e $\mathbf{A}_{T}\left(\Theta_{T},f_{k}\right)$ and
$\mathbf{A}_{R}\left(\Theta_{R},f_{k}\right)$ respectively. To address this
issue, we seek a sparse recovery technique, which is robust to the choice of
dictionary matrix and ensures convergence. In this context, the PA-BL
algorithm emerges as an ideal choice. Additionally, to enhance the estimation
accuracy and mitigate the impact of the spatial wideband effect, this work
also proposes a data-aided Bayesian learning technique. This technique
iteratively incorporates data symbol estimates along with the development of a
robust ZF based data symbol detector. Notably, other sparse signal recovery
techniques such as OMP, FOCUSS and LASSO perform poorly in the presence of the
spatial-wideband effect. Furthermore, there is no systematic approach to
formulate their data-aided extensions. Thus, the relative performance
improvement of the Bayesian learning based schemes over OMP/ FOCUSS/ LASSO is
significant in the presence of the spatial-wideband effect. In short, the
excellent CSI estimation performance, coupled with lower overheads and robust
data-detection, render the proposed frameworks eminently suitable for
application in dual-wideband THz MIMO systems.

##  VIII Conclusions

New channel estimation paradigms considering angular domain sparsity in
practical THz MIMO wideband channels have been conceived. The THz channel was
modeled rigorously using parameters derived from the HITRAN database,
considering the free-space losses in addition to the dual-wideband effect
naturally occurring in such systems. An SC-FDE based I/O model was developed
to beneficially transform the TD model to its equivalent FD representation.
Subsequently, exploiting the sparsity of the dual-wideband THz channel, novel
Bayesian learning approaches were proposed for CSI recovery. To begin with,
the PA-BL methodology was conceived for channel estimation, that uses only
pilot symbols, followed by the DA-BL approach for superior performance, that
leverages the concatenated pilot and data outputs. By explicitly deriving the
BCRLBs, we characterize the estimation performance of PA-BL and DA-BL
frameworks. A comprehensive set of simulation results was presented that
evidences the NMSE and BER performance benefits of the various schemes
formulated in this paper. In particular, the NMSE of the DA-BL paradigm was
seen to closely follow that of the Genie-aided transceiver and the BCRLB
references at high SNRs, a feature that shows the advantages of practical
deployment of the different algorithms devised.

##  Appendix A  calculation of the matrix $\boldsymbol{\Xi}^{(j)}[k]$

The optimization cost function in (78) can be simplified as

$\displaystyle\mathbb{E}\left\\{\left\|\mathbf{Y}_{d}[k]-\mathbf{H}_{\textrm{eq}}[k]\mathbf{X}_{d}[k]\right\|_{F}^{2}\right\\}=$
$\displaystyle\mathrm{Tr}\bigg{\\{}\mathbf{Y}_{d}^{H}[k]\mathbf{Y}_{d}[k]-\mathbf{Y}_{d}^{H}[k]\hat{\mathbf{H}}_{\textrm{eq}}^{(j)}[k]\mathbf{X}_{d}[k]-$
$\displaystyle\mathbf{X}_{d}^{H}[k]\left(\hat{\mathbf{H}}_{\textrm{eq}}^{(j)}[k]\right)^{H}\mathbf{Y}_{d}[k]-\mathbf{X}_{d}^{H}[k]\mathbb{E}\left\\{\mathbf{H}_{\textrm{eq}}^{H}[k]\mathbf{H}_{\textrm{eq}}[k]\right\\}\mathbf{X}_{d}[k]\bigg{\\}},$
(84)

where $\hat{\mathbf{H}}_{\textrm{eq}}^{(j)}[k]$ is given by (81). Also the
quantity
$\mathbb{E}\left\\{\mathbf{H}_{\textrm{eq}}^{H}[k]\mathbf{H}_{\textrm{eq}}[k]\right\\}$
is expressed as

$\displaystyle\mathbb{E}\left\\{\mathbf{H}_{\textrm{eq}}^{H}[k]\mathbf{H}_{\textrm{eq}}[k]\right\\}=\left(\hat{\mathbf{H}}_{\textrm{eq}}^{(j)}[k]\right)^{H}\hat{\mathbf{H}}_{\textrm{eq}}^{(j)}[k]+\boldsymbol{\Xi}^{(j)}[k].$
(85)

The procedure for computation of the matrix $\boldsymbol{\Xi}^{(j)}[k]$ is
described next. Let
$\mathbf{h}_{\textrm{eq}}[k]\in\mathbb{C}^{N_{RF}^{2}\times 1}$ be the
vectorized form of the equivalent channel matrix given by
$\mathbf{h}_{\textrm{eq}}[k]=\mathrm{vec}\left(\mathbf{H}_{\textrm{eq}}[k]\right)$
and its estimate $\hat{\mathbf{h}}_{\textrm{eq}}^{(j)}[k]$ be defined as
$\hat{\mathbf{h}}_{\textrm{eq}}^{(j)}[k]=\mathrm{vec}\left(\hat{\mathbf{H}}_{\textrm{eq}}^{(j)}[k]\right)$.
Its error covariance matrix
$\boldsymbol{\Sigma}_{\textrm{eq}}^{(j)}[k]\in\mathbb{C}^{N_{RF}^{2}\times
N_{RF}^{2}}$ can be derived by substituting the value of
$\hat{\mathbf{H}}_{\textrm{eq}}^{(j)}[k]$ from (81) given

$\displaystyle\boldsymbol{\Sigma}_{\textrm{eq}}^{(j)}[k]=\widetilde{\boldsymbol{\Phi}}[k]\boldsymbol{\Sigma}^{(j)}[k]\widetilde{\boldsymbol{\Phi}}^{H}[k],$
(86)

where
$\widetilde{\boldsymbol{\Phi}}[k]=\left[\mathbf{F}_{{RF}}^{T}\otimes\mathbf{W}_{{RF}}^{H}\right]\boldsymbol{\Psi}[k]\in\mathbb{C}^{N_{RF}^{2}\times
G_{R}G_{T}}$. Subsequently, the $(a,b)$th element of the matrix
$\boldsymbol{\Xi}^{(j)}[k]$ is defined

$\displaystyle\boldsymbol{\Xi}^{(j)}[k](a,b)=$
$\displaystyle\mathrm{Tr}\big{[}\boldsymbol{\Sigma}_{\textrm{eq}}^{(j)}[k](\tilde{a}-N_{RF}+1:\tilde{a},$
$\displaystyle\tilde{b}-N_{RF}+1:\tilde{b})\big{]},$ (87)

where $\tilde{a}=aN_{RF}$ and $\tilde{b}=bN_{RF}$.

##  Appendix B bayesian cramér-rao lower bound for data-aided learning

Let $\mathbf{J}_{\mathrm{DA}}[k]\in\mathbb{C}^{G_{R}G_{T}\times G_{R}G_{T}}$
denote the combined Bayesian FIM, which can be evaluated as

$\displaystyle\mathbf{J}_{\mathrm{DA}}[k]=\mathbf{J}_{d,\mathrm{DA}}[k]+\mathbf{J}_{b,\mathrm{DA}}[k],$
(88)

where $\mathbf{J}_{d,\mathrm{DA}}[k]\in\mathbb{C}^{G_{R}G_{T}\times
G_{R}G_{T}}$ denotes the FIM associated with the data output, while
$\mathbf{J}_{b,\mathrm{DA}}[k]$ represents the FIM corresponding to the prior
information of the beamspace CSI $\mathbf{h}_{b}[k]$. The above quantities can
be computed as

$\displaystyle\mathbf{J}_{d,\mathrm{DA}}[k]=-\mathbb{E}_{\mathbf{y}[k],\mathbf{h}_{b}[k]}\left\\{\frac{\partial^{2}\mathcal{L}(\mathbf{y}[k]|\mathbf{h}_{b}[k];\mathbf{X}_{d}[k])}{\partial\mathbf{h}_{b}[k]\partial\mathbf{h}^{H}_{b}[k]}\right\\}$
$\displaystyle\mathbf{J}_{b,\mathrm{DA}}[k]=-\mathbb{E}_{\mathbf{h}_{b}[k]}\left\\{\frac{\partial^{2}\mathcal{L}(\mathbf{h}_{b}[k];\mathbf{\Gamma}_{k})}{\partial\mathbf{h}_{b}[k]\partial\mathbf{h}^{H}_{b}[k]}\right\\}.$
(89)

Using the conditional PDF as obtained in (67), the FIM
$\mathbf{J}_{d,\mathrm{DA}}[k]$ can be obtained as

$\displaystyle\mathbf{J}_{d,\mathrm{DA}}[k]=\tilde{\mathbf{\Psi}}[k]\mathbf{R}^{-1}_{v}\tilde{\mathbf{\Psi}}[k].$
(90)

Similarly, as obtained in the PA-BL, the FIM associated with
$\mathbf{h}_{b}[k]$ is obtained using the prior PDF from equation (68), and
the corresponding expression for $\mathbf{J}_{b}[k]$ can be obtained as

$\displaystyle\mathbf{J}_{b,\mathrm{DA}}[k]=\hat{\mathbf{\Gamma}}_{k}^{-1}.$
(91)

Using the expressions in (90) and (91), the Bayesian FIM can be represented as

$\displaystyle\mathbf{J}_{\mathrm{DA}}[\mathit{k}]=\tilde{\mathbf{\Psi}}^{\mathit{H}}[\mathit{k}]\mathbf{R}_{\mathit{v}}^{-1}\tilde{\mathbf{\Psi}}[\mathit{k}]+\hat{\mathbf{\Gamma}}_{\mathit{k}}^{-1}.$
(92)

Thus, the lower bound on the MSE of the estimate
$\hat{\mathbf{h}}_{\mathit{b}}[\mathit{k}]$, can be obtained as

$\displaystyle\begin{split}\textrm{MSE}(\hat{\mathbf{h}}_{\mathit{b}}[\mathit{k}])&=\mathbb{E}\left\\{\parallel\hat{\mathbf{h}}_{\mathit{b}}[\mathit{k}]-\mathbf{h}_{\mathit{b}}[\mathit{k}]\parallel^{2}\right\\}\geq\textrm{Tr}\left\\{\mathbf{J}^{-1}_{\mathrm{DA}}[\mathit{k}]\right\\}\\\
&=\textrm{Tr}\left\\{\left(\tilde{\mathbf{\Psi}}^{\mathit{H}}[\mathit{k}]\mathbf{R}_{\mathit{v}}^{-1}\tilde{\mathbf{\Psi}}[\mathit{k}]+\hat{\mathbf{\Gamma}}_{\mathit{k}}^{-1}\right)^{-1}\right\\}.\end{split}$
(93)

Exploiting the relationship between the beamspace representation of the THz
MIMO channel and its vectorized version, the BCRLB for the estimated CSI is
finally given by

$\displaystyle\textrm{MSE}(\hat{\mathbf{H}}[\mathit{k}])\geq\textrm{Tr}\left\\{\mathbf{\Psi}[\mathit{k}]\mathbf{J}^{-1}_{\mathrm{DA}}[\mathit{k}]\mathbf{\Psi}^{\mathit{H}}[\mathit{k}]\right\\}.$
(94)

## References

  * [1] I. F. Akyildiz, C. Han, Z. Hu, S. Nie, and J. M. Jornet, “Terahertz band communication: An old problem revisited and research directions for the next decade,” _IEEE Transactions on Communications_ , vol. 70, no. 6, pp. 4250–4285, 2022.
  * [2] A. Shafie, N. Yang, C. Han, J. M. Jornet, M. Juntti, and T. Kurner, “Terahertz communications for 6G and beyond wireless networks: Challenges, key advancements, and opportunities,” _IEEE Network (Early Access)_ , pp. 1–8, 2022.
  * [3] H. Sarieddeen, N. Saeed, T. Y. Al-Naffouri, and M.-S. Alouini, “Next generation Terahertz communications: A rendezvous of sensing, imaging, and localization,” _IEEE Communications Magazine_ , vol. 58, no. 5, pp. 69–75, 2020.
  * [4] X. Yang, M. Li, Y. Wang, and Y. Zeng, “Atmospheric absorption loss estimation of Terahertz-wave band satellite ground detection,” in _Proceedings of the 18th IEEE International Conference on Communication Technology_ , 2018, pp. 741–744.
  * [5] H. Sarieddeen, M.-S. Alouini, and T. Y. Al-Naffouri, “An overview of signal processing techniques for Terahertz communications,” _Proceedings of the IEEE_ , vol. 109, no. 10, pp. 1628–1665, 2021.
  * [6] S. Srivastava, A. Tripathi, N. Varshney, A. K. Jagannatham, and L. Hanzo, “Hybrid transceiver design for Tera-hertz MIMO systems relying on Bayesian learning aided sparse channel estimation,” _IEEE Transactions on Wireless Communications_ , vol. 22, no. 4, pp. 2231–2245, 2023\.
  * [7] A. Faisal, H. Sarieddeen, H. Dahrouj, T. Y. Al-Naffouri, and M.-S. Alouini, “Ultra-massive MIMO systems at Terahertz bands: Prospects and challenges,” _IEEE Vehicular Technology Magazine_ , vol. 15, no. 4, pp. 33–42, 2020.
  * [8] Y.-P. Lin, “Hybrid MIMO-OFDM beamforming for wideband mmWave channels without instantaneous feedback,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 19, pp. 5142–5151, 2018.
  * [9] H. Li, M. Li, Q. Liu, and A. L. Swindlehurst, “Dynamic hybrid beamforming with low-resolution PSs for wideband mmWave MIMO-OFDM systems,” _IEEE Journal on Selected Areas in Communications_ , no. 9, pp. 2168–2181, 2020.
  * [10] B. Ning, Z. Tian, Z. Chen, C. Han, J. Yuan, and S. Li, “Prospective beamforming technologies for ultra-massive MIMO in Terahertz communications: A tutorial,” _arXiv preprint arXiv:2107.03032_ , 2021.
  * [11] H. Elayan, O. Amin, R. M. Shubair, and M.-S. Alouini, “Terahertz communication: The opportunities of wireless technology beyond 5G,” in _2018 International Conference on Advanced Communication Technologies and Networking (CommNet)_. IEEE, 2018, pp. 1–5.
  * [12] L. Dai, J. Tan, Z. Chen, and H. V. Poor, “Delay-phase precoding for wideband THz massive MIMO,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 9, pp. 7271–7286, 2022.
  * [13] M. Jian, F. Gao, Z. Tian, S. Jin, and S. Ma, “Angle-domain aided UL/DL channel estimation for wideband mmwave massive MIMO systems with beam squint,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 7, pp. 3515–3527, 2019.
  * [14] J. M. Jornet and I. F. Akyildiz, “Channel modeling and capacity analysis for electromagnetic wireless nanonetworks in the Terahertz band,” _IEEE Transactions on Wireless Communications_ , vol. 10, no. 10, pp. 3211–3221, 2011\.
  * [15] A. Saeed, O. Gurbuz, and M. A. Akkas, “Terahertz communications at various atmospheric altitudes,” _Physical Communication_ , vol. 41, p. 101113, 2020\.
  * [16] C. Lin and G. Y. L. Li, “Terahertz communications: An array-of-subarrays solution,” _IEEE Communications Magazine_ , vol. 54, no. 12, pp. 124–131, 2016.
  * [17] Y. Zhang, Y. Huo, D. Wang, X. Dong, and X. You, “Channel estimation and hybrid precoding for distributed phased arrays based MIMO wireless communications,” _IEEE Transactions on Vehicular Technology_ , no. 11, pp. 12 921–12 937, 2020.
  * [18] S. Srivastava, J. Nath, and A. K. Jagannatham, “Data aided quasistatic and doubly-selective CSI estimation using affine-precoded superimposed pilots in millimeter wave MIMO-OFDM systems,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 7, pp. 6983–6998, 2021.
  * [19] K. Venugopal, A. Alkhateeb, N. G. Prelcic, and R. W. Heath, “Channel estimation for hybrid architecture-based wideband millimeter wave systems,” _IEEE Journal on Selected Areas in Communications_ , vol. 35, no. 9, pp. 1996–2009, 2017.
  * [20] S. Srivastava, P. Sharma, S. Dwivedi, A. K. Jagannatham, and L. Hanzo, “Fast block LMS based estimation of angularly sparse channels for single-carrier wideband millimeter wave hybrid MIMO systems,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 1, pp. 666–681, 2021.
  * [21] X. Gao, L. Dai, S. Zhou, A. M. Sayeed, and L. Hanzo, “Wideband beamspace channel estimation for millimeter-wave MIMO systems relying on lens antenna arrays,” _IEEE Transactions on Signal Processing_ , no. 18, pp. 4809–4824, 2019.
  * [22] B. Wang, F. Gao, S. Jin, H. Lin, and G. Y. Li, “Spatial-and frequency-wideband effects in millimeter-wave massive MIMO systems,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 13, pp. 3393–3406, 2018.
  * [23] B. Wang, M. Jian, F. Gao, G. Y. Li, and H. Lin, “Beam squint and channel estimation for wideband mmWave massive MIMO-OFDM systems,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 23, pp. 5893–5908, 2019.
  * [24] K. Dovelos, M. Matthaiou, H. Q. Ngo, and B. Bellalta, “Channel estimation and hybrid combining for wideband Terahertz massive MIMO systems,” _IEEE Journal on Selected Areas in Communications_ , vol. 39, no. 6, pp. 1604–1620, 2021.
  * [25] T.-H. Chou, N. Michelusi, D. J. Love, and J. V. Krogmeier, “Dual-Wideband Time-Varying Sub-Terahertz Massive MIMO Systems: A Compressed Training Framework,” _arXiv preprint arXiv:2201.00992_ , 2022.
  * [26] J. Rodríguez-Fernández, N. González-Prelcic, K. Venugopal, and R. W. Heath, “Frequency-domain compressive channel estimation for frequency-selective hybrid millimeter wave MIMO systems,” _IEEE Transactions on Wireless Communications_ , vol. 17, no. 5, pp. 2946–2960, 2018\.
  * [27] Z. Sha and Z. Wang, “Channel estimation and equalization for Terahertz receiver with RF impairments,” _IEEE Journal on Selected Areas in Communications_ , vol. 39, no. 6, pp. 1621–1635, 2021.
  * [28] A. M. Elbir, W. Shi, K. V. Mishra, and S. Chatzinotas, “Federated Learning for THz Channel Estimation,” _arXiv preprint arXiv:2207.06017_ , 2022\.
  * [29] M. Wang, F. Gao, N. Shlezinger, M. F. Flanagan, and Y. C. Eldar, “A block sparsity based estimator for mmWave massive MIMO channels with beam squint,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 49–64, 2019\.
  * [30] S. Ju, Y. Xing, O. Kanhere, and T. S. Rappaport, “Millimeter wave and sub-terahertz spatial statistical channel model for an indoor office building,” _IEEE Journal on Selected Areas in Communications_ , vol. 39, no. 6, pp. 1561–1575, 2021.
  * [31] R. Piesiewicz, C. Jansen, D. Mittleman, T. Kleine-Ostmann, M. Koch, and T. Kurner, “Scattering analysis for the modeling of THz communication systems,” _IEEE Transactions on Antennas and Propagation_ , vol. 55, no. 11, pp. 3002–3009, 2007.
  * [32] J. Rodriguez-Fernandez and N. Gonzalez-Prelcic, “Channel estimation for frequency-selective mmwave MIMO systems with beam-squint,” in _Proceedings of the IEEE Global Communications Conference (GLOBECOM)_ , 2018, pp. 1–6.
  * [33] C. Ekanadham, D. Tranchina, and E. P. Simoncelli, “Recovery of sparse translation-invariant signals with continuous basis pursuit,” _IEEE transactions on signal processing_ , vol. 59, no. 10, pp. 4735–4744, 2011.
  * [34] S. F. Cotter, B. D. Rao, K. Engan, and K. Kreutz-Delgado, “Sparse solutions to linear inverse problems with multiple measurement vectors,” _IEEE Transactions on Signal Processing_ , vol. 53, no. 7, pp. 2477–2488, 2005.
  * [35] A. Kulkarni and T. Mohsenin, “Low overhead architectures for OMP compressive sensing reconstruction algorithm,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 64, no. 6, pp. 1468–1480, 2017.
  * [36] D. P. Wipf and B. D. Rao, “Sparse Bayesian learning for basis selection,” _IEEE Transactions on Signal processing_ , vol. 52, no. 8, pp. 2153–2164, 2004.
  * [37] A. Mishra, A. Rajoriya, A. K. Jagannatham, and G. Ascheid, “Sparse Bayesian learning-based channel estimation in millimeter wave hybrid MIMO systems,” in _Proceedings of the 18th IEEE International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_ , 2017, pp. 1–5.
  * [38] S. Alzeer, S. Almatrudi, H. A. Bukhari, Y. Han, Y. Ding, and B. D. Rao, “Millimeter wave channel estimation using data-aided DoA estimation,” in _Proceedings of the 52nd IEEE Asilomar Conference on Signals, Systems, and Computers_ , 2018, pp. 1567–1571.
  * [39] A. Garg, S. Srivastava, N. Yadav, A. K. Jagannatham, and L. Hanzo. “Technical Report: Angularly Sparse Channel Estimation in Dual-Wideband Tera-Hertz (THz) Hybrid MIMO Systems Relying on Bayesian Learning" IIT Kanpur, Tech. Rep., 2023. [Online]. Available: https://www.iitk.ac.in/mwn/documents/MWNLab_TR_THz2023.pdf
  * [40] C. Hill. (1999) HITRAN Online Interface to the HITRAN Database. [Online]. Available: https://hitran.org/

|  Abhisha Garg (Graduate Student Member, IEEE) received the B.Tech degree in
Electronics & Communication from Graphic Era University, Dehradun, India, in
$2016$, and the M.Tech degree in Signal Processing from Guru Gobind Singh