# Low-Complexity Eigenvector Prediction-based Precoding Matrix Prediction in
Massive MIMO with Mobility

Ziao Qin, Haifan Yin, , and Weidong Li Z. Qin, H. Yin and W. Li are with
School of Electronic Information and Communications, Huazhong University of
Science and Technology, 430074 Wuhan, China (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>weidongli@hust.edu.cn).This work was supported by the
National Natural Science Foundation of China under Grant 62071191. The
corresponding author is Haifan Yin.

###### Abstract

In practical massive multiple-input multiple-output (MIMO) systems, the
precoding matrix is often obtained from the eigenvectors of channel matrices
and is challenging to update in time due to finite computation resources at
the base station, especially in mobile scenarios. In order to reduce the
precoding complexity while enhancing the spectral efficiency (SE), a novel
precoding matrix prediction method based on the eigenvector prediction (EGVP)
is proposed. The basic idea is to decompose the periodic uplink channel
eigenvector samples into a linear combination of the channel state information
(CSI) and channel weights. We further prove that the channel weights can be
interpolated by an exponential model corresponding to the Doppler
characteristics of the CSI. A fast matrix pencil prediction (FMPP) method is
also devised to predict the CSI. We also prove that our scheme achieves
asymptotically error-free precoder prediction with a distinct complexity
advantage. Simulation results show that under the perfect non-delayed CSI, the
proposed EGVP method reduces floating point operations by 80% without losing
SE performance compared to the traditional full-time precoding scheme. In more
realistic cases with CSI delays, the proposed EGVP-FMPP scheme has clear SE
performance gains compared to the precoding scheme widely used in current
communication systems.

###### Index Terms:

massive MIMO, CSI, mobility, precoding matrix prediction, channel prediction.

##  I Introduction

Massive multiple-input multiple-output (MIMO) is an integral part of the fifth
generation (5G) new radio (NR) mobile communication systems due to its
outstanding performance in terms of spectral efficiency (SE) and energy
efficiency [1, 2]. One of the key feature of massive MIMO is the large number
of antennas at the base station (BS). A sufficient number of BS antennas
enables linear precoders at the BS, such as zero-forcing (ZF), to achieve
near-optimal SE, as well as error-free channel estimation and perfect noise
cancellation [3]. In order to achieve high spatial multiplexing gains, the
downlink precoders are calculated based on the prior-estimated downlink
channel state information (CSI).

The CSI acquisition dilemma in massive MIMO systems with non-negligible
mobility has recently attracted considerable attention in both academia [4, 5]
and industry [6, 7]. Usually, channel prediction is an effective way to help
the BS to timely update CSI in the fast time-varying massive MIMO systems. In
this scenario, additional methods should be introduced to mitigate the
influence of high mobility on CSI acquisition [5]. Recently, many research
works have proposed effective channel prediction methods. By exploiting the
specific angle-delay-Doppler structure of the multipath, a Prony-based
angular-delay domain (PAD) channel prediction was proposed in [8]. For an
orthogonal time frequency system (OTFS) system, the channel prediction was
achieved by a path division multiple access (PDMA) in [9]. The authors of [10]
proposed a channel prediction method based on the attention mechanism in
machine learning to avoid propagation error. These works leveraged the channel
prediction to update the CSI through a polynomial framework [8] or non-
polynomial algorithm [10].

In addition to the dilemma of timely CSI acquisition, the overwhelming
precoding complexity is one of the major challenges for the practical massive
MIMO systems. In current mobile communication systems, the UEs are equipped
with two or more antennas [11] and eigen zero-forcing (EZF) [12] is the de-
facto method for precoding [13, 14]. The EZF precoding scheme inevitability
introduces an eigenvector decomposition (EVD) or singular value decomposition
(SVD) of the CSI to select a particular eigenvector as the precoder for the
UEs. The large number of BS antennas results in a high dimensional SVD
operation. There are some existing works addressing the precoding complexity
reduction in massive MIMO. The work [15] proposed two low-complexity
eigenvector-based estimators utilizing partitioned sub-array auto-correlation
and cross-correlation combining. The authors in [16] exploited cross-carrier
correlation to interpolate the Gram matrix for precoding and reduce the
complexity. And the authors in [17] reduced the precoding complexity by
interpolating the CSI with a flag-manifold-based method.

In practical 5G systems with mobility, the timely precoding is more
challenging. Due to the limited computational ability, the BS often updates
the precoders in a periodic way [11] and it is hard to update the precoders in
every subframe even though the CSI may be predicted [8, 9, 10]. Therefore, the
major challenge lies in the timely and accurate precoding matrix prediction
with low-complexity.

To track the fast time-varying channel while reducing the precoding
complexity, we utilize several periodic eigenvectors of the CSI to predict the
rest of the precoders, called the eigenvector prediction (EGVP) method. The
key idea of EVGP is to decompose the precoder into a linear combination of the
channel weights and channels. On the one hand, we prove that the channel
weights can be estimated with a complex exponential model, which enables a
polynomial-complexity channel weight interpolation by the periodic eigenvector
samples. On the other hand, the downlink channel prediction is based on an
uplink channel parameter prediction with several uplink channel estimation
samples and the channel reciprocity in the time division duplex (TDD) system.
As a result, the precoding matrix is predicted with polynomial complexity and
the CSI delay can be handled at the same time.

The proposed EGVP-based prediction method is different from the traditional
channel prediction schemes [8, 9, 10] and precoders complexity reduction
schemes [15, 16, 17]. The traditional channel prediction methods like [8, 9,
10] often assumed a full-time precoder updating scheme, which is difficult to
implement in the practical system with limited computational capability. The
studies in [15, 16, 17] reduced the precoding complexity in one way or
another, while failing to consider the fast-varying channel in mobile
environments. To the best of our knowledge, this is the first paper to predict
the precoders through a low-complexity eigenvector prediction method in a
massive MIMO system with mobility. The main contributions of our paper are
summarized below:

  * •

We predict the precoding matrix in a massive MIMO system by a low-complexity
channel eigenvector prediction, called EGVP method. A fast matrix pencil
prediction-based EGVP (EGVP-FMPP) method is further proposed to deal with the
CSI delay. The two schemes both show great advantages in terms of the
precoding complexity and SE, especially in practical systems with limited SVD
capability at the BS.

  * •

We decompose the precoding matrix prediction problem into two sub-problems:
the channel weight interpolation problem and the channel prediction problem.
We further prove that the channel weight can be estimated by a complex
exponential model associated with the Doppler characteristics of the channels.
The relationship between the periodic length of the SVD of the channel matrix
and the maximum speed of the UEs is derived to achieve a non-overlapping
channel eigenvector sampling in EGVP.

  * •

We give asymptotic analyses of the prediction error (PE) performances of the
proposed schemes. The results show that two uplink channel samples and twice
the largest number of non-orthogonal paths among all downlink channels may
lead to an asymptotically error-free channel prediction and eigenvector
interpolation, respectively. Moreover, the numerical results show that EGVP
method still functions well in practical systems with limited number of
antennas and bandwidth.

  * •

We demonstrate the performance advantages of the EGVP method over the
benchmark schemes, such as the SE, PE and complexity of precoding. For a
common system parameter configuration, the EGVP scheme can reduce the
complexity of the traditional full-time SVD-based precoding scheme by 80%
without lossing SE performance. According to the numerical results under an
industrial channel model, the SE improvements of the EGVP-FMPP scheme over the
benchmarks vary from 20.6% to 49.1%.

The rest of the paper is divided into six sections. Sec. II formulates the
downlink channel and the corresponding eigenvector. Sec. III proposes EGVP
algorithm given perfect uplink CSI and Sec. IV introduces EGVP-FMPP method to
combat the CSI delay. Sec. V analyzes the asymptotic performance and the
complexity of our schemes. Sec. VI shows the numerical results of our proposed
methods in different scenarios. Sec. VII concludes this paper.

The bold symbol stands for a matrix or a vector. The matrix dimension is
denoted by ${\mathbb{C}^{a\times b}}$, where $a$ is the number of rows and $b$
is the number of columns. The Kronecker product symbol is $\otimes$. The
Moore-Penrose inversion, conjugate transpose, conjugation, transpose of the
matrix are denoted by the superscripts of ${\dagger}$, $H$, $*$, and $T$,
respectively. The absolute value of a complex number is $\left|x\right|$. The
expectation of a random variable $x$ over the time and the frequency is
$\mathbb{E}\left\\{{{x}}\right\\}$. The inner product of two vectors is
$\left\langle{{\bf{x}},{\bf{y}}}\right\rangle$. The convergence to infinity of
a variable is denoted by $x\to\infty$. The modulo operation is abbreviated as
$\bmod\left({x,y}\right)$.

##  II Formulation of CSI and eigenvector

We consider a TDD massive MIMO system with $N_{t}$ antennas at the BS. The
downlink channel can be estimated through uplink sounding reference signals
(SRS) and the exploitation of channel reciprocity. The $K$ UEs are randomly
distributed in the cell and each UE is equipped with $M$ antennas. The number
of UE antennas is usually limited in practical communication systems, $M\leq
4$ . The $k$-th UE moves at a certain speed $v_{k}$. The wideband system has
$N_{f}$ subcarriers and the center frequency is $f_{0}$. The duration of a
subframe is $\Delta_{t}$ and $\Delta_{f}$ is the subcarrier spacing.

###  II-A Channel model

The channel sparsity in spatial and frequency domain can be exploited by
projecting the wideband channel into the angular-delay domain [8, 11]. As a
result, we consider the wideband channel in this paper. However, the
generalization of our methods to the narrow-band case is straightforward.

Let the collection of the wideband downlink channel at all subcarriers of the
$k$-th UE be

${{\bf{H}}_{k}}\left(t\right)=\left[{\begin{array}[]{*{20}{c}}{{{\bf{h}}_{1,k}}\left(t\right)}&{{{\bf{h}}_{2,k}}\left(t\right)}&\cdots&{{{\bf{h}}_{M,k}}\left(t\right)}\end{array}}\right].$
(1)

In the following, we drop the subscript $k$ for simplicity. Each column
${\bf{h}}_{m}\in\mathbb{C}^{N_{f}N_{t}\times 1}$ represents the downlink
channel of the $m$-th antenna at all subcarriers. Assume a uniform planar
array (UPA) at the BS [18], the channel ${\bf{h}}_{m}$ is modeled as

${{\bf{h}}_{m}}\left(t\right)=\sum\limits_{p=1}^{P}{{\beta_{p,m}}{e^{j{\omega_{p,m}}t}}{{\bf{d}}_{m}}\left({{\theta_{p,m}},{\varphi_{p,m}},{\tau_{p,m}}}\right)},$
(2)

where
${{\bf{d}}_{m}}\left({{\theta_{p,m}},{\varphi_{p,m}},{\tau_{p,m}}}\right)={{\boldsymbol{\alpha}}_{m}}\left({{\theta_{p,m}},{\varphi_{p,m}}}\right)\otimes{{\boldsymbol{\tau}}_{m}}\left({{\tau_{p,m}}}\right)$
is an angle-delay signature of the $p$-th path. The large scale channel gain
is $\beta_{p,m}$. The 2D steering vector
${{\boldsymbol{\alpha}}_{m}}\left({{\theta_{p,m}},{\varphi_{p,m}}}\right)$ is
a Kronecker product of the horizontal steering and vertical steering vector
[19]. ${\theta_{p,m}}$, ${\varphi_{p,m}}$ are the zenith angle and azimuth
angle, respectively. The total number of paths is $P$. The Doppler frequency
is
${\omega_{p,m}}=2\pi{{v\cos{\phi_{p,m}}{f_{0}}}\mathord{\left/{\vphantom{{v\cos{\phi_{p,m}}{f_{0}}}c}}\right.\kern-1.2pt}c}$,
where $c$ is the speed of light and ${\phi_{p,m}}$ is the angle between the
path $p$ and the speed vector of the UE $k$. The delay vector
${{\boldsymbol{\tau}}_{m}}\left({{\tau_{p,m}}}\right)$ characterizes the delay
information ${{\tau_{p,m}}}$ at all subcarriers

${{\boldsymbol{\tau}}_{m}}\left({{\tau_{p,m}}}\right)={e^{2\pi
j{\tau_{p,m}}{f_{0}}}}\left[{\begin{array}[]{*{20}{c}}1&\cdots&{{e^{j2\pi{\tau_{p,m}}\left({{N_{f}}-1}\right){\Delta_{f}}}}}\end{array}}\right].$

According to the virtual channel representation idea [20], the channel (2) can
be projected into the angle-delay domain. Define a transform matrix as the
Kronecker product of a discrete Fourier transform (DFT) matrix
${{\mathbf{W}}_{1}}\left({{N_{f}}}\right)\in{{\mathbb{C}}^{{N_{t}}\times{N_{t}}}}$
and a DFT matrix
${{\mathbf{W}}_{2}}\left({{N_{f}}}\right)\in{{\mathbb{C}}^{{N_{f}}\times{N_{f}}}}$

${\bf{W}}\left({{N_{f}},{N_{t}}}\right){\rm{=}}{{\bf{W}}_{1}}{\left({{N_{f}}}\right)^{H}}\otimes{{\bf{W}}_{2}}\left({{N_{t}}}\right).$
(3)

Then the downlink channel in angle-delay domain ${{\bf{g}}_{m}}\left(t\right)$
is

${{\bf{g}}_{m}}\left(t\right)={\bf{W}}{\left({{N_{t}},{N_{f}}}\right)^{H}}{{\bf{h}}_{m}}\left(t\right).$
(4)

As the matrix ${\bf{W}}{\left({{N_{t}},{N_{f}}}\right)}$ is an unitary matrix,
the wideband channel (2) can be equivalently calculated by

${{\bf{h}}_{m}}\left(t\right)={\bf{W}}\left({{N_{t}},{N_{f}}}\right){{\bf{g}}_{m}}\left(t\right)=\sum\limits_{n=1}^{{N_{t}}{N_{f}}}{{\eta_{m,n}}\left(t\right){{\bf{w}}_{n}}},$
(5)

where ${{\eta_{m,n}}\left(t\right)}$ is the $n$-th element of the vector
${{\bf{g}}_{m}}\left(t\right)$. It represents the amplitude of the
corresponding angle-delay vector ${\bf{w}}_{n}$, which is the $n$-th column
vector of ${\bf{W}}{\left({{N_{t}},{N_{f}}}\right)}$.

###  II-B Periodic eigenvector acquisition

In our paper, the EZF precoding scheme is enabled by the SVD of the wideband
channel matrix. To be more specific, the downlink precoder is selected from
the eigenvectors of ${\bf{H}}\left(t\right)$. Due to the limited SVD
computation capability of most practical communication systems, the BS often
updates the precoders periodically with a cycle length $T_{\rm{svd}}$, which
is longer than the duration of a subframe.

We now discuss the eigenvector acquisition for precoding at the BS. First, the
SVD of ${\bf{H}}\left(t\right)$ is given by

${{\bf{H}}}\left(t\right)={\bf{U}}\left(t\right){\bf{\Sigma}}\left(t\right){\bf{V}}{\left(t\right)^{H}}={{\bf{U}}_{M}}\left(t\right){{\bf{\Sigma}}_{M}}\left(t\right){{\bf{V}}_{M}}{\left(t\right)^{H}},$
(6)

where ${\bf{U}}\left(t\right)\in\mathbb{C}^{N_{f}N_{t}\times N_{f}N_{t}}$ is
the left-singular matrix and ${\bf{V}}\left(t\right)\in\mathbb{C}^{M\times M}$
is the right-singular matrix.
${\bf{\Sigma}}\left(t\right)\in\mathbb{C}^{N_{f}N_{t}\times M}$ is the
singular value matrix of which elements are the singular values sorted in a
descending order. The subscript $M$ means $M$-truncated in columns. And the
$M$-truncated eigenmatrix
${{\bf{U}}_{M}}\left(t\right)\in\mathbb{C}^{N_{f}N_{t}\times M}$ is defined by

${{\bf{U}}_{M}}\left(t\right)=\left[{\begin{array}[]{*{20}{c}}{{{\bf{u}}_{1}}\left(t\right)}&{{{\bf{u}}_{2}}\left(t\right)}&\cdots&{{{\bf{u}}_{M}}\left(t\right)}\end{array}}\right].$
(7)

Note that the eigenmatrix ${{\bf{U}}_{M}}\left(t\right)$ suffers from a random
phase ambiguity due to the SVD nature. Therefore, a phase-calibration progress
of the eigenvectors (7) should be introduced. During a period of time
${t}\in\left[{{t_{{\rm{in}}}},{t_{{\rm{ed}}}}}\right]$, the eigenvectors
${{\bf{u}}_{m}}\left(t\right)$ are phase-calibrated by

$\left[{\begin{array}[]{*{20}{c}}{{{\overline{\bf{u}}}_{m}}\left({{t_{{\rm{in}}}}}\right)}\\\
{{{\overline{\bf{u}}}_{m}}\left({{t_{{\rm{in}}}}+{\Delta_{t}}}\right)}\\\
\cdots\\\
{{{\overline{\bf{u}}}_{m}}\left({{t_{{\rm{ed}}}}}\right)}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{{{\bf{u}}_{m}}\left({{t_{{\rm{in}}}}}\right)}\\\
{{\Delta_{m}}\left({{t_{{\rm{in}}}}+{\Delta_{t}}}\right){{\bf{u}}_{m}}\left({{t_{{\rm{in}}}}+{\Delta_{t}}}\right)}\\\
\cdots\\\
{{\Delta_{m}}\left({{t_{{\rm{ed}}}}}\right){{\bf{u}}_{m}}\left({{t_{{\rm{ed}}}}}\right)}\end{array}}\right],$

where the phase calibration is calculated by
${\Delta_{m}}\left({{t}}\right)=\frac{{{{\bf{u}}_{m}}{{\left({{t}}\right)}^{H}}{{\bf{u}}_{m}}\left({{t_{{\rm{in}}}}}\right)}}{{\left|{{{\bf{u}}_{m}}{{\left({{t}}\right)}^{H}}{{\bf{u}}_{m}}\left({{t_{{\rm{in}}}}}\right)}\right|}}$.
The initial subframe and ending subframe are $t_{\rm{in}}$ and $t_{\rm{ed}}$,
respectively. The current subframe is $t$.

Obviously, the phase-calibrated eigenvector
$\overline{\bf{u}}_{m}\left(t\right)$ is also the eigenvector of
${{\bf{H}}}\left(t\right)$. Learned from basic algebra knowledge, the phase-
calibration does not affect the SE.

Figure 1: Illustration of the precoding matrix prediction in the EGVP scheme
within a period of $T_{\rm{egsp}}$ subframes.

As the result of the periodic SVD at the BS, eigenvectors are sampled every
$T_{\rm{svd}}$ subframes and the rest of them need to be interpolated. In this
case, our EGVP scheme is devised to interpolate the precoders when the SVD
operation is not available. We refers these subframes as SVD gaps in the rest
of the paper. The basic idea of EGVP lies in utilizing a small amount
$N_{\rm{svd}}$ of eigenvector samples to interpolate the rest of eigenvectors
within a period of time $T_{\rm{egsp}}=N_{s}T_{\rm{svd}}$ which is referred to
as the EGVP interval in the following. The illustration of our EGVP within one
$T_{\rm{egsp}}$ is shown in Fig. 1 and the details will be discussed in the
next section.

##  III EGVP with perfect CSI

In this section, the detailed algorithm of the EGVP method is derived
considering the ideal case where the downlink CSI is known. The realistic case
where the CSI suffers from a non-negligible delay will be discussed in the
next section. Given a perfect CSI, the precoder prediction problem is
transformed into a channel weight interpolation problem. First, the calibrated
eigenmatrix $\overline{\bf{U}}_{M}\left(t\right)$ is written as a linear
combination of the channels and channel weights. It is then shown that the
channel weight can be approximated by a exponential model, which enables a
low-complexity channel weight interpolation. In the end, the precoders are
predicted by the interpolated channel weights and the CSI.

###  III-A Linear eigenvector decomposition

The following proposition is first derived.

###### Proposition 1

The eigenvectors $\overline{\bf{u}}_{m}\left(t\right)$ can be decomposed into
a linear combination of the channels ${\bf{h}}_{j}\left(t\right)$ and channel
weights $a_{m,j}\left(t\right)$

$\overline{\bf{u}}_{m}\left(t\right)=\sum\limits_{j=1}^{M}{{a_{m,j}}\left(t\right){{\bf{h}}_{j}}}\left(t\right).$
(8)

_Proof:_ Please refer to Appendix A.

Based on (31) in Appendix A, the channel weight $a_{m,j}\left(t\right)$
actually characterizes the ratio of the inner product of
${\bf{h}}_{j}\left(t\right)$ and $\overline{\bf{u}}_{m}\left(t\right)$ to the
$m$-th eigenvalue $\chi_{m}\left(t\right)$

${a_{m,j}}\left(t\right)=\frac{{\left\langle{\overline{\bf{u}}_{m}\left(t\right),{{\bf{h}}_{j}}\left(t\right)}\right\rangle}}{{{\chi_{m}}\left(t\right)}}.$
(9)

Utilizing Proposition 1, the eigenmatrix is equivalently written

$\overline{\bf{U}}_{M}\left(t\right)=\left[{\begin{array}[]{*{20}{c}}{{{\bf{h}}_{1}}\left({t}\right)}&{{{\bf{h}}_{2}}\left({t}\right)}&\cdots&{{{\bf{h}}_{M}}\left({t}\right)}\end{array}}\right]{\bf{A}}\left(t\right),$
(10)

where the channel weight matrix ${\bf{A}}\left(t\right)$ is

${\bf{A}}\left(t\right)=\left[{\begin{array}[]{*{20}{c}}{{a_{1,1}}\left(t\right)}&{{a_{2,1}}\left(t\right)}&\cdots&{{a_{M,1}}\left(t\right)}\\\
{{a_{1,2}}\left(t\right)}&{{a_{2,2}}\left(t\right)}&\cdots&{{a_{M,2}}\left(t\right)}\\\
\cdots&\cdots&\cdots&\cdots\\\
{{a_{1,M}}\left(t\right)}&{{a_{2,M}}\left(t\right)}&\cdots&{{a_{M,M}}\left(t\right)}\end{array}}\right].$
(11)

Therefore, the eigenvector interpolation problem is equivalent to the channel
weight interpolation problem.

###  III-B Channel weight estimation model

The interpolation of the channel weights is difficult to obtain due to the
non-linearity and non-singularity nature of SVD. However, it might be
asymptotically approximated by a complex exponential model. Unless otherwise
specified, the asymptotic condition refers to the case where the bandwidth
$N_{f}$ and the number of base station antennas $N_{t}$ approaches infinity.
We will propose a theorem to deal with the channel weight estimation problem.
To begin with, we introduce the following lemma to prove that the inner
product of two downlink channels
${s_{i,j}}\left(t\right)=\left\langle{{{\bf{h}}_{j}}\left(t\right),{{\bf{h}}_{i}}\left(t\right)}\right\rangle$
can be approximated by a complex exponential model.

###### Lemma 1

The channel inner product ${s_{i,j}}\left(t\right)$ can be asymptotically
approximated by a complex exponential model

$\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}{{s_{i,j}}\left(t\right)=\sum\limits_{x\in{{\cal{X}}_{i,j}}}{\rho_{i,j}^{\left(x\right)}{e^{j\omega_{i,j}^{\left(x\right)}t}}}},$
(12)

where ${\rho_{i,j}^{\left(x\right)}}$ is the complex amplitude and
${{e^{j\omega_{i,j}^{\left(x\right)}t}}}$ is the corresponding exponential.
The set ${\cal{X}}_{i,j}$ contains the indices of non-zero angle-delay vectors
and the size of it is ${N_{{\rm{no}},i,j}}$.

_Proof:_ Please refer to Appendix B.

According to Appendix B, the physical meaning of $N_{{\rm{no}},i,j}$ is the
number of non-orthogonal multi-paths between the channel
${\bf{h}}_{i}\left(t\right)$ and ${\bf{h}}_{j}\left(t\right)$. The exponential
${{e^{j\omega_{i,j}^{\left(x\right)}t}}}$ reflects the difference of Doppler
frequency shift between the non-orthogonal path $x$ of the channel
${\bf{h}}_{i}\left(t\right)$ and ${\bf{h}}_{j}\left(t\right)$. We introduce a
technical assumption below on the asymptotic value of
${s_{i,j}}\left(t\right)$.

###### Assumption 1

The asymptotic value of ${s_{i,j}}\left(t\right)$ satisfies

$\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}\frac{{{s_{i,i}}\left(t\right)}}{{{N_{t}}{N_{f}}}}=\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}\frac{{\left\|{{{\bf{h}}_{i}}\left(t\right)}\right\|_{2}^{2}}}{{{N_{t}}{N_{f}}}}={\cal{A}}_{i},$
(13)
$\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}\frac{{{s_{j,j}}\left(t\right)}}{{{N_{t}}{N_{f}}}}=\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}\frac{{\left\|{{{\bf{h}}_{j}}\left(t\right)}\right\|_{2}^{2}}}{{{N_{t}}{N_{f}}}}={\cal{B}}_{j},$
(14)
$\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}\frac{{\left|{{s_{i,j}}\left(t\right)}\right|}}{{{N_{t}}{N_{f}}}}=\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}\frac{{\left|{\left\langle{{{\bf{h}}_{i}}\left(t\right),{{\bf{h}}_{j}}\left(t\right)}\right\rangle}\right|}}{{{N_{t}}{N_{f}}}}={\cal{C}}_{i,j}\neq
0.$ (15)

Based on the law of large numbers, the equations (13) and (14) can be
reasonably assumed. The non-zero constants ${\cal{A}}_{i}$ and ${\cal{B}}_{j}$
mean that the normalized channel gain of ${\bf{h}}_{i}\left(t\right)$ and
${\bf{h}}_{j}\left(t\right)$ converge to constant values. The assumption (15)
is the result of the non-trivial correlation between the channel
${{\bf{h}}_{i}}\left(t\right)$ and ${{\bf{h}}_{j}}\left(t\right)$, which is
caused by the multi-antennas at the UE side and the finite number of multi-
paths of the channel. This assumption is valid under the channel model of 5G
NR [18].

In the end, the following theorem is proposed to prove that the channel weight
$a_{m,j}\left(t\right)$ can be asymptotically approximated by a complex
exponential model.

###### Theorem 1

The channel weight $a_{m,j}\left(t\right)$ can be approximated by a linear
combination of complex exponential functions under Assumption 1.

$\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}{a_{m,j}}\left(t\right){\rm{}}=\sum\limits_{l=1}^{{L_{i,j}}}{b_{m,j}^{\left(l\right)}{e^{j\omega_{m,j}^{\left(l\right)}t}}},$
(16)

where ${e^{j\omega_{m,j}^{\left(l\right)}t}}$ is the exponential and
${b_{m,j}^{\left(l\right)}}$ is the corresponding amplitude. The order
$L_{i,j}$ is the number of exponentials.

_Proof:_ Please refer to Appendix C.

We should note that the complex exponential model of the channel weight in
(16) is not unique. The reason lies in the plain fact of non-unique
eigenvector solution of (42). Fortunately, Appendix C proves that this does
not affect the conclusion of Theorem 1. Moreover, it gives a particular
eigenvector solution ${\overline{\bf{a}}_{m}}\left(t\right)$ related to the
channel inner product ${s_{i,j}}\left(t\right)$. In such case, the number of
exponentials $L_{i,j}$ is associated with the number of non-orthogonal multi-
paths $N_{{\rm{no}},i,j}$. The exponential
${e^{j\omega_{m,j}^{\left(l\right)}t}}$ characterizes the difference of
Doppler shifts between the channel ${\bf{h}}_{m}\left(t\right)$ and
${\bf{h}}_{j}\left(t\right)$. The exact value of $L_{i,j}$ can be computed by
a signal detecting technology [21]. However, a unified $L_{i,j}=L_{\rm{svd}}$
is also a practical option, which will be demonstrated in Sec. VI.

Theorem 1 is the basis of EGVP method. Then we discuss how to predict the
precoders by the eigenvector interpolation.

###  III-C Precoder prediction

In order to interpolate the channel weights within SVD gaps, we utilize
$N_{\rm{svd}}$ eigenvector samples to obtain the channel weight estimation
model in (16). The precoders are predicted by the interpolated eigenvectors
and the CSI.

Denote the eigenvector samples within one EGVP interval $T_{\rm{egsp}}$ by

$\overline{\bf{u}}_{m}\left({{t_{{\rm{inl}}}}}\right),\overline{\bf{u}}_{m}\left({{t_{{\rm{inl}}}}+{T_{{\rm{svd}}}}}\right),\cdots,\overline{\bf{u}}_{m}\left({{t_{{\rm{inl}}}}+\left({{N_{\rm{svd}}}-1}\right){T_{{\rm{svd}}}}}\right),$

where $t_{\rm{inl}}$ is the initial subframe. According to the physical
meaning of the channel weight model (16) in Theorem 1, the following corollary
is easily derived to ensure the Nyquist sampling frequency.

###### Corollary 1

The SVD sample length is upper bounded by

${T_{{\rm{svd}}}}\leq\frac{c}{{2{v_{\max}}{f_{0}}}},$ (17)

where the right hand side denotes the inverse of the maximum Doppler frequency
of the channel.

This corollary points out an upper bound of the SVD cycle given a maximum UE
velocity. Considering a TDD system with center frequency $f_{0}=$ 3.5 GHz and
subframe length $\Delta_{t}=$ 0.5 ms, if the maximum speed of UE is about 123
km/h, the SVD cycle should not exceed $T_{\rm{svd}}=1.5$ ms. In this case, the
eigenvector samples update every 3 subframes.

The channel weight estimation based on (16) is a classic exponential parameter
estimating problem. Many traditional algorithms, like Prony [22], maximum
likelihood [23], Matrix Pencil [24], are capable to solve this problem. Due to
the low complexity and robustness to noise, we utilize Matrix Pencil method to
estimate the complex exponential model (16) under a prediction order
constraint $N_{\rm{svd}}\geq 2L_{\rm{svd}}$. The detailed Matrix Pencil
algorithm is omitted in our paper and can be found in [24, 25]. Based on the
estimated complex exponential model (16), the interpolated channel weight
${{\hat{a}}_{m,j}}\left({{t_{p}}}\right)$ at subframe $t_{p}$ is

${{\hat{a}}_{m,j}}\left({{t_{p}}}\right)=\sum\limits_{l=1}^{{L_{{\rm{svd}}}}}{b_{m,j}^{\left(l\right)}{e^{j\frac{{\omega_{m,j}^{\left(l\right)}}}{{{T_{{\rm{svd}}}}}}{t_{p}}}}},$
(18)

where the subframe $t_{p}$ satisfies

$\bmod\left({t_{p}-{t_{{\rm{in}}}},T_{\rm{svd}}}\right)\neq
0,t_{p}\in\left[{{t_{{\rm{in}}}},{t_{{\rm{in}}}}+\left({{N_{\rm{svd}}}-1}\right){T_{{\rm{svd}}}}}\right].$

Combined with the CSI, the eigenvector is reconstructed by

$\hat{\bf{u}}_{m}\left({{t_{p}}}\right)=\sum\limits_{j=1}^{M}{{{\hat{a}}_{m,j}}}\left({{t_{p}}}\right){{\bf{h}}_{j}}\left({{t_{p}}}\right).$
(19)

The whole EGVP algorithm is summarized in Algorithm 1.

Algorithm 1 EGVP algorithm with perfect CSI in one $T_{\rm{egsp}}$

1: Initialize $T_{\rm{svd}},L_{\rm{svd}},T_{\rm{egsp}}$, $t_{\rm{in}}$,
$t_{\rm{ed}}$ and obtain ${\bf{H}}\left(t\right)$;

2: Obtain periodic eigenvector samples based on (6) and calibrate the phase as
$\overline{\bf{u}}_{m}\left(t\right)$;

3: Decompose the eigenvector $\overline{\bf{u}}_{m}\left(t\right)$ and obtain
the channel weights $a_{m,j}\left(t\right)$ with (8);

4: Generate a prediction matrix with $a_{m,j}\left(t\right)$ and estimate the
model in (16) using Matrix Pencil method;

5: for
$t\in\left[{{t_{{\rm{in}}}},{t_{{\rm{ed}}}}}\right]\cap\,\bmod\,\left({t-{t_{{\rm{in}}}},{T_{{\rm{svd}}}}}\right)\neq
0$ do

6: Interpolate the channel weight utilizing (18) and reconstruct the
eigenvector with the interpolated channel weights and the CSI based on (19);

7: Update $t=t+1$;

8: end for

9: Return the eigenvectors $\hat{\bf{u}}_{m}(t)$ and obtain the precoders.

In this section, we demonstrate the specific procedures of EGVP method where
the perfect downlink CSI is available. However, in fast time-varying channel,
the non-neglectable CSI delay problem declines the system performance. The
delayed CSI will be dealt with in EGVG-FMPP method in the next section.

##  IV EGVP-FMPP with delayed CSI

In Sec. Sec. III, the precoders are predicted where perfect and timely CSI is
assumed. The fast time-varying channel in practical systems makes it highly
challenging to obtain the timely CSI for the BS. To be more specific, the
delay of CSI includes three parts: the signal transmission delay
${T_{{\text{trs}}}^{d}}$, the SVD delay ${T_{{\text{svd}}}^{d}}$ and the
interpolation delay of EGVP ${T_{{\text{int}}}^{d}}$. The total CSI delay is
${T^{d}}=T_{{\text{trs}}}^{d}+T_{{\text{int}}}^{d}+T_{{\text{svd}}}^{d}$. The
signal transmission delay is associated with the Doppler effect. Moreover, the
SVD delay cannot be overlooked due to the limited SVD capability of the base
station.

In this section, we consider the more realistic setting with CSI delays. A
channel prediction method, i.e., FMPP method is introduced in EGVP method to
predict the precoders given the delayed CSI. It shows distinct complexity
advantage over the traditional MP method [24]. The core idea of EGVP-FMPP lies
in the CSI prediction which is based on the channel reciprocity in TDD and
uplink channel parameter prediction before applying EGVP procedure. The
predicted channel is

${\hat{\bf{h}}_{m}}\left({{t+T^{d}}}\right)=\sum\limits_{n=1}^{{N_{t}}{N_{f}}}{{\eta_{m,n}}\left({{t+T^{d}}}\right){{\bf{w}}_{n}}},$
(20)

where the complex amplitude ${{\eta_{m,n}}\left(t\right)}$ can be approximated
by a summation of complex exponential functions [25]

${\eta_{m,n}}\left(t\right)=\sum\limits_{l=1}^{{L_{{\rm{ce}}}}}{{\lambda_{m,n,l}}{e^{j{\omega_{m,n,l}}t}}},$
(21)

where ${L_{{\rm{ce}}}}$ is the number of Doppler shifts
${{e^{j{\omega_{m,n,l}}t}}}$. The complex amplitudes ${{\lambda_{m,n,l}}}$
denote the channel gains and random phases. Considering the consistency of
deployment, we still apply the idea of Matrix Pencil method here to acquire
the exponential model (21). However, the proposed FMPP method shows a
complexity advantage compared to the original Matrix Pencil method. The
details of FMPP method are elaborated below.

Within a sampling interval
$t\in\left[{{t_{\rm{1}}},{t_{N_{{\rm{ce}}}}}}\right]$, the total
${N_{{\rm{ce}}}}$ history complex amplitudes ${\eta_{m,n}}\left(t\right)$ are
utilized to construct two Hankel prediction matrices. The prediction order
${L_{{\rm{ce}}}}$ satisfies ${N_{{\rm{ce}}}}\geq 2{L_{{\rm{ce}}}}$. For
representation simplicity, the subscripts $\left({m,n}\right)$ are dropped
below. One Hankel matrix is
${{\boldsymbol{\eta}}_{1}}\in{\mathbb{C}^{\left({{N_{{\rm{ce}}}}-{L_{{\rm{ce}}}}}\right)\times{L_{{\rm{ce}}}}}}$

${{\boldsymbol{\eta}}_{1}}{\rm{=}}\left[{\begin{array}[]{*{20}{c}}{\eta\left({{t_{{L_{{\rm{ce}}}}+1}}}\right)}&{\eta\left({{t_{{L_{{\rm{ce}}}}}}}\right)}&\cdots&{\eta\left({{t_{2}}}\right)}\\\
{\eta\left({{t_{{L_{{\rm{ce}}}}+2}}}\right)}&{\eta\left({{t_{{L_{{\rm{ce}}}}+1}}}\right)}&\cdots&{\eta\left({{t_{3}}}\right)}\\\
\vdots&\vdots&\ddots&\vdots\\\
{\eta\left({{t_{{N_{{\rm{ce}}}}}}}\right)}&{\eta\left({{t_{{N_{{\rm{ce}}}}-1}}}\right)}&\cdots&{\eta\left({{t_{{N_{{\rm{ce}}}}-{L_{{\rm{ce}}}}+1}}}\right)}\end{array}}\right].$

Similarly, the other is
${{\boldsymbol{\eta}}_{0}}\in{\mathbb{C}^{\left({{N_{{\rm{ce}}}}-{L_{{\rm{ce}}}}}\right)\times{L_{{\rm{ce}}}}}}$

${{\boldsymbol{\eta}}_{0}}{\rm{=}}\left[{\begin{array}[]{*{20}{c}}{\eta\left({{t_{{L_{{\rm{ce}}}}}}}\right)}&{\eta\left({{t_{{L_{{\rm{ce}}}}-1}}}\right)}&\cdots&{\eta\left({{t_{1}}}\right)}\\\
{\eta\left({{t_{{L_{{\rm{ce}}}}+1}}}\right)}&{\eta\left({{t_{{L_{{\rm{ce}}}}}}}\right)}&\cdots&{\eta\left({{t_{2}}}\right)}\\\
\vdots&\vdots&\ddots&\vdots\\\
{\eta\left({{t_{{N_{{\rm{ce}}}}-1}}}\right)}&{\eta\left({{t_{{N_{{\rm{ce}}}}-2}}}\right)}&\cdots&{\eta\left({{t_{{N_{{\rm{ce}}}}-{L_{{\rm{ce}}}}}}}\right)}\end{array}}\right].$

The following relationship holds according to [24]

${{\boldsymbol{\eta}}_{0}}={{\bf{Z}}_{1}}{\bf{\Lambda}}{{\bf{Z}}_{2}},{{\boldsymbol{\eta}}_{1}}={{\bf{Z}}_{1}}{\bf{\Lambda}}{{\bf{Z}}_{0}}{{\bf{Z}}_{2}},$
(22)

where the matrices
${{{\bf{Z}}_{0}}}\in{\mathbb{C}^{{L_{{\rm{ce}}}}\times{L_{{\rm{ce}}}}}}$,
${{{\bf{Z}}_{1}}}\in{\mathbb{C}^{\left({{N_{{\rm{ce}}}}-{L_{{\rm{ce}}}}}\right)\times{L_{{\rm{ce}}}}}}$
and ${{{\bf{Z}}_{2}}}\in{\mathbb{C}^{{L_{{\rm{ce}}}}\times{{L_{{\rm{ce}}}}}}}$
are obtained by the exponentials ${{e^{j{\omega_{m,n,l}}t}}}$ according to the
definitions in [24]. The diagonal amplitude matrix is ${\bf{\Lambda}}$ and the
diagonal elements of it are ${{\lambda_{m,n,l}}}$. Based on the definition,
${{\bf{Z}}_{1}}$ is of full column rank and ${{\bf{Z}}_{2}}$ is of full row
rank. Hence, Eq. (22) is equivalently written as

${\bf{\Lambda}}={{\bf{Z}}_{1}}^{\dagger}{{\boldsymbol{\eta}}_{0}}{{\bf{Z}}_{2}}^{\dagger},{{\boldsymbol{\eta}}_{1}}={{\bf{Z}}_{1}}{\bf{\Lambda}}{{\bf{Z}}_{0}}{{\bf{Z}}_{2}}={{\boldsymbol{\eta}}_{0}}{{\bf{Z}}_{2}}^{\dagger}{{\bf{Z}}_{0}}{{\bf{Z}}_{2}}.$
(23)

Therefore, the following relationship can be derived

$\displaystyle{{\boldsymbol{\eta}}_{{T^{d}}+1}}$
$\displaystyle={{\boldsymbol{\eta}}_{{T^{d}}-1}}{{\boldsymbol{Z}}_{2}}^{\dagger}{{\bf{Z}}_{0}}{{\bf{Z}}_{2}}={{\boldsymbol{\eta}}_{{T^{d}}-2}}{{\bf{Z}}_{2}}^{\dagger}{{\bf{Z}}_{0}}^{2}{{\bf{Z}}_{2}}=\cdots$
(24)
$\displaystyle={{\boldsymbol{\eta}}_{0}}{{\bf{Z}}_{2}}^{\dagger}{\left({{{\bf{Z}}_{0}}}\right)^{{T^{d}}}}{{\bf{Z}}_{2}}.$

As a result, the element of the first column and last row of
${{\boldsymbol{\eta}}_{{T^{d}}+1}}$ is the predicted amplitude by $T^{d}$
subframes, i.e., ${\eta_{m,n}}\left({{t_{d}}}\right)$, where
${t_{d}}={t_{1}}+{N_{{\rm{ce}}}}+{T^{d}}$. Hence, the predicted channel
${\hat{\bf{h}}_{m}}\left({{t_{d}}}\right)$ is obtained by (20). Similar to
EGVP method in Sec. III, substitute the channels with
${\hat{\bf{h}}_{m}}\left({{t_{d}}}\right)$ in (19) and the precoders can be
predicted. In this case, the EGVP method with predicted CSI is called EGVP-
FMPP in our paper.

The complexity advantage of the FMPP over the traditional MP method [24] is
distinct. The complexity of the MP is
${\mathcal{O}}\left({{N_{t}}{N_{f}}\left({{N_{{\text{ce}}}}^{2}+{L_{{\text{ce}}}}^{2}+{L_{{\text{ce}}}}\log{T^{d}}}\right)}\right)$
whereas the complexity of the FMPP is
${\mathcal{O}}\left({{N_{t}}{N_{f}}\left({{L_{{\text{ce}}}}^{2}\left({{N_{{\text{ce}}}}-{L_{{\text{ce}}}}}\right)+{L_{{\text{ce}}}}\log{T^{d}}}\right)}\right)$.
Since ${N_{{\text{ce}}}}\geqslant 2{L_{{\text{ce}}}}$, the complexity
reduction from the MP to the FMPP is
${\mathcal{O}}\left({{N_{t}}{N_{f}}\left({{N_{{\text{ce}}}}^{2}+{L_{{\text{ce}}}}^{2}+{L_{{\text{ce}}}}^{3}-{L_{{\text{ce}}}}^{2}{N_{{\text{ce}}}}}\right)}\right)\geqslant{\mathcal{O}}\left({{N_{t}}{N_{f}}{L_{{\text{ce}}}}^{2}\left({5-{L_{{\text{ce}}}}}\right)}\right)$.
Therefore, given $L_{\rm{ce}}\leq 5$, the FMPP method shows a distinct
complexity advantage. Moreover, this condition is often valid due to the
channel sparsity in massive MIMO systems.

The proposed FMPP method is also different from our previous channel
prediction method, joint angel-delay-Doppler (JADD) in FDD [25]. To be more
specific, the downlink channel prediction in JADD scheme depends on the
feedback amplitudes and the estimated exponetials. On the contrary, FMPP
method only relies on the estimated exponetials predicted by the uplink
channel samples, which means that there is no need to compute the amplitudes.

In the next section, the theoretical performances of our two EGVP schemes are
demonstrated.

##  V Performance analysis

In this section, we first analyze how many samples are sufficient to achieve
asymptotically error-free channel prediction and eigenvector prediction. Then
a complexity analysis of our EGVP method is given to demonstrate its advantage
over the state of the art (SoTA) schemes.

###  V-A Prediction error performance

First, the asymptotic channel prediction performance is derived in the
following theorem.

###### Theorem 2

When ${N_{t}},{N_{f}}\to\infty$ and the number of channel samples
${N_{{\rm{ce}}}}=2$, the channel prediction ${\hat{\bf{h}}_{m}}\left(t\right)$
obtained by FMPP is an asymptotically error-free estimation.

$\begin{array}[]{l}\mathop{\lim}\limits_{{N_{t}},{N_{f}}\to\infty}\mathbb{E}\left\\{{\frac{{\left\|{{{\bf{h}}_{m}}\left(t\right)-{{\hat{\bf{h}}}_{m}}\left(t\right)}\right\|_{2}^{2}}}{{\left\|{{{\bf{h}}_{m}}\left(t\right)}\right\|_{2}^{2}}}}\right\\}=0.\end{array}$
(25)

_Proof:_ Please refer to Appendix D.

Theorem 2 shows that two history channel samples are enough to asymptotically
achieve an error-free channel prediction in FMPP method. Furthermore, the PE
is irrelevant with the channel delay $T^{d}$, which means that the effect of
delayed CSI can be compensated.

We should note that this conclusion is based on the assumption that the angle-
delay-Doppler structure of the channel barely varies within $T^{d}$. This
assumption usually holds in moderate mobility [8, 25], e.g., $T^{d}=5$ ms,
$v=60$ km/h.

Then, we aim to analyze the eigenvector PE of EGVP which is measured by the
normalized mean square error (NMSE) between the eigenvectors
${\overline{\bf{u}}_{m}\left(t\right)}$ and the interpolated eigenvectors
${{{\hat{\bf{u}}}_{m}}\left(t\right)}$.

###### Theorem 3

When ${N_{t}},{N_{f}}\to\infty$ and the number of eigenvector samples
$N_{\rm{svd}}\geq 2N_{\rm{no}}^{\rm{max}}$, EGVP method leads to an
asymptotically error-free eigenvector prediction

$\mathop{\lim}\limits_{{N_{t}},{N_{f}}\to\infty}\mathbb{E}\left\\{{\left\|{\frac{{{{\overline{\bf{u}}}_{m}}\left(t\right)-{{\hat{\bf{u}}}_{m}}\left(t\right)}}{{{{\overline{\bf{u}}}_{m}}\left(t\right)}}}\right\|_{2}^{2}}\right\\}=0,$
(26)

where $N_{\rm{no}}^{\rm{max}}$ denotes the maximum number of non-orthogonal
multi-paths among all downlink channels.

_Proof:_ Please refer to Appendix E.

Theorem 3 gives a minimum number of eigenvector samples to achieve an error-
free EGVP. It is associated with the maximum number of non-orthogonal multi-
paths among all downlink channels. In fact, it is exhausting to find the exact
$N_{\rm{no}}^{\rm{max}}$ in practice due to the complicated channel
environments. Even though it is possible to detect the size of
$N_{\rm{no}}^{\rm{max}}$ [26], a fixed sample quantity is more efficient.
Furthermore, a small $N_{\rm{no}}^{\rm{max}}$ means more frequent eigenvector
interpolation and is possible thanks to the channel sparsity of wideband
massive MIMO [27]. Therefore, we tend to configure a small fixed
$N_{\rm{no}}^{\rm{max}}$ and will testify the performance with numerical
results in the next section.

###  V-B Complexity analysis

One of the major motivations of EGVP scheme is to reduce the complexity of the
precoding matrix acquisition at the BS. To demonstrate the complexity
advantage of EGVP scheme, four SoTA schemes are introduced for comparison,
i.e., full-time SVD, periodic SVD, Wiener prediction and approximated gram
matrix interpolation (AGMI) scheme. The downlink CSI is prior-known by the BS
in the complexity analysis.

The full-time SVD scheme refers to the traditional scheme where the precoders
are updated in every subframes like [8, 10]. The complexity order of it is
${\mathcal{O}}\left({{N_{{\text{svd}}}}{T_{{\text{svd}}}}M{{\left({{N_{f}}{N_{t}}}\right)}^{2}}}\right)$.

The periodic SVD scheme means that a periodic precoder updating scheme is
adopted at the BS, similar to the current precoding scheme in 5G NR [11]. AGMI
scheme can be utilized to reduce the complexity of precoding in the frequency
domain [16], but here, the eigenvectors are interpolated in the time domain

${\hat{\mathbf{u}}_{1}}\left({{t_{p}}}\right)={\hat{\mathbf{u}}_{1}}\left({{t_{k}}}\right)\left({\frac{{{t_{k+1}}-{t_{p}}}}{{{T_{{\text{svd}}}}}}}\right)+{\hat{\mathbf{u}}_{1}}\left({{t_{k+1}}}\right)\left({1-\frac{{{t_{k+1}}-{t_{p}}}}{{{T_{{\text{svd}}}}}}}\right),$
(27)

where the interpolated subframe ${t_{p}}$ satisfies
${t_{k}}<{t_{p}}<{t_{k+1}}$. ${t_{k}},{t_{k+1}}$ are the subframes of the
$k$-th and $k+1$-th eigenvector samples, respectively. The complexity order of
periodic SVD and AGMI scheme are the same as
${\mathcal{O}}\left({{N_{{\text{svd}}}}M{{\left({{N_{f}}{N_{t}}}\right)}^{2}}}\right)$.

Wiener prediction method is a classic interpolation method for discrete linear
time-invariant system [28, 29]. It first utilizes the prior-known auto-
correction sequences ${{\bf{r}}_{a}}\left(\tau\right)$ of the signals to
calculate all $L_{w}$ order wiener filter coefficients through the Wiener-Hopf
equation [30]

${\bf{w}}={{\bf{R}}_{a}}{\left(\tau\right)^{\dagger}}{{\bf{r}}_{a}}\left({\tau+{T^{d}}}\right),$
(28)

where ${{\bf{R}}_{a}}\left(\tau\right)$ is the Hermitian Toeplitz matrix of
${{\bf{r}}_{a}}\left(\tau\right)$. The number of interpolated subframes is
$T^{d}$.

In this paper, the precoders are predicted with a linear combination of the
$L_{w}=4$ wiener filter coefficients and recent $L_{w}$ eigenvector samples

${{{u}}_{m}}\left({t+{T^{d}}}\right)=\sum\limits_{l=1}^{{L_{w}}}{w\left(l\right){{{u}}_{m}}\left({t-l}\right)},$
(29)

where $w\left(l\right)$ is the element of ${\bf{w}}$ and
${{u_{m}}\left(t\right)},t\in\left[{0,T-1}\right]$ is the element of the
eigenvector ${{{\bf{u}}_{m}}\left(t\right)}$. $T$ is the period of time. The
auto-correction sequences are obtained by

${{\bf{r}}_{a}}\left(\tau\right)={\mathbb{E}}\left\\{{{u_{m}}\left(t\right){u_{m}}{{\left({t-\tau}\right)}^{*}}}\right\\},\tau\in\left[{-T+1,T-1}\right].$

Therefore, the complexity order of Wiener prediction is
${\mathcal{O}}\left({{N_{{\text{svd}}}}M{{\left({{N_{f}}{N_{t}}}\right)}^{2}}}\right)+{\mathcal{O}}\left({{N_{{\text{svd}}}}^{2}{T_{{\text{svd}}}}^{2}\left({{N_{f}}{N_{t}}}\right)}\right)$.

According to Sec. III, the complexity order of EGVP is
${\mathcal{O}}\left({{N_{{\text{svd}}}}M{{\left({{N_{f}}{N_{t}}}\right)}^{2}}}\right)+{\mathcal{O}}\left({{N_{{\text{svd}}}}{T_{{\text{svd}}}}{M^{2}}\left({{N_{f}}{N_{t}}}\right)}\right)$.
Obviously, the complexity order of the full-time SVD scheme is larger than the
others by
${\mathcal{O}}\left({\left({{T_{{\text{svd}}}}-1}\right){{\left({{N_{f}}{N_{t}}}\right)}^{2}}}\right)$.
The complexity order of EGVP increases
${\mathcal{O}}\left({{N_{{\text{svd}}}}{T_{{\text{svd}}}}{M^{2}}{N_{f}}{N_{t}}}\right)$
compared to the periodic SVD and AGMI scheme. However, compared to the Wiener
prediction scheme, EGVP scheme reduces the complexity order by
${\mathcal{O}}\left({{N_{{\text{svd}}}}{T_{{\text{svd}}}}{N_{f}}{N_{t}}\left({{N_{{\text{svd}}}}{T_{{\text{svd}}}}-M}\right)}\right)$,
which is often positive under the practical configuration $M\leqslant
4,{N_{{\text{svd}}}}{T_{{\text{svd}}}}\geqslant 4$.

Figure 2: Complexity comparison between EGVP and the other schemes in FLOPs
under different configurations of the SVD cycle length and the number of
eigenvector samples.

Furthermore, the complexity within one eigenvector sample cycle
${T_{{\text{egsp}}}}={N_{{\text{svd}}}}{T_{{\text{svd}}}}$ is measured in
floating point operations (FLOPs). Fig. 2 shows the complexity of different
schemes in FLOPs under different configurations of the SVD cycle length and
the number of eigenvector samples. Under a common configuration
$\left({{N_{t}},{N_{f}},{N_{{\text{svd}}}},{T_{{\text{svd}}}}}\right)=\left({64,51,7,7}\right)$,
the complexity of EGVP decreases 79.9% compared to the full-time SVD scheme
and $0.87\%$ compared to the Wiener prediction scheme. The complexity of EGVP
only increases $0.29\%$ compared to the periodic SVD and AGMI scheme.

As a result, the complexity advantage of EGVP scheme is explicit and the
numerical results of the SE performance of EGVP are further shown in the next
section.

##  VI Simulation results

The numerical results of the proposed EGVP and EGVP-FMPP schemes are discussed
in this section. A multi-path channel model widely used in industry, cluster
delay line (CDL), is adopted in our simulations. Based on 3GPP specification
[18], CDL-A channel model is generated with parameters in Table I. In this
model, 23 clusters with 460 multi-paths are aggregated. The delay spread is
300 ns. The center frequency is 3.5 GHz. The bandwidth is 20 MHz and the
subcarrier spacing (SCS) is 30 kHz. Hence, 51 resource blocks are configured
and each subframe lasts 0.5 ms according to [31]. The BS is equipped with a
UPA configured as
$\left({{N_{v}},{N_{h}},{N_{pl}}}\right)=\left({4,8,2}\right)$, where $N_{v}$
and $N_{h}$ are the number of antennas in one column and one row of the
antenna array, respectively. The number of polarizations is denoted by
$N_{pl}$. At the UE side, the receive antennas of all 8 UEs are configured as
$\left({{N_{v}},{N_{h}},{N_{pl}}}\right)=\left({1,2,2}\right)$. We set the
length of SVD cycle as 2.5 ms, i.e., the precoding matrix updates every 5
subframes. The EZF precoding is utilized at the BS and minimum mean mquare
error-interference rejection combining (MMSE-IRC) is applied at the UEs.
Unless other specified, we configure the prediction order of the eigenvector
interpolation $L_{\rm{svd}}$ and channel prediction order $L_{\rm{ce}}$ both
as 3. The number of the channel samples and the eigenvector samples are
${N_{{\text{ce}}}}=2{L_{{\text{ce}}}}+1$ and
${N_{{\text{svd}}}}=2{L_{{\text{svd}}}}+1$, respectively. The system
performance is measured in SE or PE. The PE of the eigenvector is calculated
by (26). And the SE is computed

${R_{{\rm{se}}}}{\rm{=}}\mathbb{E}\left\\{{\sum\limits_{k=1}^{K}{{\rm{log}}\left({1+\frac{{\left\|{{{\bf{h}}_{k}}\left(t,f\right){{\bf{g}}_{k}}\left(t,f\right)}\right\|_{2}^{2}}}{{\sigma_{k}^{2}+\sum\limits_{j\neq
k}^{K}{{{\left|{{{\bf{h}}_{k}}\left(t,f\right){{\bf{g}}_{j}}\left(t,f\right)}\right|}^{2}}}}}}\right)}}\right\\},$

where ${{{\bf{g}}_{k}}\left({t,f}\right)}$ is the precoder of the $k$-th UE
calculated by EZF and $\sigma_{k}^{2}$ is the power of Gaussian noise at the
receiver side. The expectation is taken over subcarriers and subframes. The
SNR at the UE varies from 0 dB to 30 dB.

TABLE I: System Parameters in Simulations Channel model | CDL-A
---|---
Bandwidth $B$ | 20 $\rm{MHz}$
UL/DL carrier frequency $f_{0}$ | 3.5 $\rm{GHz}$
Subcarrier spacing $\Delta_{f}$ | 30 $\rm{kHz}$
Subframe duration ${\Delta_{t}}$ | 0.5 $\rm{ms}$
Resource block $N_{f}$ | 51 RB
Delay spread $\tau_{s}$ | 300 $\rm{ns}$
Number of paths $P$ | 460
| Transmit antenna
---
$\left({{N_{v}},{N_{h}},{N_{pl}}}\right)$
|  $\left({{N_{v}},{N_{h}},{N_{pl}}}\right)=\left({4,8,2}\right)$,
---
the polarization directions are ${0^{\circ}},{90^{\circ}}$
| Receive antenna
---
$\left({{N_{v}},{N_{h}},{N_{pl}}}\right)$
|  $\left({{N_{v}},{N_{h}},{N_{pl}}}\right)=\left({1,2,2}\right)$,
---
the polarization directions are $\pm{45^{\circ}}$
Number of UEs $K$ | 8
SVD cycle $T_{\rm{svd}}$ | 2.5 $\rm{ms}$
CSI delay $T_{{\text{trs}}}^{d}$ | 2.5 $\rm{ms}$
EGVP interpolation delay $T_{{\text{int}}}^{d}$ | 1.5 $\rm{ms}$
SVD delay $T_{{\text{svd}}}^{d}$ | 2.5 $\rm{ms}$
Channel estimation order $L_{\rm{ce}}$ | 3
EGVP order $L_{\rm{svd}}$ | 3

We consider two CSI acquisition scenarios including the timely CSI case and
the delayed CSI case. The CSI delay in the benchmarks is
${T^{d}}=T_{{\text{trs}}}^{d}+T_{{\text{svd}}}^{d}=5$ ms while
${T^{d}}=T_{{\text{trs}}}^{d}+T_{{\text{svd}}}^{d}+T_{{\text{int}}}^{d}=6.5$
ms in EGVP-FMPP scheme considering the channel weight interpolation algorithm.
The performance upper bound throughout our simulation refers to an ideal
scenario where SVD operates at every subframe, denoted by “Full-time SVD”. As
described in Sec. V, the periodic SVD, Wiener prediction and AGMI schemes are
evaluated for performance comparison.

The proposed schemes and benchmarks are evaluated in various scenarios,
including different UE speeds, different CSI delays, different BS antenna
configurations, different SVD cycle lengths and noisy channel sampling cases.
Our EGVP scheme is testified under the perfect CSI assumption and EGVP-FMPP
scheme indicates that a channel prediction method FMPP is applied to combat
the CSI delay $T^{d}$.

In Fig. 3, the SE performances of EGVP and EGVP-FMPP are demonstrated when $v$
= 30 km/h. When perfect CSI is available, our EGVP scheme approaches the upper
bound performance while reducing 79.9% of complexity. When the SNR equals 30
dB, the EGVP increases SE by 5.5% over the AGMI scheme and by 15.8% SE over
the periodic SVD scheme, however introducing only 0.29% additional complexity.
The EGVP scheme reduces 0.87% complexity and enhances 15.5% SE compared to
Wiener scheme. However, when the CSI is delayed, the SE of all schemes
declines. It is noteworthy that our EGVP-FMPP scheme exhibits less degradation
than the benchmarks. The SE performance improvement of EGVP-FMPP scheme varies
from 13.2% to 32.7% compared to the benchmarks. Therefore, we can conclude
that the two EGVP schemes function well in a moderate-mobility system despite
the CSI delay.

Figure 3: SE performances with noise-free channel samples, $v=30$ km/h,
$L_{\rm{ce}}=3$.

Furthermore, the SE performance of the EGVP scheme is shown in Fig. 4 when the
UE speed is $v=60$ km/h. Based on the configuration in Table I and Corollary
1, the maximum UE velocity is $v_{\rm{max}}=61.8$ km/h. When perfect CSI is
assumed, the EGVP scheme approaches 97.8% SE of the upper bound and increases
10.7% over AGMI and 22.3% SE over Wiener and the periodic SVD scheme. In case
of delayed CSI, the orders of the channel estimator for the EGVP-FMPP scheme
are configured to $L_{\rm{ce}}=3$ and $L_{\rm{ce}}=6$. Obviously, the larger
$L_{\rm{ce}}$ results in better SE performance. It can be concluded that a
higher order of channel estimator is more suitable in high mobility scenarios.
In this scenario, the EGVP-FMPP scheme demonstrates a greater SE advantage
over the benchmarks. To be more precise, when $L_{\rm{ce}}=6$, the EGVP-FMPP
scheme enhances the SE performance from 20.6% to 49.1% in comparison to the
benchmarks.

Figure 4: SE performances with noise-free channel samples, $v=60$ km/h.

Fig. 5 demonstrates the PE performances of the two EGVP scheme given different
number of the BS antennas ${N_{t}}\in\left\\{{4,8,16,64,128,256}\right\\}$ and
different speeds of the UEs. In the EGVP-FMPP scheme, the PE decreases as
$N_{t}$ increases. However, the PE of the EGVP scheme remains relatively
stable with respect to $N_{t}$. Furthermore, the results demonstrate that both
EGVP schemes are applicable for a relatively small $N_{t}$ which is more
applicable in practical systems.

Figure 5: PE performances vs the number of antennas with noise-free channel
samples, $L_{\rm{ce}}=3$.

Then, we evaluate the PE performances given different SVD cycle lengths
$T_{\rm{svd}}$ in Fig. 6. The PE performance of all schemes decreases slowly
with $T_{\rm{svd}}$. Both EGVP schemes show distinct PE advantage over all
benchmarks. No doubt a longer SVD cycle length leads to lower computation
complexity yet larger PE. Therefore, the trade-off between the SVD cycle
length and computation complexity can be tricky and the equation (17) is
helpful in this regard.

Figure 6: PE performances of the eigenvectors vs different lengths of the SVD
cycle with noise-free channel samples, $v=30$ km/h, $L_{\rm{ce}}=3$. Figure 7:
PE performances of the eigenvectors vs different CSI delays with noise-free
channel samples, $v=30$ km/h.

Fig. 7 illustrates the PE performances of EGVP-FMPP scheme and the benchmarks
under different CSI delays. The benchmarks consider only the transmission
delay $T_{{\text{trs}}}^{d}$ and the SVD delay $T_{{\text{svd}}}^{d}$ whereas
EGVP-FMPP considers the additional interpolation delay $T_{{\text{int}}}^{d}$.
It can be observed that as the CSI delay increases, the PE performance of all
schemes deteriorates. However, EGVP-FMPP demonstrates superior PE performance
compared to the benchmarks.

Figure 8: SE performances with noisy channel samples, $v$ = 30 km/h, $\gamma$
= 20 dB, 40 dB, $L_{\rm{ce}}=3$.

In the end, we should note that the above results are based on the assumption
of the noise-free channel samples. In reality, noisy sampling is inevitable
and no doubt deteriorates the performances of our schemes. Many classic noise-
canceling methods can be applied in this scenario, such as Prony [32], minimum
decription length (MDL) [33], total linear squares (TLS) [34] and Wiener
filtering [35]. The detailed discussion is beyond the scope of our paper and
the MDL method [33] is applied here to detect signals and suppress noise.
Given two different channel sampling SNR $\gamma$, the SE performances of our
schemes are illustrated in Fig. 8. It can be concluded that the channel
sampling noise has limited impact on the EGVP scheme with perfect CSI.
However, it deteriorates the SE of the EGVP-FMPP scheme, especially when the
noise power is large. Fortunately, the channel sampling SNR $\gamma$ = 40 dB
is sufficient for the EGVP-FMPP scheme to outperform the benchmarks with
noise-free sampling.

##  VII Conclusion

This paper proposed a precoding matrix prediction method called EGVP in
massive MIMO systems with mobility. Based on several periodic eigenvector
samples of the downlink channels, the precoders at the BS were predicted by a
polynomial-complexity channel weight interpolation and channel prediction. It
was enabled by the decomposition of the precoders into a linear combinations
of the CSI and the corresponding channel weights. The channel weights were
interpolated by a complex exponential model. Considering the CSI delay, we
further proposed a EGVP-FMPP method with an additional low-complexity channel
prediction procedure. The asymptotic analyses proved that an error-free
prediction can be achieved in our proposed EGVP schemes. The complexity
advantage of our scheme over the traditional channel prediction schemes and
the other precoding matrix prediction benchmarks were analyzed. Based on the
industrial channel model, our schemes showed substantial improvement in terms
of SE with low complexity in mobile environments.

##  Appendix A Proof of Proposition 1

The eigenvectors ${\overline{\bf{u}}_{m}}\left(t\right)$ obtained from the SVD
of the channel matrix ${\bf{H}}\left(t\right)$ equal to the eigenvectors of
${\bf{H}}\left(t\right){\bf{H}}{\left(t\right)^{T}}$

$\left({\sum\limits_{j=1}^{M}{{{\bf{h}}_{j}}\left(t\right){{\bf{h}}_{j}}{{\left(t\right)}^{H}}}}\right)\cdot{\overline{\bf{u}}_{m}}\left(t\right)={\lambda_{m}}\left(t\right){\overline{\bf{u}}_{m}}\left(t\right).$
(30)

The left hand side of the above equation can be written as

$\left({\sum\limits_{j=1}^{M}{{{\bf{h}}_{j}}\left(t\right){{\bf{h}}_{j}}{{\left(t\right)}^{H}}}}\right){\overline{\bf{u}}_{m}}\left(t\right)=\sum\limits_{j=1}^{M}{{{\bf{h}}_{j}}\left(t\right)\sum\limits_{n=1}^{{N_{t}}{N_{f}}}{{h_{n,j}}\left(t\right){{\bar{u}}_{n,m}}\left(t\right)}},$

where ${{h_{n,j}}\left(t\right)}$ and ${{\overline{u}}_{n,m}}\left(t\right)$
are the $n$-th element of ${{{\bf{h}}_{j}}\left(t\right)}$ and
${\overline{\bf{u}}_{m}}\left(t\right)$, respectively. Denote the channel
weight as

${a_{m,j}}\left(t\right)=\frac{1}{{{\chi_{m}}\left(t\right)}}\sum\limits_{n=1}^{{N_{t}N_{f}}}{{h_{n,j}}\left(t\right){{\overline{u}}_{n,m}}\left(t\right)}.$
(31)

Substitute the weight (31) into the equation (30) and obtain the following
equation

${\overline{\bf{u}}_{m}}\left(t\right)=\sum\limits_{j=1}^{M}{{a_{m,j}}\left(t\right){{\bf{h}}_{j}}}\left(t\right).$
(32)

Thus, Proposition 1 is proved.

##  Appendix B Proof for Lemma 1

Based on the channel sparsity in the angle-delay domain, the $n$-th amplitude
of the vector ${{\bf{g}}_{m}}\left(t\right)$, i.e.,
${{\eta_{m,n}}\left(t\right)}$ can be approximately fitted by a complex
exponential model as mentioned in [8, 25]

${\eta_{m,n}}\left(t\right)\approx\sum\limits_{q=1}^{Q_{m,n}}{{c_{m,n,q}}{e^{j{\omega_{m,n,q}}t}}},$
(33)

where $Q_{m,n}$ is the number of exponentials and ${c_{m,n,q}}$ is the complex
amplitude corresponding to the exponential ${e^{j{\omega_{m,n,q}}t}}$. The DFT
leakage problem of the projection matrix
${\bf{W}}\left({{N_{f}},{N_{t}}}\right)$ limits the approximation accuracy.
Each angle-delay vector maps a corresponding angle-delay signature
${\bf{d}}_{m}$

$\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}{{\bf{w}}_{n}}=\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}\frac{{{{\bf{d}}_{m}}\left({{\theta_{p}},{\varphi_{p}},{\tau_{p}}}\right)}}{{\sqrt{{N_{t}}{N_{f}}}}},\forall
p\in{{\cal P}},n\in{{\cal N}},$ (34)

where the set $\cal{P}$ and $\cal{N}$ are the multi-path index set and the
angle-delay vector index set, respectively. Moreover, in asymptotic condition,
each angular frequency ${{\omega_{m,n}}}$ corresponds to the Doppler frequency
${{\omega_{m,p}}}$ of the $p$-th path, i.e, $Q_{m,n}=1$ [25]. This means that
the complex amplitude $c_{m,n}$ maps the large scale channel parameter
$\beta_{m,p}$ of the $p$-th path. Therefore, the asymptotic approximation
model of (33) is

${\eta_{m,n}}\left(t\right)={c_{m,n}}{e^{j{\omega_{m,n}}t}}\buildrel\Delta\over{=}\frac{{{\beta_{m,p}}{e^{j{\omega_{m,p}}t}}}}{{\sqrt{{N_{t}}{N_{f}}}}}.$
(35)

Based on the channel in angle-delay domain equation (5), the channel inner
product ${s_{i,j}}\left(t\right)$ satisfies

$\begin{array}[]{l}{s_{i,j}}\left(t\right)=\left\langle{{{\bf{h}}_{i}}\left(t\right),{{\bf{h}}_{j}}\left(t\right)}\right\rangle\\\
=\left({\sum\limits_{n=1}^{{N_{t}}{N_{f}}}{{\eta_{j,n}}{{\left(t\right)}^{H}}{{\bf{w}}_{n}}^{H}}}\right)\left({\sum\limits_{n=1}^{{N_{t}}{N_{f}}}{{\eta_{j,n}}\left(t\right){{\bf{w}}_{n}}}}\right)\\\
=\sum\limits_{n=1}^{{N_{t}}{N_{f}}}{{\eta_{j,n}}{{\left(t\right)}^{H}}{\eta_{i,n}}\left(t\right)}\end{array}.$
(36)

Substitute ${\eta_{m,n}}\left(t\right)$ with the equation (35), then the inner
product ${s_{i,j}}\left(t\right)$ is

${s_{i,j}}\left(t\right)=\sum\limits_{n=1}^{{N_{f}}{N_{t}}}{{c_{j,n}}^{H}{e^{-j{\omega_{j,n}}t}}{c_{i,n}}{e^{j{\omega_{i,n}}t}}=\sum\limits_{x\in{{\cal{X}}_{i,j}}}{\rho_{i,j}^{\left(x\right)}{e^{j\omega_{i,j}^{\left(x\right)}t}}},}$

where the set ${{\cal{X}}_{i,j}}$ denotes the indices collection of non-zero
amplitudes. The size of ${{\cal{X}}_{i,j}}$ is denoted by $\left|{{\cal
X}_{i,j}}\right|={N_{{\rm{no}},i,j}}\leq{N_{f}}{N_{t}}$. The non-zero complex
amplitude is ${\rho_{i,j}^{\left(x\right)}}$ and the corresponding exponential
is ${{e^{j\omega_{i,j}^{\left(x\right)}t}}}$.

Therefore, Lemma 1 is proved.

##  Appendix C Proof of Theorem 1

Based on the channel weight definition (31), we rewrite it

${a_{m,j}}\left(t\right)=\frac{{\left\langle{{\overline{\bf{u}}_{m}}\left(t\right),{{\bf{h}}_{j}}\left(t\right)}\right\rangle}}{{{\chi_{m}}\left(t\right)}}=\frac{{{{\bf{h}}_{j}}{{\left(t\right)}^{H}}\sum\limits_{i=1}^{M}{{a_{m,i}}\left(t\right){{\bf{h}}_{i}}\left(t\right)}}}{{{\chi_{m}}\left(t\right)}}.$
(37)

Substitute the channel ${{{\bf{h}}_{j}}\left(t\right)}$ in (37) with (5)

$\begin{array}[]{l}{a_{m,j}}\left(t\right)=\frac{{\sum\limits_{i=1}^{M}{\left({\sum\limits_{n=1}^{{N_{t}}{N_{f}}}{{\eta_{j,n}}{{\left(t\right)}^{H}}{\eta_{i,n}}\left(t\right)}}\right){a_{m,i}}\left(t\right)}}}{{{\chi_{m}}\left(t\right)}}\\\
=\frac{{\sum\limits_{i=1}^{M}{{s_{i,j}}\left(t\right){a_{m,i}}\left(t\right)}}}{{{\chi_{m}}\left(t\right)}}\end{array}.$
(38)

Then, rewrite the above equation in a matrix form

${\bf{S}}\left(t\right)\left[{\begin{array}[]{*{20}{c}}{{a_{1,m}}\left(t\right)}\\\
{{a_{2,m}}\left(t\right)}\\\ \cdots\\\
{{a_{M,m}}\left(t\right)}\end{array}}\right]={\chi_{m}}\left(t\right)\left[{\begin{array}[]{*{20}{c}}{{a_{1,m}}\left(t\right)}\\\
{{a_{2,m}}\left(t\right)}\\\ \cdots\\\
{{a_{M,m}}\left(t\right)}\end{array}}\right],$ (39)

where ${\bf{S}}\left(t\right)$ is the collection of all channel inner products

${\bf{S}}\left(t\right)=\left[{\begin{array}[]{*{20}{c}}{{s_{1,1}}\left(t\right)}&{{s_{1,2}}\left(t\right)}&\cdots&{{s_{1,M}}\left(t\right)}\\\
{{s_{2,1}}\left(t\right)}&{{s_{2,2}}\left(t\right)}&\cdots&{{s_{2,M}}\left(t\right)}\\\
\cdots&\cdots&\cdots&\cdots\\\
{{s_{M,1}}\left(t\right)}&{{s_{M,2}}\left(t\right)}&\cdots&{{s_{M,M}}\left(t\right)}\end{array}}\right],$
(40)

where the element ${s_{i,j}}\left(t\right)$ is defined in (12). Obviously, the
equation (39) is the $m$-th eigen equation of matrix ${\bf{S}}\left(t\right)$.
It shares the same eigenvalue $\chi_{m}\left(t\right)$ with the matrix
${\bf{H}}\left(t\right)$. At the same time, the corresponding eigenvector is
the collection of the channel weights

${{\bf{a}}_{m}}\left(t\right)={\left[{\begin{array}[]{*{20}{c}}{{a_{1,m}}\left(t\right)}&{{a_{2,m}}\left(t\right)}&\cdots&{{a_{M,m}}\left(t\right)}\end{array}}\right]^{T}}.$
(41)

Due to the definition of ${s_{i,j}}\left(t\right)$ in (12), the matrix
${\bf{S}}\left(t\right)$ is a Hermitian matrix. And the eigenvalue
$\chi_{m}\left(t\right)$ is non-zero. Therefore, the matrix
${\bf{S}}\left(t\right)$ is definite and non-singular.

The eigenvalue $\chi_{m}\left(t\right)$ is obtained by solving

$\left({{\bf{S}}\left(t\right)-{\chi_{m}}\left(t\right){{\bf{I}}_{M}}}\right){{\bf{a}}_{m}}\left(t\right)=\bf{0},$
(42)

where ${{\bf{I}}_{M}}$ is a $M$-dimension unit matrix.

We first let $M=2$ and focus on obtaining the first eigenvector
${{\bf{a}}_{1}}\left(t\right)$. The other channel weight solution
${{\bf{a}}_{2}}\left(t\right)$ can be similarly proved. Rewrite the equation
(42) as

$\left[{\begin{array}[]{*{20}{c}}{{s_{1,1}}\left(t\right)-{\chi_{1}}\left(t\right)}&{{s_{1,2}}\left(t\right)}\\\
{{s_{2,1}}\left(t\right)}&{{s_{2,2}}\left(t\right)-{\chi_{1}}\left(t\right)}\end{array}}\right]{{\bf{a}}_{m}}\left(t\right)=\bf{0}.$
(43)

Using Gaussian elimination, turn the the left hand side matrix in (43) into an
upper triangular matrix $\widetilde{\bf{S}}\left(t\right)$

$\left[{\begin{array}[]{*{20}{c}}{{{{S}}_{1,1}}\left(t\right)}&{{{{S}}_{1,2}}\left(t\right)}\\\
0&{{{{S}}_{2,2}}\left(t\right)}\end{array}}\right]{{\bf{a}}_{m}}\left(t\right)={\bf{0}},$
(44)

where each element of $\widetilde{\bf{S}}\left(t\right)$ is denoted by

$\left\\{\begin{array}[]{l}{{{S}}_{1,1}}\left(t\right)={s_{1,1}}\left(t\right)-{\chi_{1}}\left(t\right),\\\
{{{S}}_{1,2}}\left(t\right)={s_{1,2}}\left(t\right),\\\
{{{S}}_{2,2}}\left(t\right)={s_{2,2}}\left(t\right)-{\chi_{1}}\left(t\right)-\frac{{{{\left|{{{{S}}_{1,2}}\left(t\right)}\right|}^{2}}}}{{{{{S}}_{1,1}}\left(t\right)}}.\end{array}\right.$
(45)

As ${\bf{S}}\left(t\right)$ is non-singular, the rank of the matrix
$\widetilde{\bf{S}}\left(t\right)$ is 1. Therefore, the eigen equation (44)
satisfies

$\left\\{\begin{array}[]{l}{a_{1,1}}\left(t\right){{{S}}_{1,1}}\left(t\right)+{a_{2,1}}\left(t\right){{{S}}_{1,2}}\left(t\right)=0,\\\
{{{S}}_{2,2}}\left(t\right)=0,{{{S}}_{1,2}}\left(t\right)\neq
0,\end{array}\right.$ (46)

which equals to

${a_{1,1}}\left(t\right)=-\frac{{{s_{1,2}}\left(t\right)}}{{{s_{1,1}}\left(t\right)-{\chi_{1}}\left(t\right)}}{a_{2,1}}\left(t\right),$
(47)
$\left({{s_{2,2}}\left(t\right)-{\chi_{1}}\left(t\right)}\right)\left({{s_{1,1}}\left(t\right)-{\chi_{1}}\left(t\right)}\right)={\left|{{s_{1,2}}\left(t\right)}\right|^{2}}.$
(48)

Solve the equation (48) and obatin the eigenvalue

${\chi_{1}}\left(t\right){\rm{=}}\frac{{{s_{1,1}}\left(t\right)+{s_{2,2}}\left(t\right)\pm\sqrt{{{\left({{s_{1,1}}\left(t\right)-{s_{2,2}}\left(t\right)}\right)}^{2}}+4{{\left|{{s_{1,2}}\left(t\right)}\right|}^{2}}}}}{2}.$
(49)

As the first eigenvalue $\chi_{1}\left(t\right)$ is the greater one, we take
the greater value of (49). Then the equation (47) becomes

$\left\\{\begin{array}[]{l}{a_{1,1}}\left(t\right)={{{\cal
K}}_{1,2}}\left(t\right){a_{2,1}}\left(t\right),\\\ {{{\cal
K}}_{1,2}}\left(t\right)=\frac{{2{s_{1,2}}\left(t\right)}}{{\sqrt{{{\left({{s_{1,1}}\left(t\right)-{s_{2,2}}\left(t\right)}\right)}^{2}}+4{{\left|{{s_{1,2}}\left(t\right)}\right|}^{2}}}+\left({{s_{2,2}}\left(t\right)-{s_{1,1}}\left(t\right)}\right)}},\end{array}\right.$
(50)

where ${\cal K}_{1,2}\left(t\right)$ is the ratio of ${a_{1,1}}\left(t\right)$
to ${a_{2,1}}\left(t\right)$.

Then, we analyze the asymptotic value of ${\cal K}_{1,2}\left(t\right)$. Let
the normalized denominator of ${\cal K}_{1,2}\left(t\right)$ be

$\begin{array}[]{l}{D_{1,2}}\left(t\right)=\sqrt{\frac{{{{\left({{s_{1,1}}\left(t\right)-{s_{2,2}}\left(t\right)}\right)}^{2}}+4{{\left|{{s_{1,2}}\left(t\right)}\right|}^{2}}}}{{{N_{t}}^{2}{N_{f}}^{2}}}}\\\
+\frac{{\left({{s_{2,2}}\left(t\right)-{s_{1,1}}\left(t\right)}\right)}}{{{N_{t}}{N_{f}}}}.\end{array}$
(51)

Based on Assumption 1, we have

$\begin{array}[]{l}\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}{{\cal{K}}_{1,2}}\left(t\right)=\frac{{\frac{{{s_{1,2}}\left(t\right)}}{{{N_{t}}{N_{f}}}}}}{{{D_{1,2}}\left(t\right)}}\\\
=\frac{{\frac{{{s_{1,2}}\left(t\right)}}{{{N_{t}}{N_{f}}}}}}{{\sqrt{{{\left({{\cal{A}}-{\cal{B}}}\right)}^{2}}+{{\cal{C}}^{2}}}+{\cal{A}}-{\cal{B}}}}\\\
\end{array}.$ (52)

And the absolute value of ${\cal K}_{1,2}\left(t\right)$ converges to

$\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}\left|{{{\cal
K}}_{1,2}\left(t\right)}\right|=\frac{{\cal{C}}}{{\sqrt{{{\left({{\cal{A}}-{\cal{B}}}\right)}^{2}}+{{\cal{C}}^{2}}}+{\cal{A}}-{\cal{B}}}}.$
(53)

Learned from (52) and (53), ${\cal K}_{1,2}\left(t\right)$ can be modeled by
${s_{1,2}}\left(t\right)$ multiplied by a time-invariant scaling factor.

However, the eigenvector ${{\bf{a}}_{1}}\left(t\right)$ is non-unique. Without
losing generality, we first consider a particular solution

$\overline{a}_{2,1}\left(t\right)={{\cal U}}{s_{1,2}}\left(t\right),$ (54)

where $\cal U$ is a time-invariant scaling factor. Learned from (50) and Lemma
1, the particular solution $\overline{a}_{1,1}\left(t\right)$ is

$\begin{array}[]{l}\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}{\overline{a}_{1,1}}\left(t\right)={{\cal{K}}_{1,2}}\left(t\right){\cal{U}}{s_{1,2}}\left(t\right)\\\
=\widetilde{{\cal
U}}{s_{1,2}}{\left(t\right)^{2}}\buildrel\Delta\over{=}\sum\limits_{l=1}^{L_{i,j}}{b_{1,1}^{\left(l\right)}{e^{j\omega_{1,1}^{\left(l\right)}t}}}\end{array},$
(55)

where, $\widetilde{{\cal U}}$ is a time-invariant scaling factor. Thus, the
particular solution of channel weight ${{\bf{\bar{a}}}_{1,1}}\left(t\right)$
can be modeled by a complex exponential model associated with
${s_{1,2}}\left(t\right)$.

Likewise, the general solution of the channel weight $a_{2,1}\left(t\right)$
is also estimated by a complex exponential model

$\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}{a_{2,1}}\left(t\right)\buildrel\Delta\over{=}\sum\limits_{l=1}^{L_{i,j}}{b_{2,1}^{\left(l\right)}{e^{j\omega_{2,1}^{\left(l\right)}t}}}.$
(56)

According to (50), both $a_{2,1}\left(t\right)$ and $a_{1,1}\left(t\right)$
can also be modeled by a complex exponential model which may be irrelevant
with ${s_{1,2}}\left(t\right)$.

Even though the above derivation is made when $M=2$, this limitation can be
relaxed to $M>2$. In this case, an iterative form of channel weight is also
proven to be asymptotically estimated by a complex exponential model.

When $M>2$, perform an elementary row operation and turn the matrix
${{\bf{S}}\left(t\right)-{\chi_{1}}\left(t\right){{\bf{I}}_{M}}}$ to an upper
triangular matrix $\widetilde{\bf{S}}\left(t\right)$. First, define an
iterative matrix
$\widetilde{\bf{S}}^{I}\left(t\right),I\in\left\\{0,1,2,\cdots,M-1\right\\}$
as the result of the $I$-th elementary row operation of
${{\bf{S}}\left(t\right)-{\chi_{1}}\left(t\right){{\bf{I}}_{M}}}$. And the
initial matrix is
$\widetilde{\bf{S}}^{0}\left(t\right)={{\bf{S}}\left(t\right)-{\chi_{1}}\left(t\right){{\bf{I}}_{M}}}$.
Obviously, the last iterative matrix $\widetilde{\bf{S}}^{M-1}\left(t\right)$
equals to the matrix $\widetilde{\bf{S}}\left(t\right)$

$\left[{\begin{array}[]{*{20}{c}}{S_{1,1}^{M{\rm{-}}1}\left(t\right)}&{S_{1,2}^{M{\rm{-}}1}\left(t\right)}&\cdots&{S_{1,M}^{M{\rm{-}}1}\left(t\right)}\\\
0&{S_{2,2}^{M{\rm{-}}1}\left(t\right)}&\cdots&{S_{2,M}^{M{\rm{-}}1}\left(t\right)}\\\
\cdots&\cdots&\cdots&\cdots\\\
0&\cdots&{S_{M{\rm{-}}1,M{\rm{-}}1}^{M{\rm{-}}1}\left(t\right)}&{S_{M{\rm{-}}1,M}^{M{\rm{-}}1}\left(t\right)}\\\
0&0&\cdots&{S_{M,M}^{M{\rm{-}}1}\left(t\right)}\end{array}}\right].$ (57)

Then, we aim to derive the iterative form of the elements of (57) instead of
the closed-form like (50). The elements of the initial matrix
$S_{i,j}^{0}\left(t\right)$ satisfy

$S_{i,j}^{0}\left(t\right)=\left\\{\begin{array}[]{l}{s_{i,j}}\left(t\right),i\neq
j,\\\ {s_{i,j}}\left(t\right)-{\chi_{1}}\left(t\right),i=j.\end{array}\right.$
(58)

When $I=1$, $\widetilde{\bf{S}}^{1}\left(t\right)$ becomes

$\left[{\begin{array}[]{*{20}{c}}{S_{1,1}^{1}\left(t\right)}&{S_{1,2}^{1}\left(t\right)}&\cdots&{S_{1,M}^{1}\left(t\right)}\\\
0&{S_{2,2}^{1}\left(t\right)}&\cdots&{S_{2,M}^{1}\left(t\right)}\\\
\cdots&\cdots&\cdots&\cdots\\\
0&\cdots&{S_{M{\rm{-}}1,M{\rm{-}}1}^{1}\left(t\right)}&{S_{M{\rm{-}}1,M}^{1}\left(t\right)}\\\
0&{S_{M,2}^{1}\left(t\right)}&\cdots&{{S_{M,M}}\left(t\right)}\end{array}}\right],$
(59)

where the elements satisfy

$S_{i,j}^{1}\left(t\right)=\left\\{\begin{array}[]{l}S_{i,j}^{0}\left(t\right),i\leq
1,\\\
S_{i,j}^{0}\left(t\right)-\frac{{S_{i,1}^{0}\left(t\right)}}{{S_{1,1}^{0}\left(t\right)}}S_{1,j}^{0}\left(t\right),1<i\leq
j,\\\ 0,i>j.\end{array}\right.$ (60)

The indices $i$ and $j$ denote the $i$-th row and $j$-column of
$\widetilde{\bf{S}}\left(t\right)$, respectively. When $I\geq 1$, the general
form of the elements of $\widetilde{\bf{S}}^{I}\left(t\right)$ are

$S_{i,j}^{I}\left(t\right)=\left\\{\begin{array}[]{l}S_{i,j}^{I-1}\left(t\right),i\leq
I,\\\
S_{i,j}^{I-1}\left(t\right)-\frac{{S_{i,I}^{I-1}\left(t\right)}}{{S_{I,I}^{I-1}\left(t\right)}}S_{I,j}^{I-1}\left(t\right),i>I,j>I,\\\
0,j\leq I<i.\end{array}\right.$ (61)

Based on the iterative equation (61) and initial condition (58), the eigen
equation
$\widetilde{\bf{S}}\left(t\right){{\bf{a}}_{m}}\left(t\right)={\bf{0}}$ is
equally rewritten as

$\sum\limits_{i=n}^{M}{S_{n,i}^{M-1}\left(t\right){a_{i,m}}\left(t\right)=0,n\in
1\cdots M},$ (62)

where $n$ denotes the $n$-th row of $\widetilde{\bf{S}}\left(t\right)$. The
solution to the series of equations (62) is

${a_{n,m}}\left(t\right)=-\sum\limits_{i=n}^{M-1}{{{\cal
K}}_{n,i+1}^{M-1}\left(t\right){a_{i+1,m}}\left(t\right),n\in 1\cdots M-1},$
(63)

where the ratio is defined by ${{\cal
K}}_{i,j}^{M-1}\left(t\right)=\frac{{S_{i,j}^{M-1}\left(t\right)}}{{S_{i,i}^{M-1}\left(t\right)}}$.
According to Appendix C, the solution ${{\bf{a}}_{m}}\left(t\right)$ is non-
unique and ${a_{M,m}}\left(t\right)$ can be any possible value.

The ratio ${{\cal K}}_{i,j}^{I-1}\left(t\right)$ is first analyzed in order to
prove the iterative property of the channel weight $a_{n,m}\left(t\right)$. We
have proved that the ratio ${{{\cal
K}}_{1,2}}\left(t\right)=\frac{{{S_{1,2}}\left(t\right)}}{{{S_{1,1}}\left(t\right)}}$
can be modeled by ${s_{1,2}}\left(t\right)$ multiplied by a time-invariant
scaling factor. Likewise, when $I=1$, the ratio ${{\cal
K}}_{i,j}^{0}\left(t\right)=\frac{{S_{i,j}^{0}\left(t\right)}}{{S_{i,i}^{0}\left(t\right)}}$
shares the same property. When $I>1$, based on the iterative equation

$S_{i,j}^{I}\left(t\right)=S_{i,j}^{I-1}\left(t\right)-{{\cal
K}}_{I,j}^{I-1}\left(t\right)S_{i,I}^{I-1}\left(t\right),i>I,j>I,$ (64)

the next iteration $S_{i,j}^{1}\left(t\right)$ can be modeled by a complex
exponential model like $S_{i,j}^{0}\left(t\right)$. As a result, the ratio
${\cal
K}_{i,j}^{1}\left(t\right)=\frac{{S_{i,j}^{1}\left(t\right)}}{{S_{i,i}^{1}\left(t\right)}}$
can be also estimated by a complex exponential model, which means that all
ratios ${{\cal K}}_{i,j}^{I-1}\left(t\right)$ hold this conclusion.

In the end, a series of solutions ${{\bf{a}}_{m}}\left(t\right)$ can be
obtained by solving (63) given the particular solution
${{\overline{a}}_{M,m}}\left(t\right)={{\cal U}}{s_{M,m}}\left(t\right)$ or
the general solution ${{a}_{M,m}}\left(t\right)$. Therefore, the weight
${a_{M,m}}\left(t\right)$ can be also modeled by a complex exponential model
when $M>2$. Based on the above discussion about cases when $M=2$ and $M>2$,
Theorem 1 is proved.

##  Appendix D Proof of Theorem 2

Given ${N_{{\rm{ce}}}}=2$, the prediction order of FMPP is $L_{\rm{ce}}=1$.
When ${N_{t}},{N_{f}}\to\infty$, Appendix B has proved that each exponential
${{e^{j{\omega_{m,n}}}}}$ of the complex amplitude
${\eta_{m,n}}\left(t\right)$ corresponds to the Doppler frequency shift
${{e^{j{\omega_{m,p}}}}}$ of the channel ${\bf{h}}_{m}\left(t\right)$. Based
on (35), the asymptotic channel PE is

$\begin{array}[]{l}\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}\mathbb{E}\left\\{{\frac{{\left\|{{{\bf{h}}_{m}}\left(t\right)-{{\hat{\bf{h}}}_{m}}\left(t\right)}\right\|_{2}^{2}}}{{\left\|{{{\bf{h}}_{m}}\left(t\right)}\right\|_{2}^{2}}}}\right\\}\\\
=\mathbb{E}\left\\{{\frac{{\left\|{\sum\limits_{p=1}^{P}{{\beta_{p,m}}{e^{j{\omega_{p,m}}t}}{{\bf{d}}_{m}}}-\sum\limits_{n=1}^{{N_{t}}{N_{f}}}{{c_{m,n}}{e^{j{\omega_{m,n}}t}}{{\bf{w}}_{n}}}}\right\|_{2}^{2}}}{{\left\|{{{\bf{h}}_{m}}\left(t\right)}\right\|_{2}^{2}}}}\right\\}\\\
=\mathbb{E}\left\\{{\frac{{\left\|{\sum\limits_{n\in{{{\cal
N}}_{{\rm{nz}}}}}{\left({\frac{{{\beta_{m,n}}{e^{j{\omega_{m,n}}t}}}}{{\sqrt{{N_{f}}{N_{t}}}}}-{c_{m,n}}{e^{j{\omega_{m,n}}t}}}\right){{\bf{w}}_{n}}}}\right\|_{2}^{2}}}{{\left\|{{{\bf{h}}_{m}}\left(t\right)}\right\|_{2}^{2}}}}\right\\}=0,\end{array}$
(65)

where ${{{{\cal N}}_{{\rm{nz}}}}}$ is the collection of non-zero indices of
the angle-delay vector ${\bf{w}}_{n}$. Thus, Theorem 2 is proved.

##  Appendix E Proof of Theorem 3

The eigenvector ${\overline{\bf{u}}_{1}\left(t\right)}$ corresponding to the
maximum eigenvalue is analyzed first. Based on (8) and (19), the eigenvector
PE is equally rewritten as

$\begin{array}[]{l}\overline{\bf{u}}_{1}\left(t\right)-{\hat{\bf{u}}_{1}}\left(t\right)=\sum\limits_{j=1}^{M}{{{\bf{h}}_{j}}\left(t\right)\left({{a_{1,j}}\left(t\right)-{{\hat{a}}_{1,j}}\left(t\right)}\right)}\\\
={\bf{H}}{\left(t\right)^{T}}\left({{{\bf{a}}_{1}}\left(t\right)-{{\hat{\bf{a}}}_{1}}\left(t\right)}\right).\end{array}$
(66)

The particular solution of the predicted channel weight is given in Appendix
C. The closed form of the channel weight is first (54) proved, i.e.,
${\hat{\bf{a}}_{1}}\left(t\right)\buildrel\Delta\over{=}{\overline{\bf{a}}_{1}}\left(t\right)$.
Then the iterative form of the channel weight (63) can be similarly derived.

According to Lemma 1 and (12), the number of the exponentials in (16) is
$N_{{\rm{no}}}^{{\rm{max}}}$. Therefore, at least $2N_{\rm{no}}^{\rm{max}}$
samples ensure a successful prediction of the
${\hat{\bf{a}}_{1}}\left(t\right)$

$\begin{array}[]{l}\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}{\hat{\bf{a}}_{1}}\left(t\right)={\left[{\begin{array}[]{*{20}{c}}{{{\hat{a}}_{1,1}}\left(t\right)}&{{{\hat{a}}_{1,2}}\left(t\right)}\end{array}}\right]^{T}}\\\
{\left[{\begin{array}[]{*{20}{c}}{\sum\limits_{l=1}^{N_{{\rm{no}}}^{\max}}{b_{1,1}^{\left(l\right)}{e^{j\omega_{1,1}^{\left(l\right)}t}}}}&{\sum\limits_{l=1}^{N_{{\rm{no}}}^{\max}}{b_{1,2}^{\left(l\right)}{e^{j\omega_{1,2}^{\left(l\right)}t}}}}\end{array}}\right]^{T}}\end{array}.$
(67)

According to Theorem 1, each amplitude ${b_{1,2}^{\left(l\right)}}$ of
${{\hat{a}}_{1,2}}\left(t\right)$ maps the amplitude ${{\rho_{1,2,x}}}$ of
${s_{1,2}}\left(t\right)$. And the poles
${{e^{j\omega_{1,2}^{\left(l\right)}}}}$ of ${{\hat{a}}_{1,2}}\left(t\right)$
maps the poles ${e^{j{\omega_{1,2,x}}}}$ of ${s_{1,2}}\left(t\right)$. Lemma 1
shows that $s_{1,2}\left(t\right)$ can be approximated by a
$N_{\rm{no}}^{\rm{max}}$-order complex exponential model (12). Therefore, the
interpolated channel weight satisfies
$\mathop{\lim}\limits_{{N_{t}}{N_{f}}{\to\infty}}{\hat{a}_{1,2}}\left(t\right){\rm{}}={{\cal
U}}{s_{1,2}}\left(t\right)$. In terms of the estimation of
${{{\hat{a}}_{1,1}}\left(t\right)}$, based on Appendix B, the poles
${{e^{j{\omega_{1,2,x}}}}}$ of $s_{1,2}\left(t\right)$ denotes the product of
the non-orthogonal Doppler frequencies of channel ${\bf{h}}_{1}\left(t\right)$
and channel ${\bf{h}}_{2}\left(t\right)$. Thus the order of ${\widetilde{{\cal
U}}{s_{1,2}}{{\left(t\right)}^{2}}}$ is also $N_{\rm{no}}^{\rm{max}}$. Then
the real channel weight ${\widetilde{{\cal
U}}{s_{1,2}}{{\left(t\right)}^{2}}}$ can also be modeled by an
$N_{\rm{no}}^{\rm{max}}$-order complex exponential model (12), which means
that the estimated channel weight
$\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}{\hat{a}_{1,1}}\left(t\right)=\sum\limits_{l=1}^{N_{{\rm{no}}}^{\max}}{b_{1,1}^{\left(l\right)}{e^{j\omega_{1,1}^{\left(l\right)}t}}}$
can map the real channel weight ${\widetilde{{\cal
U}}{s_{1,2}}{{\left(t\right)}^{2}}}$. Therefore, the error-free channel weight
estimation is achieved
$\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}{\hat{\bf{a}}_{1}}\left(t\right){\rm{}}={\overline{\bf{a}}_{1}}\left(t\right)$
as well as the eigenvector estimation
$\mathop{\lim}\limits_{{N_{t}}{N_{f}}\to\infty}{\widehat{\bf{u}}_{1}}\left(t\right)={\overline{\bf{u}}_{1}}\left(t\right)$.

Likewise, the above derivation is also valid for the other eigenvectors. In
the end, Theorem 3 is proved.

## References

  * [1] T. L. Marzetta, “Noncooperative cellular wireless with unlimited numbers of base station antennas,” _IEEE Trans. Wireless Commun._ , vol. 9, no. 11, pp. 3590–3600, 2010.
  * [2] H. Q. Ngo, E. G. Larsson, and T. L. Marzetta, “Energy and spectral efficiency of very large multiuser MIMO systems,” _IEEE Trans. Commun._ , vol. 61, no. 4, pp. 1436–1449, 2013.
  * [3] J. Hoydis, S. ten Brink, and M. Debbah, “Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?” _IEEE J. Sel. Areas Commun._ , vol. 31, no. 2, pp. 160–171, 2013.
  * [4] W. Zhou, J. Wu, and P. Fan, “High mobility wireless communications with doppler diversity: Fundamental performance limits,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 12, pp. 6981–6992, 2015.
  * [5] J. Wu and P. Fan, “A survey on high mobility wireless communications: Challenges, opportunities and solutions,” _IEEE Access_ , vol. 4, pp. 450–476, 2016.
  * [6] Fraunhofer IIS and Fraunhofer HHI, “RP-191951: Mobility enhancements for MIMO,” in _3GPP TSG RAN WG#85_ , January 2019, Newport Beach California, USA.
  * [7] Joern Krause, “RP-213565: New WID on Further NR mobility enhancements,” in _3GPP TSG RAN WG#94-e_ , Online, December 2021.
  * [8] H. Yin, H. Wang, Y. Liu _et al._ , “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, 2020.
  * [9] M. Li, S. Zhang, F. Gao _et al._ , “A new path division multiple access for the massive MIMO-OTFS networks,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 4, pp. 903–918, 2021.
  * [10] H. Jiang, M. Cui, D. W. K. Ng _et al._ , “Accurate channel prediction based on transformer: Making mobility negligible,” _IEEE J. Sel. Areas Commun._ , vol. 40, no. 9, pp. 2717–2732, 2022.
  * [11] 3GPP, _NR; Physical layer procedures for data (Release 17)_. Technical Report TR 38.214, available: http://www.3gpp.org, 2022.
  * [12] 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, 2010.
  * [13] L. Dong, H. Zhao, Y. Chen _et al._ , “Introduction on IMT-2020 5G trials in china,” _IEEE J. Sel. Areas Commun._ , vol. 35, no. 8, pp. 1849–1866, 2017.
  * [14] T. Okuyama, S. Suyama, N. Nonaka _et al._ , “Outdoor experimental trials of millimeter-wave base station cooperation with digital beamforming in high-mobility environments for 5G evolution,” in _2020 IEEE 92nd Vehicular Technology Conference (VTC2020-Fall)_ , 2020, pp. 1–5.
  * [15] Y. Chen, X. Zhan, F. Shu _et al._ , “Two low-complexity DOA estimators for massive/ultra-massive MIMO receive array,” _IEEE Wireless Commun. Lett._ , vol. 11, no. 11, pp. 2385–2389, 2022.
  * [16] C. Jeon, Z. Li, and C. Studer, “Approximate gram-matrix interpolation for wideband massive MU-MIMO systems,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 5, pp. 4677–4688, 2020.
  * [17] S. Nijhawan, A. Gupta, K. Appaiah, R. Vaze, and N. Karamchandani, “Flag manifold-based precoder interpolation techniques for MIMO-OFDM systems,” _IEEE Trans. Wireless Commun._ , vol. 69, no. 7, pp. 4347–4359, 2021.
  * [18] 3GPP, _Study on channel model for frequencies from 0.5 to 100 GHz (Release 17)_. Technical Report TR 38.901, available: http://www.3gpp.org, 2022.
  * [19] Y.-H. Nam, B. L. Ng, K. Sayana _et al._ , “Full-dimension MIMO (FD-MIMO) for next generation cellular technology,” _IEEE Commun. Mag._ , vol. 51, no. 6, pp. 172–179, 2013.
  * [20] A. M. Sayeed, “Deconstructing multiantenna fading channels,” _IEEE Trans. Signal Process._ , vol. 50, no. 10, pp. 2563–2579, 2002.
  * [21] W. Chen, K. M. Wong, and J. P. Reilly, “Detection of the number of signals: A predicted eigen-threshold approach,” _IEEE Trans. Signal Process._ , vol. 39, no. 5, pp. 1088–1098, 1991.
  * [22] D. W. Tufts and R. Kumaresan, “Estimation of frequencies of multiple sinusoids: Making linear prediction perform like maximum likelihood,” _Proc. IEEE_ , vol. 70, no. 9, pp. 975–989, 1982.
  * [23] Y. Bresler and A. Macovski, “Exact maximum likelihood parameter estimation of superimposed exponential signals in noise,” _IEEE Trans. Acoust., Speech, Signal Process._ , vol. 34, no. 5, pp. 1081–1089, 1986.
  * [24] Y. Hua and T. K. Sarkar, “Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise,” _IEEE Trans. Acoust., Speech, Signal Process._ , vol. 38, no. 5, pp. 814–824, 1990.
  * [25] Z. Qin, H. Yin, Y. Cao _et al._ , “A partial reciprocity-based channel prediction framework for FDD massive MIMO with high mobility,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 11, pp. 9638–9652, 2022.
  * [26] R. R. Nadakuditi and A. Edelman, “Sample eigenvalue based detection of high-dimensional signals in white noise using relatively few samples,” _IEEE Trans. Signal Process._ , vol. 56, no. 7, pp. 2625–2638, 2008.
  * [27] H. Yin, D. Gesbert, M. Filippou _et al._ , “A coordinated approach to channel estimation in large-scale multiple-antenna systems,” _IEEE J. Sel. Areas Commun._ , vol. 31, no. 2, pp. 264–273, 2013.
  * [28] P. Hoeher, S. Kaiser, and P. Robertson, “Two-dimensional pilot-symbol-aided channel estimation by wiener filtering,” in _1997 IEEE International Conference on Acoustics, Speech, and Signal Processing_ , vol. 3, 1997, pp. 1845–1848 vol.3.
  * [29] K. P. Schober, Renaud-Alexandre, and R. Wichman, “Improved user-specific channel estimation using geodesical interpolation at the transmitter,” _IEEE Wireless Commun. Lett._ , vol. 4, no. 2, pp. 165–168, 2015.
  * [30] M. H. Hayes, _Statistical digital signal processing and modeling_. John Wiley & Sons, 1996.
  * [31] 3GPP, _NR; Physical layer procedures for data (Release 17)_. Technical Report TR 38.104, available: http://www.3gpp.org, 2022.
  * [32] R. Kumaresan, D. Tufts, and L. Scharf, “A prony method for noisy data: Choosing the signal components and selecting the order in exponential signal models,” _Proc. IEEE_ , vol. 72, no. 2, pp. 230–233, 1984.
  * [33] T. Roos, P. Myllymaki, and J. Rissanen, “MDL denoising revisited,” _IEEE Trans. Signal Process._ , vol. 57, no. 9, pp. 3347–3360, 2009.
  * [34] M. Rahman and K.-B. Yu, “Total least squares approach for frequency estimation using linear prediction,” _IEEE Trans. Acoust., Speech, Signal Process._ , vol. 35, no. 10, pp. 1440–1454, 1987.
  * [35] J. Chen, J. Benesty, Y. Huang _et al._ , “New insights into the noise reduction wiener filter,” _IEEE/ACM Trans. Audio, Speech, Language Process._ , vol. 14, no. 4, pp. 1218–1234, 2006.