Time Varying Channel Tracking with Spatial and Temporal BEM for Massive MIMO Systems
Jianwei Zhao Hongxiang Xie Feifei Gao Weimin Jia Shi Jin Hai Lin




Time Varying Channel Tracking with Spatial and Temporal BEM for Massive MIMO Systems
    Jianwei Zhao, Hongxiang Xie,
        Feifei Gao,
         Weimin Jia, Shi Jin, and Hai Lin 
This work was presented in part at the 2017 IEEE ICC <cit.>.
J. Zhao, H. Xie and F. Gao are with  Tsinghua National Laboratory for Information Science and Technology (TNList) Beijing 100084, P. R. China, (e-mail:zhaojw15@mails.tsinghua.edu.cn, xiehx14@mails.tsinghua.edu.cn, feifeigao@ieee.org). J. Zhao is also with
High-Tech Institute of Xian, Xian, Shaanxi 710025, China.
  W. Jia is with High-Tech Institute of Xian, Xian, Shaanxi 710025, China (e-mail: jwm602@163.com).
  S. Jin is with the National Communications Research Laboratory, Southeast University, Nanjing 210096, P. R. China (email: jinshi@seu.edu.cn). 
 H. Lin is with the Department of Electrical and Information Systems,
Graduate School of Engineering, Osaka Prefecture University, Sakai,
Osaka, Japan (e-mail: hai.lin@ieee.org).


    December 30, 2023
================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================


In this paper, we propose a channel tracking method for massive MIMO systems under both time-varying and spatial-varying circumstance. Exploiting the characteristics of massive antenna array, a spatial-temporal basis expansion model (ST-BEM) is designed to reduce the effective dimensions of uplink/downlink channel, which decomposes channel state information (CSI) into the time-varying spatial information and gain information. We firstly model the users' movements as a one-order unknown Markov process, which is blindly learned by the expectation and maximization (EM) approach. Then, the uplink time-varying spatial information can be blindly tracked by unscented Kalman filter (UKF) and  Taylor series expansion of the steering vector, while the rest uplink channel gain information can be trained by only a few pilot symbols. Due to angle reciprocity (spatial reciprocity), the spatial information of the downlink channel can be immediately obtained from the uplink counterpart, which greatly reduces the complexity of downlink channel tracking. Various numerical results are provided to demonstrate the effectiveness of the proposed method.





Massive MIMO, channel tracking,  spatial and temporal BEM (ST-BEM), DOA, angle reciprocity, unscented Kalman Filter (UKF).




§ INTRODUCTION

 Massive multiple-input multiple-output (MIMO), as one of the most important techniques in 5G communications, has attracted enormous attention from both academy and industry <cit.>. It applies hundreds or even thousands of antennas at the base station (BS) to simultaneously serve tens of terminals in the same time-frequency resource, making the communications much more efficient, secure, and robust<cit.>.

The potential gains of massive MIMO depend on the channel state information (CSI), and numerous works have been proposed to solve channel estimation issues <cit.>. In most cases, the massive array antennas are adopted at base station (BS), which makes the elements of channel possess high correlation in the spatial domain, and the corresponding channel covariance appears to be low rank. Based on this fact,  <cit.> proposed to use dominant eigenvectors to span the channel vector for the time-invariant circumstance and reduce the effective channel parameters. The compressive sensing (CS) technique was proposed to reduce the training as well as the feedback overhead in <cit.>, while a novel pilot decontamination approach was presented in <cit.> that offers a powerful way of discriminating across interfering users with even strongly correlated pilot sequences. Recently, the angle domain massive MIMO channel estimation scheme were proposed in <cit.>, where the angle information is exploited to separate users, and the array signal processing method is utilized to simplify the subsequent design, also named as angle division multiple access (ADMA). The works in <cit.> all utilized the fact that the angle information could be distinguished precisely for massive MIMO systems with large number of antennas, by contrast, which would be not valid for the traditional MIMO systems, such as the smart antenna. Therefore, massive MIMO is an effective juncture of the wireless communications and the array signal processing.

Nevertheless, a practical propagation of wireless signals would face the time-varying environment. The first work considering  channel aging effect on massive MIMO was <cit.>, where the channel variation was characterized as a function of different system parameters, based on which a channel prediction method  was designed. In <cit.>, the authors investigated the impact of the general channel aging conditions on the downlink performance of massive MIMO systems, while reference <cit.> discussed the effect of channel aging on the sum rate of uplink massive MIMO systems. However, most of these works   focus on the performance analysis rather than providing a concrete approach  for  dynamic channel tracking.

Conventionally, three approaches can be used for channel aging problem: (1) the Gauss-Markov model <cit.>, which captures the channel temporal variation through the symbol-by-symbol updating; (2) the basis expansion model (BEM) <cit.>, which decomposes the channels into the superposition of the time-varying basis functions weighted by the time-invariant coefficients; (3) the temporal channel covariance matrix method, where the dominant eigenvectors of the temporal covariance matrix serve as the basis vectors to span the time-varying channels. However, there are no related channel tracking methods for large-scale massive MIMO systems to the best of our knowledge.

In this paper, we propose a channel tracking method for massive MIMO systems under both the time-varying and the spatial-varying circumstance. A spatial-temporal basis expansion model (ST-BEM) is designed to reduce the effective dimensions of uplink/downlink channels, which decomposes channel state information (CSI) into the time-varying spatial information and gain information. Besides, the spatial information can be further determined by the central direction of arrival (DOA) and angular spread (AS) of the incoming signal. We firstly model the users' movements as a one-order Markov process whose unknown parameters are estimated by expectation and maximization (EM) learning. Then, the central DOAs can be blindly tracked by unscented Kalman Filter (UKF), while the AS can also be blindly obtained through the Taylor series expansion of the steering
vector. With the tracked spatial information, the uplink channel gain can be obtained by only a few pilot symbols. Besides, due to angle reciprocity, the downlink channel spatial information can be immediately obtained from uplink counterpart, which greatly reduces the complexity of downlink channel tracking. Finally, various numerical results are provided to demonstrate the effectiveness
of the proposed method. Compared to the conference version <cit.>, this journal version has more contributions as:

  
  * The users' movements are modeled as a one-order unknown Markov process which is a more practical case of mobile users.
  
  * The DOAs of the incident signals are assumed as time varying and are tracked by EM learning based UKF.
  
  * The ADMA scheme is utilized to decrease the user interference and improve the spectrum efficiency.
  
  * Various numerical results are provided to verify the effectiveness of the proposed method.


The rest of the paper is organized as follows. In section <ref>, the system model and the channel model are given.
The problem formulation of channel tracking is described in section <ref>. Spatial signature tracking is discussed in section <ref>, while channel gain estimation is presented in section <ref>. Numerical simulations and results are displayed in section <ref>. Finally, conclusions are drawn in section <ref>.

Notations: Vector is denoted by boldface small letter, while matrix is expressed by
capital letter;  𝐀^T, 𝐀^*, 𝐀^H, 𝐀^-1 and 𝐀^† represents the transpose, complex conjugate, hermitian,
inverse, and pseudo-inverse of the matrix 𝐀 respectively; the trace of 𝐀 is tr(𝐀); the (i,j)th entry of 𝐀 is [𝐀]_ij; a diagonal matrix
with the diagonal elements constructed from 𝐚 is denoted by diag{𝐚}, while diag{𝐀}
is a vector with elements extracted from the diagonal components of matrices 𝐀;
𝔼{·} is the statistical
expectation and 𝐡 is the Euclidean norm of 𝐡; ⌈·⌉ and ⌊·⌋ denote the integer ceiling and integer floor, respectively; [·]_:,ℬ_k and [·]_ℬ_k,: represent the
sub-matrices by collecting the related ℬ_k columns or rows, respectively.



§ SYSTEM AND CHANNEL MODEL



Consider a multiuser massive MIMO system shown in Fig. <ref>, where base station (BS) is equipped with M antennas in the form of uniform linear array (ULA), and K users with single antenna are randomly distributed in the cell.  The baseband channels between users and BS are assumed to be time-selective flat-fading in each block ζ that contains N information symbols, and the corresponding M× 1 uplink channel between user k and BS at time slot n of block ζ can be expressed as <cit.>

    𝐡_k(n)=1/√(P)∑_p=1^P α_kpe^-j(2π f_d nT_scosφ_kp+ϕ_kp)𝐚[θ_kp(ζ)], 0≤ n≤ N-1,

where P is the number of the multi-incoming rays, α_kp is the complex gain of the p-th ray, f_d and T_s are the maximum Doppler frequency and the system sampling period respectively, φ_kp is the angle between uplink transmitted signal and the motion direction of user-k, and ϕ_kp represents  the initial phase, which is uniformly distributed in [0,2π]. Moreover, 𝐚[θ_kp(ζ)]∈ℂ^M× 1 is the steering vector and its m-th element can be expressed as

    {𝐚[θ_kp(ζ)]}_m=e^j2π md/λsinθ_kp(ζ),m=0,…,M-1,

where d represents the antenna spacing, λ is the signal carrier wavelength, and θ_kp(ζ) represents DOA of the p-th ray of user k seen by BS in block ζ.
Denote the central DOA of user k in block ζ as θ̅_k(ζ), and then each ray of user k can be expressed as

    θ_kp(ζ)=θ̅_k(ζ)+θ̃_kp(ζ),

where θ̃_kp(ζ) is the corresponding random AS satisfying |θ̃_kp(ζ)|≤Δθ_k(ζ), and Δθ_k(ζ) is the maximum AS of user k, which is small for massive MIMO systems.


Moreover, θ̃_kp(ζ) could be assumed as independent and identically distributed random variables satisfying <cit.>

    𝔼{θ̃_kp(ζ)θ̃_k̃p̃(̃ζ̃)̃}=σ_k^2δ(k-k̃)δ(p-p̃)δ(ζ-ζ̃),

where σ_k^2 is the variance of θ̃_kp(ζ). Since the spatial location of the users changes on the order of seconds, the DOA information of users seen by the BS can be viewed as unchanged in each data block ζ, but  might change from block to block.







§ PROBLEM FORMULATION
 
For array with large number of antennas, channel vector would exhibit many new features in both spatial domain and temporal  domain, which can be used for channel tracking.


 §.§ Spatial Domain Channel Representation with SBEM



Let us define the M-point DFT of the channel vector 𝐡̃_k(n) as

    𝐡̃_k(n)=𝐅𝐡_k(n),         k=1,…,K,

where 𝐅 is the M× M DFT matrix with (p,q)-th element [𝐅]_pq=1/√(M)e^-j2π/Mpq.

According to (<ref>), the q-th component of 𝐡̃_k(n) can be computed as

    [𝐡̃_k(n) ]_q   =
    1/√(MP)∑_p=1^P∑_m=0^M-1α_kpe^ϑ_kpe^-j[2π/Mmq-2π/λmd sinθ_kp(ζ)]
       =1/√(MP)∑_p=1^Pα_kpe^(ϑ_kp-jM-1/2η_kp)sin(Mη_kp/2)/sin(η_kp/2),

where ϑ_kp=-j(2π f_d nT_scosφ_kp+ϕ_kp) and η_kp=2π/Mq-2π d/λsinθ_kp(ζ).

Due to both the narrow AS and the large number of antennas, 𝐡̃_k(n) is a highly sparse vector and most power is contained in a small set, which is called as  spatial signature index (SSI) set and is denoted as ℬ_k(ζ). Besides, the left bound of ℬ_k(ζ) is determined by the leftmost ray θ̅_k(ζ)-Δθ_k(ζ), while the right bound of ℬ_k(ζ) is determined by the rightmost ray θ̅_k(ζ)+Δθ_k(ζ). When M→∞, the size of ℬ_k(ζ) is determined by the central DOA θ̅_k(ζ) and AS Δθ_k(ζ)<cit.> as

    B_k=  |ℬ_k(ζ)|   ≈⌈ Md/λsin[θ̅_k(ζ)+Δθ_k(ζ)] ⌉-⌊ Md/λsin[θ̅_k(ζ)-Δθ_k(ζ)] ⌋
       ≈⌈2Md/λ|cosθ̅_k(ζ)|Δθ_k(ζ)⌉,

Moreover, the central index of  ℬ_k(ζ), denoted as  q_c,k(ζ),  is determined purely by the central DOA θ̅_k(ζ) as

    q_c,k(ζ)≈⌊M d/λsinθ̅_k(ζ)⌉.



An example of with AS [27^∘,29^∘] is given in Fig. <ref> for M=128, where the DFT of entire multi-rays is illustrated. The SSI  ℬ_k(ζ) containing η=95% of the total power is 11, which is very small in comparison with  M.



Therefore,  the channel vector 𝐡_k(n) could be expanded by a limited number of orthogonal basis as

    𝐡_k(n)   =𝐅^H𝐡̃_k(n)
      ≈[𝐅^H]_:,ℬ_k(ζ)[𝐡̃_k(n)]_ℬ_k(ζ),:
        =∑_q∈ℬ_k(ζ)h̃_k,q(n) 𝐟_q,

where h̃_k,q(n)≜ [𝐡̃_k(n)]_q,  and 𝐟_q
is the orthogonal basis coming from the q-th column of 𝐅^H.  Equation (<ref>) is named as spatial basis expansion model (SBEM) <cit.>.




 §.§ Time Domain Channel Representation with CE-BEM


The components of 𝐡̃_k(n), h̃_k,q(n), q=0,…,M-1 are band-limited, and the maximum bandwidth of the power spectra is equal to the maximum Doppler frequency f_d for massive MIMO system.



Let us define the time domain discrete correlation matrix of 𝐡_k(n) as 𝐑_k(m)≜𝔼{𝐡_k(n)𝐡_k^H(n+m)}, and the (i,l)-th component of 𝐑_k(m) is given as

    [𝐑_k(m)]_i,l   =𝔼{h_k,i(n)h^*_k,l(n+m)}
       =1/P∑_p=1^P 𝔼{|α_kp|^2}𝔼{e^-j2π f_d mT_s cosφ_kp}·𝔼{e^j2π d/λ(i-l)sinθ_kp}.



The proof process is different from the conventional Clarke's reference model <cit.>, where the incoming angles of user seen by BS are assumed to be uniformly distributed in [0,2π], the AS of user k keeps within a small range, i.e., [θ_k(ζ)-Δθ_k(ζ),θ_k(ζ)+Δθ_k(ζ)], while φ_kp is assumed to be randomly distributed at [0,2π] for the random mobility of users. Consequently, the expectation in (<ref>) is mainly focused on φ_kp, and then it holds that

    [𝐑_k(m)]_i,l   =J_0(2π f_d m T_s)· g(2π d/λ(l-i)),

where J_0(x) is the first kind zero-order Bessel functions, which is given by

    J_0(x)=1/2π∫_-π^π e^-jxcos ydy,

while g(x) is defined as

    g(x)=1/2Δθ_k∫_θ_k-Δθ_k^θ_k+Δθ_ke^-jxsin ydy.

Note that g(x) does not affect the bandwidth of (<ref>). Besides, according to <cit.>, the power spectrum of J_0(2π f_d m T_s) is the “U-shape” function,
  
    S_J_0(f)=1/π f_d√(1-f^2/f_d^2),   f∈[-f_d,f_d],

and hence the bandwidth of (<ref>) is a constant f_d. Interestingly,  the range of the AS (spatial property) will not affect time-domain bandwidth of h_k,q(n) (frequency property).

Since h̃_k,q(n)=∑_i=0^M-1h_k,i(n)e^-j2π/Miq, we have

    𝔼{h̃_k,q(n)[h̃_k,q(n+m)]^*}
    =∑_i=0^M-1∑_l=0^M-1𝔼{h_k,i(n)[h_k,l(n+m)]^*}e^-j2π/Mq(i-l),

where 𝔼{h̃_k,q(n)[h̃_k,q(n+m)]^*} is the superposition of multiple band-limited signals with the same bandwidth f_d. Therefore, the components h̃_k,q(n), q=0,…,M-1 are all band-limited and upper bounded by f_d.


From Property <ref>, we know h̃_k,q(n) could be timely expanded by μ≪ M limited orthogonal time basis to capture the rapid variation of h̃_k,q(n) for massive MIMO systems <cit.> as

    h̃_k,q(n)=∑_r=0^μγ_k,q^re^j2π(r-μ/2)n/N,   0≤ n≤ N-1,

where γ_k,q^r's are the invariant coefficients, while the order μ is a function of the
channel bandwidth f_d and the length of sampling interval NT_s. Equation (<ref>) is also known as the complex exponential basis expansion model (CE-BEM).


 §.§ Channel  Dimension Reduction with ST-BEM



According to previous discussions, the overall time-varying channel vector 𝐡_k(n) can be jointly expanded as

    𝐡_k(n)    ≈∑_q∈ℬ_k(ζ)h̃_k,q(n)𝐟_q=
      ∑_q∈ℬ_k(ζ)∑_r=0^μγ_k,q^r e^j2π(r-μ/2)n/N𝐟_q
       =∑_q∈ℬ_k(ζ)γ_k,q^T 𝐜_n 𝐟_q,     n=0,…,N-1,

where γ_k,q=[γ_k,q^0,…,γ_k,q^μ]^T and
𝐜_n=[e^-j2π n/Nμ/2,…,e^j2π n/Nμ/2]^T. Moreover, (<ref>) can be named as spatial-temporal basis expansion model (ST-BEM), since it represents the channel as an combination of the SBEM and CE-BEM.

It is worth noticing that μ should satisfy  μ≥ 2⌈ f_dNT_s⌉ to offer sufficient degrees of freedom <cit.>. In order to provide an explicit illustration, an example of h̃_k,2(n) is offered in Fig. <ref>, where the simulation parameters are chosen as M=128,d=λ/2, f_d=200Hz, T_s=0.1 ms and N=100. It can be found that when μ≥ 2⌈ f_dNT_s⌉=4, the approximation of h̃_k,2(n) is well, while the ambiguous estimation happens for μ=2 due to the lack of the sufficient sampling degrees of freedom.




§ TRACKING THE SPATIAL INFORMATION WITH PARAMETER LEARNING
 
In this section, we provide a blind uplink spatial information tracking method based on the uplink received symbols. With (<ref>), channel tracking problem is transmitted to tracking SSI set ℬ_k(ζ) and estimating the CE-BEM coefficients γ_k,q.
Furthermore, since ℬ_k(ζ) is determined by central DOA  θ̅_k(ζ) and AS Δθ_k(ζ) according to (<ref>) and (<ref>), tracking the spatial information can be decomposed into tracking θ̅_k(ζ) and Δθ_k(ζ) separately.


 §.§ Central DOA θ̅_k(ζ) Modelling 

According to <cit.>, users that can be simultaneously scheduled should have distinct spatial information to decrease user interference. Therefore, users can be allocated into the same group if their SSI sets do not overlap and also guarded by a interval Ω, i.e.,

    ℬ_k(ζ)∩ℬ_l(ζ)=∅,  dist[ℬ_k(ζ),ℬ_l(ζ)]≥Ω,

Note that the users are scheduled according to their DOA information, and hence the corresponding transmission scheme can be named as angle division multiple access (ADMA). One can find much more detailed descriptions and methodologies
of ADMA in <cit.>. The initial DOA information of each user can be derived from method in  <cit.>. Interested readers could refer to <cit.> for more details.

Assuming that all K users are divided into G groups, and the g-th group is denoted as 𝒰_g with K_g users. When users in group 𝒰_g are scheduled, the received signals at BS in time slot n of block ζ can be expressed as

    𝐱(n)=∑_k=1^K_g𝐡_k(n)s_k(n)+𝐰(n)=∑_k=1^K_g s_k(n)/√(P)∑_p=1^P α_kpe^ξ_kp𝐚[θ_kp(ζ)]+𝐰(n),

where s_k(n) denotes the transmitted signal of user k, and 𝐰(n) is the noise vector whose elements are i.i.d. with distribution 𝒞𝒩(0,σ_n^2).

Since the data transmission also possesses the  Vandermonde structure of the channel <cit.>, the DFT of 𝐱̃(n) can still tell the spatial signature of the instantaneous 𝐡_k(n). According to (<ref>), it holds that

    [𝐱̃(n)]_ℬ_k(ζ)=[s_k(n)𝐅𝐡_k(n)]_ℬ_k(ζ).


Note that [𝐱̃(n)]_ℬ_k(ζ) serves a humble observation of the spatial signature. Then a straightforward way to track the SSI set ℬ_k(ζ) is to obtain the size of ℬ_k(ζ) from the size of [𝐱̃(n)]_ℬ_k(ζ), while to derive central SSI q_c,k(ζ) of ℬ_k(ζ) from the central SSI of [𝐱̃(n)]_ℬ_k(ζ).  The corresponding method can be named as DFT searching. However, when M is not infinite, power leakage effect <cit.> would actually distort the central position as well as the size of ℬ_k(ζ). Therefore, we will propose a more effective way to track the spatial information.

From (<ref>), we can derive the observed central SSI q_c,k(ζ), which provides us a practical measurement of central DOA  θ̅_k(ζ). According to (<ref>), the relationship between q_c,k(ζ) and θ̅_k(ζ) can be represented by

    q_c,k(ζ) =dM/λsinθ̅_k(ζ)+u_k(ζ),

where u_k(ζ) is the measurement noise meeting i.i.d. Gaussian distribution 𝒞𝒩(0,Q_u_k). Equation (<ref>) can be viewed as the measurement equation.

Besides, the one-order Markov process <cit.> could be utilized to describe the variation of central DOAs as

    θ̅_k(ζ) =θ̅_k(ζ-1)+ω_k(ζ-1),

where θ̅_k(ζ-1) are the central DOA in block ζ-1. Moreover, ω_k(ζ-1) is the system noise, which meets 𝔼[ω_k(ζ-1)^Hω_k(ζ-1)]=Q_ω_k. Equation (<ref>) can be viewed as the system equation.



 §.§ EM based Parameter Learning

Before proceeding, we need first derive the unknown parameters Q_ω_k and Q_u_k. Denote  β_k=[ Q_ω_k,Q_u_k] and 𝐪_k=[q_c,k(0),q_c,k(1),…,q_c,k(ς-1)]^T as the unknown parameter vector and the observation vector, respectively, where ς is the dimension of the observation vector for parameter learning. A data-aided maximum likelihood (DA-ML) estimator for β_k can be formulated as

    β̂_k=   max p(𝐪_k|β_k).

Nevertheless, the direct ML estimation of β_k is not feasible for its large computation overhead. An effective solution is to search for the ML solution iteratively via the EM algorithm, which consists of two steps, namely the expectation step and the maximization step. Besides, it has been verified that β̂_k^(l) can converge to one stationary point of the likelihood function p(ℛ|β_k) under fairly general conditions <cit.>.

Let us denote θ̅_k=[θ̅_k(0),…,θ̅_k(ς-1)]. Using the EM algorithm, the vector 𝐪_k forms
the incomplete data set ℛ, while the vectors  𝐪_k together with θ̅_k form
the complete data set 𝒵. During the l-th iteration, we firstly compute the objective function L(β_k,β̂_̂k̂^(l-1)) in the expectation step as

    L(β_k,β̂_̂k̂^(l-1))=∫_𝒵p(𝒵|ℛ,β̂_̂k̂^(l-1))ln p(𝒵|β_k) d𝒵.


Note that L(β_k,β̂_̂k̂^(l-1)) depends on the estimation β̂^(l-1) in the (l-1)th iteration and the trial value β_k. Since the set ℛ is known, the function L(β_k,β̂_̂k̂^(l-1))  can be
rewritten as

    L(β_k,β̂_̂k̂^(l-1))=   ∫_θ̅_k
    p(θ̅_k|ℛ,β̂_k^(l-1))ln[p(ℛ|θ̅_k,β_k)p(θ̅_k|β_k)]dθ̅_k
    
    =   ∫_θ̅_k
    p(θ̅_k|ℛ,β̂_k^(l-1))[
    ln p(ℛ|θ̅_k,β_k)
    +ln p(θ̅_k|β_k)]dθ̅_k
    
    
    
    
    =   𝔼{ln p(ℛ|θ̅_k,β_k)|
    ℛ,β̂_k^(l-1)}+
    𝔼{ln p(θ̅_k|β_k)|
    ℛ,β̂_k^(l-1)}.


According to the system equation (<ref>) and measurement equation (<ref>), we have


    ln p(𝐪_k|θ̅_k,β_k)   =∑_ζ=0^ς-1ln p(𝐪_𝐤_ζ|θ̅_kζ,β_k)
       =∑_ζ=0^ς-1ln1/√(2π Q_u_k)exp-(q_c,k(ζ) -dM/λsinθ̅_k(ζ))^2/2 Q_u_k,
    ln p(θ̅_k|β_k)   =-1/2ςln 2π Q_ω_k-∑_ζ=0^ς-1(θ̅_k(ζ) -θ̅_k(ζ-1))^2/2 Q_ω_k.

Therefore, it can be further derived that

    𝔼{ln p(𝐪_k|θ̅_k,β_k)|𝐪_k,β̂_k^(l-1)}=-1/2ςln 2π Q_u_k^(l-1)-∑_ζ=0^ς-1[q_c,k(ζ)]^2/2Q_u_k^(l-1)+
       ∑_ζ=0^ς-1[q_c,k(ζ)]dM/λ𝔼{sinθ̅_k(ζ)| 𝐪_k,β̂_k^(l-1)}/Q_u_k^(l-1)-∑_ζ=0^ς-1[dM/λ]^2E{sinθ̅_k^2(ζ)| 𝐪_k,β̂_k^(l-1)}/2Q_u_k^(l-1),
       𝔼{p(θ̅_k|β_k)|𝐪_k,β̂_k^(l-1)}
    =-1/2ςln 2π Q_ω_k^(l-1)-∑_ζ=0^ς-1𝔼{θ̅_k^2(ζ)| 𝐪_k,β̂_k^(l-1)}/2 Q_ω_k^(l-1)
       -∑_ζ=0^ς-1𝔼{θ̅_k^2(ζ-1)| 𝐪_k,β̂_k^(l-1)}/2 Q_ω_k^(l-1)+∑_ζ=0^ς-1𝔼{θ̅_k(ζ)θ̅_k(ζ-1)| 𝐪_k,β̂_k^(l-1)}/Q_ω_k^(l-1).


Secondly, we turn to the maximization of L(β_k,β̂_̂k̂^(l-1)) as

    β̂_k^(l)=max_β_k{ L(β_k,β̂_̂k̂^(l-1))}.


 The solution of (<ref>) can be derived from ∂ L(β_k,β̂_̂k̂^(l-1))/∂ Q_ω_k=0 and ∂ L(β_k,β̂_̂k̂^(l-1))/∂ Q_u_k=0 as

    Q_u_k=   1/ς∑_ζ=0^ς-1[q_c,k(ζ)]^2+1/ς∑_ζ=0^ς-1[dM/λ]^2E{sinθ̅_k^2(ζ)| 𝐪_k,β̂_k^(l-1)}
       -2/ς∑_ζ=0^ς-1q_c,k(ζ)dM/λ𝔼{sinθ̅_k(ζ)| 𝐪_k,β̂_k^(l-1)},
    Q_ω_k=   1/ς∑_ζ=0^ς-1𝔼{θ̅_k^2(ζ)| 𝐪_k,β̂_k^(l-1)}+1/ς∑_ζ=0^ς-1𝔼{θ̅_k^2(ζ-1)| 𝐪_k,β̂_k^(l-1)}
       -2/ς∑_ζ=0^ς-1𝔼{θ̅_k(ζ)θ̅_k(ζ-1)| 𝐪_k,β̂_k^(l-1)}.


When we can derive the statistics of the system noise and measurement noise in (<ref>) and (<ref>), the central DOA can be tracked by recursive predicting and updating with UKF.



 §.§ Tracking Central DOA θ̅_k(ζ) with UKF and URTSS

Based on Bayesian filtering and smoothing framework <cit.>, the tracking of the central DOA can then be partitioned into two phases: the forward tracking and the backward smoothing. During the former phase, Bayesian filter is adopted to sequentially estimate the current state, while during the latter phase, the corresponding Bayesian smoother is utilized to reconstruct the system states. Since the observation vector 𝐪_k is not linear with respect to the central DOA vector θ̅_k, the dynamic state space is not a linear Gaussian state space either. Thus,
the conventional Kalman filter (KF) and Rauch-Tung-Striebel smoother (RTSS) can not be utilized. Though extended Kalman filter (EKF) is a widely adopted nonlinear filter but cannot be used here since Taylor series expansion of EKF has large truncation error, which would severely degrade the channel tracking performance. In the following, we would utilize unscented Kalman filter (UKF) to implement forward tracking and unscented Rauch-Tung-Striebel smoother (URTSS) to achieve the backward smoothing respectively for each iteration  l <cit.>.

Let us  define the following posteriori statistics to describe the distribution of the system states:
    θ̅_k^s,(l)(ζ)=   𝔼{θ̅_k(ζ)|𝐪_k,β_k^(l-1)},
    
    P_k^s,(l)(ζ)=   𝔼{[θ̅_k(ζ)-θ̅_k^s,(l)(ζ)][θ̅_k(ζ)-θ̅_k^s,(l)(ζ)]^H|𝐪_k,β_k^(l-1)},
     C_k,ζ-1,ζ^s,(l)=   𝔼{[θ̅_k(ζ-1)-θ̅_k^s,(l)(ζ-1)][θ̅_k(ζ)-θ̅_k^s,(l)(ζ)]^H|𝐪_k,β_k^(l-1)}.


Besides, it can be further derived that

    𝔼{θ̅_k(ζ)θ̅_k^H(ζ)|𝐪_k,β_k^(l-1)}
    =   θ̅_k^s,(l)(ζ)[θ̅_k^s,(l)(ζ)]^H+
     P_k^s,(l)(ζ)
    ,
    𝔼{θ̅_k(ζ-1)θ̅_k^H(ζ)|𝐪_k,β_k^(l-1)}
    =   θ̅_k^s,(l)(ζ-1)[θ̅_k^s,(l)(ζ)]^H+
     C_k,ζ-1,ζ^s,(l)
    .


Now, let us start estimating the posteriori statistics
θ̅_k^s,(l)(ζ), P_k^s,(l)(ζ), and C_k,ζ-1,ζ^s,(l)
conditioned on 𝐪_k and θ̅_k^(l-1).



  §.§.§ Forward Tracking by UKF


 According to the procedure of UKF, a total of 2R+1 sigma points χ_k^(i) are calculated from  √( P_k(ζ-1)) and θ̅_k(ζ-1) in block ζ-1 to describe the central DOA distribution of each user for massive MIMO systems.

    {[                           χ^(0)_k= θ̅_k(ζ-1),i = 0,;    χ^(i)_k= θ̅_k(ζ-1)+[√(R+ε)√( P_k(ζ-1))],i=1,…,R,; χ^(i)_k= θ̅_k(ζ-1)-[√(R+ε)√( P_k(ζ-1))],i=R+1,…,2R, ].

where ε=α^2(R+κ)-R, α is a scaling parameter controlling the spread of sigma points around θ̅_k(ζ-1), R is the number of the system states, and κ is a secondary scaling parameter, which is usually set as 3-R.

Next, the system equation is applied to all the sigma points χ^(i)_k, i=0,…,2R to yield a cloud of the transformed posterior
points as

    ι_k^(i) = f(χ_k^(i)),

where f is the system state transformed function determined by (<ref>). Besides, the transformed points  ι_k^(i) can be viewed as predicted sigma points of the central DOA.

Then, the statistics of the posterior points ι_k^(i)  are calculated to form the transformed mean θ̅_k^ -(ζ)  and covariance P_k^-(ζ) for central DOA tracking as

    θ̅_k^ -(ζ)  = ∑_i = 0^2R W_m^(i)ι_k^(i),
       P_k^-(ζ) = ∑_i=0^2RW_c^(i)[ι_k^(i)- θ̅_k^ -(ζ)][ι_k^(i)- θ̅_k^ -(ζ)]^H + Q_ω_k,

where W_m^(i) and W_c^(i) are the weights for the mean and covariance respectively. Moreover, these weights can be derived as

    {[         W_m^(0)=ε/R+ε,i = 0,;   W_m^(i)=ε/2(R+ε),i=1,…,2R,; W_c^(0)=ε/R+ε+1-α^2+β,i = 0,;   W_c^(i)=ε/2(R+ε),i=1,…,2R, ].

where β is used to incorporate the prior knowledge of the distribution of the system states, and the optimal choice of β for Gaussian distribution is normally taken as 2 <cit.>.

Next, we turn to the prediction of the measurement states. Similarly, unscented transformation is used to approximate a probability distribution with the predicted θ̅_k^-(ζ) and P_k^-(ζ), and the predicted sigma points can be derived as

    {[                             χ^(0)-_k= θ̅_k^-(ζ-1),i = 0,;    χ^(i)-_k= θ̅_k^-(ζ-1)+[√(R+ε)√( P_k^-(ζ-1))],i=1,…,R,; χ^(i)-_k= θ̅_k^-(ζ-1)-[√(R+ε)√( P_k^-(ζ-1))],i=R+1,…,2R, ].

such that their mean and covariance are θ̅_k^-(ζ) and P_k^-(ζ) of the predicted central DOA.

Then, the measurement equation is applied to χ_k^(i) to form measurement sigma points ξ_k^(i) as

    ξ_k^(i) = ⌈dM/λsinχ_k^(i)-⌋.

Since ξ_k^(i) represents the possible candidates of the predicted central SSI for the central DOA tracking, we can use the mean of ξ_k^(i) to represent the predicted central SSI as

    y_k^ -  = ∑_i = 0^2R+1 W_m^(i)ξ_k^(i).


Then, the forward Kalman gain for UKF can be derived as

    K_k = P_xyP_yy^ - 1,

where P_xy and P_yy are the predicted state-measurement cross covariance and measurement covariance
respectively with the expressions:

    P_xy = ∑_i = 0^2RW_c^(i)[ι_k^(i)- θ̅_k^ -(ζ)][ ξ_k^(i) -y_k^- ]^T,
       P_yy = ∑_i = 0^2RW_c^(i)[ξ_k^(i) - y_k^ - ]( ξ_k^(i)-y_k^-)^T + Q_u_k.


At the last step of UKF, the filtered system states and the corresponding error covariance for central DOA tracking can be respectively obtained by

    θ̅_k(ζ)= θ̅_k^-(ζ)+K_k[q_c,k(ζ) - y_k^- ],
       P_k(ζ) =P_k^-(ζ)-K_k P_yyK_k^T.



The complexity of the UKF based forward DOA tracking can be evaluated by the calculation of the K_g × K_g square roots of P_k^-(ζ), P_k(ζ) in (<ref>)-(<ref>) through the Cholesky factorization operation, whose complexity  is  𝒪(K_g^3/6). In comparison, the covariance based channel tracking method needs to perform singular value decomposition, and the corresponding complexity is 𝒪(M^3) (K_g≪ M). Therefore, the proposed  channel tracking method has a relatively lower complexity.




  §.§.§ Backward Smoothing with URTSS

 UKF can only derive the system states based on the previous measurement, but not the whole trajectory of the measurements. When we get the forward estimated values of θ̅_k(ζ) and P_k(ζ), we could use URTSS to obtain a more desirable estimation and the posteriori statistic of the central DOA  θ̅_k^s(ζ), ζ=ς,…,0.

The 2R+1 sigma points χ^s,(i)_k can be derived from  the forward UKF as

    {[                         χ^s,(0)_k=θ̅_k(ζ),i = 0,;    χ^s,(i)_k=θ̅_k(ζ)+[√(R+ε)√( P_k(ζ))],i=1,…,R,; χ^s,(i)_k=θ̅_k(ζ)-[√(R+ε)√( P_k(ζ))],i=R+1,…,2R. ].


Then, the transformed sigma points ι_k^s,(i) and the predicted mean y_k^ s- of the central DOA can be expressed as

    ι_k^s,(i)= f(χ^s,(i)_k),
       y_k^ s-  = ∑_i = 0^2R+1 W_m^(i)ι_k^s,(i).


Similar with the procedure of UKF (<ref>)-(<ref>), the URTSS gain could be computed as

    K_k^s =C_k,ζ-1,ζ^sP_yy^s^-1,

where the predicted cross covariance C_k,ζ-1,ζ^s and the predicted covariance P_yy^s can be separately derived as

    C_k,ζ-1,ζ^s = ∑_i = 0^2R W_c^(i)[χ^s,(i)_k -θ̅_k(ζ-1) ]( ι_k^s,(i) - y_ζ^ s- )^T,
       P_yy^s = ∑_i = 0^2RW_c^(i)( ι_k^s,(i) - y_ζ^ s- )( ι_k^s,(i) - y_ζ^ s- )^T + Q_ω_k.


Afterwards, the smoothed central DOA and the corresponding error covariance of the central DOA can be acquired as

    θ̅_k^s(ζ)=θ̅_k(ζ)+ K_k^s[θ̅_k^s(ζ+1) - y_k^ s- ],
        P_k^s(ζ)= P_k(ζ)+ K_k^s[P_k^s(ζ+1)-P_yy^s]K_k^s^T.



The most computationally expensive operation of UKF and URTSS is to form the sigma points in Eq. (<ref>), (<ref>) and (<ref>), which needs the Cholesky decomposition of the state covariance matrices. The complexity can be decreased by intaking the square-root UKF and the square-root URTSS, which propagates the Cholesky factors rather than the state covariance matrices to avoid the decomposition.



 §.§  AS Δθ_k(ζ) Tracking

 Let us take the first order Taylor series expansion of 𝐚[θ_kp(ζ)] around central DOA of θ̅_k(ζ)  as
    
    𝐚[θ_kp(ζ)]   =𝐚[θ̅_k(ζ)+θ̃_kp(ζ)]≈𝐚[θ̅_k(ζ)]+θ̃_kp(ζ)∂𝐚[θ̅_k(ζ)]/∂θ̅_k(ζ),

where the high order terms of the series are neglected since the value of θ̃_kp(ζ) is small for massive MIMO systems <cit.>.

The uplink received signal can be rewritten as

    𝐱(n)=∑_k=1^K[𝐚[θ̅_k(ζ)]z_k1(n) + ∂𝐚[θ̅_k(ζ)]/∂θ̅_k(ζ)z_k2(n)]+𝐰(n),

where z_k1(n)=s_k(n)∑_p=1^Pe^ϑ_kpα_kp,
and z_k2(n)=s_k(n)∑_p=1^P e^ϑ_kpα_kpθ̃_kp(n).
Moreover, it can be computed that E[z_k1(n)z_k1^H(n)]=S_kσ_kα^2,
E[z_k2(n)z_k2^H(n)]=S_kσ_k^2σ_kα^2, where σ_kα is the covariance of α_kp and S_k=E|s_k(n)^2|.

Define
    𝐀=[ 𝐚[θ̅_1(ζ)],…, 𝐚[θ̅_K(ζ)], ∂𝐚[θ̅_1(ζ)]/∂θ̅_1(ζ),…,∂𝐚[θ̅_K(ζ)]/∂θ̅_K(ζ)],
       𝐳(n)=[z_11(n),…,z_K1(n),z_12(n),…,z_K2(n) ]^T.

The uplink received signal can be rewritten as

    𝐱(n)=𝐀𝐳(n)+𝐰(n).

The covariance matrix in each block ζ can be computed as

    𝐑_𝐱=E[𝐱(n)𝐱^H(n)]=𝐀Σ𝐀^H+σ_n^2𝐈_M,  n=0,…,N-1,

where Σ= ρdiag( S_1σ_1α^2,…,S_Kσ_Kα^2,S_1σ_1^2σ_1α^2,…,S_1σ_K^2σ_Kα^2).

From (<ref>),  Σ  can be obtained as

    Σ̂=𝐀̂^†(𝐑_𝐱-σ_n^2𝐈_M)(𝐀̂^H)^†,

where 𝐀̂ is the estimation of 𝐀, which can be obtained by replacing the DOA with the tracked θ̅_k.

According to the matrix structure of Σ̂, the variance σ_k^2 can be estimated as

    σ̂_k^2=[Σ̂]_K+k,K+k/[Σ̂]_k,k.


Since θ_kp is uniformly distributed within [θ̅_k(ζ)-Δθ_k(ζ), θ̅_k(ζ)+Δθ_k(ζ)],
we can compute Δθ_k(ζ)=√(3σ̂_k^2) for the current block. Then, according to (<ref>), the SSI set of ℬ_k(ζ) for the uplink channel can be derived as

    ℬ_k(ζ)=[⌊ Md/λsin(θ̅_k(ζ)-Δθ_k(ζ))⌋, ⌈ Md/λsin(θ̅_k(ζ)+Δθ_k(ζ))⌉].



For different deployment circumstances, we can  assume other system equations. For example, we can consider the fixed trajectories for the high-speed railway and unmanned aerial vehicle (UAV) applications. This new information could definitely be utilized to improve the tracking accuracy.




§ TRACKING THE CHANNEL GAINS 

In this section, we will use a small amount of pilot symbols to track the time-varying channel gains for  both uplink and downlink channels with ST-BEM.



 §.§  Tracking Uplink Channel Coefficients

The pilot symbol aided modulation  technique <cit.> is used to derive the time-varying channels, where the pilots are inserted among the information symbols in each interval of NT_s. Define 𝒯_t={n_0,n_1,…,n_T-1}⊂{0,…,N-1}
as the time index set for pilot symbols. Since channels of different users in the same group could be distinguished by their spatial
signatures, we could assign the same pilot sequence to  each user in the same group to save
the orthogonal training resources.
Denote the received training sequences at BS as
𝐘=[𝐲(n_0),𝐲(n_1),…, 𝐲(n_T-1)],
and the common pilot sequences for the g-th group
as 𝐒_g={s_g(n_0),s_g(n_1),…,s_g(n_T-1)} with ∑_i=0^T-1|s_g(n_i)|^2=1. Then, it holds that

    𝐘   =∑_g=1^G∑_k∈𝒰_g[𝐡_k(n_0), …,𝐡_k(n_T-1)]√(P_k^)𝐒_g+𝐍=∑_g=1^G∑_k∈𝒰_g√(P_k^)𝐅^HΓ_k 𝐂𝐒_g +𝐍
       =𝐅^H(∑_k∈𝒰_1√(P_k^)Γ_k,…,∑_k∈𝒰_G√(P_k^)Γ_k)  [(𝐂𝐒_1)^H,…, (𝐂𝐒_G)^H]^H+𝐍
       =𝐅^HΓ[(𝐂𝐒_1)^H,…, (𝐂𝐒_G)^H]^H+𝐍,

where P_k^ is the uplink power constraint of user k, Γ_k=[γ_k,0,γ_k,1,…,γ_k,M-1]^T denotes the CE-BEM coefficients for user-k, 𝐂=[𝐜_n_0,𝐜_n_1,…,𝐜_n_T-1],
 and 𝐍 is the noise matrix whose elements are  i.i.d. 𝒞𝒩(0,σ_n^2).

When T≥ G(μ+1), there will be adequate observations to estimate all unknowns parameters in Γ, and the standard least square (LS) estimator yields <cit.>

    Γ̂ = 𝐅𝐘{[(𝐂𝐒_1)^H,…, (𝐂𝐒_G)^H]^H}^†.

The mean square error (MSE) of Γ̂ can be computed as

    𝔼{Γ̂-Γ_F^2}=Mσ_n^2 {([(𝐂𝐒_1)^H,…, (𝐂𝐒_K)^H]^H[𝐂𝐒_1,…,𝐂𝐒_K])^-1}.


According to <cit.>, we know that
the optimal pilot symbols to minimize the MSE in (<ref>) should satisfy the following constraints.

    𝐂𝐒_g𝐒_g^H𝐂^H=𝐈_μ+1,𝐂𝐒_g𝐒_g'^H𝐂^H=0_μ+1,  ∀ g≠ g'.

Then, the optimal pilot symbols for different groups are proved to be equi-powered, equi-spaced over {0,…,N-1},
and should be phase shift orthogonal  <cit.>. One typical example of this kind of
pilot sequence is

    s_g(n_i)=√(1/T)e^j2π i (g-1)(μ+1)/T, g=1,…,G,  i=0,…,T-1.


Next, we focus on users in group g and derive that

    Γ̂_k = Γ_k+∑_l∈{𝒰_g\ k}√(P_l^/P_k^)Γ_l+
      1/√(P_k^/σ_n^2)𝐍_k,

where 𝐍_k∈ℂ^M×μ+1 has the i.i.d. 𝒞𝒩(0,1) elements.

Considering the distinct spatial signatures of users in the same group, we can  extract

    [Γ_k]_ℬ_k(ζ),:=[Γ̂_k ]_ℬ_k(ζ),: =[Γ_k]_ℬ_k(ζ),:
    +∑_l∈{𝒰_g\ k}√(P_l^/P_k^)[Γ_l]_ℬ_k(ζ),:
    +[𝐍_k]_ℬ_k(ζ),:/√(P_k^/σ_n^2),

where the second term ∑_l∈{𝒰_g\ k}√(P_l^/P_k^)[Γ_l]_ℬ_k(ζ),: is the pilot contamination term caused by reusing the same pilot in one group. Since ℬ_l(ζ) and ℬ_k(ζ) are kept away from each other, the entries of [Γ_l]_ℬ_k(ζ),: in (<ref>) are very small and negligible. Therefore, Γ_k can be approximated as

    Γ_k=[0^T [Γ_k]_ℬ_k(ζ),:^H 0^T]^H.




While T=K(μ+1) pilots are needed for traditional channel estimation scheme <cit.>, the total pilot overheads of uplink channel tracking is significantly reduced to T=G(μ+1) (G≪ K) through the proposed method, which greatly decreases the uplink training overhead.




 §.§  Tracking Downlink Channel Coefficients
 
Denote the downlink channel from BS to user-k as 𝐠_k(n)∈ℂ^M× 1. Similar to (<ref>), 𝐠_k(n)∈ℂ^M× 1 can be modeled as
 
    𝐠_k(n)≈∑_q∈ℬ_kg̃_k,q(n)𝐟_q= ∑_q∈ℬ_k(ζ)γ_k,q^d^T𝐜_n 𝐟_q,n=0,…,N-1,

where γ_k,q^d denotes the downlink channel CE-BEM coefficients and
ℬ_k(ζ) is the downlink spatial signatures. Consistent with the uplink channel tracking, the downlink channel tracking can also be simplified to tracking the downlink SSI set ℬ_k(ζ) and estimating those remaining unknown coefficients γ_k,q^d.

Traditionally, in the mode of TDD, the downlink channel state information can be directly derived from channel reciprocity. However, in the
mode of FDD where the channel reciprocity would not hold, the training overhead for the high
dimension channel would be prohibitive, which is the main challenge of massive MIMO
system. Interestingly, electromagnetic characteristics do not change in several dozens of GHz
<cit.>, and the propagation path of electromagnetic wave is reciprocal, i.e., uplink DOA and downlink DOD are reciprocal. Therefore, we could utilize the angle reciprocity
to simplify the complexity of downlink channel tracking in FDD mode.


With angle reciprocity, the downlink spatial signature, denoted as  ℬ_k (ζ), can be directly obtained from the tracked uplink SSI set ℬ_k(ζ). With the angle reciprocity, there is

    sinθ_kp(ζ) = qλ_1/Md=q'λ_2/Md,  q∈ℬ_k(ζ),q'∈ℬ_k(ζ),

where λ_1 and λ_2 denote the uplink and downlink carrier wavelength respectively.
Then the downlink SSI set ℬ_k(ζ) can
be expressed as

    ℬ_k(ζ)=[⌊λ_1/λ_2q_min⌋,⌈λ_1/λ_2q_max⌉],

where q_min≤ q≤ q_max, ∀ q∈ℬ_k(ζ).

Denote the size of ℬ_k(ζ) as τ, i.e., the effective dimension of downlink channel.  According to (<ref>), T=τ(μ+1) pilot symbols are required for downlink channel parameter estimation. Adopting the uplink user grouping strategy directly, the received 𝐲^d_k=[y_k(n_0),…,y_k(n_T-1)]^T at user-k in 𝒰_g can be expressed as

    𝐲^d_k=   ∑_l∈𝒰_g∑_1≤ i≤τ, q_i∈ℬ_l(ζ)√(P_l^/τ)𝐒^d_i𝐂^Tγ_k,q_i^d𝐟_q_i+𝐧_k
    
      =   ∑_l∈𝒰_g[√(P_1^/τ)𝐒^d_1𝐂^T,…,√(P_τ^/τ)𝐒^d_τ𝐂^T]
      ([Γ_k^d]_ℬ_l(ζ)^T)
      𝐅_ℬ_l(ζ),:+𝐧_k,

where Γ_k^d is the downlink channel parameters, 𝐧_k∼𝒞𝒩(0,σ_n^2𝐈_T) is the noise vector, P_l^ is the total downlink power constraint for user l and diagonal matrices 𝐒^d_i∈ℂ^T× T,i=1,…,τ denote the transmitted pilot sequences for users in 𝒰_g.



Similar to (<ref>), the optimal pilot sequences 𝐒^d_i∈ℂ^T× T,i=1,…,τ
for downlink channel tracking are also equi-powered, equi-spaced and
phase shift orthogonal. Then, the downlink channels can be recovered by

    [Γ_k^d]_ℬ_k(ζ)=[Γ_k^d]_ℬ_k(ζ)
    +∑_l∈{𝒰_g\ k}√(P_l^/P_k^)[Γ_l^d]_ℬ_k(ζ)
    +1/√(P_k^/σ_n^2)𝐧_k^',

where 𝐧_k^'∈ℂ^τ× 1 has the i.i.d. 𝒞𝒩(0,1) elements.

It can be found that user k does not need the knowledge of spatial SSI set ℬ_k(ζ) to perform the estimation of [Γ_k^d]_ℬ_k(ζ). Each user only needs to feed back τ components [Γ_k^d]_ℬ_l(ζ)
to BS to fulfill the downlink channel tracking, which removes the necessity of feedback from BS to the user. This is a key advantage that makes the proposed downlink channel tracking strategy suitable for fast-fading environments. The overall channel tracking mechanism is illustrated  in Fig. <ref>.



§ SIMULATIONS 





 In this section, numerical results are presented to demonstrate the effectiveness of the proposed methods. We set M=128,  d=λ/2 and  K=12. Besides, all K users are gathered into 4 spatially distributed clusters, and users can be divided into G=3 groups according to the DOA distribution. The channel vector between the users and the BS is generated according to (<ref>) with AS Δθ=2^∘.
 The performance metric of the channel estimation is taken as the normalized MSE, i.e.

    MSE≜1/N∑_n=1^N∑_k=1^K𝐡_k(n)-𝐡̂_k(n)^2/𝐡_k(n)^2.



In the first example, the performance of the central DOA tracking is shown in Fig. <ref>, where the real trajectory, central DOA from DFT searching, and central DOA tracking without EM leaning are also displayed for comparison. Meanwhile, the DOA tracking performance as a function of SNR is also displayed in Fig. <ref>. The users moves with maximum speed 80 Km/h, and the variation of DOA is generated according to the system equation (<ref>). It can be found from Fig. <ref> that the angle variation is small for massive MIMO systems, and both the DFT searching method and the central DOA tracking without EM learning are approximately consistent with the trajectory but have larger error than the DOA tracking with EM learning. The reason is that the one-order Markov process can well describe the movements of the users, and EM learning can precisely derive the unknown parameters from the system equation and measurement equation. The MSE performance of central DOA tracking with different SNR in Fig. <ref> further verifies the aforementioned analysis. Besides, we can also find that there are error floors for all the displayed method, which arises from the observation error.



Fig. <ref> displays the performances of AS tracking from DFT searching and Taylor series expansion, respectively. Meanwhile, the reference curve, which represents ℬ_k(ζ) that contains 98% power of the channel, is also displayed as comparison. It is seen that the result from DFT searching is not desirable and yields large fluctuation errors, while the one from Taylor series expansion method is much better and closer to the reference curve. The reason is that when the number of the antennas M is not infinite, the power leakage effect and the measurement noise will affect the precision of SSI set tracking. 



Fig. <ref> displays the uplink channel tracking performances as a function of SNR and MSE, where the aging channel without spatial information tracking and the channel with tracked ℬ_k(ζ) are also displayed as comparison. Besides, since the communications standards normally regulate a constant number of system parameters rather than a dynamic quantity for practical deployments, the performances of the proposed method with a predefined size of SSI set, denoted as Υ, are also displayed as comparison. It can be seen that as the SNR increases, there are error floors for the ST-BEM based methods. This phenomenon arises from the truncation error of BEM from the real channel and can also be observed in <cit.>. Besides, we can also see that the performance of channel tracking behaves better with the increase of Υ, and the MSE of the proposed method with the tracked ℬ_k(ζ) is the best in the displayed curves.



Fig. <ref> compares the downlink MSE performance of the proposed method, where κ =M(μ+1) pilot symbols are utilized for the conventional LS method, while  T=κ/8, κ/4, κ/2 are allocated for the ST-BEM method separately. The power constraints are the same ∑_k=1^KP^_k=KTρ for both methods for fairness. It can be found that the proposed method is much superior to the conventional LS method, not only for the higher estimation accuracy in low SNR regions, but also due to the less training overheads, which is consistent with the analysis in the Section <ref>.



Lastly, the bit error rate (BER) performance under QPSK modulation is shown for the downlink data transmission in Fig. <ref>. The results  of perfect CSI, CSI from the proposed method, CSI from the SBEM method <cit.>, and CSI from the conventional LS method are displayed under the matched filter precoding scheme. The overall training power is set as the same for each method to keep the comparison fairness. It is easy to find that the method in <cit.> is not efficient for highly time-varying channel. The proposed channel tracking method performs better than the LS method at the low SNR region, and the reasons can be found from (<ref>) where the proposed method only involves τ components of the noise vector, while the conventional LS method includes the whole noise power.
Moreover, it can be seen that the BER achieved by ST-BEM has about 0.5 dB gap from that of perfect CSI, which corroborates the effectiveness of the proposed ST-BEM channel tracking method.



§ CONCLUSION

In this paper, the problem of time-varying channel tracking was considered for massive MIMO systems. We proposed an ST-BEM to reduce the effective dimensions of the channels, where the spatial channel is decomposed into the time-varying spatial information and the time-varying gain information. Since the spatial information is characterized by the central DOA and AS of the incident signal, we use the one-order Markov process to describe the user's movements and UKF to blindly track the central DOAs, where the unknown system parameters in the system equation and measurement equation are estimated by EM learning. In addition, AS is also blindly tracked through Taylor series expansion of the steering vector. With the derived SSI set, the time-varying gain information of ST-BEM can be derived through a few pilot symbols. Simulation results show that the proposed method offers a practical and effective way for massive MIMO channel tracking.






1.5
99

ICC
J. Zhao, F. Gao, W. Jia, J. Zhao, and W. Zhang, “Channel tracking for massive MIMO systems with spatial-temporal basis expansion model,” in Proc. of IEEE ICC'17, Pairs, France, July 2017.

marzetta2010noncooperative
T. L. Marzetta, “Noncooperative cellular wireless with unlimited numbers of
  base station antennas,” IEEE Trans. Wireless Commun., vol. 9, no. 11,
  pp. 3590–3600, Nov. 2010.

chen1
S. Chen, and J. Zhao, “The requirements, challenges, and technologies for 5G of terrestrial mobile telecommunication,”  IEEE Commun. Mag., vol. 52, no. 5, pp. 36–43, May  2014.
chen2
S. Chen, F. Qin, B. Hu, X. Li, and Z. Chen, “User-centric ultra-dense networks for 5G: challenges, methodologies, and directions,”  IEEE Wireless Commun., vol. 23, no. 2, pp. 78–85, May  2016.

five
F. Boccardi, R. W. Heath, A. Lozano, T. L. Marzetta, and P. Popovski, “Five
  disruptive technology directions for 5G,”  IEEE Commun. Mag., vol. 52, no. 2, pp. 74–80, Feb. 2014.

han(9)
S. Hur, T. Kim, D. J. Love, J. V. Krogmeier, T. A. Thomas, and A. Ghosh, “Millimeter wave beamforming for wireless backhaul and access in small cell networks,” IEEE Trans. Commun., vol. 61, no. 10, pp. 4391–4403, Oct. 2013.


scalingupMIMO
F. Rusek, D. Persson, B. K. Lau, E. G. Larsson, T. L. Marzetta, O. Edfors, and
  F. Tufvesson, “Scaling up MIMO: Opportunities and challenges with very large
  arrays,” IEEE Signal Process. Mag., vol. 30, no. 1, pp. 40–60, Jan. 2013.




Caire
A. Adhikary, J. Nam, J.-Y. Ahn, and G. Caire, “Joint spatial division and
  multiplexing–the large-scale array regime,” IEEE Trans. Inf. Theory,
  vol. 59, no. 10, pp. 6441–6463, Oct. 2013.

KLT
X. Rao, and V. K. Lau, “Distributed compressive CSIT estimation and feedback
  for FDD multi-user massive MIMO systems,” IEEE Trans. Signal
  Process., vol. 62, no. 12, pp. 3261–3271, June. 2014.

CCS
W. U. Bajwa, J. Haupt, A. M. Sayeed, and R. Nowak, “Compressed channel sensing:
A new approach to esitmating sparse multipath channels,” Proc. IEEE, vol. 98, no. 6, pp. 1058–1076, June 2010.

yin
H. Yin, D. Gesbert, M. Filippou, and Y. Liu, “A coordinated approach to
  channel estimation in large-scale multiple-antenna systems,” IEEE J.
  Sel. Areas Commun., vol. 31, no. 2, pp. 264–273, Feb. 2013.


xieCR
H. Xie, B. Wang, F. Gao, and S. Jin, “A full-space spectrum-sharing strategy for massive MIMO cognitive radio,”  IEEE J. Select. Areas Commun., vol. 34, no. 10, pp. 2537–2549, Oct. 2016.

xie
H. Xie, F. Gao, S. Zhang, and S. Jin, “An unified transmission strategy for TDD/FDD massive MIMO systems with spatial basis expansion model,”  IEEE Trans. Veh. Technol.,  vol. 66, no. 4, pp. 3170–3184, Apr. 2017.

overview
H. Xie, F. Gao, and S. Jin, “An overview of low-rank channel estimation for massive MIMO systems,”  IEEE Access, vol. 4, pp. 7313–7321, Nov.  2016.

fandian
D. Fan, F. Gao, G. Wang, Z. Zhong, and A. Nallanathan, “Angle domain signal processing aided channel estimation for indoor 60GHz TDD/FDD massive MIMO systems,”  IEEE J. Select. Areas Commun., vol. 35, no. 9, pp. 1948–1961, Sep. 2017.

linhai
H. Lin, F. Gao, S. Jin, and G. Ye Li, “A new view of multi-user hybrid massive MIMO: non-orthogonal angle division multiple access,”  IEEE J. Select. Areas Commun., vol. 35, no. 10, pp. 2268–2280, Oct. 2017.

zhao
J. Zhao, F. Gao, W. Jia, S. Zhang, S. Jin, and H. Lin, “Angle domain hybrid precoding and channel tracking for mmWave massive MIMO systems,” IEEE Trans. Wireless Commun., vol. 16, no. 10, pp. 6868–6880, Oct. 2017.

channelaging1
K.T. Truong, and R.W. Jr, “Effects of channel aging in massive MIMO systems,” J. Commun. networks, vol. 50, no. 12, pp. 338–351, Dec. 2013.

channelaging2
A. Papazafeiropoulos, “Impact of general channel aging conditions on the downlink performance of massive MIMO,” IEEE Trans. Veh. Technol., vol. 50, no. 12, pp. 338-351, Dec. 2016.

channelaging3
C. Kong, “Sum-rate and power scaling of massive MIMO systems with channel aging,” IEEE Trans. Commun, vol. 63, no. 12, pp. 4879-4893, Dec. 2015.

gao_relay14
G. Arredondo, W. H. Chriss, and E. H. Walker, “A multipath fading
simulator for mobile radio,” IEEE Trans. Commun., vol. 21, no. 11,
pp. 1325–1328, Jan. 2003.

globalcom
H. Xie, F. Gao, S. Zhang, and S. Jin, “Spatial-temporal BEM and channel estimation strategy for massive MIMO time-varying systems,” in Proc. of IEEE GLOBECOM'16, Washington DC, USA, Dec. 2016.

GiannakisBEM
G. B. Giannakis, and C. Tepedelenlioglu, “Basis expansion models and diversity
  techniques for blind identification and equalization of time-varying
  channels,” Proc. IEEE, vol. 86, no. 10, pp. 1969–1986, Oct. 1998.

jinshi
C. Sun, X. Gao, S. Jin, M. Matthaiou, Z. Ding, and C. Xiao, “Beam division
  multiple access transmission for massive MIMO communications,”
  IEEE Trans. Commun., vol. 63, no. 6, pp. 2170–2184, June 2015.

2D-ESPRIT
A. Hu, T. Lv, and H. Gao, “An ESPRIT-based approach for 2-D localization of incoherently distributed sources in massive MIMO systems,” IEEE J. Sel. Topics Signal Process., vol. 8, no. 5, pp. 996–1011, Oct. 2014.

beamspace1
J. Brady, N. Behdad, and A. M. Sayeed, “Beamspace MIMO
for millimeter-wave communications: System architecture, modeling,
analysis, and measurements,”  IEEE Trans. Antennas Propag., vol. 61, no. 7,
  pp. 3814–3827, July 2013.
Robert
A. Alkhateeb, O. Ayach, G. Leus, and R. Heath, “Channel estimation
and hybrid precoding for millimeter wave cellular systems,” IEEE J. Sel.
Areas Signal Process., vol. 8, no. 5,
  pp. 831–846, Oct. 2014.









chengshanxiao
Y. R. Zheng, and C. Xiao, “Simulation models with correct statistical
properties for rayleigh fading channels,” IEEE Trans. Commun., vol. 51, no. 6, pp. 920–928, June 2003.


Clarke
R. H. Clarke, “A statistical theory of mobile-radio reception,” Bell System Technical Journal,
vol. 47, no. 6, pp. 957–1000, July 1968.

gao_relay
G. Wang, F. Gao, W. Chen, and C. Tellambura, “
Channel estimation and training design for two-way relay networks in time-selective fading environments,”
IEEE Trans. Wireless Commun., vol. 10, no. 8, pp. 2687–2691, Aug. 2011.


temporalBEM1
G. B. Giannakis, and C. Tepedelenlioglu, “Basis expansion models and diversity
  techniques for blind identification and equalization of time-varying
  channels,” Proc. IEEE, vol. 86, no. 10, pp. 1969–1986, Oct. 1998.

data-aided
J. Ma, and L. Ping, “Data-aided channel estimation in large antenna systems,”
in Proc. of IEEE ICC'14, Sydney, NSW, pp. 4626–4631, June 2014.


mmdoa
C. Zhang, D. Guo, and P. Fan, “Tracking angles of departure and arrival in a mobile millimeter wave channel,” in Proc. of IEEE ICC'16,  Kuala Lumpur, pp. 1–6, May 2016.

 liulingjia0
A. Wang, L. Liu, and J. Zhang, “Low complexity direction of arprival (DoA) esitmation for 2D massive MIMO systems,” in IEEE GLOBECOM Workshop on Emerging Technologies for LTE-Advanced
and Beyond-4G, Anaheim, CA, pp. 703–707, Dec. 2012.




 EM
S. Zhang, F. Gao, J. Li, and H. Li, “Time varying channel estimation for DSTC-based
relay networks: tracking, smoothing and BCRBs,” IEEE Trans. Wireless Commun., vol. 14, no. 9, pp. 5022–5037, Sep. 2015.

 UKF
S. J. Julier, and J. K. Uhlmann, “Unscented filtering and nonlinear estimation,”
Proc. IEEE, vol. 92, no. 3, pp. 401–422, Mar. 2004.

 URTSS
S. Zhang, F. Gao, and H. Li, “Time varying individual channel estimation for one-way relay networks with UKF and URTSS,” in Proc of IEEE Globecom'14, pp. 3567–3572, Dec. 2014.

PSAM
J. K. Cavers, “An analysis of pilot symbol assisted modulation for
Rayleigh fading channels,” IEEE Trans. Veh. Technol., vol. 40, no. 4, pp. 686–693, Nov. 1991.








DOAreciprocity1
A. Kuchar, M.l Tangemann, and E. Bonek, “A real-time DOA-based smart antenna processor,”
IEEE Trans. Veh. Technol., vol. 51, no. 6, pp. 1279–1293, Nov. 2002.

DOAreciprocity2
Y. Han, J. Ni, and G. Du, “The potential approaches to achieve channel reciprocity in FDD system with frequency correction algorithms,” 5th International ICST Conference on Communications and Networking in China (CHINACOM), Beijing, China, pp. 1–5, Aug. 2010.

DOAreciprocity3
K. Hugl, K. Kalliola, and J. Laurila, “Spatial reciprocity of uplink and downlink radio channels in FDD systems,” in COST 273 TD(02)066, 2002.



