Estimation of Time-Varying Channels in Virtual Angular Domain for Massive MIMO Systems

This paper is concerned with the problem of channel estimation in virtual angular domain in massive MIMO orthogonal frequency division multiplexing (OFDM) systems, in the presence of time varying channels. The proposed channel estimator is based on dichotomous coordinate decent joint sparse recovery (DCD-JSR) algorithm, which was well suited for massive MIMO systems with channels static over several OFDM symbols, but it is not accurate for time-varying channels. The estimator proposed in this paper exploits the basis expansion model (BEM) for approximation of the time-variation of the channel and the sparsity of the channel in the virtual angular domain with a common support over OFDM subcarriers and OFDM symbols for joint estimation of the BEM expansion coefficients. It is shown that, when using the Legendre polynomials as the BEM and depending on the normalized Doppler frequency, only a small number of expansion coefficients is required to provide accurate channel estimation.


I. INTRODUCTION
Massive multiple-input multiple-output (MIMO) has been proposed for the next generation communication systems, since it can provide high data rates by employing the wireless transmitter with a high number of antennas and provide high spectral efficiency by utilising the additional degree of freedom in the spatial domain [1], [2]. The massive MIMO faces challenges such as practical system modeling, channel estimation and so on [3]. Especially, accurate channel state information (CSI) serves as the key information for techniques such as coherent detection, therefore, an accurate channel estimation is necessary for practical applications of massive MIMO systems [4].
In [5], experimental results have shown that massive MIMO channels exhibit sparsity, due to the relatively small angle spread perceived from a base station (BS) between user and BS. According to [6], [7], and [8], the common sparsity is shared by different subcarriers due to the spatial The associate editor coordinating the review of this manuscript and approving it for publication was Wei Feng .
propagation property of the wireless channel, such as the number of scatterers is remaining nearly unchanged over the system bandwidth, which is referred to as the spatially common sparsity over multiple subcarriers. Besides, since the angle variation from the user to the BS is relatively slow, we can consider the channel support is static over a sequence of OFDM symbols [7], [8]. Thus, in this paper, we consider that the channel exhibits sparsity, and shares common support over several OFDM symbols.
Sparse recovery techniques are attractive for channel estimation [9], [10], [11]. The channel estimation deals with estimation of complex-valued parameters. In [12], the dichotomous coordinate descent (DCD) algorithm is proposed for solving complex-valued sparse problems. It has been indicated in [12] and [13] that the computational complexity of the algorithm is dominated by the computational complexity of a small number of successful iterations, while most of the operations of the DCD algorithm are additions and bit-shifts. The algorithm provides accurate channel estimation and is computationally efficient [14], however, it does not suit massive MIMO channel estimation in the VOLUME 11, 2023 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ virtual angular domain. In [15], by employing the DCD algorithm and considering that the channel is statistic over several OFDM symbols, the DCD-JSR algorithm is proposed for the channel estimation in massive MIMO in the virtual angular domain; the simulation results have shown that the DCD-JSR algorithm outperformes the distributed sparsity adaptive matching pursuit (DSAMP) algorithm [7], however, it cannot provide accurate channel estimation when the channel is time-varying. In [16] and [17] and many other publications, it has been shown that, the time-varying channel can be approximated accurately by employing the basis expansion model (BEM). Consequently, estimation of a realization of random process describing the time-varying channel is transformed into estimation of time-invariant expansion coefficients [18]. In [19], the Karhunen-Loeve BEM has been proposed for estimating the time-varying channel, however, it is very sensitive to the variations of the channel statistics. In [20] and [21], algebraic polynomial BEMs have been employed to estimate the timevarying channel, where the channel vectors can be approximated as a linear combination of a set of polynomials. In [17], the experimental results have indicated that, by employing the Legendre polynomials as the BEM, the channel variation could be accurately approximated. Thus, in this paper, we consider employing the Legendre polynomials to approximate the time-variation of the channel in the virtual angular domain.
In this paper, by combining the DCD-JSR algorithm [15] and the BEM, we show that the modified DCD-JSR algorithm can estimate the channel in OFDM system operating over frequency selective and highly mobile wireless timevarying channels. Simulation results show that, compared to the original DCD-JSR algorithm, the modified DCD-JSR algorithm could provide better MSE performance when estimating time-varying channels.
The paper is organized as follows. Section II describes the system model. In Section III, channel model for time varying channel is introduced. In Section IV, the processes of employing the basis functions, and estimating the time-varying channel are described. In Section V, numerical examples are analysed and, finally, Section VI presents the conclusion.
In this paper, capital and small bold fonts are used to denote matrices and vectors, respectively, (x) n denotes the nth element of the vector x, R q denotes the qth column of the matrix R, and R n denotes the nth row of the matrix R, R m,n denotes an element of the matrix R. The transpose operator is given by (.) T , (.) * denotes the conjugate operator, (.) † denotes the Moore-Penrose inversion, and (.) H denotes the Hermitian transpose operator. The 0 -norm and 2 -norm are represented by ||. 0 and ||. 2 , respectively. We use I to denote a support, |I| is the cardinality of the support I, I c is the complement of I, R I is a matrix obtained from R, and which only contains rows corresponding to support I. R I,I is an |I| × |I| matrix obtained from R by collecting elements from columns and rows corresponding to I, and x I is the subset of x that includes non-zero elements from x corresponding to I. We use h to denote a channel vector andh to denote the channel vector in the virtual angular domain,h n denotes the channel vector corresponding to the nth subcarrier. R denotes the real part of a complex number, and j = √ −1.

II. SYSTEM MODEL A. CHANNEL MODEL
In a massive MIMO system, consider a time interval consisting of J OFDM symbols. M antennas are employed at the BS to serve K single-antenna users simultaneously, where M k. At the tth OFDM symbol, 1 ≤ t ≤ J, for the nth subcarrier, 1 ≤ n ≤ N, the received signal model for the kth user, 1 ≤ k ≤ K, in the frequency domain is given by where h t k,n ∈ C M×1 represents the downlink channel between the M BS antennas and the kth user, x t n ∈ C M×1 is the vector of transmitted symbols (data or pilot symbols) and w t k,n is the corresponding additive white Gaussian noise (AWGN). For a single user, we can drop the index k, thus we can write Matrix A B is used to modify the channel vector h t n into a channel vectorh t n in the virtual angular domain, and it is determined by the geometric structure of the antenna array. We consider a uniform linear array with the antenna spacing λ/2, where λ is the wavelength, then A B becomes the discrete Fourier transform (DFT) matrix. Thus we obtain where h t n T = h t n T A * B . The channel vector in the angular domain divides the covering area of the BS into angular intervals, and the mth element ofh t n corresponds to the mth virtual angle, where 1 ≤ m ≤ M.
In massive MIMO systems, the BS is often positioned at a high elevation with a small number of scatterers (in contrast to the number of antennas) [5], [22], but the scatterers on the user side are relatively rich. Thus, we can consider for the kth user, in the virtual angular domain, only a small number |I| of multipath arrivals contains the majority of the channel energy. Hence, we have |I| ≤ M, and the channel vector exhibits sparsity in the virtual angular domain.
We consider that the channel is static over one OFDM symbol, furthermore, as indicated in [7], [8], and [15], since the spatial propagation characteristics such as scatterers are almost unchanged over the system bandwidth, as shown in Fig.1, we can consider that in the same OFDM symbol, for different subcarriers, the subchannel exhibits common sparsity.
According to [23], in time-varying scenarios, the variation of the arrival angles is usually much slower than that of channel gain, thus, we can consider that during J successive OFDM symbols, the channel exhibits common sparsity. Thus, we can obtain where I j n is the support set for the Jth OFDM symbol at the nth subcarrier.
In this paper, we consider the pilot-aided channel estimation. For the tth OFDM symbol, a part of subcarriers is used for transmitting pilot symbols s t p ∈ C M×1 , and the received signal at the pilot subcarrier n(p) is given by: where θ t,m,p are independent random numbers uniformly distributed in (0, 2π].

B. CHANNEL ESTIMATION APPROACH
The conventional method to acquire the channel state information (CSI) in frequency-division-duplexing (FDD) systems is as follows: the BS transmits downlink pilot symbols to a user, so the user can estimate the downlink CSI locally and then feed it back to the BS via an uplink channel [24]. If we are employing conventional CSI estimation techniques (such as the minimum mean square error (MMSE) estimator), since the number of pilot symbols required at the BS has to scale linearly with the number of transmit antennas [22], it would cause prohibitively large overhead for both pilot training (downlink) and CSI feedback (uplink). Hence, to solve the overhead issue, as suggested in [7], the channel estimation is performed at the BS. The channel estimation scheme is summarized as follows.
1 . In each OFDM symbol, every BS antenna broadcasts pilot symbols to users, the kth user receives the signal y k ∈ C M×1 and feeds it back to the BS. The BS recovers the CSI for each user based on the feedback signals y k , k = 1, . . . , K. As shown in Fig. 2, each OFDM symbol contains N subcarriers, while P subcarriers are used to transmit pilot symbols. The user feeds back the received signal to the BS without performing downlink channel estimation. 2 . At the BS, by employing the BEM to approximate the time variation of the channel, the DCD-JSR algorithm can jointly estimate the common support I for multiple sparse virtual angular domain channels. The least squares (LS) algorithm [25] is then employed to acquire the CSI based on an estimate of the common support I.

C. RECEIVED SIGNAL
For the tth OFDM symbol, at the pth pilot subcarrier, the signal received at the BS is given by: Here, φ t p = s t p T A * B T ∈ C 1×M is the sensing vector defined by the DFT matrix and the pilot symbols.h t n(p) ∈ C M×1 is the sparse channel vector for the n(p)th subcarrier, and v t p is the corresponding noise, which contains both downlink and uplink channel noise.

III. TIME VARYING CHANNEL
In practice, due to the user mobility, the propagation of wireless signals would face the time-varying environment [1]. Due to the simple implementation and the orthogonality between columns of the basis-expansion matrix, using the Legendre polynomials as basis functions on the investigated time interval has been considered in the literature [17]. In this paper, we consider that for the pth pilot subcarrier, the time-varying channel vector can be approximated by N b basis functions:ĥ where b i (t) is the tth element of a vector b i ∈ C J×1 representing samples of the basis function b i (t), c i,p ∈ C M×1 are expansion coefficients for the ith basis function at the pth pilot VOLUME 11, 2023 subcarrier. This approximation is valid if the channel variation are slow enough over one OFDM symbol, so that we can ignore the intercarrier interference between the OFDM subcarriers. Thus, we mostly deal with variation of the channel over a sequence of J OFDM symbols. In such scenarios, the channel approximation in (8) can be made arbitrary accurate by choosing a large enough number of basis functions N b . In our numerical investigation below, we consider scenarios with a limited product f d T of the Doppler frequency f d by the (orthogonality) duration of the OFDM symbol T, specifically, f d T ≤ 0.05. In this paper, we consider employing the Legendre polynomials as the basis functions [26], defined as: For the tth OFDM symbol at the pth pilot subcarrier, by substituting (8) into (7), we can obtain: Sinceh t n(p) ∈ C M×1 exhibit common sparsity, the expansion coefficient vectors c i,p ∈ C M×1 also exhibit common sparsity. Thus, the task of estimating JM channel coefficients is transformed into estimating only N b |I| expansion coefficients with usually N b J and |I| ≤ M. We collect the received signal samples r t p , 1 ≤ t ≤ J, in a vector r p = r 1 p , r 2 p , . . . , r J p T ∈ C J×1 , then we have: . . , F t i,p ] ∈ C J×M is a matrix whose tth row is F t i,p , and v p = v 1 p , v 2 p , . . . , v J p T ∈ C J×1 is the noise vector. Since the expansion coefficients exhibit common sparsity, we can firstly estimate the common support and then find the expansion coefficients.

IV. DCD-JSR ALGORITHM FOR TIME-VARYING CHANNEL
Here, the homotopy DCD algorithm [12] is used to estimate the support of the expansion coefficients, as shown in Table 1. First, we apply the homotopy DCD algorithm to the 2 0 optimization problem of minimizing: Here, we solve the optimization problem for the pth pilot subcarrier of the ith expansion coefficient, and τ ∈ [0, 1) is a regularization parameter. The second term in (12) makes it a non-convex problem and the solution of it is NP-hard.
To solve the problem, we initially assign the support set I p = ∅, and by following the proposition in [12] we can add new elements into the support or remove elements from  the support in several iterations, thus, the estimated expansion coefficientsc i,p can be obtained. Therefore we need to assign initially a high value to the regularization parameter τ = τ max , so that the second term in (12) dominates the cost function to provide an empty support I p = ∅. In the homotopy iterations, by gradually reducing value of τ as τ ← γ τ , where γ ∈ [0, 1), new elements can be added to the support or removed from the support [12]. The algorithm stops when τ < τ min , where τ min = µ τ τ max and µ τ ∈ [0, 1) is a predefined parameter.
To reduce the computational complexity [12], instead of solving the LS problem in Table 1, step 3, we employ the DCD iterations to solve the LS problem, as shown in Table 2, where N u is the number of successful DCD iterations, and a successful DCD iteration means that the solution is updated.
Following is an example of how we estimate the common support for the expansion coefficients. For the  simulation scenario, we consider a massive MIMO system, SNR=20 dB, M = 128, P = 64, J = 100, the normalized Doppler frequency f d T = 0.05 and N b = 3, |I| = 8. Fig. 3 shows the normalized energy of the elements of the estimated expansion coefficient c 1,1 , for the first basis function (zero-order Legendre polynomial) at the first pilot subcarrier in the angular intervals. It is easy to see that, due to the large variance of the energy of the elements in the angular intervals, we can not clearly identify the support.
Since all expansion coefficients share a common support, for the pth pilot subcarrier, we can compute another vector with contribution from all the expansion coefficients: Here, as shown in Fig. 4,q p is a vector that contains normalized energy of elements for all expansion coefficients at the pth pilot subcarrier in the angular intervals. The new plot shows clearly 4 of 8 non-zero directions. However, the variance is still large and we can not estimate reliably the support at this step. As indicated in the previous section, since the expansion coefficients c i,p share the common support among P subcarriers, we can compute a new vector, which takes this into account: q = q/(max(q)N b P).
As shown in Fig.5, here,q ∈ C M×1 is a sparse vector with elements averaging contribution from all pilot subcarriers and all expansion coefficients. We can acquire now the common supportĨ by using the hard thresholding where ξ is a predefined thresholding parameter. Based on the support estimateĨ, the MMSE approach [25] is employed to estimate the expansion coefficientsc p : Here, w is the identity matrix with size of |Ĩ| × |Ĩ|,c p ∈ C 1×MN b is a vector containingc i,p for N b basis functions, and F p ∈ C J×MN b is the matrix containing N b vectors F i,p , and σ 2 is the noise variance, which is assumed to be known. The estimation of noise variance is a well known problem (e.g., see [27], [28], [29]) and it is not considered here.

V. SIMULATION RESULTS
The mean squared error (MSE) of the channel estimation will be used to asses the algorithm performance. First, we compute the MSE for the tth OFDM symbol at the pth subcarrier: ||h t n(p) || 2 2 = (h t n(p) ) H (h t n(p) ), whereĥ t p is the estimated channel vector obtained from (8), andh t p is the true channel vector. Then the overall MSE for the channel estimation is computed by: The MSE in (22) is further averaged over the simulation trials. We consider simulation scenarios corresponding to a MIMO system with a uniform linear array. For the massive MIMO system, in most simulation scenarios, we consider the number of antennas M = 128, the number of pilot subcarriers P = 64, the sampling frequency f s = 15.36 MHz, the time interval for one OFDM symbol T = 66.7 µs, the carrier frequency f c = 2.5 GHz and the number of simulation trials N s = 500. The performance of the oracle LS algorithm [27] with known support is adopted as the performance bound. In most scenarios, we consider two cases, SNR = 10 dB and SNR = 20 dB, and the user mobility with v = 120 km h and v = 300 km h . The Doppler frequency f d can be obtained by using: where v c = 3 × 10 8 m/s is the light speed.
The channel estimation performance is investigated in several ways. First, we investigate the MSE performance of the proposed algorithm with different number of employed OFDM symbols, then we compare the MSE performance against the number of basis functions, with normalized Doppler frequencies f d T = 0.02, f d T = 0.05, where the Doppler frequency is approximately f d = 300 Hz, f d = 700 Hz, respectively. After that, we compare the MSE performance for different SNR, considering the normalized Doppler frequency f d T = 0.05. The MSE performance against the number of DCD-iterations is investigated to show the convergence of the proposed algorithm, and the MSE performance against the number of employed antennas is also investigated. Furthermore, we compare the MSE performance of the DCD-JSR algorithm and the distributed sparsity adaptive matching pursuit (DSAMP) algorithm from [7] against the number of non-zero virtual angles, the oracle LS algorithm with known support is adopted as the performance bound [25]. At last, the computational complexity of the DCD-JSR algorithm is analyzed.
In Fig. 6, we show the MSE performance of the DCD-JSR algorithm for different number of employed OFDM symbols, for SNR = 20 dB, |I| = 3, f d T = 0.02 and f d T = 0.05. It can be seen that, as the number of OFDM symbols increases, the better MSE performance is provided by employing more basis functions. Thus we can conclude that, with longer data packets, more basis functions are required for the DCD-JSR algorithm to provide the best MSE performance, and the minimum MSE is also reduced.
In Fig. 7, we compare the MSE performance for different normalized Doppler frequencies, f d T = 0.02 and f d T = 0.05. The number of employed OFDM symbols is set to J = 100, and the number of non-zero virtual angles |I| = 3. It can be seen that, for the higher normalized Doppler frequency, we need to employ more basis functions to obtain the best  MSE performance. In other words, to provide the best MSE performance, the number of basis functions to be employed increases with the normalized Doppler frequency, which means that with higher user mobility, the more basis functions is required to provide accurate channel estimation.
In Fig. 8, we investigate the number of basis functions required to provide the best MSE performance under different SNR scenarios, for the case f d T = 0.02, J = 100, |I| = 3. It can be seen that, as the SNR increases, the number of basis functions required to provide the best MSE performance is increased. This is because when the SNR is low, the main issue for channel estimation is the noise, a small number of basis functions is required to approximate the channel. When the SNR is high, the algorithm can focus more on the time variation of the channel, thus the number of basis functions required will be larger. Hence, we can say that, the number of basis functions required to approximate the time-varying channel should be higher for higher SNR. To provide the best MSE performance, the thresholding factor ξ needs to be properly adjusted. In Fig. 9, we investigate the MSE performance of the DCD-JSR algorithm against the hard thresholding factor ξ , for the case N b = 3, f d T = 0.05, J = 40. It is clear that, for both cases SNR = 10 dB and SNR = 20 dB, as the number of non-zero virtual angles |I| increases, the range for the thresholding factor ξ which can provide the best MSE performance decreases. However, in all the cases, the thresholding factor can be chosen in the interval [0.30, 0.55] to provide the minimum MSE. Fig. 10 shows the MSE performance of the DCD-JSR algorithm in scenarios with different number of non-zero virtual angles against the number of DCD iterations, for the case N b = 2, f d T = 0.02, J = 40. It can be seen that, in all these scenarios, after a few DCD iterations, the algorithm converges to the best MSE. However, the smaller number of non-zero angles, the faster is the convergence. For |I| ≤ 9, a single DCD iteration is enough for the convergence.
In Fig. 11(a), we show the MSE performance for different number of employed antennas, for the case N b = 3, f d T = 0.05, J = 20, |I| = 4. It can be seen that, with a small number of antennas, the MSE performance of the DCD-JSR algorithm is poor. For SNR = 10 dB, it requires M = 56 to approach the oracle performance, for SNR = 20dB, it requires at least M = 32. In Fig. 11 (b), we show the probability of perfect support estimation against the number of employed antennas, where a perfect support estimation means the estimated support is exactly the same as the true support, for the case |I| = 4, f d T = 0.05, J = 20. It can be seen that, at SNR = 10 dB, with a small number of employed antennas, we cannot estimate the support correctly until M = 64; this explains why the MSE performance is poor with small number of antennas. We have run our simulations up to M = 512 and observed that the MSE performance does not change.
In Fig. 12(a) and Fig. 12(b), for the DCD-JSR algorithm with N b = 1 and N b = 2, and the distributed sparsity  adaptive matching pursuit algorithm (DSAMP) [7], we show the MSE performance for different number of non-zero virtual angles, and the probability of perfect support estimation, respectively; here, the DCD-JSR algorithm with N b = 1 corresponds to the version of the DCD-JSR algorithm previously proposed in [15] for time-invariant channels. For simulation scenario, we consider the normalized Doppler frequency f d T = 0.02, J = 40, SNR = 20 dB. It can be seen that, in Fig. 12(a), when f d T = 0.02, since the time variation of the channel is slow, we can estimate the channel quite well using only one basis function in the DCD-JSR algorithm or using the DSAMP algorithm, while both of them shows the MSE performance close to the oracle performance. The DCD-JSR algorithm with N b = 2 also shows close to the oralce MSE performance. In Fig. 12(b), it is seen that for the DCD-JSR algorithm with N b = 2, the support estimation is slightly better than that for the DCD-JSR algorithm with N b = 1 and DSAMP algorithm, which explains why the DCD-JSR algorithm with N b = 2 can provide a better MSE performance in this case.
In Fig. 13, a lower SNR is considered compared to Fig. 12, SNR = 10 dB. It can be seen that, when the noise level is higher, the MSE performance provided by the DCD-JSR algorithm with N b = 2 still shows close to the oracle performance and provides the perfect support estimation, whereas the DCD-JSR algorithm with a single basis function (N b = 1) and the DSAMP algorithm, both developed for time-invariant channels, show inferior MSE performance and support estimation, although the DCD-JSR algorithm with N b = 1 still outperforms the DSAMP algorithm.
In Fig. 14(a) and Fig. 14(b), for the DCD-JSR algorithm, and the DSAMP algorithm, we show the MSE performance for different number of non-zero virtual angles, and the probability of perfect support estimation, respectively. For simulation scenario, we consider the normalized Doppler frequency f d T = 0.05, J = 40, and SNR = 20 dB. It can be seen in Fig. 14(a) that the DCD-JSR algorithm shows close to the oracle MSE performance, while the DSAMP algorithm shows a poor performance. In Fig. 14(b), it is seen that the DCD-JSR algorithm with N b = 3 always provides the perfect support estimation, while the DCD-JSR algorithm with only one basis function N b = 1 shows inferior performance, and the DSAMP algorithm cannot estimate the support accurately. This is because the DSAMP algorithm is developed for static channel, and when f d T = 0.05, i.e. the time variation of the channel is fast, the algorithm is incapable of providing a high estimation performance.
In Fig. 15, results are shown for a higher noise level compared to Fig. 14, we set SNR = 10 dB. It can be seen in Fig. 15(a), that the DSAMP algorithm has again a poor MSE performance, while the DCD-JSR algorithm with N b = 1 and N b = 3 shows close to the oracle performance. In Fig. 15(b), it is clear that the DSAMP algorithm cannot provide an accurate support estimation in this case. The probability of perfect support estimation provided by the DCD-JSR algorithm with N b = 1 decreses as the number of non-zero virtual angels increases, while the DCD-JSR algorithm with N b = 3 can always provides the perfect support estimation. This is because as the time variation of the channel becomes faster, more basis functions is required to accurately approximate the channel.
Hence, from Fig. 12 to Fig. 15, we can conclude that the DCD-JSR algorithm outperforms the DSAMP algorithm. The improvement in the performance provided by the DCD-JSR algorithm against the DSAMP algorithm is more  significant in time-varying channels. For faster time varying channels, by employing more basis functions,we can significantly improve the performance of the DCD-JSR algorithm compared to the case N b = 1 developed in [15] for static channels. Fig. 16 shows the computational complexity of the DCD-JSR algorithm and DSAMP algorithm against the number of non-zero virtual angles, obtained for the case f d T = 0.02, J = 20, SNR = 20 dB. It can be seen that the DCD-JSR algorithm has significantly lower computational complexity than the DSAMP algorithm. The DCD-JSR algorithm, when N b = 2, has slightly higher computational complexity than the DCD-JSR algorithm with N b = 1, while the increase in the number of basis functions provides a significantly better MSE performance.

VI. CONCLUSION
In this paper, by combining the BEM approach and the DCD-JSR algorithm, an efficient algorithm for estimation of the fast time-varying channels in virtual angular domain for massive MIMO systems is proposed. Simulation results have shown that compared to the previously proposed algorithms for channel estimation in virtual angular domain designed for time-invariant channels, the proposed DCD-JSR algorithm could provide significantly better MSE performance in timevarying channels.