{'fraction_non_alphanumeric': 0.04278860699312789, 'fraction_numerical': 0.008838656697065307, 'mean_word_length': 4.4268793244986515, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Angular-Domain Selective Channel Tracking and Doppler Compensation for High-Mobility mmWave Massive MIMO

In this paper, we consider a mmWave massive multiple-input multiple-output (MIMO) communication system with one static base station (BS) serving a fast-moving user, both equipped with a very large array. The transmitted signal arrives at the user through multiple paths, each with a different angle of-arrival and hence Doppler frequency offset, thus resulting in a fast time-varying multipath fading MIMO channel. In order to mitigate the Doppler-induced channel aging for reduced pilot overhead, we propose a new angular-domain selective channel tracking and Doppler compensation scheme at the user side. We formulate the joint estimation of partial angular-domain channel and DFO parameters as a dynamic compressive sensing problem. Then we propose a Doppler-aware-dynamic variational Bayesian inference (DD-VBI) algorithm to solve this problem efficiently. Finally, we propose a practical DFO compensation scheme which selects the dominant paths of the fast time-varying channel for DFO compensation and thereby converts it into a slow time-varying effective channel. Compared with the existing methods, the proposed scheme can enjoy the huge array gain provided by the massive MIMO and also balance the tradeoff between the CSI signaling overhead and spatial multiplexing gain. Simulation results verify the advantages of the proposed scheme over various baseline schemes.


I. INTRODUCTION
With the development of the Fifth Generation (5G) wireless communication systems, high-mobility scenarios such as highspeed rail and Vehicle-To-everything (V2x) communications have gained increasingly more interest. Due to the highspeed relative motion between the transmitter and receiver, the transmitted signal, propagating through multiple different paths, arrives at the receiver with different Doppler frequency offsets (DFOs), thus resulting in a fast time-varying multipath fading channel. In this case, the link performance such as achievable data rate will be degraded significantly due to the Doppler-induced channel aging effect [1]. To overcome this challenge, several methods have been proposed in previous works, as elaborated below.
Direct Channel Estimation/Prediction: For line-of-sight (LOS) channels with a single LOS path, it is relatively easy to compensate the Doppler effect and resolve the channel aging issue by estimating the DFO parameter of the LOS path. However, when there are multiple different paths due to rich scattering, it is challenging to compensate the Doppler effect because different paths with different DFOs are mixed together in the received signal. Nevertheless, some works have proposed to directly estimate the fast time-varying channels in the time/frequency domain [2], [3]. Existing channel estimators can be classified into two types. The first type approximates time-varying channels using a linearly time-varying (LTV) channel model [4], [5]. For example, a hybrid frequency/timedomain channel estimation algorithm is proposed in [4] based on the LTV model and two methods are introduced to mitigate the Doppler effect. However, this algorithm introduces a processing delay of at least one Orthogonal Frequency Division Multiplexing (OFDM) symbol. The second type of estimators adopt the basis expansion model (BEM) [6] to convert the problem of estimating the channel impulse response (CIR) to that of estimating the basis function weights [7]. For example, in [7], the channel estimation and Doppler mitigation are jointly considered by exploiting the correlations in time and frequency domains, and the basis function coefficients are estimated via the linear minimum mean squared error (LMMSE) approach. However, accurate knowledge on the maximum DFO is required to determine the minimum order of basis function and the computational burden is also heavy for multi-antenna systems [8]. Moreover, the BEM inevitably introduces approximation error to channel estimation due to the imperfect model assumed.
Orthogonal Time Frequency Space (OTFS) Modulation: OTFS modulation [9] is an emerging technique which is able to handle the fast time-varying channels. This method modulates transmitted symbols in the delay-Doppler domain instead of time/frequency domain as in traditional modulation techniques such as OFDM.. The idea is to transform the timevarying channel into a time-invariant channel in the delay-Doppler domain. Early works on OTFS modulation focused on the single-input single-output (SISO) systems [9]. Later, OTFS is extended to multiple-input multiple-output (MIMO) systems by transmitting consecutive impulses with proper guard time between two adjacent ones to distinguish different base station (BS) antennas [10]. However, the channel estimation method in [10] cannot be directly applied to massive MIMO system since a large number of antennas are required to be distinguished by transmitting such impulses, which will lead to large pilot overhead.
Angular-Domain DFO Estimation and Compensation: Since the different DFOs of multiple paths are resulted from their different angles of arrival (AoAs)/angles of departure (AoDs) at the receiver, they can be separated in the angle domain via spatial processing. As such, angular-domain DFO estimation and compensation is another popular approach to address the Doppler-induced channel aging issue [1], [11]. For MIMO systems, prior work [12] first pointed out that channel time-variation can be slowed down through beamforming with a large number of transmit/receive antennas. Motivated by this, a small-scale uniform circular antenna-array (UCA) is adopted in [13], [14] to separate multiple DFOs via array beamforming. However, the DFO compensation methods in [13], [14] only apply to scenarios with very sparse channels due to the limited spatial resolution of small-scale MIMO. Recently, some works have exploited high-spatial resolution provided by massive MIMO to address the high-mobility induced challenges [15], [16], [17], [18], [19]. For example, the authors of [8] propose to separate the DFOs in angular-domain by beamforming with a large-scale uniform linear antenna-array (ULA) at the mobile user side. After estimating and compensating the DFO in each angle, the resultant quasi time-invariant channel can be estimated more efficiently. However, the signaling overhead for the maximum likelihood (ML) based joint estimation of the massive MIMO channel matrix and DFO parameters is extremely high. Moreover, only a single data stream is transmitted from the BS to the mobile user via all possible channel directions, and thus it cannot enjoy the spatial multiplexing gain as well as the huge array signal-to-noise ratio (SNR) gain provided by massive MIMO.
The above works have focused on channel estimation. It is also possible to directly search for the best beamforming vectors without explicit channel estimation. For example, in [20], the authors propose an exhaustive search (ES) scheme, which examines all beam pairs in the codebook and determines the best pair that maximizes a given performance metric (e.g., beamforming gain). To reduce the training overhead of the ES scheme, the hierarchical search (HS) scheme proposed in [21] utilizes hierarchical codebooks and has a favorable performance at low SNRs. However, these codebook-based beam search methods suffer from the quantization error caused by the codebook and the channel aging effect.
In this paper, we consider high-mobility mmWave massive MIMO systems, where both the BS and users are equipped with massive antenna arrays to facilitate Doppler compensation, improve the spectrum and energy efficiency, as well as overcome the large path loss at high frequency band. As such, combining the mmWave and massive MIMO technologies has the potential to significantly improve the capacity and reliability of high-mobility wireless communications. However, it is very challenging to design an efficient channel estimation scheme. For example, since both ends have massive MIMO, the dimension of the channel matrix is huge and conventional downlink/uplink channel estimation or codebookbased beam search will lead to large signaling overhead. Although various pilot overhead reduction methods such as those based on compressive sensing (CS) have been proposed for the estimation of slow massive MIMO fading channels [22][23], they do not consider the Doppler effect and thus cannot be applied to the high mobility scenario. To overcome this challenge, we propose a novel angular-domain selective channel tracking and Doppler compensation scheme, which exploits the dynamic sparsity of the mmWave massive MIMO channel as well as precoded training in both downlink and uplink to significantly reduce the signaling overhead. The main new contributions of our paper are given as follows.
• Angular-Domain Selective Channel Tracking: We propose a selective channel tracking scheme to only estimate partial angular-domain channel parameters at the user side that are sufficient for Doppler effect compensation to significantly reduce the pilot overhead. Moreover, we propose an efficient downlink training vector design at the BS side to strike a balance between exploitation of most promising channel directions for array (SNR) gain and exploration of unknown channel directions. Compared to the conventional random training vector design, the proposed design can exploit the massive MIMO array gain to further enhance the channel tracking performance. • Angular-Domain Selective Doppler Compensation: We propose an angular-domain selective DFO compensation scheme at the user side which selectively converts the dominant paths of the fast time-varying channel into a slow time-varying effective channel. Compared to the non-selective DFO compensation scheme in [18], the proposed scheme can enjoy the huge array (SNR) gain provided by the massive MIMO and also balance the tradeoff between the CSI signaling overhead reduction and spatial multiplexing gain maximization. • Channel Tracking Algorithm based on Dynamic VBI: To further reduce the pilot overhead, the proposed partial channel tracking design is formulated as a dynamic CS problem with unknown DFO parameters in the measurement matrix. Then, we adopt a three-layer hierarchical Markov model to capture the dynamic sparsity of the partial angular-domain channel. The existing methods, such as Variational Bayesian inference (VBI) [24] and Sparse Bayesian Learning (SBL) [25], cannot be directly applied to this three-layer hierarchical prior. To address this challenge, we propose a Doppler-aware-dynamic Variational Bayesian inference (DD-VBI) algorithm, which combines the VBI and message-passing approaches to achieve superior channel tracking performance.
The rest of this paper is organized as follows. In Section II, we describe the system model and frame structure. In Section III, we give a brief introduction of the proposed angulardomain selective channel tracking and Doppler compensation scheme. In Sections IV and V, we present the three-layer hierarchical Markov model for partial angular channel vector and the proposed DD-VBI algorithm. The simulation results and conclusions are given in Sections VI and VII. Number of propagation paths αt,q The complex path gain of the q-th path f d,t The maximum DFO ξt,q(ϑt,q) The AoD (AoA) of the q-th path ηt Rotation angle of user's antenna array The AoD/AoA off-grid vector N b Number of RF chains at the user N Number of AoA grid Table I: The key notations used in the paper.
to denote the j-th column of a matrix X. The key notations are summarized in Table I.

A. System Architecture and Frame Structure
Consider a time-division duplexing (TDD) mmWave massive MIMO system with one static BS serving a fast-moving user 1 . The BS is equipped with M ≫ 1 antennas. The user is equipped with N ≫ 1 antennas. The time is divided into frames, with each frame containing a downlink subframe and an uplink subframe, as illustrated in Fig. 1. Each subframe contains a large number of symbol durations. In the t-th downlink subframe, there are N p uniformly distributed training vectors, which are set to be the same vector denoted as v t . For convenience, we use N p to denote the symbol index set for the N p training vectors. Note that inserting N p identical training vectors uniformly in the downlink subframe facilitates the estimation of AoAs and Doppler parameters at the user, as will be explained later. Based on the estimated AoAs and Doppler parameters, the user applies a Doppler compensation matrix to mitigate the Doppler effect and essentially converts the fast time-varying channel into a slow time-varying effective channel. In the uplink subframe, there are two sets of N u p training vectors at the beginning and end of the uplink subframe, respectively. The two sets of uplink training vectors are used to estimate the slow time-varying effective channel after Doppler compensation. Specifically, the uplink transmission (e.g., beamforming and power allocation) in the t-th uplink subframe is optimized based on the slow time-varying effective channel estimated at the t-th uplink subframe. On the other hand, by making use of the channel reciprocity, the downlink transmission in the t-th downlink subframe is optimized based on the slow time-varying effective channel estimated at the end of the (t − 1)-th uplink subframe. Since the effective channel after Doppler compensation changes slowly compared to the subframe duration, such a design can effectively overcome the channel aging issue caused by the Doppler effect.

B. Doppler Multipath Channel Model
For clarity, we focus on the case when both the BS and mobile user are equipped with a half-wavelength space ULA and the channel is flat fading. To incorporate the DFO with conventional mmWave channels, the downlink channel model for the antenna pair {n t , n r } is given by [23] h nr ,ntt,i = Lt q=1 α t,q e j[2πf d,t icos(ϑt,q+ηt)+ψn t (ξt,q )+ψn r (ϑt,q)] , (1) where t stands for the frame index, i stands for the symbol index, L t is the total number of propagation paths, α t,q is the random complex path gain associated with the qth propagation path, f d,t is the maximum DFO of the tth frame, ξ t,q and ϑ t,q are the AoD and AoA of the qth path, respectively, and η t is rotation angle of the user's antenna array with respect to the moving direction in the tth frame. Here, ψ nt (ξ t,q ) and ψ nr (ϑ t,q ) represent the phase shifts induced at the n t -th transmit antenna and the n r -th receive antenna, respectively, which depend on the antenna structure, position, and the direction of the path. Note that in (1), we have implicitly assumed that the channel parameters L t , α t,q , ξ t,q , ϑ t,q , f d,t , η t are fixed within each frame but may change over different frames, which is usually true even for high-speed users [18]. However, the channel h nr,nt,t,i itself may change over different symbols at a much faster timescale due to the fast changing phase term 2πf d,t icos(ϑ t,q + η t ) caused by the Doppler effect.

C. Angular Domain Channel Representation
To obtain the angular domain channel representation, we introduce a uniform grid ofM AoDs andÑ AoAs over [0, 2π) In practice, the true AoDs/AOAs usually do not lie exactly on the grid points. In this case, there will be mismatches between the true AoDs/AOAs and the nearest grid point. To overcome this issue, we introduce an off-grid basis for the angular domain channel representation, as in [25]. Specifically, let θ T,mt,q and θ R,nt,q denote the nearest grid point to ξ t,q and ϑ t,q , respectively. We introduce the AoD off-grid vector For half-wavelength space ULAs, the array response vectors at the BS and user side are given by a T (θ) = Furthermore, defineX t ∈ CÑ ×M as the angular domain channel matrix with the (n, m)-th element given bỹ x t,n,m = α t,q , (n, m) = (n t,q , m t,q ), q = 1, 2, ..., L t 0, otherwise .
Then, for given AoA off-grid, DFO parameter , rotation angle pair ϕ t = β R,t , f d,t , η t , and AoD off-grid vector β T,t , H t,i can be expressed in a compact form as where t e jψn t (ξt,q) and e jψn r (ϑt,q) in (1) are implicitly contained in the array response matrices A R,i (ϕ t ) and A H T (β T,t ). Note that we can also define the angular domain representation for more general 2-dimensional (2D) antenna arrays. In this case, the array response vector a T (θ, φ) (or a R (θ, φ)) can be expressed as a function of the azimuth angle θ and elevation angle φ. Please refer to [26] for the details.

III. ANGULAR-DOMAIN SELECTIVE CHANNEL TRACKING
AND DOPPLER COMPENSATION In this section, we propose an efficient angular-domain selective channel tracking and Doppler compensation scheme at the user side. The proposed scheme can exploit both dynamic sparsity of mmWave massive MIMO channel and high resolution of AoA at multi-antenna mobile users to accurately estimate the downlink AoAs and maximum DFO. Using these estimated parameters, a Doppler compensation matrix is applied at the user to convert the fast time-varying channel into a slow time-varying effective channel, based on which efficient downlink/uplink transmissions can be achieved. The proposed scheme includes four key components, namely the Angular-Domain Selective Channel Tracking, Selective Doppler Compensation, Slow Time-Varying Effective Channel Estimation, and Downlink Training Vector Design. Fig. 2 illustrates a toplevel diagram of the proposed scheme and the details of each component are elaborated below. The frame index t will be omitted when there is no ambiguity.

A. Outline of Angular-Domain Selective Channel Tracking at the User
This component is used to estimate the downlink AoAs, rotation angle and maximum DFO based on the N p downlink training vectors. Thanks to the high-spatial resolution provided by the large array at the user side, the user can distinguish DFOs associated with different AOAs from multiple active paths. However, since both the BS and the user are equipped with the large array, the parameter space can be very large if we attempt to estimate the full angular-domain channel parameters (i.e., the full angular domain channel matrixX, rotation angle η and maximum DFO f d ). Since the DFO only occurs at the mobile user side, we propose to only estimate partial angular-domain channel parameters that are just sufficient to obtain AoAs, rotation angle and maximum DFO for Doppler compensation.
Specifically, the product of channel H i and downlink training vector v can be expressed as: )v are called partial angular channel coefficients since x only contain partial information about the full angular channelX. Specifically, a non-zero |x n | 2 with value larger than the noise floor indicates that there is an active path to the n-th AoA direction at the user side. Therefore, we only need to estimateÑ partial channel parameters x, the AoA off-grid vector β R , rotation angle η and maximum DFO f d , which are much less than the originalÑM full channel parameters, the off-grid vector , rotation angle and maximum DFO. Note that if the N p training vectors are different, there will beÑ N p partial angular channel coefficients, leading to a larger parameter space to be estimated. Moreover, with uniformly distributed training vectors, the phase rotation due to the Doppler term e j2πf d icos(θR,n+βR,t,n+ηt) is larger compared to the case when the N p training vectors are squeezed in the beginning of the downlink subframe, leading to a better estimation performance for the Doppler parameter f d . Therefore, such a selective channel tracking design based on N p uniformly distributed and identical training vectors can significantly reduce the number of downlink training vectors N p required to achieve accurate estimation of AoAs, rotation angle and Doppler parameters. The received baseband pilot signal is given by where v ∈ C M is the training vector for downlink channel tracking, and n i is the additive white Gaussian noise (AWGN) with each element having zero mean and variance σ 2 , respectively. The exact choice of v is postponed to Section III-D. The aggregate received pilot signal (channel measurements) of all the N p downlink pilot symbols (training vectors) in the t-th frame can be expressed in a compact form as Based on the received downlink training vectors, the user obtains the estimated partial channel parametersx,β R ,η and f d using a selective channel tracking algorithm. The detailed problem formulation and algorithm design for selective channel tracking scheme are postponed to Section IV.

B. Angular-Domain Selective Doppler Compensation at the User
This component is used to convert the fast time-varying channel into a slow time-varying effective channel after obtaining the estimated partial channel parametersx,β R ,η t andf d . We first select a set of N d most significant AoA directions with the largest energy, where the energy of the n-th AoA direction θ R,n + β R,n is defined as |x n | 2 . The parameter N d ≤ N is used to control the tradeoff between the spatial multiplexing gain and the effective CSI signaling overhead (i.e., the CSI signaling overhead required to obtain effective channel H s i in (6)). Let N d ⊆ {1, ..., N } denote the index set of the selected N d most significant AoA directions. Then, in order to mitigate the Doppler effect and perform per-AoA DFO compensation for each selected AoA direction, a DFO compensation matrix where In the following, we explain why the Doppler effect can be alleviated by applying the above DFO compensation matrix to obtain an effective channel H s i . For half-wavelength space ULA, if there is no estimation error for the partial channel parameters, we have as N → ∞ [27]. From (7), H s i is constant within a frame if we ignore the small order term O 1 √ N , i.e., the Doppler effect can be completely eliminated for sufficiently large N . This observation is also consistent with the results in [18].

C. Slow Time-Varying Effective Channel Estimation at the BS
The user can simply transmit N d orthogonal pilots in the uplink training stage. Then the conventional Least Squares (LS) based channel estimation method can be used at the BS to obtain the estimated slow time-varying effective channel H s i . Based onĤ s i , the BS can optimize the precoder for both uplink and downlink transmissions. Note that the optimization of MIMO precoder is a standard problem and there are many existing solutions with different performance and complexity tradeoff. Then, the optimized uplink precoder is fed back to the user for uplink transmission. Since N d can be much less than N , the feedback overhead for the optimized uplink precoder is acceptable for practice.

D. Training Vector Design at the BS
The training vector v t at the BS is designed according to the slow time-varying effective channelĤ s t−1 estimated at the end of the (t − 1)-th uplink subframe. The basic idea for training vector design is to strike a balance between exploitation of known channel directions (i.e., transmitting training signal over the most promising channel directions with large effective channel energy to achieve beamforming gain) and exploration of unknown channel directions (i.e., transmitting training signal over other channel directions to detect unknown channel directions). Since the effective channel H s i changes slowly, the effective channelĤ s t−1 estimated at the end of the (t − 1)-th uplink subframe is expected to provide valuable information for the most promising channel directions. On the other hand, the information about the most promising channel directions extracted fromĤ s t−1 may not be perfect due to the estimation error and CSI delay. In addition, some new direction may arise in the next frame. Therefore, the other channel directions should also be incorporated into the training vector to facilitate the detection of unknown channel directions.
Specifically, we first project the estimated effective channel where µ is a threshold which is chosen to be closed to 1. In other words, the most promising channel directions contain µ fraction of the total effective channel energy. Let N s = |M * |. Finally, the training vector is given by where the first term in (9) exploites the information about the most promising N s channel directions extracted fromĤ s t−1 , ρ is a system parameter which determines the proportion of transmit power used to exploit the most promising channel directions, the second term is used to detect the other unknown channel directions, θ s m is randomly generated from [0, 2π].  In Fig. 3, we illustrate how the achievable data rate is affected by changing the parameters ρ and µ to achieve different tradeoffs between the exploration and exploitation. It can be seen that setting µ = 0.9 and ρ = 0.5 can strike a good balance between exploitation and exploration.

A. Three-layer Hierarchical Markov Model for Partial Angular Channel Vector
The dynamic sparsity of the partial angular channel coefficients x t is captured using a three-layer hierarchical Markov model, as illustrated in Fig. 4. The first layer of random variable is the channel support vector s t ∈ {0, 1}Ñ , whose n-th element, denoted by s t,n , indicates whether the channel coefficient x t,n is active (s t,n = 1) or not (s t,n = 0). The second layer of random variable is the precision vector γ t = [γ t,1 , · · · , γ t,Ñ ] T , where γ t,n represents the precision (inverse of the variance) of x t,n . The third layer of random variables are partial angular channel coefficients x t . For convenience, denote a time series of vectors {x τ } t τ =1 as x 1:t (same for γ 1:t , s 1:t , f d,1:t ,β R,1:t ). Then the three-layer hierarchical Markov prior distribution (joint distribution of x 1:t , γ 1:t and s 1:t ) is given by where p (s 1 |s 0 ) p (s 1 ), the conditional probability p (x τ |γ τ ) has a product form p (x τ |γ τ ) = Ñ n=1 p (x τ,n |γ τ,n ) and each is modeled as a Gaussian prior distribution The conditional prior of precision vector γ τ is given by Γ (γ; a γ , b γ ) is a Gamma hyperprior. a τ , b τ are the shape and rate parameters of the channel precision γ τ,n conditioned on s τ,n = 1 and they should be chosen such that aτ bτ = E[γ τ,n ] = Θ(1), since the variance γ −1 τ,n of x τ,n is Θ(1) when it is active (s τ,n = 1). a τ , b τ are the shape and rate parameters , conditioned on the opposite event (i.e., s τ,n = 0). In this case, the shape and rate parameters a τ , b τ of the precision γ τ,n should be chosen such that aτ bτ = E[γ τ,n ] ≫ 1, since the variance γ −1 τ,n of x τ,n is close to zero when it is inactive. Note that the exact channel distribution is usually unknown in practice. In this case, it is reasonable to choose a prior distribution such that the derived algorithm can promote sparsity with low complexity and achieve robust performance to different channel distributions. By controlling the parameters in the Gamma distribution of γ τ,n , one can easily promote sparsity based on the knowledge of channel support s τ , as explained above. Moreover, Since the Gamma distribution for γ τ,n is the conjugate probability distribution of the Gaussian distribution for x τ,n , the above hierarchical prior for x τ,n and γ τ,n facilitates low-complexity VBI algorithm design with closed-form update equations [24]. Finally, the VBItype algorithm derived from the such a hierarchical prior is well known to be insensitive to the true distribution of the sparse signals [24], [28]. As a result, similar hierarchical prior distribution has been widely adopted in sparse Bayesian learning [29].
Due to the slowly changing propagation environment, the channel supports often change slowly over time, which implies that s τ,n depends on s τ −1,n , e.g., if s τ −1,n = 1, then there is a higher probability that s τ,n is also 1. Such dynamic sparsity of support vectors can be naturally modeled as a temporal Markov model with an initial prior distribution p(s 1 ) and a transition probability: where the transition probability is given by p(s τ,n = 1|s τ −1,n = 0) = ρ 0,1 , and p(s τ,n = 0|s τ −1,n = 1) = ρ 1,0 . The Markov parameters {ρ 1,0 , ρ 0,1 } characterize the degree of temporal correlation of the channel support. Specifically, smaller ρ 1,0 or ρ 0,1 lead to highly correlated supports across time, which means the propagation environment between the user and BS is changing slowly. Larger ρ 1,0 or ρ 0,1 can allow support to change substantially across time, which means the propagation environment is changing significantly. Moreover, the statistic parameters {ρ 1,0 , ρ 0,1 } could be automatically learned based on the EM framework during the recovery process [30], as detailed in Appendix A. The initial distribution p(s 1,n ), ∀n is set to be the steady state distribution of the Markov chain in (13), i.e., λ p(s 1,n ) = ρ 0,1 ρ 0,1 + ρ 1,0 .

B. Selective Channel Tracking Formulation
Using the angular domain channel representation, the receive signal at t-th frame y t ∈ C N Np can be rewritten as a CS model with an unknown AoA off-grid and DFO parameter pair ϕ t = β R,t , η t , f d,t in the measurement matrix as where the measurement matrix is given by In each frame t, the user needs to estimate the partial channel parameters x t , the AoA off-grid and DFO parameter pair ϕ t = β R,t , η t , f d,t , given the observations up to t frame y 1:t in model (14), the estimated AoA off-grid and DFO parameter pairsφ 1:t−1 = β R,1:t−1 ,η 1:t−1 ,f d,1:t−1 up to (t − 1) frame. In particular, for given ϕ t , we are interested in computing minimum mean-squared error (MMSE) estimates of x t,n ,x t,n = E x t,n |y 1:t ;φ 1:t−1 , ϕ t , where the expectation is over the marginal posterior: where v t = {x t , s t , γ t , κ t }, −x t,n denotes the vector collections integration over v t except for the element x t,n and ∝ denotes equality after scaling. On the other hand, the optimal ϕ t at the t-th frame is obtained by ML as follows [25]: Once we obtain the ML estimate ofφ t , and the associated conditional marginal posterior p(x t,n |y 1:t ;φ 1:t−1 , ϕ t ), we can obtain the MMSE estimates of {x t,n }.
One challenge in computing the MMSE estimate is the calculation of the exact posterior in (15) whose factor graph has loops. In the next subsection, we propose a Doppler-awaredynamic-VBI (DD-VBI) algorithm to approximately calculate the marginal posteriors p(x t,n |y 1:t ;φ 1:t−1 , ϕ t ) by combining the message passing and VBI approaches, and use the inexact majorization-minimization (MM) method (which is a generalization of the EM method) [25] to find an approximate solution for (16). The proposed DD-VBI algorithm is shown in the simulations to achieve a good performance.

A. Decomposition and Approximation of Joint Probability Distribution
This section is to decompose and approximate the joint probability distribution in (16) such that the joint probability distribution at the t-th frame only involves the probability density function (PDF) of the current hidden variables v t , the current observation y t , and the messagesp(s t |y 1:t−1 ,φ 1:t−1 ) passed from the previous frame, based on which a more efficient algorithm can be designed.
For simplicity, we definê In the rest of this section, we will omitφ 1:t−1 in the PDFs when there is no ambiguity.

B. Outline of the Doppler-Aware-Dynamic-VBI algorithm in Frame t
The basic idea of the DD-VBI algorithm is that, at every frame t, simultaneously approximates the marginal posterior {p(x t,n |y 1:t ; ϕ t )} and maximizes the log-likelihood ln p(y 1:t ; ϕ t ) with respect to ϕ t , based on the noisy measurements of the t-th frame and the messagesp(s t |y 1:t−1 ) passed from the previous frame. In summary, for every frame, the DD-VBI algorithm performs iterations between the following two major steps until convergence, as shown in Fig. 5.
• DD-VBI-E Step: Given ϕ t at t-th frame and messagesp(s t |y 1:t−1 ;φ 1:t−1 ) passed from the previous frame, calculate the approximate marginal posterior of p(v t |y 1:t ; ϕ t ), denoted as q(v t |y 1:t ; ϕ t ), using the sparse VBI approach, as elaborated in subsection V-D. • DD-VBI-M Step: Given q(v t |y 1:t ; ϕ t ) ≈ p(v t |y 1:t ; ϕ t ), construct a surrogate function for the objective function ln p(y 1:t ; ϕ t ), then maximize the surrogate function with respect to ϕ t as elaborated in subsection V-C. After convergence, the messagesp(s t+1 |y 1:t ;φ 1:t ) are calculated based on q(s t |y 1:t ; ϕ 1:t ) and passed to the next frame. In the following, we first elaborate the M step, which is a variation of the in-exact MM method in [25]. After that, we will elaborate how to approximately calculate the posterior p(v t |y 1:t ; ϕ t ) ≈ q(v t |y 1:t ; ϕ t ) in the E step, which is required to construct the surrogate function in the M step. Figure 5: Iinteraction between the two modules of the DD-VBI algorithm within a frame.

C. DD-VBI-M Step
It is difficult to directly maximize the log-likelihood function ln p(y 1:t ; ϕ t ), because there is no closed-form expression due to the multi-dimensional integration over v t as in (16). To make the problem tractable, in the DD-VBI-M Step, we adopt an in-exact MM method in [31], [25], which maximizes a surrogate function of ln p(y 1:t ; ϕ t ) with respect to ϕ t , to find an approximate solution of (16). Specifically, let u(ϕ t ;φ t ) be the surrogate function constructed at some fixed pointφ t , which satisfies the following properties: Inspired by the EM algorithm [31], we use the following surrogate function: where q(v t |y 1:t ; ϕ t ) is a tractable approximation of p(v t |y 1:t ; ϕ t ). In subsection V-A, we have approximated the joint probability distribution in (19) usingp(y 1:t , v t ; ϕ t ). Therefore, the surrogate function can be approximated aŝ When there is no approximation error for the associated PDFs, i.e., q(v t |y 1:t ; ϕ t ) = p(v t |y 1:t ; ϕ t ) andp(y 1:t , v t ; ϕ t ) = p(y 1:t , v t ; ϕ t ), it can be verified thatû(ϕ t ;φ t ) satisfies the properties in (18).
In the M step of the j-th iteration, we update ϕ t as where (·) j stands for the j-th iteration. In our problem,û(ϕ t ;φ t ) is a non-convex function and it is difficult to find its optimal solution. Therefore, we use a simple gradient update as in [25], i.e.

DD-VBI-E
Step performs the sparse VBI to approximate the conditional marginal posteriors p(v t |y 1:t ; ϕ t ) based on the following joint prior distribution: The corresponding approximate posterior distributions q (v t ) obtained by the sparse VBI will be given by (30)-(37).

1) Outline of Sparse VBI within a Frame:
For convenience, we use v t,n to denote an individual variable in v t . Let H = {n|∀v t,n ∈ v t }. Moreover, we use q (v t ) as a simplified notation for q (v t |y 1:t ; ϕ t ) when there is no ambiguity. The approximate conditional marginal posterior q (v t ) could be calculated by minimizing the Kullback-Leibler divergence (KLD) between p (v t |y 1:t ; ϕ t ) and q (v t ) subject to a factorized form constraint on q (v t ) as q (v t,n ) dv t,n = 1, ∀n ∈ H.
Problem A VBI is non-convex, we aim at finding a stationary solution (denoted by q * (v t )) of A VBI , as defined below.
is called a stationary solution of Problem A VBI if it satisfies all the constraints in A VBI and ∀n ∈ H, By finding a stationary solution q * (v t ) of A VBI , we could obtain the approximate posterior q * (v t,n ) ≈ p (v t,n |y 1:t ; ϕ t ) , ∀n ∈ H.
A stationary solution of A VBI can be obtained via alternately optimizing each individual density q (v t,n ) , n ∈ H. For given q (v t,l ) , ∀l = n, the optimal q (v t,n ) that minimizes the KLD in A VBI is given by (17), (28) can be approximated as Based on (29), the update equations of all variables are given in the subsequent subsections. The detailed derivation can be found in Appendix B. Note that the operator · v t,l is equivalent to · q(v t,l ) and the expectation f (v t,l ) q(v t,l ) w.r.t. its own approximate posterior is simplified as f (v t,l ) .

2) Initialization of Sparse VBI:
In order to trigger the alternating optimization (AO) algorithm, we use the following initializations for the distribution functions q (x t ) , q(γ t ) in the first iteration of every frame and q (s 1 ) in the first iteration of first frame. In the rest iterations, we initialize q (x t ),q (s t ),q(γ t ) to the (approximate) posterior calculated in the previous frame.

17: end for
Finally, the messagesp(s t+1 |y 1:t ;φ 1:t ), as the prior to the channel support, is passed to the next frame.
The overall DD-VBI algorithm is summarized in Algorithm 1. Note that in the t-th frame of DD-VBI, the contribution of the previous observations y 1:t−1 on the estimation of v t and ϕ t is summarized in the messagesp(s t |y 1:t−1 ;φ 1:t−1 ) passed from frame t−1. Therefore, in the t-th frame, there is no need to store all the observations y 1:t−1 up to frame t − 1.
Remark 2. In practical mmWave massive MIMO systems, the number of RF chains can be less than the number of antennas at both the BS and user sides to reduce the hardware cost and power consumption. The proposed scheme can be easily extended to the case with limited RF chains. For example, suppose there are only M b < M RF chains at the BS and N b < N RF chains at the user. In this case, when the BS transmits the training vector v = gF ∈ C M for downlink channel tracking in the i-th symbol duration, the user employs U i = W i G i ∈ C N ×N b as a combining matrix to combine the received signal into N b baseband channel measurements, where F ∈ C M×M b and g ∈ C M b are the RF training matrix and baseband training vector at the BS, respectively, and W i ∈ C N ×N b and G i ∈ C N b ×N b are the RF and baseband combining matrix at the mobile user in the i-th symbol duration, respectively. We can still write the received pilot signal as a CS model as in (14), but with

F. Complexity and Signaling Overhead Comparison
The computational complexity of the proposed algorithm is dominated by the update of q(x t ). Assuming the arithmetic with individual elements has complexity O(1), the computational complexity of matrix inversion in (32) Complexity order Typical complexity order Supposing the algorithm executes I iterations, the total complexity order of the proposed method is O I(N b N p ) 2Ñ , considering thatÑ is usually larger than N b N p . In Table II and III, we assume the mmWave channel has L D dominant paths and compare the complexity and signaling overhead of the proposed algorithm with the following baseline algorithms: • Baseline 1 (ML) [8]: This is the ML based joint DFO and channel estimation algorithm proposed in [8]. • Baseline 2 (ES) [20]: This is the codebook based exhaustive beam search algorithm in [20]. K E denotes the number of beamforming vectors of the codebook in ES. • Baseline 3 (HS) [21]: This is the codebook based hierarchical beam search algorithm in [21]. S denotes the total level of hierarchical codebook. K H denotes the number of beamforming vectors of each codebook level in HS. As seen from Table II [21]. In this typical case, the complexity of the proposed scheme and the baseline schemes are shown in the third column of Table  II. On the other hand, the proposed algorithm has the lowest signaling overhead for both the general case and the typical case, as shown in Table III.

VI. SIMULATION RESULTS
In this section, we compare the performance of the proposed algorithm with the baseline algorithms described in Section V-F. For the ML baseline, we also consider the case when only partial channel parameters are estimated as in the proposed scheme (ML-Partial). For the proposed scheme, we also consider the case when the training vector is generated randomly (DD-VBI-Random). The channel parameters are based on the millimeter-wave statistical spatial channel model (mm-SSCM) as specified in [32], which was developed according to the 28and 73-GHz ultrawideband propagation measurements in New York City. The signal bandwidth is 50M Hz and the frame duration is set as T b = 0.5ms. The carrier frequency is 28GHz . The mobile user employs a ULA of M = 128 antennas and the inter antenna spacing is λ/2, and the BS also employs a ULA. The user velocity is assumed to be 380km/h, which translates to f d,t ≈ 10KHz.
In the simulation, we will consider both cases when the user is equipped with a full set of RF chains and limited RF chains. For the case with limited RF chains, the number of RF chains is set to be N b = 16. The MSE for DFO estimation and the uplink achievable data rate are adopted as the performance metrics. The frequency MSE and the channel MSE is defined as x t−xt 2 xt 2 , respectively. The parameter N d is chosen to be equal to the number of dominant AoA directions at the user side (an AoA direction is called a dominant AoA direction if its energy is no less than 10% of the most significant AoA direction) and N d data streams are transmitted over the N d dominant AoA directions in the uplink with equal power allocation. Since the uplink transmission is only designed based on the dominant AoA directions which are known at the user, there is no need for the BS to feed back the slow time-varying effective channel. The Doppler frequency and channel MSE performance of different algorithms versus the pilot number, SNR and number of BS antennas are shown in Fig. 6 and Fig. 7. It can be seen that the proposed DD-VBI algorithm achieves large performance gain over all the baseline algorithms, under both full RF chains and limited RF chains. By using the proposed training vector design to strike a balance between exploitation of known channel directions and exploration of unknown channel directions, the DD-VBI algorithm could further improve the MSE performance compared to the case with a random training vector. This demonstrates that the proposed algorithm can effectively estimate the maximum DFO by selective channel tracking and efficient training vector design. Note that as the number of BS antennas increases, the parameters to be estimated increase significantly and the spatial resolution and array gain will also increase. Since the ML-Full method does not exploit the selective channel estimation method or channel sparsity to reduce the number of free parameters, its performance may degrade with the number of BS antennas. However, the performance gain of the proposed algorithm improves because in this case, the number of parameters to be estimated does not increase with the number of BS antennas. This demonstrates that the proposed algorithm is a powerful method for accurate DFO estimation, even when both BS and mobile user are equipped with massive MIMO, and the number of RF chains at the user side is limited.

B. Achievable Data Rate Performance
The achievable data rates of different algorithms versus the pilot number, SNR and number of BS antennas are shown in Fig. 8. It can be observed that the performance of all algorithms increases with the number of pilots, SNR and number of BS antennas. The proposed DD-VBI algorithm can achieve large performance gain over various baselines under both full RF chains and limited RF chains. Moreover, by using the proposed training vector design, the proposed DD-VBI algorithm could further improve the achievable data rate. This verifies that the proposed selective channel tracking and Doppler compensation scheme can also enhance the achievable data rate with low pilot overhead, under different SNRs and numbers of BS antennas.

C. Complexity versus Realized Gain
In practice, we can control the tradeoff between the complexity and realized gain of the proposed algorithm by adjusting the number of iterations. In Fig. 9, we plot the achievable rate versus the CPU time. It can be seen that the proposed algorithm can achieve a better performance than the baseline algorithms for the same CPU time. Moreover, the performance gain increases with the CPU time. Therefore, the proposed algorithm provides a better and more flexible tradeoff between the performance and computational power.

VII. CONCLUSION
We propose an angular-domain selective channel tracking and Doppler compensation scheme for high-mobility massive MIMO systems. Firstly, we propose a selective channel tracking scheme and the associated Doppler-aware-dynamic-VBI algorithm to accurately estimate the DFO and partial angulardomain channel parameters with reduced pilot overhead. Then, we propose an angular-domain selective DFO compensation Achievable data rate bps/Hz DD-VBI ML-Partial ES HS Figure 9: The tradeoff between the complexity and realized gain of the proposed algorithm.
scheme to convert the dominant paths of the fast time-varying channel into a slow time-varying effective channel, based on which efficient uplink and downlink transmissions can be achieved. Simulations verify that the proposed scheme not only can mitigate the Doppler and channel aging effect with much less pilots than existing schemes, but also can achieve a good tradeoff between the CSI signaling overhead and spatial multiplexing/array gain.

APPENDIX
A. EM Estimation of the Parameters ρ 1,0 and ρ 0,1 Denote the s t as the channel support vector, as defined in SubSec.IV-A. We apply the EM method in [31] to obtain an estimate of the parameter ρ 1,0 by perform the following iterations until convergence: ρ i+1 1,0 = arg max ρ1,0 Epi (s1:t) ln p(s 1:t , y 1:t ; ρ 1,0 )|ρ i 1,0 , (41) where ρ i 1,0 stands for the i-th iteration, Epi (s1:t) {·} denotes expectation over the posterior distributionp i (s 1:t ) = p(s 1:t |y 1:t , ρ i 1,0 ) conditioned on y 1:t and ρ i 1,0 . Solving the problem (41), we can get a closed-form expression for updating the ρ 1,0 as Since the channel support vector s 1:t forms a Markov chain, the posterior distributionp i (s 1:t ) can be approximately calculated using the sum-product message passing over the factor graph of the associated Markov chain with the priors of st ,n 's given by the output of the DD-VBI algorithm at each framet [30]. After the iterations in (42) convergence, we can obtain the EM estimation of the parameter ρ 1,0 in the t-th frame. The EM estimation of the parameter ρ 0,1 can be obtained in a similar way.