# Hybrid mmWave MIMO Systems under Hardware Impairments and Beam Squint:
Channel Model and Dictionary Learning-aided Configuration

Hongxiang Xie, Joan Palacios and Nuria González-Prelcic∗ ∗Hongxiang Xie was
with the Department of Electrical and Computer Engineering at the University
of Texas at Austin when this work was performed. Nuria González-Prelcic and
Joan Palacios are with the Electrical and Computer Engineering Department,
North Carolina State University, USA. The corresponding author is Nuria
González-Prelcic<EMAIL_ADDRESS>
###### Abstract

Low overhead channel estimation based on compressive sensing (CS) has been
widely investigated for hybrid wideband millimeter wave (mmWave) multiple-
input multiple-output (MIMO) systems. The channel sparsifying dictionaries
used in prior work are built from ideal array response vectors evaluated on
discrete angles of arrival/departure. In addition, these dictionaries are
assumed to be the same for all subcarriers, without considering the impacts of
hardware impairments and beam squint. In this manuscript, we derive a general
channel and signal model that explicitly incorporates the impacts of hardware
impairments, practical pulse shaping functions, and beam squint, overcoming
the limitations of mmWave MIMO channel and signal models commonly used in
previous work. Then, we propose a dictionary learning (DL) algorithm to obtain
the sparsifying dictionaries embedding hardware impairments, by considering
the effect of beam squint without introducing it into the learning process. We
also design a novel CS channel estimation algorithm under beam squint and
hardware impairments, where the channel structures at different subcarriers
are exploited to enable channel parameter estimation with low complexity and
high accuracy. Numerical results demonstrate the effectiveness of the proposed
DL and channel estimation strategy when applied to realistic mmWave channels.

###### Index Terms:

Dictionary learning, compressive sensing, millimeter wave (mmWave), massive
MIMO, beam squint, spatial wideband effect, hardware impairments, mutual
coupling, antenna spacing error, channel estimation, sparse coding, dictionary
update.

##  I Introduction

The acquisition of channel state information (CSI) is crucial for mmWave link
configuration, and challenging when operating with hybrid beamforming
architectures. To reduce the training overhead associated with CSI
acquisition, prior work has made full use of the sparse nature of mmWave
channels in the angular or delay domains [1, 2, 3, 4, 5]. Nevertheless, some
relevant practical aspects have not been fully considered in previous
compressive channel models and estimation algorithms: the beam squint effect,
calibration errors, and hardware impairments. Specifically, the channel
sparsifying dictionaries used in prior work are typically assumed to be
(overcomplete) discrete Fourier transform (DFT) matrices, or constructed from
the ideal array response matrices (IARM) evaluated on discrete grids of
quantized angles of arrivals and departures (AoAs/AoDs) [2, 4]. These
assumptions are valid, however, only when the beam squint effect is negligible
and no hardware impairments or calibration errors exist. In this paper, we
show that under hardware impairments such as mutual coupling or antenna
separation disturbances, the array response vectors will no longer be the
Vandermonde vectors, and that different array response vectors should be
considered at every frequency for channel modeling under beam squint. In other
words, the assumptions and modeling of wideband mmWave MIMO channels in prior
work are not valid, and therefore, the prior CSI acquisition strategies are
not effective when beam squint and hardware impairments are considered.

The impact of beamsquint has been analyzed in prior work. As shown in [6, 7,
8], the time delay of the same data symbol across the antenna array aperture
is non-negligible in the large-scale MIMO configurations and/or the wideband
systems. Due to the spatial delay difference of each data symbol at different
antennas, the array steering vectors will have different responses at each
frequency, what leads to the beam squint effect. The work presented in [6]
includes the derivation of a channel model for the large-scale MIMO system
under beam squint, showing that the array response vectors for channel
modeling have to be frequency-dependent. Aware of this beam squint impact, [9]
and [10] also considered the frequency domain channel models by using explicit
frequency-dependent array response vectors for different subcarriers.
Meanwhile, to avoid the impacts of beam squint on channel estimation at
different subcarriers, [11, 12, 13] proposed to estimate the channel
parameters from the perspective of angle and delay domains, together with user
scheduling to alleviate inter-user interference. The works in [14] and [15]
consider a massive MIMO setting at mmWave or THz bands with single antenna
users, and propose solutions for beam tracking and channel estimation,
respectively, considering beam squint, but they do not introduce the filtering
effect or hardware impairments in the channel model. Generally, these prior
works on channel estimation/tracking under beam squint do not consider beam
squint combined with filtering effects or other hardware imperfections. In
this manuscript, we will show that beam squint not only leads to the
frequency-dependence on array steering vectors, but also yields additional
distortions at different antennas across all subcarriers, especially on those
at the band edge, when combined with the filtering operations at the
transceivers. The new channel model that results from these considerations has
not been derived in previous work. The pioneering work in [16] showed the need
and impact of considering the pulse shaping function in the MIMO channel
model, but beam squint was neglected. In summary, prior work did not consider
the combined effect of beam squint and the filtering effects, and also
neglected the impact of other hardware impairments such as calibration error,
mutual coupling or antenna separation disturbances when developing the signal
model and algorithms for a mmWave system operating with a hybrid MIMO
architecture.

The inclusion of hardware impairments on MIMO channel models has been
investigated in [17, 18, 19, 20]. There are many different hardware
impairments in the practical radio frequency (RF) chains, although three of
them dominate the effects in the resulting channel. First, due to the
manufacture and calibration errors, the antenna array will generate unexpected
radiation patterns, including both gain and phase errors on each antenna
element. Second, any perturbation on antenna locations or inter-element
spacing between antenna elements will result in irregular linear arrays rather
than perfect half-wavelength uniform linear arrays (ULAs). Finally, the
antenna spacing disturbance also creates the mutual coupling effect between
antenna elements. Taking into account all these hardware impairments, it is
apparent that the aforementioned sparsifying dictionaries constructed from
IARM evaluated at quantized angles are no longer the best choice for
exploiting channel sparsity. Previous work has shown that dictionary learning
(DL) is an effective technique to capture the underlying structure of mmWave
MIMO channels associated to specific types of sites [21] or with various
hardware impairments [18, 19, 20]. Specifically, [21] exploited the K-SVD
algorithm to find a dictionary to represent a collection of observed channel
realizations. Following this idea, [18] proposed a joint uplink/downlink
sparsifying DL algorithm for narrow band massive MIMO systems operating at
lower frequencies. Our previous work [19, 20] further investigated the DL and
channel estimation strategy for hybrid wideband mmWave MIMO systems under low
SNR conditions. Nevertheless, the beam squint effect was not incorporated in
previous work, and it cannot be ignored under certain relationships between
the carrier frequency and the bandwidth [7].

Motivated by these limitations, we propose a DL-based channel estimation
strategy for hybrid mmWave MIMO systems under the impacts of both hardware
impairments and beam squint. The main contributions of this manuscript are
summarized as follows:

  * •

We derive a general wideband mmWave MIMO channel model under both hardware
impairments and beam squint. The new model not only incorporates the hardware
impairments of antenna spacing disturbances, gain/phase errors, and array
mutual coupling, but also explicitly considers the impacts of combined pulse
shaping, filtering and beam squint, showing the limitations of existing MIMO
channel models with beam squint and the associated channel estimation schemes.
The derivation of the combined effect of beam squint and pulse
shaping/filtering is of particular interest. Previous literature does not
provide an alternative channel model including all these effects.

  * •

We propose a DL algorithm for finding sparsifying dictionaries that embed
hardware impairments, which fully exploits the channel properties at different
subcarriers under beam squint. This algorithm does not learn the impact of
beam squint but rather considers its model. Comparing to existing DL
strategies which has to learn a general dictionary for both hardware
impairments and beam squint, the newly proposed DL scheme enables better
adaptation to the impacts of hardware impairments and facilitates the
management of the beam squint impact at different subcarriers.

  * •

We design a novel orthogonal matching pursuit (OMP)-based algorithm for
compressive channel estimation under beam squint, which exploits the simple
structures of the new channel model at central subcarriers to obtain initial
parameter estimates with low complexity. Then it compensates the additional
distortions at side subcarriers induced by beam squint. In this way, the
measurements at all subcarriers can be used to achieve higher parameter
estimation accuracy with lower complexity.

  * •

We evaluate the proposed DL and channel estimation algorithms via numerical
simulations. Results show that the training overhead of channel estimation
with learned dictionaries can be significantly reduced compared to traditional
dictionaries without considering hardware impairment or beam squint. This
validates the developed channel models and corroborates the effectiveness of
the proposed DL algorithms for hybrid wideband mmWave MIMO systems under both
hardware impairments and beam squint effects.

Notations: Vectors and matrices are denoted by boldface small and capital
letters; the transpose, conjugate, Hermitian (conjugate transpose), inverse,
and pseudo-inverse of the matrix ${\mathbf{A}}$ are denoted by
${\mathbf{A}}^{T}$, ${\overline{{\mathbf{A}}}}$, ${\mathbf{A}}^{*}$,
${\mathbf{A}}^{-1}$ and ${\mathbf{A}}^{\dagger}$; ${{\mathbf{I}}_{M}}$ is an
$M\times M$ identity matrix; $[{\mathbf{a}}]_{n}$ denotes the $n$-th element
of ${\mathbf{a}}$ and $[{\mathbf{A}}]_{:,j}$ denotes the $j$-th column vector
of ${\mathbf{A}}$; $[{\mathbf{A}}]_{m,n}$ denotes the $(m,n)$-th element of
${\mathbf{A}}$; $\triangleq$ represents new definitions;
$\mathcal{I}(N)\triangleq\\{0,1,\ldots,N-1\\}$ denotes the index set of
cardinality $N$; ${\mathbb{C}}$ and ${\mathbb{R}}$ denote the sets of complex
and real numbers; $\text{tr}\\{{\mathbf{A}}\\}$ is the trace of
${\mathbf{A}}$; $\text{diag}\\{{\mathbf{a}}\\}$ denotes a diagonal matrix with
its diagonal elements given in ${\mathbf{a}}$ and
$\text{diag}\\{{\mathbf{A}}\\}$ formulates a vector by extracting the diagonal
elements of ${\mathbf{A}}$; $\otimes$, $\odot$ and $\star$ denote the
Kronecker, Hadamard and Khatri-Rao product between vectors/matrices;
$\text{j}=\sqrt{-1}$ denotes the imaginary unit; $|\mathcal{K}|$ denotes the
cardinality of a set; $\lfloor x\rfloor$ denotes the largest integer less than
or equal to $x$ and $\lceil x\rceil$ denotes the smallest integer great than
or equal to $x$; $\mathcal{F}(\cdot)[f]$ denotes the Fourier transform
evaluated at frequency $f$; and $\|{\mathbf{a}}\|_{0}$ denotes the $\ell_{0}$
norm of vector, i.e., the number of nonzero elements of ${\mathbf{a}}$.

##  II New Signal and Channel Models Incorporating Hardware Impairments and
Beam Squint

In this section, we derive the new general time domain signal model and
frequency domain channel model for the hybrid mmWave MIMO systems under the
impact of both hardware impairments and beam squint. We will prove that the
channel models with beam squint assumed in previous works [6, 10, 9, 11, 12,
13] are not complete or valid in some cases.

We consider a fully connected hybrid mmWave MIMO system.The transmitter (TX)
is equipped with ${N_{\text{T}}}$ antennas and ${L_{\text{T}}}$ RF chains, and
the receiver (RX) has ${N_{\text{R}}}$ antennas and ${L_{\text{R}}}$ RF
chains. The channel between the TX and the RX is frequency-selective. Pulse
shaped orthogonal frequency division multiplexing (OFDM) with $K$ subcarriers
is considered to simultaneously transmit
${N_{\text{s}}}~{}(\leq\min({L_{\text{T}}},{L_{\text{R}}}))$ data streams. The
system sampling period is denoted by ${T_{\text{s}}}$. We also consider a
pulse shaping function with a roll-off factor $\beta$ and an overall system
bandwidth of $(1+\beta)/{T_{\text{s}}}$. The center carrier frequency and
wavelength are denoted by ${f_{\text{c}}}$ and
${\lambda_{\text{c}}}=c/{f_{\text{c}}}$ (with $c$ the speed of light). We use
the index $k$ to denote the frequency domain subcarriers, with

$\displaystyle f_{k}$ $\displaystyle=f_{\text{c}}+\Delta
f_{k}=f_{\text{c}}-\frac{1}{2{T_{\text{s}}}}+\frac{k}{K{T_{\text{s}}}},\quad\text{for}\
k=0,\ldots,K-1.$ (1)

the corresponding subcarrier frequency value [22].

Let us start by defining a model for practical antenna arrays. As shown in
[23, 17, 18, 19, 24], various impairments exist in practical implementations
of antenna arrays, such as the gain/phase errors on each antenna element, the
errors at each antenna element location, and the coupling effects between
antenna elements. Specifically, we denote
$\mathbf{C}_{\text{R}}\in{\mathbb{C}}^{{N_{\text{R}}}\times{N_{\text{R}}}}$ as
the symmetric mutual coupling matrix for the RX antenna array, representing
the unwanted interchange of energy between elements in the arrays, and denote
$\bm{\gamma}_{\text{R}}\in\mathbb{C}^{N_{\text{R}}\times 1}$ with
$[\bm{\gamma}_{\text{R}}]_{n_{\text{R}}}=g_{\text{R},n_{\text{R}}}e^{\text{j}\nu_{\text{R},n_{\text{R}}}}$
as the antenna gain and phase errors, in which $g_{\text{R},n_{\text{R}}}$ and
$\nu_{\text{R},n_{\text{R}}}$ are the receive gain error normalized to a
reference amplitude, and the additional receive phase error for the
$n_{\text{R}}$-th antenna element. Moreover, let
$\bm{\epsilon}_{\text{R}}\in{\mathbb{R}}^{{N_{\text{R}}}\times 1}$ with
$[\bm{\epsilon}_{\text{R}}]_{n_{\text{R}}}=\epsilon_{\text{R},n_{\text{R}}}$
be the vector of antenna location errors at all RX antenna elements, where the
location error of the first antenna element is normalized to
$\epsilon_{\text{R},0}=0$. The variables defining the hardware impairments for
the TX antenna array are defined in an analogous way.

We define now the signal model for the transmitted signal during training. To
sound the channel, the TX sends $Q$ symbols over a ${L_{\text{T}}}\times 1$
signal vector ${\bf s}(t)$ at time instant $t$, where the $l_{\text{T}}$-th
($l_{\text{T}}=1,\ldots,{L_{\text{T}}}$) element of the complex exponential
representation of the transmitted signal ${\bf s}(t)$ is defined as

$\displaystyle s_{l_{\text{T}}}(t)\triangleq[{\bf
s}(t)]_{l_{\text{T}}}=\sum_{q=1}^{Q}[{\bf F}_{\rm BB}{\bf
S}]_{l_{\text{T}},q}p_{\rm T}(t-(q-1){T_{\text{s}}})e^{j2\pi f_{\rm c}t},$ (2)

where $p_{\rm T}(.)$ is the transmit pulse, $e^{j2\pi f_{\rm c}t}$ models the
RF upconversion stage, ${\bf S}\in\mathbb{C}^{{L_{\text{T}}}\times Q}$
collects all the time domain signal symbols to be transmitted and ${\bf
F}_{\rm BB}\in\mathbb{C}^{{L_{\text{T}}}\times{L_{\text{T}}}}$ is the digital
precoder. Moreover, ${\bf s}(t)$ satisfies ${\mathbb{E}}\\{{\bf s}(t){\bf
s}(t)^{*}\\}=\frac{P_{\text{T}}}{{L_{\text{T}}}}{\mathbf{I}}_{{L_{\text{T}}}}$,
with $P_{\text{T}}$ the transmit power constraint. Using an analog precoder
${\mathbf{F}}_{\rm RF}\in\mathbb{C}^{{N_{\text{T}}}\times{L_{\text{T}}}}$, the
ideal transmitted signal at the $n_{\rm T}$-th antenna element, $n_{\rm
T}=1,\ldots,{N_{\text{T}}}$, can be computed as

$\displaystyle i_{n_{\rm
T}}(t)=\sum_{l_{\text{T}}=1}^{{L_{\text{T}}}}[{\mathbf{F}}_{\rm RF}]_{n_{\rm
T},l_{\text{T}}}[{\bf s}(t)]_{l_{\text{T}}}.$ (3)

Due to hardware imperfections, the signal at the ${n_{\text{T}}}$-th transmit
antenna is multiplied by a phase and amplitude perturbation
$[\bm{\gamma}_{\text{T}}]_{n_{\text{T}}}^{*}$, and then leaks into the signal
at the ${n^{\prime}_{\text{T}}}$-th transmit antenna due to mutual coupling by
a factor $[{\bf C}_{\text{T}}]_{n_{\text{T}},n^{\prime}_{\text{T}}}^{*}$,
leading to the expression of the actual transmitted signal

$\displaystyle x_{n_{\rm T}}(t)=\sum_{n^{\prime}_{\rm
T}=1}^{{N_{\text{T}}}}[{\bf
C}_{\text{T}}]_{n^{\prime}_{\text{T}},n_{\text{T}}}^{*}[\bm{\gamma}_{\text{T}}]_{n^{\prime}_{\text{T}}}^{*}i_{n^{\prime}_{\rm
T}}(t)=\sum_{n^{\prime}_{\rm
T}=1}^{{N_{\text{T}}}}[\left(\text{diag}(\bm{\gamma}_{\text{T}}){\bf
C}_{\text{T}}\right)^{*}]_{n_{\rm T},n^{\prime}_{\rm T}}i_{n^{\prime}_{\rm
T}}(t).$ (4)

Assuming a passband geometric channel model with $L$ paths, the ideal received
signal at the $n_{\rm R}$-th antenna element, $n_{\rm
R}=1,\ldots,{N_{\text{R}}}$, can be written as

$\displaystyle i_{n_{\rm
R}}(t)=\sum_{\ell=1}^{L}\sum_{n_{\text{T}}=1}^{{N_{\text{T}}}}\alpha_{l}^{\rm
P}\cdot x_{n_{\rm T}}(t-\tau_{\ell,n_{\text{R}},n_{\text{T}}})+z_{n_{\rm
R}}(t),$ (5)

where $z_{n_{\rm R}}(t)$ is the noise term, $\alpha_{l}^{\rm
P}\in{\mathbb{C}}$ is the path gain of the passband channel including phase
changes, and $\tau_{\ell,n_{\text{R}},n_{\text{T}}}$ is the delay of $\ell$-th
path between the $n_{\text{R}}$-th receive antenna and the $n_{\text{T}}$-th
transmit antenna. Similarly to what happens at the transmit array, owing to
hardware imperfections, the signal at each receive antenna ${n_{\text{R}}}$
leaks into the signal at each antenna ${n^{\prime}_{\text{R}}}$ through mutual
coupling by a factor $[{\bf
C}_{\text{R}}]_{n^{\prime}_{\text{R}},n_{\text{R}}}$, and then receives a
gain/phase weight $[\bm{\gamma}_{\text{R}}]_{n^{\prime}_{\text{R}}}$ to
generate the actual receive signal

$\displaystyle y_{n_{\rm R}}(t)=\sum_{n^{\prime}_{\rm
R}=1}^{{N_{\text{R}}}}[\bm{\gamma}_{\text{R}}]_{n_{\text{T}}}[{\bf
C}_{\text{R}}]_{n_{\text{R}},n^{\prime}_{\text{R}}}i_{n^{\prime}_{\rm
R}}(t)=\sum_{n^{\prime}_{\rm
R}=1}^{{N_{\text{R}}}}[\text{diag}(\bm{\gamma}_{\text{R}}){\bf
C}_{\text{R}}]_{n_{\text{R}},n^{\prime}_{\text{R}}}i_{n^{\prime}_{\rm R}}(t).$
(6)

After an analog combining stage modeled by the matrix ${{\mathbf{W}}}_{\rm
RF}\in\mathbb{C}^{N_{\rm R}\times L_{\rm R}}$, we obtain the received signal
vector $\bm{y}^{\prime}(t)\in{\mathbb{C}}^{{L_{\text{R}}}\times 1}$, with the
$l_{\text{R}}$-th element given as

$\displaystyle[\bm{y}^{\prime}(t)]_{l_{\rm R}}=\sum_{n_{\rm
R}=1}^{{N_{\text{R}}}}[{\bf W}_{\rm RF}]_{n_{\rm R},l_{\rm R}}^{*}y_{n_{\rm
R}}(t).$ (7)

To obtain the downconverted signal we get rid of the complex representation of
the signal and multiply by $2e^{-j2\pi f_{\rm c}t}$ to get
$[\bm{y}^{\prime\prime}]_{l^{\prime}_{\rm
R}}(t)=([\bm{y}^{\prime}]_{l^{\prime}_{\rm
R}}(t)+[\bm{y}^{\prime}]_{l^{\prime}_{\rm R}}^{*}(t))e^{-j2\pi f_{\rm c}t}$.
After matched filtering and sampling with sampling period $T_{\rm s}$, we have

$\displaystyle\big{[}{\bm{y}}_{\rm BB}[d]\big{]}_{{l_{\rm R}}}=(p_{\rm
R}*[\bm{y}^{\prime\prime}]_{l^{\prime}_{\rm R}})((d-1){T_{\text{s}}}),\
\text{for}\ d=0,1,\ldots$ (8)

where $p_{\rm R}$ is the matched filter. Since the bandwidths of $p_{\rm T}$
and $p_{\rm R}$, namely $B_{\rm T}$ and $B_{\rm R}$, satisfy $B_{\rm T}+B_{\rm
R}\leq 2f_{\rm c}$ and $[\bm{y}^{\prime}]_{l^{\prime}_{\rm R}}^{*}(t)e^{-j2\pi
f_{\rm c}t}$ is centered at $-2f_{\rm c}$, the convolution between the matched
filter and the conjugate signal is zero, i.e. $p_{\rm
R}(t)*([\bm{y}^{\prime}]_{l^{\prime}_{\rm R}}^{*}(t)e^{-j2\pi f_{\rm c}t})=0$.
After the digital combining stage represented by ${\bf W}_{\rm BB}$, the
discrete time received signal can be written as

$\displaystyle\big{[}{\bm{y}}[d]\big{]}_{{l_{\rm R}}}=\sum_{l^{\prime}_{\rm
R}}[{\bf W}_{\rm BB}]_{l^{\prime}_{\rm R},l_{\rm R}}^{*}\big{[}{\bm{y}}_{\rm
BB}[d]\big{]}_{{l_{\rm R}}},\ \text{for}\ d=0,1,\ldots$ (9)

By developing the previous steps, we reach an equivalent expression for the
received signal

$\displaystyle\big{[}{\bm{y}}[d]\big{]}_{{l_{\rm R}}}=\sum_{\ell,n_{\rm
R},n_{\rm T},l_{\rm T}}\alpha_{l}^{\rm P}[\breve{\bf W}]_{n_{\rm R},l_{\rm
R}}^{*}[\breve{\bf F}]_{n_{\rm T},l_{\rm T}}\sum_{q=1}^{Q}p((d-q)T_{\rm
s}-{\tau}_{\ell,n_{\rm R},n_{\rm T}})e^{-j2\pi f_{\rm c}{\tau}_{\ell,n_{\rm
R},n_{\rm T}}}[{\mathbf{s}}[q]]_{{l_{\rm T}}}+[\breve{\bf W}]_{n_{\rm
R},l_{\rm R}}^{*}z_{n_{\rm R}}(dT_{\rm s}),$ (10)

where $\breve{\bf W}\triangleq\left(\text{diag}(\bm{\gamma}_{\text{R}}){\bf
C}_{\text{R}}\right)^{*}{\bf W}$ and $\breve{\bf
F}\triangleq\left(\text{diag}(\bm{\gamma}_{\text{T}}){\bf
C}_{\text{T}}\right)^{*}{\bf F}$, being ${\bf W}={\bf W}_{\rm RF}{\bf W}_{\rm
BB}$ and ${\bf F}={\bf F}_{\rm RF}{\bf F}_{\rm BB}$ the hybrid combiner and
precoder, respectively, while $p=p_{\rm R}*p_{\rm T}$ combines the effects of
transmit and receive filters.

(a) Propagation of the $\ell$-th path between RX and TX with AoA
${\phi_{\ell}}$ and AoD ${\theta_{\ell}}$.

(b) Received signals at RX antennas, where the delay of each path varies
across array geometry.

Figure 1: Illustration of the spatial wideband effect (from the perspective of
time domain) or the beam squint effect (from the perspective of frequency
domain) under large-scale antenna array regimes.

Next, following the derivation in the Appendix I, the DFT of (10) is computed
as

$\displaystyle[{\bf y}[k]]_{{l_{\rm R}}}=\sum_{\ell,n_{\rm R},n_{\rm T},l_{\rm
T}}\alpha_{l}^{\rm P}[\breve{\bf W}]_{n_{\rm R},l_{\rm R}}^{*}[\breve{\bf
F}]_{n_{\rm T},l_{\rm T}}g(k,{\tau}_{\ell,n_{\rm R},n_{\rm
T}})e^{-\text{j}2\pi(f_{\rm c}+\Delta f_{k}){\tau}_{\ell,n_{\rm R},n_{\rm
T}}}[{\mathbf{s}}[k]]_{{l_{\rm T}}}+[\breve{\bf W}]_{n_{\rm R},l_{\rm
R}}^{*}[{\mathbf{z}}[k]]_{n_{{\rm R}}},$ (11)

where $g(k,{\tau}_{\ell,n_{\rm R},n_{\rm T}})$ is defined in (53) in the
Appendix I, representing the frequency response at the $k$-th subcarrier of
the combined filter $p$ for the $\ell$-th path between the $n_{\text{R}}$-th
RX antenna and the $n_{\text{T}}$-th TX antenna. Moreover,
$[{\mathbf{s}}[k]]_{{l_{\rm T}}}$ is the frequency domain version of the
transmitted signal defined in (2), while $[{\mathbf{z}}[k]]_{n_{\rm
R}}\triangleq\sum_{d=0}^{K-1}z_{n_{\rm R}}(dT_{\rm s})e^{-\text{j}\frac{2\pi
dk}{K}}$ is the noise term in the frequency domain. Regarding the delay
${\tau}_{\ell,n_{\rm R},n_{\rm T}}$, as shown in Fig. 1, it can be expressed
as

$\displaystyle{\tau}_{\ell,n_{\rm R},n_{\rm T}}={\tau_{\ell}}+(n_{\rm
R}d_{\text{R}}+\epsilon_{\text{R},n_{\rm R}})\sin({\phi_{\ell}})/c-(n_{\rm
T}d_{\text{T}}+\epsilon_{\text{T},n_{\rm T}})\sin({\theta_{\ell}})/c.$ (12)

Now, by defining the frequency-dependent array steering vectors

$\displaystyle\breve{{\mathbf{a}}}_{\text{R},k}(\phi)=\text{diag}(\bm{\gamma}_{\text{R}}){\bf
C}_{\text{R}}\check{{\mathbf{a}}}_{\text{R},k}(\phi),$ (13)
$\displaystyle\breve{{\mathbf{a}}}_{\text{T},k}(\theta)=\text{diag}(\bm{\gamma}_{\text{T}}){\bf
C}_{\text{T}}\check{{\mathbf{a}}}_{\text{T},k}(\theta),$ (14)

with

$\displaystyle[\check{{\mathbf{a}}}_{\text{R},k}(\phi)]_{n}=\frac{1}{\sqrt{N_{\rm
R}}}e^{-\text{j}2\pi
f_{k}(n{d_{\text{R}}}+\epsilon_{\text{R},n})\sin(\phi)/c}$ (15)
$\displaystyle[\check{{\mathbf{a}}}_{\text{T},k}(\theta)]_{n}=\frac{1}{\sqrt{N_{\rm
T}}}e^{-\text{j}2\pi
f_{k}(n{d_{\text{T}}}+\epsilon_{\text{T},n})\sin(\theta)/c},$ (16)

denoting the distortion matrix by

$\displaystyle[{\bf
G}_{k}({\tau_{\ell}},{\phi_{\ell}},{\theta_{\ell}})]_{n_{\rm R},n_{\rm
T}}=g(k,{\tau}_{\ell,n_{\rm R},n_{\rm T}}),$ (17)

and defining the quivalent complex gain $\alpha_{l}=\alpha_{l}^{\rm
P}e^{-\text{j}2\pi f_{\rm c}\tau_{l}}$, we can rewrite (11) as

$\displaystyle{\bf y}[k]={\bf
W}^{*}\sum_{\ell=1}^{L}\left({\alpha_{\ell}}e^{-\text{j}2\pi\Delta
f_{k}{\tau_{\ell}}}{\bf
G}_{k}({\tau_{\ell}},{\phi_{\ell}},{\theta_{\ell}})\odot(\breve{{\mathbf{a}}}_{\text{R},k}({\phi_{\ell}})\breve{{\mathbf{a}}}_{\text{T},k}^{*}({\theta_{\ell}}))\right){\bf
F}\hat{\bf s}[k]+\breve{\bf W}^{*}{\bf z}[k].$ (18)

Based on this expression, we can define the general frequency domain channel
matrix as

$\displaystyle{\mathbf{H}}[k]=\sum_{\ell=1}^{L}{\alpha_{\ell}}e^{-\text{j}2\pi\Delta
f_{k}{\tau_{\ell}}}{\bf
G}_{k}({\tau_{\ell}},{\phi_{\ell}},{\theta_{\ell}})\odot(\breve{{\mathbf{a}}}_{\text{R},k}({\phi_{\ell}})\breve{{\mathbf{a}}}_{\text{T},k}^{*}({\theta_{\ell}})),$
(19)

and then, the final expression of the received signal under these definitions
becomes

$\displaystyle{\bf y}[k]={\bf W}^{*}{\mathbf{H}}[k]{\bf F}{\bf
s}[k]+\breve{\bf W}^{*}{\mathbf{z}}[k].$ (20)

Note that the equivalent combiner that impacts the noise term includes the
hardware impairments as well, and thereby the noise covariance matrix is
${\mathbf{C}}_{{\mathbf{z}}[k]}=\sigma^{2}\breve{\bf W}^{*}\breve{\bf W}$.

As per the definition of $g(k,\tau)$ in the Appendix I, i.e.,

$\displaystyle
g(k,\tau)=\left\\{\begin{array}[]{cl}\frac{\mathcal{F}(p)[\Delta
f_{k}]}{T_{\rm s}}+\frac{\mathcal{F}(p)[\Delta f_{k}-1/T_{\rm s}]}{T_{\rm
s}}e^{\text{j}2\pi\tau/T_{\rm s}}&\text{if}\ \ 2\Delta f_{k}>(1-\beta)/T_{\rm
s},\\\ \frac{\mathcal{F}(p)[\Delta f_{k}]}{T_{\rm s}}&\text{if}\ \ 2|\Delta
f_{k}|\leq(1-\beta)/T_{\rm s},\\\ \frac{\mathcal{F}(p)[\Delta f_{k}]}{T_{\rm
s}}+\frac{\mathcal{F}(p)[\Delta f_{k}+1/T_{\rm s}]}{T_{\rm
s}}e^{-\text{j}2\pi\tau/T_{\rm s}}&\text{if}\ \ 2\Delta
f_{k}<-(1-\beta)/T_{\rm s},\end{array}\right.$ (24)

where $\mathcal{F}(p)[\Delta f_{k}]$ is the Fourier transform of the combined
filter $p(t)$ evaluated for the frequency difference $\Delta f_{k}$, the
general channel matrix in (19) can be rewritten as

$\displaystyle{\mathbf{H}}[k]=\begin{cases}\sum_{\ell=1}^{L}{\alpha_{\ell}}e^{-\text{j}2\pi\Delta
f_{k}{\tau_{\ell}}}{\bf
G}_{k}({\tau_{\ell}},{\phi_{\ell}},{\theta_{\ell}})\odot(\breve{{\mathbf{a}}}_{\text{R},k}({\phi_{\ell}})\breve{{\mathbf{a}}}_{\text{T},k}^{*}({\theta_{\ell}})),\quad\forall|\Delta
f_{k}|>(1-\beta)/2{T_{\text{s}}}\\\
\sum_{\ell=1}^{L}{\alpha_{\ell}}\frac{\mathcal{F}(p)[\Delta
f_{k}]}{{T_{\text{s}}}}e^{-\text{j}2\pi\Delta
f_{k}{\tau_{\ell}}}\breve{{\mathbf{a}}}_{\text{R},k}({\phi_{\ell}})\breve{{\mathbf{a}}}_{\text{T},k}^{*}({\theta_{\ell}}),\quad\forall|\Delta
f_{k}|\leq(1-\beta)/2{T_{\text{s}}}\end{cases}$ (25)

where the second case is due to the fact that when $|\Delta
f_{k}|\leq(1-\beta)/2{T_{\text{s}}}$, the factor $g(k,{\tau}_{\ell,n_{\rm
R},n_{\rm T}})=\frac{\mathcal{F}(p)[\Delta f_{k}]}{{T_{\text{s}}}}$ is
independent of the delay ${\tau}_{\ell,n_{\rm R},n_{\rm T}}$ at different
antenna indices, and thus enables the simplification. Note that (25) is valid
for any pulse shaping function (as shown in the derivation in Appendix I), and
therefore, this does not introduce any new assumption on the channel model.

######  Remark 1.

We have derived the new general model for MIMO channels under both hardware
impairments and beam squint. This new model shows that the beam squint not
only induces the frequency-dependence on array steering vectors, but also
yields additional distortions at different antennas across all subcarriers,
especially on side subcarriers. Comparing (25) to the existing MIMO channel
models with beam squint in [6, 10, 9, 11, 12, 13], it is obvious that those
models only considered the second case in (25), and assumed the same
expression for all subcarriers. In other words, the prior work only considered
the frequency-dependence impact of beam squint, but ignored the additional
distortions. This is due to the fact that the model was derived in the
continuous time domain, and did not take into account the impact of the extra
bandwidth of the combined filter (as shown in Appendix I). For ease of
subsequent exposition, we will denote the set of central subcarriers by
${\mathcal{K}}_{\text{cen}}=\\{k\big{|}|\Delta
f_{k}|\leq(1-\beta)/2{T_{\text{s}}}\\}$ and the set of side subcarriers by
${\mathcal{K}}_{\text{side}}=\\{k\big{|}|\Delta
f_{k}|>(1-\beta)/2{T_{\text{s}}}\\}$. As shown in Fig. 2, there are
approximately $|{\mathcal{K}}_{\text{cen}}|\approx\lfloor(1-\beta)K\rfloor$
central subcarriers inside ${\mathcal{K}}_{\text{cen}}$ and
$|{\mathcal{K}}_{\text{side}}|\approx\lceil\beta K\rceil$ side subcarriers
inside ${\mathcal{K}}_{\text{side}}$.

Figure 2: Frequency response of the pulse shaping function with roll-off
factor $\beta$, center frequency ${f_{\text{c}}}$ and
${f_{\text{s}}}=1/{T_{\text{s}}}$. For OFDM systems with $K$ subcarriers,
there are approximately $\lfloor(1-\beta)K\rfloor$ central subcarriers inside
${\mathcal{K}}_{\text{cen}}$, and $\lceil\beta K\rceil$ side subcarriers
inside ${\mathcal{K}}_{\text{side}}$.

###  II-A Frequency domain signal model for training

With the aforementioned system and channel model, the received signal model
for training is described as follows. During the training phase, the TX will
send pilot signals over several OFDM frames, so the RX can collect
measurements for initial channel estimation and DL. Specifically, during the
$m$-th ($m=1,2,\ldots$) training OFDM symbol, the TX will send an
${L_{\text{T}}}\times 1$ pilot signal ${\mathbf{s}}_{m}[k]$ at the $k$-th
subcarrier using a frequency-flat precoder
${{\mathbf{F}}}_{m}\in{\mathbb{C}}^{{N_{\text{T}}}\times{L_{\text{T}}}}$,
while the RX employs a frequency-flat combiner
${{\mathbf{W}}}_{m}\in{\mathbb{C}}^{{N_{\text{R}}}\times{L_{\text{R}}}}$. The
received signal is then given as

$\displaystyle{\mathbf{y}}_{m}[k]={\mathbf{W}}_{m}^{*}{\mathbf{H}}[k]{\mathbf{F}}_{m}{{\mathbf{s}}}_{m}[k]+{{\mathbf{n}}}_{m}[k],$
(26)

where ${\mathbf{n}}_{m}[k]\in{\mathbb{C}}^{{L_{\text{R}}}\times 1}$ is the
additive Gaussian noise vector, distributed as $\mathcal{C}\mathcal{N}(\bf
0,\sigma^{2}\breve{{\mathbf{W}}}^{*}_{m}\breve{{\mathbf{W}}}_{m})$. Generally,
we can decompose the transmitted signal as
${\mathbf{s}}_{m}[k]\triangleq{\mathbf{q}}_{m}s[k]$, with
${\mathbf{q}}_{m}\in{\mathbb{C}}^{{L_{\text{T}}}\times 1}$ a frequency-flat
training vector and $s[k]$ a frequency-dependent training symbol. In doing so,
we can multiply the received signal ${\mathbf{y}}_{m}[k]$ by $(s[k])^{-1}$ and
get a frequency-flat observation matrix at the RX as follows

$\displaystyle\tilde{{\mathbf{y}}}_{m}[k]\triangleq\mathrm{vec}\left({\mathbf{y}}_{m}[k](s[k])^{-1}\right)=({\mathbf{q}}_{m}^{T}{\mathbf{F}}_{m}^{T}\otimes{\mathbf{W}}_{m}^{*})\mathrm{vec}({\mathbf{H}}[k])+\tilde{{\mathbf{n}}}_{m}[k]={\mathbf{\Phi}}_{m}\mathrm{vec}({\mathbf{H}}[k])+\tilde{{\mathbf{n}}}_{m}[k],$
(27)

where
$\tilde{{\mathbf{n}}}_{m}[k]\triangleq\mathrm{vec}\big{(}{{\mathbf{n}}}_{m}[k](s[k])^{-1}\big{)}$
and ${\mathbf{\Phi}}_{m}$ is defined accordingly, representing the sensing
matrix for the $m$-th OFDM symbol. To get a higher effective SNR for the
received measurements, we will use training spreading to average out the noise
[19]. Moreover, to enable the initial channel estimate, the measurements over
$M$ OFDM symbols are stacked together such that

$\displaystyle\underbrace{\left[\tilde{{\mathbf{y}}}_{1}[k]^{T},\ldots,\tilde{{\mathbf{y}}}_{M}[k]^{T}\right]^{T}}_{\tilde{{\mathbf{y}}}[k]}=\underbrace{\left[{\mathbf{\Phi}}_{1}^{T},\ldots,{\mathbf{\Phi}}_{M}^{T}\right]^{T}}_{{\mathbf{\Phi}}}{\mathbf{h}}[k]+\underbrace{\left[\tilde{{\mathbf{n}}}_{1}[k]^{T},\ldots,\tilde{{\mathbf{n}}}_{M}[k]^{T}\right]^{T}}_{\tilde{{\mathbf{n}}}[k]},$
(28)

where $\tilde{{\mathbf{y}}}[k]\in{\mathbb{C}}^{M{L_{\text{R}}}\times 1}$,
${\mathbf{\Phi}}\in{\mathbb{C}}^{M{L_{\text{R}}}\times{N_{\text{R}}}{N_{\text{T}}}}$
and $\tilde{{\mathbf{n}}}[k]\in{\mathbb{C}}^{M{L_{\text{R}}}\times 1}$ are
defined accordingly. Therefore, we can obtain the initial least squares (LS)
channel estimate as
$\hat{{\mathbf{h}}}[k]={\mathbf{\Phi}}^{{\dagger}}\tilde{{\mathbf{y}}}[k]$,
while the corresponding channel matrix
$\hat{{\mathbf{H}}}[k]=\text{unvec}(\hat{{\mathbf{h}}}[k])$.

During the training phase, it becomes necessary to collect initial channel
measurements at different locations across the coverage area as the RX moves
around. This not only creates a large training data set for learning the
hardware imperfections, but also ensures the data set is diverse for the
environment. In doing so, the learned sparsifying dictionary is not dedicated
for a specific location, but adapted to the fixed hardware impairments. As in
[19], this procedure of collecting measurements at multiple locations can be
done at the stage of network setup, like in the case of an indoor WiFi
scenario. Therefore, we assume the measurements are collected at
${N_{\text{sa}}}$ locations and the initial channel estimate at the $u$-th
location is then represented as $\hat{{\mathbf{h}}}^{(u)}[k]$ for
$u\in\mathcal{I}({N_{\text{sa}}})$. In the following, these initial channel
estimates obtained at the training phases will be used for learning the
hardware impairments related dictionaries.

##  III Dictionary Learning for Hardware Impairments under beam squint
(DLHWBS)

In this section, we present the problem formulation and optimization for
learning the TX and RX dictionaries adapted to hardware impairments.
Generally, the overall dictionaries are partitioned into two parts, i.e., one
frequency-flat part adapting to the hardware impairments and one frequency-
dependent part accounting for the beam squint effect. The learning of hardware
impairment related dictionaries will take into consideration of the channel
properties under beam squint at different subcarriers.

###  III-A Formulation of the dictionary learning problem

As shown in (25), the frequency domain channel model can be partitioned into
two sets. We treat each one separately here. For central subcarriers inside
${\mathcal{K}}_{\text{cen}}$, the channel model can be simplified so that it
is more convenient and straightforward to separate the hardware impairments
and beam squint impact on the channel matrix. Let us first rewrite the channel
matrix for $k\in{\mathcal{K}}_{\text{cen}}$ as follows

$\displaystyle{\bf H}[k]$
$\displaystyle=\sum_{\ell=1}^{L}{\alpha_{\ell}}\frac{\mathcal{F}(p)[\Delta
f_{k}]}{{T_{\text{s}}}}e^{-\text{j}2\pi\Delta
f_{k}{\tau_{\ell}}}\breve{{\mathbf{a}}}_{\text{R},k}({\phi_{\ell}})\breve{{\mathbf{a}}}^{*}_{\text{T},k}({\theta_{\ell}})=\breve{{\mathbf{A}}}_{\text{R},k}(\bm{\phi})\bm{\Lambda}[k]\breve{{\mathbf{A}}}_{\text{T},k}^{*}(\bm{\theta})$
$\displaystyle\approx\breve{{\mathbf{A}}}^{\rm
v}_{\text{R},k}\bm{\Omega}[k]\breve{{\mathbf{A}}}_{\text{T},k}^{\rm v*},$ (29)

where
$\bm{\Lambda}[k]\triangleq\text{diag}\big{\\{}[\alpha_{1}\frac{\mathcal{F}(p)[\Delta
f_{k}]}{{T_{\text{s}}}}e^{-\text{j}2\pi\Delta
f_{k}\tau_{1}},\ldots,\alpha_{L}\frac{\mathcal{F}(p)[\Delta
f_{k}]}{{T_{\text{s}}}}e^{-\text{j}2\pi\Delta f_{k}\tau_{L}}]^{T}\big{\\}}$ is
a diagonal matrix containing the frequency domain path gains, and
$\breve{{\mathbf{A}}}_{\text{R},k}(\bm{\phi})\triangleq\left[\breve{{\mathbf{a}}}_{\text{R},k}(\phi_{1}),\ldots,\breve{{\mathbf{a}}}_{\text{R},k}(\phi_{L})\right]\in{\mathbb{C}}^{{N_{\text{R}}}\times
L}$ collects the receive array response vectors at all paths. Moreover, the
approximation is obtained by discretizing the AoA/AoD spaces with on-grid
angles, i.e., $\breve{{\mathbf{A}}}^{\rm
v}_{\text{R},k}\triangleq[\breve{\mathbf{a}}_{\text{R},k}(\phi^{\rm
v}_{1}),\ldots,\allowbreak\breve{\mathbf{a}}_{\text{R},k}(\phi^{\rm
v}_{{K_{\text{R}}}})]\in{\mathbb{C}}^{{N_{\text{R}}}\times{K_{\text{R}}}}$ and
$\breve{{\mathbf{A}}}_{\text{T},k}^{\rm
v}\triangleq[\breve{\mathbf{a}}_{\text{T},k}(\theta^{\rm
v}_{1}),\ldots,\allowbreak\breve{\mathbf{a}}_{\text{T},k}(\theta^{\rm
v}_{{K_{\text{T}}}})]\in{\mathbb{C}}^{{N_{\text{T}}}\times{K_{\text{T}}}}$
collect the virtual receive and transmit array response vectors evaluated on
${K_{\text{R}}}$ quantized angles $\bm{\phi}^{\rm v}=\\{\phi^{\rm
v}_{i}\\}_{i=1}^{{K_{\text{R}}}}$ for AoAs and ${K_{\text{T}}}$ quantized
angles $\bm{\theta}^{\rm v}=\\{\theta^{\rm v}_{i}\\}_{i=1}^{{K_{\text{T}}}}$
for AoDs, and
$\bm{\Omega}[k]\in{\mathbb{C}}^{{K_{\text{R}}}\times{K_{\text{T}}}}$ is a
sparse matrix, containing the path gains of these discrete quantized AoAs/AoDs
at its non-zero elements. To separate the impact of hardware impairments and
beam squint on channel matrix, we can rewrite the array steering vector in
(15) as

$\displaystyle[\check{{\mathbf{a}}}_{\text{R},k}(\phi)]_{n_{\text{R}}}$
$\displaystyle=\frac{1}{\sqrt{{N_{\text{R}}}}}e^{-\text{j}2\pi
f_{k}(n_{\text{R}}{d_{\text{R}}}+\epsilon_{\text{R},n_{\text{R}}})\sin(\phi)/c}$
$\displaystyle=\frac{1}{\sqrt{{N_{\text{R}}}}}e^{-\text{j}2\pi(f_{\rm
c}\epsilon_{\text{R},n_{\text{R}}}+\Delta
f_{k}\epsilon_{\text{R},n_{\text{R}}})\sin(\phi)/c}e^{-\text{j}2\pi(f_{k}n_{\text{R}}{d_{\text{R}}})\sin(\phi)/c}$
$\displaystyle\approx
e^{-\text{j}2\pi({f_{\text{c}}}\epsilon_{\text{R},n_{\text{R}}})\sin(\phi)/c}\cdot\frac{1}{\sqrt{{N_{\text{R}}}}}e^{-\text{j}2\pi(f_{k}n_{\text{R}}{d_{\text{R}}})\sin(\phi)/c},$
(30) $\displaystyle\check{{\mathbf{a}}}_{\text{R},k}(\phi)$
$\displaystyle\approx{\bf e}_{\text{R}}(\phi)\odot{\mathbf{a}}_{{\rm
R},k}(\phi),$ (31)

where $[{\bf e}_{\text{R}}(\phi)]_{n_{\text{R}}}=e^{-\text{j}2\pi(f_{\rm
c}\epsilon_{\text{R},n_{\text{R}}})\sin(\phi)/c}$ comprises the hardware
antenna element location imperfection effects. Therefore, the expression of
$\breve{{\mathbf{A}}}^{\rm v}_{\text{R},k}$ in (III-A) can be rewritten as

$\displaystyle\breve{{\mathbf{A}}}^{\rm v}_{\text{R},k}$
$\displaystyle\triangleq[\breve{\mathbf{a}}_{\text{R},k}(\phi^{\rm
v}_{1}),\ldots,\breve{\mathbf{a}}_{\text{R},k}(\phi^{\rm
v}_{{K_{\text{R}}}})]$
$\displaystyle=\underbrace{\text{diag}(\bm{\gamma}_{\text{R}}){\bf
C}_{\text{R}}}_{{\mathbf{D}}_{\text{R},1}}\bigg{(}\underbrace{\left[{\bf
e}_{\text{R}}(\phi^{\rm v}_{1}),\ldots,{\bf e}_{\text{R}}(\phi^{\rm
v}_{{K_{\text{R}}}})\right]}_{{\mathbf{D}}_{\text{R},2}}\odot\underbrace{\left[{\mathbf{a}}_{\text{R},k}(\phi^{\rm
v}_{1}),\ldots,{\mathbf{a}}_{\text{R},k}(\phi^{\rm
v}_{{K_{\text{R}}}})\right]}_{{\mathbf{A}}_{\text{R},k}^{\rm v}}\bigg{)}$
$\displaystyle={\mathbf{D}}_{\text{R},1}\left({\mathbf{D}}_{\text{R},2}\odot{\mathbf{A}}_{\text{R},k}^{\rm
v}\right),$ (32)

where the general RX dictionary $\breve{{\mathbf{A}}}^{\rm v}_{\text{R},k}$ is
partitioned into three parts, i.e., ${\mathbf{D}}_{\text{R},1}$ includes the
impacts of antenna coupling and gain/phase errors, ${\mathbf{D}}_{\text{R},2}$
accounts for the antenna location perturbations, and
${\mathbf{A}}_{\text{R},k}^{\rm
v}\in{\mathbb{C}}^{{N_{\text{R}}}\times{K_{\text{R}}}}$ handles the beam
squint effect at each subcarriers. Note that ${\mathbf{D}}_{\text{R},1}$ and
${\mathbf{D}}_{\text{R},2}$ are hardware impairments related dictionary
components to be learned while ${\mathbf{A}}_{\text{R},k}^{\rm v}$ is known
for all subcarriers. Moreover, recalling the notation of antenna location
error vector $\bm{\epsilon}_{\text{R}}$ and the definition of $\bf
e_{\text{R}}(\phi)$ in (III-A), we can further express
${\mathbf{D}}_{\text{R},2}$ as a function of the antenna location error vector
$\bm{\epsilon}_{\text{R}}$, i.e.,
${\mathbf{D}}_{\text{R},2}=f(\bm{\epsilon}_{\text{R}})=e^{-\text{j}2\pi{f_{\text{c}}}\bm{\epsilon}_{\text{R}}\cdot\sin(\bm{\phi}^{\rm
v})^{T}/c}$. Similarly, the general TX dictionary $\breve{{\mathbf{A}}}^{\rm
v}_{\text{T},k}$ can be expressed as $\breve{{\mathbf{A}}}^{\rm
v}_{\text{T},k}={\mathbf{D}}_{\text{T},1}\left({\mathbf{D}}_{\text{T},2}\odot{\mathbf{A}}_{\text{T},k}^{\rm
v}\right)$ and ${\mathbf{D}}_{\text{T},2}$ is a function of the antenna
location error vector $\bm{\epsilon}_{\text{T}}$, i.e.,
${\mathbf{D}}_{\text{T},2}=f(\bm{\epsilon}_{\text{T}})=e^{-\text{j}2\pi{f_{\text{c}}}\bm{\epsilon}_{\text{T}}\cdot\sin(\bm{\theta}^{\rm
v})^{T}/c}$. Then the vectorization of
${\mathbf{H}}[k],k\in{\mathcal{K}}_{\text{cen}}$ in (III-A) can be expressed
as

$\displaystyle{\mathbf{h}}[k]$
$\displaystyle=\mathrm{vec}({\mathbf{H}}[k])\approx\left(\overline{\breve{{\mathbf{A}}}_{\text{T},k}^{\rm
v}}\otimes\breve{{\mathbf{A}}}^{\rm
v}_{\text{R},k}\right)\mathrm{vec}(\bm{\Omega}[k])$
$\displaystyle=\left(\overline{{\mathbf{D}}_{\text{T},1}}\otimes{\mathbf{D}}_{\text{R},1}\right)\left[\left(\overline{{\mathbf{D}}_{\text{T},2}}\otimes{\mathbf{D}}_{\text{R},2}\right)\odot\left(\overline{{\mathbf{A}}_{\text{T},k}^{\rm
v}}\otimes{\mathbf{A}}_{\text{R},k}^{\rm
v}\right)\right]\mathrm{vec}(\bm{\Omega}[k]).$ (33)

The channel model for side subcarriers $k\in{\mathcal{K}}_{\text{side}}$ in
(25) can be represented and approximated as follows

$\displaystyle{\mathbf{H}}[k]$
$\displaystyle=\sum_{\ell=1}^{L}{\alpha_{\ell}}e^{-\text{j}2\pi\Delta
f_{k}{\tau_{\ell}}}{\bf
G}_{k}({\tau_{\ell}},{\phi_{\ell}},{\theta_{\ell}})\odot(\breve{{\mathbf{a}}}_{\text{R},k}({\phi_{\ell}})\breve{{\mathbf{a}}}^{*}_{\text{T},k}({\theta_{\ell}}))$
$\displaystyle\approx\sum_{l=1}^{L_{\rm
v}}\sum_{i=1}^{{K_{\text{R}}}}\sum_{j=1}^{{K_{\text{T}}}}b_{l,i,j}e^{-\text{j}2\pi\Delta
f_{k}\tau_{l}^{\rm v}}{\bf G}_{k}(\tau_{l}^{\rm v},\phi_{i}^{\rm
v},\theta_{j}^{\rm v})\odot(\breve{{\mathbf{a}}}_{\text{R},k}(\phi_{i}^{\rm
v})\breve{{\mathbf{a}}}^{*}_{\text{T},k}(\theta_{j}^{\rm v}))$ (34)

where the approximation is obtained by first discretizing the AoA/AoD spaces
as did in (III-A) and then discretizing the delay space with $L_{\rm v}$ on-
grid delays $\\{\tau_{l}^{\rm v}\\}_{l=1}^{L_{\rm v}}$. Moreover, $b_{l,i,j}$
is the corresponding path gain on each pair of discretized AoA/AoD/delay
grids, which theoretically is nonzero only at the discrete AoA/AoD/delay bin
$(\phi^{\text{v}}_{i},\theta^{\text{v}}_{j},\tau_{l}^{\rm v})$ corresponding
to the true AoA/AoD/delay $({\phi_{\ell}},{\theta_{\ell}},{\tau_{\ell}})$.
Then the vectorization of ${\mathbf{H}}[k],k\in{\mathcal{K}}_{\text{side}}$
can be given as

$\displaystyle{\mathbf{h}}[k]$
$\displaystyle=\mathrm{vec}({\mathbf{H}}[k])=\sum_{\ell=1}^{L}{\alpha_{\ell}}e^{-\text{j}2\pi\Delta
f_{k}{\tau_{\ell}}}\mathrm{vec}({\bf
G}_{k}({\tau_{\ell}},{\phi_{\ell}},{\theta_{\ell}}))\odot(\overline{\breve{{\mathbf{a}}}_{\text{T},k}({\theta_{\ell}})}\otimes\breve{{\mathbf{a}}}_{\text{R},k}({\phi_{\ell}}))$
$\displaystyle\approx\sum_{l=1}^{L_{\rm
v}}\Big{[}{\bm{\Psi}}_{k}(\tau_{l}^{\rm
v},\bm{\phi}^{\text{v}},\bm{\theta}^{\text{v}})\odot\left(\overline{\breve{{\mathbf{A}}}_{\text{T},k}^{\rm
v}}\otimes\breve{{\mathbf{A}}}^{\rm v}_{\text{R},k}\right)\Big{]}\bm{b}_{l},$
(35)

where ${\bm{\Psi}}_{k}(\tau_{l}^{\rm
v},\bm{\phi}^{\text{v}},\bm{\theta}^{\text{v}})\in{\mathbb{C}}^{{N_{\text{R}}}{N_{\text{T}}}\times{K_{\text{R}}}{K_{\text{T}}}}$
and its columns are defined as

$\displaystyle\big{[}{\bm{\Psi}}_{k}(\tau_{l}^{\rm
v},\bm{\phi}^{\text{v}},\bm{\theta}^{\text{v}})\big{]}_{:,(j-1){K_{\text{T}}}+i}=e^{-\text{j}2\pi\Delta
f_{k}\tau_{l}^{\rm v}}\mathrm{vec}({\bf G}_{k}(\tau_{l}^{\rm
v},\phi^{\text{v}}_{i},\theta^{\text{v}}_{j})),\ \forall
i=1,\ldots,{K_{\text{R}}},j=1,\ldots,{K_{\text{T}}},$ (36)

and $\bm{b}_{l}\in{\mathbb{C}}^{{K_{\text{R}}}{K_{\text{T}}}\times 1}$
collects $b_{l,i,j}$, satisfying $\|\bm{b}_{l}\|_{0}\leq 1$ and
$\sum_{l=1}^{L_{\rm v}}\|\bm{b}_{l}\|_{0}\leq L$.

######  Remark 2.

Comparing the channel approximation expressions for central subcarriers in
(III-A) and for side subcarriers in (III-A), it is clear that the additional
distortion matrix
${\mathbf{G}}_{k}({\tau_{\ell}},{\phi_{\ell}},{\theta_{\ell}})$ for side
subcarriers entangles the three parameters of delay, AoA and AoD, such that a
combination of discretized delay, AoA and AoD grids is needed. This will
increase the sparsifying dictionary dimension significantly and induce
overwhelming computational complexity during the CS recovery of these
parameters. Instead, for the central subcarriers, the discretization of delay,
AoA and AoD is decoupled without the additional distortion matrix and thus
even one-dimensional search can be done sequentially for each parameter, which
will help reduce the computational complexity to a large extent.

Next, as per the approximate channel models in (III-A)-(III-A) as well as
(III-A)-(III-A), we can formulate the final DL problem for hardware
impairments. Specifically, stacking the initial channel estimates at all
subcarriers from all locations $u\in\mathcal{I}({N_{\text{sa}}})$, the problem
formulation of DLHW can be expressed as

$\displaystyle\min_{\begin{subarray}{c}{\mathbf{D}}_{\text{T},1},{\mathbf{D}}_{\text{T},2},{\mathbf{D}}_{\text{R},1},{\mathbf{D}}_{\text{R},2},\\\
\bm{\Omega}^{(u)}[k],\bm{b}^{(u)}_{l}\end{subarray}}\
\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{cen}}}\left\|\hat{{\mathbf{h}}}^{(u)}[k]-\left(\overline{\breve{{\mathbf{A}}}_{\text{T},k}^{\rm
v}}\otimes\breve{{\mathbf{A}}}^{\rm
v}_{\text{R},k}\right)\mathrm{vec}(\bm{\Omega}^{(u)}[k])\right\|_{F}^{2}$
$\displaystyle\kern
75.0pt+\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{side}}}\Big{\|}\hat{{\mathbf{h}}}^{(u)}[k]-\sum_{l=1}^{L_{\rm
v}}\Big{[}{\bm{\Psi}}_{k}(\tau_{l}^{\rm
v},\bm{\phi}^{\text{v}},\bm{\theta}^{\text{v}})\odot\left(\overline{\breve{{\mathbf{A}}}_{\text{T},k}^{\rm
v}}\otimes\breve{{\mathbf{A}}}^{\rm
v}_{\text{R},k}\right)\Big{]}\bm{b}_{l}^{(u)}\Big{\|}^{2}_{F},$
$\displaystyle\kern 20.0pt\text{subject
to}\quad\quad\big{\|}\mathrm{vec}(\bm{\Omega}^{(u)}[k])\big{\|}_{0}\leq
L^{(u)};\quad\big{\|}\bm{b}_{l}^{(u)}\big{\|}_{0}\leq
1,\quad\sum_{l=1}^{L_{\rm v}}\big{\|}\bm{b}_{l}^{(u)}\big{\|}_{0}\leq
L^{(u)}.$ (37)

Note that the two sum terms in (III-A) are format-consistent as
$\big{[}{\bm{\Psi}}_{k}(\tau_{l}^{\rm
v},\bm{\phi}^{\text{v}},\bm{\theta}^{\text{v}})\big{]}_{:,(j-1){K_{\text{T}}}+i}=\frac{\mathcal{F}(p)[f_{k}]}{{T_{\text{s}}}}e^{-\text{j}2\pi
f_{k}\tau_{l}^{\rm v}}\bm{1}_{{N_{\text{R}}}{N_{\text{T}}}}$ for central
subcarriers $k\in{\mathcal{K}}_{\text{cen}}$. We will see that the simplified
expression for central subcarriers can help facilitate the derivations of DL
algorithms.

###  III-B Speeding up the dictionary learning with a new sparse coding
algorithm

The optimization problem in (III-A) is not jointly convex with respect to the
variables ${\mathbf{D}}_{\text{T},1}$, ${\mathbf{D}}_{\text{T},2}$,
${\mathbf{D}}_{\text{R},1}$, ${\mathbf{D}}_{\text{R},2}$,
$\bm{\Omega}^{(u)}[k]$, and $\bm{b}^{(u)}_{l}$, but it can still be solved by
the alternating optimization techniques. As in the typical DL problems [25],
the optimization of (III-A) is split into two stages: sparse coding and
dictionary update.

####  III-B1 Sparse coding stage

In this stage, we fix all dictionary parts and update the channel coefficients
$\bm{\Omega}^{(u)}[k],\forall
k\in{\mathcal{K}}_{\text{cen}},u\in\mathcal{I}({N_{\text{sa}}})$ and
$\bm{b}_{l}^{(u)},\forall
k\in{\mathcal{K}}_{\text{side}},u\in\mathcal{I}({N_{\text{sa}}})$.
Specifically, for each $u\in\mathcal{I}({N_{\text{sa}}})$, the optimization
problem of (III-A) is reduced to (omitting the superscript (u) for simplicity)

$\displaystyle\min_{\bm{\Omega}[k],\bm{b}_{l}}\quad\sum_{k\in{\mathcal{K}}_{\text{cen}}}\left\|\hat{{\mathbf{h}}}[k]-\left(\overline{\breve{{\mathbf{A}}}_{\text{T},k}^{\rm
v}}\otimes\breve{{\mathbf{A}}}^{\rm
v}_{\text{R},k}\right)\mathrm{vec}(\bm{\Omega}[k])\right\|_{F}^{2}$
$\displaystyle\kern
29.0pt+\sum_{k\in{\mathcal{K}}_{\text{side}}}\Big{\|}\hat{{\mathbf{h}}}[k]-\sum_{l=1}^{L_{\rm
v}}\Big{[}{\bm{\Psi}}_{k}(\tau_{l}^{\rm
v},\bm{\phi}^{\text{v}},\bm{\theta}^{\text{v}})\odot\left(\overline{\breve{{\mathbf{A}}}_{\text{T},k}^{\rm
v}}\otimes\breve{{\mathbf{A}}}^{\rm
v}_{\text{R},k}\right)\Big{]}\bm{b}_{l}\Big{\|}^{2}_{F},$
$\displaystyle\text{subject
to}\quad\quad\big{\|}\mathrm{vec}(\bm{\Omega}[k])\big{\|}_{0}\leq
L;\quad\big{\|}\bm{b}_{l}\big{\|}_{0}\leq 1,\quad\sum_{l=1}^{L_{\rm
v}}\big{\|}\bm{b}_{l}\big{\|}_{0}\leq L.$ (38)

Generally, this problem can be optimized by various CS techniques, such as
orthogonal matching pursuit (OMP) [26], simultaneous OMP (SOMP) [27], or
simultaneous weighted OMP (SWOMP) [2], to name a few. In this paper, we
propose a new sparse coding algorithm to solve (III-B1), named Dictionary
Adaptive OMP under beam squint (DA-OMP-BS) as shown in Algorithm 1. This
algorithm exploits the following two important properties of the channel
models under beam squint, which will enable low-complexity and high-accuracy
recovery of the channel coefficients and parameters:

  * •

First, the common sparsity support property of channel vectors across
subcarriers will still be considered but with frequency-dependent dictionaries
accounting for beam squint impacts. For SOMP and OMP in previous work, the
common sparsity support of channel coefficients between subcarriers is assumed
under the same dictionary for all subcarriers. This is an approximate result
for the case without considering beam squint effect. When the steering vector
at the center frequency is used for all subcarriers, there exist approximation
errors no matter how significant the beam squint effect is. Under beam squint
circumstances, we can exploit the true property of common sparsity support
between subcarriers but with different frequency-dependent dictionaries
accounting for the beam squint impacts. The principle behind this argument is
that the physical AoAs and AoDs associated with propagation paths are constant
and independent of subcarriers, and thereby when the channel vector at each
subcarrier is projected to the corresponding sparsifying dictionary, only a
few bins out of the ${K_{\text{R}}}$ or ${K_{\text{T}}}$ virtual angular bins
corresponding to the physical AoAs and AoDs are nonzero.

  * •

Second, as mentioned above, the additional distortion matrix at side
subcarriers tangles the three parameters of delay, AoA and AoD. Then if
traditional SOMP or OMP is directly applied for (III-B1), the overall
sparsifying dictionary will be a three-dimensional (3D) dictionary [4] and the
dimension of this 3D dictionary $\allowbreak[{\bm{\Psi}}_{k}(\tau_{1}^{\rm
v},\bm{\phi}^{\text{v}},\bm{\theta}^{\text{v}})\odot\big{(}\overline{\breve{{\mathbf{A}}}_{\text{T},k}^{\rm
v}}\otimes\breve{{\mathbf{A}}}^{\rm v}_{\text{R},k}\big{)},\ldots,\\\
{\bm{\Psi}}_{k}(\tau_{L_{\rm v}}^{\rm
v},\bm{\phi}^{\text{v}},\bm{\theta}^{\text{v}})\odot\big{(}\overline{\breve{{\mathbf{A}}}_{\text{T},k}^{\rm
v}}\otimes\breve{{\mathbf{A}}}^{\rm v}_{\text{R},k}\big{)}]$ will be
${N_{\text{R}}}{N_{\text{T}}}\times L_{\rm v}{K_{\text{R}}}{K_{\text{T}}}$,
proportional to the product of the numbers of delay/AoA/AoD grids. This will
induce overwhelming complexity. As for the newly proposed DA-OMP-BS, it will
first exploit the channel models at central subcarriers to obtain initial
estimates of the delay/AoA/AoD parameters. Without the impact of additional
distortion matrix at central subcarriers, this can be done even by one-
dimensional search over each parameter space iteratively. These initial
estimates of parameters are then used to reduce the effective dictionary
dimension at side subcarriers and compensate the distortions induced by beam
squint, so that all the subcarriers can be collected to improve the estimates
of delays/AoAs/AoDs again. In doing so, the sparse coding of (III-B1) can be
solved with much lower complexity and higher accuracy.

In line of these ideas, the DA-OMP-BS is expected to outperform SOMP and OMP
without considering beam squint. We summarize the procedure of DA-OMP-BS in
Algorithm 1. Once the sparse coding stage of (III-B1) is done, the
coefficients $\bm{\Omega}^{(u)}[k],k\in{\mathcal{K}}_{\text{cen}}$ and
$\bm{b}_{l}^{(u)},k\in{\mathcal{K}}_{\text{side}}$ as well as the estimates of
path gains, delays, AoAs, and AoDs
$\hat{\bm{\varphi}}_{\hat{l}}^{(u)}\triangleq\\{\hat{\alpha}_{\hat{l}}^{(u)},\hat{\phi}_{\hat{l}}^{(u)},\hat{\theta}_{\hat{l}}^{(u)},\hat{\tau}_{\hat{l}}^{(u)}\\}$
for $\hat{l}=1,\ldots,\hat{L}^{(u)},\forall u\in\mathcal{I}({N_{\text{sa}}})$
can be obtained, which will be used to update the hardware impairments related
dictionaries in the next subsections.

Algorithm 1 : DA-OMP-BS algorithm for sparse coding

1:procedure DA-OMP-BS (${{\mathbf{y}}}^{(u)}[k]$, ${\mathbf{\Phi}}$,
${\mathbf{D}}_{\text{T},1}$, ${\mathbf{D}}_{\text{T},2}$,
${\mathbf{D}}_{\text{R},1}$,
${\mathbf{D}}_{\text{R},2},{K_{\text{R}}},{K_{\text{T}}},L_{\rm v}$)

2: Construct $\breve{{\mathbf{A}}}_{\text{R},k}^{\rm v}$ and
$\breve{{\mathbf{A}}}_{\text{T},k}^{\rm v}$ at each subcarrier with input
dictionaries

3: Compute the whitened equivalent sensing matrix and received signals

4: ${\mathbf{L}}_{\rm
w}=\text{blkdiag}\left\\{(\breve{\mathbf{W}}_{1}^{*}\breve{\mathbf{W}}_{1})^{-1/2},\ldots,(\breve{\mathbf{W}}_{M}^{*}\breve{\mathbf{W}}_{M})^{-1/2}\right\\}$
and ${\mathbf{\Phi}}_{\rm w}={\mathbf{L}}_{\rm w}{\mathbf{\Phi}}$

5: ${\mathbf{y}}_{\text{w}}^{(u)}[k]={\mathbf{L}}_{\rm
w}{\mathbf{y}}^{(u)}[k],\ \text{for}\ k=0,\ldots,K-1,\
u\in\mathcal{I}({N_{\text{sa}}})$

6: for $u=0,\ldots,{N_{\text{sa}}}-1$ do

7: Initialize the residual vectors and set of estimated parameters

8: ${\mathbf{r}}[k]={\mathbf{\Phi}}_{\rm
w}^{*}{\mathbf{y}}_{\text{w}}^{(u)}[k],\ \forall k=0,\ldots,K-1$ and
$\hat{L}^{(u)}=0,\ \hat{\bm{\varphi}}^{(u)}_{l}=\emptyset$

9: while $\text{MSE}>\epsilon$ do:

10: Define: $O(\tau,\phi,\theta,\mathcal{K}_{\rm
set})=\frac{\big{\|}\sum_{k\in\mathcal{K}_{\rm
set}}(\overline{\breve{{\mathbf{a}}}_{\text{T},k}^{\rm
v}(\theta)}\otimes\breve{{\mathbf{a}}}^{\rm
v}_{\text{R},k}(\phi))^{*}({\mathbf{r}}[k]e^{\text{j}2\pi\Delta
f_{k}\tau})\big{\|}_{F}^{2}}{\big{\|}\sum_{k\in\mathcal{K}_{\rm
set}}(\overline{\breve{{\mathbf{a}}}_{\text{T},k}^{\rm
v}(\theta)}\otimes\breve{{\mathbf{a}}}^{\rm
v}_{\text{R},k}(\phi))^{*}({\mathbf{\Phi}}_{\rm w}^{*}e^{\text{j}2\pi\Delta
f_{k}\tau})\big{\|}_{F}^{2}}$

11: Initial parameter estimation using central subcarriers

12: delay: $\hat{\tau}=\operatorname*{argmax}_{\tau_{l}^{\rm v}}\
\big{\|}\sum_{k\in{\mathcal{K}}_{\text{cen}}}{\mathbf{r}}[k]e^{\text{j}2\pi\Delta
f_{k}\tau_{l}^{\rm v}}\big{\|}_{F}^{2}$

13: AoA: $\hat{\phi}=\operatorname*{argmax}_{\phi_{i}^{\rm v}}\
\frac{\big{\|}\sum_{k\in{\mathcal{K}}_{\text{cen}}}({\mathbf{I}}_{{N_{\text{T}}}}\otimes\breve{{\mathbf{a}}}^{\rm
v}_{\text{R},k}(\phi_{i}^{\rm v}))^{*}({\mathbf{r}}[k]e^{\text{j}2\pi\Delta
f_{k}\hat{\tau}})\big{\|}_{F}^{2}}{\big{\|}\sum_{k\in{\mathcal{K}}_{\text{cen}}}({\mathbf{I}}_{{N_{\text{T}}}}\otimes\breve{{\mathbf{a}}}^{\rm
v}_{\text{R},k}(\phi_{i}^{\rm v}))^{*}({\mathbf{\Phi}}_{\rm
w}^{*}e^{\text{j}2\pi\Delta f_{k}\hat{\tau}})\big{\|}_{F}^{2}}$

14: AoD: $\hat{\theta}=\operatorname*{argmax}_{\theta_{j}^{\rm v}}\
O(\hat{\tau},\hat{\phi},\theta_{j}^{\rm v},\mathcal{K}_{\rm cen})$

15: Refine the initial parameter estimates

16: delay: $\hat{\tau}=\operatorname*{argmax}_{\tau_{l}^{\rm v}}\
O(\tau_{l}^{\rm v},\hat{\phi},\hat{\theta},\mathcal{K}_{\rm cen})$

17: AoA: $\hat{\phi}=\operatorname*{argmax}_{\phi_{i}^{\rm v}}\
O(\hat{\tau},\phi_{i}^{\rm v},\hat{\theta},\mathcal{K}_{\rm cen})$

18: AoD: $\hat{\theta}=\operatorname*{argmax}_{\theta_{j}^{\rm v}}\
O(\hat{\tau},\hat{\phi},\theta_{j}^{\rm v},\mathcal{K}_{\rm cen})$

19: Distortion compensation at side subcarriers

20:
${\mathbf{r}}[k]=\text{diag}\big{\\{}\mathrm{vec}({\mathbf{G}}_{k}(\hat{\tau},\hat{\phi},\hat{\theta}))\big{\\}}^{-1}{\mathbf{r}}[k],\
\forall k\in{\mathcal{K}}_{\text{side}}$

21: Update parameter estimates using all subcarriers

22: delay: $\hat{\tau}=\operatorname*{argmax}_{\tau_{l}^{\rm v}}\
O(\tau_{l}^{\rm v},\hat{\phi},\hat{\theta},\mathcal{K}_{\rm
cen}\cup\mathcal{K}_{\rm side})$

23: AoA: $\hat{\phi}=\operatorname*{argmax}_{\phi_{i}^{\rm v}}\
O(\hat{\tau},\phi_{i}^{\rm v},\hat{\theta},\mathcal{K}_{\rm
cen}\cup\mathcal{K}_{\rm side})$

24: AoD: $\hat{\theta}=\operatorname*{argmax}_{\theta_{j}^{\rm v}}\
O(\hat{\tau},\hat{\phi},\theta_{j}^{\rm v},\mathcal{K}_{\rm
cen}\cup\mathcal{K}_{\rm side})$

25: Update the set of estimated parameters: $\hat{L}^{(u)}=\hat{L}^{(u)}+1,\
\hat{\bm{\varphi}}^{(u)}_{l}=\\{\hat{\tau},\hat{\phi},\hat{\theta}\\}$

Algorithm 2 : DA-OMP-BS algorithm for sparse coding (continued)

26: Update the path gains $\hat{\alpha}_{l}$ and coefficient vectors by
minimizing

27:
$\sum_{k=0}^{K-1}\big{\|}{\mathbf{y}}_{\text{w}}^{(u)}[k]{-}{\mathbf{\Phi}}_{\rm
w}\sum_{l=1}^{\hat{L}^{(u)}}\alpha_{l}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{l}}\mathrm{vec}({\bf
G}_{k}(\hat{\tau}_{l},\hat{\phi}_{l},\hat{\theta}_{l}))\odot(\overline{\breve{{\mathbf{a}}}_{\text{T},k}(\hat{\theta}_{l})}\otimes\breve{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{l}))\big{\|}_{2}^{2}$

28: Update the residual for each subcarrier

29: ${\mathbf{r}}[k]{=}{\mathbf{\Phi}}_{\rm
w}^{*}\big{(}{\mathbf{y}}_{\text{w}}^{(u)}[k]{-}{\mathbf{\Phi}}_{\rm
w}\sum_{l=1}^{\hat{L}^{(u)}}\hat{\alpha}_{l}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{l}}\mathrm{vec}({\bf
G}_{k}(\hat{\tau}_{l},\hat{\phi}_{l},\hat{\theta}_{l}))\odot(\overline{\breve{{\mathbf{a}}}_{\text{T},k}(\hat{\theta}_{l})}\otimes\breve{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{l}))\big{)}$

30: Update the current MSE:
$\text{MSE}{=}\frac{1}{M{L_{\text{R}}}K}\sum_{k=0}^{K-1}\|({\mathbf{\Phi}}^{*}_{\rm
w})^{\dagger}{\mathbf{r}}[k]\|_{2}^{2}$

31: end while

32: end for

33: Output:
$\hat{\bm{\varphi}}^{(u)}_{l}{=}\\{\hat{\alpha}_{l},\hat{\tau}_{l},\hat{\phi}_{l},\hat{\theta}_{l}\\}$,
$\mathrm{vec}(\hat{\bm{\Omega}}^{(u)}[k])$ and $\bm{b}_{l}^{(u)}$, for
$u\in\mathcal{I}({N_{\text{sa}}})$

####  III-B2 Dictionary update stage

Next, we fix the channel coefficients and path parameters in (III-A) in
preparation of updating the hardware impairments related dictionaries
${\mathbf{D}}_{\text{T},1},{\mathbf{D}}_{\text{T},2},{\mathbf{D}}_{\text{R},1},{\mathbf{D}}_{\text{R},2}$.
Note that due to the special structure in this problem that the TX and RX
dictionaries are entangled with the beam squint effect at different
subcarriers, the typical dictionary update algorithms, like the method of
optimal directions (MOD) [28] or K-SVD [25], cannot be directly used.
Therefore, we apply alternating optimization in this sub-stage as well to
subsequently update the four dictionaries.

For the update of ${\mathbf{D}}_{\text{R},1}$, the problem of (III-A) can be
reduced to

$\displaystyle\min_{{\mathbf{D}}_{\text{R},1}}\quad\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{cen}}}\bigg{\|}\hat{{\mathbf{H}}}^{(u)}[k]-{\mathbf{D}}_{\text{R},1}\underbrace{\left({\mathbf{D}}_{\text{R},2}\odot{\mathbf{A}}_{\text{R},k}^{\rm
v}\right)\bm{\Omega}^{(u)}[k]\left({\mathbf{D}}_{\text{T},2}\odot{\mathbf{A}}_{\text{T},k}^{\rm
v}\right)^{*}{\mathbf{D}}_{\text{T},1}^{*}}_{{\mathbf{X}}_{\text{R},1}^{(u)}[k]}\bigg{\|}^{2}_{F}$
$\displaystyle+\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{side}}}\bigg{\|}\hat{{\mathbf{H}}}^{(u)}[k]-\sum_{\hat{l}=1}^{\hat{L}^{(u)}}\hat{\alpha}^{(u)}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}^{(u)}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}^{(u)}_{\hat{l}},\hat{\phi}^{(u)}_{\hat{l}},\hat{\theta}^{(u)}_{\hat{l}})\odot\big{(}\breve{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}^{(u)}_{\hat{l}})\breve{{\mathbf{a}}}_{\text{T},k}^{*}(\hat{\theta}^{(u)}_{\hat{l}})\big{)}\bigg{\|}_{F}^{2}.$
(39)

which is equivalent to

$\displaystyle\min_{{\mathbf{D}}_{\text{R},1}}\quad\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{cen}}}\bigg{\|}\hat{{\mathbf{H}}}^{(u)}[k]-{\mathbf{D}}_{\text{R},1}\underbrace{\left({\mathbf{D}}_{\text{R},2}\odot{\mathbf{A}}_{\text{R},k}^{\rm
v}\right)\bm{\Omega}^{(u)}[k]\left({\mathbf{D}}_{\text{T},2}\odot{\mathbf{A}}_{\text{T},k}^{\rm
v}\right)^{*}{\mathbf{D}}_{\text{T},1}^{*}}_{{\mathbf{X}}_{\text{R},1}^{(u)}[k]}\bigg{\|}^{2}_{F}$
$\displaystyle\kern
20.0pt+\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{side}}}\bigg{\|}\hat{{\mathbf{H}}}^{(u)}[k]-\sum_{\hat{l}=1}^{\hat{L}^{(u)}}\hat{\alpha}^{(u)}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}^{(u)}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}^{(u)}_{\hat{l}},\hat{\phi}^{(u)}_{\hat{l}},\hat{\theta}^{(u)}_{\hat{l}})$
$\displaystyle\kern
100.0pt\odot\Big{(}{\mathbf{D}}_{\text{R},1}\underbrace{\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}^{(u)}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}^{(u)}_{\hat{l}})\big{)}\big{(}{\mathbf{e}}_{\text{T}}(\hat{\theta}^{(u)}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{T},k}(\hat{\theta}^{(u)}_{\hat{l}})\big{)}^{*}{\mathbf{D}}_{\text{T},1}^{*}}_{{\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k]}\Big{)}\bigg{\|}_{F}^{2},$
(40)

where ${\mathbf{X}}_{\text{R},1}^{(u)}[k]$ and
${\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k]$ are defined accordingly for ease
of expression. To update ${\mathbf{D}}_{\text{R},1}$, we need to calculate the
derivative of the objective function with respect to
${\mathbf{D}}_{\text{R},1}$, i.e., $\displaystyle\frac{\partial
J}{\partial{\mathbf{D}}_{\text{R},1}}$, which is expressed as follows (A proof
is provided in Appendix II.)

$\displaystyle\frac{\partial
J}{\partial{\mathbf{D}}_{\text{R},1}}=-\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{cen}}}\overline{(\hat{{\mathbf{H}}}^{(u)}[k]-{\mathbf{D}}_{\text{R},1}{\mathbf{X}}_{\text{R},1}^{(u)}[k])}({\mathbf{X}}_{\text{R},1}^{(u)}[k])^{T}-\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{side}}}\sum_{\hat{l}=1}^{\hat{L}^{(u)}}\hat{\alpha}^{(u)}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}^{(u)}_{\hat{l}}}$ $\displaystyle\kern
40.0pt\cdot\bigg{[}\Big{(}\overline{\hat{{\mathbf{H}}}^{(u)}[k]}-\sum_{\hat{l}=1}^{\hat{L}^{(u)}}\overline{\hat{\alpha}^{(u)}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}^{(u)}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}^{(u)}_{\hat{l}},\hat{\phi}^{(u)}_{\hat{l}},\hat{\theta}^{(u)}_{\hat{l}})\odot{\mathbf{D}}_{\text{R},1}{\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k]}\Big{)}$
$\displaystyle\kern 40.0pt\odot{\bf
G}_{k}(\hat{\tau}^{(u)}_{\hat{l}},\hat{\phi}^{(u)}_{\hat{l}},\hat{\theta}^{(u)}_{\hat{l}})\bigg{]}\cdot\big{(}{\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k]\big{)}^{T}.$
(41)

Therefore, the update of ${\mathbf{D}}_{\text{R},1}$ can be obtained by
(stochastic) gradient decent as

$\displaystyle\left[{\mathbf{D}}_{\text{R},1}\right]^{\rm
new}=\left[{\mathbf{D}}_{\text{R},1}\right]^{\rm old}-\eta\frac{\partial
J}{\partial{{\mathbf{D}}_{\text{R},1}}},$ (42)

where $\eta$ is the step-size of gradient descent and can be determined by
backtracking line search. Similarly, the update of ${\mathbf{D}}_{\text{T},1}$
can be obtained.

For the update of
${\mathbf{D}}_{\text{R},2}=f(\bm{\epsilon}_{\text{R}})\triangleq
e^{-\text{j}2\pi{f_{\text{c}}}\bm{\epsilon}_{\text{R}}\cdot\sin(\bm{\phi}^{\rm
v})^{T}/c}$, it is equivalent to updating $\bm{\epsilon}_{\text{R}}$. By
stacking all subcarriers, we have the objective of updating
$\bm{\epsilon}_{\text{R}}$ as follows

$\displaystyle\min_{\bm{\epsilon}_{\text{R}}}\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{cen}}}\bigg{\|}\underbrace{{\mathbf{D}}_{\text{R},1}^{\dagger}\hat{{\mathbf{H}}}^{(u)}[k]{\mathbf{D}}_{\text{T},1}^{*{\dagger}}}_{{\mathbf{Y}}_{\text{R},2}^{(u)}[k]}-\left({\mathbf{D}}_{\text{R},2}\odot{\mathbf{A}}_{\text{R},k}^{\rm
v}\right)\underbrace{\bm{\Omega}^{(u)}[k]\left({\mathbf{D}}_{\text{T},2}\odot{\mathbf{A}}_{\text{T},k}^{\rm
v}\right)^{*}}_{{\mathbf{X}}_{\text{R},2}^{(u)}[k]}\bigg{\|}^{2}_{F}$
$\displaystyle\kern
10.0pt+\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{side}}}\bigg{\|}\hat{{\mathbf{H}}}^{(u)}[k]-\sum_{\hat{l}=1}^{\hat{L}^{(u)}}\hat{\alpha}^{(u)}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}^{(u)}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}^{(u)}_{\hat{l}},\hat{\phi}^{(u)}_{\hat{l}},\hat{\theta}^{(u)}_{\hat{l}})$
$\displaystyle\kern
70.0pt\odot\Big{(}{\mathbf{D}}_{\text{R},1}\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}^{(u)}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}^{(u)}_{\hat{l}})\big{)}\underbrace{\big{(}{\mathbf{e}}_{\text{T}}(\hat{\theta}^{(u)}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{T},k}(\hat{\theta}^{(u)}_{\hat{l}})\big{)}^{*}{\mathbf{D}}_{\text{T},1}^{*}}_{{\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k]}\Big{)}\bigg{\|}_{F}^{2},$
(43)

where ${\mathbf{X}}_{\text{R},2}^{(u)}[k]$ and
${\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k]$ are defined accordingly for ease
of expression. To update $\bm{\epsilon}_{\text{R}}$, we need to calculate the
derivative of the objective function with respect to
$\bm{\epsilon}_{\text{R}}$, which can be expressed as (A proof is given in
Appendix III)

$\displaystyle\frac{\partial
J}{\partial\bm{\epsilon}_{\text{R}}}=\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{cen}}}2\mathcal{R}\left\\{\left\\{\frac{\partial
J_{k}^{(u)}}{\partial{\mathbf{D}}_{\text{R},2}}\odot{\mathbf{D}}_{\text{R},2}\right\\}\cdot\frac{-\text{j}2\pi{f_{\text{c}}}\sin(\bm{\phi}^{\rm
v})}{c}\right\\}$
$\displaystyle-\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{side}}}2\mathcal{R}\Bigg{\\{}\sum_{\hat{l}=1}^{\hat{L}^{(u)}}\hat{\alpha}^{(u)}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}^{(u)}_{\hat{l}}}\cdot\frac{-\text{j}2\pi{f_{\text{c}}}\sin(\hat{\phi}^{(u)}_{\hat{l}})}{c}\cdot\bigg{[}{\mathbf{D}}_{\text{R},1}^{T}$
$\displaystyle\cdot\bigg{(}\Big{(}\overline{\hat{{\mathbf{H}}}^{(u)}[k]}-\sum_{\hat{l}=1}^{\hat{L}^{(u)}}\overline{\hat{\alpha}^{(u)}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}^{(u)}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}^{(u)}_{\hat{l}},\hat{\phi}^{(u)}_{\hat{l}},\hat{\theta}^{(u)}_{\hat{l}})\odot\big{(}{\mathbf{D}}_{\text{R},1}\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}^{(u)}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}^{(u)}_{\hat{l}})\big{)}{\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k]\big{)}}\Big{)}$
$\displaystyle\odot{\bf
G}_{k}(\hat{\tau}^{(u)}_{\hat{l}},\hat{\phi}^{(u)}_{\hat{l}},\hat{\theta}^{(u)}_{\hat{l}})\bigg{)}({\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k])^{T}\bigg{]}\odot\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}^{(u)}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}^{(u)}_{\hat{l}})\big{)}\Bigg{\\}}.$
(44)

Therefore, the update of ${\mathbf{D}}_{\text{R},2}$ is given as

$\displaystyle\left[{\mathbf{D}}_{\text{R},2}\right]^{\rm
new}=f([\bm{\epsilon}_{\text{R}}]^{\rm
new})=f\left([\bm{\epsilon}_{\text{R}}]^{\rm old}-\eta\frac{\partial
J}{\partial\bm{\epsilon}_{\text{R}}}\right).$ (45)

Similarly, the update of ${\mathbf{D}}_{\text{T},2}$ can be obtained.

Algorithm 3 : Dictionary learning for hardware impairments under beam squint
(DLHWBS)

  * •

Input: Initial channel estimates $\hat{{\mathbf{h}}}^{(u)}[k],\forall
k=0,\ldots,K-1,u\in\mathcal{I}({N_{\text{sa}}})$.

  * •

Initialization: Set the dictionary matrices
${\mathbf{D}}_{\text{R},1}\in{\mathbb{C}}^{{N_{\text{R}}}\times{N_{\text{R}}}}$
and
${\mathbf{D}}_{\text{T},1}\in{\mathbb{C}}^{{N_{\text{T}}}\times{N_{\text{T}}}}$
using measurement data based on DIA algorithm [29], and set the dictionary
matrices ${\mathbf{D}}_{\text{T},2}$ and ${\mathbf{D}}_{\text{R},2}$ as all-
one matrices.

  * •

While not converge do

  * 1\. Sparse coding stage: Fixing all dictionaries, solve (III-B1) using Algorithm 1 to update channel coefficients $\bm{\Omega}^{(u)}[k]$ and $\bm{b}^{(u)}_{l}$, as well as the path parameters $\hat{\bm{\varphi}}_{\hat{l}}^{(u)}$.

  * 2\. Dictionary update stage: Fixing coefficients, update dictionaries as follows

  * While not converge do

  * Update ${\mathbf{D}}_{\text{R},1}$ using (42) and update ${\mathbf{D}}_{\text{T},1}$ similarly,

  * Update ${\mathbf{D}}_{\text{R},2}$ using (45) and update ${\mathbf{D}}_{\text{T},2}$ similarly.

  * end while

  * end while

  * •

Output: The optimal dictionaries
${\mathbf{D}}_{\text{T},1},{\mathbf{D}}_{\text{T},2},{\mathbf{D}}_{\text{R},1},{\mathbf{D}}_{\text{R},2}$.

The overall procedure of the proposed DL scheme is summarized in Algorithm 3,
and the flow diagram is represented in Fig. 3.

Figure 3: Flow diagram of the dictionary learning stage.

Once the hardware impairments related dictionaries
${\mathbf{D}}_{\text{R},1},{\mathbf{D}}_{\text{R},2},{\mathbf{D}}_{\text{T},1},{\mathbf{D}}_{\text{T},2}$
are obtained, the overall RX and TX sparsifying dictionaries
$\breve{{\mathbf{A}}}^{\rm v}_{\text{R},k}$ and $\breve{{\mathbf{A}}}^{\rm
v}_{\text{T},k}$ can be constructed at each subcarrier to incorporate the beam
squint impact as in (III-A), and then used for subsequent online compressive
channel estimation, which is expected to help reduce the training overhead
significantly.

###  III-C Convergence and complexity analysis

The iterative refinement between the sparse coding and the dictionary update
stages is a decreasing process of the objective function in (III-A), which
takes positive values. In other words, we are minimizing a function bounded by
zero and its domain is closed. Therefore, the convergence is guaranteed, and
the system will be able to learn the hardware impairments.

As for the complexity, since the DL phase can be implemented offline, the
complexity involved in the dictionary update stage does not increase the
overall complexity of the online sparse coding stage. Therefore, we compare
the computational overhead of the proposed DA-OMP-BS with those of TD-OMP [4],
in terms of complex multiplication operations for each iteration. For ease of
comparison, we set the numbers of AoA/AoD/delay grids to be the same
${K_{\text{R}}}={K_{\text{T}}}=L_{\rm v}=S$. For sparse coding with DA-OMP-BS,
the parameters can be estimated by iterative one-dimensional search over each
parameter space, and thus the complexity order per iteration is
$\mathcal{O}(({N_{\text{R}}}+{N_{\text{T}}}+{N_{\text{c}}})S{N_{\text{sa}}})$.
For TD-OMP, the joint search over the 3D parameter space induces a much higher
complexity order of
$\mathcal{O}\big{(}{N_{\text{R}}}{N_{\text{T}}}+S^{3}{N_{\text{c}}}{N_{\text{sa}}}\big{)}$.

##  IV Numerical Results

In this section, we present numerical results to demonstrate the effectiveness
of the proposed dictionary learning and channel estimation algorithms for a
realistic hybrid wideband mmWave MIMO system under both hardware impairments
and beam squint.

The TX and RX are equipped with ULAs with half-wavelength spacing, i.e.,
${d_{\text{R}}}={d_{\text{T}}}={\lambda_{\text{c}}}/2$, and
${N_{\text{T}}}=32$, ${N_{\text{R}}}=8$. Regarding the RF chains,
${L_{\text{T}}}={L_{\text{R}}}={N_{\text{s}}}=2$. The number of OFDM
subcarriers is ${N_{\text{c}}}=64$ and ${T_{\text{s}}}=1~{}ns$.

The channels are generated based on (25) with $L=6$ multipath components, as
in typical indoor scenarios. The pulse shaping function $p(t)$ is assumed to
be a raised-cosine filter with roll-off factor of $\beta=0.25$. The angle of
arrival $\phi_{l}$ and angle of departure $\theta_{l}$ are uniformly
distributed in $[-\pi,\pi]$, while the delays are uniformly distributed in
$[0,16T_{\rm s}]$. The gains are generated following a complex Gaussian
distribution with the variance adjusted to achieve a specific SNR for the
evaluation of the channel estimation strategies. For the parameters of
hardware impairments, as in [17, 19, 18], the maximum gain and phase error
variances for each antenna element are set as $5\%$ and
$20^{\circ}\pi/180^{\circ}$ respectively. Finally, the mutual coupling
coefficients among antennas are within $[0.01,0.4]$ and the antenna location
errors are assumed to be uniformly distributed between
$[-0.1{\lambda_{\text{c}}},0.1{\lambda_{\text{c}}}]$.

We evaluate the performance of our proposed wideband channel estimation
algorithm, DA-OMP-BS, in addition to two baseline algorithms described in
prior work, denoted as TD-OMP [4] and WB-ADMM [1]. The sizes for the discrete
angle and delay grids are set as ${K_{\text{T}}}=2{N_{\text{T}}}$,
${K_{\text{R}}}=2{N_{\text{R}}}$ and $L_{\rm v}=2{N_{\text{c}}}$ for all the
different algorithms evaluated in the simulations. Note, however, that the
proposed Algorithm 1 is of much lower complexity, so that finer angle and
delay grids could have been used to improve parameter estimation accuracy. The
performance of the different channel estimation strategies is evaluated in
conjunction with different strategies for building the sparsifying
dictionaries. In particular, we will consider the dictionaries obtained with
our proposed approach DLHWBS, a dictionary constructed from overcomplete ideal
array response matrices (IARM) ${\mathbf{A}}_{\text{R},k}^{\rm v}$ and
${\mathbf{A}}_{\text{T},k}^{\rm v}$ assuming no any hardware impairments, the
general SeDL algorithm proposed in [19] which learns a frequency-flat combined
dictionary for hardware impairments, and the ideal dictionary built from the
known impairments and the corresponding array response vectors.

During the DL phase, the number of OFDM symbols for training is set as $500$,
and a spreading factor $10$ is used to increase the effective SNR by 10 dBs.
For the optimization of DLHWBS, a revised version of the dictionary
initialization algorithm (DIA) in [29] is used for initialization of the
hardware impairments. Specifically, the original DIA algorithm [29] is first
applied to initialize the combined dictionary
${\mathbf{D}}_{\text{R}}^{\text{init}}\triangleq{\mathbf{D}}_{\text{R},1}({\mathbf{D}}_{\text{R},2}\odot{\mathbf{A}}_{\text{R},k}^{\text{v}})$.
Then we set the initial values of the antenna spacing errors as zeros, i.e.,
$\bm{\epsilon}_{\text{R}}={\mathbf{0}}$ and
${\mathbf{D}}_{\text{R},2}=\bm{1}_{{N_{\text{R}}}\times{K_{\text{R}}}}$, and
thus the initialization of gain, phase and coupling matrix is obtained by
${\mathbf{D}}_{\text{R},1}={\mathbf{D}}_{\text{R}}^{\text{init}}({\mathbf{A}}_{\text{R},k}^{\text{v}})^{{\dagger}}$.

First we evaluate the average computational complexity time for each channel
estimation algorithm considering 100 channel realizations. The results are
shown in Fig. 4, where it can be seen that with the selected grid sizes, fixed
for all the algorithms, our method is around 500 times faster then TD-OMP, and
around 100 times faster than WD-ADMM. In the next simulations we will show
that despite this significant reduction in complexity, our method always
outperforms WD-ADMM for any selection of the dictionary, and, depending on the
system parameters, it slightly outperforms or performs similarly to TD-OMP.

Figure 4: Average computation complexity time of each evaluated channel
estimation method.

Next we evaluate the normalized mean squared error (NMSE) for the different
combinations of channel estimation algorithms and sparsifying dictionaries, as
a function of the number of training symbols $M$ and the SNR considering 100
channel realizations. Fig. 5(a) shows the NMSE when the SNR=0dB and the number
of training symbols varies from 20 to 120. It can be observed that the NMSE
reduction provided by DA-OMP-BS when exploiting the ideal dictionary based on
known impairments varies from 1 to 1.5 dB. DA-OMP-BS outperforms any other
channel estimation algorithm independently of the considered dictionary. Fig.
5(b) shows the NMSE results as function of the SNR when the number of training
symbols $M$ is set to 60. For SNR=0dB, both DA-OMP-BS leads to an NMSE value
of -10 dB as TD-OMP but with a cost in complexity 500 times lower. Leaving
aside the ideal dictionary built from known impairments, the best performing
dictionary is DLHWBS, as it can be seen in both Fig. 5(a) and Fig. 5(b).
DLHWBS outperforms SeDL because the latter attempts to learn a frequency-flat
dictionary for all subcarriers under both hardware impairments and beam squint
effect, and assumes all the channels follow the second case of (25). In other
words, the distortions on the channel models at side subcarriers are ignored,
and thus there already exist some modeling errors when SeDL is applied. These
results confirm the effectivenes of our proposed approaches both for channel
estimation and dictionary learning.

Figure 5: Comparisons of NMSE performance: (a) as a function of the number of
training OFDM symbols and SNR set to $0{\rm dB}$; (b) as a function of the SNR
for $60$ training OFDM symbols.

In Fig. 6, we also compare the spectral efficiency (SE) performances
corresponding to the various sparsifying dictionaries and sparse coding
algorithms. As in previous work [19, 2], the SE is computed by assuming fully-
digital precoding and combining using estimates for the ${N_{\text{s}}}$
dominant left and right singular vectors of the channel estimates. To be
clear, SE is defined as

$\displaystyle\text{SE}=\frac{1}{K}\sum_{k=0}^{K-1}\sum_{n=1}^{{N_{\text{s}}}}\log\left(1+\frac{\text{SNR}}{{N_{\text{s}}}}\lambda_{n}({\mathbf{H}}_{\text{eff}}[k])^{2}\right),$

where ${\mathbf{H}}_{\text{eff}}[k]$ is the effective channel after
precoding/combining and $\lambda_{n}({\mathbf{H}}_{\text{eff}}[k])$ takes the
singular values of ${\mathbf{H}}_{\text{eff}}[k]$. For both DA-OMP-BS and
conventional TD-OMP, the SE can be significantly increased when exploiting the
dictionaries learned with DLHWBS instead of IARM dictionaries. Moreover, the
performance gap between the proposed DLHWBS algorithm and the case of ideal
hardware impairment knowledge is small.

Figure 6: Comparison of SE performance: (a) as a function of the number of
training OFDM symbols and SNR set to $0{\rm dB}$; (b) as a function of the SNR
for $60$ training OFDM symbols.

Finally, we compute the BER when considering the quadrature phase shift keying
(QPSK) modulation and minimum mean squared error (MMSE) detection. Fig. 7.
shows how better channel estimates translate into lower BER for DA-OMP-BS
combined with DLHWBS. The average BER performance gap between the proposed
DLHWBS and the ideal case of known impairments is less than 0.5 dB. This
reaffirms the effectiveness of our proposed algorithm for hybrid wideband
channel under hardware impairments and beam squint.

Figure 7: Comparisons of BER performance: (a) as a function of the number of
training OFDM symbols and SNR set to $0{\rm dB}$; (b) as a function of the SNR
for $60$ training OFDM symbols.

##  V Conclusions

In this manuscript, we derived a general channel model for MIMO systems by
explicitly considering the combined impact of hardware impairments, pulse
shaping/filtering effects and beam squint. The resulting model is an extension
of existing MIMO channel models with beam squint. Based on this general
channel model, we formulated a DL problem to obtain the sparsifying
dictionaries for channel representation that account for hardware impairments.
The effect of beam squint is considered, but it is not part of the learning
process, since it can be mathematically modeled. We also proposed a novel
compressive channel estimation algorithm under beam squint, which exploits the
structure of the beam squint at different subcarriers to facilitate the
estimation of the channel parameters at a much lower complexity. Numerical
results have demonstrated the effectiveness of the proposed DL and channel
estimation strategies and the significant gains with respect to the strategies
proposed in prior work.

## Appendix I

## Proof of equation (11)

Let the delayed filter be

$\displaystyle r(t)=p(t-\tau).$ (46)

If we have a set of measurements ${\bf
r}=[r(0),r({T_{\text{s}}}),\ldots,r((K-1){T_{\text{s}}})]^{T}$, the $k$-th
entry after DFT can be computed as

$\displaystyle\hat{r}[k]=\sum_{d=0}^{K-1}[{\mathbf{r}}]_{d}\cdot
e^{-\text{j}\frac{2\pi kd}{K}}.$ (47)

Equivalently, we can use the continuous interpretation of the DFT to represent
this as

$\displaystyle\hat{r}[k]=\mathcal{F}\Big{(}\sum_{d}\delta_{dT_{\rm
s}}(p*\delta_{\tau})\Big{)}[\Delta f_{k}],$ (48)

where $\sum_{d}\delta_{dT_{\rm
s}}=\sum_{d=-\infty}^{\infty}\delta(t-d{T_{\text{s}}})$ is the sampling
function with period ${T_{\text{s}}}$ and
$\delta_{\tau}\triangleq\delta(t-\tau)$. Moreover, $*$ denotes convolution
operation and $\mathcal{F}(\cdot)[\Delta f_{k}]$ is the Fourier transform
evaluated for the frequency difference $\Delta f_{k}$. Next we can make use of
Fourier product and convolution identities to get to

$\displaystyle\hat{r}[k]=\Big{(}\mathcal{F}\big{(}\sum_{d}\delta_{dT_{\rm
s}}\big{)}*\big{(}\mathcal{F}(p)\mathcal{F}(\delta_{\tau})\big{)}\Big{)}[\Delta
f_{k}].$ (49)

Using the Dirac comb Fourier identity formula we reach

$\displaystyle\hat{r}[k]=\Big{(}\big{(}\frac{1}{T_{\rm
s}}\sum_{i}\delta_{i/T_{\rm
s}}\big{)}*\big{(}\mathcal{F}(p)\mathcal{F}(\delta_{\tau})\big{)}\Big{)}[\Delta
f_{k}].$ (50)

Now we can express the convolution in terms of an integral

$\displaystyle\hat{r}[k]=\frac{1}{T_{\rm
s}}\int\mathcal{F}(p)[f]{e^{-\text{j}2\pi f\tau}}\sum_{i}\delta_{i/T_{\rm
s}}(\Delta f_{k}-f)df,$ (51)

and this can be solved by evaluating the Dirac comb as follows

$\displaystyle\hat{r}[k]=\sum_{i}\frac{\mathcal{F}(p)[\Delta f_{k}+i/T_{\rm
s}]}{T_{\rm s}}e^{-\text{j}2\pi(\Delta f_{k}+i/T_{\rm
s})\tau}=\Big{(}\sum_{i}\frac{\mathcal{F}(p)[\Delta f_{k}+i/T_{\rm s}]}{T_{\rm
s}}e^{-\text{j}2\pi i\tau/T_{\rm s}}\Big{)}e^{-\text{j}2\pi\Delta f_{k}\tau}.$
(52)

Note that since $\mathcal{F}(p)$ is bounded, this sum only includes a few
terms. Let us define the delay-frequency distortion as

$\displaystyle g(k,\tau)=\sum_{k^{\prime}}\frac{\mathcal{F}(p)[\Delta
f_{k}+k^{\prime}/T_{\rm s}]}{T_{\rm s}}e^{-\text{j}2\pi i\tau/T_{\rm s}},$
(53)

and then (52) can be simplified to

$\displaystyle\hat{r}[k]=g(k,\tau)e^{-\text{j}2\pi\Delta f_{k}\tau}.$ (54)

If we assume the domain of $\mathcal{F}(p)$ to be in
$[-\frac{1+\beta}{2T_{s}},+\frac{1+\beta}{2T_{s}}]$ with $\beta\in[0,1]$, then
we have the expression for $g(k,\tau)$ as

$\displaystyle
g(k,\tau)=\left\\{\begin{array}[]{cl}\frac{\mathcal{F}(p)[\Delta
f_{k}]}{T_{\rm s}}+\frac{\mathcal{F}(p)[\Delta f_{k}-1/T_{\rm s}]}{T_{\rm
s}}e^{\text{j}2\pi\tau/T_{\rm s}}&\text{if}\ \ 2\Delta f_{k}>(1-\beta)/T_{\rm
s}\\\ \frac{\mathcal{F}(p)[\Delta f_{k}]}{T_{\rm s}}&\text{if}\ \ 2|\Delta
f_{k}|\leq(1-\beta)/T_{\rm s}\\\ \frac{\mathcal{F}(p)[\Delta f_{k}]}{T_{\rm
s}}+\frac{\mathcal{F}(p)[\Delta f_{k}+1/T_{\rm s}]}{T_{\rm
s}}e^{-\text{j}2\pi\tau/T_{\rm s}}&\text{if}\ \ 2\Delta
f_{k}<-(1-\beta)/T_{\rm s}\end{array}\right.$ (58)

In the case of a raised cosine filter with parameters $T=T_{\rm s}$ and
$\beta\in[0,1]$, it is straightforward to prove that

$\displaystyle
g(k,\tau)=\left\\{\begin{array}[]{cl}\frac{1}{2}(1+e^{\text{j}\frac{2\pi\tau}{T_{\rm
s}}}+(1-e^{\text{j}\frac{2\pi\tau}{T_{\rm s}}})\cos(\frac{\pi T_{\rm
s}}{\beta}(|\Delta f_{k}|-\frac{1-\beta}{2T_{\rm s}})))&\text{if}\ \ 2\Delta
f_{k}>(1-\beta)/T_{\rm s}\\\ 1&\text{if}\ \ 2|\Delta
f_{k}|\leq(1-\beta)/T_{\rm s}\\\
\frac{1}{2}(1+e^{-\text{j}\frac{2\pi\tau}{T_{\rm
s}}}+(1-e^{-\text{j}\frac{2\pi\tau}{T_{\rm s}}})\cos(\frac{\pi T_{\rm
s}}{\beta}(|\Delta f_{k}|-\frac{1-\beta}{2T_{\rm s}})))&\text{if}\ \ 2\Delta
f_{k}<-(1-\beta)/T_{\rm s}\end{array}\right.$ (62)

## Appendix II

## Derivative of the objective function in (III-B2) with respect to
${\mathbf{D}}_{\rm R,1}$

First, we calculate the derivative of the first sum term in (III-B2) with
respect to ${\mathbf{D}}_{\rm R,1}$. For
$u\in\mathcal{I}({N_{\text{sa}}}),k\in{\mathcal{K}}_{\text{cen}}$, we have

$\displaystyle\frac{\partial
J_{k}^{(u)}}{\partial{\mathbf{D}}_{\text{R},1}}=-\overline{(\hat{{\mathbf{H}}}^{(u)}[k]-{\mathbf{D}}_{\text{R},1}{\mathbf{X}}_{\text{R},1}^{(u)}[k])}({\mathbf{X}}_{\text{R},1}^{(u)}[k])^{T}.$
(63)

Next, we calculate the derivative of the second sum term in (III-B2) with
respect to ${\mathbf{D}}_{\rm R,1}$. Recalling the chain rule, we can express
the Jacobian matrix of the second sum term with respect to ${\mathbf{D}}_{\rm
R,1}$ as, $\forall
u\in\mathcal{I}({N_{\text{sa}}}),k\in{\mathcal{K}}_{\text{side}}$,

$\displaystyle\frac{\partial
J_{k}^{(u)}}{\partial\mathrm{vec}({\mathbf{D}}_{\text{R},1})^{T}}$
$\displaystyle=\frac{\partial
J_{k}^{(u)}}{\partial\mathrm{vec}(\sum_{\hat{l}=1}^{\hat{L}}\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\odot{\mathbf{D}}_{\text{R},1}{\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k])^{T}}$
$\displaystyle\cdot\frac{\partial\mathrm{vec}(\sum_{\hat{l}=1}^{\hat{L}}\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\odot{\mathbf{D}}_{\text{R},1}{\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k])}{\partial\mathrm{vec}({\mathbf{D}}_{\text{R},1}{\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k])^{T}}\cdot\frac{\partial\mathrm{vec}\big{(}{\mathbf{D}}_{\text{R},1}{\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k]\big{)}}{\partial\mathrm{vec}({\mathbf{D}}_{\text{R},1})^{T}}$
$\displaystyle=-\mathrm{vec}\Big{(}\overline{\hat{{\mathbf{H}}}^{(u)}[k]}-\sum_{\hat{l}=1}^{\hat{L}}\overline{\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\odot{\mathbf{D}}_{\text{R},1}{\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k]}\Big{)}^{T}$
$\displaystyle\cdot\sum_{\hat{l}=1}^{\hat{L}}\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}\text{diag}\big{\\{}\mathrm{vec}({\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}}))\big{\\}}\big{(}({\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k])^{T}\otimes{\mathbf{I}}_{{N_{\text{R}}}}\big{)}.$
(64)

With the Jocobian matrix, we can expressive the derivative of the second sum
term in (III-B2) with respect to ${\mathbf{D}}_{\rm R,1}$ as follows

$\displaystyle\frac{\partial
J_{k}^{(u)}}{\partial{\mathbf{D}}_{\text{R},1}}=\text{unvec}\Big{\\{}\frac{\partial
J_{k}^{(u)}}{\partial\mathrm{vec}({\mathbf{D}}_{\text{R},1})}\Big{\\}}=\text{unvec}\Big{\\{}-\sum_{\hat{l}=1}^{\hat{L}}\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}\big{(}{\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k]\otimes{\mathbf{I}}_{{N_{\text{R}}}}\big{)}$
$\displaystyle\kern 10.0pt\cdot\text{diag}\big{\\{}\mathrm{vec}({\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}}))\big{\\}}\cdot\mathrm{vec}\Big{(}\overline{\hat{{\mathbf{H}}}^{(u)}[k]}-\sum_{\hat{l}=1}^{\hat{L}}\overline{\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\odot{\mathbf{D}}_{\text{R},1}{\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k]}\Big{)}\Big{\\}}$
$\displaystyle=-\sum_{\hat{l}=1}^{\hat{L}}\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}\Big{[}\Big{(}\overline{\hat{{\mathbf{H}}}^{(u)}[k]}-\sum_{\hat{l}=1}^{\hat{L}}\overline{\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi
f_{k}\hat{\tau}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\odot{\mathbf{D}}_{\text{R},1}{\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k]}\Big{)}\odot{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\Big{]}$
$\displaystyle\kern
40.0pt\cdot\big{(}{\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k]\big{)}^{T}.$ (65)

Then combining the derivatives in $\eqref{equ:derivative_DR1_1}$ for
$k\in{\mathcal{K}}_{\text{cen}}$ and the derivatives in
$\eqref{equ:derivative_DR1_2}$ for $k\in{\mathcal{K}}_{\text{side}}$, we can
obtain the final derivative of the objective function in (III-B2) with respect
to ${\mathbf{D}}_{\rm R,1}$, i.e.,

$\displaystyle\frac{\partial
J}{\partial{\mathbf{D}}_{\text{R},1}}=-\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{cen}}}\overline{(\hat{{\mathbf{H}}}^{(u)}[k]-{\mathbf{D}}_{\text{R},1}{\mathbf{X}}_{\text{R},1}^{(u)}[k])}({\mathbf{X}}_{\text{R},1}^{(u)}[k])^{T}-\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{side}}}\sum_{\hat{l}=1}^{\hat{L}}\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}$
$\displaystyle\cdot\Big{[}\Big{(}\overline{\hat{{\mathbf{H}}}^{(u)}[k]}-\sum_{\hat{l}=1}^{\hat{L}}\overline{\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\odot{\mathbf{D}}_{\text{R},1}{\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k]}\Big{)}\odot{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\Big{]}\big{(}{\mathbf{X}}_{\text{R},1,\hat{l}}^{(u)}[k]\big{)}^{T}.$
(66)

## Appendix III

## Derivative of the objective function in (III-B2) with respect to
$\bm{\epsilon}_{\rm R}$

First, we calculate the derivative of the first sum term in (III-B2) with
respect to $\bm{\epsilon}_{\text{R}}$. Note that the gradient of any element
of ${\mathbf{D}}_{\text{R},2}$ with respect to the antenna location error
$\epsilon_{\text{R},m}$ can be expressed as

$\displaystyle\frac{\partial[{\mathbf{D}}_{\text{R},2}]_{m,n}}{\partial\epsilon_{\text{R},m}}=\frac{\partial
e^{-\text{j}2\pi{f_{\text{c}}}\epsilon_{\text{R},m}\cdot\sin(\phi_{n}^{\rm
v})/c}}{\partial\epsilon_{\text{R},m}}=[{\mathbf{D}}_{\text{R},2}]_{m,n}\cdot\frac{-\text{j}2\pi{f_{\text{c}}}\sin(\phi_{n}^{\rm
v})}{c}.$ (67)

Therefore, we have the gradient of $J_{k}^{(u)}$ with respect to
$\bm{\epsilon}_{\text{R}}$ as

$\displaystyle\frac{\partial
J_{k}^{(u)}}{\partial\bm{\epsilon}_{\text{R}}}=2\mathcal{R}\left\\{\left\\{\frac{\partial
J_{k}^{(u)}}{\partial{\mathbf{D}}_{\text{R},2}}\odot{\mathbf{D}}_{\text{R},2}\right\\}\cdot\frac{-\text{j}2\pi{f_{\text{c}}}\sin(\bm{\phi}^{\rm
v})}{c}\right\\},$ (68)

We next calculate the derivative of the second sum term in (III-B2) with
respect to $\bm{\epsilon}_{\text{R}}$. Similar to Appendix II, using the chain
rule of Jacobian matrix, we have, $\forall
u\in\mathcal{I}({N_{\text{sa}}}),k\in{\mathcal{K}}_{\text{side}}$,

$\displaystyle\frac{\partial
J_{k}^{(u)}}{\partial\bm{\epsilon}_{\text{R}}^{T}}=2\mathcal{R}\Bigg{\\{}\frac{\partial
J_{k}^{(u)}}{\partial\mathrm{vec}\Big{(}\sum_{\hat{l}=1}^{\hat{L}}\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\odot\big{(}{\mathbf{D}}_{\text{R},1}\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{\hat{l}})\big{)}{\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k]\big{)}\Big{)}^{T}}$
$\displaystyle\cdot\frac{\partial\mathrm{vec}\Big{(}\sum_{\hat{l}=1}^{\hat{L}}\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\odot\big{(}{\mathbf{D}}_{\text{R},1}\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{\hat{l}})\big{)}{\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k]\big{)}\Big{)}}{\partial\mathrm{vec}\big{(}{\mathbf{D}}_{\text{R},1}\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{\hat{l}})\big{)}{\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k]\big{)}^{T}}$
$\displaystyle\cdot\frac{\partial\mathrm{vec}\big{(}{\mathbf{D}}_{\text{R},1}\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{\hat{l}})\big{)}{\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k]\big{)}}{\partial\mathrm{vec}\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{\hat{l}})\big{)}^{T}}\cdot\frac{\partial\mathrm{vec}\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{\hat{l}})\big{)}}{\partial\mathrm{vec}\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\big{)}^{T}}\cdot\frac{\partial\mathrm{vec}\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\big{)}}{\partial\bm{\epsilon}_{\text{R}}^{T}}\Bigg{\\}}$
$\displaystyle=2\mathcal{R}\Bigg{\\{}-\mathrm{vec}\Big{(}\overline{\hat{{\mathbf{H}}}^{(u)}[k]}-\sum_{\hat{l}=1}^{\hat{L}}\overline{\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\odot\big{(}{\mathbf{D}}_{\text{R},1}\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{\hat{l}})\big{)}{\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k]\big{)}}\Big{)}^{T}$
$\displaystyle\cdot\sum_{\hat{l}=1}^{\hat{L}}\alpha_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}\text{diag}\big{\\{}\mathrm{vec}({\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}}))\big{\\}}$
$\displaystyle\cdot\big{(}({\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k])^{T}\otimes{\mathbf{D}}_{\text{R},1}\big{)}\cdot\text{diag}\big{\\{}\mathrm{vec}({{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{\hat{l}}))\big{\\}}\cdot\text{diag}\big{\\{}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\big{\\}}\cdot\frac{-\text{j}2\pi{f_{\text{c}}}\sin(\hat{\phi}_{\hat{l}})}{c}\Bigg{\\}}.$
(69)

Therefore, we have

$\displaystyle\frac{\partial J_{k}^{(u)}}{\partial\bm{\epsilon}_{\text{R}}}$
$\displaystyle=2\mathcal{R}\Bigg{\\{}-\sum_{\hat{l}=1}^{\hat{L}}\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}\cdot\frac{-\text{j}2\pi{f_{\text{c}}}\sin(\hat{\phi}_{\hat{l}})}{c}\cdot\text{diag}\big{\\{}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{\hat{l}})\big{\\}}$
$\displaystyle\cdot\big{(}({\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k])\otimes{\mathbf{D}}_{\text{R},1}^{T}\big{)}\cdot\text{diag}\big{\\{}\mathrm{vec}({\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}}))\big{\\}}$
$\displaystyle\cdot\mathrm{vec}\Big{(}\overline{\hat{{\mathbf{H}}}^{(u)}[k]}-\sum_{\hat{l}=1}^{\hat{L}}\overline{\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\odot\big{(}{\mathbf{D}}_{\text{R},1}\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{\hat{l}})\big{)}{\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k]\big{)}}\Big{)}\Bigg{\\}}$
$\displaystyle=2\mathcal{R}\Bigg{\\{}-\sum_{\hat{l}=1}^{\hat{L}}\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}\cdot\frac{-\text{j}2\pi{f_{\text{c}}}\sin(\hat{\phi}_{\hat{l}})}{c}\cdot\bigg{[}{\mathbf{D}}_{\text{R},1}^{T}$
$\displaystyle\cdot\bigg{(}\Big{(}\overline{\hat{{\mathbf{H}}}^{(u)}[k]}-\sum_{\hat{l}=1}^{\hat{L}}\overline{\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\odot\big{(}{\mathbf{D}}_{\text{R},1}\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{\hat{l}})\big{)}{\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k]\big{)}}\Big{)}$
$\displaystyle\odot{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\bigg{)}({\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k])^{T}\bigg{]}\odot\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{\hat{l}})\big{)}\Bigg{\\}}.$
(70)

Then combining the derivatives in $\eqref{equ:epsilon_update1}$ for
$k\in{\mathcal{K}}_{\text{cen}}$ and the derivatives in
$\eqref{equ:epsilon_update2}$ for $k\in{\mathcal{K}}_{\text{side}}$, we can
obtain the final derivative of the objective function in (III-B2) with respect
to $\bm{\epsilon}_{\text{R}}$, i.e.,

$\displaystyle\frac{\partial J}{\partial\bm{\epsilon}_{\text{R}}}$
$\displaystyle=\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{cen}}}2\mathcal{R}\left\\{\left\\{\frac{\partial
J_{k}^{(u)}}{\partial{\mathbf{D}}_{\text{R},2}}\odot{\mathbf{D}}_{\text{R},2}\right\\}\cdot\frac{-\text{j}2\pi{f_{\text{c}}}\sin(\bm{\phi}^{\rm
v})}{c}\right\\}$
$\displaystyle-\sum_{u\in\mathcal{I}({N_{\text{sa}}})}\sum_{k\in{\mathcal{K}}_{\text{side}}}2\mathcal{R}\Bigg{\\{}\sum_{\hat{l}=1}^{\hat{L}}\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}\cdot\frac{-\text{j}2\pi{f_{\text{c}}}\sin(\hat{\phi}_{\hat{l}})}{c}\cdot\bigg{[}{\mathbf{D}}_{\text{R},1}^{T}$
$\displaystyle\cdot\bigg{(}\Big{(}\overline{\hat{{\mathbf{H}}}^{(u)}[k]}-\sum_{\hat{l}=1}^{\hat{L}}\overline{\hat{\alpha}_{\hat{l}}e^{-\text{j}2\pi\Delta
f_{k}\hat{\tau}_{\hat{l}}}{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\odot\big{(}{\mathbf{D}}_{\text{R},1}\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{\hat{l}})\big{)}{\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k]\big{)}}\Big{)}$
$\displaystyle\odot{\bf
G}_{k}(\hat{\tau}_{\hat{l}},\hat{\phi}_{\hat{l}},\hat{\theta}_{\hat{l}})\bigg{)}({\mathbf{X}}_{\text{R},2,\hat{l}}^{(u)}[k])^{T}\bigg{]}\odot\big{(}{\mathbf{e}}_{\text{R}}(\hat{\phi}_{\hat{l}})\odot{{\mathbf{a}}}_{\text{R},k}(\hat{\phi}_{\hat{l}})\big{)}\Bigg{\\}}.$
(71)

## References

  * [1] E. Vlachos, G. C. Alexandropoulos, and J. Thompson, “Wideband MIMO Channel Estimation for Hybrid Beamforming Millimeter Wave Systems via Random Spatial Sampling,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 13, no. 5, pp. 1136–1150, 2019.
  * [2] J. Rodríguez-Fernández, N. González-Prelcic, K. Venugopal, and R. W. Heath, “Frequency-domain compressive channel estimation for frequency-selective hybrid millimeter wave MIMO systems,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 5, pp. 2946–2960, May 2018.
  * [3] Y. Han and J. Lee, “Two-stage compressed sensing for millimeter wave channel estimation,” in _Proc. of IEEE ISIT’16_ , July 2016, pp. 860–864.
  * [4] K. Venugopal, A. Alkhateeb, N. González-Prelcic, and R. W. Heath, “Channel estimation for hybrid architecture-based wideband millimeter wave systems,” _IEEE J. Sel. Areas Commun._ , vol. 35, no. 9, pp. 1996–2009, Sept. 2017\.
  * [5] J. Lee, G. Gil, and Y. H. Lee, “Exploiting spatial sparsity for estimating channels of hybrid MIMO systems in millimeter wave communications,” in _Proc. of IEEE GLOBECOM’14_ , Dec. 2014, pp. 3326–3331.
  * [6] B. Wang, F. Gao, S. Jin, H. Lin, and G. Y. Li, “Spatial- and frequency-wideband effects in millimeter-wave massive MIMO systems,” _IEEE Trans. Signal Process._ , vol. 66, no. 13, pp. 3393–3406, July 2018\.
  * [7] J. H. Brady and A. M. Sayeed, “Wideband communication with high-dimensional arrays: New results and transceiver architectures,” in _Proc. of IEEE ICCW’15_ , June 2015, pp. 1042–1047.
  * [8] T. Chu and H. Hashemi, “True-time-delay-based multi-beam arrays,” _IEEE Trans. Microw. Theory Tech._ , vol. 61, no. 8, pp. 3072–3082, Aug. 2013.
  * [9] J. Rodriguez-Fernandez and N. González-Prelcic, “Channel estimation for frequency-selective mmWave MIMO systems with beam-squint,” in _Proc. of IEEE GLOBECOM’18_ , 2018, pp. 1–6.
  * [10] J. P. González-Coma, W. Utschick, and L. Castedo, “Hybrid LISA for wideband multiuser millimeter wave communication systems under beam squint,” _CoRR_ , vol. abs/1804.09223, 2018. [Online]. Available: http://arxiv.org/abs/1804.09223
  * [11] 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 Trans. Signal Process._ , vol. 68, pp. 49–64, 2019.
  * [12] 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 Trans. Wireless Commun._ , vol. 18, no. 7, pp. 3515–3527, 2019.
  * [13] 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 Trans. Signal Process._ , vol. 67, no. 23, pp. 5893–5908, 2019.
  * [14] J. Tan and L. Dai, “Wideband Beam Tracking in THz Massive MIMO Systems,” _IEEE Journal on Selected Areas in Communications_ , vol. 39, no. 6, pp. 1693–1710, 2021.
  * [15] ——, “Wideband channel estimation for THz massive MIMO,” _China Communications_ , vol. 18, no. 5, pp. 66–80, 2021.
  * [16] P. Schniter and A. Sayeed, “Channel estimation and precoder design for millimeter-wave communications: The sparse way,” in _2014 48th Asilomar Conference on Signals, Systems and Computers_ , 2014, pp. 273–277.
  * [17] M. Eberhardt, P. Eschlwech, and E. Biebl, “Investigations on antenna array calibration algorithms for direction-of-arrival estimation,” _Advances in Radio Science: ARS_ , vol. 14, pp. 181–190, 2016.
  * [18] Y. Ding and B. D. Rao, “Dictionary learning-based sparse channel representation and estimation for FDD massive MIMO systems,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 8, pp. 5437–5451, 2018.
  * [19] H. Xie and N. González-Prelcic, “Dictionary learning for channel estimation in hybrid frequency-selective mmWave MIMO systems,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 11, pp. 7407–7422, 2020.
  * [20] J. Palacios, N. González-Prelcic, and J. Widmer, “Managing hardware impairments in hybrid millimeter wave MIMO systems: A dictionary learning-based approach,” in _Proc. of IEEE Asilomar’19_ , 2019, pp. 168–172.
  * [21] T. Wiese, L. Weiland, and W. Utschick, “Low-rank approximations for spatial channel models,” in _Proc. of IEEE WSA’16_ , Mar. 2016, pp. 1–5.
  * [22] P. Tan and N. C. Beaulieu, “Reduced ICI in OFDM systems using the” better than” raised-cosine pulse,” _IEEE Communications Letters_ , vol. 8, no. 3, pp. 135–137, 2004.
  * [23] B. Friedlander and A. J. Weiss, “Direction finding in the presence of mutual coupling,” _IEEE Trans. Antennas Propag._ , vol. 39, no. 3, pp. 273–284, 1991.
  * [24] K. V. Stavropoulos and A. Manikas, “Array calibration in the presence of unknown sensor characteristics and mutual coupling,” in _Proc. of IEEE EUSIPCO’00_ , 2000, pp. 1–4.
  * [25] M. Aharon, M. Elad, and A. Bruckstein, “K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation,” _IEEE Trans. Signal Process._ , vol. 54, no. 11, pp. 4311–4322, Nov. 2006.
  * [26] J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” _IEEE Trans. Inf. Theory_ , vol. 53, no. 12, pp. 4655–4666, Dec. 2007.
  * [27] J. A. Tropp, A. C. Gilbert, and M. J. Strauss, “Simultaneous sparse approximation via greedy pursuit,” in _Proc. of IEEE ICASSP’05_ , vol. 5, March 2005, pp. 721–724.
  * [28] K. Engan, S. O. Aase, and J. H. Husoy, “Method of optimal directions for frame design,” in _Proc. of IEEE ICASSP’99_ , vol. 5, 1999, pp. 2443–2446.
  * [29] C. Rusu and B. Dumitrescu, “An initialization strategy for the dictionary learning problem,” in _Proc. of IEEE ICASSP’14_ , May 2014, pp. 6731–6735.