# Channel Tracking and Prediction for IRS-aided Wireless Communications

Yi Wei, Ming-Min Zhao, An Liu, and Min-Jian Zhao The authors are with the
College of Information Science and Electronic Engineering, Zhejiang
University, Hangzhou 310027, China (email: {21731133, zmmblack, anliu,
mjzhao}@zju.edu.cn).

###### Abstract

For intelligent reflecting surface (IRS)-aided wireless communications,
channel estimation is essential and usually requires excessive channel
training overhead when the number of IRS reflecting elements is large. The
acquisition of accurate channel state information (CSI) becomes more
challenging when the channel is not quasi-static due to the mobility of the
transmitter and/or receiver. In this work, we study an IRS-aided wireless
communication system with a time-varying channel model and propose an
innovative two-stage transmission protocol. In the first stage, we send pilot
symbols and track the direct/reflected channels based on the received signal,
and then data signals are transmitted. In the second stage, instead of sending
pilot symbols first, we directly predict the direct/reflected channels and all
the time slots are used for data transmission. Based on the proposed
transmission protocol, we propose a two-stage channel tracking and prediction
(2SCTP) scheme to obtain the direct and reflected channels with low channel
training overhead, which is achieved by exploiting the temporal correlation of
the time-varying channels. Specifically, we first consider a special case
where the IRS-access point (AP) channel is assumed to be static, for which a
Kalman filter (KF)-based algorithm and a long short-term memory (LSTM)-based
neural network are proposed for channel tracking and prediction, respectively.
Then, for the more general case where the IRS-AP, user-IRS and user-AP
channels are all assumed to be time-varying, we present a generalized KF
(GKF)-based channel tracking algorithm, where proper approximations are
employed to handle the underlying non-Gaussian random variables. Numerical
simulations are provided to verify the effectiveness of our proposed
transmission protocol and channel tracking/prediction algorithms as compared
to existing ones.

###### Index Terms:

Intelligent reflecting surface (IRS), channel tracking, channel prediction,
Kalman filter, deep learning.

##  I Introduction

As a promising cost-effective technology for enhancing the spectral and energy
efficiency of wireless communication systems, intelligent reflecting surface
(IRS), also known as reconfigurable intelligent surface (RIS), has drawn
significant attention in both academy and industry [1]. The IRS is composed of
a massive number of low-cost passive reflecting elements, which can be smartly
controlled to dynamically configurate the wireless communication environment
for transmission enhancement and interference suppression. Due to its passive
nature, IRS requires lower hardware cost and energy consumption as compared to
the traditional active relays. As such, IRS can be densely deployed in
wireless communication systems to flexibly reconfigurate the propagation
environment, achieving improved communication capacity and reliability [2].

Several current research activities focus on how to implement IRSs [3, 4], as
well as the potentials and challenges of IRS-aided wireless communications [5,
6]. In most of the existing works, the acquisition of channel state
information (CSI) is essential since the performance gain provided by IRS is
heavily dependent on the channel estimation accuracy. However, since IRS is
passive in general and can neither send nor receive pilot symbols, the IRS-
user and access point (AP)-IRS channels cannot be estimated separately.
Therefore, the channel estimation problem in IRS-aided wireless communication
systems is much more challenging as compared to those in conventional systems
without IRS. Fortunately, the knowledge of the cascaded user-IRS-AP channel,
also known as the _reflected channel_ , is sufficient for signal detection and
beamforming design [7]. As such, most of the existing channel estimation
related contributions in the IRS literature focused on the reflected channel
estimation problem [8, 9, 10, 11, 12, 13, 14, 15, 16]. Specifically, in [8],
the authors proposed to estimate the cascaded channel coefficients one-by-one
by switching only one IRS reflecting element on at each time. In [9], a
discrete Fourier transform (DFT) based IRS phase-shift matrix (also named as
_reflection pattern_) was proposed, where all IRS reflecting elements are
designed to be active and the estimation variance is reduced as compared to
the scheme in [8]. Moreover, the works [10] and [11] focused on the channel
estimation problem in IRS-aided multi-user MISO systems. Both of them
exploited the correlation among the reflected channels to reduce the channel
training overhead, and the latter further improved the scheme in [10] by
alleviating the negative effects caused by error propagation. In addition,
several works designed channel estimation algorithms based on some special
properties of the channels [12, 13, 14]. In particular, the works [12] and
[13] formulated the channel estimation problem as a sparse channel matrix
recovery problem using the compressive sensing (CS) technique. The work [14]
exploited the low-rank structure of the channels in massive multi-input multi-
output (MIMO) systems and formulated the reflected channel estimation problem
as a combined sparse matrix factorization problem. Furthermore, when the full
instantaneous CSI is not available, the statistical CSI can be utilized to
design advanced active and passive beamforming algorithms, which usually
incurs much less channel estimation overhead [15, 16].

In most of the aforementioned works, the quasi-static channel model is
assumed, i.e., the channels are assumed to be approximately constant over a
relatively long coherence time, such that accurate estimation of the
instantaneous CSI is possible. However, in practice, when the users are with
mobility, the channel coefficients are more likely to vary and be temporally
correlated, which can be utilized to further reduce the channel training
overhead. For such time-varying channels, the existing channel estimation
algorithms are not efficient in general and may even be inapplicable. Instead,
efficient channel tracking methods are usually required to be designed to
obtain the time-varying CSI. Channel tracking has been widely studied in the
literature for traditional wireless communication systems without IRS [17, 18,
19]. Specifically, the work [17] studied the channel tracking and equalization
problem for time-varying frequency-selective MIMO channels by approximating
the MIMO channel variation using a low-order autoregressive model and tracking
the approximated channel via a Kalman filter (KF). In orthogonal frequency
division multiplexing (OFDM) systems, when time-varying frequency selective
channels are considered, the work [18] proposed to first successively track
the delay-subspace by KF and then track the channel impulse response. The
authors in [19] considered an FDD massive MIMO system with limited scattering
around the base station (BS) and proposed a two-dimensional (2D) Markov model
to capture the 2D dynamic sparsity of massive MIMO channels. An efficient
message passing algorithm was derived to recursively track the dynamic
channel. For IRS-aided communication systems, the work [20] considered a time-
varying channel model and proposed to track the downlink direct and reflected
channel using two independent KFs. It was assumed in [20] that the reflected
channel, i.e., the cascaded user-IRS-AP channel, follows the Gauss-Markov
model. However, this assumption may not hold in practice. Although the
individual user-IRS, IRS-AP and user-AP channels can be well approximately as
Gauss-Markov models as reported in many related works, the cascaded user-IRS-
AP channel is no longer Gaussian distributed.

In this work, we advance the abovementioned works by studying the channel
tracking and prediction (CTP) problem in an IRS-aided wireless communication
system and considering a more general time-varying channel model.
Specifically, the user-IRS, IRS-AP and user-AP channels are assumed to be
independent with each other and all follow stationary steady-state Gauss-
Markov processes. A two-stage transmission protocol is proposed which is able
to reduce the channel training overhead by exploiting the temporal correlation
of the channels. Based on the proposed transmission protocol, we design an
innovative two-stage CTP (2SCTP) scheme, where a KF-based channel tracking
algorithm and a deep learning (DL)-enabled channel prediction algorithm are
integrated for performance enhancement. The main contributions of this work
are summarized as follows:

  * •

First, we propose a two-stage transmission protocol and a 2SCTP scheme.
Specifically, the first stage contains two phases, in the channel training
phase, the user sends pilot symbols and the AP tracks the channel based on the
received pilot signals, while in the data transmitting phase, data signals are
sent at the AP and the user aims to recover the data based on the tracked
channel obtained in the channel training phase. In the second stage, the AP
obtains the channel coefficients via prediction and all time slots in this
stage are allocated for data transmission. Adopting the proposed two-stage
transmission protocol is able to reduce the channel training overhead
significantly, especially in the second stage.

  * •

Second, for the special case where the IRS-AP channel is assumed to be
approximately constant, we propose a KF-based algorithm for channel tracking
and a long short-term memory (LSTM)-based neural network for channel
prediction.111Considering this special case is meaningful since the locations
of the IRS and AP are usually fixed in practice and thus the IRS-AP channel
can remain constant for a relatively long time interval. Specifically, the KF-
based channel tracking algorithm is derived by modeling the channel varying
model as a linear state equation and using historical channel information
(i.e., the covariance matrix of the previously estimated channels). It is able
to recover the high-dimensional channel vector from received pilot signals
(i.e., observations), whose dimension is much lower. In the second stage, the
LSTM-based neural network, namely the _observation (OB)-LSTM_ , is designed to
predict some imaginary observations, instead of directly predicting the
channels, and then the KF-based algorithm is employed again to predict the
channel coefficients based on these imaginary observations. Note that
predicting the imaginary observations and then employing the KF-based
algorithm is more efficient than directly predicting the channels using LSTM
since the dimensions of the channel vectors are usually much larger than the
observations, especially when the number of reflecting elements is large,
which makes directly predicting the channel vectors very difficult.

  * •

Third, we study the general case where the user-IRS, IRS-AP and user-AP
channels are all assumed to be time-varying. In this case, the KF-based
channel tracking algorithm cannot be directly applied, since not all the
random variables (i.e., the channel coefficients) involved in the state
equation are complex Gaussian distributed. To tackle this difficulty, we
propose a novel approximation method to approximate the underlying non-
Gaussian variables with Gaussian ones, based on which the state equation is
transformed into a linear and Gaussian system. Then, a generalized KF
(GKF)-based channel tracking algorithm is proposed, and by combining it with
the LSTM-based channel prediction algorithm, we show that the 2SCTP scheme is
also effective in the considered more general case.

  * •

Finally, extensive numerical simulations are presented to demonstrate the
effectiveness of the proposed 2SCTP scheme. We show that by exploiting the
temporal correlation of the channel, much lower channel training overhead can
be achieved by the proposed 2SCTP scheme as compared to the existing channel
estimation methods.

The rest of the paper is organized as follows. Section II presents the system
model and the framework of the proposed 2SCTP scheme. Then, the proposed CTP
algorithms for the special and general cases are respectively given in Section
III and Section IV. Numerical results are presented in Section V, and finally
Section VI concludes the paper.

Notations: Scalars, vectors and matrices are respectively denoted by lower
(upper) case, boldface lower case and boldface upper case letters. For a
matrix $\mathbf{X}$ of arbitrary size, $\mathbf{X}^{T}$, $\mathbf{X}^{*}$,
$\mathbf{X}^{H}$ denote the transpose, conjugate and conjugate transpose of
$\mathbf{X}$, respectively. $[\mathbf{x}]_{n}$ represents the $n$-th element
of the vector $\mathbf{x}$, and $\mathbf{X}_{n}$ denotes the $n$-th column of
the matrix $\mathbf{X}$. The symbol $||\cdot||$ denotes the Euclidean norm of
a complex vector, and $|\cdot|$ is the absolute value of a complex scalar.
$\text{diag}(x_{1},\cdots,x_{N})$ denotes a diagonal matrix whose diagonal
elements are set as $x_{1},\cdots,x_{N}$. ${\mathbb{C}}^{m\times n}$ denotes
the space of $m\times n$ complex matrices. The notations
$\mathbb{E}\big{[}\cdot\big{]}$ and $\text{Var}\big{[}\cdot\big{]}$ represent
the expectation and variance of a random variable, and
$\text{Cov}\big{[}\cdot,\cdot\big{]}$ denote the covariance of two random
variables. The remainder operation is denoted by %, i.e.,
$x_{1}\%x_{2}=x_{1}-x_{2}\lfloor\frac{x_{1}}{x_{2}}\rfloor$ with $\lfloor
x\rfloor$ denoting the maximum integer that is smaller than $x$. The symbol
$j$ is used to represent $\sqrt{-1}$. Finally, we define the complex Gaussian
distribution with mean $\mu$ and variance ${\sigma}^{2}$ as
$\mathcal{C}\mathcal{N}(\mu,{\sigma}^{2})$.

##  II System Model and Transmission Protocol

###  II-A System Model

Figure 1:  System model of the considered IRS-aided communication system.
Figure 2:  Illustration of the proposed frame structure.

As shown Fig. 1, we consider an IRS-aided time-division duplexing (TDD) system
where a single-antenna AP communicates with a single-antenna user via an IRS
with $N$ reflecting elements. The CSI is obtained via uplink pilot
transmission in the considered system with the assumption of channel
reciprocity. Let $\mathbf{g}\in\mathbb{C}^{N\times 1}$,
$\mathbf{h}_{r}\in\mathbb{C}^{N\times 1}$ and $h_{d}$ denote the baseband
equivalent IRS-AP, user-IRS and user-AP channels, respectively. Assume that
the channel statistics will stay unchanged for a long period of time called a
super frame, which includes a number of frames. Each frame is divided into $T$
time intervals and each time interval further consists of $\tau$ time slots,
as illustrated in Fig. 2. It is assumed that all the channels remain
approximately constant in each time interval and vary from the current time
interval to the next. Due to the insufficient angular spread of the scattering
environment and closely spaced antennas/reflecting elements, both line-of-
sight (LoS) and non-LoS (NLoS) components may exist in practical channels. As
a result, the IRS-AP, user-IRS and user-AP channels in the $t$-th time
interval can be respectively modeled as follows:

$\mathbf{g}(t)=\sqrt{\frac{l_{\text{IA}}\beta_{\text{IA}}}{1+\beta_{\text{IA}}}}\mathbf{g}^{\text{LoS}}+\sqrt{\frac{l_{\text{IA}}}{1+\beta_{\text{IA}}}}\mathbf{g}^{\text{NLoS}},$
(1)
$\mathbf{h}_{r}(t)=\sqrt{\frac{l_{\text{UI}}\beta_{\text{UI}}}{1+\beta_{\text{UI}}}}\mathbf{h}_{r}^{\text{LoS}}+\sqrt{\frac{l_{\text{UI}}}{1+\beta_{\text{UI}}}}\mathbf{h}_{r}^{\text{NLoS}},$
(2)
${h}_{d}(t)=\sqrt{\frac{l_{\text{UA}}\beta_{\text{UA}}}{1+\beta_{\text{UA}}}}{h}_{d}^{\text{LoS}}+\sqrt{\frac{l_{\text{UA}}}{1+\beta_{\text{UA}}}}{h}_{d}^{\text{NLoS}},$
(3)

where $\mathbf{g}^{\text{LoS}}$, $\mathbf{h}_{r}^{\text{LoS}}$ and
${h}_{d}^{\text{LoS}}$ represent the LoS components;
$\mathbf{g}^{\text{NLoS}}$, $\mathbf{h}_{r}^{\text{NLoS}}$ and
${h}_{d}^{\text{NLoS}}$ denote the NLoS components; $\beta_{\text{IA}}$,
$\beta_{\text{UI}}$ and $\beta_{\text{UA}}$ are the Rician factors of the IRS-
AP, user-IRS and user-AP channels, respectively. $l_{\text{IA}}$,
$l_{\text{UI}}$ and $l_{\text{UA}}$ are the corresponding path losses, which
are given by
$l_{\text{IA}}=l_{0}\left({d_{\text{IA}}}/{d_{0}}\right)^{-\gamma^{\text{IA}}}$,
$l_{\text{UI}}=l_{0}\left({d_{\text{UI}}}/{d_{0}}\right)^{-\gamma^{\text{UI}}}$
and
$l_{\text{UA}}=l_{0}\left({d_{\text{UA}}}/{d_{0}}\right)^{-\gamma^{\text{UA}}}$,
respectively, where $d_{0}$ represents the reference distance and $l_{0}$ is
the path loss at the reference distance, ${d_{\text{IA}}}$, ${d_{\text{UI}}}$
and ${d_{\text{UA}}}$ denote the link distances from the IRS to the AP, from
the user to the IRS, and from the user to the AP, respectively;
$\gamma^{\text{IA}}$, $\gamma^{\text{UI}}$ and $\gamma^{\text{UA}}$ denote the
path-loss exponents. Equivalently, we have
$\mathbf{g}(t)\sim\mathcal{CN}(\bar{\mathbf{g}},\mathbf{C}_{\text{IA}})$,
$\mathbf{h}_{r}(t)\sim\mathcal{CN}(\bar{\mathbf{h}}_{r},\mathbf{C}_{\text{UI}})$
and ${h}_{d}(t)\sim\mathcal{CN}(\bar{{h}}_{d},{C}_{\text{UA}})$, where
$\mathbf{C}_{\text{IA}}=\frac{l_{\text{IA}}}{1+\beta_{\text{IA}}}\mathbb{E}\big{[}\mathbf{g}^{\text{NLoS}}(\mathbf{g}^{\text{NLoS}})^{H}\big{]}$,
$\mathbf{C}_{\text{UI}}=\frac{l_{\text{UI}}}{1+\beta_{\text{UI}}}\mathbb{E}\big{[}\mathbf{h}_{r}^{\text{NLoS}}(\mathbf{h}_{r}^{\text{NLoS}})^{H}\big{]}$
and
${C}_{\text{UA}}=\frac{l_{\text{UA}}}{1+\beta_{\text{UA}}}\mathbb{E}\big{[}{h}_{d}^{\text{NLoS}}({h}_{d}^{\text{NLoS}})^{*}\big{]}$.

Due to the mobility of the AP and/or user, the channels are usually time
varying and exhibit a high degree of temporal correlation[21], [22].
Therefore, we employ the independent stationary steady-state Gauss-Markov
process [23], [24] to model the temporal evolution of the channel parameters,
which are given as follows:

$\mathbf{g}(t)=\sqrt{1-\alpha_{\text{IA}}}\big{(}\mathbf{g}(t-1)-\bar{\mathbf{g}}\big{)}+\sqrt{\alpha_{\text{IA}}}\mathbf{u}_{\textrm{IA}}(t)+\bar{\mathbf{g}},$
(4a)
$\mathbf{h}_{r}(t)=\sqrt{1-\alpha_{\text{UI}}}\big{(}\mathbf{h}_{r}(t-1)-\bar{\mathbf{h}}_{r}\big{)}+\sqrt{\alpha_{\text{UI}}}\mathbf{u}_{\textrm{UI}}(t)+\bar{\mathbf{h}}_{r},$
(4b)
${h}_{d}(t)=\sqrt{1-\alpha_{\text{UA}}}\big{(}{h}_{d}(t-1)-\bar{h}_{d}\big{)}+\sqrt{\alpha_{\text{UA}}}{u}_{\textrm{UA}}(t)+\bar{h}_{d},$
(4c)

where
$\mathbf{u}_{\text{IA}}(t)\sim\mathcal{CN}(\mathbf{0},\sigma^{2}_{\text{IA}}\mathbf{I})$,
$\mathbf{u}_{\text{UI}}(t)\sim\mathcal{CN}(\mathbf{0},\sigma^{2}_{\text{UI}}\mathbf{I})$
and
${u}_{\text{UA}}(t)\sim\mathcal{CN}({0},{\sigma}^{2}_{\text{UA}}\mathbf{I})$
are the perturbation terms in the IRS-AP, user-IRS and user-AP channels,
respectively; $\alpha_{\text{IA}}$, $\alpha_{\text{UI}}$ and
$\alpha_{\text{UA}}$ represent the temporal correlation coefficients of the
IRS-AP, user-IRS and user-AP channels, respectively. Note that the channel
realizations $\mathbf{g}(t-1)$, $\mathbf{h}_{r}(t-1)$ and ${h}_{d}(t-1)$ are
statistically independent of $\mathbf{u}_{\text{IA}}(t)$,
$\mathbf{u}_{\text{UI}}(t)$ and $u_{\text{UA}}(t)$.

In the $t$-th time interval, the received signal at the AP can be expressed as

$y(t)=\sqrt{p}[\mathbf{h}_{r}^{H}(t)\bm{\Theta}(t)\mathbf{g}(t)+h_{d}(t)]s(t)+z(t),$
(5)

where $s(t)$ denotes the transmit symbol, $p$ represents the transmit power,
$z(t)$ denotes the additive white Gaussian noise (AWGN) with zero-mean and
variance $\sigma^{2}$, $\bm{\Theta}(t)$ represents the reflection pattern at
the IRS. Let us define
$\mathbf{h}\triangleq[h_{d};\text{diag}(\mathbf{h}_{r}^{H})\mathbf{g}]$ as the
user-AP equivalent channel, then the received signal in (5) can be
equivalently rewritten as

$y(t)=\sqrt{p}\mathbf{v}^{H}(t)\mathbf{h}(t)s(t)+z(t),$ (6)

where $\mathbf{v}(t)\triangleq[1,\theta_{1}(t),\cdots,\theta_{N}(t)]^{H}$ with
$\theta_{n}(t)$ denoting the $n$-th diagonal element of the reflection pattern
$\bm{\Theta}(t)$. Notice that when the number of IRS reflecting elements is
large, the acquisition of the instantaneous CSI may cause considerable channel
training/estimation overhead, which leads to reduced user transmission rate
due to the limited time left for data transmission. To address this issue, we
propose to exploit the temporal correlation of the channel in this work and
present a 2SCTP scheme, where a KF-based channel tracking algorithm and a DL-
enabled channel prediction algorithm are integrated to reduce the channel
training overhead and thus improve the transmission rate.

###  II-B Transmission Protocol

The proposed two-stage transmission protocol is shown in Fig. 2. As can be
seen, in the first stage that contains $T_{1}$ time intervals, each time
interval is divided into two phases, i.e., the channel training phase
(consists of $\tau_{1}$ time slots) and the data transmission phase (consists
of $\tau_{2}=\tau-\tau_{1}$ time slots). More specifically, in the first
$\tau_{1}$ time slots of the first stage, the user sends pilot symbols that
are known at the AP such that the AP can track the channel according to the
received pilot signals, then the remaining $\tau-\tau_{1}$ time slots of each
time interval are used for data transmission from the AP to the user. In the
second stage which includes $T_{2}=T-T_{1}$ time intervals, there is no
channel training phase and the channels are predicted such that all the time
slots in these time intervals are utilized for data transmission.

As compared with the existing channel estimation methods, the proposed 2SCTP
scheme is able to reduce the channel training overhead in the following two
ways:

  * •

In the first stage, we exploit the temporal correlation of the channels to
reduce the channel training overhead in each time interval, thus $\tau_{1}$
can be much smaller than $N+1$.

  * •

In the second stage, we predict the channels based on the statistical channel
information extracted from the CSI collected in the first stage, hence the
channel training overhead can be completely removed.

Based on the proposed transmission protocol, we propose the 2SCTP scheme,
whose details are given in the following two sections. Note that at the
beginning of each super frame, we need to rerun the proposed 2SCTP scheme to
track the channel since the channel statistics change. Therefore, in the rest
of the paper, we focus on the CTP algorithm design within one frame. Besides,
the temporal correlation coefficients $\alpha_{\text{IA}}$,
$\alpha_{\text{UI}}$ and $\alpha_{\text{UA}}$, the variances of the
perturbation terms $\sigma^{2}_{\text{IA}}$, $\sigma^{2}_{\text{UI}}$ and
${\sigma}^{2}_{\text{UA}}$, and the channel means $\bar{\mathbf{g}}$,
$\bar{\mathbf{h}}_{r}$ and $\bar{h}_{d}$ are regarded as prior knowledge of
the proposed 2SCTP scheme.222 These parameters can be obtained by efficient
parameter estimation methods, e.g., the expectation-maximization (EM)
algorithm [25], and further investigation is left for future work.

##  III 2SCTP Scheme for the Special Case

In this section, we focus on the special case where the IRS-AP channel is
assumed to change much slower than the user-IRS and user-AP channels, i.e.,
$\alpha_{\text{IA}}\ll\alpha_{\text{UI}}$ and
$\alpha_{\text{IA}}\ll\alpha_{\text{UA}}$. Furthermore, for simplicity, we
assume that the user-AP, IRS-AP and user-IRS channels follow the Rayleigh
channel model, extension to the more general Rician channel model will be
discussed later. Based on these assumptions, the user-AP equivalent channel
$\mathbf{h}(t)$ can be modeled as

$\mathbf{h}(t)=\mathbf{A}_{h}\mathbf{h}(t-1)+\mathbf{B}_{h}\mathbf{u}(t),$ (7)

where

$\mathbf{A}_{h}=\left[\begin{array}[]{cc}\sqrt{1-\alpha_{\text{UA}}}&\mathbf{0}_{1\times
N}\\\ \mathbf{0}_{N\times
1}&\sqrt{1-\alpha_{\text{UI}}}\mathbf{I}\end{array}\right],\mathbf{B}_{h}=\left[\begin{array}[]{cc}\sqrt{\alpha_{\text{UA}}}&\mathbf{0}_{1\times
N}\\\ \mathbf{0}_{N\times
1}&\sqrt{\alpha_{\text{UI}}}\mathbf{I}\end{array}\right],$ (8)

$\mathbf{u}(t)\sim\mathcal{CN}(\mathbf{0},\mathbf{C}_{h})$ and
$\mathbf{C}_{h}$ is given by

$\mathbf{C}_{h}=\left[\begin{array}[]{cc}\sigma^{2}_{\text{UA}}&\mathbf{0}_{1\times
N}\\\ \mathbf{0}_{N\times
1}&\sigma^{2}_{\text{UI}}\text{diag}(\mathbf{g}\mathbf{g}^{H})\end{array}\right].$
(9)

In this special case, although the reflected channel
$\mathbf{g}^{H}(t)\text{diag}\big{(}\mathbf{h}_{r}(t)\big{)}$ is related with
the IRS-AP and user-IRS channels, it can be assumed to change according to one
temporal correlation coefficient $\alpha_{\text{UI}}$ since
$\alpha_{\text{IA}}\ll\alpha_{\text{UI}}$. At one extreme, if
$\alpha_{\text{UA}}=\alpha_{\text{UI}}=0$, the reflected channel parameters
are totally correlated, (i.e., $\mathbf{h}(t)=\mathbf{h}(t-1)$), while at the
other extreme when $\alpha_{\text{UA}}=\alpha_{\text{UI}}=1$, the reflected
channel coefficients evolve according to an uncorrelated Gaussian random
process.

###  III-A KF-based Channel Tracking in Both Stages

As a well-known and powerful variable estimation method, the KF method has
been widely applied in time series analysis, such as signal processing and
econometrics [26, 27]. The KF method keeps track of the estimated state of the
system and the uncertainty of the estimate, and only the estimated state from
the previous time step and the current measurement are required to obtain the
current state estimate and the corresponding uncertainty. Inspired by the
superiority of the KF method in time series processing, we propose a KF-based
channel tracking algorithm to estimate the time-varying channels for the
considered special case with limited channel training overhead. Notice that
the proposed channel tracking algorithm is employed in both stages.

First of all, let us introduce the problem setting for a general KF method.
The KF model assumes that the state of a system at time step $t$ evolved from
the prior state at time step $t-1$ according to the linear state equation
${x}_{t}={F}_{t}{x}_{t-1}+{w}_{t},$ where $F_{t}$ is the state transition
model, the state variable $x_{t}$ and state noise variable $w_{t}$ follow
independent complex Gaussian distributions. Then, the observation at time step
$t$ is obtained according to the linear observation equation
${z}_{t}={H}_{t}{x}_{t-1}+{r}_{t},$ where $H_{t}$ is the measurement matrix,
the observation variable $z_{t}$ and observation noise variable $r_{t}$ also
follow complex Gaussian distributions. The KF method aims to predict the state
variable by utilizing a series of observations.

To apply the KF method to the considered channel tracking problem, we can
regard (7) as the state equation. By assuming that the pilot symbol is set as
$s=1$ without loss of generality, the observation functions in the first and
second stages are respectively given by

$\mathbf{y}(t)=\sqrt{p}\mathbf{V}(t)\mathbf{h}(t)+\mathbf{z}(t),\;\;\;t=1,\cdots,T_{1},$
(10)
${\mathbf{y}}_{i}(t)=\sqrt{p}\mathbf{V}(t)\mathbf{h}(t)+\mathbf{z}(t),\;\;\;t=T_{1}+1,\cdots,T,$
(11)

where the received signal $\mathbf{y}(t)\in\mathbb{C}^{\tau_{1}\times 1}$ and
${\mathbf{y}}_{i}(t)\in\mathbb{C}^{\tau_{1}\times 1}$ serve as the real and
imaginary observations, respectively,
$\mathbf{V}(t)\triangleq[\mathbf{v}_{1}(t),\cdots,\mathbf{v}_{\tau_{1}}(t)]^{T}\in\mathbb{C}^{\tau_{1}\times(N+1)}$
is regarded as the measurement matrix, and $\mathbf{v}_{i}(t)$ represents the
reflection pattern in the $i$-th time slot of the $t$-th time interval. Note
that the _imaginary observation_ is defined as the received signal if we
imagine that the user sends the pilot symbol $1$ with the known reflection
patterns for $t=T_{1}+1,...,T$, and the details of how to obtain (predict) the
imaginary observations in $t=T_{1}+1,...,T$ from the real observations in
$t=1,...,T_{1}$ using the LSTM-based neural network will be elaborated in
Section III-B. Furthermore, the value of $\tau_{1}$, i.e., the dimension of
the received signal $\mathbf{y}(t)$, will affect the convergence of the
proposed KF-based channel tracking algorithm and the required channel training
overhead, thus it should be carefully selected. Besides, to facilitate the
prediction of the imaginary observations using the LSTM-based neural network
and improve the channel tracking performance, we employ a periodic measurement
matrix which satisfies

$\mathbf{V}(t)=[\mathbf{Q}_{(\tau_{1}(t-1)+1)\%(N+1)},\cdots,\mathbf{Q}_{(\tau_{1}t)\%(N+1)}]^{T},$
(12)

where each column is chosen from a full-rank reference matrix
$\mathbf{Q}\in\mathbb{C}^{(N+1)\times(N+1)}$. The design of the measurement
matrix will be discussed in Section III-B, and the effect of choosing
different $\tau_{1}$ and employing different reference matrices will be
investigated in Section V. For clarity, we introduce

$\tilde{\mathbf{y}}(t)=\left\\{\begin{array}[]{ll}\mathbf{y}(t),&t=1,\cdots,T_{1},\\\
\mathbf{y}_{i}(t),&t=T_{1}+1,\cdots,T.\end{array}\right.$ (13)

Let us define
$\mathbf{y}_{t}\triangleq[\tilde{\mathbf{y}}^{T}(t),\tilde{\mathbf{y}}^{T}(t-1),\cdots,\tilde{\mathbf{y}}^{T}(1)]^{T}$
as the vector of observations and define
$\hat{\mathbf{y}}(t)\triangleq\tilde{\mathbf{y}}(t)-\mathbf{V}(t)\hat{\mathbf{h}}(t)$
as the innovation with
$\hat{\mathbf{h}}(t)=\mathbb{E}\big{[}\mathbf{h}(t)|\mathbf{y}_{t}\big{]}$
denoting the predicted channel vector. Then, the minimum mean square error
(MMSE) estimate of $\mathbf{h}(t)$ can be recursively obtained by the KF
equations [28], where each iteration works in the following two-step process.

  * •

In the prediction step, the estimate of the current state variable
$\hat{\mathbf{h}}(t)$ along with its covariance matrix $\mathbf{M}(t)$ are
predicted using the state estimated from the previous time interval, i.e.,

$\hat{\mathbf{h}}(t)=\mathbf{A}_{h}\mathbf{h}_{\text{KF}}(t-1),$ (14)
$\mathbf{M}(t)=\mathbb{E}\big{[}\big{(}\mathbf{h}(t)-\hat{\mathbf{h}}(t)\big{)}\big{(}\mathbf{h}(t)-\hat{\mathbf{h}}(t)\big{)}^{H}\big{]}=\mathbf{A}_{h}\mathbf{M}_{\text{KF}}(t-1)\mathbf{A}_{h}^{T}+\mathbf{B}_{h}\mathbf{C}_{h}\mathbf{B}_{h}^{T},$
(15)

where $\mathbf{h}_{\text{KF}}(t)$ represents the refined state estimate in the
$t$-th time interval (also named as _correction_), which will be introduced
later, and
$\mathbf{M}_{\text{KF}}(t)=\mathbb{E}\big{[}\big{(}\mathbf{h}(t)-\mathbf{h}_{\text{KF}}(t-1)\big{)}\big{(}\mathbf{h}(t)-\mathbf{h}_{\text{KF}}(t-1)\big{)}^{H}\big{]}$
denotes the estimation covariance matrix.

  * •

In the update step, we first update the Kalman gain
$\mathbf{G}_{\text{KF}}(t)$ as

$\mathbf{G}_{\text{KF}}(t)=\mathbf{M}(t)\mathbf{V}^{H}(t)\big{(}\mathbf{V}(t)\mathbf{M}(t)\mathbf{V}^{H}(t)+\sigma^{2}\mathbf{I}\big{)}^{-1}.$
(16)

Then, multiplying the innovation $\hat{\mathbf{y}}(t)$ by the Kalman gain
$\mathbf{G}_{\text{KF}}(t)$ and combining $\hat{\mathbf{y}}(t)$ with the state
estimate $\hat{\mathbf{h}}(t)$, the correction $\mathbf{h}_{\text{KF}}(t)$ can
be obtained via

$\mathbf{h}_{\text{KF}}(t)=\mathbb{E}\big{[}\mathbf{h}(t)|\mathbf{y}_{t-1}\big{]}+\mathbb{E}\big{[}\mathbf{h}(t)|\hat{\mathbf{y}}(t)\big{]}=\hat{\mathbf{h}}(t)+\mathbf{G}_{\text{KF}}(t)\hat{\mathbf{y}}(t).$
(17)

Finally, the estimation covariance matrix $\mathbf{M}_{\text{KF}}(t)$ required
to calculate the prediction covariance matrix $\mathbf{M}(t)$ in (16), is
updated by

$\mathbf{M}_{\text{KF}}(t)=\mathbb{E}\big{[}\big{(}\mathbf{h}(t)-\mathbf{h}_{\text{KF}}(t-1)\big{)}\big{(}\mathbf{h}(t)-\mathbf{h}_{\text{KF}}(t-1)\big{)}^{H}\big{]}=\mathbf{M}(t)-\mathbf{G}_{\text{KF}}(t)\mathbf{V}^{H}(t)\mathbf{M}(t).$
(18)

To summarize, the proposed KF-based channel tracking algorithm is shown in
Algorithm 1, where we set $\mathbf{h}_{\text{KF}}(-1)=\mathbf{0}$ and
$\mathbf{M}_{\text{MF}}=\mathbb{E}\big{[}\mathbf{h}\mathbf{h}^{H}\big{]}$ to
initialize the KF iterations. Besides, Algorithm 1 can be readily extended to
the Rician channel model. Specifically, let $\bar{\mathbf{h}}$ represent the
mean of $\mathbf{h}(t)$, and by substituting
$\tilde{\mathbf{h}}(t)=\mathbf{h}(t)-\bar{\mathbf{h}}$ into (7), we obtain

$\tilde{\mathbf{h}}(t)=\mathbf{A}_{h}\tilde{\mathbf{h}}(t-1)+\mathbf{B}_{h}\mathbf{u}(t).$
(19)

Regarding (19) as the state equation, the estimation of the state variable in
the $t$-th time interval, i.e., $\tilde{\mathbf{h}}_{\text{KF}}(t)$, can be
obtained by directly applying Algorithm 1, after which we can acquire the
refined state estimate in the $t$-th time interval as
$\mathbf{h}_{\text{KF}}(t)=\tilde{\mathbf{h}}_{\text{KF}}(t)+\bar{\mathbf{h}}$.

Algorithm 1 KF-based Channel Tracking Algorithm

Input: covariance matrix
$\mathbf{C}_{h}=\mathbb{E}\big{[}\mathbf{h}\mathbf{h}^{H}\big{]}$ and
observations $\tilde{\mathbf{y}}(t)$
Output: Estimated channel vectors $\mathbf{h}_{\text{KF}}(t)$

1:Initialization: $t=0$, $\mathbf{h}_{\text{KF}}(-1)=\mathbf{0}$,
$\mathbf{M}_{\text{MF}}=\mathbb{E}\big{[}\mathbf{h}\mathbf{h}^{H}\big{]}$

2:while $t>=1$ do

3:  _Prediction step_ : Calculate the prediction $\hat{\mathbf{h}}(t)$
according to (14), Calculate the prediction covariance matrix $\mathbf{M}(t)$
according to (15),

4:  _Update step_ : Update the Kalman gain vector $\mathbf{G}_{\text{KF}}(t)$
according to (16), Calculate the correction $\mathbf{h}_{\text{KF}}(t)$
according to (17), Calculate the estimation covariance matrix
$\mathbf{M}_{\text{KF}}(t)$ according to (18),

5: $t=t+1$,

6:end while

###  III-B DL-enabled Channel Prediction in the Second Stage

Solving the channel prediction problem can be interpreted as a time series
processing task. Recently, DL techniques designed for time series processing
have been widely used to solve channel prediction problems, e.g., using
recurrent neural network (RNN) for frequency-domain channel prediction [29],
deploying LSTM networks to predict fading channels [30], etc. However,
predicting a high-dimensional channel vector in massive MIMO systems or IRS-
aided communication systems is not easy, and it is challenging to achieve high
prediction accuracy. Moreover, in order to process the high-dimensional data,
the prediction neural network requires a huge number of neurons/layers to
ensure adequate network capacity. This will lead to excessive computing
resources and memory costs. Assume that the varying reflection patterns
$\mathbf{v}(t),t=1,\cdots,T$, are known at the AP and user, we can obtain the
imaginary observation $\mathbf{y}_{i}(t)$ in the $t$-th time interval
according to (11). Then, to tackle the abovementioned challenges, we propose
to design a novel OB-LSTM network to predict the low-dimensional imaginary
observations, instead of the high-dimensional channel vector. These predicted
imaginary observations are then regarded as the input of the KF-based channel
tracking algorithm, which in turn outputs the current channel vector. Based on
the above strategy, no pilot symbols are required to be transmitted in the
second stage, thus the channel training overhead is zero in this stage and the
transmission rate of the user is maximized.

The detailed structure of the proposed OB-LSTM network is provided as follows.
First, let $L_{I}$ and $L_{P}$ denote the lengths of the input observations
and predicted imaginary observations, respectively, where $L_{I}\leq T_{1}$
and $L_{P}\leq T_{2}$. Then, in the $t$-th time interval for
$t=T_{1},\cdots,T$, we feed $L_{I}$ observations
$\\{\tilde{\mathbf{y}}(t-L_{I}+1),\cdots,\tilde{\mathbf{y}}(t)\\}$ into the
OB-LSTM network and train it to predict the next $L_{P}$ imaginary
observations, i.e.,
$\\{{\mathbf{y}}_{i}(t+1),\cdots,{\mathbf{y}}_{i}(t+L_{P})\\}$. Since the
observations are related to both the unknown channels and the known reflection
matrix $\mathbf{V}(t)$, we also treat $\mathbf{V}(t)$ as the part of the
network input. Besides, the historical observations are pre-processed by
normalization to improve the performance of the multilayer perceptron models
in the proposed OB-LSTM network. After data normalization, the overall input
data sequence is given by
$\\{\mathbf{x}(1),\cdots,\mathbf{x}(t),\cdots,\mathbf{x}(L_{I})\\}$, where
$\mathbf{x}(t)=[F_{n}(\tilde{\mathbf{y}}^{T}(t)),F_{n}([\mathbf{v}_{1}^{T}(t),\cdots,$
$\mathbf{v}_{\tau_{1}}^{T}(t)])]^{T}\in\mathbb{C}^{D_{I}\times 1}$ with
$D_{I}\triangleq\tau_{1}(N+1)$ and $F_{n}(\cdot)$ denoting the dimension of
the input and the normalization function, respectively. Accordingly, we
collect a set of real observations and construct the training data-label
sample set as
$\big{\\{}\\{\mathbf{x}(j+1),\cdots,\mathbf{x}(j+L_{I})\\},\\{\mathbf{y}(j+L_{I}+1),\cdots,\mathbf{y}(j+L_{I}+L_{P})\\}\big{\\}}_{j=0}^{J-1}$.

Different from classical deep learning applications, e.g., computer vision and
natural language processing, complex data format is usually considered in
wireless communications. To handle the complex input data, the proposed OB-
LSTM network is designed to contain two sub-networks with the same structure,
which are called the R sub-network and the I sub-network,
respectively.333Alternatively, we can combine the real and imaginary parts of
the input observations into one vector as the network input, such that no sub-
networks are needed. However, the proposed sub-network structure can achieve
similar performance with less network parameters, thus is more favorable here.
In the following, we only focus on the structure of the R sub-network since
the I sub-network can be similarly designed. Specifically, each sub-network
contains one input fully-connected layer, $K$ LSTM layers and one output
fully-connected layer. The input fully-connected layer is designed to augment
the input dimension of the LSTM units, which can help to extract the high-
dimensional features of the input data and improve the prediction performance.
Let $\epsilon$ represent an expansion factor, then the output of the input
fully-connected layer, denoted by $\mathbf{x}_{h}(t)\in\mathbb{R}^{\epsilon
D_{I}\times 1}$, can be expressed as

$\mathbf{x}_{h}(t)=\text{ReLU}\big{(}\mathbf{W}_{e}\mathbf{x}(t)+\mathbf{b}_{e}\big{)},$
(20)

where $\mathbf{W}_{e}\in{\mathbb{C}}^{\epsilon D_{I}\times D_{I}}$ and
$\mathbf{b}_{e}\in\mathbb{C}^{\epsilon D_{I}\times 1}$ denote the weight and
bias (i.e., the learnable network parameters), respectively, and
$\text{ReLU}(\cdot)$ is employed as the activation function. The LSTM layers
are the core of the proposed network, since they are developed to deal with
the vanishing gradient problem that might be encountered when training
traditional RNNs and are applicable to tasks such as classifying, processing
and predicting based on time series data [31]. In this work, one LSTM layer
consists of $L_{I}$ cascading LSTM units. Each LSTM unit is composed of a
cell, an input gate, an output gate and a forget gate, where the cell
remembers information over arbitrary time interval and the three gates
regulate the flow of information into and out of the cell. As shown in Fig. 3,
in the $(t,k)$-LSTM unit (i.e., the $t$-th LSTM unit in the $k$-th LSTM
layer), the forget gate $\mathbf{f}_{t,k}$, the input gate $\mathbf{i}_{t,k}$,
the output gate $\mathbf{o}_{t,k}$, the cell state $\mathbf{c}_{t,k}$ and the
output $\mathbf{y}_{t}$ are calculated as follows:

Figure 3: The structure of one LSTM unit.

$\displaystyle\mathbf{f}_{t,k}=\sigma\left(\mathbf{W}_{t,k}^{fx}\mathbf{x}_{t,k}+\mathbf{W}_{t,k}^{fy}\mathbf{y}_{t-1,k}+\mathbf{b}_{t,k}^{f}\right),$
(21)
$\displaystyle\mathbf{i}_{t,k}=\sigma\left(\mathbf{W}_{t,k}^{ix}\mathbf{x}_{t,k}+\mathbf{W}_{t,k}^{iy}\mathbf{y}_{t-1,k}+\mathbf{b}_{t,k}^{i}\right),$
$\displaystyle\mathbf{o}_{t,k}=\sigma\left(\mathbf{W}_{t,k}^{ox}\mathbf{x}_{t,k}+\mathbf{W}_{t,k}^{oy}\mathbf{y}_{t-1,k}+\mathbf{b}_{t,k}^{o}\right),$
$\displaystyle\mathbf{c}_{t,k}=\mathbf{f}_{t,k}\otimes\mathbf{c}_{t-1,k}+\mathbf{i}_{t,k}\otimes\tanh\left(\mathbf{W}_{t,k}^{cx}\mathbf{x}_{t,k}+\mathbf{W}_{t,k}^{cy}\mathbf{y}_{t-1,k}+\mathbf{b}_{t,k}^{c}\right),$
$\displaystyle\mathbf{y}_{t,k}=\mathbf{o}_{t,k}\otimes\tanh\left(\mathbf{c}_{t,k}\right),$

where $\mathbf{x}_{t,k}$ denotes the input of this LSTM unit;
$\mathbf{W}_{t,k}^{\kappa}\in\mathbb{R}^{\epsilon D_{I}\times
1},\kappa\in\\{fx,fy,ix,iy,ox,oy,cx,cy\\}$ and
$\mathbf{b}_{t,k}^{\nu}\in\mathbb{R}^{\epsilon D_{I}\times
1},\nu\in\\{f,i,o,c\\}$ represent the corresponding weights and biases,
respectively; the sigmoid function $\sigma(\cdot)$ and the tanh function
$\text{tanh}(\cdot)$ are considered as the activation functions of the gates
and the cell, respectively. Note that the output of $\sigma(\cdot)$ varies in
$[0,1]$, hence this activation function describes how much information can
pass through this gate. Finally, the output fully-connected layer transforms
the $L_{I}$ outputs of the last LSTM layer, i.e.,
$\mathbf{y}_{1,K},\cdots,\mathbf{y}_{L_{I},K}$, into the prediction of the
$L_{P}$ imaginary observations, i.e.,

$\mathbf{y}_{p}=\mathbf{W}_{p}[\mathbf{y}_{1,K}^{T},\cdots,\mathbf{y}_{L_{I},K}^{T}]^{T}+\mathbf{b}_{p},$
(22)

where $\mathbf{y}_{p}\in\mathbb{R}^{\tau_{1}L_{P}\times 1}$ denotes the
output, $\mathbf{W}_{p}\in\mathbb{R}^{\tau_{1}L_{P}\times\epsilon D_{I}L_{I}}$
and $\mathbf{b}_{p}\in\mathbb{R}^{\tau_{1}L_{P}\times 1}$ represent the
corresponding weight and bias. The outputs of the R and I sub-networks are
regarded as the real and imaginary parts, which are then utilized to construct
the final $L_{P}$ predicted imaginary observations
$\hat{\mathbf{y}}(L_{I}+1),\cdots,\hat{\mathbf{y}}(L_{I}+L_{P})$.

To summarize, the output of the proposed OB-LSTM network can be expressed as

$\\{\hat{\mathbf{y}}(j+L_{I}+1),\cdots,\hat{\mathbf{y}}(j+L_{I}+L_{P})\\}=\text{OB-
LSTM}\big{(}\\{\mathbf{x}(j+1),\cdots,\mathbf{x}(j+L_{I})\\};\\{\bm{\Theta}_{R},\bm{\Theta}_{I}\\}\big{)},$
(23)

where $\bm{\Theta}_{R}$ $(\bm{\Theta}_{I})$ represents the collection of
network parameters included in the R (I) sub-network.

Figure 4: Overall procedure of the proposed channel prediction method.

The overall procedure of the proposed channel prediction method is depicted in
Fig. 4, where the structure of the proposed OB-LSTM network is also
illustrated.

In this work, we adopt the supervised learning strategy and construct the
following loss function which is based on the $\ell_{2}$-norm of the
difference between the future observations and their estimates

$\mathcal{L}_{2}(\bm{\Theta}_{R},\bm{\Theta}_{I})=\frac{1}{J}\sum_{j=0}^{J-1}\sum_{k=1}^{T_{P}}\|\hat{\mathbf{y}}(j+L_{I}+k)-\mathbf{y}(j+L_{I}+k)\|_{2}^{2}.$
(24)

All the network parameters in the proposed OB-LSTM network are updated by
minimizing this loss function using the Adam optimizer [32]. It is noteworthy
that the parameters of the OB-LSTM network are learned through offline
training. When the channel statistics change, the user should send pilot
symbols in the first several frames of a super frame, as such a set of real
observations can be obtained, based on which we can further finetune the
network parameters to alleviate the performance loss caused by model mismatch.

When the training is finished, the OB-LSTM network can be flexibly integrated
into the proposed 2SCPT scheme. We take the case of $T_{1}=6$ and $T_{2}=6$ as
a toy example to elaborate this. The most simple strategy (also called
_Strategy A_) is to let $L_{I}=T_{1}=6$ and $L_{P}=T_{2}=6$, then $6$ real
observations, i.e., $\mathbf{y}(1),\cdots,\mathbf{y}(6)$, collected in the
first stage can be directly input into the OB-LSTM network to predict $6$
imaginary observations, i.e., $\mathbf{y}_{i}(7),\cdots,\mathbf{y}_{i}(12)$.
However, in this strategy, if the value of $T_{1}$ or $T_{2}$ is large, the
computational complexity of the OB-LSTM network will be high since the
input/output dimensions and the number of neurons/layers are large. Thus, it
would be difficult to train such a network with good prediction performance.
An alternative strategy (also called _Strategy B_) is to set $L_{P}<L_{I}\leq
T_{1}$, e.g., $L_{I}=6$ and $L_{P}=1$, such that the predicted imaginary
observations can be reused as the network input. Specifically, in the first
time interval of the second stage, the OB-LSTM network predicts the next
imaginary observation $\mathbf{y}_{i}(7)$ based on the last 6 real
observations $\mathbf{y}(1),\cdots,\mathbf{y}(6)$ collected in the first
stage, while in the subsequent time intervals, we can construct the input data
using 1 imaginary observation and the last $5$ real observations, i.e.,
$\mathbf{y}(2),\cdots,\mathbf{y}(6),\mathbf{y}_{i}(7)$, to obtain
$\mathbf{y}_{i}(8)$. Note that the hyper-parameters of the OB-LSTM network,
e.g., $L_{I}$, $L_{P}$ and $\epsilon$, and the employed strategies should be
carefully chosen to balance the prediction performance, the computational
complexity of the network and the required channel training overhead.
Investigation into the impacts of different hyper-parameters and strategies
will be shown in Section V.

######  Remark 1.

When we consider the case where the AP is deployed with multiple antennas, the
user-AP equivalent channel matrices, i.e., the reflected channel matrices and
the direct channel vectors, are required to be tracked and predicted based on
the observations. By vectorizing the channel matrices and constructing the
observation functions as in (10) and (11), the proposed 2SCTP scheme can be
readily extended to this more general case.

##  IV 2SCTP Scheme For the General Case

In this section, we extend the proposed 2SCTP scheme to the more general case
where the IRS-AP, user-IRS and user-AP channels are all assumed to be time-
varying. In this case, the user-AP equivalent channel $\mathbf{h}(t)$ can be
written as

$\mathbf{h}(t)=[\sqrt{1-\alpha_{\text{UA}}}{h}_{d}(t-1)+\sqrt{\alpha_{\text{UA}}}{u}_{\textrm{UA}}(t),\tilde{{h}}_{1}(t),\cdots,\tilde{{h}}_{N}(t)]^{T},$
(25)

where $\tilde{h}_{n}(t),n\in[1,N]$ consists of four random variables and is
given by

$\displaystyle\tilde{{h}}_{n}(t)=$
$\displaystyle\sqrt{(1-\alpha_{\text{IA}})(1-\alpha_{\text{UI}})}[\mathbf{g}(t-1)]_{n}[\mathbf{h}_{r}(t-1)]_{n}+\sqrt{\alpha_{\text{IA}}(1-\alpha_{\text{UI}})}[\mathbf{u}_{\textrm{IA}}(t)]_{n}[\mathbf{h}_{r}(t-1)]_{n}$
(26)
$\displaystyle+\sqrt{\alpha_{\text{UI}}(1-\alpha_{\text{IA}})}[\mathbf{g}(t-1)]_{n}[\mathbf{u}_{\textrm{UI}}(t)]_{n}+\sqrt{\alpha_{\text{UI}}\alpha_{\text{IA}}}[\mathbf{u}_{\textrm{IA}}(t)]_{n}[\mathbf{u}_{\textrm{UI}}(t)]_{n}.$

Then, we can simplify the state equation (25) as

$\mathbf{h}(t)=\mathbf{A}_{g}\mathbf{h}(t-1)+\mathbf{u}_{g}(t),$ (27)

where

$\mathbf{A}_{g}=\left[\begin{array}[]{cc}\sqrt{1-\alpha_{\text{UA}}}&\mathbf{0}_{1\times
N}\\\ \mathbf{0}_{N\times
1}&\sqrt{(1-\alpha_{\text{IA}})(1-\alpha_{\text{UI}})}\mathbf{I}\end{array}\right],$
(28)

$\mathbf{u}_{g}(t)$ represents the composite noise variable that satisfies
$[\mathbf{u}_{g}(t)]_{1}={u}_{\textrm{UA}}(t)$ and
$[\mathbf{u}_{g}(t)]_{n}=\sqrt{\alpha_{\text{IA}}(1-\alpha_{\text{UI}})}$
$[\mathbf{u}_{\textrm{IA}}(t)]_{n}[\mathbf{h}_{r}(t-1)]_{n}+\sqrt{\alpha_{\text{UI}}(1-\alpha_{\text{IA}})}[\mathbf{g}(t-1)]_{n}[\mathbf{u}_{\textrm{UI}}(t)]_{n}+\sqrt{\alpha_{\text{UI}}\alpha_{\text{IA}}}[\mathbf{u}_{\textrm{IA}}(t)]_{n}$
$[\mathbf{u}_{\textrm{UI}}(t)]_{n},n\in[2,N+1]$. The covariance matrix of
$\mathbf{u}_{g}(t)$ is denoted as
$\mathbf{C}_{g}=\mathbb{E}\big{[}\mathbf{u}_{g}\mathbf{u}_{g}^{H}\big{]}$. It
is readily seen that although the state equation (27) is linear, the elements
in $\mathbf{h}(t)$ do not follow complex Gaussian distributions in general,
which implies that the state variable $\mathbf{h}(t)$ and noise variable
$\mathbf{u}_{g}(t)$ are not Gaussian vectors. This is quite different from the
state equation (7) discussed in the special case, and the KF-based channel
tracking algorithm cannot be directly applied to this general case. Note that
for the non-Gaussian scenario, the particle filter (PF) method is widely-used
[33], such as in target tracking, signal processing and automatic control,
etc. The PF method uses a set of particles, i.e., observations, to represent
the posterior distribution of the state variable which follows an arbitrary
distribution. However, collecting these particles can be quite time-consuming,
thus applying the PF method for the general case can be very inefficient.

To tackle this difficulty, we propose in this paper a GKF-based channel
tracking algorithm for the general case.444The proposed GKF-based channel
tracking algorithm can be extended to handle other channel models as well as
long as the state equation is linear. Our idea is inspired by the extended KF
(EKF) and unscented KF (UKF) methods [34], both of which are designed for non-
linear systems. Their difference is that the EKF method employs multivariate
Taylor series expansions to linearize the state equation at the working point,
while the UKF method tries to approximate the probability density function of
the output of the non-linear function included in the system by a number of
deterministic sampling points which represent the underlying distribution as a
Gaussian distribution. Then, the original KF method can work on these modified
or linearized systems. In the considered problem, our state equation and
observation function are both linear, but the composite noise variable
$\mathbf{u}_{g}(t)$ in the state equation is not Gaussian. Therefore, some
necessary modifications should be made on $\mathbf{u}_{g}(t)$ such that the
resulting problem can be solved by the KF method. Specifically, to tackle the
difficulty caused by the non-Gaussian noise variable $\mathbf{u}_{g}(t)$, we
propose to use a complex Gaussian distribution to approximate the distribution
of $\mathbf{u}_{g}(t)$ in each time interval. Then, the predicted state
variable $\mathbf{h}(t)$, obtained from the linear state equation (27), also
follows a complex Gaussian distribution. With the help of such an
approximation, we can apply the KF method to address the resulting channel
tracking problem. In the following, we will first present the proposed
approximation method and then introduce the GKF-based channel tracking
algorithm.

###  IV-A Channel Distribution Analysis

In this subsection, we aim to analyse the distribution of the noise variable
$\mathbf{u}_{g}(t)$ included in the state function (27). First, it is readily
seen that $[\mathbf{u}_{g}(t)]_{1}$ follows the complex Gaussian distribution,
i.e.,
$[\mathbf{u}_{g}(t)]_{1}\sim\mathcal{CN}(0,\alpha_{\text{UA}}\beta_{\text{UA}})$.
The other elements in $\mathbf{u}_{g}(t)$ can be expressed as the sum of three
random variables, i.e.,

$[\mathbf{u}_{g}(t)]_{n}=[\mathbf{u}_{g,1}(t)]_{n}+[\mathbf{u}_{g,2}(t)]_{n}+[\mathbf{u}_{g,3}(t)]_{n},$
(29)

where
$[\mathbf{u}_{g,1}(t)]_{n}\triangleq\sqrt{\alpha_{\text{IA}}(1-\alpha_{\text{UI}})}[\mathbf{h}_{r}(t-1)(t)]_{n}[\mathbf{u}_{\textrm{IA}}(t)]_{n}$,
$[\mathbf{u}_{g,2}(t)]_{n}\triangleq\sqrt{\alpha_{\text{UI}}(1-\alpha_{\text{IA}})}[\mathbf{g}(t-1)]_{n}[\mathbf{u}_{\textrm{UI}}(t)]_{n}$
and
$[\mathbf{u}_{g,3}(t)]_{n}\triangleq\sqrt{\alpha_{\text{IA}}\alpha_{\text{UI}}}[\mathbf{u}_{\textrm{IA}}(t)]_{n}[\mathbf{u}_{\textrm{UI}}(t)]_{n}$.
We can see that $[\mathbf{u}_{g,1}(t)]_{n}$ and $[\mathbf{u}_{g,2}(t)]_{n}$
both follow complex Gaussian distributions, i.e.,
$[\mathbf{u}_{g,1}(t)]_{n}\sim\mathcal{CN}(0,\alpha_{\text{IA}}(1-\alpha_{\text{UI}})\sigma^{2}_{\text{IA}}|[\mathbf{h}_{r}(t-1)]_{n}|^{2})$,
$[\mathbf{u}_{g,2}(t)]_{n}\sim\mathcal{CN}(0,\alpha_{\text{UI}}(1-\alpha_{\text{IA}})|[\mathbf{g}_{r}(t-1)]_{n}|^{2}\sigma^{2}_{\text{UI}})$,
while the third random variable $[\mathbf{u}_{g,3}(t)]_{n}$ is the product of
two independent complex Gaussian random variables, i.e.,
$\sqrt{\alpha_{\text{IA}}}[\mathbf{u}_{\textrm{IA}}(t)]_{n}$ and
$\sqrt{\alpha_{\text{UI}}}[\mathbf{u}_{\textrm{UI}}(t)]_{n}$.

Figure 5: Marginal probability distribution of the real and imaginary parts of
$X,\\{X_{1},X_{2}\\}$ and $\tilde{X}$.

Then, we focus on the analysis of $[\mathbf{u}_{g,3}(t)]_{n}$. For clarity, we
refer to the distribution of this type of random variables as the product
complex Gaussian distribution in the following. Note that the product of two
real Gaussian distributions, named as product real Gaussian distribution, has
been studied before in [35, 36], where the authors proved that it can be
expressed in terms of Meijer G-functions. Consider a complex random variable
$X$ that satisfies $X=X_{1}X_{2}$, where
$X_{1},X_{2}\sim\mathcal{CN}(0,\frac{1}{4})$. Since it is difficult to
directly analyse the distribution of $X$, we first investigate it via Monte-
Carlo simulation using $10^{5}$ samples. Fig 5 demonstrates the marginal
probability distributions of the real and imaginary parts of $X$ and
$X_{1}/X_{2}$. It is observed that the joint (or marginal) distribution
resembles a shaper and slimmer complex Gaussian (or Gaussian) distribution,
which inspires us to approximate the product complex Gaussian distribution by
a simple complex Gaussian distribution with the same mean and variance.
However, since $[\mathbf{u}_{g,3}(t)]_{n}$ is in essence a component of the
system noise variable $[\mathbf{u}_{g}(t)]_{n}$, we cannot sample this random
variable independently and calculate its mean and variance. Hence, the ideas
employed in the PF and UKF methods, i.e., sampling a set of particles to
represent the posterior distribution of the state variable, cannot be applied
here to approximate the distribution of $[\mathbf{u}_{g,3}(t)]_{n}$. To
overcome this challenge, we resort to the following result about the mean and
variance of the random variable $X$, i.e., $\mu_{P}$ and $\sigma^{2}_{P}$,

$\displaystyle\mu_{P}=\mathbb{E}\big{[}X_{1}\big{]}\mathbb{E}\big{[}X_{2}\big{]}=0,$
(30)
$\displaystyle\sigma^{2}_{P}=\mathbb{E}\big{[}XX^{*}\big{]}=\mathbb{E}\big{[}X_{1}X_{2}X_{2}^{*}X_{1}^{*}\big{]}=\mathbb{E}\big{[}X_{1}X_{1}^{*}X_{2}X_{2}^{*}\big{]}=\mathbb{E}\big{[}X_{1}X_{1}^{*}\big{]}\mathbb{E}\big{[}X_{2}X_{2}^{*}\big{]}=\sigma^{2}_{1}\sigma^{2}_{2}.$

As such, the mean and variance of $[\mathbf{u}_{g,3}(t)]_{n}$ can be directly
obtained from those of
$\sqrt{\alpha_{\text{IA}}}[\mathbf{u}_{\textrm{IA}}(t)]_{n}$ and
$\sqrt{\alpha_{\text{UI}}}[\mathbf{u}_{\textrm{UI}}(t)]_{n}$. Then, based on
(30) and the numerical results shown in Fig. 5, a random product complex
Gaussian variable $X$ with mean $\mu_{P}$ and variance $\sigma^{2}_{P}$ can be
approximated by a complex Gaussian variable
$\tilde{X}\sim\mathcal{CN}(\mu_{P},\sigma^{2}_{P})$. As can be seen, although
the distribution of $\tilde{X}$ is similar to that of $X$, it is still not
sharp enough as compared to the true distribution. However, since only one
component in $[\mathbf{u}_{g}(t)]_{n}$, i.e., $[\mathbf{u}_{g,3}(t)]_{n}$,
requires to be approximated and the scaling coefficient of
$[\mathbf{u}_{g,3}(t)]_{n}$, i.e.,
$\sqrt{\alpha_{\text{IA}}\alpha_{\text{UI}}}$, is smaller than those of
$[\mathbf{u}_{g,1}(t)]_{n}$ and $[\mathbf{u}_{g,2}(t)]_{n}$, i.e.,
$\sqrt{\alpha_{\text{IA}}(1-\alpha_{\text{UI}})}$ and
$\sqrt{(1-\alpha_{\text{IA}})\alpha_{\text{UI}}}$, we can infer that the
negative effect caused by the approximation error should be limited. In
Section V, we will show that the proposed GKF-based channel tracking algorithm
can achieve good performance although under the employed approximation.

###  IV-B Complex Gaussian Approximation

Next, we focus on approximating the distribution of $[\mathbf{u}_{g}(t)]_{n}$
based on the result in Section IV-A. According to the analysis in the previous
subsection, $[\mathbf{u}_{g,3}(t)]_{n}$ can be approximated by
$[\hat{\mathbf{u}}_{g,3}(t)]_{n}$ which satisfies
$[\hat{\mathbf{u}}_{g,3}(t)]_{n}\sim\mathcal{CN}(0,\alpha_{\text{IA}}\alpha_{\text{UI}}\sigma^{2}_{\text{IA}}\sigma^{2}_{\text{UI}})$,
thus $[\mathbf{u}_{g}(t)]_{n}$ can be approximated by
$[\hat{\mathbf{u}}_{g}(t)]_{n}\triangleq[\mathbf{u}_{g,1}(t)]_{n}+[\mathbf{u}_{g,2}(t)]_{n}+[\hat{\mathbf{u}}_{g,3}(t)]_{n},n\in[2,N+1]$.
It is noteworthy that the temporal correlation coefficients, i.e.,
$\alpha_{\text{IA}}$, $\alpha_{\text{UI}}$, $\alpha_{\text{UA}}$, and the
variances of $[\mathbf{u}_{\text{IA}}]_{n}$, $[\mathbf{u}_{\text{UI}}]_{n}$
and $[\mathbf{u}_{\text{UA}}]_{n}$, i.e., $\sigma^{2}_{\text{IA}}$,
$\sigma^{2}_{\text{UA}}$, $\sigma^{2}_{\text{UI}}$, are prior knowledges that
are known. Since the sum of Gaussian distributions is still Gaussian, the
covariance matrix of $\mathbf{u}_{g}(t)$, i.e., $\mathbf{C}_{g}(t)$, can be
approximated by

$\hat{\mathbf{C}}_{g,\text{I}}(t)=\text{diag}\big{(}\alpha_{\text{UA}}\sigma^{2}_{\text{UA}},\sigma^{2}_{g,1}(t),\cdots,\sigma^{2}_{g,N}(t)\big{)},$
(31)

where

$\sigma^{2}_{g,n}(t)=\alpha_{\text{IA}}(1-\alpha_{\text{UI}})|[\mathbf{h}_{r}(t-1)]_{n}|^{2}\sigma^{2}_{\text{IA}}+\alpha_{\text{UI}}(1-\alpha_{\text{IA}})|[\mathbf{g}(t-1)]_{n}|^{2}\sigma^{2}_{\text{UI}}+\alpha_{\text{IA}}\alpha_{\text{UI}}\sigma^{2}_{\text{IA}}\sigma^{2}_{\text{UI}}.$
(32)

This is referred to as _type one complex Gaussian approximation (CGA-I)_ in
the following. It can be seen that the amplitudes of
$[\mathbf{h}_{r}(t-1)]_{n}$ and $[\mathbf{g}(t-1)]_{n}$ are required in order
to calculate $\sigma^{2}_{g,n}(t)$, which is however difficult to obtain in
the tracking process due to the passive nature of IRS, i.e., the CSI of the
IRS-AP and user-IRS channels are difficult to obtain. To address this issue,
we propose to approximate
$\alpha_{\text{IA}}(1-\alpha_{\text{UI}})|[\mathbf{h}_{r}(t-1)]_{n}|^{2}\sigma^{2}_{\text{IA}}+\alpha_{\text{UI}}(1-\alpha_{\text{IA}})|[\mathbf{g}(t-1)]_{n}|^{2}\sigma^{2}_{\text{UI}}$
using its lower bound, which is based on the following theorem.

######  Theorem 1.

Assuming that $a=a_{r}+ja_{i}$ and $b=b_{r}+jb_{i}$ are two arbitrary complex
scalars, then we have

$|a|^{2}+|b|^{2}\geq\big{|}\Re(ab)\big{|}+\big{|}\Im(ab)\big{|},$ (33)

###### Proof.

Please relegate to Appendix. ∎

Based on Theorem 1, we employ a random variable
$[\tilde{\mathbf{u}}_{g}(t)]_{n}\sim\mathcal{CN}(0,\tilde{\sigma}^{2}_{g,n}(t))$
to approximate $[{\mathbf{u}}_{g}(t)]_{n}$, where
$\tilde{\sigma}^{2}_{g,n}(t)\leq\text{Var}\big{[}[{\mathbf{u}}_{g}(t)]_{n}\big{]}$.555
Note that in the KF method, the covariance matrix of the noise variable, i.e.,
$\mathbf{u}_{g}(t)$, is used to measure the estimation error of the state
variable, i.e., $\mathbf{h}(t)$. Therefore, approximating
$\text{Var}\big{[}[{\mathbf{u}}_{g}(t)]_{n}\big{]}$ using the lower bound
$\tilde{\sigma}^{2}_{g,n}(t)$ will make the estimated variable more stable.
Specifically, the value of $\tilde{\sigma}^{2}_{g,n}(t)$ is given by

$\tilde{\sigma}^{2}_{g,n}(t)=\delta_{1}\Big{(}\big{|}\Re([\mathbf{h}(t-1)]_{n})\big{|}+\big{|}\Im([\mathbf{h}(t-1)]_{n})\big{|}\Big{)}+\delta_{2},$
(34)

where
$\delta_{1}\triangleq\sqrt{\alpha_{\text{IA}}(1-\alpha_{\text{IA}})\alpha_{\text{UI}}(1-\alpha_{\text{UI}})}\sigma_{\text{IA}}\sigma_{\text{UI}}$
and
$\delta_{2}\triangleq\alpha_{\text{IA}}\alpha_{\text{UI}}\sigma^{2}_{\text{IA}}\sigma^{2}_{\text{UI}}$.
Then, the covariance matrix $\mathbf{C}_{g}(t)$ can be approximated by

$\hat{\mathbf{C}}_{g,\text{II}}(t)=\left[\begin{array}[]{cc}\alpha_{\text{UA}}\sigma^{2}_{\text{UA}}&\mathbf{0}_{1\times
N}\\\ \mathbf{0}_{N\times
1}&\delta_{1}\text{diag}\big{(}\big{|}\Re([\mathbf{h}(t-1)]_{2:N+1})\big{|}+\big{|}\Im([\mathbf{h}(t-1)]_{2:N+1})\big{|}\big{)}+\delta_{2}\mathbf{I}_{N}\end{array}\right],$
(35)

where
$[\mathbf{h}(t-1)]_{2:N+1}\triangleq\big{[}[\mathbf{h}_{2,N}(t-1)]_{2},\cdots,[\mathbf{h}(t-1)]_{N+1}\big{]}^{T}\in\mathbb{C}^{N\times
1}$. This is named as _type two complex Gaussian approximation (CGA-II)_ in
the following. Note that CGA-I and CGA-II are different in the sense that the
IRS-AP and user-IRS channel coefficients are required in CGA-I, while they are
not required in CGA-II. However, the performance achieved by CGA-I and CGA-II
is quite close, as will be shown in the simulation results.

###  IV-C GKF-based Channel Tracking

Finally, based on the analysis presented in Section IV-A and IV-B, we extend
the KF-based channel tracking algorithm (i.e., Algorithm 1) to the general
case and propose in Algorithm 2 the GKF-based channel tracking algorithm.
Different from Algorithm 1 where only the Kalman gain
$\mathbf{G}_{\text{KF}}(t)$, the correction $\mathbf{h}_{\text{KF}}(t)$ and
the estimation covariance $\mathbf{M}_{\text{KF}}(t)$ are updated in each time
interval, the covariance matrix $\hat{\mathbf{C}}_{gl}(t)$ is also required to
be updated in Algorithm 2 since it varies with time in the general case.
According to (31), the estimation of $\hat{\mathbf{C}}_{gl}(t)$ can be
obtained by

$\hat{\mathbf{C}}_{gl}(t)=\left[\begin{array}[]{cc}\alpha_{\text{UA}}\sigma^{2}_{\text{UA}}&\mathbf{0}_{1\times
N}\\\ \mathbf{0}_{N\times
1}&\delta_{1}\text{diag}\big{(}\big{|}\Re([\mathbf{h}_{\text{KF}}(t-1)]_{2:N+1})\big{|}+\big{|}\Im([\mathbf{h}_{\text{KF}}(t-1)]_{2:N+1})\big{|}\big{)}+\delta_{2}\mathbf{I}_{N}\end{array}\right].$
(36)

In the second stage of the proposed transmission protocol, imaginary
observations can be similarly obtained by the OB-LSTM network as in Section
III, which serve as the input of Algorithm 2, then the channels in this stage
can be estimated and no channel training overhead is required.

Algorithm 2 GKF-based Channel Tracking Algorithm

Input: covariance matrix
$\mathbf{C}_{h}=\mathbb{E}\big{[}\mathbf{h}\mathbf{h}^{H}\big{]}$ and
observations $\tilde{\mathbf{y}}(t)$
Output: Estimated channel vector $\mathbf{h}_{\text{KF}}(t)$

1:Initialization: $t=0$, $\mathbf{h}_{\text{KF}}(-1)=\mathbf{0}$,
$\mathbf{M}_{\text{MF}}=\mathbb{E}\big{[}\mathbf{h}\mathbf{h}^{H}\big{]}$,
$\hat{\mathbf{C}}_{gl}(0)=\mathbb{E}\big{[}\mathbf{h}\mathbf{h}^{H}\big{]}$

2:while $t\geq 1$ do _Prediction step_ :

3: Calculate the prediction according to (14),

4: Calculate the prediction covariance matrix according to
$\mathbf{M}(t)=\mathbb{E}\big{[}\big{(}\mathbf{h}(t)-\hat{\mathbf{h}}(t)\big{)}\big{(}\mathbf{h}(t)-\hat{\mathbf{h}}(t)\big{)}^{H}\big{]}=\mathbf{A}_{h}\mathbf{M}_{\text{KF}}(t-1)\mathbf{A}_{h}^{T}+\mathbf{B}_{h}\hat{\mathbf{C}}_{gl}(t)\mathbf{B}_{h}^{T}$,
_Update step_ :

5: Update the Kalman gain vector according to (16),

6: Update the correction according to (17),

7: Update the estimation covariance matrix according to (18),

8: Update the approximated covariance matrix of $\mathbf{u}_{g}$,
$\hat{\mathbf{C}}_{gl}(t+1)$, according to (36),

9: $t=t+1$,

10:end while

##  V Simulation Results

In this section, we present numerical results to validate the effectiveness of
the proposed two-stage transmission protocol and the 2SCTP scheme. We consider
an IRS-aided wireless communication system as shown in Fig. 6, where a three-
dimensional coordinate system is assumed and the AP and IRS are located on the
$x$-axis and $y-z$ plane, respectively. The reference antenna/element at the
AP/IRS are located at $(3$ $\text{m},0,0)$ and $(0,50$ $\text{m},2$
$\text{m})$, and the user is located at $(2$ $\text{m},50$ $\text{m},0)$. In
the simulations, the IRS is equipped with a $5\times 7$ uniform rectangular
array.666If the number of IRS reflecting elements is large, the channel
training overhead required to obtain the observations and the number of
network parameters included in the OB-LSTM network will also be increased. In
this case, we can employ the grouping and partition method [37] to achieve a
good trade-off between channel tracking performance and channel training
overhead/network complexity. The reference distance and the pass loss at the
reference distance are set as $d_{0}=1$ m and $l_{0}=-30$ dB. The IRS-AP,
user-IRS and user-AP path-loss exponents are fixed to
$\gamma^{\text{IA}}=2.2$, $\gamma^{\text{UI}}=2.2$ and
$\gamma^{\text{UA}}=3.6$, respectively. The transmit power and noise variance
are set to $p=26$ dBm and $\sigma^{2}=-80$ dBm. Furthermore, the variances of
the elements of $\mathbf{u}_{\text{IA}}$, $\mathbf{u}_{\text{UI}}$ and
${u}_{\text{UA}}$ are respectively set to
$\sigma^{2}_{\text{IA}}=l_{\text{IA}}$, $\sigma^{2}_{\text{UA}}=l_{\text{UA}}$
and $\sigma^{2}_{\text{UI}}=l_{\text{UI}}$. In all our simulations, the
definitions of the normalized mean square error (NMSE) and average NMSE
(ANMSE) in the $t$-th time interval are given by
$\text{NMSE}(t)=\frac{\|\mathbf{h}_{\text{KF}}(t)-{\mathbf{h}}(t)\|^{2}}{\|\mathbf{h}(t)\|^{2}},\;\;\;\;\text{ANMSE}(t)=\frac{1}{t}\sum_{i=1}^{t}\frac{\|\mathbf{h}_{\text{KF}}(i)-{\mathbf{h}}(i)\|^{2}}{\|\mathbf{h}(i)\|^{2}}.$

Figure 6: Simulation setup.

###  V-A Special Case

In this subsection, we consider the special case where the IRS-AP channel is
assumed to stay unchanged, i.e., $\alpha^{\text{IA}}=0$, and the temporal
correlation coefficients of the user-IRS and user-AP channels are set as
$\alpha^{\text{UI}}=0.01$ and $\alpha^{\text{UA}}=0.01$, respectively.

First, we investigate the channel tracking performance in the first stage.
Fig. 8 shows the NMSE performance of the KF-based channel tracking algorithm
(i.e., Algorithm 1) with different numbers of time intervals. The number of
time slots in each time interval for transmitting pilot symbols is assumed to
be $\tau_{1}\in\\{2,6,10\\}$. For comparison, we also provide the performance
achieved by the channel estimation (CE) scheme proposed in [9] as the
benchmark, i.e., in each time interval, $N+1$ time slots are allocated to
transmit pilot symbols and the channels are estimated by using the DFT
reflection pattern and the MMSE channel estimation algorithm. As can be seen,
the performance of Algorithm 1 improves with the increasing of $\tau_{1}$ and
in the meantime, higher convergence speed can be achieved. For example, when
$\tau_{1}=10$, the number of time intervals required by Algorithm 1 to achieve
convergence is less than 10, while when $\tau_{1}=2$, this number increases to
20. Moreover, with larger $\tau_{1}$, Algorithm 1 shows more stable NMSE
performance after convergence, e.g., the fluctuation of the NMSE curve when
$\tau_{1}=10$ is smaller than that when $\tau_{1}=2$. Furthermore, as compared
to the benchmark, Algorithm 1 can achieve much lower (150 time slots versus
525 time slots if 15 time intervals and $\tau_{1}=10$ are considered) channel
training overhead with only minor NMSE performance loss.

Figure 7: Channel tracking performance in terms of NMSE when different
$\tau_{1}$ are assumed in the first stage (special case).

Figure 8: Channel tracking performance in terms of ANMSE when two different
measurement matrices are considered (special case).

In Fig. 8, we investigate the ANMSE performance of Algorithm 1 in the first
stage when two different measurement matrices $\mathbf{V}(t)$ are employed,
i.e., the random measurement matrix and the DFT measurement matrix. In the
random measurement matrix, each element of the reference matrix $\mathbf{Q}$
in (12) is generated from a complex Gaussian distribution $\mathcal{CN}(0,1)$,
while in the DFT measurement matrix, $\mathbf{Q}$ is an $(N+1)$-DFT matrix.
$\tau_{1}$ is fixed to 6. One can see that as compared to using random
measurement matrix, Algorithm 1 with the DFT measurement matrix exhibits
better convergence. This is due to the fact that all the columns of the DFT
measurement matrix are independent with each other, which forms a better
sensing matrix.

Second, we investigate the channel prediction performance in the second stage.
In our simulations, the number of LSTM layers in the OB-LSTM network is set to
$K=4$. The DFT measurement matrix is used due to its superiority demonstrated
in Fig. 8. The number of training samples is set to $10^{4}$, and the
batchsize and learning rate are fixed to 10 and 0.0001, respectively. The
training process is conducted offline on a Windows server with Intel Xeon Gold
6230 CPU and an Nvidia 2080Ti GPU, and the proposed networks are implemented
in Python using the TensorFlow library with the Adam optimizer.

We present in Fig. 9 the convergence of the proposed OB-LSTM network during
training, where the $\mathcal{L}_{2}$ loss in (24) when feeding the validation
data set into the OB-LSTM network is regarded as the performance metric. We
consider four different configurations of the network structure, i.e., the
scaling factor of the LSTM units $\epsilon$ is set to be varying in
$[1,3,5,10]$. Besides, two different configurations of the lengths of the
input observations and predicted imaginary observations are considered: (a)
$L_{I}=6,L_{P}=1$, and (b) $L_{I}=6,L_{P}=6$. From Fig. 9, we can observe that
the $\mathcal{L}_{2}$ loss achieved by the OB-LSTM network is able to converge
in about 20000 iterations. Besides, the performance of the OB-LSTM network
improves with the increasing of the scaling factor $\epsilon$. However, larger
$\epsilon$ also means higher computational complexity. Therefore, the hyper-
parameter $\epsilon$ should be carefully chosen before training to strike a
good balance between performance and computational complexity. Furthermore, by
comparing Fig. 9 (a) and (b), we can see that with the same value of
$\epsilon$, the OB-LSTM network in case (a) outperforms that in case (b). This
implies that the prediction performance of the OB-LSTM network deteriorates as
the length of the predicted imaginary observations $L_{P}$ increases, and
there is a tradeoff between channel training overhead (i.e., prediction
length) and prediction performance.

Figure 9: Convergence of the proposed OB-LSTM network during training.

Next, we exhibit in TABLE I the NMSE of the predicted imaginary observations
(OB-NMSE) achieved by the proposed OB-LSTM network, when the number of time
intervals allocated for the two stages are fixed to $T_{1}=6$ and $T_{2}=6$,
respectively. We consider _Strategy A_ and _Strategy B_ discussed in Section
III-B which meet the requirement of $T_{1}=6$ and $T_{2}=6$. As shown in TABLE
I, the OB-NMSE achieved by _Strategy A_ is quite stable over all the
considered time intervals, while that by _Strategy B_ gradually deteriorates
over time. Besides, _Strategy B_ outperforms _Strategy A_ in the first three
time intervals, while _Strategy A_ shows better performance in the other time
intervals. The former is due to the non-ideal processing introduced by the OB-
LSTM network when the output dimension $L_{P}$ is larger ($L_{P}=6$) in
_Strategy A_. The latter is because the OB-LSTM network in _Strategy B_
outputs the following imaginary observations one by one based on both real and
imaginary observations, as such the prediction errors in the previous time
intervals will deteriorate the prediction performance in the subsequent time
intervals and thus _Strategy B_ becomes worse than _Strategy A_ in the later
time intervals.

TABLE I: Observation prediction performance in the second stage. OB-NMSE | Time interval index in the second stage
---|---
1 | 2 | 3 | 4 | 5 | 6
_Strategy A_ | 0.1058 | 0.1263 | 0.1208 | 0.1098 | 0.1026 | 0.1541
_Strategy B_ | 0.0297 | 0.0579 | 0.0748 | 0.1161 | 0.1235 | 0.1449
Figure 10: NMSE performance of the proposed 2SCTP scheme in different
scenarios.

Then, in Fig. 10, we investigate the overall NMSE performance of the proposed
2SCTP scheme in two different scenarios, i.e., Scenario A and Scenario B.
Specifically, in Scenario A, we set $T_{1}=6$, $T_{2}=6$, $L_{I}=6$, $L_{P}=6$
and _Strategy A_ is employed to generate the predicted imaginary observations;
while in Scenario B, we set $T_{1}=6$, $T_{2}=3$, $L_{I}=6$, $L_{P}=1$ and
_Strategy B_ is employed. For comparison, we consider a pure channel tracking
(CT) scheme, where there is no channel prediction, i.e., all the time
intervals are allocated to the first stage. Note that the ratio
$\frac{T_{2}}{T_{1}}$ of Scenario A is larger than that of Scenario B, thus
lower channel training overhead can be achieved by Scenario A. As shown in
Fig. 10, in the first stage, the NMSE performance of the 2SCTP scheme in
Scenario B is the same to that in Scenario A, however in the second stage, the
performance in Scenario B is better than that in Scenario A and it is quite
close to that of the benchmark. Therefore, for the proposed 2SCTP scheme,
adopting the parameters in Scenario B can effectively reduce the channel
training overhead with almost no performance loss.

At last, Table II compares the channel training overhead required by the 2SCTP
scheme in Scenarios A and B when the number of time intervals is fixed to
$T=3600$. The amount of channel training overhead is measured by the number of
time slots allocated for pilot transmission. The overheads required by the CE
and CT schemes are regarded as benchmarks. It can be observed from Table II
that in Scenario B, the 2SCTP scheme requires less channel training overhead
than that in Scenario A, however the NMSE performance is worse in Scenario B
(see Fig. 10). Moreover, the proposed 2SCTP scheme is able to reduce the
channel training overhead significantly as compared to the CE and CT schemes.
In particular, the overhead required in Scenario B is only 8.3% and 50.0% of
those by the CE and CT schemes, respectively.

.

TABLE II: Channel Training Overhead Comparison. Scheme | 2SCTP scheme (Scenario A) | 2SCTP scheme (Scenario B)
---|---|---
Channel training overhead | $1.44\times 10^{4}$ | $1.08\times 10^{4}$
Scheme | CT scheme | CE scheme
Channel training overhead | $2.16\times 10^{4}$ | $1.30\times 10^{5}$

###  V-B General Case

Figure 11: Marginal distributions of the random variable
$\Re([\mathbf{u}_{g}]_{2})$.

Figure 12: Channel tracking process in terms of NMSE where different
approximated covariance matrices of $\mathbf{u}_{g}$ are assumed (general
case).

In this subsection, we consider the general case and set the temporal
correlation coefficients as
$\alpha^{\text{IA}}=\alpha^{\text{UI}}=\alpha^{\text{UA}}=0.01$.

First, we show in Fig. 12 the real distribution of $\mathbf{u}_{g}$ and its
approximations using CGA-I in (31) and CGA-II in (35). For simplicity, we take
the marginal distribution of the real and imaginary parts of
$[\mathbf{u}_{g}]_{2}$ as an example. Besides, since the marginal
distributions of $\Re([\mathbf{u}_{g}]_{2})$ and $\Im([\mathbf{u}_{g}]_{2})$
are the same, we only plot the distribution of $\Re([\mathbf{u}_{g}]_{2})$ and
its approximation in Fig. 12. It is observed that the distribution achieved by
CGA-I is almost the same with the real distribution, while the distribution
achieved by CGA-II has a smaller variance than the real distribution.

Then, Fig. 12 shows the performance achieved by the proposed GKF-based channel
tracking algorithm (i.e., Algorithm 2), when two different approximated
covariance matrices of $\mathbf{u}_{g}$ are used, i.e.,
$\hat{\mathbf{C}}_{g,\text{I}}$ given in (31) and $\hat{\mathbf{C}}_{gl}$
given in (36), which are obtained via CGA-I and CGA-II, respectively. It is
noteworthy that $\mathbf{h}_{r}(t-1)$ and $\mathbf{g}(t-1)$ are needed to
calculate $\hat{\mathbf{C}}_{g,\text{I}}(t)$, while only
$\mathbf{h}_{\text{KF}}(t-1)$ is required to obtain
$\hat{\mathbf{C}}_{gl}(t)$. To employ $\hat{\mathbf{C}}_{g,\text{I}}$ in
Algorithm 2, we replace (36) with (31) in step 8 to update the approximated
covariance matrix of $\mathbf{u}_{g}$, and $\mathbf{h}_{r}(t-1)$ and
$\mathbf{g}(t-1)$ are assumed to be known.777Due to the passive nature of IRS,
the IRS-AP channel $\mathbf{g}$ and the user-IRS channel $\mathbf{h}_{r}$ are
difficult to be estimated. Here, to show the performance of the proposed CGA-I
method, we assume that $\mathbf{g}$ and $\mathbf{h}_{r}$ are known. The number
of time slots allocated for the channel training phase is set to $\tau_{1}=6$.
As shown in Fig. 12, by using $\hat{\mathbf{C}}_{gl}$ and
$\hat{\mathbf{C}}_{g,\text{I}}$, Algorithm 2 can achieve similar NMSE
performance, and the converged NMSE can reach 0.05 after about 10 time
intervals.

Next, in Fig. 14, we investigate the NMSE performance achieved by Algorithm 2
with different values of $\tau_{1}$, i.e., $\tau_{1}\in[2,4,6,8]$. Similar to
Fig. 8, we regard the NMSE performance of the CE scheme as the benchmark. It
can be observed that when $\tau_{1}$ is small, i.e., $\tau_{1}=2$, Algorithm 2
almost fails to track the varying channel. When $\tau_{1}=4,6,8$, the NMSE
performance achieved by Algorithm 2 can achieve convergence within $10$ time
intervals. Besides, as $\tau_{1}$ increases, the convergence will become
faster and more stable. However, there is a performance gap (though not large)
between Algorithm 2 and the benchmark, which is because the number of pilots
required by Algorithm 2 is much less than that of the CE scheme.

Figure 13: Channel tracking process in terms of NMSE when different $\tau_{1}$
are assumed in the first stage (general case).

Figure 14: NMSE performance of the proposed 2SCTP scheme in different
scenarios (the general case).

Finally, Fig. 14 presents the NMSE performance of the proposed 2SCTP scheme in
the general case by combining Algorithm 2 and the OB-LSTM network. Different
from that in Fig. 10, we consider two scenarios with the same amount of
channel training overhead, i.e., in Scenario A, we set $T_{1}=6$, $T_{2}=3$,
$L_{I}=6$, $L_{P}=3$ and employ _Strategy A_ to generate the predicted
imaginary observations; while in Scenario B, we set $T_{1}=6$, $T_{2}=3$,
$L_{I}=6$, $L_{P}=1$ and use _Strategy B_. It can be seen that in the general
case, the proposed 2SCTP scheme can also achieve very similar NMSE performance
with the CT scheme, yet with much lower channel training overhead.

##  VI Conclusion

In this paper, we investigated the CTP problem in an IRS-aided wireless
communication system with time-varying channel and designed an innovative two-
stage transmission protocol. Based on the proposed transmission protocol, we
proposed a novel 2SCTP scheme to track and predict the channels with low
channel training overhead. By exploiting the temporal correlation of the
channels, we proposed a KF-based channel tracking algorithm for the special
case when the IRS-AP channel is static, while for the general case, we
developed a GKF-based channel tracking algorithm by devising a simple Gaussian
approximation method. Furthermore, we presented an LSTM-based neural network,
namely the _OB-LSTM_ network, to predict the imaginary observations based on
which the channels can be estimated by applying the KF/GKF-based channel
tracking algorithms. Numerical results showed that the proposed 2SCTP scheme
is able to outperform the existing CE and CT schemes significantly in terms of
channel training overhead.

###### Proof.

First, the quadratic sum of $a$ and $b$, i.e., $|a|^{2}+|b|^{2}$, can be
rewritten as

$|a|^{2}+|b|^{2}=\frac{1}{2}\big{|}(a_{r}^{2}+b_{r}^{2})+(a_{i}^{2}+b_{i}^{2})\big{|}+\frac{1}{2}\big{|}(a_{r}^{2}+b_{i}^{2})+(a_{i}^{2}+b_{r}^{2})\big{|}.$
(37)

According to the arithmetic and geometric (AM-GM) inequality [38],
$x^{2}+y^{2}\geq 2xy$ and $x^{2}+y^{2}\geq-2xy$ hold for all
$x,y\in\mathbb{R}$. As such, we have

$\displaystyle\left.\begin{aligned}
&(a_{r}^{2}+b_{r}^{2})+(a_{i}^{2}+b_{i}^{2})\geq 2a_{r}b_{r}-2a_{i}b_{i}\\\
&(a_{r}^{2}+b_{r}^{2})+(a_{i}^{2}+b_{i}^{2})\geq-2a_{r}b_{r}+2a_{i}b_{i}\end{aligned}\right\\}\Rightarrow\big{|}(a_{r}^{2}+b_{r}^{2})+(a_{i}^{2}+b_{i}^{2})\big{|}\geq
2\big{|}a_{r}b_{r}-a_{i}b_{i}\big{|},$ (38)
$\displaystyle\left.\begin{aligned}
&(a_{r}^{2}+b_{i}^{2})+(a_{i}^{2}+b_{r}^{2})\geq 2a_{r}b_{i}+2a_{i}b_{r}\\\
&(a_{r}^{2}+b_{i}^{2})+(a_{i}^{2}+b_{r}^{2})\geq-2a_{r}b_{i}-2a_{i}b_{r}\end{aligned}\right\\}\Rightarrow\big{|}(a_{r}^{2}+b_{i}^{2})+(a_{i}^{2}+b_{r}^{2})\big{|}\geq
2\big{|}a_{r}b_{i}+a_{i}b_{r}\big{|}.$

Then, it is easy to see that $|a|^{2}+|b|^{2}$ can be lower bounded by

$|a|^{2}+|b|^{2}=\frac{1}{2}\big{|}(a_{r}^{2}+b_{r}^{2})+(a_{i}^{2}+b_{i}^{2})\big{|}+\frac{1}{2}\big{|}(a_{r}^{2}+b_{i}^{2})+(a_{i}^{2}+b_{r}^{2})\big{|}\geq|a_{r}b_{r}-a_{i}b_{i}|+|a_{r}b_{i}+a_{i}b_{r}|.$
(39)

Besides, since $ab=(a_{r}a_{i}-b_{r}b_{i})+j(a_{r}b_{i}+a_{i}b_{r})$ holds,
(39) can be written as
$|a|^{2}+|b|^{2}\geq\big{|}\Re(ab)\big{|}+\big{|}\Im(ab)\big{|}$. This thus
completes the proof. ∎

## References

  * [1] Q. Wu, S. Zhang, B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface-aided wireless communications: A tutorial,” _IEEE Trans. Commun._ , vol. 69, no. 5, pp. 3313–3351, May 2021.
  * [2] B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface assisted multi-user OFDMA: Channel estimation and training design,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 12, pp. 8315–8329, Dec. 2020.
  * [3] X. Tan, Z. Sun, D. Koutsonikolas, and J. M. Jornet, “Enabling indoor mobile millimeter-wave networks based on smart reflect-arrays,” in _IEEE INFOCOM_ , 2018, pp. 270–278.
  * [4] L. Subrt and P. Pechac, “Controlling propagation environments using intelligent walls,” in _EUCAP_ , 2012, pp. 1–5.
  * [5] M. M. Zhao, Q. Wu, M. J. Zhao, and R. Zhang, “Exploiting amplitude control in intelligent reflecting surface aided wireless communication with imperfect CSI,” _IEEE Trans. Commun._ , vol. 69, no. 6, pp. 4216–4231, Jun. 2021\.
  * [6] S. Gong, X. Lu, D. T. Hoang, D. Niyato, L. Shu, D. I. Kim, and Y. C. Liang, “Towards smart radio environment for wireless communications via intelligent reflecting surfaces: A comprehensive survey,” _IEEE Commun. Surveys Tuts._ , vol. 22, no. 4, pp. 2283–2314, 2020.
  * [7] S. Abeywickrama, R. Zhang, Q. Wu, and C. Yuen, “Intelligent reflecting surface: Practical phase shift model and beamforming optimization,” _IEEE Trans. Commun._ , vol. 68, no. 9, pp. 5849–5863, Jun. 2020.
  * [8] D. Mishra and H. Johansson, “Channel estimation and low-complexity beamforming design for passive intelligent surface assisted MISO wireless energy transfer,” in _IEEE ICASSP_ , 2019, pp. 4659–4663.
  * [9] T. L. Jensen and E. D. Carvalho, “An optimal channel estimation scheme for intelligent reflecting surfaces based on a minimum variance unbiased estimator,” in _IEEE ICASSP_ , May 2020, pp. 5000–5004.
  * [10] Z. Wang, L. Liu, and S. Cui, “Channel estimation for intelligent reflecting surface assisted multiuser communications: Framework, algorithms, and analysis,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 10, pp. 6607–6620, Jun. 2020.
  * [11] Y. Wei, M. M. Zhao, M. J. Zhao, and Y. Cai, “Channel estimation for IRS-aided multiuser communications with reduced error propagation,” _IEEE Trans. Wireless Commun._ , DOI:10.1109/TWC.2021.3115161, 2021.
  * [12] J. Chen, Y. Liang, H. V. Cheng, and W. Yu, “Channel estimation for reconfigurable intelligent surface aided multi-user MIMO systems,” _arXiv:1912.03619v1_ , 2019.
  * [13] A. Taha, M. Alrabeiah, and A. Alkhateeb, “Enabling large intelligent surfaces with compressive sensing and deep learning,” _IEEE Access_ , vol. 9, pp. 44 304–44 321, Mar. 2021.
  * [14] Z. He and X. Yuan, “Cascaded channel estimation for large intelligent metasurface assisted massive MIMO,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 2, pp. 210–214, Oct. 2020.
  * [15] M. M. Zhao, A. Liu, Y. Wan, and R. Zhang, “Two-timescale beamforming optimization for intelligent reflecting surface aided multiuser communication with QoS constraints,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 9, pp. 6179–6194, Sept. 2021.
  * [16] M. M. Zhao, Q. Wu, M. J. Zhao, and R. Zhang, “Intelligent reflecting surface enhanced wireless networks: Two-timescale beamforming optimization,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 1, pp. 2–17, Jan. 2021.
  * [17] C. Komninakis, C. Fragouli, A. Sayed, and R. Wesel, “Multi-input multi-output fading channel tracking and equalization using Kalman estimation,” _IEEE Trans. Signal Proces._ , vol. 50, no. 5, pp. 1065–1076, Dec. 2002.
  * [18] Z. Jellali and L. Najjar Atallah, “Fast fading channel estimation by Kalman filtering and CIR support tracking,” _IEEE Trans. Broadcast._ , vol. 63, no. 4, pp. 635–643, Dec. 2017.
  * [19] L. Lian, A. Liu, and V. K. N. Lau, “Exploiting dynamic sparsity for downlink FDD-massive MIMO channel tracking,” _IEEE Trans. Signal Proces._ , vol. 67, no. 8, pp. 2007–2021, Apr. 2019.
  * [20] P. Cai, J. Zong, X. Luo, Y. Zhou, S. Chen, and H. Qian, “Downlink channel tracking for intelligent reflecting surface-aided FDD MIMO systems,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 4, pp. 3341–3353, Apr. 2021.
  * [21] I. E. Telatar and D. N. C. Tse, “Capacity and mutual information of wideband multipath fading channels,” _IEEE Trans. Inf. Theory_ , vol. 46, no. 4, pp. 1384–1400, Jul. 2000.
  * [22] S. Borade and L. Zheng, “Writing on fading paper, dirty tape with little ink: Wideband limits for causal transmitter CSI,” _IEEE Trans. Inf. Theory_ , vol. 58, no. 8, pp. 5388–5397, Aug. 2012.
  * [23] J. Ziniel and P. Schniter, “Dynamic compressive sensing of time-varying signals via approximate message passing,” _IEEE Trans. Signal Process._ , vol. 61, no. 21, pp. 5270–5284, Aug. 2013.
  * [24] L. Lian, A. Liu, and V. K. N. Lau, “Exploiting dynamic sparsity for downlink FDD-massive MIMO channel tracking,” _IEEE Trans. Signal Process._ , vol. 67, no. 8, pp. 2007–2021, Apr. 2019.
  * [25] G. J. McLachlan and T. Krishnan, _The EM algorithm and extensions_. John Wiley & Sons, 2007, vol. 382.
  * [26] P. Zarchan and H. Musoff, _Fundamentals of Kalman Filtering: A Practical Approach_. USA: American Institute of Aeronautics and Astronautics, 2000.
  * [27] E. Ghysels and M. Marcellino, _Applied Economic Forecasting using Time Series Methods_. New York, USA: Oxford University Press., 2018.
  * [28] S. M. Kay, _Fundamentals of Statistical Signal Processing_. USA: Prentice-Hall, 1993.
  * [29] W. Jiang and H. D. Schotten, “Recurrent neural network-based frequency-domain channel prediction for wideband communications,” in _VTC2019-Spring_ , 2019, pp. 1–6.
  * [30] W. Jiang and H. D. Schotten, “Deep learning for fading channel prediction,” _IEEE OJ-COMS_ , vol. 1, pp. 320–332, Mar. 2020.
  * [31] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural Computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [32] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _Proc. Int. Conf. Learn. Repres._ , 2015.
  * [33] M. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking,” _IEEE Trans. Signal Process._ , vol. 50, no. 2, pp. 174–188, Feb. 2002.
  * [34] S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear estimation,” _Proceedings of the IEEE_ , vol. 92, no. 3, pp. 401–422, Mar. 2004.
  * [35] M. D. Springer and W. E. Thompson, “The distribution of products of beta, gamma and gaussian random variables,” _SIAM J. Appl. Math._ , vol. 18, p. 721, 1970.
  * [36] Z. Stojanac, D.Suess, and M. Kliesch, “On products of gaussian random variables,” _arXiv1711.10516_ , 2018.
  * [37] C. You, B. Zheng, and R. Zhang, “Channel estimation and passive beamforming for intelligent reflecting surface: Discrete phase shift and progressive refinement,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2604–2620, Nov. 2020.
  * [38] J. M. Steele, _The Cauchy-Schwarz Master Class: An Introduction to the Art of Mathematical Inequalities_. USA: Cambridge University Press., 2004.