# A Multi-Dimensional Matrix Pencil-Based
Channel Prediction Method for
Massive MIMO with Mobility

Weidong Li, Haifan Yin, Ziao Qin, Yandi Cao, and Mérouane Debbah W. Li, H.
Yin, Z. Qin and Y. Cao are with Huazhong University of Science and Technology,
430074 Wuhan, China (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>yandicao@hust.edu.cn).M. Debbah is with the Technology
Innovation Institute, and also with the Mohamed Bin Zayed University of
Artificial Intelligence, 9639 Masdar City, Abu Dhabi, United Arab Emirates
(e-mail: merouane.debbah@tii.ae).A part of this work [1] was presented in the
IEEE International Conference on Communications (IEEE ICC 2022).This work is
supported by the National Natural Science Foundation of China under Grant
62071191.

###### Abstract

This paper addresses the mobility problem in massive multiple-input multiple-
output systems, which leads to significant performance losses in the practical
deployment of the fifth generation mobile communication networks. We propose a
novel channel prediction method based on multi-dimensional matrix pencil
(MDMP), which estimates the path parameters by exploiting the angular-
frequency-domain and angular-time-domain structures of the wideband channel.
The MDMP method also entails a novel path pairing scheme to pair the delay and
Doppler, based on the super-resolution property of the angle estimation. Our
method is able to deal with the realistic constraint of time-varying path
delays introduced by user movements, which has not been considered so far in
the literature. We prove theoretically that in the scenario with time-varying
path delays, the prediction error converges to zero with the increasing number
of the base station (BS) antennas, providing that only two arbitrary channel
samples are known. We also derive a lower-bound of the number of the BS
antennas to achieve a satisfactory performance. Simulation results under the
industrial channel model of 3GPP demonstrate that our proposed MDMP method
approaches the performance of the stationary scenario even when the users’
velocity reaches 120 km/h and the latency of the channel state information is
as large as 16 ms.

###### Index Terms:

Massive MIMO, mobility, channel prediction, CSI delay, matrix pencil, MDMP
prediction method, channel structure.

##  I Introduction

Massive multiple-input multiple-output (MIMO) technology is playing a key role
in enhancing the spectral efficiency of the fifth generation (5G) mobile
communication systems [2]. Compared to the conventional MIMO, massive MIMO
greatly improves the capacity and reliability of mobile communication system
by deploying more antenna elements at the base station (BS) [3].

Theoretically, the performance of massive MIMO is governed by the accuracy of
the channel state information (CSI). In practice, there are several causes of
inaccurate CSI, including pilot contamination [4], imperfect CSI feedback in
frequency division duplexing (FDD) mode [2], the mobility problem [5] (or the
“curse of mobility”), etc. Although a rich body of literature has addressed
the problems of pilot contamination and CSI feedback in FDD, e.g., in [6, 7,
8, 9], the mobility problem has received relatively little attention so far,
and is still a challenging problem to be solved. Such a problem mainly results
from the mobility of user equipment (UE) and the delay of CSI [5]. The
movement of the UE makes the estimation of CSI outdated and unusable for
multi-user precoding, particularly in high-mobility scenario with large CSI
delay. Some papers have focused on the effect of CSI delay and prove that it
is harmful for the spectral efficiency performance of massive MIMO [10, 11,
12, 13].

One way to overcome this problem is by channel prediction, which has been
investigated in the literature. The work in [14] proposes a spatial-temproal
basis expansion model (ST-BEM) method to predict the downlink (DL) channel.
The authors of [15] propose a compressed sensing channel prediction method. An
efficient approximated maximum likelihood estimator is proposed in [16].
However, the aforementioned literature generally assumes the channel
parameters are time-invariant, which might be questionable in practical
mobility scenarios. The time-varying delay in Wi-Fi communication scenario is
studied in [17], where the authors propose an iterative channel estimator
basing on a two-dimensional (2-D) subspace spanned by the discrete prolate
spheroidal (DPS) sequence. The DPS sequence has less spectrum leakage than the
discrete Fourier transform (DFT) sequence. However, the DPS sequence needs to
meet two specific requirements in applications: ${\omega_{\max}}{T}\ll 1$ and
${\tau_{\max}}\Delta f\ll 1$, where ${\omega_{\max}}$ is the maximum Doppler,
${\tau_{\max}}$ denotes the maximum delay, ${T}$ is the duration of sample,
and $\Delta f$ represents the sub-carrier spacing. These requirements might
not be easily full-filled in practice.

Recently, machine-learning (ML) and neural network (NN) algorithms, such as
recurrent NN, conventional NN and deep NN, have been applied to predict the
channel through the trained network with the historical data [18, 19, 20, 21].
The authors in [21] propose a data-driven channel prediction method by
adopting complex-valued NN algorithm. However, these algorithms usually need
enormous data, which cannot be collected easily in realistic systems.
Additionally, these algorithms may perform well in static and low-mobility
scenarios, yet they might not achieve the expected performance in high-
mobility scenario due to the long training time. Moreover, the trained network
may not perform ideally in the varying environment, due to the challenge of
network generalization.

By transforming the wideband channel into angular-delay domain, the paper [5]
proposes a Prony-based angular-delay domain (PAD) prediction method and proves
that it can overcome the problem of CSI aging. However, it does not take the
effect of time-varying path delay into consideration and assumes the CSI delay
to be an integral multiple of the pilot interval. To the best of our
knowledge, the realistic effect of time-varying path delay is overlooked in
most of the existing literature. Such an effect makes the path delay difficult
to track and leads to the distortion of the traditional time-domain multipath
channel structure.

In order to break the limitations above and tackle the mobility problem, we
propose a novel matrix pencil (MP) based super-resolution channel prediction
method. The super-resolution property of our method comes from the super-
accurate parameters estimated by MP algorithm. In the literature, some papers
adopt MP algorithm to estimate channel parameters. In [22] the authors adopt
2-D MP algorithm to estimate 2-D frequencies. The work in [23] estimates the
path angle and delay in the IEEE 802.11ac context. The elevation angle of
departure (EOD), azimuth angle of departure (AOD) and delay are estimated in
[24]. Compared to other traditional super-resolution algorithms, e.g.,
multiple signal classification (MUSIC) [25] and estimation of signal
parameters via rational invariance techniques (ESPRIT) [26], MP algorithm has
some distinct advantages: it does not search and find peak value in space, and
it needs less channel samples. However, the traditional MP algorithm cannot be
directly applied to solve the mobility problem of massive MIMO, because of the
time-varying path delay and the multi-dimensional structure of the wideband
massive MIMO channel.

In order to achieve more accurate channel predictions in this paper, we extend
the traditional MP method to multi-dimensional algorithm that extracts the
EOD, AOD, time-varying delay and Doppler simultaneously. Specifically, we
first introduce a three-dimensional (3-D) MP algorithm to estimate the EOD,
AOD, and delay on the angular-frequency domain, and then put forward another
3-D MP algorithm to estimate the EOD, AOD and Doppler on the angular-time
domain. Based on the super-resolution property of the estimated EOD and AOD,
we propose a pairing algorithm to pair the path delay and Doppler. The future
CSI is reconstructed with the estimated parameters. To the best of our
knowledge, our proposed multi-dimensional matrix pencil (MDMP) method is the
first attempt to estimate the EODs, AODs, Doppler and the time-varying delays
of multiple paths simultaneously.

The contributions of this paper are summarized as follows:

  * •

We propose a super-resolution MDMP prediction method to address the mobility
problem with time-varying path delays, which estimates the multipath EODs,
AODs, delays and Doppler simultaneously through angular-time-domain and
angular-frequency-domain structures of the channel. Compared to the
traditional methods, significant gains of our proposed method are confirmed in
simulations.

  * •

We prove that the prediction error of the MDMP method converges to zero in
mobility scenario with arbitrary CSI delay, when the number of the BS antennas
and the bandwidth are large enough. Our proposed method breaks the limitation
of the PAD method whose CSI delay is an integral multiple of the pilot
interval.

  * •

We also prove that the prediction error converges to zero providing that only
two arbitrary samples are known, when the number of the BS antennas is large
enough. Our proposed MDMP method does not necessarily require the samples to
be neighboring ones, yet it is a common assumption in most existing works.

  * •

We derive a lower-bound of the number of the BS antennas in a wideband channel
to give a satisfactory performance of the MDMP method. The lower-bound is
correlated with the number of samples. As the number of samples increases, the
lower-bound decreases.

This paper is organized as follows: Sec. II introduces the channel model. In
Sec. III, MDMP prediction method is proposed. The performance of MDMP method
is analyzed in Sec. IV. The simulation results are shown in Sec. V, and Sec.
VI concludes the paper.

Notations: We use boldface to represent vectors and matrices, where
${{\bf{0}}_{{M_{1}}\times{N_{1}}}}$, ${{\bf{I}}_{{M_{2}}}}$,
${\Upsilon_{{M_{3}}}}$ and ${{\bf{1}}_{{M_{4}}\times{N_{5}}}}$ denote
${M_{1}}\times{N_{1}}$ zero matrix, ${M_{2}}\times{M_{2}}$ identity matrix,
${M_{3}}\times{M_{3}}$ anti-identity matrix and ${M_{4}}\times{N_{5}}$ all-
ones matrix. Let ${({\bf{X}})^{T}}$, ${({\bf{X}})^{*}}$, ${({\bf{X}})^{H}}$,
${({\bf{X}})^{-1}}$ and ${({\bf{X}})^{\dagger}}$ denote the transpose,
conjugate, conjugate transpose, inverse and Moore-Penrose pseudoinverse of the
matrix ${\bf{X}}$, respectively. ${\rm{card}}(\cdot)$ denotes the number of
the elements in a set. $r\\{\cdot\\}$ denotes the rank of a matrix.
${\left\|\cdot\right\|_{F}}$ stands for the Frobenius norm. ${\mathop{\rm
Re}\nolimits}\left\\{\cdot\right\\}$ and ${\mathop{\rm
Im}\nolimits}\left\\{\cdot\right\\}$ take the real and imaginary components of
a complex number. $E\\{\cdot\\}$ is the expectation operation.
${\rm{diag}}\\{\cdot\\}$ denotes the diagonal operation of a matrix.
$[{\bf{X}}:{\bf{Y}}]$ is the extending matrix of ${\bf{X}}$ and ${\bf{Y}}$.
${\bf{X}}\odot{\bf{Y}}$ and ${\bf{X}}\otimes{\bf{Y}}$ denote the Hadamard
product and Kronecker product of ${\bf{X}}$ and ${\bf{Y}}$.

##  II Channel Model

We consider a wideband time division duplexed (TDD) massive MIMO system, where
a BS serves multiple UEs. In such a system, the BS estimates the CSI from the
uplink (UL) pilot, and the DL CSI is acquired based on channel reciprocity.

Figure 1: The multipath channel between the BS and the UE.

Fig. 1 illustrates the channel between the BS and the UE. In most 5G
commercial systems, the BS is equipped with a uniform planar array (UPA),
which consists of $N_{h}^{{\rm{BS}}}$ columns and $N_{v}^{{\rm{BS}}}$ rows of
antenna elements. The numbers of the BS antennas and the UE antennas are
respectively $N_{t}={N_{h}^{{\rm{BS}}}}{N_{v}^{{\rm{BS}}}}$ and
$N_{r}^{{\rm{UE}}}$. There are ${M_{\rm{cluster}}}$ different clusters, which
are denoted by ${S_{i}}$, $i=1,2,\cdots,{M_{\rm{cluster}}}$. Each cluster
contains many scattering rays. The UL and DL channels share the same
bandwidth, which is composed of ${N_{f}}$ subcarriers with spacing $\Delta f$.

The DL MIMO channel between the BS and the UE is expressed as
${\bf{H}}(t,f)={\left[{{h_{u,{s_{h}},{s_{v}}}}(t,f)}\right]_{N_{r}^{{\rm{UE}}}\times
N_{h}^{{\rm{BS}}}N_{v}^{{\rm{BS}}}}}$, where ${h_{u,{s_{h}},{s_{v}}}}(t,f)$
denotes the channel frequency response between the ${s_{h}}$-th column and the
${s_{v}}$-th row of the BS antenna array and the $u$-th antenna of the UE, and
is modeled as [27]

$\begin{array}[]{l}\\!\\!\\!{h_{u,{s_{h}},{s_{v}}}}(t,f)\\\
\\!\\!\\!=\sum\limits_{p=1}^{P}{}{\beta_{p}}{e^{\frac{{j2\pi({\bf{\hat{r}}}{{{}_{p}^{{\rm{rx}}}})^{T}}{\bf{\bar{d}}}_{u}^{{\rm{rx}}}}}{{{\lambda_{0}}}}}}{e^{\frac{{j2\pi({\bf{\hat{r}}}{{{}_{p}^{{\rm{tx}}}})^{T}}{\bf{\bar{d}}}_{{s_{h}},{s_{v}}}^{{\rm{tx}}}}}{{{\lambda_{0}}}}}}{e^{j2\pi{\omega_{p}}t}}{e^{-j2\pi
f{\tau_{p}}(t)}}\end{array},\\!\\!\\!$ (1)

where $P$ denotes the number of paths, and ${\beta_{p}}$ is the complex
amplitude of the $p$-th path. Also, ${\lambda_{0}}=\frac{{\rm{c}}}{{{f_{c}}}}$
is the wavelength, where ${f_{c}}$ is the central carrier frequency and
${\rm{c}}$ is the speed of light. Furthermore,
${\bf{\hat{r}}}_{p}^{{\rm{rx}}}$ and ${\bf{\hat{r}}}_{p}^{{\rm{tx}}}$ denote
the spherical unit vectors of the UE and the BS, and are expressed as

$\begin{array}[]{l}{\bf{\hat{r}}}_{p}^{{\rm{rx}}}=\left[\begin{array}[]{l}\cos{\theta_{p,{\rm{EOA}}}}\cos{\phi_{p,{\rm{AOA}}}}\\\
\cos{\theta_{p,{\rm{EOA}}}}\sin{\phi_{p,{\rm{AOA}}}}\\\ \ \ \ \ \ \ \
\sin{\theta_{p,{\rm{EOA}}}}\end{array}\right],\end{array}\\!$ (2)
$\begin{array}[]{l}{\bf{\hat{r}}}_{p}^{{\rm{tx}}}=\left[\begin{array}[]{l}\cos{\theta_{p,{\rm{EOD}}}}\cos{\phi_{p,{\rm{AOD}}}}\\\
\cos{\theta_{p,{\rm{EOD}}}}\sin{\phi_{p,{\rm{AOD}}}}\\\ \ \ \ \ \ \ \
\sin{\theta_{p,{\rm{EOD}}}}\end{array}\right],\end{array}\\!$ (3)

where ${\theta_{p,{\rm{EOA}}}}$, ${\phi_{p,{\rm{AOA}}}}$,
${\theta_{p,{\rm{EOD}}}}$ and ${\phi_{p,{\rm{AOD}}}}$ are the elevation angle
of arrival (EOA), azimuth angle of arrival (AOA), EOD and AOD, respectively.
Additionally, ${\bf{\bar{d}}}_{{s_{h}}{s_{v}}}^{{\rm{tx}}}$ and
${\bf{\bar{d}}}_{u}^{{\rm{rx}}}$ are the location vectors of the BS and the UE
antennas:

${\bf{\bar{d}}}_{{s_{h}},{s_{v}}}^{{\rm{tx}}}={\left[{\begin{array}[]{*{20}{c}}{0,}&{d_{h}^{{\rm{tx}}}({s_{h}}-1),}&{d_{v}^{{\rm{tx}}}({s_{v}}-1)}\end{array}}\right]^{T}},\\!$
(4)

where $d_{h}^{{\rm{tx}}}$ is the horizontal antenna spacing, and
$d_{v}^{{\rm{tx}}}$ is the vertical antenna spacing. Moreover,
${\omega_{p}}=\frac{{({\bf{\hat{r}}}{{{}_{p}^{{\rm{rx}}}})^{T}}{\bf{v}}}}{{{\lambda_{0}}}}$
represents the Doppler, and ${\bf{v}}$ is the velocity vector of the UE:

${\bf{v}}=v{\left[{\begin{array}[]{*{20}{c}}{\cos{\theta_{v}}\cos{\phi_{v}},}&{\cos{\theta_{v}}\sin{\phi_{v}},}&{\sin{\theta_{v}}}\end{array}}\right]^{T}},\\!$
(5)

where $v$ is the speed, ${\theta_{v}}$ is the elevation angle of velocity, and
${\phi_{v}}$ is the azimuth angle of velocity. The $p$-th path delay
${\tau_{p}}(t)$ is time-varying and modeled as [27]

${\tau_{p}}(t)={\tau_{p,0}}+{k_{{\tau_{p}}}}t={\tau_{p,0}}-\frac{{({\bf{\hat{r}}}{{{}_{p}^{{\rm{rx}}}})^{T}}{\bf{v}}}}{{\rm{c}}}t={\tau_{p,0}}-\frac{{{\omega_{p}}}}{{{f_{c}}}}t,\\!$
(6)

where ${\tau_{p,0}}$ is the initial value, and ${k_{{\tau_{p}}}}$ denotes the
changing rate of the path delay. Notice that, besides the time domain, Doppler
also has an effect on the frequency domain. The 3-D steering vector
${{\bf{a}}^{{\rm{tx}}}}({\theta_{p}},{\phi_{p}})$ is expressed as

$\begin{array}[]{l}{{\bf{a}}^{{\rm{tx}}}}({\theta_{p}},{\phi_{p}})={\left[{\begin{array}[]{*{20}{c}}{{(a_{h}^{{\rm{tx}}})}^{0},}&{{(a_{h}^{{\rm{tx}}})}^{1},}&{\cdots,}&{{(a_{h}^{{\rm{tx}}})}^{(N_{h}^{{\rm{BS}}}-1)}}\end{array}}\right]^{T}}\\\
\ \ \ \ \ \ \ \ \ \ \ \ \ \
\otimes{\left[{\begin{array}[]{*{20}{c}}{{(a_{v}^{{\rm{tx}}})}^{0},}&{{(a_{v}^{{\rm{tx}}})}^{1},}&{\cdots,}&{{(a_{v}^{{\rm{tx}}})}^{(N_{v}^{{\rm{BS}}}-1)}}\end{array}}\right]^{T}},\end{array}\\!$
(7)

where ${\theta_{p}}$ and ${\phi_{p}}$ are ${\theta_{p,{\rm{EOD}}}}$ and
${\phi_{p,{\rm{AOD}}}}$ for simplicity. Furthermore, $a_{h}^{{\rm{tx}}}$ and
$a_{v}^{{\rm{tx}}}$ denote the spatial signatures in the directions of
$\cos({\theta_{p}})\sin({\phi_{p}})$ and $\sin({\theta_{p}})$, which are
expressed as

$a_{h}^{{\rm{tx}}}={\rm{exp}}({\frac{{j2\pi
d_{h}^{{\rm{tx}}}\cos({\theta_{p}})\sin({\phi_{p}})}}{{{\lambda_{0}}}}}),\\!$
(8) $a_{v}^{{\rm{tx}}}={\rm{exp}}({\frac{{j2\pi
d_{v}^{{\rm{tx}}}\sin({\theta_{p}})}}{{{\lambda_{0}}}}}).\\!$ (9)

Let ${{\bf{h}}_{u}}(t,f)$ denote the channel between all BS antennas and the
$u$-th UE antenna at time $t$ and frequency $f$:

$\begin{array}[]{l}\\!\\!\\!\\!\\!\\!\\!\\!\\!{{\bf{h}}_{u}}(t,f)\\!=\\!\\!\sum\limits_{p=1}^{P}{}\\!\\!{\beta_{p}}{e^{\frac{{j2\pi({\bf{\hat{r}}}{{{}_{p}^{{\rm{rx}}}}})^{T}{\bf{\bar{d}}}_{u}^{{\rm{rx}}}}}{{{\lambda_{0}}}}}}\\!\\!{e^{j2\pi{\omega_{p}}t}}{e^{-j2\pi
f{\tau_{p}}(t)}}{{\bf{a}}^{{\rm{tx}}}}({\theta_{p}},{\phi_{p}}).\\!\\!\\!\\!\\!\end{array}\\!\\!\\!$
(10)

In the time domain, the channels at all subcarriers are expressed as

$\begin{array}[]{l}\\!\\!\\!\\!\\!\\!\\!\\!\\!{{\bf{H}}_{u}}(t)=[{\begin{array}[]{*{20}{c}}{{\bf{h}}_{u}^{T}(t,{f_{1}}),}&{{\bf{h}}_{u}^{T}(t,{f_{2}}),}&{\cdots,}&{{\bf{h}}_{u}^{T}(t,{f_{{N_{f}}}})}\end{array}}],\\!\\!\\!\end{array}\\!\\!\\!$
(11)

where ${f_{{n_{f}}}}$ is the frequency of the ${n_{f}}$-th subcarrier:

${f_{{n_{f}}}}={f_{1}}+({n_{f}}-1)\Delta
f,{\rm{}}{n_{f}}=1,2,\cdots{N_{f}}.\\!$ (12)

Based on the angular-frequency-domain channel structure, ${{\bf{H}}_{u}}(t)$
is rewritten as

${{\bf{H}}_{u}}(t)={\bf{A}}_{u}^{{\rm{tx}}}{{\bf{C}}_{u}}{{\bf{B}}_{u}},\\!$
(13)

where ${\bf{A}}_{u}^{{\rm{tx}}}$ contains the 3-D steering vectors of all
paths:

$\begin{array}[]{l}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!{\bf{A}}_{u}^{{\rm{tx}}}=[{\begin{array}[]{*{20}{c}}{{{\bf{a}}^{{\rm{tx}}}}({\theta_{1}},\\!{\phi_{1}}),}&{\cdots,}&{{{\bf{a}}^{{\rm{tx}}}}({\theta_{P}},\\!{\phi_{P}})}\end{array}}].\\!\\!\\!\\!\\!\\!\\!\end{array}\\!\\!\\!\\!\\!\\!$
(14)

The frequency-domain matrix ${{\bf{B}}_{u}}$ is expressed as

${{\bf{B}}_{u}}=\left[{\begin{array}[]{*{20}{c}}{{\bf{b}}({f_{1}}),}&{{\bf{b}}({f_{2}}),}&{\cdots,}&{{\bf{b}}({f_{{N_{f}}}})}\end{array}}\right],\\!$
(15)

where ${\bf{b}}({f_{{n_{f}}}})$ denotes the delay response vector at the
${n_{f}}$-th subcarrier:

${\bf{b}}({f_{{n_{f}}}})={\left[{\begin{array}[]{*{20}{c}}{e^{-j2\pi{f_{{n_{f}}}}{\tau_{1}}(t)},}&{\cdots,}&{{e^{-j2\pi{f_{{n_{f}}}}{\tau_{P}}(t)}}}\end{array}}\right]^{T}}.\\!\\!\\!\\!$
(16)

The diagonal matrix ${{\bf{C}}_{u}}$ is expressed as

$\begin{array}[]{l}\\!\\!\\!\\!\\!\\!\\!\\!{{\bf{C}}_{u}}\\\
\\!\\!\\!\\!\\!\\!\\!=\\!{\rm{diag}}\\!\\!\left\\{\\!{{\beta_{1}}{e^{\frac{{j2\pi({\bf{\hat{r}}}{{{}_{1}^{{\rm{rx}}}}})^{T}{\bf{\bar{d}}}_{u}^{{\rm{rx}}}}}{{{\lambda_{0}}}}}}\\!{e^{j2\pi{\omega_{1}}t}},\\!\cdots\\!,{\beta_{P}}{e^{\frac{{j2\pi({\bf{\hat{r}}}{{{}_{P}^{{\rm{rx}}}}})^{T}{\bf{\bar{d}}}_{u}^{{\rm{rx}}}}}{{{\lambda_{0}}}}}}\\!{e^{j2\pi{\omega_{P}}t}}}\\!\right\\}\\!\\!.\\!\\!\end{array}\\!\\!\\!\\!\\!\\!$
(17)

According to the angular-time-domain channel structure, ${{\bf{H}}_{u}}(t)$ is
also rewritten as

${{\bf{H}}_{u}}(t)={\bf{A}}_{u}^{{\rm{tx}}}{{\bf{D}}_{u}}({{\bf{E}}_{u}}\odot{{\bf{F}}_{u}}),\\!\\!$
(18)

where ${{\bf{D}}_{u}}$ is a diagonal matrix:

${{\bf{D}}_{u}}={\rm{diag}}\left\\{{\begin{array}[]{*{20}{c}}{{\beta_{1}}{e^{\frac{{j2\pi({\bf{\hat{r}}}{{{}_{1}^{{\rm{rx}}}}})^{T}{\bf{\bar{d}}}_{u}^{{\rm{rx}}}}}{{{\lambda_{0}}}}}},}&{\cdots,}&{{\beta_{P}}{e^{\frac{{j2\pi({\bf{\hat{r}}}{{{}_{P}^{{\rm{rx}}}}})^{T}{\bf{\bar{d}}}_{u}^{{\rm{rx}}}}}{{{\lambda_{0}}}}}}}\end{array}}\right\\}.\\!\\!$
(19)

The matrix ${{\bf{E}}_{u}}$ is composed of the initial delay response vectors
at all subcarriers:

${{\bf{E}}_{u}}=\left[{\begin{array}[]{*{20}{c}}{{\bf{e}}({f_{1}}),}&{{\bf{e}}({f_{2}}),}&{\cdots,}&{{\bf{e}}({f_{{N_{f}}}})}\end{array}}\right],\\!\\!$
(20)

where ${\bf{e}}({f_{{n_{f}}}})$ is for the ${n_{f}}$-th subcarrier:

${\bf{e}}({f_{{n_{f}}}})={\left[{\begin{array}[]{*{20}{c}}{e^{-j2\pi{f_{{n_{f}}}}{\tau_{1,0}}},}&{\cdots,}&{e^{-j2\pi{f_{{n_{f}}}}{\tau_{P,0}}}}\end{array}}\right]^{T}}.\\!\\!$
(21)

The time-domain matrix ${{\bf{F}}_{u}}$ is expressed as

${{\bf{F}}_{u}}={\left[{\begin{array}[]{*{20}{c}}{{\bf{f}}({f_{1}}),}&{{\bf{f}}({f_{2}}),}&{\cdots,}&{{\bf{f}}({f_{{N_{f}}}})}\end{array}}\right]^{T}},\\!\\!$
(22)

where ${\bf{f}}({f_{{n_{f}}}}),1\leq n_{f}\leq N_{f}$, is defined as:

$\\!{\bf{f}}({f_{{n_{f}}}})={\left[{\begin{array}[]{*{20}{c}}{e^{j2\pi({\omega_{1}}-{f_{{n_{f}}}}{k_{{\tau_{1}}}})t},}&{\cdots,}&{e^{j2\pi({\omega_{P}}-{f_{{n_{f}}}}{k_{{\tau_{P}}}})t}}\end{array}}\right]^{T}}\\!\\!\\!.\\!\\!\\!\\!\\!\\!$
(23)

Notice that Eq. (18) is different from the result in [5], as we simultaneously
consider the effect of Doppler on time domain and frequency domain. The
multipath delays, angles and Doppler will be estimated jointly by our proposed
MDMP method in the next section.

##  III The Proposed MDMP Prediction Method

In this section, we introduce our proposed MDMP super-resolution method for
channel prediction. In general, our approach first estimates the EOD, AOD,
delay, and Doppler information by exploiting the angular-frequency-domain and
angular-time-domain structures of the channel, and then performs a path
pairing procedure to determine the corresponding parameters of the paths.
Finally the BS predicts the CSI by reconstructing the future channel with the
estimated parameters. To ease the exposition, we briefly introduce here the
derivations in Sec. III-A and Sec. III-B. Since Sec. III-A and Sec. III-B have
similar derivations, we take Sec. III-A for example. The derivation mainly
contains five steps, i.e., generating the 3-D MP matrix from the channels,
transforming the 3-D complex MP matrix to a real matrix, determining the
number of paths, estimating the path delays, and estimating the EODs and the
AODs.

###  III-A The 3-D Angle-Delay Estimation

We first generate a 3-D MP matrix. For ease of exposition, we start with the
one-dimensional (1-D) setting of the matrix pencil method, and then move on to
the 2-D and 3-D cases. More specifically, we first generate a 1-D MP matrix
only containing the information of $\cos{\theta_{p}}\sin{\phi_{p}}$, by
windowing the BS antenna panel in the horizontal direction. Then, we generate
a 2-D MP matrix by adding the second window in the vertical direction of the
BS antenna panel, which also contains the information of $\sin{\theta_{p}}$.
Finally, based on the 2-D MP matrix, we generate a 3-D MP matrix by adding the
third window in the frequency domain, which contains the information of path
delay. More details will be shown below.

Based on the angular-frequency-domain channel structure, a 1-D MP matrix is
generated by sliding a horizontal window within the BS antenna panel. Such a
MP matrix
${{\bf{G}}_{u,r}}(t,{n_{f}})\in{{\mathbb{C}}^{L\times(N_{h}^{{\rm{BS}}}-L+1)}}$
is defined as

${{\bf{G}}_{u,r}}(t,{n_{f}})\\!=\\!\\!\left[\\!\\!\begin{array}[]{l}{h_{u,1,r}}(t,{n_{f}}),\cdots\\!,{h_{u,(N_{h}^{{\rm{BS}}}\\!-\\!L\\!+\\!1),r}}(t,{n_{f}})\\\
\ \ \ \ \ \ \vdots\ \ \ \ \ \ \ \ \ \ddots\ \ \ \ \ \ \ \ \vdots\\\
{h_{u,L,r}}(t,{n_{f}}),\cdots,\ \
{h_{u,N_{h}^{{\rm{BS}}},r}}(t,{n_{f}})\end{array}\\!\\!\\!\right]\\!\\!,\\!\\!\\!\\!$
(24)

which is composed of the channels between the $u$-th UE antenna and all
horizontal antennas in the $r$-th row of the BS antenna panel at the
${n_{f}}$-th subcarrier, and $L$ is the pencil size that satisfies
$P<L<N_{h}^{{\rm{BS}}}-P+1$. Based on Eq. (13), the 1-D MP matrix is rewritten
as

${{\bf{G}}_{u,r}}(t,{n_{f}})={{\bf{X}}_{\theta,\phi,1}}{{\bf{Z}}_{\theta}^{r}}{\bf{Y}}{{\bf{Z}}_{\tau}^{{n_{f}}-1}}{{\bf{X}}_{\theta,\phi,2}},\\!$
(25)

where the $p$-th column of
${{\bf{X}}_{\theta,\phi,1}}\in{{\mathbb{C}}^{L\times P}}$ and the $p$-th row
of
${{\bf{X}}_{\theta,\phi,2}}\in{{\mathbb{C}}^{P\times(N_{h}^{{\rm{BS}}}-L+1)}}$,
are the steering vectors of two subsets of antenna array with size $1\times L$
and $1\times(N_{h}^{{\rm{BS}}}-L+1)$ in the BS antenna panel, respectively.
The two matrices ${{\bf{X}}_{\theta,\phi,1}}$ and ${{\bf{X}}_{\theta,\phi,2}}$
are defined as

${{\bf{X}}_{\theta,\phi,1}}=\left[{\begin{array}[]{*{20}{c}}{{\bf{a}}_{h,L-1}^{{\rm{tx}}}({\theta_{1}},{\phi_{1}}),}&{\cdots,}&{{\bf{a}}_{h,L-1}^{{\rm{tx}}}({\theta_{P}}\\!,\\!{\phi_{P}})}\end{array}}\right],\\!\\!\\!\\!$
(26)
$\\!\\!\\!\\!\\!\\!{{\bf{X}}_{\theta,\phi,2}}={\left[{\begin{array}[]{*{20}{c}}{{\bf{a}}_{h,N_{h}^{{\rm{BS}}}-L}^{{\rm{tx}}}({\theta_{1}},{\phi_{1}}),}&{\cdots,}&{{\bf{a}}_{h,N_{h}^{{\rm{BS}}}-L}^{{\rm{tx}}}({\theta_{P}},{\phi_{P}})}\end{array}}\\!\\!\right]^{T}}\\!\\!,\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
(27)

where
${\bf{a}}_{h,{l_{1}}}^{{\rm{tx}}}({\theta_{p}},{\phi_{p}})={\left[{{{(a_{h}^{{\rm{tx}}})}^{0}},{{(a_{h}^{{\rm{tx}}})}^{1}},\cdots,{{(a_{h}^{{\rm{tx}}})}^{{l_{1}}}}}\right]^{T}}$,
${l_{1}}\in\\{L-1,N_{h}^{{\rm{BS}}}-L\\}$. The matrix ${{\bf{Z}}_{\tau}}$
contains the phase differences between two neighboring subcarriers for the
paths:

${{\bf{Z}}_{\tau}}={\rm{diag}}\left\\{{\begin{array}[]{*{20}{c}}{e^{-j2\pi\Delta
f{\tau_{1}}(t)},}&{\cdots,}&{e^{-j2\pi\Delta
f{\tau_{P}}(t)}}\end{array}}\right\\}.\\!$ (28)

Likewise, ${{\bf{Z}}_{\theta}}$ denotes the matrix of the phase differences
between two neighboring vertical BS antennas for the paths:

${{\bf{Z}}_{\theta}}={\rm{diag}}\left\\{{\begin{array}[]{*{20}{c}}{e^{\frac{{j2\pi
d_{v}^{{\rm{tx}}}\sin({\theta_{1}})}}{{{\lambda_{0}}}}},}&{\cdots,}&{e^{\frac{{j2\pi
d_{v}^{{\rm{tx}}}\sin({\theta_{P}})}}{{{\lambda_{0}}}}}}\end{array}}\right\\}.\\!$
(29)

The matrix ${\bf{Y}}$ is defined as

${\bf{Y}}={\rm{diag}}\left\\{{{{\bf{C}}_{u}}{\bf{b}}({f_{1}})}\right\\},\\!$
(30)

where ${\bf{b}}({f_{1}})$ and ${{\bf{C}}_{u}}$ are defined in Eq. (16) and Eq.
(17).

By adding a new window along the vertical direction of the BS antenna panel, a
2-D MP matrix
${{\bf{G}}_{u}}(t,{n_{f}})\in{{\mathbb{C}}^{LR\times(N_{h}^{{\rm{BS}}}-L+1)(N_{v}^{{\rm{BS}}}-R+1)}}$
is introduced

${{\bf{G}}_{u}}(t,{n_{f}})\\!=\\!\\!\\!\left[\\!\\!\begin{array}[]{l}{{\bf{G}}_{u,1}}(t,{n_{f}}),\cdots\\!,\
{{\bf{G}}_{u,(N_{v}^{{\rm{BS}}}\\!-\\!R+1)}}(t,{n_{f}})\\\ \ \ \ \ \vdots\ \ \
\ \ \ \ \ \ \ \ddots\ \ \ \ \ \ \ \ \ \ \ \vdots\\\
{{\bf{G}}_{u,R}}(t,{n_{f}}),\cdots,\ \ \ \
{{\bf{G}}_{u,N_{v}^{{\rm{BS}}}}}(t,{n_{f}})\end{array}\\!\right]\\!,\\!\\!\\!$
(31)

which is composed of the channels between the $u$-th UE antenna and all BS
antennas at the ${n_{f}}$-th subcarrier, and contains the angular information
of the paths. The corresponding pencil size is denoted by $R$ that satisfies
$P<R<N_{v}^{{\rm{BS}}}-P+1$. Substituting Eq. (25) into Eq. (31),
${{\bf{G}}_{u}}(t,{n_{f}})$ is rewritten as

${{\bf{G}}_{u}}(t,{n_{f}})={{\check{\bf{E}}}_{1}}{\bf{Y}}{{\bf{Z}}_{\tau}^{{n_{f}}-1}}{{\check{\bf{F}}}_{1}},\\!$
(32)

where the $p$-th column of ${{\check{\bf{E}}}_{1}}\in{{\mathbb{C}}^{LR\times
P}}$ and the $p$-th row of
${{\check{\bf{F}}}_{1}}\in{{\mathbb{C}}^{P\times(N_{h}^{{\rm{BS}}}-L+1)(N_{v}^{{\rm{BS}}}-R+1)}}$
are the 3-D steering vectors of two subsets of antenna array with size
$R\times L$ and $(N_{v}^{{\rm{BS}}}-R+1)\times(N_{h}^{{\rm{BS}}}-L+1)$ on the
BS antenna panel, respectively. They are defined as

${{\check{\bf{E}}}_{1}}={\left[{\begin{array}[]{*{20}{c}}{({{\bf{X}}_{\theta,\phi,1}}{{\bf{Z}}_{\theta}^{0}})^{T},}&{\cdots,}&{({{\bf{X}}_{\theta,\phi,1}}{{\bf{Z}}_{\theta}^{R-1}})^{T}}\end{array}}\right]^{T}},\\!$
(33)
${{\check{\bf{F}}}_{1}}=\left[{\begin{array}[]{*{20}{c}}{{{\bf{Z}}_{\theta}^{0}}{{\bf{X}}_{\theta,\phi,2}},}&{\cdots,}&{{{\bf{Z}}_{\theta}^{N_{v}^{{\rm{BS}}}-R}}{{\bf{X}}_{\theta,\phi,2}}}\end{array}}\right].\\!$
(34)

Adding the third window in the frequency domain, a 3-D MP matrix
${{\bf{G}}_{u}}(t)\in{{\mathbb{C}}^{LRK\times(N_{h}^{{\rm{BS}}}-L+1)(N_{v}^{{\rm{BS}}}-R+1)({N_{f}}-K+1)}}$
is generated from ${{\bf{G}}_{u}}(t,{n_{f}})$ and expressed as

${{\bf{G}}_{u}}(t)\\!=\\!\\!\left[\\!\\!\begin{array}[]{l}{{\bf{G}}_{u}}(t,1),\
\cdots,{{\bf{G}}_{u}}(t,{N_{f}}\\!-\\!K\\!+\\!1)\\\ \ \ \ \ \ \vdots\ \ \ \ \
\ \ddots\ \ \ \ \ \ \ \ \ \ \ \vdots\\\ {{\bf{G}}_{u}}(t,K),\cdots,\ \ \ \
{{\bf{G}}_{u}}(t,{N_{f}})\end{array}\right]\\!\\!,\\!\\!\\!\\!\\!$ (35)

where $K$ is the pencil size that satisfies $P<K<{N_{f}}-P+1$. The matrix
${{\bf{G}}_{u}}(t)$ is composed of the channels between the $u$-th UE antenna
and all BS antennas at all subcarriers, and also contains the angular and
delay information of the paths. For notational simplicity, we define
$\mu_{1}=LRK$ and
$\mu_{2}=(N_{h}^{{\rm{BS}}}-L+1)(N_{v}^{{\rm{BS}}}-R+1)({N_{f}}-K+1)$.
According to the channel angular-frequency-domain structure in Eq. (13),
${{\bf{G}}_{u}}(t)$ is rewritten as

${{\bf{G}}_{u}}(t)={{\check{\bf{E}}}_{2}}{\bf{Y}}{{\check{\bf{F}}}_{2}},\\!$
(36)

where ${{\check{\bf{E}}}_{2}}\in{{\mathbb{C}}^{\mu_{1}\times P}}$ and
${{\check{\bf{F}}}_{2}}\in{{\mathbb{C}}^{P\times\mu_{2}}}$ are defined as

${{\check{\bf{E}}}_{2}}={\left[{\begin{array}[]{*{20}{c}}{({{\check{\bf{E}}}_{1}}{{\bf{Z}}_{\tau}^{0}})^{T},}&{\cdots,}&{({{\check{\bf{E}}}_{1}}{{\bf{Z}}_{\tau}^{K-1}})^{T}}\end{array}}\right]^{T}},\\!$
(37)
${{\check{\bf{F}}}_{2}}=\left[{\begin{array}[]{*{20}{c}}{{{\bf{Z}}_{\tau}^{0}}{{\check{\bf{F}}}_{1}},}&{\cdots,}&{{{\bf{Z}}_{\tau}^{{N_{f}}-K}}{{\check{\bf{F}}}_{1}}}\end{array}}\right].\\!$
(38)

The physical meaning of
${{\check{\bf{E}}}_{1}}{{\bf{Z}}_{\tau}^{{k_{1}}}},{k_{1}}=0,\cdots,K-1$ is
the EOD-AOD steering vectors of a subset of antenna array with size $R\times
L$ on the BS antenna panel at the $({k_{1}}+1)$-th subcarrier of the paths.
Likewise,
${{\bf{Z}}_{\tau}^{{k_{2}}}}{{\check{\bf{F}}}_{1}},{k_{2}}=0,\cdots,{N_{f}}-K$
is the EOD-AOD steering vectors of a subset of antenna array with size
$(N_{v}^{{\rm{BS}}}-R+1)\times(N_{h}^{{\rm{BS}}}-L+1)$ on the BS antenna panel
at the $({k_{2}}+1)$-th subcarrier for the paths. The 3-D MP matrix
${{\bf{G}}_{u}}(t)$ contains the information of multipaths, if its rank
satisfies $r({{\bf{G}}_{u}}(t))=P$. In this case, the three pencil sizes yield

$\left\\{\begin{array}[]{l}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ LR(K-1)\geq
P\\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ LK(R-1)\geq P\\\ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ RK(L-1)\geq P\\\
(N_{h}^{{\rm{BS}}}-L+1)(N_{v}^{{\rm{BS}}}-R+1)({N_{f}}-K+1)\geq
P\end{array}\\!\\!\\!\\!\right..\\!\\!\\!$ (39)

Then, we transform the 3-D MP matrix ${{\bf{G}}_{u}}(t)$ to a real matrix
${{\bf{G}}_{\rm{re}}}(t)\in{{\mathbb{C}}^{{\mu_{1}}\times 2{\mu_{2}}}}$ by the
unitary matrix pencil (UMP) method [28], which reduces the computational
complexity without losing accuracy. The real matrix ${{\bf{G}}_{\rm{re}}}(t)$
is transformed by

${{\bf{G}}_{\rm{re}}}(t)={\bf{Q}}_{\mu_{1}}^{H}{{\bf{G}}_{\rm{ex}}}(t){{\bf{Q}}_{2{\mu_{2}}}},\\!\\!\\!$
(40)

where

${{\bf{G}}_{\rm{ex}}}(t)=\left[{{{\bf{G}}_{u}}(t):{\Upsilon_{\mu_{1}}}{{\bf{G}}_{u}^{*}(t)}{\Upsilon_{{\mu_{2}}}}}\right].\\!\\!\\!\\!\\!$
(41)

We define ${{\bf{Q}}_{\mu_{1}}}$ and ${{\bf{Q}}_{2{\mu_{2}}}}$ as two unitary
matrices depending on the size of $\mu_{1}$ and $2{\mu_{2}}$. Taking
${{\bf{Q}}_{\mu_{1}}}$ for example, if $\mu_{1}$ is even,

${{\bf{Q}}_{\mu_{1}}}=\frac{1}{{\sqrt{2}}}\left[\begin{array}[]{l}{{\bf{I}}_{\frac{{\mu_{1}}}{2}}}\
\ \ \ \ j{{\bf{I}}_{\frac{{\mu_{1}}}{2}}}\\\ {\Upsilon_{\frac{{\mu_{1}}}{2}}}\
-j{\Upsilon_{\frac{{\mu_{1}}}{2}}}\end{array}\right],\\!\\!\\!\\!\\!$ (42)

and if it is odd,

${{\bf{Q}}_{\mu_{1}}}=\frac{1}{{\sqrt{2}}}\left[\begin{array}[]{l}\
{{\bf{I}}_{\frac{{{\mu_{1}}-1}}{2}}}\ \ \ \
{{\bf{0}}_{(\frac{{{\mu_{1}}-1}}{2})\times 1}}\ \ \
j{{\bf{I}}_{\frac{{{\mu_{1}}-1}}{2}}}\\\
{{\bf{0}}_{1\times(\frac{{{\mu_{1}}-1}}{2})}}\ \ {\sqrt{2}}\ \ \ \ \ \ \ \
{{\bf{0}}_{1\times(\frac{{{\mu_{1}}-1}}{2})}}\\\ \
{\Upsilon_{\frac{{{\mu_{1}}-1}}{2}}}\ \ \
{{\bf{0}}_{(\frac{{{\mu_{1}}-1}}{2})\times 1}}\
-j{\Upsilon_{\frac{{{\mu_{1}}-1}}{2}}}\end{array}\right].\\!\\!\\!\\!\\!$ (43)

Next we will determine the number of paths by the singular value decomposition
(SVD) of ${{\bf{G}}_{\rm{re}}}(t)$:
${{\bf{G}}_{\rm{re}}}(t)={{\bf{U}}_{\tau}}{{\bf{S}}_{\tau}}{\bf{V}}_{\tau}^{H}$,
where ${{\bf{U}}_{\tau}}=[{{\bf{u}}_{1}},\cdots,{{\bf{u}}_{\mu_{1}}}]$,
${{\bf{S}}_{\tau}}={\rm{diag}}\\{{s_{1}},\cdots,{s_{{M_{\tau}}}}\\}$, and
${M_{\tau}}=\min(\mu_{1},2\mu_{2})$. Without loss of generality, we assume the
diagonal elements of ${{\bf{S}}_{\tau}}$ are in non-increasing order. To
determine the number of non-negligible paths, we define a set

$\begin{array}[]{l}{\mathfrak{M}}=\\{{m_{\tau}}|\lvert{{s_{{m_{\tau}}}}}\rvert\geq{\gamma_{1}},1\leq{m_{\tau}}\leq{M_{\tau}}\\}\\\
\ \ \ \
=\\{{\begin{array}[]{*{20}{c}}{m_{1},}&{\cdots,}&{m_{P}}\end{array}}\\},\end{array}\\!\\!\\!\\!\\!$
(44)

where the positive threshold ${\gamma_{1}}$ is close to zero. The number of
non-negligible paths $P={\rm{card}}({\mathfrak{M}})$.

Hereinafter, we estimate the path delay, which contains two substeps, i.e.,
generating a real matrix, and the EVD. We first define a matrix
${{\bf{U}}_{s}}=[{{\bf{u}}_{{m_{1}}}},\cdots,{{\bf{u}}_{{m_{P}}}}]$, which is
generated by selecting the $P$ columns from ${{\bf{U}}_{\tau}}$. A real matrix
related to the path delays is defined as

${{\bf{\Psi}}_{\tau}}={\left[{{\mathop{\rm
Re}\nolimits}({\bf{Q}}_{\mu_{3}}^{H}{{\bf{J}}_{1}}{\bf{Q}}_{\mu_{1}}){{\bf{U}}_{s}}}\right]^{\dagger}}{\mathop{\rm
Im}\nolimits}({\bf{Q}}_{\mu_{3}}^{H}{{\bf{J}}_{1}}{\bf{Q}}_{\mu_{1}}){{\bf{U}}_{s}},\\!$
(45)

where ${\mu_{3}}=KRL-RL$ and
${{\bf{J}}_{1}}=\left[{{{\bf{I}}_{\mu_{3}}}:{{\bf{0}}_{{\mu_{3}}\times
RL}}}\right]$. The unitary matrix ${{\bf{Q}}_{\mu_{3}}}$ only depends on the
size ${\mu_{3}}$, and is expressed as Eq. (42) or Eq. (43). We also introduce
a matrix ${{\bf{\hat{Z}}}_{\tau}}$, which contains the delay information:

${{\bf{\hat{Z}}}_{\tau}}={\rm{diag}}\left\\{\\!\\!{\begin{array}[]{*{20}{c}}{\tan\left({\pi\Delta
f{\tau_{1}}(t)}\right),}\\!\\!&\\!\\!{\cdots,}\\!\\!&\\!\\!{\tan\left({\pi\Delta
f{\tau_{{P}}}(t)}\right)}\end{array}}\\!\\!\right\\}.\\!$ (46)

The two matrices ${{\bf{\Psi}}_{\tau}}$ and ${{\bf{\hat{Z}}}_{\tau}}$ are
similarity matrices that share the same eigenvalues [28]. Then, the path delay
will be estimated by the EVD of ${{\bf{\Psi}}_{\tau}}$:
${{\bf{\Psi}}_{\tau}}={{\bf{W}}_{\tau}}{{\bf{\hat{Z}}}_{\tau}}{{\bf{W}}_{\tau}^{-1}}$,
where ${{\bf{W}}_{\tau}}$ is an eigenvectors collection of
${{\bf{\Psi}}_{\tau}}$. Since ${{\bf{\Psi}}_{\tau}}$ and
${{\bf{\hat{Z}}}_{\tau}}$ are similarity matrices, ${{\bf{\hat{Z}}}_{\tau}}$
can be calculated by
${{\bf{\hat{Z}}}_{\tau}}={\bf{W}}_{\tau}^{-1}{{\bf{\Psi}}_{\tau}}{{\bf{W}}_{\tau}}$.
According to the angular-frequency-domain channel structure in Eq. (13), the
estimated diagonal elements in ${{\bf{\hat{Z}}}_{\tau}}$ may reflect the
multipath delay phase differences between two neighboring subcarriers. Thus,
the $p$-th path delay ${\hat{\tau}_{p}}(t)$ is estimated as

${\hat{\tau}_{p}}(t)=\frac{{{{\tan}^{-1}}\left({{{\hat{z}}_{\tau,p}}}\right)}}{{\pi\Delta
f}},\\!$ (47)

where ${\hat{z}_{\tau,p}}$ denotes the $p$-th estimated diagonal element of
${{\bf{\hat{Z}}}_{\tau}}$. Moreover, the delay changing rate
${\hat{k}_{{\tau_{p}}}}$ is estimated with two different time samples:

${\hat{k}_{{\tau_{p}}}}=\frac{{{{\hat{\tau}}_{p}}({t_{2}})-{{\hat{\tau}}_{p}}({t_{1}})}}{{{t_{2}}-{t_{1}}}}.\\!$
(48)

Up to now, the multipath delays have been estimated. Next, we estimate the
multipath EODs and AODs mainly by two substeps, i.e., introducing some
shuffling matrices, generating two real matrices related to the EODs and the
AODs. Since ${{\bf{G}}_{u}}(t)$ contains the angular information of the paths,
we only need to permute ${{\bf{G}}_{u}}(t)$ to generate two new 3-D MP
matrices ${\bf{G}}_{h}(t)$ and ${\bf{G}}_{v}(t)$. First, four shuffling
matrices ${{\bf{S}}_{{\rm{left}},h}}$, ${{\bf{S}}_{{\rm{right}},h}}$,
${{\bf{S}}_{{\rm{left}},v}}$ and ${{\bf{S}}_{{\rm{right}},v}}$ are introduced
as

$\begin{array}[]{l}\\!\\!\\!\\!\\!\\!\\!\\!{{\bf{S}}_{{\rm{left}},h}}\\!=\\![{\bf{s}}(1),\\!\cdots\\!,{\bf{s}}(1\\!+\\!\mu_{1}\\!-\\!L),{\bf{s}}(2),\cdots,\\\
\ \
{\bf{s}}(2\\!+\\!\mu_{1}\\!-\\!L),\\!\cdots\\!,{\bf{s}}(L),\\!\cdots\\!{\bf{s}}(\mu_{1}){]^{T}},\end{array}\\!\\!\\!\\!$
(49)
$\begin{array}[]{l}\\!\\!\\!\\!\\!{{\bf{S}}_{{\rm{right}},h}}\\!=\\![{\bf{s}}(1),\\!\cdots\\!,{\bf{s}}(\mu_{2}\\!-\\!(N_{h}^{{\rm{BS}}}-L)),{\bf{s}}(2),\cdots,\\\
\ \ \ \ \ \ \
{\bf{s}}(\mu_{2}\\!-\\!(N_{h}^{{\rm{BS}}}\\!-\\!L\\!-\\!1)),\\!\cdots\\!,{\bf{s}}(N_{h}^{{\rm{BS}}}\\!-\\!L\\!+\\!1),\\!\cdots\\!{\bf{s}}(\mu_{2}){]^{T}},\end{array}\\!\\!\\!\\!$
(50)
$\begin{array}[]{l}\\!\\!\\!\\!\\!\\!{{\bf{S}}_{{\rm{left}},v}}\\!=\\!\\![{\bf{s}}(1),\\!\cdots\\!,\\!{\bf{s}}(L),\\!\cdots\\!,{\bf{s}}(1\\!\\!+\\!\\!\mu_{1}\\!\\!-\\!\\!LR),\\!\cdots\\!,{\bf{s}}(L\\!\\!+\\!\\!\mu_{1}\\!\\!-\\!\\!LR),\\!\cdots\\!,\\\
\ \ \ \ \
{\bf{s}}(1\\!\\!+\\!\\!(R\\!\\!-\\!\\!1)L),\\!\cdots\\!,{\bf{s}}(L\\!\\!+\\!\\!(R\\!\\!-\\!\\!1)L),\\!\cdots\\!,\\\
\ \ \ \ \
{\bf{s}}(1\\!\\!+\\!\\!(R\\!\\!-\\!\\!1)L\\!\\!+\\!\\!(K\\!\\!-\\!\\!1)LR),\\!\cdots\\!,{\bf{s}}(\mu_{1}){]^{T}},\end{array}\\!\\!\\!\\!$
(51)
$\begin{array}[]{l}\\!\\!\\!\\!\\!\\!{{\bf{S}}_{{\rm{right}},v}}\\!=\\!\\![{\bf{s}}(1),\\!\cdots\\!,\\!{\bf{s}}(N_{h}^{{\rm{BS}}}\\!\\!-\\!\\!L\\!\\!+\\!\\!1),\\!\cdots\\!,\\\
\ \ \ \ \
{\bf{s}}((N_{v}^{{\rm{BS}}}\\!\\!-\\!\\!R\\!\\!+\\!\\!1)(N_{h}^{{\rm{BS}}}\\!\\!-\\!\\!L\\!\\!+\\!\\!1)\\!\\!-\\!\\!(N_{h}^{{\rm{BS}}}\\!\\!-\\!\\!L)),\\!\cdots\\!,\\\
\ \ \ \ \
{\bf{s}}((N_{v}^{{\rm{BS}}}-R+1)(N_{h}^{{\rm{BS}}}-L+1)),\\!\cdots\\!,\\\ \ \
\ \ \
{\bf{s}}(1\\!\\!-\\!\\!(N_{h}^{{\rm{BS}}}\\!\\!-\\!\\!L\\!\\!+\\!\\!1)\\!\\!+\\!\\!\mu_{2}),\\!\cdots\\!,{\bf{s}}(\mu_{2}){]^{T}},\end{array}\\!\\!\\!\\!$
(52)

where ${\bf{s}}(m),m=1,\cdots,LRK$, is an $LRK\times 1$ unit vector with the
$m$-th element being one. Permuted from ${{\bf{G}}_{u}}(t)$, the two matrices
${\bf{G}}_{h}(t)$ and ${\bf{G}}_{v}(t)$ are expressed as
${\bf{G}}_{h}(t)={{\bf{S}}_{{\rm{left}},h}}{{\bf{G}}_{u}}(t){{{\bf{S}}_{{\rm{right}},h}^{H}}}$
and
${\bf{G}}_{v}(t)={{\bf{S}}_{{\rm{left}},v}}{{\bf{G}}_{u}}(t){{{\bf{S}}_{{\rm{right}},v}^{H}}}$,
respectively. Notice that in practice, the two 3-D MP matrices
${\bf{G}}_{h}(t)$ and ${\bf{G}}_{v}(t)$, as well as the two shuffling matrices
${{\bf{S}}_{{\rm{right}},h}}$ and ${{\bf{S}}_{{\rm{right}},v}}$, are not
needed. In fact, the EODs and the AODs might be estimated by only introducing
${{\bf{S}}_{{\rm{left}},h}}$ and ${{\bf{S}}_{{\rm{left}},v}}$.

Then, we will focus on the generation of the two real matrices related to the
EODs and the AODs. By following a similar procedure to generate
${{\bf{\Psi}}_{\tau}}$ from ${{\bf{G}}_{u}}(t)$, i.e., from Eq. (40) to Eq.
(45), we may obtain two real matrices ${{\bf{\Psi}}_{\phi}}$ and
${{\bf{\Psi}}_{\theta}}$, which are generated from ${\bf{G}}_{h}(t)$ and
${\bf{G}}_{v}(t)$ respectively. The matrices ${{\bf{\Psi}}_{\phi}}$ and
${{\bf{\Psi}}_{\theta}}$ are defined by

$\begin{array}[]{l}{{\bf{\Psi}}_{\phi}}={\left({{\mathop{\rm
Re}\nolimits}({\bf{Q}}_{\mu_{4}}^{H}{{\bf{J}}_{2}}{\bf{Q}}_{\mu_{1}}){\bf{Q}}_{\mu_{1}}{{\bf{S}}_{{\rm{left}},h}}{\bf{Q}}_{\mu_{1}}^{H}{{\bf{U}}_{s}}}\right)^{\dagger}}\\\
\ \ \ \ \ \ \ \ \ {\mathop{\rm
Im}\nolimits}({\bf{Q}}_{\mu_{4}}^{H}{{\bf{J}}_{2}}{\bf{Q}}_{\mu_{1}}){\bf{Q}}_{\mu_{1}}{{\bf{S}}_{{\rm{left}},h}}{\bf{Q}}_{\mu_{1}}^{H}{{\bf{U}}_{s}},\end{array}\\!\\!\\!\\!$
(53) $\begin{array}[]{l}{{\bf{\Psi}}_{\theta}}={\left({{\mathop{\rm
Re}\nolimits}({\bf{Q}}_{\mu_{5}}^{H}{{\bf{J}}_{3}}{\bf{Q}}_{\mu_{1}}){\bf{Q}}_{\mu_{1}}{{\bf{S}}_{{\rm{left}},v}}{\bf{Q}}_{\mu_{1}}^{H}{{\bf{U}}_{s}}}\right)^{\dagger}}\\\
\ \ \ \ \ \ \ {\mathop{\rm
Im}\nolimits}({\bf{Q}}_{\mu_{5}}^{H}{{\bf{J}}_{3}}{\bf{Q}}_{\mu_{1}}){\bf{Q}}_{\mu_{1}}{{\bf{S}}_{{\rm{left}},v}}{\bf{Q}}_{\mu_{1}}^{H}{{\bf{U}}_{s}},\end{array}\\!\\!\\!\\!$
(54)

where ${\mu_{4}}=KRL-KR$, ${\mu_{5}}=KRL-KL$,
${{\bf{J}}_{2}}=\left[{{{\bf{I}}_{\mu_{4}}}:{{\bf{0}}_{{\mu_{4}}\times
KR}}}\right]$, and
${{\bf{J}}_{3}}=\left[{{{\bf{I}}_{\mu_{5}}}:{{\bf{0}}_{{\mu_{5}}\times
KL}}}\right]$. The two unitary matrices ${{\bf{Q}}_{\mu_{4}}}$ and
${{\bf{Q}}_{\mu_{5}}}$ are defined in Eq. (42) or Eq. (43) according to their
sizes ${\mu_{4}}$ and ${\mu_{5}}$.

Define two real diagonal matrices ${{\bf{\hat{Z}}}_{\phi}}$ and
${{\bf{\hat{Z}}}_{\theta}}$, which contain the angular information

$\begin{array}[]{l}{{{\bf{\hat{Z}}}}_{\phi}}={\rm{diag}}\\{\tan(\frac{{\pi
d_{h}^{{\rm{tx}}}\cos({\theta_{1}})\sin({\phi_{1}})}}{{{\lambda_{0}}}}),\\\ \
\ \ \ \ \ \ \ \cdots,\tan(\frac{{\pi
d_{h}^{{\rm{tx}}}\cos({\theta_{{P}}})\sin({\phi_{{P}}})}}{{{\lambda_{0}}}})\\},\end{array}\\!\\!\\!\\!$
(55)
$\begin{array}[]{l}\\!\\!\\!\\!\\!\\!\\!\\!{{\bf{\hat{Z}}}_{\theta}}={\rm{diag}}\\{\\!{\begin{array}[]{*{20}{c}}{\tan(\frac{{\pi
d_{v}^{{\rm{tx}}}\sin({\theta_{1}})}}{{{\lambda_{0}}}}),}\\!&\\!{\cdots,}\\!&\\!{\tan(\frac{{\pi
d_{v}^{{\rm{tx}}}\sin({\theta_{{P}}})}}{{{\lambda_{0}}}})}\end{array}}\\!\\}.\end{array}\\!\\!\\!\\!$
(56)

The diagonal elements in ${{\bf{\hat{Z}}}_{\phi}}$ reflect the phase
differences between the two neighboring antenna elements in the horizontal
direction of the BS antenna panel of the paths. Likewise,
${{\bf{\hat{Z}}}_{\theta}}$ reflects the phase differences between the two
neighboring antenna elements, in the vertical direction of the BS antenna
panel for the paths.

Finally, we will estimate the EODs and the AODs as follows: The two real
matrices ${{\bf{\Psi}}_{\phi}}$ and ${{\bf{\Psi}}_{\theta}}$ are similar to
${{\bf{\hat{Z}}}_{\phi}}$ and ${{\bf{\hat{Z}}}_{\theta}}$, respectively, i.e.,
${{\bf{\Psi}}_{\phi}}$ and ${{\bf{\hat{Z}}}_{\phi}}$ share the same
eigenvalues, ${{\bf{\Psi}}_{\theta}}$ and ${{\bf{\hat{Z}}}_{\theta}}$ also
share the same eigenvalues. Thus, we may obtain ${{\bf{\hat{Z}}}_{\phi}}$ and
${{\bf{\hat{Z}}}_{\theta}}$ via the EVDs of ${{\bf{\Psi}}_{\phi}}$ and
${{\bf{\Psi}}_{\theta}}$. Note that ${{\bf{\Psi}}_{\phi}}$,
${{\bf{\Psi}}_{\theta}}$, and ${{\bf{\Psi}}_{\tau}}$ share the same
eigenvectors, since ${\bf{G}}_{h}(t)$ and ${\bf{G}}_{v}(t)$ are simply
permuted versions of ${{\bf{G}}_{u}}(t)$. More specifically, we may denote the
EVDs by
${{\bf{\Psi}}_{\phi}}={{\bf{W}}_{\tau}}{{\bf{\hat{Z}}}_{\phi}}{{\bf{W}}_{\tau}^{-1}}$
and
${{\bf{\Psi}}_{\theta}}={{\bf{W}}_{\tau}}{{\bf{\hat{Z}}}_{\theta}}{{\bf{W}}_{\tau}^{-1}}$,
where ${\bf{W}}_{\tau}$ contains the eigenvectors of ${{\bf{\Psi}}_{\tau}}$.
As a result, ${{\bf{\hat{Z}}}_{\phi}}$ and ${{\bf{\hat{Z}}}_{\theta}}$ are
calculated as
${{\bf{\hat{Z}}}_{\phi}}={\bf{W}}_{\tau}^{-1}{{\bf{\Psi}}_{\phi}}{{\bf{W}}_{\tau}}$
and
${{\bf{\hat{Z}}}_{\theta}}={\bf{W}}_{\tau}^{-1}{{\bf{\Psi}}_{\theta}}{{\bf{W}}_{\tau}}$.
The EOD and the AOD are estimated by

${\hat{\theta}_{p}}={\sin^{-1}}\left({\frac{{{\tan^{-1}}({{\hat{z}}_{\theta,p}}){\lambda_{0}}}}{{\pi
d_{v}^{{\rm{tx}}}}}}\right),\\!\\!\\!\\!$ (57)
${\hat{\phi}_{p}}={\sin^{-1}}\left({\frac{{{\tan^{-1}}({{\hat{z}}_{\phi,p}}){\lambda_{0}}}}{{\pi
d_{h}^{{\rm{tx}}}\cos\left({{{\sin}^{-1}}\left({\frac{{{\tan^{-1}}({{\hat{z}}_{\theta,p}}){\lambda_{0}}}}{{\pi
d_{v}^{{\rm{tx}}}}}}\right)}\right)}}}\right),\\!\\!\\!\\!$ (58)

where ${\hat{z}_{\phi,p}}$ and ${\hat{z}_{\theta,p}}$ denote the $p$-th
estimated diagonal elements of ${{\bf{\hat{Z}}}_{\phi}}$ and
${{\bf{\hat{Z}}}_{\theta}}$, respectively. So far, the EODs, AODs and delays
have been estimated by our proposed 3-D angle-delay estimation scheme, which
is summarized in Algorithm 1. Notice that ${{\bf{H}}_{u}}(t)$ as an input of
Algorithm 1 is a historical sample obtained from the UL sounding reference
signal (SRS). With Algorithm 1, we will estimate the channel parameters, which
are used for reconstructing and predicting the future DL channel.

We now analyze the computational complexity of this algorithm in terms of time
and memory. The time complexity is dominated by step 3 - step 8. The
complexity of step 3 is
$\mathcal{O}(2\mu_{1}^{2}{\mu_{2}})+\mathcal{O}(4{\mu_{1}}\mu_{2}^{2})$. The
main complexity of step 4 is the SVD, i.e.,
$\mathcal{O}(4{\mu_{1}}\mu_{2}^{2})$. Step 5 contains the matrix inversion and
multiplications, and has a complexity order of
$\mathcal{O}(\mu_{3}^{2}{\mu_{1}})+\mathcal{O}({\mu_{3}}\mu_{1}^{2})+\mathcal{O}({P^{2.37}})$.
In step 6, the computation complexity is dominated by the EVD, i.e.,
$\mathcal{O}({P^{3}})$. The complexity of step 7 is
$\mathcal{O}(\mu_{4}^{2}{\mu_{1}})+\mathcal{O}({\mu_{4}}\mu_{1}^{2})+\mathcal{O}(\mu_{1}^{3})+\mathcal{O}({P^{2.37}})$,
which mainly contains the matrix inversion and multiplication. Similar to step
6, step 8 has a complexity order of $\mathcal{O}({P^{3}})$. According to the
conditions of pencil sizes in Eq. (39), we may obtain the global time
complexity of Algorithm 1 as
$\mathcal{O}(2{\mu_{1}^{2}}{\mu_{2}})+\mathcal{O}(4{\mu_{1}}{\mu_{2}^{2}})+\mathcal{O}({\mu_{1}^{3}})$.
Then, the memory complexity of Algorithm 1 is dominated by step 3, i.e., the
two matrices ${{\bf{Q}}_{\mu_{1}}}$ and ${{\bf{Q}}_{\mu_{2}}}$ in step 3 lead
to the total memory complexity order of
$\mathcal{S}(\mu_{1}^{2})+\mathcal{S}(4\mu_{2}^{2})$.

Algorithm 1 Proposed 3-D angle-delay estimation scheme.

0: ${{\bf{H}}_{u}}(t)$, $L$, $R$, $K$, ${\gamma_{1}}$;

1: for $t=\left\\{{{t_{1}},{t_{2}}}\right\\}$ do

2: Generate a 3-D MP matrix ${{\bf{G}}_{u}}(t)$ following the procedure
between Eq. (24) and Eq. (36);

3: Convert ${{\bf{G}}_{u}}(t)$ to 3-D real matrix ${{\bf{G}}_{\rm{re}}}(t)$
according to Eq. (40) and Eq. (41);

4: Determine the number of paths by Eq. (44);

5: Generate the real matrix ${{\bf{\Psi}}_{\tau}}$ with Eq. (45);

6: Estimate the delay ${\hat{\tau}_{p}}(t)$ as Eq. (47) by the EVD of
${{\bf{\Psi}}_{\tau}}$;

7: Define two shuffling matrices as Eq. (49) and Eq. (51), and generate
${{\bf{\Psi}}_{\phi}}$ and ${{\bf{\Psi}}_{\theta}}$ with Eq. (53) and Eq.
(54);

8: Estimate the EOD ${\hat{\theta}_{p}}$ and the AOD ${\hat{\phi}_{p}}$ by Eq.
(57) and Eq. (58), respectively;

9: end for

10: Estimate the parameter ${\hat{k}_{{\tau_{p}}}}$ by Eq. (48);

10: ${\hat{\tau}_{p}}(t)$, ${\hat{\theta}_{p}}$, ${\hat{\phi}_{p}}$,
${\hat{k}_{{\tau_{p}}}}$, $p\in\left\\{{1,\cdots,{P}}\right\\}$

###  III-B The 3-D Angle-Doppler Estimation

With the angular and delay information obtained in Sec. III-A, we now estimate
the multipath angles and Doppler in this section. Since the duration of the
channel sample is $T$, the ${n_{s}}$-th sample time is denoted by
$t={n_{s}}{T}$. Similar to Sec. III-A, we start by defining a 3-D MP matrix
${{\bm{\mathcal{G}}}_{u}}({f_{1}})\in{{\mathbb{C}}^{LRQ\times(N_{h}^{{\rm{BS}}}-L+1)(N_{v}^{{\rm{BS}}}-R+1)({N_{s}}-Q+1)}}$
as

${{\bm{\mathcal{G}}}_{u}}({f_{1}})=\left[\begin{array}[]{l}{{\bm{\mathcal{G}}}_{u}}(1,{f_{1}}),\
\cdots,{{\bm{\mathcal{G}}}_{u}}({N_{s}}\\!-\\!Q\\!+\\!1,{f_{1}})\\\ \ \ \ \
\vdots\ \ \ \ \ \ \ \ \ddots\ \ \ \ \ \ \ \ \ \vdots\\\
{{\bm{\mathcal{G}}}_{u}}(Q,{f_{1}}),\cdots,\ \ \ \
{{\bm{\mathcal{G}}}_{u}}({N_{s}},{f_{1}})\end{array}\right],\\!\\!\\!\\!$ (59)

where ${N_{s}}$ is the number of samples, and $Q$ is the pencil size
satisfying $P<Q<{N_{s}}-P+1$. The matrix ${{\bm{\mathcal{G}}}_{u}}({f_{1}})$
is composed of the channels between the $u$-th UE antenna and all BS antennas
at the first subcarrier and all sample times. For simplicity, we define
${{\omega}_{\mu_{1}}}=LRQ$ and
${{\omega}_{\mu_{2}}}=(N_{h}^{{\rm{BS}}}-L+1)(N_{v}^{{\rm{BS}}}-R+1)({N_{s}}-Q+1)$.
The 2-D MP matrix
${{\bm{\mathcal{G}}}_{u}}({n_{s}},{f_{1}})\in{{\mathbb{C}}^{LR\times(N_{h}^{{\rm{BS}}}-L+1)(N_{v}^{{\rm{BS}}}-R+1)}}$,
${n_{s}}=1,\cdots,{N_{s}}$ consists of the channels between the $u$-th UE
antenna and all BS antennas at the first subcarrier and the ${n_{s}}$-th
sample. We rewrite ${{\bm{\mathcal{G}}}_{u}}({n_{s}},{f_{1}})$ as

${{\bm{\mathcal{G}}}_{u}}({n_{s}},{f_{1}})={{\check{\bf{E}}}_{1}}{{\bf{Y}}_{w}}{\bf{Z}}_{w_{\tau}}^{{n_{s}}}{{\check{\bf{F}}}_{1}},\\!\\!\\!\\!$
(60)

where ${\bf{Z}}_{\omega_{\tau}}$ denotes the phase differences between two
neighboring samples of the paths:

${\bf{Z}}_{\omega_{\tau}}={\rm{diag}}\left\\{{\begin{array}[]{*{20}{c}}{e^{j2\pi{\omega_{\tau_{1}}}{T}},}&{\cdots,}&{e^{j2\pi{\omega_{\tau_{P}}}{T}}}\end{array}}\right\\},\\!\\!\\!\\!$
(61)

where

${\omega_{\tau_{p}}}={\omega_{p}}-{f_{1}}{k_{{\tau_{p}}}}.\\!\\!\\!\\!$ (62)

The matrix ${{\bf{Y}}_{\omega}}$ is defined as

${{\bf{Y}}_{\omega}}={\rm{diag}}\left\\{{{{\bf{D}}_{u}}{\bf{e}}({f_{1}})}\right\\},\\!\\!\\!\\!$
(63)

where ${{\bf{D}}_{u}}$ and ${\bf{e}}({f_{1}})$ are shown in Eq. (19) and Eq.
(21). According to the angular-time-domain channel structure in Eq. (18),
${{\bm{\mathcal{G}}}_{u}}({f_{1}})$ is rewritten as

${{\bm{\mathcal{G}}}_{u}}({f_{1}})={{\check{\bf{E}}}_{\omega,2}}{{\bf{Y}}_{\omega}}{{\check{\bf{F}}}_{\omega,2}},\\!\\!\\!\\!$
(64)

where
${{\check{\bf{E}}}_{\omega,2}}\in{{\mathbb{C}}^{{{\omega}_{\mu_{1}}}\times
P}}$ and
${{\check{\bf{F}}}_{\omega,2}}\in{{\mathbb{C}}^{P\times{{\omega}_{\mu_{2}}}}}$
are defined as

${{\check{\bf{E}}}_{\omega,2}}={\left[\\!{\begin{array}[]{*{20}{c}}{({{\check{\bf{E}}}_{1}}{\bf{Z}}_{\omega_{\tau}}^{0})^{T},}\\!&\\!{\cdots,}\\!&\\!{({{\check{\bf{E}}}_{1}}{\bf{Z}}_{\omega_{\tau}}^{Q-1})^{T}}\end{array}}\\!\right]^{T}},\\!\\!\\!\\!$
(65)
${{\check{\bf{F}}}_{\omega,2}}=\left[{\begin{array}[]{*{20}{c}}{{\bf{Z}}_{\omega_{\tau}}^{0}{{\check{\bf{F}}}_{1}},}&{\cdots,}&{{\bf{Z}}_{\omega_{\tau}}^{{N_{s}}-Q}{{\check{\bf{F}}}_{1}}}\end{array}}\right].\\!\\!\\!\\!$
(66)

The physical meaning of
${{\check{\bf{E}}}_{1}}{\bf{Z}}_{\omega_{\tau}}^{q_{1}}$,
${q_{1}}=0,\cdots,Q-1$ is the 3-D steering vectors of a subset of antenna
array with size $R\times L$ on the BS antenna panel at the $({q_{1}}+1)$-th
sample. Likewise, ${\bf{Z}}_{\omega_{\tau}}^{q_{2}}{{\check{\bf{F}}}_{1}}$,
${q_{2}}=0,\cdots,{N_{s}}-Q$ is the 3-D steering vectors of a subset of
antenna array with size $(N_{v}^{{\rm{BS}}}-R+1)\times(N_{h}^{{\rm{BS}}}-L+1)$
on the BS antenna panel at the $({q_{2}}+1)$-th sample.

The three pencil sizes in ${{\bm{\mathcal{G}}}_{u}}({f_{1}})$ yield

$\left\\{\begin{array}[]{l}\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ LR(Q-1)\geq
P\\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ LQ(R-1)\geq P\\\ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ RQ(L-1)\geq P\\\
(N_{h}^{{\rm{BS}}}-L+1)(N_{v}^{{\rm{BS}}}-R+1)({N_{s}}-Q+1)\geq
P\end{array}\right.\\!\\!\\!\\!\\!\\!.\\!\\!\\!\\!$ (67)

By following the similar estimation procedures of Sec. III-A, we then
transform ${{\bm{\mathcal{G}}}_{u}}({f_{1}})$ to a real matrix as
${{\bm{\mathcal{G}}}_{\rm{re}}}({f_{1}})$. The step of determining the number
of paths is not needed, because this work has been finished in Sec. III-A.
Next, we will estimate the parameter ${\hat{\omega}_{\tau_{p}}}$ by generating
a real matrix and the EVD. More specifically, a real matrix related to Doppler
is given by

$\ \\!\\!{{\bf{\Psi}}_{\omega_{\tau}}}\\!\\!=\\!\\!{\left(\\!{{\mathop{\rm
Re}\nolimits}({\bf{Q}}_{{\omega}_{\mu_{3}}}^{H}{{\bf{J}}_{\omega,1}}{{\bf{Q}}_{{\omega}_{\mu_{1}}}}){{\bf{U}}_{\omega,s}}}\\!\right)^{\dagger}}\\!{\rm
Im}({\bf{Q}}_{{\omega}_{\mu_{3}}}^{H}{{\bf{J}}_{\omega,1}}{{\bf{Q}}_{{\omega}_{\mu_{1}}}}){{\bf{U}}_{\omega,s}},\\!\\!\\!\\!$
(68)

where ${{\omega}_{\mu_{3}}}=QRL-RL$, and
${{\bf{J}}_{\omega,1}}=\left[{{{\bf{I}}_{{\omega}_{\mu_{3}}}}:{{\bf{0}}_{{{\omega}_{\mu_{3}}}\times
RL}}}\right]$. The two unitary matrices ${{\bf{Q}}_{{\omega}_{\mu_{1}}}}$ and
${{\bf{Q}}_{{\omega}_{\mu_{3}}}}$ depending on the size of
${{\omega}_{\mu_{1}}}$ and ${{\omega}_{\mu_{3}}}$, are expressed as Eq. (42)
or Eq. (43). Perform the SVD of ${{\bm{\mathcal{G}}}_{\rm{re}}}({f_{1}})$:
${{\bm{\mathcal{G}}}_{\rm{re}}}({f_{1}})={{\bf{U}}_{\omega}}{{\bf{S}}_{\omega}}{\bf{V}}_{\omega}^{H}$.
We select the ${P}$ columns from ${{\bf{U}}_{\omega}}$ as
${{\bf{U}}_{\omega,s}}$, which follows the similar procedure of
${{\bf{U}}_{s}}$ in Sec. III-A. The matrix ${{\bf{\Psi}}_{\omega_{\tau}}}$ and
a real diagonal matrix ${{\bf{\hat{Z}}}_{\omega_{\tau}}}$ are similarity
matrices, and share the same eigenvalues. The matrix
${{\bf{\hat{Z}}}_{\omega_{\tau}}}$ containing the Doppler information is
defined as

${{{\bf{\hat{Z}}}}_{{\omega_{\tau}}}}={\rm{diag}}\\{\\!{\begin{array}[]{*{20}{c}}{-\tan({\omega_{{\tau_{1}}}}\pi
T),}\\!&\\!{\cdots,}\\!&\\!{-\tan({\omega_{{\tau_{{P}}}}}\pi
T)}\end{array}}\\!\\}.\\!\\!\\!\\!$ (69)

The diagonal elements in ${{{\bf{\hat{Z}}}}_{{\omega_{\tau}}}}$ reflect the
phase differences between two samples of the paths.

Until now, the real matrix ${{\bf{\Psi}}_{\omega_{\tau}}}$ has been generated,
and the following step is the EVD. We perform the EVD of
${{\bf{\Psi}}_{\omega_{\tau}}}$:
${{\bf{\Psi}}_{\omega_{\tau}}}={{\bf{W}}_{\omega_{\tau}}}{{\bf{\hat{Z}}}_{\omega_{\tau}}}{{\bf{W}}_{\omega_{\tau}}^{-1}}$,
where ${{\bf{W}}_{\omega_{\tau}}}$ is composed of the eigenvectors of
${{\bf{\Psi}}_{\omega_{\tau}}}$. The matrix ${{\bf{\hat{Z}}}_{\omega_{\tau}}}$
is calculated as
${{\bf{\hat{Z}}}_{\omega_{\tau}}}={\bf{W}}_{\omega_{\tau}}^{-1}{{\bf{\Psi}}_{\omega_{\tau}}}{{\bf{W}}_{\omega_{\tau}}}$.
The parameter ${\hat{\omega}_{\tau_{p}}}$ is thus estimated by

${\hat{\omega}_{\tau_{p}}}=-\frac{{{{\tan}^{-1}}({{\hat{z}}_{\omega_{\tau},p}})}}{{\pi{T}}},\\!\\!\\!\\!$
(70)

where ${\hat{z}_{\omega_{\tau},p}}$ denotes the $p$-th estimated diagonal
element of ${{\bf{\hat{Z}}}_{\omega_{\tau}}}$. Likewise, we will estimate the
EODs and the AODs as follows: The two real matrices related to the EODs and
the AODs are given by

$\begin{array}[]{l}\\!\\!\\!\\!\\!\\!\\!\\!{{\bf{\Psi}}_{\omega_{\tau},\phi}}\\!=\\!\\!{\left(\\!{{\mathop{\rm
Re}\nolimits}({\bf{Q}}_{{\omega}_{\mu_{4}}}^{H}\\!{{\bf{J}}_{\omega,2}}{{\bf{Q}}_{{\omega}_{\mu_{1}}}}\\!){\bf{Q}}_{{\omega}_{\mu_{1}}}\\!{{\bf{S}}_{{\rm{left}},\omega,h}}{\bf{Q}}_{{\omega}_{\mu_{1}}}^{H}\\!\\!{{\bf{U}}_{\omega,s}}}\\!\right)^{\dagger}}\\\
\ \ \ \ \ \ \ {\mathop{\rm
Im}\nolimits}({\bf{Q}}_{{\omega}_{\mu_{4}}}^{H}\\!{{\bf{J}}_{\omega,2}}{\bf{Q}}_{{\omega}_{\mu_{1}}}\\!){\bf{Q}}_{{\omega}_{\mu_{1}}}\\!\\!{{\bf{S}}_{{\rm{left}},\omega,h}}{\bf{Q}}_{{\omega}_{\mu_{1}}}^{H}{{\bf{U}}_{\omega,s}},\end{array}\\!\\!$
(71)
$\begin{array}[]{l}\\!\\!\\!\\!\\!\\!\\!{{\bf{\Psi}}_{\omega_{\tau},\theta}}\\!=\\!\\!{\left(\\!{{\mathop{\rm
Re}\nolimits}({\bf{Q}}_{{\omega}_{\mu_{5}}}^{H}{{\bf{J}}_{\omega,3}}{{\bf{Q}}_{{\omega}_{\mu_{1}}}}\\!){\bf{Q}}_{{\omega}_{\mu_{1}}}{{\bf{S}}_{{\rm{left}},\omega,v}}{\bf{Q}}_{{\omega}_{\mu_{1}}}^{H}{{\bf{U}}_{\omega,s}}}\\!\right)^{\dagger}}\\\
\ \ \ \ \ \ {\mathop{\rm
Im}\nolimits}({\bf{Q}}_{{\omega}_{\mu_{5}}}^{H}{{\bf{J}}_{\omega,3}}{\bf{Q}}_{{\omega}_{\mu_{1}}}\\!){\bf{Q}}_{{\omega}_{\mu_{1}}}{{\bf{S}}_{{\rm{left}},\omega,v}}{\bf{Q}}_{{\omega}_{\mu_{1}}}^{H}{{\bf{U}}_{\omega,s}},\end{array}\\!\\!$
(72)

where ${{\omega}_{\mu_{4}}}=QRL-QR$, ${{\omega}_{\mu_{5}}}=QRL-QL$,
${{\bf{J}}_{\omega,2}}=\left[{{{\bf{I}}_{{\omega}_{\mu_{4}}}}:{{\bf{0}}_{{{\omega}_{\mu_{4}}}\times
QR}}}\right]$, and
${{\bf{J}}_{\omega,3}}=\left[{{{\bf{I}}_{{\omega}_{\mu_{5}}}}:{{\bf{0}}_{{{\omega}_{\mu_{5}}}\times
QL}}}\right]$. Define ${{\bf{S}}_{{\rm{left}},\omega,h}}$ and
${{\bf{S}}_{{\rm{left}},\omega,v}}$ as two new shuffling matrices that share
the similar expressions with ${{\bf{S}}_{{\rm{left}},h}}$ and
${{\bf{S}}_{{\rm{left}},v}}$ in Eq. (49) and Eq. (51). The two shuffling
matrices ${{\bf{S}}_{{\rm{left}},\omega,h}}$ and
${{\bf{S}}_{{\rm{left}},\omega,v}}$ are generated by replacing the parameter
$K$ with $Q$ in Eq. (49) and Eq. (51). The two matrices
${{\bf{\Psi}}_{\omega_{\tau},\phi}}$ and
${{\bf{\Psi}}_{\omega_{\tau},\theta}}$ are similar to two real diagonal
matrices ${{\bf{\hat{Z}}}_{\omega_{\tau},\phi}}$ and
${{\bf{\hat{Z}}}_{\omega_{\tau},\theta}}$:

$\begin{array}[]{l}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!{{{\bf{\hat{Z}}}}_{{\omega_{\tau}},\theta}}\\!=\\!{\rm{diag}}\\{\tan(\frac{{\pi
d_{v}^{{\rm{tx}}}\sin({\theta_{\omega,1}})}}{{{\lambda_{0}}}}),\cdots\\!,\tan(\frac{{\pi
d_{v}^{{\rm{tx}}}\sin({\theta_{\omega,{{P}}}})}}{{{\lambda_{0}}}})\\},\end{array}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
(73)
$\begin{array}[]{l}{{{\bf{\hat{Z}}}}_{{\omega_{\tau}},\phi}}={\rm{diag}}\\{\tan(\frac{{\pi
d_{h}^{{\rm{tx}}}\cos({\theta_{\omega,1}}){\rm{sin}}({\phi_{\omega,1}})}}{{{\lambda_{0}}}}),\\\
\ \ \ \ \ \ \ \ \ \ \cdots,\tan(\frac{{\pi
d_{h}^{{\rm{tx}}}\cos({\theta_{\omega,{{P}}}}){\rm{sin}}({\phi_{\omega,{{P}}}})}}{{{\lambda_{0}}}})\\},\end{array}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
(74)

which are given by
${{\bf{\hat{Z}}}_{\omega_{\tau},\phi}}={\bf{W}}_{\omega_{\tau}}^{-1}{{\bf{\Psi}}_{\omega_{\tau},\phi}}{{\bf{W}}_{\omega_{\tau}}}$
and
${{\bf{\hat{Z}}}_{\omega_{\tau},\theta}}={\bf{W}}_{\omega_{\tau}}^{-1}{{\bf{\Psi}}_{\omega_{\tau},\theta}}{{\bf{W}}_{\omega_{\tau}}}$.
The EOD and the AOD are estimated by

${\hat{\theta}_{\omega,p}}={\sin^{-1}}\left({\frac{{{\tan^{-1}}({{\hat{z}}_{\omega_{\tau},\phi,p}}){\lambda_{0}}}}{{\pi
d_{v}^{{\rm{tx}}}}}}\right),\\!\\!\\!\\!$ (75) $\
{\hat{\phi}_{\omega,p}}\\!=\\!{\sin^{-1}}\\!\left({\frac{{{\tan^{-1}}({{\hat{z}}_{\omega_{\tau},\theta,p}}){\lambda_{0}}}}{{\pi
d_{h}^{{\rm{tx}}}\cos\left({\sin{{}^{-1}}\left({\frac{{{\tan^{-1}}({{\hat{z}}_{\omega_{\tau},\phi,p}}){\lambda_{0}}}}{{\pi
d_{v}^{{\rm{tx}}}}}}\right)}\right)}}}\right),\\!\\!\\!\\!$ (76)

where ${\hat{z}_{\omega_{\tau},\phi,p}}$ and
${\hat{z}_{\omega_{\tau},\theta,p}}$ denote the $p$-th estimated diagonal
elements of ${{\bf{\hat{Z}}}_{\omega_{\tau},\phi}}$ and
${{\bf{\hat{Z}}}_{\omega_{\tau},\theta}}$.

The procedures of our proposed 3-D angle-Doppler estimation scheme are
summarized in Algorithm 2. Similar to the time and memory complexity analysis
of Algorithm 1, we may obtain the time complexity of Algorithm 2 as
$\mathcal{O}(2{\omega_{{\mu_{1}}}^{2}}{\omega_{{\mu_{2}}}})+\mathcal{O}(4{\omega_{{\mu_{1}}}}{\omega_{{\mu_{2}}}^{2}})+\mathcal{O}({\omega_{{\mu_{1}}}^{3}})$,
and the memory complexity as
$\mathcal{S}({\omega_{{\mu_{1}}}^{2}})+\mathcal{S}(4{\omega_{{\mu_{2}}}^{2}})$.

Algorithm 2 Proposed 3-D angle-Doppler estimation scheme.

0: ${{\bf{H}}_{u}}(t)$, $L$, $R$, $Q$, ${N_{s}}$, ${P}$;

1: Generate ${{\bm{\mathcal{G}}}_{u}}({f_{1}})$ according to the procedure
between Eq. (59) and Eq. (64);

2: Convert ${{\bm{\mathcal{G}}}_{u}}({f_{1}})$ to a real matrix
${{\bm{\mathcal{G}}}_{\rm{re}}}({f_{1}})$;

3: Generate the matrix ${{\bf{\Psi}}_{\omega_{\tau}}}$ with Eq. (68);

4: Estimate ${\hat{\omega}_{\tau_{p}}}$ by Eq. (70) with the EVD of
${{\bf{\Psi}}_{\omega_{\tau}}}$;

5: Generate two matrices ${{\bf{\Psi}}_{\omega_{\tau},\phi}}$ and
${{\bf{\Psi}}_{\omega_{\tau},\theta}}$ with Eq. (71) and Eq. (72);

6: Estimate the EOD ${\hat{\theta}_{\omega,p}}$ and the AOD
${\hat{\phi}_{\omega,p}}$ by Eq. (75) and Eq. (76);

6: ${\hat{\omega}_{\tau_{p}}}$, ${\hat{\theta}_{\omega,p}}$,
${\hat{\phi}_{\omega,p}}$, $p\in\left\\{{1,\cdots,{P}}\right\\}$

Although we have obtained the angular-delay information of the paths in
Algorithm 1, and the angular-Doppler information of the paths in Algorithm 2,
the angles, delay, and Doppler of the paths have to be paired so as to
reconstruct the channel.

###  III-C The Proposed Angle-Delay-Doppler Pairing Scheme

We propose to pair the angles, delay, and Doppler of the paths, with the
super-resolution property of the angle estimation. Let some ${P}\times 1$
vectors $\bm{\hat{\theta}}$, $\bm{\hat{\phi}}$, $\bm{\hat{\tau}(t)}$,
$\bm{\hat{k}_{\tau}}$, $\bm{\hat{\theta}_{\omega}}$,
$\bm{\hat{\phi}_{\omega}}$ and $\bm{\hat{\omega}_{\tau}}$ denote the
corresponding estimated parameters of multipaths. A pairing matrix
${[{{\bf{s}}_{{\rm{pair}},1}},\cdots,{{\bf{s}}_{{\rm{pair}},{P}}}]}$ is
defined to map the relationship between $[\bm{\hat{\theta}},\bm{\hat{\phi}}]$
and $[\bm{\hat{\theta}_{\omega}},\bm{\hat{\phi}_{\omega}}]$ as

$[\bm{\hat{\theta}},\bm{\hat{\phi}}]={[{{\bf{s}}_{{\rm{pair}},1}},\cdots,{{\bf{s}}_{{\rm{pair}},{P}}}]}[\bm{\hat{\theta}_{\omega}},\bm{\hat{\phi}_{\omega}}],\\!\\!\\!\\!$
(77)

where ${{\bf{s}}_{{\rm{pair}},p}}=[{s_{p,1}},\cdots,{s_{p,{P}}}]^{T}$,
$p=1,\cdots,{P}$ is a ${P}\times 1$ unit vector with only one element being
one. Denote the $p$-th row of
$[{\bm{\hat{\theta}_{w}}},{\bm{\hat{\phi}_{w}}}]$ as
$[{\hat{\theta}_{w,p}},{\hat{\phi}_{w,p}}]$. By solving the minimization
problem

${{\bf{s}}_{{\rm{pair}},p}}=\arg\mathop{\min}\limits_{{{\bf{s}}_{{\rm{pair}},p}}}(\lvert{{{\hat{\theta}}_{w,p}}{{\bf{1}}_{{P}\times{1}}}-\bm{\hat{\theta}}}\rvert+\lvert{{{\hat{\phi}}_{w,p}}{{\bf{1}}_{{P}\times{1}}}-\bm{\hat{\phi}}}\rvert),\\!\\!\\!\\!$
(78)

we may obtain ${{\bf{s}}_{\text{pair},p}}$. Specifically,
$\lvert{{{\hat{\theta}}_{w,p}}{{\bf{1}}_{{P}\times{1}}}-\bm{\hat{\theta}}}\rvert$
is the error vector between ${\hat{\theta}_{w,p}}$ and $\bm{\hat{\theta}}$,
and
$\lvert{{{\hat{\phi}}_{w,p}}{{\bf{1}}_{{P}\times{1}}}-\bm{\hat{\phi}}}\rvert$
is the error vector between ${\hat{\phi}_{w,p}}$ and $\bm{\hat{\phi}}$. The
row index of $[{\bm{\hat{\theta}_{w}}},{\bm{\hat{\phi}_{w}}}]$ satisfying the
minimize value of
$\lvert{{{\hat{\theta}}_{w,p}}{{\bf{1}}_{{P}\times{1}}}-\bm{\hat{\theta}}}\rvert+\lvert{{{\hat{\phi}}_{w,p}}{{\bf{1}}_{{P}\times{1}}}-\bm{\hat{\phi}}}\rvert$
is the row index of the only non-zero entry in ${{\bf{s}}_{\text{pair},p}}$.

By reordering the entries of $\bm{\hat{\omega}_{\tau}}$ with the pairing
result, the Doppler $\bm{\hat{\omega}}$ is obtained as

$\bm{\hat{\omega}}={({[{{\bf{s}}_{{\rm{pair}},1}},\cdots,{{\bf{s}}_{{\rm{pair}},{P}}}]})}\bm{\hat{\omega}_{\tau}}+{f_{1}}{\bm{\hat{k}_{\tau}}}.\\!\\!\\!\\!$
(79)

Now, $\bm{\hat{\tau}(t)}$ and $\bm{\hat{\omega}}$ are paired correctly and
associated with the corresponding paths. With the estimated parameters, the
future channel at time $t$ can be reconstructed as ${{\bf{\hat{H}}}_{u}}(t)$.

Algorithm 3 Proposed MDMP channel prediction method.

1: Estimate $\bm{\hat{\theta}}$, $\bm{\hat{\phi}}$, $\bm{\hat{\tau}(t)}$ and
$\bm{\hat{k}_{\tau}}$ with Algorithm 1;

2: Estimate $\bm{\hat{\theta}_{\omega}}$, $\bm{\hat{\phi}_{\omega}}$ and
$\bm{\hat{\omega}_{\tau_{p}}}$ with Algorithm 2;

3: for $p=1:{P}$ do

4: Find the pairing matrix
${[{{\bf{s}}_{{\rm{pair}},1}},\cdots,{{\bf{s}}_{{\rm{pair}},{P}}}]^{T}}$ by
solving the minimization problem of Eq. (78);

5: end for

6: Calculate the Doppler $\bm{\hat{\omega}}$ by Eq. (79);

7: Predict future channel ${{\bf{\hat{H}}}_{u}}(t)$ with the estimated
parameters;

Our proposed MDMP method is summarized in Algorithm 3. The step 3 - step 6 of
Algorithm 3 has a time complexity order of $\mathcal{O}({P^{2}})$, and a
memory complexity order of $\mathcal{S}({P^{2}})$. The time complexity of
channel reconstruction in step 7 is
$\mathcal{O}({N_{h}^{{\rm{BS}}}}{N_{v}^{{\rm{BS}}}}{N_{f}}{P})$. The memory
complexity of step 7 is
$\mathcal{S}({N_{h}^{{\rm{BS}}}}{N_{v}^{{\rm{BS}}}}{N_{f}})$. The time and
memory complexity of the MDMP method are primarily determined by Algorithm 1
and Algorithm 2. Therefore, the global time complexity of the MDMP method is
$\mathcal{O}(2{\mu_{1}^{2}}{\mu_{2}})+\mathcal{O}(4{\mu_{1}}{\mu_{2}^{2}})+\mathcal{O}({\mu_{1}^{3}})+\mathcal{O}(2{\omega_{{\mu_{1}}}^{2}}{\omega_{{\mu_{2}}}})+\mathcal{O}(4{\omega_{{\mu_{1}}}}{\omega_{{\mu_{2}}}^{2}})+\mathcal{O}({\omega_{{\mu_{1}}}^{3}})$.
The MDMP method has a global memory complexity order of
$\max(\mathcal{S}(\mu_{1}^{2})+\mathcal{S}(4\mu_{2}^{2}),\mathcal{S}({\omega_{{\mu_{1}}}^{2}})+\mathcal{S}(4{\omega_{{\mu_{2}}}^{2}}))$.
The PAD method in [5] has a complexity order of
$\mathcal{O}(N{N_{t}}{N_{f}}\log({N_{t}}{N_{f}}))+\mathcal{O}({N_{p}}{N^{2.37}})+\mathcal{O}({N_{d}}{N_{p}}N)$,
where ${N_{p}}\ll{N_{t}}{N_{f}}$, $N$ and ${N_{d}}$ denote the prediction
order of the PAD method and the number of the predicted samples. The MDMP
method may have a larger time complexity compared to the PAD method. However,
the MDMP method has better performance and is applicable in more general
settings, i.e., the CSI delay does not have to be an interval multiple of the
pilot interval, the historical channel samples are not necessarily neighboring
ones. The detailed proofs will be shown in the next section.

##  IV Performance Analysis of the MDMP Prediction Method

In this section, we show the performance analysis of our proposed MDMP
prediction method. Denote the observation sample at time $t$ by
${{{\bf{\tilde{H}}}}_{u}}(t)$. The vectorized forms of ${{\bf{H}}_{u}}(t)$,
${{{\bf{\tilde{H}}}}_{u}}(t)$ and ${{\bf{\hat{H}}}_{u}}(t)$ are denoted
respectively by ${{\bf{h}}_{u}}(t)$, ${{{{\bf{\tilde{h}}}}_{u}}(t)}$, and
${{{{\bf{\hat{h}}}}_{u}}(t)}$, i.e.,
${{\bf{h}}_{u}}(t)\\!=\\!{\rm{vec}}({{\bf{H}}_{u}}(t))$,
${{{\bf{\tilde{h}}}}_{u}}(t)={\rm{vec(}}{{{\bf{\tilde{H}}}}_{u}}(t))$ and
${{\bf{\hat{h}}}_{u}}(t)={\rm{vec(}}{{{\bf{\hat{H}}}}_{u}}(t))$. Before the
asymptotic performance analysis, we introduce an assumption.

######  Assumption 1.

The observation sample yields

${{\bf{\tilde{h}}}_{u}}(t)={{\bf{h}}_{u}}(t)+{\bf{n}},\\!\\!\\!\\!$ (80)

where
${\bf{n}}={[{n_{1}},\cdots,{n_{N_{h}^{{\rm{BS}}}N_{v}^{{\rm{BS}}}{N_{f}}}}]^{T}}\in{{\mathbb{C}}^{N_{h}^{{\rm{BS}}}N_{v}^{{\rm{BS}}}{N_{f}}\times
1}}$ is the independent identically distributed (i.i.d.) Gaussian white noise
with zero mean and element-wise variance ${\sigma^{2}}$. As
$N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty$, the variance
${\sigma^{2}}$ converges to zero, such that:

$\mathop{\lim}\limits_{N_{v}^{{\rm{BS}}},N_{h}^{{\rm{BS}}},{N_{f}}\to\infty}\frac{{\left\|{{{{\bf{\tilde{h}}}}_{u}}(t)-{{\bf{h}}_{u}}(t)}\right\|_{2}^{2}}}{{\left\|{{{\bf{h}}_{u}}(t)}\right\|_{2}^{2}}}=0.\\!\\!$
(81)

Remarks: This technical assumption means the normalized channel sample error
converges to zero when the number of the BS antennas and the bandwidth
increase. In fact, the condition of Eq. (81) can be achieved even in the
multi-user multi-cell scenario with pilot contamination, by some non-linear
signal processing technologies [29].

For simplicity, we also introduce a vector in Definition 1.

######  Definition 1.

Define a vector containing the parameters of multipaths as

${\bf{\Omega}}(t)={\left[{{{\bf{\Omega}}_{1}}(t),\cdots,{{\bf{\Omega}}_{P}}(t)}\right]^{T}},\\!\\!\\!\\!$
(82)

where
${{\bf{\Omega}}_{p}}(t)=\left[{\cos{\theta_{p}}\sin{\phi_{p}},\sin{\theta_{p}},{\omega_{p}},{\tau_{p}}(t)}\right]$,
$p=1,\cdots,P$.

######  Theorem 1.

Under Assumption 1, for an arbitrary CSI delay ${t_{\tau}}$, the asymptotic
performance of MDMP prediction method yields:

$\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}\\!\\!\\!\\!\\!\\!\frac{{\left\|{{{{{\bf{\hat{h}}}}_{u}}(t+{t_{\tau}})}\\!-\\!{{{\bf{h}}_{u}}(t+{t_{\tau}})}}\right\|_{2}^{2}}}{{\left\|{{{{\bf{h}}_{u}}(t+{t_{\tau}})}}\right\|_{2}^{2}}}=0,\\!\\!$
(83)

providing that the pencil sizes satisfy $N_{h}^{{\rm{BS}}}-P+1>L>P$,
$N_{v}^{{\rm{BS}}}-P+1>R>P$, ${N_{f}}-P+1>K>P$, and ${N_{s}}-P+1>Q>P$.

###### Proof.

The proof can be found in Appendix A.

Remarks: Theorem 1 indicates that the channel prediction error converges to
zero when the number of the BS antennas and the bandwidth are large. Note that
the CSI delay in Theorem 1 does not have to be an integral multiple of the
pilot interval, which indicates our proposed MDMP method is more general than
the PAD method in [5] because the PAD achieves asymptotically error-free
performance when the CSI delay is an integral multiple of the channel sampling
interval.

If $Q=P+1$, the least number of samples in Theorem 1 should satisfy
${N_{s}}=2P+1$. The number of paths $P$ might be large in rich scattering
environment. In such cases, Theorem 1 indicates that we may need a large
number of samples $N_{s}$ to achieve asymptotically error-free prediction. In
the following, we will show the asymptotic performance of the prediction error
with only two arbitrary samples known.

######  Theorem 2.

Under Assumption 1, if the EOD and the AOD satisfy
$({\theta_{p}},{\phi_{p}})\neq({\theta_{q}},{\phi_{q}})$, $\forall p\neq q$,
the number of subcarriers satisfies ${N_{f}}\geq 2$, and the configuration of
the BS antennas satisfies $N_{h}^{{\rm{BS}}}-P+1>L>P$ and
$N_{v}^{{\rm{BS}}}-P+1>R>P$, then with only two arbitrary samples
${{{\bf{\tilde{h}}}}_{u}}({n_{{s_{1}}}}{T})$ and
${{{\bf{\tilde{h}}}}_{u}}({n_{{s_{2}}}}{T})$ known, the asymptotic performance
of the MDMP prediction method yields:

$\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}}\to\infty}\frac{{\left\|{{{{{\bf{\hat{h}}}}_{u}}(t+{t_{\tau}})}\\!-\\!{{{\bf{h}}_{u}}(t+{t_{\tau}})}}\right\|_{2}^{2}}}{{\left\|{{{{\bf{h}}_{u}}(t+{t_{\tau}})}}\right\|_{2}^{2}}}=0.\\!\\!\\!\\!$
(84)

###### Proof.

The proof can be found in Appendix B.

Remarks: Theorem 2 requires that the angles of any two paths are different,
which is a stronger assumption than in Theorem 1. However a better result is
obtained since only two samples are needed. Note that the two samples are not
necessarily neighboring ones, which is also more general than the PAD method.

In Theorem 1 and Theorem 2, the conditions of the BS antennas, e.g.,
$N_{h}^{{\rm{BS}}}-P+1>L>P$ and $N_{v}^{{\rm{BS}}}-P+1>R>P$, ensure the rank
of the 3-D MP matrices satisfy
$r({{{\bf{G}}}_{{\rm{re}}}}(t+{t_{\tau}}))=r({{\bm{\mathcal{G}}}_{u}}({f_{1}}))=P$.
If $L=P+1$ and $R=P+1$, the least number of the BS antennas is
${N_{t}}=(2P+1)^{2}$, which seems to be large. Next, by assuming
$r({{{\bf{G}}}_{{\rm{re}}}}(t+{t_{\tau}}))=r({{\bm{\mathcal{G}}}_{u}}({f_{1}}))=P$
known, we derive a lower-bound of the BS antennas to give a satisfactory
performance. In order to do so, we introduce two functions.

######  Definition 2.

Define two functions ${F_{1}}(L,R)$ and ${F_{2}}(R)$.

${F_{1}}(L,R)=\frac{{\cal{N}}_{s}}{{N_{v}^{{\rm{BS}}}-R+1}}+L-1,\\!\\!\\!\\!$
(85)

where $L$ and $R$ satisfy

$\left\\{\begin{array}[]{l}LR\geq{\cal Q}\\\ L\geq 2\\\ N_{v}^{{\rm{BS}}}\geq
R\geq 2\end{array}\right.,\\!\\!\\!\\!$ (86)

${\cal{N}}_{s}=\max(\frac{Q}{{{N_{s}}-Q+1}},1)$ and ${\cal
Q}=\max(\frac{P}{{Q-1}},4)$. The other function ${F_{2}}(R)$ is defined as

${F_{2}}(R)\\!=\\!\frac{P}{{({N_{s}}\\!-\\!Q\\!+\\!1)(N_{v}^{{\rm{BS}}}\\!-\\!R\\!+\\!1)}}\\!+\\!\frac{P}{{(R\\!-\\!1)Q}}\\!-\\!1,\\!\\!\\!\\!$
(87)

with

$\left\\{\begin{array}[]{l}LQ\geq{{\max\left({\frac{P}{{R-1}},4}\right)}}\\\
N_{v}^{{\rm{BS}}}\geq R\geq 2\end{array}\right..\\!\\!\\!\\!$ (88)

Next, we derive a lower-bound of the number of the BS antennas in Proposition
1 to achieve a satisfactory performance. For an arbitrary given
$N_{v}^{{\rm{BS}}}$, we derive the bound of the total number of the BS
antennas $N_{t}$. In a wideband system, ${N_{f}}$ and ${N_{s}}$ usually
satisfy ${N_{f}}>{N_{s}}$. According to Eq. (39) and Eq. (67), the lower-bound
of the number of the BS antennas is determined by ${N_{s}}$ in Eq. (67). The
lower-bound of $N_{t}$ consists of two sub-bounds
${f_{N_{v}^{{\rm{BS}}},{N_{s}},1}}$ and ${f_{N_{v}^{{\rm{BS}}},{N_{s}},2}}$,
where ${f_{N_{v}^{{\rm{BS}}},{N_{s}},1}}$ is derived from the 3-D MP matrix in
estimating Doppler, and ${f_{N_{v}^{{\rm{BS}}},{N_{s}},2}}$ is derived from
the 3-D MP matrix in estimating angles.

######  Proposition 1.

In a wideband channel, for a given ${N_{v}^{{\rm{BS}}}}$, a lower-bound of the
number of the BS antennas $N_{t}$ is given as

$N_{t}\geq\max({f_{N_{v}^{{\rm{BS}}},{N_{s}},1}},{f_{N_{v}^{{\rm{BS}}},{N_{s}},2}}),\\!\\!\\!\\!$
(89)

where

$\begin{array}[]{l}\\!\\!\\!\\!\\!{f_{N_{v}^{{\rm{BS}}},{N_{s}},1}}=\\!\\!\max(\\!{N_{v}^{{\rm{BS}}}}{F_{1}}({L_{1}},\\!{R_{1}}),\\!{N_{v}^{{\rm{BS}}}}{F_{1}}({L_{2}},\\!{R_{2}}),\\\
\\!{N_{v}^{{\rm{BS}}}}{F_{1}}({L_{3}},\\!{R_{3}}),{N_{v}^{{\rm{BS}}}}{F_{1}}({L_{4}},\\!{R_{4}}),\\!{N_{v}^{{\rm{BS}}}}{F_{1}}({L_{5}},\\!{R_{5}})),\\!\\!\\!\end{array}\\!\\!\\!\\!$
(90)

and $({L_{{n_{1}}}},{R_{{n_{1}}}})$, ${n_{1}}=1,2,3,4,5$ are defined as

$\begin{array}[]{l}{R_{1}}=\frac{{\sqrt{{\cal
Q}}}}{{\sqrt{{\cal{N}}_{s}}+\sqrt{{\cal Q}}}}(N_{v}^{{\rm{BS}}}+1)\\\
{L_{1}}=\frac{{\sqrt{{\cal Q}}(\sqrt{{\cal{N}}_{s}}+\sqrt{{\cal
Q}})}}{{N_{v}^{{\rm{BS}}}+1}}\end{array},\\!\\!\\!\\!$ (91)
$\begin{array}[]{l}{R_{2}}=\frac{{\sqrt{{\cal Q}}}}{{\sqrt{{\cal
Q}}-\sqrt{{\cal{N}}_{s}}}}(N_{v}^{{\rm{BS}}}+1)\\\ {L_{2}}=\frac{{\sqrt{{\cal
Q}}(\sqrt{{\cal
Q}}-\sqrt{{\cal{N}}_{s}})}}{{N_{v}^{{\rm{BS}}}+1}}\end{array},\\!\\!\\!\\!$
(92)

and $({R_{3}}=\frac{{{\cal Q}}}{2},{L_{3}}=2)$,
$({R_{4}}=2,{L_{4}}=\frac{{{\cal Q}}}{2})$, and
$({R_{5}}=N_{v}^{{\rm{BS}}},{L_{5}}=\frac{{{\cal Q}}}{N_{v}^{{\rm{BS}}}})$.
The value ${f_{N_{v}^{{\rm{BS}}},{N_{s}},2}}$ is defined as

$\begin{array}[]{l}\\!\\!{f_{N_{v}^{{\rm{BS}}},{N_{s}},2}}\\!\\!=\\\
\\!\\!\left\\{\begin{array}[]{l}\\!\\!\\!\\!\\!{N_{v}^{{\rm{BS}}}}(\frac{1}{{({N_{s}}-1)}}+1),\
\ \\!\\!\\!\\!\\!{\rm{if}}\
\\!N_{v}^{{\rm{BS}}}\\!-\\!P\\!+\\!1\\!\geq\\!R\\!\geq\\!\frac{P}{4}\\!+\\!1\\\
\\!\\!\\!\\!\\!{N_{v}^{{\rm{BS}}}}(\frac{1}{{({N_{s}}-Q+1)}}\\!+\\!\frac{P}{Q}\\!-\\!1),\\!{\rm{if}}\
\\!\min(N_{v}^{{\rm{BS}}}\\!-\\!\\!P\\!\\!+\\!1,\\!\frac{P}{4}\\!\\!+\\!1)\\!\geq\\!R\\!\geq\\!2\\\
\\!\\!\\!\\!\\!{N_{v}^{{\rm{BS}}}}(\frac{P}{{({N_{s}}-1)}}+1),\\!{\rm{if}}\
\\!N_{v}^{{\rm{BS}}}\\!\\!\geq\\!\\!R\\!\\!\geq\\!\\!\max(N_{v}^{{\rm{BS}}}\\!\\!-\\!P\\!\\!+\\!1,\frac{P}{4}\\!+\\!\\!1,2)\\\
\\!\\!\\!\\!\\!{N_{v}^{{\rm{BS}}}}{F_{2,\rm{max}}}(R),\\!{\rm{if}}\
\\!P\\!>\\!R\\!\geq\\!\max(N_{v}^{{\rm{BS}}}\\!-\\!P\\!+\\!1,2)\end{array}\right.\\!\\!\\!\\!\\!\\!,\end{array}\\!\\!\\!\\!\\!$
(93)

where

$\begin{array}[]{l}\\!\\!{F_{2,\rm{max}}}(R)\\!\\!=\\\
\\!\\!\left\\{\begin{array}[]{l}\\!\\!\\!\\!\max(\\!{F_{2}}(\max(2,\\!{R_{6}})),\\!{F_{2}}(\max(2,\\!{R_{7}})),\\\
\\!{F_{2}}(N_{v}^{{\rm{BS}}})),\;\\!{\rm{if}}\;Q<\frac{{{N_{s}}+1}}{2}\\\
\\!\\!\\!\\!\max({F_{2}}(\max(2,{R_{6}})),{F_{2}}(\max(2,{R_{7}}))),{\rm{if}}\;Q>\frac{{{N_{s}}+1}}{2}\\\
\\!\\!\\!\\!{F_{2}}({\frac{(N_{v}^{{\rm{BS}}}+1)}{2}}),\;{\rm{if}}\
\\!Q=\frac{{{N_{s}}+1}}{2}\end{array}\right.\\!,\end{array}\\!\\!\\!\\!\\!$
(94)

and

${R_{6}}=\frac{{\sqrt{{N_{s}}-Q+1}(N_{v}^{{\rm{BS}}}+1)}}{{\sqrt{{N_{s}}-Q+1}+\sqrt{Q}}},\\!\\!\\!\\!$
(95)

${R_{7}}=\frac{{\sqrt{{N_{s}}-Q+1}(N_{v}^{{\rm{BS}}}+1)}}{{\sqrt{{N_{s}}-Q+1}-\sqrt{Q}}}.\\!\\!\\!\\!$
(96)

###### Proof.

The proof can be found in Appendix C.

Remarks: In Proposition 1, for a given $N_{v}^{{\rm{BS}}}$, we derive the
lower-bound of $N_{t}$. Such a bound means that the MDMP method does not work
if the BS antenna configuration cannot satisfy the lower-bound. As a remedy,
the lower-bound can be made smaller to satisfy the antenna configuration, by
increasing the number of channel samples $N_{s}$. Note that if
$N_{h}^{{\rm{BS}}}$ is given, the lower-bound of $N_{t}$ can also be derived.
The detailed derivation is omitted.

##  V Numerical Results

In this section, the simulation results of our proposed method are shown. The
clustered delay line (CDL) channel model of 3GPP is adopted. The number of
scattering clusters is 9. Each cluster contains 20 rays and the total number
of propagation paths is 180. According to [27], the $p$-th path delay is
modeled as
${\tau_{p}}(t)=\tau{{}^{\prime}_{p,0}}-\min\left\\{{\tau{{}^{\prime}_{1,0}},\cdots,\tau{{}^{\prime}_{P,0}}}\right\\}+\tau^{\prime}{{}^{\prime}_{p}}(t)$,
where $\tau^{\prime}{{}^{\prime}_{p}}(t)$ is the cluster delay, and
$\tau{{}^{\prime}_{p,0}}-\min\left\\{{\tau{{}^{\prime}_{1,0}},\cdots,\tau{{}^{\prime}_{P,0}}}\right\\}$
denotes the initial delay and is modeled as a spatially random variable
related to the correlation distance. The root mean square (RMS) angular
spreads of AOD, EOD, AOA and EOA are $87.1^{\circ}$, $56.4^{\circ}$,
$102.1^{\circ}$ and $65.3^{\circ}$, respectively. The main simulation
parameters are listed in Table I. Considering the 3-D Urban Macro (3-D UMa)
scenario, the UEs have certain velocity, e.g., 60 km/h and 120 km/h. The
central frequency is 3.5 GHz and the bandwidth is 100 MHz, which is composed
of 273 resource blocks (RBs) in frequency domain. The duration of a slot is
0.5 ms, which contains 14 OFDM symbols. One channel sample is available for
each slot.

TABLE I: The main simulation parameters. Scenario |  3D Urban Macro (3D UMa)
---|---
Carrier frequency (GHz) |  3.5
Bandwidth (MHz) |  100 (273 RBs)
Subcarrier spacing (kHz) |  30
Number of UEs |  16
BS antenna configuration |  $(\underline{M},\underline{N},\underline{P},\underline{M}_{g},\underline{N}_{g})=(8,8,1,1,1)$, $(d_{h}^{{\rm{tx}}},d_{v}^{{\rm{tx}}})=(0.5\lambda,0.8\lambda)$
UE antenna configuration |  $(\underline{M},\underline{N},\underline{P},\underline{M}_{g},\underline{N}_{g})=(1,1,2,1,1)$, the polarization angles are $0^{\circ}$ and $90^{\circ}$
Delay Spread (ns) |  616
CSI delay (ms) |  12, 16
UEs speed (km/h) |  60, 120

The antenna configuration is denoted by
$(\underline{M},\underline{N},\underline{P},\underline{M}_{g},\underline{N}_{g})$,
where each antenna panel has $\underline{M}$ rows and $\underline{N}$ columns
of antenna elements; $\underline{P}$ is the number of polarizations;
$\underline{M}_{g}$ is the number of panels in a row and $\underline{N}_{g}$
is the number of panels in a column. In the horizontal and vertical
directions, the antenna spacings are $0.5\lambda$ and $0.8\lambda$
respectively. The pencil sizes such as $L$, $R$, $K$ and $Q$, are set as 6, 5,
137 and 15, respectively. The DL precoder is eigen-based zero-forcing (EZF)
[30]. Two performance metrics are evaluated, which are the DL spectral
efficiency and the prediction error. The DL spectral efficiency is calculated
by
$\sum\limits_{u=1}^{N_{\text{UE}}}{E\left\\{{{{\log}_{2}}(1+{\rm{SINR}}_{u})}\right\\}}$,
where $N_{\text{UE}}$ is the number of UEs, ${\rm{SINR}}_{u}$ is the signal-
to-noise ratio of the $u$-th UE, and the expectation is taken over frequency
and time. The DL prediction error is defined by
$10\log\left\\{{E\left\\{{\frac{{\left\|{{{{\bf{\hat{H}}}}_{u}}-{{\bf{H}}_{u}}}\right\|_{F}^{2}}}{{\left\|{{{\bf{H}}_{u}}}\right\|_{F}^{2}}}}\right\\}}\right\\}$,
where ${{{\bf{\hat{H}}}}_{u}}$ and ${{\bf{H}}_{u}}$ are the predicted and
exact channels, respectively, and the expectation is taken over time,
frequency and UEs.

Figure 2: (a) The spectral efficiency versus SNR and (b) the prediction error
versus the CSI delay, the UEs move at 120 km/h, the BS has 64 antennas.

Let CDL-TID and CDL-TVD denote the simulated CDL model with time-invariant and
time-varying path delay, respectively. In Fig. 2(a) and (b), we show the
spectral efficiencies and prediction errors of different methods with the UEs
moving at the speed of 120 km/h. In Fig. 2(a), the CSI delay is relatively
large, e.g., 12 ms and 16 ms, while in Fig. 2(b), the CSI delay is within the
range of 4 to 20 ms. The ideal case of the stationary setting is also shown as
a reference curve labeled as “Stationary channel”, which is the upper-bound of
the performance. The curves labeled with “No prediction” mean channel
predictions are not carried out. In Fig. 2(a), the different performances of
the PAD method in CDL-TID and CDL-TVD show that time-varying delay brings an
observable decrease of the spectral efficiency. However, our proposed MDMP
method is able to deal with the time-varying path delay, as it approaches the
ideal case of the stationary setting in the scenario with a large CSI delay of
16 ms and high velocity level of 120 km/h.

Figure 3: The spectral efficiency versus SNR, 16 ms of CSI delay, the BS has
64 antennas. Figure 4: The spectral efficiency versus SNR, 16 ms of CSI delay,
the BS equipped with 64 antennas, multiple velocity levels of UEs, i.e., four
at 30 km/h, four at 60 km/h, four at 90 km/h and four at 120 km/h.

Fig. 3 compares the spectral efficiencies of different methods with UEs moving
at 60 km/h and 120 km/h, respectively. The curves labeled as “APF-RNS”, are
the prediction performance of an adaptive and parameter-free recurrent neural
structure (APF-RNS) [31], which uses 19 historical sequential channel samples
to predict the channel at the next moment. In CDL-TVD, it can be observed that
the performance of PAD method will decrease more obviously at higher velocity
because of the faster varying delay and Doppler. However, our proposed MDMP
method is close to the upper-bound of the stationary setting even in high-
mobility scenarios.

Fig. 4 shows the spectral efficiencies of the MDMP and PAD methods with
multiple-speed UEs, where every four UEs move at 30 km/h, 60 km/h, 90 km/h and
120 km/h, respectively. One may observe that our proposed MDMP method still
performs well in this setting.

Figure 5: The prediction error versus the numbers of BS antennas, UEs move at
120 km/h, 16 ms of CSI delay.

Fig. 5 gives the prediction errors of MDMP and PAD methods with different
numbers of the BS antennas, e.g., 64, 128, 192, 256, 320 and 384. The
corresponding antenna configurations are $(8,8,1,1,1)$, $(8,8,1,1,2)$,
$(8,8,1,1,3)$, $(8,8,1,1,4)$, $(8,8,1,1,5)$ and $(8,8,1,1,6)$, respectively.
One may notice that the prediction error of our MDMP method decreases
obviously with the increasing numbers of the BS antennas. One may notice the
error floor of the MDMP method. This is because the bandwith is fixed and not
large enough.

Figure 6: The spectral efficiency versus SNR, UEs move at 120 km/h, the BS
equipped with 64 antennas.

By adding a line-of-sight (LOS) path, CDL-TVD is changed to a new model named
CDL-LOSTVD, which contains 181 propagation paths. The RMS angular spreads of
AOD, EOD, AOA and EOA are updated to $47.9^{\circ}$, $82.9^{\circ}$,
$89.9^{\circ}$ and $84.6^{\circ}$, respectively. The spectral efficiencies of
different methods in CDL-LOSTVD are shown in Fig. 6. It is clear that our
proposed MDMP method is capable of dealing with the effect of time-varying
delay in this setting, and nearly approaches the upper-bound of the stationary
setting.

##  VI Conclusion

In this paper, we addressed the challenge of mobility under the industrial
channel model with time-varying path delay in massive MIMO. We proposed a
novel multi-dimension matrix pencil channel prediction method, which estimates
the angles, delay and Doppler simultaneously by exploiting the angular-time-
domain and angular-frequency-domain structures of the wideband massive MIMO
channel. We also proposed a pairing method to associate the paths with the
estimated angles, delay and Doppler by exploiting the super-resolution
property of the angle estimation. In the theoretical analysis, we proved that
the proposed MDMP method asymptotically achieves the error-free performance
with arbitrary CSI delay, providing that only two arbitrary samples are known.
We also derived a lower-bound of the number of the BS antennas for this method
to work. By increasing the number of channel samples, the bound can be made
lower. Simulation results demonstrated that our proposed MDMP method in high-
mobility scenarios with a large CSI delay, is very close to the performance of
the stationary setting.

###  -A Proof of Theorem 1

For ease of exposition, we will first calculate the estimations of parameters
with the observation channel samples, and then we prove that the estimations
of parameters converge to the exact ones.

Firstly, the channel parameters are estimated as follows: After an arbitrary
CSI delay ${t_{\tau}}$, $t$ is updated as $t+{t_{\tau}}$. The channel in Eq.
(13) is rewritten as

${{\bf{\tilde{H}}}_{u}}(t+{t_{\tau}})={{\bf{A}}_{u}}{{\bf{C}}_{u}}{{\bf{\bar{C}}}_{u}}\left({{{\bf{\bar{B}}}_{u}}\odot{{\bf{B}}_{u}}}\right)+{\bf{N}},\\!\\!\\!\\!$
(97)

where ${\bf{N}}$ is the sample noise matrix. The matrix ${{\bf{\bar{C}}}_{u}}$
is introduced by UE movement during ${t_{\tau}}$, and is expressed as

${{\bf{\bar{C}}}_{u}}={\rm{diag}}\left\\{{\begin{array}[]{*{20}{c}}{e^{j2\pi{\omega_{1}}{t_{\tau}}},}&{\cdots,}&{e^{j2\pi{\omega_{P}}{t_{\tau}}}}\end{array}}\right\\}.\\!\\!\\!\\!$
(98)

The matrix ${{\bf{\bar{B}}}_{u}}$ contains multipath delay response vectors
and is given by

${{\bf{\bar{B}}}_{u}}=\left[{\begin{array}[]{*{20}{c}}{{\bf{\bar{b}}}({f_{1}}),}&{\cdots,}&{{\bf{\bar{b}}}({f_{{N_{f}}}})}\end{array}}\right],\\!\\!\\!\\!$
(99)

where

${\bf{\bar{b}}}({f_{n_{f}}})={\left[{\begin{array}[]{*{20}{c}}{e^{-j2\pi{f_{n_{f}}}{k_{{\tau_{1}}}}{t_{\tau}}},}&{\cdots,}&{e^{-j2\pi{f_{n_{f}}}{k_{{\tau_{P}}}}{t_{\tau}}}}\end{array}}\right]^{T}}\\!\\!.\\!\\!\\!\\!$
(100)

Then, the 3-D MP matrix is expressed as

$\\!{{\bf{\tilde{G}}}_{u}}(t+{t_{\tau}})\\!=\\!{{\check{\bf{E}}}_{2}}(t+{t_{\tau}}){\bf{Y}}(t+{t_{\tau}}){{\check{\bf{F}}}_{2}}(t+{t_{\tau}})+{\bf{N}_{{\bf{\tilde{G}}}_{u}}},\\!\\!\\!\\!$
(101)

where
${\bf{N}_{{\bf{\tilde{G}}}_{u}}}\in{{\mathbb{C}}^{{\mu_{1}}\times{\mu_{2}}}}$
is generated by ${\bf{N}}$, and

${\bf{Y}}(t+{t_{\tau}})={\bf{Y}}{{\bf{\bar{C}}}_{u}}{\rm{diag}}\left\\{{{\bf{\bar{b}}}({f_{1}})}\right\\},\\!\\!\\!\\!$
(102) $\ \ \ \ \ \ \ \ \
{{\check{\bf{E}}}_{2}}(t+{t_{\tau}})={\left[\\!{\begin{array}[]{*{20}{c}}{({{\check{\bf{E}}}_{1}})^{T},}\\!\\!&\\!{\cdots,}\\!\\!&\\!{({{\check{\bf{E}}}_{1}}{{\bf{Z}}_{{t_{\tau}}}^{K-1}})^{T}}\end{array}}\\!\right]^{T}}\\!,\\!\\!\\!\\!$
(103)
$\\!{{\check{\bf{F}}}_{2}}(t+{t_{\tau}})=\left[\\!{\begin{array}[]{*{20}{c}}{{\check{\bf{F}}}_{1}}\\!&\\!{\cdots,}\\!&\\!{{{\bf{Z}}_{{t_{\tau}}}^{{N_{f}}-K}}{{\check{\bf{F}}}_{1}}}\end{array}}\\!\right],\\!\\!\\!\\!$
(104)

and

${{\bf{Z}}_{t_{\tau}}}={{\bf{Z}}_{\tau}}{\rm{diag}}\left\\{{{\bf{\bar{b}}}({f_{{n_{f}}+1}})}\right\\}{\left({{\rm{diag}}\left\\{{{\bf{\bar{b}}}({f_{n_{f}}})}\right\\}}\right)^{*}}.\\!\\!\\!\\!$
(105)

As $N_{h}^{{\rm{BS}}}-P+1>L>P$, $N_{v}^{{\rm{BS}}}-P+1>R>P$ and
${N_{f}}-P+1>K>P$, $r({{\bf{\tilde{G}}}_{u}}(t+{t_{\tau}}))\geq P$ and
${{\bf{\tilde{G}}}_{u}}(t+{t_{\tau}})$ contains the angular and path delay
information of the paths. The SVD of ${{\bf{\tilde{G}}}_{u}}(t+{t_{\tau}})$ is
${{\bf{\tilde{G}}}_{u}}(t+{t_{\tau}})={{\bf{\tilde{U}}}_{u,t_{\tau}}}{{\bf{\tilde{S}}}_{u,t_{\tau}}}{\bf{\tilde{V}}}_{u,t_{\tau}}^{H}$,
where
${{\bf{\tilde{U}}}_{u,{t_{\tau}}}}=\left[{{{\bf{\tilde{u}}}_{u,{t_{\tau}},1}},\cdots,{{\bf{\tilde{u}}}_{u,{t_{\tau}},{M_{\tau}}}}}\right]$,
${{\bf{\tilde{S}}}_{u,{t_{\tau}}}}={\rm{diag}}\\{{{\tilde{s}}_{u,{t_{\tau}},1}},\cdots,{{\tilde{s}}_{u,{t_{\tau}},{M_{\tau}}}}\\}$,
and the diagonal elements of ${{\bf{\tilde{S}}}_{u,{t_{\tau}}}}$ are in non-
increasing order. Without loss of generality, we let ${\mu_{1}}<{\mu_{2}}$,
and ${M_{\tau}}={\mu_{1}}$. We set a new threshold ${\gamma_{2}}$ to determine
the number of the non-negligible paths as $P_{1}\leq P$. Define the
corresponding $P_{1}$ columns of ${{\bf{\tilde{U}}}_{u,{t_{\tau}}}}$ as
${{\bf{\tilde{U}}}_{u,s,t_{\tau}}}$.

We select $P_{1}$ columns from ${{\check{\bf{E}}}_{2}}(t+{t_{\tau}})$ to form
${{\check{\bf{E}}}_{s,2}}(t+{t_{\tau}})$. Since
$r({{\bf{\tilde{U}}}_{u,s,t_{\tau}}})=r({{\check{\bf{E}}}_{s,2}}(t+{t_{\tau}}))={P_{1}}$,
they may satisfy a mapping relationship as
${{\bf{\tilde{U}}}_{u,s,t_{\tau}}}={{\check{\bf{E}}}_{s,2}}(t+{t_{\tau}}){\bf{T}}_{1}$,
where ${\bf{T}}_{1}\in{{\mathbb{C}}^{{P_{1}}\times{P_{1}}}}$ is a full-rank
matrix. Define the first $KRL-RL$ rows and the last $KRL-RL$ rows from
${{\bf{\tilde{U}}}_{u,s,t_{\tau}}}$ as ${{\bf{\tilde{U}}}_{u,s1}}$ and
${{\bf{\tilde{U}}}_{u,s2}}$, respectively, which satisfy

$\\!{{\bf{\tilde{U}}}_{u,s2}}\\!-\\!{{\bf{\tilde{U}}}_{u,s1}}{{\bm{\lambda}}_{t_{\tau}}}\\!\\!=\\!{{\bf{J}}_{1}}{{\check{\bf{E}}}_{s,2}}(t\\!+\\!{t_{\tau}})({{\bf{T}}_{1}^{-1}}{{\bf{Z}}_{t_{\tau}}}{\bf{T}}_{1}\\!-\\!{{\bm{\lambda}}_{t_{\tau}}})\\!=\\!{\bf{0}},\\!\\!\\!\\!$
(106)

where
${{\bm{\lambda}}_{t_{\tau}}}={{\bf{\tilde{U}}}_{u,s1}}^{\dagger}{{\bf{\tilde{U}}}_{u,s2}}={{\bf{T}}_{1}^{-1}}{{\bf{Z}}_{t_{\tau}}}{\bf{T}}_{1}$.
The matrix ${{\bm{\lambda}}_{t_{\tau}}}$ and the diagonal matrix
${{\bf{Z}}_{t_{\tau}}}$ are similarity matrices, and share the same
eigenvalues. With the UMP method, ${{\bf{\tilde{G}}}_{u}}(t+{t_{\tau}})$,
${{\bf{Z}}_{t_{\tau}}}$, and ${{\bm{\lambda}}_{t_{\tau}}}$ are transformed to
${{{\bf{\tilde{G}}}}_{{\rm{re}}}}(t+{t_{\tau}})$,
${{{\bf{\hat{Z}}}}_{t_{\tau}}}$ and ${{\bf{\tilde{\Psi}}}_{{t_{\tau}}}}$,
respectively, where ${{{\bf{\tilde{G}}}}_{{\rm{re}}}}(t+{t_{\tau}})$ is
defined as

$\begin{array}[]{l}{{{\bf{\tilde{G}}}}_{{\rm{re}}}}(t+{t_{\tau}})={{{\bf{G}}}_{{\rm{re}}}}(t+{t_{\tau}})+{{{\bf{N}}}_{{\rm{re}}}},\end{array}\\!\\!\\!\\!$
(107)

and ${{{\bf{N}}}_{{\rm{re}}}}$ is generated from
${\bf{N}_{{\bf{\tilde{G}}}_{u}}}$. The diagonal matrix
${{{\bf{\hat{Z}}}}_{t_{\tau}}}$ and the matrix
${{\bf{\tilde{\Psi}}}_{{t_{\tau}}}}$ are introduced by

$\begin{array}[]{l}{{{\bf{\hat{Z}}}}_{t_{\tau}}}={\rm{diag}}\\{\tan(\pi\Delta
f({\tau_{1}}(t+{t_{\tau}}))),\cdots,\\\ \ \ \ \ \ \ \ \tan(\pi\Delta
f({\tau_{{P_{1}}}}(t+{t_{\tau}})))\\},\end{array}\\!\\!\\!\\!$ (108)
${{\bf{\tilde{\Psi}}}_{{t_{\tau}}}}={\left[{{\mathop{\rm
Re}\nolimits}({\bf{Q}}_{\mu_{3}}^{H}{{\bf{J}}_{1}}{\bf{Q}}_{\mu_{1}}){{\bf{\tilde{U}}}_{s,t_{\tau}}}}\right]^{\dagger}}{\mathop{\rm
Im}\nolimits}({\bf{Q}}_{\mu_{3}}^{H}{{\bf{J}}_{1}}{\bf{Q}}_{\mu_{1}}){{\bf{\tilde{U}}}_{s,t_{\tau}}},\\!$
(109)

where ${{\bf{\tilde{U}}}_{s,t_{\tau}}}$ is the $P_{1}$ columns from
${{\bf{\tilde{U}}}_{t_{\tau}}}$, and ${{\bf{\tilde{U}}}_{t_{\tau}}}$ is
calculated by the SVD of ${{{\bf{\tilde{G}}}}_{{\rm{re}}}}(t+{t_{\tau}})$:

${{{\bf{\tilde{G}}}}_{{\rm{re}}}}(t+{t_{\tau}})={{\bf{\tilde{U}}}_{t_{\tau}}}{{\bf{\tilde{S}}}_{t_{\tau}}}{\bf{\tilde{V}}}_{t_{\tau}}^{H}={{\bf{U}}_{t_{\tau}}}{{\bf{S}}_{t_{\tau}}}{\bf{V}}_{t_{\tau}}^{H}+{{\bf{U}}_{{\bf{N}}_{\rm{re}}}}{{\bf{S}}_{{\bf{N}}_{\rm{re}}}}{\bf{V}}_{{\bf{N}}_{\rm{re}}}^{H},\\!$
(110)

where
${{\bf{\tilde{U}}}_{t_{\tau}}}=\left[{{{\bf{\tilde{u}}}_{{t_{\tau}},1}},\cdots,{{\bf{\tilde{u}}}_{{t_{\tau}},{\mu_{1}}}}}\right]$,
${{\bf{\tilde{S}}}_{{t_{\tau}}}}={\rm{diag}}\\{{{\tilde{s}}_{{t_{\tau}},1}},\cdots,{{\tilde{s}}_{{t_{\tau}},{\mu_{1}}}}\\}$,
${{\bf{U}}_{t_{\tau}}}=\left[{{{\bf{u}}_{{t_{\tau}},1}},\cdots,{{\bf{u}}_{{t_{\tau}},{\mu_{1}}}}}\right]$,
and
${{\bf{S}}_{{t_{\tau}}}}={\rm{diag}}\\{{{s}_{{t_{\tau}},1}},\cdots,{{s}_{{t_{\tau}},{\mu_{1}}}}\\}$.
Obviously, ${{\bf{\tilde{\Psi}}}_{{t_{\tau}}}}$ and
${{{\bf{\hat{Z}}}}_{t_{\tau}}}$ are similarity matrices. With the EVD of
${{\bf{\tilde{\Psi}}}_{{t_{\tau}}}}$, we may estimate
${{{\bf{\hat{Z}}}}_{t_{\tau}}}$ by
${{{\bf{\hat{Z}}}}_{t_{\tau}}}={\bf{W}}_{{t_{\tau}}}^{-1}{{\bf{\tilde{\Psi}}}_{{t_{\tau}}}}{{\bf{W}}_{{t_{\tau}}}}$,
where ${\bf{W}}_{{t_{\tau}}}$ contains the eigenvectors of
${{\bf{\tilde{\Psi}}}_{{t_{\tau}}}}$. After CSI delay ${t_{\tau}}$, the path
delay ${\hat{\tau}_{p}}(t+{t_{\tau}})$ is estimated as

$\begin{array}[]{l}{\hat{\tau}_{p}}(t+{t_{\tau}})=\frac{{{{\tan}^{-1}}\left({{{\hat{z}}_{t_{\tau},p}}(t+{t_{\tau}})}\right)}}{{\pi\Delta
f}},\end{array}\\!\\!\\!\\!$ (111)

where ${\hat{z}_{t_{\tau},p}}(t+{t_{\tau}})$ is the $p$-th estimated diagonal
element of ${{\bf{\hat{Z}}}_{t_{\tau}}}$. The changing rate of path delay
${\hat{k}_{{\tau_{p}}}}$ is estimated by two different samples as

${\hat{k}_{{\tau_{p}}}}=\frac{{{{\tan}^{-1}}\left({{{\hat{z}}_{t_{\tau},p}}({t_{2}}+{t_{\tau}})}\right)-{{\tan}^{-1}}\left({{{\hat{z}}_{t_{\tau},p}}({t_{1}}+{t_{\tau}})}\right)}}{{\pi\Delta
f({t_{2}}-{t_{1}})}}\\!.\\!\\!\\!\\!$ (112)

The other parameters, i.e., Doppler, EODs and AODs can also be estimated by
the similar estimation procedure of ${\hat{\tau}_{p}}(t+{t_{\tau}})$, and the
details are omitted. So far, the parameters have been estimated with the
observation samples. Then, we derive the asymptotic performance of parameter
estimations under Assumption 1.

When the number of the BS antennas and the bandwidth are large, the
correlation matrix of ${{{\bf{\tilde{G}}}}_{{\rm{re}}}}(t+{t_{\tau}})$ is
calculated by

$\begin{array}[]{l}\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}{{\bm{R}}_{{{{\bf{\tilde{G}}}}_{{\rm{re}}}}}}=E\\{{{{\bf{\tilde{G}}}}_{{\rm{re}}}}(t+{t_{\tau}}){{{{\bf{\tilde{G}}}}_{{\rm{re}}}}^{H}}(t+{t_{\tau}})\\}\\\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
={{\bm{R}}_{{{\bf{G}}_{{\rm{re}}}}}}+{\sigma^{2}}{{\bf{I}}_{{\mu_{1}}}},\end{array}\\!\\!\\!\\!$
(113)

where
${{\bm{R}}_{{{\bf{G}}_{{\rm{re}}}}}}=E\\{{{\bf{G}}_{{\rm{re}}}}(t+{t_{\tau}}){{\bf{G}}_{{\rm{re}}}^{H}}(t+{t_{\tau}})\\}$,
${\sigma^{2}}{{\bf{I}}_{{\mu_{1}}}}=E\\{{{\bf{N}}_{{\rm{re}}}}{\bf{N}}_{{\rm{re}}}^{H}\\}$,
and $E\\{\cdot\\}$ is the expectation over the BS antennas and bandwidth. By
performing the SVD of ${{\bm{R}}_{{{{\bf{\tilde{G}}}}_{{\rm{re}}}}}}$, we
obtain

$\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}\\!\\!\\!\\!\\!\\!\\!{{\bf{\tilde{U}}}_{{t_{\tau}}}}{{\bf{\tilde{\Sigma}}}_{{t_{\tau}}}}{\bf{\tilde{U}}}_{{t_{\tau}}}^{H}\\!=\\!{{\bf{U}}_{{t_{\tau}}}}{{\bf{\Sigma}}_{{t_{\tau}}}}{\bf{U}}_{{t_{\tau}}}^{H}\\!+\\!{{\bf{U}}_{{t_{\tau}}}}{{\bf{\Sigma}}_{{\bf{N}}_{\rm{re}}}}{\bf{U}}_{{t_{\tau}}}^{H},$
(114)

where
${{\bf{\Sigma}}_{{\bf{N}}_{\rm{re}}}}={\sigma^{2}}{{\bf{I}}_{{\mu_{1}}}}$ and

${{\bf{\tilde{\Sigma}}}_{{t_{\tau}}}}={\rm{diag}}\\{{\begin{array}[]{*{20}{c}}{\lvert{{\tilde{s}}_{{t_{\tau}},1}}\rvert^{2},}&{\cdots,}&{\lvert{{\tilde{s}}_{{t_{\tau}},{\mu_{1}}}}\rvert^{2}}\end{array}}\\},\\!\\!\\!\\!$
(115)
${{\bf{\Sigma}}_{{t_{\tau}}}}={\rm{diag}}\\{{\begin{array}[]{*{20}{c}}{\lvert{s_{{t_{\tau}},1}}\rvert^{2},}&{\cdots,}&{\lvert{s_{{t_{\tau}},{\mu_{1}}}}\rvert^{2}}\end{array}}\\}.\\!\\!\\!\\!$
(116)

Under Assumption 1, we may obtain

$\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}{{{\bf{\tilde{u}}}}_{{t_{\tau}},n}}={{\bf{u}}_{{t_{\tau}},n}}{e^{j{\vartheta_{n}}}},\\!\\!\\!\\!$
(117)

where ${\vartheta_{n}}\in[0,2\pi]$, and

$\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}\\!{{\bf{\tilde{\Sigma}}}_{{t_{\tau}}}}\\!=\\!\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}({{\bf{\Sigma}}_{{t_{\tau}}}}\\!+\\!{{\bf{\Sigma}}_{{\bf{N}}_{\rm{re}}}})={{\bf{\Sigma}}_{{t_{\tau}}}}.\\!\\!\\!\\!$
(118)

In other words, the $n$-th diagonal elements in
${{\bf{\tilde{\Sigma}}}_{{t_{\tau}}}}$ and ${{\bf{\Sigma}}_{{t_{\tau}}}}$
satisfy:
$\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}{\lvert{{\tilde{s}}_{{t_{\tau}},{n}}}\rvert^{2}}={\lvert{{s}_{{t_{\tau}},{n}}}\rvert^{2}}$.
According to Eq. (110), we may obtain
$\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}{{\bf{\tilde{S}}}_{t_{\tau}}}={{\bf{S}}_{t_{\tau}}}$
and
$\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}{P_{1}}=P$.
According to Eq. (110), we denote the $P$ columns of
${{{\bf{U}}}_{{t_{\tau}}}}$ as ${{{\bf{U}}}_{{s,t_{\tau}}}}$ and may obtain

$\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}{{{\bf{\tilde{U}}}}_{{s,t_{\tau}}}}={{{\bf{U}}}_{{s,t_{\tau}}}}{\bf{\Gamma}},\\!\\!\\!\\!$
(119)

where
${\bf{\Gamma}}=\left[{{{\bf{\chi}}_{1}},\cdots,{{\bf{\chi}}_{{P}}}}\right]\in{{\mathbb{C}}^{P\times
P}}$ and ${{\bf{\chi}}_{{n}}}$ is a ${P\times 1}$ unitary vector with the
$n$-th element being ${e^{j{\vartheta_{n}}}}$. The real matrix related to path
delays asymptotically converges to

$\begin{array}[]{l}\\!\\!\\!\\!\\!\\!\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}\\!\\!\\!\\!{{\bf{\tilde{\Psi}}}_{{t_{\tau}}}}\\!=\\!\\!\\!\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}\\!{\left[{{\rm{Re}}({\bf{Q}}_{\mu_{3}}^{H}{{\bf{J}}_{1}}{\bf{Q}}_{\mu_{1}}){{{\bf{\tilde{U}}}}_{s,{t_{\tau}}}}}\\!\right]^{\dagger}}\\\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
{\rm{Im}}({\bf{Q}}_{\mu_{3}}^{H}{{\bf{J}}_{1}}{\bf{Q}}_{\mu_{1}}){{{\bf{\tilde{U}}}}_{s,{t_{\tau}}}}\\\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
=\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}{\bf{\Gamma}}^{-1}{{\bf{\hat{\Psi}}}_{t_{\tau}}}{{\bf{\Gamma}}},\end{array}\\!$
(120)

where ${{\bf{\hat{\Psi}}}_{t_{\tau}}}$ is generated by the sample without
noise:

$\begin{array}[]{l}{{\bf{\hat{\Psi}}}_{t_{\tau}}}={\left[{{\rm{Re}}({\bf{Q}}_{\mu_{3}}^{H}{{\bf{J}}_{1}}{\bf{Q}}_{\mu_{1}}){{{\bf{U}}}_{s,{t_{\tau}}}}}\right]^{\dagger}}{\rm{Im}}({\bf{Q}}_{\mu_{3}}^{H}{{\bf{J}}_{1}}{\bf{Q}}_{\mu_{1}}){{{\bf{U}}}_{s,{t_{\tau}}}}\end{array}.\\!$
(121)

Obviously, ${{{\bf{\tilde{\Psi}}}}_{t_{\tau}}}$ and
${{\bf{\hat{\Psi}}}_{t_{\tau}}}$ are similar and share the same eigenvalues.
In other words, despite the observation noise, the estimation of the $p$-th
path delay ${\hat{\tau}_{p}}(t+{t_{\tau}})$ asymptotically converges to the
exact value, i.e.,

$\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}{\hat{\tau}_{p}}(t+{t_{\tau}})={\tau_{p}}(t+{t_{\tau}}).\\!\\!\\!\\!$
(122)

Likewise,

$\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}{{\bf{\hat{\Omega}}}(t+{t_{\tau}})}={{\bf{\Omega}}(t+{t_{\tau}})}.\\!\\!\\!\\!$
(123)

With the estimated parameters in ${{\bf{\hat{\Omega}}}(t+{t_{\tau}})}$, the
channel after ${t_{\tau}}$ is reconstructed as
${{\bf{\hat{H}}}_{u}}(t+{t_{\tau}})$. Then, we may obtain

$\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}},{N_{f}}\to\infty}\\!\\!\\!\\!\\!\\!\frac{{\left\|{{{{{\bf{\hat{h}}}}_{u}}(t+{t_{\tau}})}\\!-\\!{{{\bf{h}}_{u}}(t+{t_{\tau}})}}\right\|_{2}^{2}}}{{\left\|{{{{\bf{h}}_{u}}(t+{t_{\tau}})}}\right\|_{2}^{2}}}=0.\\!\\!$
(124)

Thus, Theorem 1 is proved. $\square$

###  -B Proof of Theorem 2

We first derive the least number of samples below. The 2-D MP matrix
${{{\bm{\tilde{\mathcal{G}}}}_{u}}({n_{s}},{f_{1}})}$ generated by the
observation samples, is defined as

${{{\bm{\tilde{\mathcal{G}}}}_{u}}({n_{s}},{f_{1}})}={{{\bm{\mathcal{G}}}_{u}}({n_{s}},{f_{1}})}+{{\bf{N}}_{{\bm{\mathcal{G}}},{n_{s}}}},\\!\\!\\!\\!$
(125)

where
${{\bf{N}}_{{\bm{\mathcal{G}}},{n_{s}}}}\in{{\mathbb{C}}^{LR\times(N_{h}^{{\rm{BS}}}-L+1)(N_{v}^{{\rm{BS}}}-R+1)}}$
is the noise matrix. The 3-D MP matrix
${{{\bm{\tilde{\mathcal{G}}}}_{u}}({f_{1}})}$ is introduced by

${{{\bm{\tilde{\mathcal{G}}}}_{u}}({f_{1}})}={{{\bm{\mathcal{G}}}_{u}}({f_{1}})}+{{\bf{N}}_{{\bm{\mathcal{G}}}}},\\!\\!\\!\\!$
(126)

where
${{\bf{N}}_{{\bm{\mathcal{G}}}}}\in{{\mathbb{C}}^{{{\omega}_{\mu_{1}}}\times{{\omega}_{\mu_{2}}}}}$.
Since $N_{h}^{{\rm{BS}}}-P+1>L>P$ and $N_{v}^{{\rm{BS}}}-P+1>R>P$,
$r({{{\bm{\mathcal{G}}}_{u}}({n_{s}},{f_{1}})})=P$. Therefore,
${{{\bm{\mathcal{G}}}_{u}}({n_{s}},{f_{1}})}$ and
${{{\bm{\tilde{\mathcal{G}}}}_{u}}({n_{s}},{f_{1}})}$ contain the angular
information. Based on the expression of
${{{\bm{\mathcal{G}}}_{u}}({n_{s}},{f_{1}})}$ in Eq. (60),
$r({{\check{\bf{E}}}_{1}})=r({{\bf{Y}}_{\omega}})=r({\bf{Z}}_{\omega_{\tau}})=r({{\check{\bf{F}}}_{1}})=P$.

Based on Eq. (65), by applying the inequalities

$r(A)\\!+\\!r(B)\\!-\\!{n_{1}}\\!\leq\\!r(AB)\\!\leq\\!\min(r(A),\\!r(B)),\\!\\!\\!\\!$
(127) $r(A+B)\leq r(A,B)\leq r(A)+r(B),\\!\\!\\!\\!$ (128)

with $A\in{{\mathbb{C}}^{{m_{1}}\times{n_{1}}}}$ and
$B\in{{\mathbb{C}}^{{n_{1}}\times{m_{2}}}}$, the rank of
${{\check{\bf{E}}}_{\omega,2}}$ satisfies

$\begin{array}[]{l}r({{\check{\bf{E}}}_{\omega,2}})\geq
r({{\check{\bf{E}}}_{1}}+\cdots+{{\check{\bf{E}}}_{1}}{\bf{Z}}_{\omega_{\tau}}^{Q-1})\\\
\ \ \ \ \ \ \ \ \ \ \geq
r({{\check{\bf{E}}}_{1}})+r((\sum\limits_{q=1}^{Q}{{\bf{Z}}_{\omega_{\tau}}^{q-1}}))-P\\\
\ \ \ \ \ \ \ \ \ \ =P,\end{array}\\!\\!\\!\\!$ (129)

and $r({{\check{\bf{E}}}_{\omega,2}})\leq\min(LRQ,P)=P$, since
${{\check{\bf{E}}}_{\omega,2}}\in{{\mathbb{C}}^{LRQ\times P}}$. Consequently,
$r({{\check{\bf{E}}}_{\omega,2}})=P$. Likewise,
$r({{\check{\bf{F}}}_{\omega,2}})=r({{\bf{Y}}_{\omega}})=P$. Based on Eq.
(127), $r({{\bm{\mathcal{G}}}_{u}}({f_{1}}))$ satisfies:

$\begin{array}[]{l}r({{\bm{\mathcal{G}}}_{u}}({f_{1}}))\\!\leq\\!\min(r({{\check{\bf{E}}}_{\omega,2}}),r({{\bf{Y}}_{\omega}}),r({{\check{\bf{F}}}_{\omega,2}}))\\!=\\!P,\end{array}\\!\\!\\!\\!$
(130) $\ \ \ \ \
r({{\bm{\mathcal{G}}}_{u}}({f_{1}}))\\!\geq\\!r({{\check{\bf{E}}}_{\omega,2}})\\!+\\!r({{\bf{Y}}_{\omega}})\\!+\\!r({{\check{\bf{F}}}_{\omega,2}})\\!-\\!2P\\!\\!=\\!P.\\!\\!\\!\\!$
(131)

Thus, $r({{\bm{\mathcal{G}}}_{u}}({f_{1}}))=P$, which still holds if
${N_{s}}=Q=2$.

Likewise, we may derive that the smallest number of subcarriers satisfies
${N_{f}}=K=2$. The estimations of parameters can be calculated by setting
${N_{s}}=Q=2$ and ${N_{f}}=K=2$. The details are omitted. Following the
similar proof procedure in Appendix -A between Eq. (113) and Eq. (122), we may
prove

$\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}}\to\infty}{{\bf{\hat{\Omega}}}(t+{t_{\tau}})}={{\bf{\Omega}}(t+{t_{\tau}})}.\\!\\!\\!\\!$
(132)

In other words, we may obtain

$\mathop{\lim}\limits_{N_{h}^{{\rm{BS}}},N_{v}^{{\rm{BS}}}\to\infty}\\!\\!\\!\\!\\!\\!\frac{{\left\|{{{{{\bf{\hat{h}}}}_{u}}(t+{t_{\tau}})}\\!-\\!{{{\bf{h}}_{u}}(t+{t_{\tau}})}}\right\|_{2}^{2}}}{{\left\|{{{{\bf{h}}_{u}}(t+{t_{\tau}})}}\right\|_{2}^{2}}}=0.\\!\\!$
(133)

Thus, Theorem 2 is proved. $\square$

###  -C Proof of Proposition 1

During the procedure of Doppler estimation, according to Eq. (67),
$N_{h}^{{\rm{BS}}}$ should satisfy $N_{h}^{{\rm{BS}}}\geq{F_{1}}(L,R)$. Under
the condition $LR={\cal Q}$, we define a Lagrange function as

$F(L,R,\lambda)={F_{1}}(L,R)+\lambda(LR-{\cal Q}).\\!\\!\\!\\!$ (134)

Letting $\frac{{\partial F(L,R,\lambda)}}{{\partial L}}=\frac{{\partial
F(L,R,\lambda)}}{{\partial R}}=\frac{{\partial
F(L,R,\lambda)}}{{\partial\lambda}}=0$, we obtain the two extreme points of
$R$ and $L$: $(R_{1},L_{1})$ and $(R_{2},L_{2})$, which are shown in Eq. (91)
and Eq. (92). If ${\cal Q}={\cal{N}}_{s}$, ${R_{2}}$ and ${L_{2}}$ are
ignored. In addition, the condition of $L$ and $R$ in Eq. (86) covers three
extra extreme points e.g., $({R_{3}}=\frac{{{\cal Q}}}{2},{L_{3}}=2)$,
$({R_{4}}=2,{L_{4}}=\frac{{{\cal Q}}}{2})$, and
$({R_{5}}=N_{v}^{{\rm{BS}}},{L_{5}}=\frac{{{\cal Q}}}{N_{v}^{{\rm{BS}}}})$. By
substituting these extreme points into ${F_{1}}(L,R)$, $N_{t}$ satisfies
$N_{t}\geq{f_{N_{v}^{{\rm{BS}}},{N_{s}},1}}$, where
${f_{N_{v}^{{\rm{BS}}},{N_{s}},1}}$ is shown in Eq. (90).

Up to now, the sub-bound of Doppler estimation is derived. Next, we derive the
sub-bound of angle estimation.

In the procedure of angle estimation, $N_{h}^{{\rm{BS}}}$ and ${N_{s}}$
satisfy

$\begin{array}[]{l}(N_{h}^{{\rm{BS}}}\\!-\\!L\\!+\\!1)({N_{s}}\\!-\\!Q\\!+\\!1)\\!\geq\\!\max\left(\\!{\frac{P}{{N_{v}^{{\rm{BS}}}-R+1}},1}\\!\right).\end{array}\\!\\!\\!\\!$
(135)

Based on the different sizes of $R$, the condition of $R$ in Eq. (88)
($LQ\geq\max\left({\frac{P}{{R-1}},4}\right)$, and $N_{v}^{{\rm{BS}}}\geq
R\geq 2$) is divided into four different conditions, i.e., i)
$N_{v}^{{\rm{BS}}}-P+1\geq R\geq\frac{P}{4}+1$; ii)
$\min(N_{v}^{{\rm{BS}}}-P+1,\frac{P}{4}+1)\geq R\geq 2$; iii)
$N_{v}^{{\rm{BS}}}\geq R\geq\max(N_{v}^{{\rm{BS}}}-P+1,\frac{P}{4}+1,2)$; iv)
$P>R\geq\max(N_{v}^{{\rm{BS}}}-P+1,2)$, which lead to
${f_{N_{v}^{{\rm{BS}}},{N_{s}},2}}$ as in Eq. (93). If
$\frac{{\partial{F_{2}}(R)}}{{\partial R}}=0$, we may obtain $R_{6}$ and
$R_{7}$, shown in Eq. (95) and Eq. (96). Depending on the value of $Q$, i.e.,
$Q<\frac{{{N_{s}}+1}}{2}$, $Q>\frac{{{N_{s}}+1}}{2}$ and
$Q=\frac{{{N_{s}}+1}}{2}$, ${F_{2,\rm{max}}}(R)$ is obtained as in Eq. (94).

Finally, the lower-bound of $N_{t}$ is determined by
$N_{t}\geq\max({f_{N_{v}^{{\rm{BS}}},{N_{s}},1}},{f_{N_{v}^{{\rm{BS}}},{N_{s}},2}})$,
and Proposition 1 is proved. $\square$

## References

  * [1] W. Li, H. Yin, and M. Debbah, “A super-resolution channel prediction approach based on extended matrix pencil method,” in _Proc. IEEE Int. Conf. Commun. (IEEE ICC)_ , Seoul, South Korea, May 2022.
  * [2] T. L. Marzetta, “Noncooperative cellular wireless with unlimited numbers of base station antennas,” _IEEE Trans. Wireless Commun._ , vol. 9, no. 11, pp. 3590–3600, Nov. 2010.
  * [3] F. Boccardi, R. W. Heath, A. Lozano, T. L. Marzetta, and P. Popovski, “Five disruptive technology directions for 5G,” _IEEE Commun. Mag._ , vol. 52, no. 2, pp. 74–80, Feb. 2014.
  * [4] J. Jose, A. Ashikhmin, T. L. Marzetta, and S. Vishwanath, “Pilot contamination and precoding in multi-cell TDD systems,” _IEEE Trans. Wireless Commun._ , vol. 10, no. 8, pp. 2640–2651, Aug. 2011.
  * [5] H. Yin, H. Wang, Y. Liu, and D. Gesbert, “Addressing the curse of mobility in massive MIMO with prony-based angular-delay domain channel predictions,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 12, pp. 2903–2917, Dec. 2020.
  * [6] H. Yin, D. Gesbert, M. Filippou, and Y. Liu, “A coordinated approach to channel estimation in large-scale multiple-antenna systems,” _IEEE J. Sel. Areas Commun._ , vol. 31, no. 2, pp. 264–273, Feb. 2013.
  * [7] R. R. Müller, L. Cottatellucci, and M. Vehkaperä, “Blind pilot decontamination,” _IEEE J. Sel. Topics Signal Process._ , vol. 8, no. 5, pp. 773–786, Oct. 2014.
  * [8] A. Adhikary, J. Nam, J. Y. Ahn, and G. Caire, “Joint spatial division and multiplexing the large-scale array regime,” _IEEE Trans. Inf. Theory_ , vol. 59, no. 10, pp. 6441–6463, Oct. 2013.
  * [9] Z. Jiang, A. F. Molisch, G. Caire, and Z. Niu, “Achievable rates of FDD massive MIMO systems with spatial channel correlation,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 5, pp. 2868–2882, May 2015.
  * [10] K. T. Truong and R. W. Heath, “Effects of channel aging in massive MIMO systems,” _Journal of Communications and Networks_ , vol. 15, no. 4, pp. 338–351, Aug. 2013.
  * [11] J. Zheng, J. Zhang, E. Björnson, and B. Ai, “Impact of channel aging on cell-free massive MIMO over spatially correlated channels,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 10, pp. 6451–6466, Oct. 2021.
  * [12] R. Chopra, C. R. Murthy, H. A. Suraweera, and E. G. Larsson, “Performance analysis of FDD massive MIMO systems under channel aging,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 2, pp. 1094–1108, Feb. 2018\.
  * [13] C. Kong, C. Zhong, A. K. Papazafeiropoulos, M. Matthaiou, and Z. Zhang, “Sum-rate and power scaling of massive MIMO systems with channel aging,” _IEEE Trans. Commun._ , vol. 63, no. 12, pp. 4879–4893, Dec. 2015.
  * [14] J. Zhao, H. Xie, F. Gao, W. Jia, S. Jin, and H. Lin, “Time varying channel tracking with spatial and temporal BEM for massive MIMO systems,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 8, pp. 5653–5666, Aug. 2018.
  * [15] S. Uehashi, Y. Ogawa, T. Nishimura, and T. Ohgane, “Prediction of time-varying multi-user MIMO channels based on DOA estimation using compressed sensing,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 1, pp. 565–577, Jan. 2019.
  * [16] L. Gaudio, M. Kobayashi, G. Caire, and G. Colavolpe, “On the effectiveness of OTFS for joint radar parameter estimation and communication,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 9, pp. 5951–5965, Sept. 2020.
  * [17] T. Zemen, L. Bernado, N. Czink, and A. F. Molisch, “Iterative time-variant channel estimation for 802.11p using generalized discrete prolate spheroidal sequences,” _IEEE Trans. Veh. Technol._ , vol. 61, no. 3, pp. 1222–1233, Mar. 2012.
  * [18] Y. Zhu, X. Dong, and T. Lu, “An adaptive and parameter-free recurrent neural structure for wireless channel prediction,” _IEEE Trans. Commun._ , vol. 67, no. 11, pp. 8086–8096, Nov. 2019.
  * [19] J. Yuan, H. Q. Ngo, and M. Matthaiou, “Machine learning-based channel prediction in massive MIMO with channel aging,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 5, pp. 2960–2973, May 2020.
  * [20] M. Najla, Z. Becvar, P. Mach, and D. Gesbert, “Predicting device-to-device channels from cellular channel measurements: a learning approach,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 11, pp. 7124–7138, Nov. 2020.
  * [21] C. Wu, X. Yi, Y. Zhu, W. Wang, L. You, and X. Gao, “Channel prediction in high-mobility massive MIMO: from spatio-temporal autoregression to deep learning,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 7, pp. 1915–1930, July 2021.
  * [22] Y. Hua, “Estimating two-dimensional frequencies by matrix enhancement and matrix pencil,” _IEEE Trans. Signal Process._ , vol. 40, no. 9, pp. 2267–2280, Sept. 1992.
  * [23] A. Gaber and A. Omar, “A study of wireless indoor positioning based on joint TDOA and DOA estimation using 2-D matrix pencil algorithms and IEEE 802.11ac,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 5, pp. 2440–2454, May 2015.
  * [24] ——, “Joint estimation of time delays and directions of arrival using proper set of antenna elements of a high-order antenna array,” _Digit. Signal Process._ , vol. 94, pp. 114–124, Nov. 2019.
  * [25] R. Schmidt, “Multiple emitter location and signal parameter estimation,” _IEEE Trans. Antennas Propag._ , vol. 34, no. 3, pp. 276–280, Mar. 1986.
  * [26] R. Roy and T. Kailath, “ESPRIT-estimation of signal parameters via rotational invariance techniques,” _IEEE Trans. Acoust., Speech, Signal Process._ , vol. 37, no. 7, pp. 984–995, July 1989.
  * [27] 3GPP, _Study on channel model for frequencies from 0.5 to 100 GHz (Release 16)_. Technical Report TR 38.901, available: http://www.3gpp.org.
  * [28] K. C. Huarng and C. C. Yeh, “A unitary transformation method for angle-of-arrival estimation,” _IEEE Trans. Signal Process._ , vol. 39, no. 4, pp. 975–977, Apr. 1991.
  * [29] H. Yin, L. Cottatellucci, D. Gesbert, R. R. Müller, and G. He, “Robust pilot decontamination based on joint angle and power domain discrimination,” _IEEE Trans. Signal Process._ , vol. 64, no. 11, pp. 2990–3003, June 2016.
  * [30] L. Sun and M. R. McKay, “Eigen-based transceivers for the MIMO broadcast channel with semi-orthogonal user selection,” _IEEE Trans. Signal Process._ , vol. 58, no. 10, pp. 5246–5261, Oct. 2010.
  * [31] Y. Zhu, X. Dong, and T. Lu, “An adaptive and parameter-free recurrent neural structure for wireless channel prediction,” _IEEE Trans. Commun._ , vol. 67, no. 11, pp. 8086–8096, Nov. 2019.