# Grouping-Based Channel Estimation and Tracking for Millimeter Wave Massive MIMO Systems

Authors: Rui Yin; Xin Zhou; Wei Qi; Celimuge Wu; Yunlong Cai
Journal: Wireless Communications and Mobile Computing
 (2021)

---

## Abstract

Although the millimeter wave (mmWave) massive multiple-input and multiple-output (MIMO) system can potentially boost the network capacity for future communications, the pilot overhead of the system in practice will greatly increase, which causes a significant decrease in system performance. In this paper, we propose a novel grouping-based channel estimation and tracking approach to reduce the pilot overhead and computational complexity while improving the estimation accuracy. Specifically, we design a low-complexity iterative channel estimation and tracking algorithm by fully exploiting the sparsity of mmWave massive MIMO channels, where the signal eigenvectors are estimated and tracked based on the received signals at the base station (BS). With the recovered signal eigenvectors, the celebrated multiple-signal classification (MUSIC) algorithm can be employed to estimate the direction of arrival (DoA) angles and the path amplitude for the user terminals (UTs). To improve the estimation accuracy and accelerate the tracking speed, we develop a closed-form solution for updating the step-size in the proposed iterative algorithm. Furthermore, a grouping method is proposed to reduce the number of sharing pilots in the scenario of multiple UTs to shorten the pilot overhead. The computational complexity of the proposed algorithm is analyzed. Simulation results are provided to verify the effectiveness of the proposed schemes in terms of the estimation accuracy, tracking speed, and overhead reduction.


---

## Body

## 1. Introduction

The fifth generation (5G) wireless communication system has been proposed to meet the increasing demands on spectrum efficiency (SE), energy efficiency (EE), and network efficiency (NE). Millimeter wave (mmWave) communications operating in the frequency band spanning from 30 to 300 GHz have received great attention for 5G and beyond 5G (B5G) systems, since they are readily combined with massive multiple-input multiple-output (MIMO) [1, 2] techniques for dramatically improving the system SE, EE, and latency performance [3, 4].The channel state information (CSI) accuracy is very important for the design of mmWave massive MIMO systems, and it plays a key role to achieve high SE and EE [5, 6]. Particularly, obtaining accurate CSI is the premise to leverage the full benefits brought by the massive MIMO technique on the mmWave frequency band [7]. Conventionally, CSI is commonly obtained with the assistance of the pilot sequences. Attributed to the orthogonality of pilot sequences, the CSI estimation with high accuracy can be obtained. Nonetheless, there is a positive linear relationship between the pilot overhead and the number of antennas. Therefore, when hundreds of antennas are deployed at the BS, the pilot overhead increases linearly, which may deteriorate the system performance significantly [8] due to the pilot contamination.To reduce the pilot overhead, proper allocation of pilot sequences is one of the main solutions. The authors in [9] have designed a pilot reuse algorithm for massive MIMO with in-phase and quadrature-phase imbalances (IQI). Moreover, a pilot reuse method has been proposed in [8], where the channel power distribution in angular domain is derived via the channel covariance matrix analysis. Accordingly, a pilot reuse pattern has been designed to reduce pilot overhead. In [10, 11], a novel dynamic pilot allocation scheme has been proposed with low complexity. A fractional pilot reuse scheme based on threshold optimization has been developed in [12], where different pilots are allocated to the users classified based on the fluid model.In fact, the pilot reuse scheme unavoidably affects the accuracy of CSI estimation, which motivates the blind CSI estimation methods developed to reduce the estimation error without using pilot sequences, such as the subspace method (SM) and monte carlo methods (MCMs). However, these conventional methods have extremely high computational complexity due to the matrix decomposition and are not suitable to implement in practice. To reduce the computational complexity, a spatial basis expansion model (SBEM) has been proposed to model the mmWave massive MIMO channel with a uniform linear array (ULA) of antennas in [13]. Then, a unified transmission strategy for the mmWave massive MIMO systems in the spatial domain has been proposed, which includes low complex CSI estimation and user scheduling. The same channel model has been extended to the uniform rectangular array- (URA-) based massive MIMO systems in [14]. By combining two-dimensional fast Fourier transformation (2D-FFT) and spatial spectrum analysis, the CSI can be estimated with low-complexity based on 2D-FFT operation in [15].By exploiting the sparse characteristic of mmWave channel, an open-loop channel estimator based on orthogonal matching pursuit has been proposed for mmWave hybrid MIMO system [16]. Based on the channel reciprocity characteristics in TDD systems, the Arnoldi iteration has been used to estimate the left and right singular subspaces of the mmWave MIMO channel in [17]. To reduce the computational complexity, a sequential and adaptive subspace estimation (SASE) method has been proposed to estimate the column and row subspaces of the mmWave MIMO channel in [18].Most of the aforementioned works cannot track the variation on the CSI, which needs to be reimplemented after the channel coherence time to catch up with the varying channel coefficients. Therefore, tracking and simultaneously estimating the CSI can further reduce the computational complexity to achieve accurate CSI. In [19], a modified unscented Kalman filter has been proposed to track the DoA angles with known channel power gain as a priori information. In [20], the predetermined user motion model has been exploited to predict the change in the CSI of massive MIMO systems. The assumptions in both works hinder their practical operation. A sparse complex-valued neural network has been applied to approximate the uplink-to-downlink mapping function in [21]. However, the overhead required for the uplink channel estimation is also intolerable. In [22], a principal subspace analysis- (PSA-) based iterative method has been applied to estimate and track the CSI. However, the computational complexity of this method is still very high due to the eigenvalue decomposition (EVD), and the step size of the iterative algorithm is not optimized.To reduce the computational complexity and the pilot overhead, we propose a novel channel estimation and tracking algorithm for the mmWave massive MIMO systems with URA antennas at the BS. First, the channel sparsity of the mmWave frequency channel is fully exploited. Based on the principal component analysis (PCA) method, a weighted iterative method with the optimized step size is proposed to estimate and track the signal eigenvectors. Second, the obtained eigenvectors are used to estimate DoA angles and path amplitudes via the classic multiple-signal classification (MUSIC) algorithm. Then, a standard inner product is used to match the elevation DoA angles with the right azimuth DoA angles on the same path. Moreover, a user terminal (UT) grouping scheme is proposed to realize the training sequence reuse in the scenario of multiple UTs, which substantially reduces the number of required pilots and the pilot contamination. The main contributions of this work can be summarized as follows.(1)
Based on the PCA, a weighted iterative method is proposed to estimate and track the eigenvectors of the signal subspace, which are used to recover the CSI. Since the complex matrix decomposition operation is avoided, the computational complexity decreases significantly(2)
To accelerate the convergence speed for the proposed iterative scheme, we develop a closed-from solution for updating the step-size in the proposed iterative algorithm(3)
To reduce the pilot overhead, a UT grouping method based on the spatial separability is proposed to divide the UTs into different groups. The spatial separability implies that the UTs in the same group do not have overlapping paths so that they can share the same pilot sequence and transmit simultaneously without causing interferenceThis paper is organized as follows. In Section2, a mmWave massive MIMO system is described. The singular value decomposition (SVD) method and the method to estimate the CSI based on the signal eigenvectors are introduced in Section 3. In Section 4, a novel iterative method is proposed. The simulation results are demonstrated in Section 6. Finally, conclusions are drawn in Section 7.


## 2. System Model

In this section, we describe the wireless channel model of the investigated mmWave massive MIMO system. As shown in Figure1, the BS equipped with a URA of M×N antennas serves UTs, each of which is equipped with a single antenna. Moreover, the wavelength is denoted as λ, and the interantenna spacing of both elevation and horizontal directions is denoted as ds which is set to be half of the λ.Figure 1
Massive MIMO system with URA of antennas at the BS.We assume that the signal transmitted from each UT,i∈I≜1,⋯,U where U is the number of UTs, to the BS arrives along Pi clusters consisting of a number of independent resolvable paths, and each path has a corresponding elevation DoA angle and azimuth DoA angle. Moreover, the mean and angle spread (AS) of the elevation DoA angle of the p-th cluster are denoted by θp and Δθp, respectively. Similarly, the mean and AS of the azimuth DoA angle of the p-th cluster are denoted by ϕp and Δϕp, respectively. When the signal propagates along the pth cluster, the channel can be written by integrating over the angular region as(1)Hp=∫θp−Δθpθp+Δθp∫ϕp−Δϕpϕp+Δϕpβθ,ϕaθbTθ,ϕdθdϕ,where Hp denotes the channel matrix corresponding to the pth cluster, βθ,ϕ=αθ,ϕexp−jΦθ,ϕ, αθ,ϕ is the attenuation and Φθ,ϕ is the initial phase, and aθ and bTθ,ϕ are the array manifold vectors (AMVs), which are written as(2)aθ=1,ejχcosθ,⋯,ejχM−1cosθT,(3)bθ,ϕ=1,ejχsinθcosϕ,⋯,ejχN−1sinθcosϕT,where χ≜2πd/λ.Then, we define a steering matrix denoted asAθ,ϕ=aθbTθ,ϕ∈ℂM×N on the basis of AMVs, which can be written as(4)1⋯ejχN−1sinθcosϕejχcosθ⋯ejχcosθ+N−1sinθcosϕ⋮⋮⋮ejχM−1cosθ⋯ejχM−1cosθ+N−1sinθcosϕ.Based on (4), the channel for the p-th cluster can be rewritten as(5)Hp=∫θp−Δθpθp+Δθp∫ϕp−Δϕpϕp+Δϕpβθ,ϕAθ,ϕdθdϕ.Based on (5), we can observe that the channel matrix, Hp, mainly depends on βθ,ϕ and Aθ,ϕ, where the former is directly related to the channel attenuation and initial phase, and the latter is determined by the DoA angles. Note that if Aθ,ϕ and βθ,ϕ are estimated with high accuracy, the CSI matrix would be reconstructed with high quality. To recover the CSI, there are six parameters to be estimated, i.e., θp, Δθp, ϕp, Δϕp, αθ,ϕ, and Φθ,ϕ.Actually, (5) can be simplified when a mmWave frequency band is used in the massive MIMO system. From [2], Δθp and Δϕp tend to be zero in the scenario where the massive MIMO system operates over the mmWave frequency band. Therefore, each cluster can be seen as an independent path, and each path is formed by the line-of-sight (LOS) and the first-order reflection components. Consequently, (5) can be rewritten as(6)Hp=βθp,ϕpAθp,ϕp.Then, the received signals at the BS is expressed as(7)Yk=∑i=1U∑p=1Piβθi,p,ϕi,pAθi,p,ϕi,pxi,k+W=∑i=1U∑p=1PiHi,pxi,k+W,where k is the index of the received signal snapshot, xi,k is the transmitted signal from UT i, W is the Additive White Gaussian Noise (AWGN) matrix, and Hi,p∈ℂM×N is the channel matrix for UT i on the pth path. From (7), we can observe that to recover the channel matrix, only four parameters need to be estimated for each path p, which are θi,p, ϕi,p, αθi,p,ϕi,p, and Φθi,p,ϕi,p. Therefore, by comparing with estimating entries in the original channel matrix, the computational complexity of recovering the parameters related to the paths is significantly reduced.


## 3. PCA-Based Eigenvectors Estimation and Tracking

In this section, a PCA-based iterative method to estimate and track the signal eigenvectors is proposed.


### 3.1. Channel Decomposition

Conventionally, SVD operation can be applied to estimate the eigenvectors related to the channel, which can be used to recover the channel. To explain the related notations and theories clearly, we assume that there is one UT in the system withP distinguishable paths first.According to the SVD operation, the received signal,Yk, can be written as Yk=SVDH, where S is a unitary matrix S∈ℂM×M, D is a unitary matrix with D∈ℂN×N, the column vectors of S and D are all referred to as singular vectors, and V is a diagonal matrix, V∈ℂM×N, which is named as singular value matrix and sorted in descending order.Besides, when using the EVD to decompose the covariance matrices,CL=EYkYkH and CR=EYkHYk, the column vectors of S and D are the eigenvectors of CL and CR, respectively. Based on the array signal processing knowledge, the first P column vectors of S and D are mainly dependent on elevation and azimuth DoA angles of P resolvable paths, which are denoted as SP and DP, respectively. Accordingly, the DoA angles and amplitudes of resolvable paths can be estimated after recovering the SP and DP.Based on above investigation, by applying the EVD operation, the covariance matrix,CL, can be written as(8)CL=SPVsSPH+SzVzLSzH,where matrixSP∈ℂM×P, and the column vectors of SP are eigenvectors, Vs is a diagonal matrix and Vs∈ℂP×P with entries from λ1 to λP, SZ is a matrix ∈ℂM×M−P, and the column vectors of SZ are eigenvectors, and VzL is a diagonal matrix VzL∈ℂM−P×M−P with entries of eigenvalues from λP+1 to λM.Similarly, applying EVD operation toCR, we can obtain(9)CR=DPVdDPH+DzVzRDzH,where DP∈ℂN×P, and the column vectors in DP are eigenvectors, Vd is a diagonal matrix, Vd∈ℂP×P, with entries from γ1 to γP, DZ∈ℂM×N−P, and the column vectors of SZ are eigenvectors, VzR is a diagonal matrix, and VzR∈ℂN−P×N−P, with entries of eigenvalues from γP+1 to γM.Even thoughSP and DP can be recovered via the SVD or EVD operation as analyzed above, each of them has very high computational complexity, which is OM3 and OM3+N3, respectively. Therefore, they are not suitable for channel estimation in practice. Moreover, it cannot trace the change on the channel matrix caused by the mobility of UTs.


### 3.2. Adaptive PCA-Based Iterative Method

In order to estimate and track the eigenvectors with low computational complexity, an adaptive PCA-based iterative method is proposed in this subsection. First, we take the estimation and tracking ofSP for instance. Define a matrix T1=SPU∈ℂM×P, where U is a unitary matrix with P×P dimension. The estimation on T1 can be formulated as the following optimization problem(10)minT1JT1=EtrYk−T1T1HYk2=trCL−2trT1HCLT1+trT1HCLT1T1HT1,which is to minimize the mean square error (MSE) of the estimation. A simple iterative gradient descent method can be used to solve the optimization problem in (10). The gradient of JT1 with respect to T1 is updated by(11)T1,j=T1,j−1−η∂JT1∂T1T1=T1,j−1,where j denotes the iteration index, η>0 is a step size, and ∂JT1/∂T1 is to represent the gradient of JT1 with respect to T1. To obtain the partial derivative of JT1 with respect to T1, we first calculate the differential of JT1 with respect to T1 via the chain rule, which is given by(12)dJT1=−2trT1HCLH+CLdT1+trdT1HCLT1T1HT1+trT1HCLdT1T1HT1+trT1HCLT1dT1HT1+trT1HCLT1T1HdT1=−4trT1HCLdT1+2trT1HT1T1HCLdT1+2trT1HCLT1T1HdT1.Based on the matrix analysis, the gradient ofJT1 is given by(13)∂JT1∂T1=−22CL−CLT1T1H−T1T1HCLT1.After substituting (13) into (11), the update rule on T1 can be rewritten as(14)T1,j=T1,j−1+2η2CL−CLT1,j−1T1,j−1H−T1,j−1T1,j−1HCLT1,j−1.With Equation (14), T1 will converge to the principal subspace, SPU, which is not equal to SP due to the unitary matrix U.


### 3.3. Optimization Problem

To estimate the corresponding eigenvectors, we reconstruct the optimization problem based on the weighted method proposed in [23]. The new optimization problem is formulated as(15)minT1JT1=trCL−2trΩT1HCLT1+trΩT1HCLT1T1HT1,where Ω is a diagonal matrix whose diagonal item is wi,i. If the diagonal items of Ω are set with different and positive values and w1,1>w2,2>⋯>wP,P, the estimation of T1 will converge to SP [23]. Based on (15), the gradient of JT1 can be rewritten as(16)∂JT1∂T1=−4CLT1Ω−CLT1ΩT1HT1−T1T1HCLT1Ω−CLT1T1HT1Ω−T1ΩT1HCLT1.With (16), matrix T1 in (14) is updated by(17)T1,j=T1,j−1+η4CLT1,k−1Ω−CLT1,j−1ΩT1HT1,j−1−T1,j−1T1,j−1HCLT1,j−1Ω−CLT1,j−1T1,j−1HT1,j−1Ω−T1,j−1ΩT1,j−1HCLT1,j−1.It has been proven in [23] that T1 converges to SP when j tends to infinity. From (17), we can observe that the covariance matrix CL is essential to update T1, but the accurate CL is difficult to be obtained due to CL is the expectation of YkYkH. To obtain CL, an approximated method has been proposed in [24], which uses a sampling covariance matrix C^L,j. The C^L,j is given by(18)C^L,j=1j∑k=0jτj−iYkYkH=τj−1jC^L,j−1+1jYjYjH,where τ∈0,1 is the forgetting factor, the variable k is used to indicate the kth time slot, and j is to indicate the current time slot. The initial matrix C^L,1 is set to be equal to Y1Y1H. Moreover, it is noteworthy that the parameter τ can be adjusted to trace channel variance. If the channel is modeled as a stationary process, (18) can be regarded as the calculation on the average value of YkYkH, and τ is equal to 1. However, the τ is taken in the range of (0,1) due to the massive MIMO channel is not stationary in practice.With (18), the updating rule for T1 can be rewritten as(19)T1,k=T1,k−1+η4C^1,k−1T1,k−1Ω−C^L,k−1T1,k−1ΩT1HT1,k−1−T1,k−1T1,k−1HC^1,k−1T1,k−1Ω−C^1,k−1T1,k−1T1,k−1HT1,k−1Ω−T1,k−1ΩT1,k−1HC^1,kT1,k−1.It is noteworthy that the step size,η, is essential to obtain T1 with high accuracy and convergence speed. However, deriving an appropriate step size η is nontrivial. On the one hand, if η is small, the convergence speed will be slow. On the other hand, the accuracy will be impaired if η is high. Therefore, in the next subsection, the optimal η is derived to accelerate the convergence speed while keeping the accuracy.


### 3.4. Optimal Step Size

To find the optimal step sizeηopt, we consider the objection function JT1 at the k+1th iteration step which is given by(20)JT1,k+1=trR1,k+1−2trΩT1,k+1HR1,k+1T1,k+1+trΩT1,k+1HR1,k+1T1,k+1T1,k+1HT1,k+1.To simplify the equation, we replace∂JT1/∂T1 with ∇. Then, combining (11) with (20), (20) can be rewritten as(21)JT1,k+1,η=trC^L,k+1−2trΩT1,k−η∇HC^L,k+1T1,k−η∇+trΩT1,k−η∇HC^L,k+1T1,k−η∇T1,k−η∇HT1,k−η∇.SinceT1,k+1HT1,k+1=UHSPHSPU=I, (21) can be written as(22)JT1,k+1,η=trC^L,k+1−trΩT1,k−η∇HC^L,k+1T1,k−η∇.Therefore,JT1,k+1 has a global maximum with the optimal step size, ηopt. Based on (22), we can differentiate JT1,k+1,η with respect to η, which is given by(23)∂JT1,k+1,η∂η=∂∂η−2trΩT1,k−η∇HC^L,k+1T1,k−η∇+trΩT1,k−η∇HC^L,k+1T1,k−η∇T1,k−η∇HT1,k−η∇=−2trΩ−∇HC^1,k+1T1,k−η∇−2trΩT1,k−∇HC^L,k+1−η∇=2trΩ∇HC^L,k+1T1,k−2ηtrΩ∇HC^L,k+1∇+2trΩT1,kHC^L,k+1∇−trΩ∇HC^L,k+1∇.When∂JT1,k+1/∂η=0, JT1,k+1 has a global maximum. As a result, the optimal step size, ηopt, is the solution of ∂JT1,k+1/∂η=0, which is given by(24)ηopt=trΩ∇HC^L,k+1T1,k+trΩT1,kHC^L,k+1∇2trΩ∇HC^L,k+1∇.Based on (24), we can achieve the optimal updating rule to estimate and track the signal eigenvector matrix, SP. In the same way, we can estimate and track the DP.


## 3.1. Channel Decomposition

Conventionally, SVD operation can be applied to estimate the eigenvectors related to the channel, which can be used to recover the channel. To explain the related notations and theories clearly, we assume that there is one UT in the system withP distinguishable paths first.According to the SVD operation, the received signal,Yk, can be written as Yk=SVDH, where S is a unitary matrix S∈ℂM×M, D is a unitary matrix with D∈ℂN×N, the column vectors of S and D are all referred to as singular vectors, and V is a diagonal matrix, V∈ℂM×N, which is named as singular value matrix and sorted in descending order.Besides, when using the EVD to decompose the covariance matrices,CL=EYkYkH and CR=EYkHYk, the column vectors of S and D are the eigenvectors of CL and CR, respectively. Based on the array signal processing knowledge, the first P column vectors of S and D are mainly dependent on elevation and azimuth DoA angles of P resolvable paths, which are denoted as SP and DP, respectively. Accordingly, the DoA angles and amplitudes of resolvable paths can be estimated after recovering the SP and DP.Based on above investigation, by applying the EVD operation, the covariance matrix,CL, can be written as(8)CL=SPVsSPH+SzVzLSzH,where matrixSP∈ℂM×P, and the column vectors of SP are eigenvectors, Vs is a diagonal matrix and Vs∈ℂP×P with entries from λ1 to λP, SZ is a matrix ∈ℂM×M−P, and the column vectors of SZ are eigenvectors, and VzL is a diagonal matrix VzL∈ℂM−P×M−P with entries of eigenvalues from λP+1 to λM.Similarly, applying EVD operation toCR, we can obtain(9)CR=DPVdDPH+DzVzRDzH,where DP∈ℂN×P, and the column vectors in DP are eigenvectors, Vd is a diagonal matrix, Vd∈ℂP×P, with entries from γ1 to γP, DZ∈ℂM×N−P, and the column vectors of SZ are eigenvectors, VzR is a diagonal matrix, and VzR∈ℂN−P×N−P, with entries of eigenvalues from γP+1 to γM.Even thoughSP and DP can be recovered via the SVD or EVD operation as analyzed above, each of them has very high computational complexity, which is OM3 and OM3+N3, respectively. Therefore, they are not suitable for channel estimation in practice. Moreover, it cannot trace the change on the channel matrix caused by the mobility of UTs.


## 3.2. Adaptive PCA-Based Iterative Method

In order to estimate and track the eigenvectors with low computational complexity, an adaptive PCA-based iterative method is proposed in this subsection. First, we take the estimation and tracking ofSP for instance. Define a matrix T1=SPU∈ℂM×P, where U is a unitary matrix with P×P dimension. The estimation on T1 can be formulated as the following optimization problem(10)minT1JT1=EtrYk−T1T1HYk2=trCL−2trT1HCLT1+trT1HCLT1T1HT1,which is to minimize the mean square error (MSE) of the estimation. A simple iterative gradient descent method can be used to solve the optimization problem in (10). The gradient of JT1 with respect to T1 is updated by(11)T1,j=T1,j−1−η∂JT1∂T1T1=T1,j−1,where j denotes the iteration index, η>0 is a step size, and ∂JT1/∂T1 is to represent the gradient of JT1 with respect to T1. To obtain the partial derivative of JT1 with respect to T1, we first calculate the differential of JT1 with respect to T1 via the chain rule, which is given by(12)dJT1=−2trT1HCLH+CLdT1+trdT1HCLT1T1HT1+trT1HCLdT1T1HT1+trT1HCLT1dT1HT1+trT1HCLT1T1HdT1=−4trT1HCLdT1+2trT1HT1T1HCLdT1+2trT1HCLT1T1HdT1.Based on the matrix analysis, the gradient ofJT1 is given by(13)∂JT1∂T1=−22CL−CLT1T1H−T1T1HCLT1.After substituting (13) into (11), the update rule on T1 can be rewritten as(14)T1,j=T1,j−1+2η2CL−CLT1,j−1T1,j−1H−T1,j−1T1,j−1HCLT1,j−1.With Equation (14), T1 will converge to the principal subspace, SPU, which is not equal to SP due to the unitary matrix U.


## 3.3. Optimization Problem

To estimate the corresponding eigenvectors, we reconstruct the optimization problem based on the weighted method proposed in [23]. The new optimization problem is formulated as(15)minT1JT1=trCL−2trΩT1HCLT1+trΩT1HCLT1T1HT1,where Ω is a diagonal matrix whose diagonal item is wi,i. If the diagonal items of Ω are set with different and positive values and w1,1>w2,2>⋯>wP,P, the estimation of T1 will converge to SP [23]. Based on (15), the gradient of JT1 can be rewritten as(16)∂JT1∂T1=−4CLT1Ω−CLT1ΩT1HT1−T1T1HCLT1Ω−CLT1T1HT1Ω−T1ΩT1HCLT1.With (16), matrix T1 in (14) is updated by(17)T1,j=T1,j−1+η4CLT1,k−1Ω−CLT1,j−1ΩT1HT1,j−1−T1,j−1T1,j−1HCLT1,j−1Ω−CLT1,j−1T1,j−1HT1,j−1Ω−T1,j−1ΩT1,j−1HCLT1,j−1.It has been proven in [23] that T1 converges to SP when j tends to infinity. From (17), we can observe that the covariance matrix CL is essential to update T1, but the accurate CL is difficult to be obtained due to CL is the expectation of YkYkH. To obtain CL, an approximated method has been proposed in [24], which uses a sampling covariance matrix C^L,j. The C^L,j is given by(18)C^L,j=1j∑k=0jτj−iYkYkH=τj−1jC^L,j−1+1jYjYjH,where τ∈0,1 is the forgetting factor, the variable k is used to indicate the kth time slot, and j is to indicate the current time slot. The initial matrix C^L,1 is set to be equal to Y1Y1H. Moreover, it is noteworthy that the parameter τ can be adjusted to trace channel variance. If the channel is modeled as a stationary process, (18) can be regarded as the calculation on the average value of YkYkH, and τ is equal to 1. However, the τ is taken in the range of (0,1) due to the massive MIMO channel is not stationary in practice.With (18), the updating rule for T1 can be rewritten as(19)T1,k=T1,k−1+η4C^1,k−1T1,k−1Ω−C^L,k−1T1,k−1ΩT1HT1,k−1−T1,k−1T1,k−1HC^1,k−1T1,k−1Ω−C^1,k−1T1,k−1T1,k−1HT1,k−1Ω−T1,k−1ΩT1,k−1HC^1,kT1,k−1.It is noteworthy that the step size,η, is essential to obtain T1 with high accuracy and convergence speed. However, deriving an appropriate step size η is nontrivial. On the one hand, if η is small, the convergence speed will be slow. On the other hand, the accuracy will be impaired if η is high. Therefore, in the next subsection, the optimal η is derived to accelerate the convergence speed while keeping the accuracy.


## 3.4. Optimal Step Size

To find the optimal step sizeηopt, we consider the objection function JT1 at the k+1th iteration step which is given by(20)JT1,k+1=trR1,k+1−2trΩT1,k+1HR1,k+1T1,k+1+trΩT1,k+1HR1,k+1T1,k+1T1,k+1HT1,k+1.To simplify the equation, we replace∂JT1/∂T1 with ∇. Then, combining (11) with (20), (20) can be rewritten as(21)JT1,k+1,η=trC^L,k+1−2trΩT1,k−η∇HC^L,k+1T1,k−η∇+trΩT1,k−η∇HC^L,k+1T1,k−η∇T1,k−η∇HT1,k−η∇.SinceT1,k+1HT1,k+1=UHSPHSPU=I, (21) can be written as(22)JT1,k+1,η=trC^L,k+1−trΩT1,k−η∇HC^L,k+1T1,k−η∇.Therefore,JT1,k+1 has a global maximum with the optimal step size, ηopt. Based on (22), we can differentiate JT1,k+1,η with respect to η, which is given by(23)∂JT1,k+1,η∂η=∂∂η−2trΩT1,k−η∇HC^L,k+1T1,k−η∇+trΩT1,k−η∇HC^L,k+1T1,k−η∇T1,k−η∇HT1,k−η∇=−2trΩ−∇HC^1,k+1T1,k−η∇−2trΩT1,k−∇HC^L,k+1−η∇=2trΩ∇HC^L,k+1T1,k−2ηtrΩ∇HC^L,k+1∇+2trΩT1,kHC^L,k+1∇−trΩ∇HC^L,k+1∇.When∂JT1,k+1/∂η=0, JT1,k+1 has a global maximum. As a result, the optimal step size, ηopt, is the solution of ∂JT1,k+1/∂η=0, which is given by(24)ηopt=trΩ∇HC^L,k+1T1,k+trΩT1,kHC^L,k+1∇2trΩ∇HC^L,k+1∇.Based on (24), we can achieve the optimal updating rule to estimate and track the signal eigenvector matrix, SP. In the same way, we can estimate and track the DP.


## 4. Estimation on the Channel Parameters

After estimatingSP and DP in Section 3, the DoA angles of the resolvable paths on both elevation and azimuth directions can be recovered based on the classic MUSIC algorithm. Define up and vp as(25)up=χcosθp,(26)vp=χsinθpcosϕp.With these two variables, the steering matrix,Aθ,ϕ=aθbTθ,ϕ, can be rewritten as Aup,vp=aupbTvp. Then, the MUSIC method is used to estimate up and vp. Then, the DoA angles can be recovered according to (25) and (26).


### 4.1. Estimation onup and vp

MultiplyingCL=EYkYkH and (8) with Sz, we obtain(27)CLSz=∑p=1PHpψiHpHSz+σ2Sz,(28)CLSz=SPΛsSPHSz+SzΛz1SzHSz=σ2Sz.where ψi=ExixiH.Then, subtracting (27) to (28), there is(29)HpHSz=0.Sincehp=βθp,ϕpAup,vp, we can derive the following equation from (29)(30)Aup,vpHSz=0.Based on matrix manipulation, from (30), we can obtain(31)aHupsz,k=∑m=0M−1ejmupsz,k,m=0,where sz,k is the kth column of matrix Sz and sz,k,m is the mth entry of sz,k.Based on (31), the pseudospectrum with respect to u can be derived as(32)fu=1∑k=P+1MaupHsz,k2=1pHζSzSzHpζ,wherepζ=1,ζ,⋯,ζM−1T. According to (32), if ζ=eju, the roots of the denominator would locate on the unit circle and be the estimation on u since the steering vector, aup,1≤p≤P, is orthogonal to the noise subspace, Sz. Moreover, SzSzH can be derived based on the estimation on SP in Section 2, which is given by(33)SzSzH=I−SPSPH.Then, inserting (33) into (32), up, ∀1≤i≤P, can be estimated by finding the peaks of the pseudospectrum in (32). By using the same method, we can estimate vp, ∀1≤p≤P, with DP.


### 4.2. Matching onup and vp

After estimatingup and vp in the above subsection, it needs to match up with the right vp on the same path since each path has a corresponding pair of elevation and azimuth DoA angles. To do so, a standard inner product is defined to match up with vp, which is given by(34)Pr=<J,K>=trJHK,where J and K are two matrices. The inner product is to calculate the energy of the received signals on certain direction. Accordingly, with the estimation of up and vp on the above subsection, the inner products of received signals matrix Y with the steering matrix Aui,vj is given by(35)Pr=<Aui,vj,Y>=trAHui,vjY=trbvj∗aHuiY=aHuiYbvj∗,where i,j=1,2,3⋯P. Pr will be the maximum value when i=j. The received signals from P paths make the most contribution to the received energy at the BS. Then, we can match up with vp by searching the maximum value of Pr. After that, based on (25) and (26), the DoA angles of each path can be estimated by the following formulas(36)θp=arccosupχ,∀1≤p≤P,(37)ϕp=arccosvpχ1−uP/χ2,∀1≤p≤P.


### 4.3. Estimation onβθ,ϕ

Based on (7) and SVD operation, the received signals Yk is rewritten as(38)Yk=SΣDH=∑p=1Pβθp,ϕpAup,vpxk+W=∑p=1PρpspdpH.where ρp is the pth eigenvalue of Yk; sp and dp are the pth eigenvectors of SP and DP, respectively. Then, we set Gp=spdpH and construct an orthonormal space as follows(39)<Gi,Gj>=trGiHGj=0,ifi≠j,1,ifi=j.According to (39), the pth eigenvalue, ρp, can be written as(40)ρp=<Gp,Yk>=<Gp,∑p=1Pβθp,ϕpAup,vpx+W>=∑p=1Pβθp,ϕpspHAup,vpdpx+spHWdp=∑p=1Pβθp,ϕpspHAup,vpdpx+σ2.Based on (40), define(41)O=dets1HAu1,v1d1s1HAu2,v2d1⋯s1HAuP,vPd1s2HAu1,v1d2s2HAu2,v2d2⋯s2HAuP,vPd2⋮⋮⋱⋮sPHAu1,v1dPsPHAu2,v2dP⋯sPHAuP,vPdP,(42)Op=dets1HAu1,v1d1⋯ρ1−σ⋯s1HAuP,vPd1s2HAu1,v1d2⋯ρ2−σ⋯s2HAPuP,vPd2⋮⋮⋮⋱⋮sPHAu1,v1dP⋯ρP−σ⋯sPHAuP,vPdP.Based on the Cramer’s Rule, the amplitude of pathi can be calculated as(43)βθp,ϕp=OpO,1≤p≤P.Then, the complete algorithm flow of the proposed scheme is given in Algorithm1.Algorithm 1: Channel estimation and tracking scheme
1: Receive thekth signal snapshot in the kth time slot, Yk. Then, the BS updates the covariance matrix, C^L,k, based on (18);2: Estimate and track the eigenvector matrix of the received signal,SP, based on the PCA algorithm in (19). Recover DP in the same way;3: Estimateup and vp,1≤i≤P, based on (32);4: Pairup,1≤p≤P, with vp,1≤p≤P, based on (35);5: Estimateθp,1≤p≤P and ϕp,1≤p≤P based on (36) and (37);6: Estimate the amplitude corresponding to each path based on (43);7: Recover the channel based on (6);


## 4.1. Estimation onup and vp

MultiplyingCL=EYkYkH and (8) with Sz, we obtain(27)CLSz=∑p=1PHpψiHpHSz+σ2Sz,(28)CLSz=SPΛsSPHSz+SzΛz1SzHSz=σ2Sz.where ψi=ExixiH.Then, subtracting (27) to (28), there is(29)HpHSz=0.Sincehp=βθp,ϕpAup,vp, we can derive the following equation from (29)(30)Aup,vpHSz=0.Based on matrix manipulation, from (30), we can obtain(31)aHupsz,k=∑m=0M−1ejmupsz,k,m=0,where sz,k is the kth column of matrix Sz and sz,k,m is the mth entry of sz,k.Based on (31), the pseudospectrum with respect to u can be derived as(32)fu=1∑k=P+1MaupHsz,k2=1pHζSzSzHpζ,wherepζ=1,ζ,⋯,ζM−1T. According to (32), if ζ=eju, the roots of the denominator would locate on the unit circle and be the estimation on u since the steering vector, aup,1≤p≤P, is orthogonal to the noise subspace, Sz. Moreover, SzSzH can be derived based on the estimation on SP in Section 2, which is given by(33)SzSzH=I−SPSPH.Then, inserting (33) into (32), up, ∀1≤i≤P, can be estimated by finding the peaks of the pseudospectrum in (32). By using the same method, we can estimate vp, ∀1≤p≤P, with DP.


## 4.2. Matching onup and vp

After estimatingup and vp in the above subsection, it needs to match up with the right vp on the same path since each path has a corresponding pair of elevation and azimuth DoA angles. To do so, a standard inner product is defined to match up with vp, which is given by(34)Pr=<J,K>=trJHK,where J and K are two matrices. The inner product is to calculate the energy of the received signals on certain direction. Accordingly, with the estimation of up and vp on the above subsection, the inner products of received signals matrix Y with the steering matrix Aui,vj is given by(35)Pr=<Aui,vj,Y>=trAHui,vjY=trbvj∗aHuiY=aHuiYbvj∗,where i,j=1,2,3⋯P. Pr will be the maximum value when i=j. The received signals from P paths make the most contribution to the received energy at the BS. Then, we can match up with vp by searching the maximum value of Pr. After that, based on (25) and (26), the DoA angles of each path can be estimated by the following formulas(36)θp=arccosupχ,∀1≤p≤P,(37)ϕp=arccosvpχ1−uP/χ2,∀1≤p≤P.


## 4.3. Estimation onβθ,ϕ

Based on (7) and SVD operation, the received signals Yk is rewritten as(38)Yk=SΣDH=∑p=1Pβθp,ϕpAup,vpxk+W=∑p=1PρpspdpH.where ρp is the pth eigenvalue of Yk; sp and dp are the pth eigenvectors of SP and DP, respectively. Then, we set Gp=spdpH and construct an orthonormal space as follows(39)<Gi,Gj>=trGiHGj=0,ifi≠j,1,ifi=j.According to (39), the pth eigenvalue, ρp, can be written as(40)ρp=<Gp,Yk>=<Gp,∑p=1Pβθp,ϕpAup,vpx+W>=∑p=1Pβθp,ϕpspHAup,vpdpx+spHWdp=∑p=1Pβθp,ϕpspHAup,vpdpx+σ2.Based on (40), define(41)O=dets1HAu1,v1d1s1HAu2,v2d1⋯s1HAuP,vPd1s2HAu1,v1d2s2HAu2,v2d2⋯s2HAuP,vPd2⋮⋮⋱⋮sPHAu1,v1dPsPHAu2,v2dP⋯sPHAuP,vPdP,(42)Op=dets1HAu1,v1d1⋯ρ1−σ⋯s1HAuP,vPd1s2HAu1,v1d2⋯ρ2−σ⋯s2HAPuP,vPd2⋮⋮⋮⋱⋮sPHAu1,v1dP⋯ρP−σ⋯sPHAuP,vPdP.Based on the Cramer’s Rule, the amplitude of pathi can be calculated as(43)βθp,ϕp=OpO,1≤p≤P.Then, the complete algorithm flow of the proposed scheme is given in Algorithm1.Algorithm 1: Channel estimation and tracking scheme
1: Receive thekth signal snapshot in the kth time slot, Yk. Then, the BS updates the covariance matrix, C^L,k, based on (18);2: Estimate and track the eigenvector matrix of the received signal,SP, based on the PCA algorithm in (19). Recover DP in the same way;3: Estimateup and vp,1≤i≤P, based on (32);4: Pairup,1≤p≤P, with vp,1≤p≤P, based on (35);5: Estimateθp,1≤p≤P and ϕp,1≤p≤P based on (36) and (37);6: Estimate the amplitude corresponding to each path based on (43);7: Recover the channel based on (6);


## 5. Extension for Multiple UTs

In this section, the proposed method is extended to the multiple UT scenario. When the UTs have the same paths, their signals cannot be distinguished at the BS. To solve this problem and reduce the pilot overhead, a user grouping scheme is proposed in this section.


### 5.1. Initial DoA Angle Estimation

The initial DoA angles of different UTs need to be known at the BS before dividing them into different groups. To reach this goal, the training process is necessary. Assume that there areI orthogonal training sequences, denoted as r1,⋯,ri,⋯,rI, where ri is an L-dimensional vector with ri2=L. During the training stage, each UT is assigned to a training sequence, and the received signals at BS can be expressed as(44)Y¯=∑i∈UH¯iri+W¯,whereY¯∈ℂM×N×L, and H¯i∈ℂM×N×1 is the uplink channel between the ith UT and the BS; U is the set of UTs. W¯∈ℂM×N×L is the additive white Gaussian noise. The orthogonality of the training sequences can be used to distinguish the received signals from different UTs. First, multiply (44) with riH to get(45)Y¯i=1LY¯riH=1L∑j∈UH¯jrj+W¯riH=H¯i+1LW¯riH.Based on (45), the transmitted signal from different UTs can be separated. With the detached signal Y¯i, the proposed channel estimation method can be used to estimate the initial DoA angles associated with the ith UT. However, the number of UTs may be more than the number of orthogonal training sequences, i.e., U>I. To solve this problem, the pilot sequences can be reused in a time division multiple access (TDMA) manner. UTs are divided into U/I groups, and the UTs in each group can reuse the same I orthogonal pilot sequences in different time slots. Accordingly, the BS can estimate the initial DoA angles of all UTs.


### 5.2. Grouping Scheme

After the initial DoA angle estimation, the UTs can be partitioned into different groups, and the UTs in the same group can use the same training sequence and transmit simultaneously.The principle of the scheme is to guarantee that the UTs in the same group do not have overlapping paths. For example, there are two UTs,u1 and u2 with the DoA angles of θu1, θu2 and ϕu1, ϕu2 estimated in the initial stage. If the minimum difference on the DoA angles between u1 and u2 is bigger than a threshold denoted as λ, minθu1−θu2,ϕu1−ϕu2≥λ, they would be put into the same group. Based on this rule, the UTs can be designated into different groups, and the UTs in the same group have different distinguishable paths. As a result, the UTs in the same group can use the same training sequence, and their CSI can be estimated and tracked by Algorithm 1.Algorithm 2: Grouping Scheme
1: Initialization: divide UTs intoU/I groups. Different UTs in the same group are allocated with different pilot sequences;2: UTs in the same group transmit the pilots in the same time slot. The BS distinguishes the received signals from different UTs;3: The BS estimates the initial DoA angles of different UTs;4: If the minimum difference on the DoA angles between any two UTs is bigger than a threshold denoted asλ, they are put into the same group;5: The UTs in the same group will use the same pilot sequence.With the grouping scheme proposed in Algorithm2, the channel estimation scheme proposed in Algorithm 1 can be extended to the multi-UT scenario. First, the UT grouping scheme is executed at the BS. Then, the UTs in the same group will use the same pilot sequence and transmit simultaneously. After the BS receives signals from the UTs in the same group, it updates the covariance matrix based on (18). Then, it estimates the channel state information based on Algorithm 1 proposed in Section 4. It is noteworthy that the BS treats the UTs in the same group as one “virtual UT.” Since there exists a certain gap in DoA angles among UTs in the same group, the paths associated with their channels can be distinguished at the BS. Based on the initial DoA estimation, the BS can identify to which the UT paths belong. When the BS cannot track the CSI accurately, it needs to rerun the grouping scheme in the algorithm. To guarantee the channel estimation and tracking accuracy, the BS may run Algorithm 2 after certain time slots in practice.


## 5.1. Initial DoA Angle Estimation

The initial DoA angles of different UTs need to be known at the BS before dividing them into different groups. To reach this goal, the training process is necessary. Assume that there areI orthogonal training sequences, denoted as r1,⋯,ri,⋯,rI, where ri is an L-dimensional vector with ri2=L. During the training stage, each UT is assigned to a training sequence, and the received signals at BS can be expressed as(44)Y¯=∑i∈UH¯iri+W¯,whereY¯∈ℂM×N×L, and H¯i∈ℂM×N×1 is the uplink channel between the ith UT and the BS; U is the set of UTs. W¯∈ℂM×N×L is the additive white Gaussian noise. The orthogonality of the training sequences can be used to distinguish the received signals from different UTs. First, multiply (44) with riH to get(45)Y¯i=1LY¯riH=1L∑j∈UH¯jrj+W¯riH=H¯i+1LW¯riH.Based on (45), the transmitted signal from different UTs can be separated. With the detached signal Y¯i, the proposed channel estimation method can be used to estimate the initial DoA angles associated with the ith UT. However, the number of UTs may be more than the number of orthogonal training sequences, i.e., U>I. To solve this problem, the pilot sequences can be reused in a time division multiple access (TDMA) manner. UTs are divided into U/I groups, and the UTs in each group can reuse the same I orthogonal pilot sequences in different time slots. Accordingly, the BS can estimate the initial DoA angles of all UTs.


## 5.2. Grouping Scheme

After the initial DoA angle estimation, the UTs can be partitioned into different groups, and the UTs in the same group can use the same training sequence and transmit simultaneously.The principle of the scheme is to guarantee that the UTs in the same group do not have overlapping paths. For example, there are two UTs,u1 and u2 with the DoA angles of θu1, θu2 and ϕu1, ϕu2 estimated in the initial stage. If the minimum difference on the DoA angles between u1 and u2 is bigger than a threshold denoted as λ, minθu1−θu2,ϕu1−ϕu2≥λ, they would be put into the same group. Based on this rule, the UTs can be designated into different groups, and the UTs in the same group have different distinguishable paths. As a result, the UTs in the same group can use the same training sequence, and their CSI can be estimated and tracked by Algorithm 1.Algorithm 2: Grouping Scheme
1: Initialization: divide UTs intoU/I groups. Different UTs in the same group are allocated with different pilot sequences;2: UTs in the same group transmit the pilots in the same time slot. The BS distinguishes the received signals from different UTs;3: The BS estimates the initial DoA angles of different UTs;4: If the minimum difference on the DoA angles between any two UTs is bigger than a threshold denoted asλ, they are put into the same group;5: The UTs in the same group will use the same pilot sequence.With the grouping scheme proposed in Algorithm2, the channel estimation scheme proposed in Algorithm 1 can be extended to the multi-UT scenario. First, the UT grouping scheme is executed at the BS. Then, the UTs in the same group will use the same pilot sequence and transmit simultaneously. After the BS receives signals from the UTs in the same group, it updates the covariance matrix based on (18). Then, it estimates the channel state information based on Algorithm 1 proposed in Section 4. It is noteworthy that the BS treats the UTs in the same group as one “virtual UT.” Since there exists a certain gap in DoA angles among UTs in the same group, the paths associated with their channels can be distinguished at the BS. Based on the initial DoA estimation, the BS can identify to which the UT paths belong. When the BS cannot track the CSI accurately, it needs to rerun the grouping scheme in the algorithm. To guarantee the channel estimation and tracking accuracy, the BS may run Algorithm 2 after certain time slots in practice.


## 6. Numerical Results

In this section, we provide the numerical results to demonstrate the superiority of the proposed scheme. In the simulation setup, the number of antennas deployed at the BS on both elevation and azimuth directions is 64 (M=N=64), and the interantenna spacing, d, is set to be λ/2. In the first scenario, we assume that there is a single UT under coverage of the BS and the received SNR is 10 dB. Moreover, we assume that the propagation from the UT to BS goes three resolvable paths. The elevation and azimuth angles corresponding to three paths are (25∘, 45∘, 80∘) and (20∘, 30∘, 70∘), respectively.


### 6.1. The Estimation Accuracy and Convergence Speed

First, to evaluate the estimation accuracy, a Tanimoto coefficient is defined which is given by(46)εk=1−1P∑i=1Pt1,i,k·sit1,i,k2+si2−t1,i,k·si,where si is the ith eigenvector of matrix SP, and t1,i,k is the ith eigenvector of the estimation matrix SP at the kth iteration step. Tanimoto coefficient defined in (46) implies that if ε tends to be 0, the estimation approaches to the real SP. To demonstrate the superiority of the optimal step size derived in Section 3.4, the Backtracking Line Search (BLS) method proposed in [25] is used for comparison. From Figure 2, we can observe that ε achieved by the PCA-based algorithm and BLS algorithm converges to 0, which indicates the availability of these algorithm on estimation of eigenvector matrix. Meanwhile, during the first 400 iterations, the ε achieved by the PCA-based algorithm is less than the BLS algorithm. The result indicates that the proposed scheme has the high accuracy under the same iterations. From another point of view, the proposed scheme has the fast convergence speed on estimation of SP.Figure 2
The Tanimoto coefficient comparison.After recovering the eigenvectors, the intermediate variables,up and vp, are estimated and paired according to the proposed scheme. To demonstrate the match process clearly, the inner products of the received signal matrix on the steering matrix are listed in Table 1, where the steering vectors are constructed by the estimated up and vp.Table 1
Matching onup and vp.
<Aup,vp∣Y>v10.6121v20.3991v30.3363u10.906118.903624.552518.9036u20.707163.173815.03032.6220u30.17416.313419.411422.7863From Table1, we can match up with the right vp by finding the maximum value. Consequently, u1, u2, and u3 should be paired with v2, v1, and v3, respectively. After pairing, the estimation of DoA angles and path amplitudes can be obtained through (36), (37), and (43), respectively, which is listed in Table 2, where the index p is to indicate the pth path.Table 2
Estimation onup, vp, and βθp,ϕp.
p123Realθp25.0000°45.0000°80.0000°Estimation ofθp25.0282°45.0005°79.9737°Realϕp20.0000°30.0000°70.0000°Estimation ofϕp19.3768°30.0451°70.0308°Realβθp,ϕp0.60000.80001.0000Estimation ofβθp,ϕp0.59930.80451.0011From Table2, we can observe that the estimations on both θp, ϕp, and βθp,ϕp are close to the real values, which verifies the effectiveness of the proposed scheme. Moreover, since the ϕp is determined by both up and vp, the estimation error on up and vp would affect the estimation accuracy on ϕp. Therefore, the estimation accuracy of the azimuth angles, ϕp, is less than that on the elevation angles, θp.


### 6.2. The Impact of Received SNR on the DoA Estimation

To reveal estimation performance with different SNRs, the normalized mean square error (MSE) on DoA angle estimation is demonstrated in this section. The MSE is defined as(47)MSEθ=1P∑p=1P∣θ¯p−θp∣θp,(48)MSEϕ=1P∑p=1Pϕ¯p−ϕpϕp.Figure3 illustrates the MSE performance on DoA angle estimation. It can be observed that the estimation MSE on DoA angles reduces as the received SNR at the BS increases. Figure 3 also compares the proposed scheme with the BLS algorithm with respect to DoA angle estimation MSE. We can observe that the proposed scheme outperforms the BLS algorithm on the estimation accuracy. Since the estimation on azimuth DoA angle depends on both the estimated up and vp, the estimation accuracy on the elevation DoA angle is higher than that on the azimuth DoA angle.Figure 3
The MSE on DoA estimation.


### 6.3. Tracking Performance

To verify the tracking ability of the proposed scheme, the SVD-based channel estimation is simulated for comparison, where the SVD operation is executed periodically to update the channel estimation. Assuming that the elevation angle changes from 45° to 53°, and azimuth angles change from 30° to 22°, with the changing rate of 0.02°/stepsize.As shown in Figure4, the SVD-based channel estimation method cannot track the change on the channel even though it may estimate channel with high accuracy in some times. On the other hand, the channel estimation based on Algorithm 1 is very close to the perfect one all the time even when the channel changes. Therefore, this simulation result confirms the tracking ability on the channel estimation of Algorithm 1.Figure 4
Tracking on the DoA angles.


### 6.4. The Effectiveness of Grouping Scheme

For the sake of illustration on the effect of the UT grouping scheme, a multiple UT scenario is considered in this section. The elevation and azimuth angles corresponding to the distinguish paths for 10 different UTs are listed in Table3. Assuming that there are only 2 orthogonal pilot sequences and set λ=5°. First, the UTs are assigned to 5 groups. UTs in 5 groups will reuse the 2 orthogonal pilot sequences in a TDMA manner, and the BS can obtain the initial DoA angle estimations of UTs. Afterwards, based on the grouping scheme (GS) in Algorithm 2, UT1, UT2, UT3, UT4, and UT7 should be put into one group, and the remaining UTs should be put into another group. To better demonstrate the effectiveness of GS, another random grouping scheme is adopted for comparison, where the UT1 to UT5 are put into the same group, and the UT6 to UT10 are put into another group. Accordingly, the estimation on DoA angles is shown in Figure 5.Table 3
Add caption.
The number of UT12345678910The corresponding elevation and azimuth angles (θp,ϕp)(25°, 20°)(45°, 30°)(80°, 70°)(37°, 82°)(60°, 38°)(75°, 43°)(16°, 8°)(86°, 53°)(17°, 32°)(53°, 28°)(55°, 60°)(65°, 10°)(14°, 50°)(50°, 25°)(25°, 19.5°)(31°, 48°)(70°, 55°)(69.5°, 60°)(48°, 14°)(27°, 66°)Figure 5
(a) When the proposed grouping scheme in Algorithm2 is adopted, all paths from different UTs can be distinguished at the BS. (b) When the random grouping scheme is adopted, the overlapping paths from different UTs cannot differentiated.
(a)
Using GS(b)
Without using GSIn Figure5, the pentagram is used to depict the estimation on DoA angles, the red circle is used to represent the real DoA angles, and the pink dotted circle is used as an indicator. From Figure 5(a), we can observe that the BS can distinguish all paths from different UTs when the proposed GS in Algorithm 2 is adopted. Even the real DoA angles from the two UTs are very close, they can be differentiated by the BS. On the contrary, in Figure 5(b), the pink circle with the dotted line shows that two overlapped circles only contain one corresponding pentagram, which means the system is failed to distinguish the paths between UT1 and UT5. Therefore, when multiple UTs are present in the system, the GS is necessary for the BS to achieve accurate CSI via the proposed channel estimation and tracking scheme.


### 6.5. Multiple UTs Tracking

Due to the mobility of UTs, tracking the CSI accurately needs to run the GS at the BS periodically. Therefore, the periodicity to execute the GS at the BS is the key to guarantee tracking accuracy. Without loss of generality, the scenario with two UTs is simulated in this section. In the simulation setup, we assume that the signal transmission from each UT to the BS will go through the same path. Moreover, different initial DoA angles and step size are simulated for comparison.Accordingly, the period to execute Algorithm2 is denoted as k time slots. To evaluate the effect of the period on the tracking accuracy, a coefficient, denoted as τ, is defined as(49)τk=∑i=1U∑j=1Jθ^ki,j−θki,j+ϕ^ki,j−ϕki,jJ,whereθ^ki,j and ϕ^ki,j are the estimations of elevation and azimuth angles of the ith UT in the jth time slot, respectively; J is the total running time slots, which is set to J=600. τk implies the average estimation error on the DoA angles on both elevation and azimuth directions.The tracking accuracy with different periods is illustrated in Figures6 and 7, which indicates that optimal GS period is directly related to the DoA angle changing rate and the initial DoA angle difference on elevation direction associated to different UTs. We can observe that the tracking accuracy decreases when the period, k, increases. Optimal period for GS is inversely proportional to the changing rate and proportional to the initial DoA angle difference on elevation direction. Moreover, it is noteworthy that the curve with the diamond marker is stationary even if the difference in the azimuth DoA angles between two UTs is small. Actually, the difference in the azimuth angles has a little effect on the tracking accuracy when the elevation angles of UTs are distinguishable.Figure 6
The impact of changing rate on the tracking accuracy.Figure 7
The impact of initial DoA angle on the tracking accuracy.


## 6.1. The Estimation Accuracy and Convergence Speed

First, to evaluate the estimation accuracy, a Tanimoto coefficient is defined which is given by(46)εk=1−1P∑i=1Pt1,i,k·sit1,i,k2+si2−t1,i,k·si,where si is the ith eigenvector of matrix SP, and t1,i,k is the ith eigenvector of the estimation matrix SP at the kth iteration step. Tanimoto coefficient defined in (46) implies that if ε tends to be 0, the estimation approaches to the real SP. To demonstrate the superiority of the optimal step size derived in Section 3.4, the Backtracking Line Search (BLS) method proposed in [25] is used for comparison. From Figure 2, we can observe that ε achieved by the PCA-based algorithm and BLS algorithm converges to 0, which indicates the availability of these algorithm on estimation of eigenvector matrix. Meanwhile, during the first 400 iterations, the ε achieved by the PCA-based algorithm is less than the BLS algorithm. The result indicates that the proposed scheme has the high accuracy under the same iterations. From another point of view, the proposed scheme has the fast convergence speed on estimation of SP.Figure 2
The Tanimoto coefficient comparison.After recovering the eigenvectors, the intermediate variables,up and vp, are estimated and paired according to the proposed scheme. To demonstrate the match process clearly, the inner products of the received signal matrix on the steering matrix are listed in Table 1, where the steering vectors are constructed by the estimated up and vp.Table 1
Matching onup and vp.
<Aup,vp∣Y>v10.6121v20.3991v30.3363u10.906118.903624.552518.9036u20.707163.173815.03032.6220u30.17416.313419.411422.7863From Table1, we can match up with the right vp by finding the maximum value. Consequently, u1, u2, and u3 should be paired with v2, v1, and v3, respectively. After pairing, the estimation of DoA angles and path amplitudes can be obtained through (36), (37), and (43), respectively, which is listed in Table 2, where the index p is to indicate the pth path.Table 2
Estimation onup, vp, and βθp,ϕp.
p123Realθp25.0000°45.0000°80.0000°Estimation ofθp25.0282°45.0005°79.9737°Realϕp20.0000°30.0000°70.0000°Estimation ofϕp19.3768°30.0451°70.0308°Realβθp,ϕp0.60000.80001.0000Estimation ofβθp,ϕp0.59930.80451.0011From Table2, we can observe that the estimations on both θp, ϕp, and βθp,ϕp are close to the real values, which verifies the effectiveness of the proposed scheme. Moreover, since the ϕp is determined by both up and vp, the estimation error on up and vp would affect the estimation accuracy on ϕp. Therefore, the estimation accuracy of the azimuth angles, ϕp, is less than that on the elevation angles, θp.


## 6.2. The Impact of Received SNR on the DoA Estimation

To reveal estimation performance with different SNRs, the normalized mean square error (MSE) on DoA angle estimation is demonstrated in this section. The MSE is defined as(47)MSEθ=1P∑p=1P∣θ¯p−θp∣θp,(48)MSEϕ=1P∑p=1Pϕ¯p−ϕpϕp.Figure3 illustrates the MSE performance on DoA angle estimation. It can be observed that the estimation MSE on DoA angles reduces as the received SNR at the BS increases. Figure 3 also compares the proposed scheme with the BLS algorithm with respect to DoA angle estimation MSE. We can observe that the proposed scheme outperforms the BLS algorithm on the estimation accuracy. Since the estimation on azimuth DoA angle depends on both the estimated up and vp, the estimation accuracy on the elevation DoA angle is higher than that on the azimuth DoA angle.Figure 3
The MSE on DoA estimation.


## 6.3. Tracking Performance

To verify the tracking ability of the proposed scheme, the SVD-based channel estimation is simulated for comparison, where the SVD operation is executed periodically to update the channel estimation. Assuming that the elevation angle changes from 45° to 53°, and azimuth angles change from 30° to 22°, with the changing rate of 0.02°/stepsize.As shown in Figure4, the SVD-based channel estimation method cannot track the change on the channel even though it may estimate channel with high accuracy in some times. On the other hand, the channel estimation based on Algorithm 1 is very close to the perfect one all the time even when the channel changes. Therefore, this simulation result confirms the tracking ability on the channel estimation of Algorithm 1.Figure 4
Tracking on the DoA angles.


## 6.4. The Effectiveness of Grouping Scheme

For the sake of illustration on the effect of the UT grouping scheme, a multiple UT scenario is considered in this section. The elevation and azimuth angles corresponding to the distinguish paths for 10 different UTs are listed in Table3. Assuming that there are only 2 orthogonal pilot sequences and set λ=5°. First, the UTs are assigned to 5 groups. UTs in 5 groups will reuse the 2 orthogonal pilot sequences in a TDMA manner, and the BS can obtain the initial DoA angle estimations of UTs. Afterwards, based on the grouping scheme (GS) in Algorithm 2, UT1, UT2, UT3, UT4, and UT7 should be put into one group, and the remaining UTs should be put into another group. To better demonstrate the effectiveness of GS, another random grouping scheme is adopted for comparison, where the UT1 to UT5 are put into the same group, and the UT6 to UT10 are put into another group. Accordingly, the estimation on DoA angles is shown in Figure 5.Table 3
Add caption.
The number of UT12345678910The corresponding elevation and azimuth angles (θp,ϕp)(25°, 20°)(45°, 30°)(80°, 70°)(37°, 82°)(60°, 38°)(75°, 43°)(16°, 8°)(86°, 53°)(17°, 32°)(53°, 28°)(55°, 60°)(65°, 10°)(14°, 50°)(50°, 25°)(25°, 19.5°)(31°, 48°)(70°, 55°)(69.5°, 60°)(48°, 14°)(27°, 66°)Figure 5
(a) When the proposed grouping scheme in Algorithm2 is adopted, all paths from different UTs can be distinguished at the BS. (b) When the random grouping scheme is adopted, the overlapping paths from different UTs cannot differentiated.
(a)
Using GS(b)
Without using GSIn Figure5, the pentagram is used to depict the estimation on DoA angles, the red circle is used to represent the real DoA angles, and the pink dotted circle is used as an indicator. From Figure 5(a), we can observe that the BS can distinguish all paths from different UTs when the proposed GS in Algorithm 2 is adopted. Even the real DoA angles from the two UTs are very close, they can be differentiated by the BS. On the contrary, in Figure 5(b), the pink circle with the dotted line shows that two overlapped circles only contain one corresponding pentagram, which means the system is failed to distinguish the paths between UT1 and UT5. Therefore, when multiple UTs are present in the system, the GS is necessary for the BS to achieve accurate CSI via the proposed channel estimation and tracking scheme.


## 6.5. Multiple UTs Tracking

Due to the mobility of UTs, tracking the CSI accurately needs to run the GS at the BS periodically. Therefore, the periodicity to execute the GS at the BS is the key to guarantee tracking accuracy. Without loss of generality, the scenario with two UTs is simulated in this section. In the simulation setup, we assume that the signal transmission from each UT to the BS will go through the same path. Moreover, different initial DoA angles and step size are simulated for comparison.Accordingly, the period to execute Algorithm2 is denoted as k time slots. To evaluate the effect of the period on the tracking accuracy, a coefficient, denoted as τ, is defined as(49)τk=∑i=1U∑j=1Jθ^ki,j−θki,j+ϕ^ki,j−ϕki,jJ,whereθ^ki,j and ϕ^ki,j are the estimations of elevation and azimuth angles of the ith UT in the jth time slot, respectively; J is the total running time slots, which is set to J=600. τk implies the average estimation error on the DoA angles on both elevation and azimuth directions.The tracking accuracy with different periods is illustrated in Figures6 and 7, which indicates that optimal GS period is directly related to the DoA angle changing rate and the initial DoA angle difference on elevation direction associated to different UTs. We can observe that the tracking accuracy decreases when the period, k, increases. Optimal period for GS is inversely proportional to the changing rate and proportional to the initial DoA angle difference on elevation direction. Moreover, it is noteworthy that the curve with the diamond marker is stationary even if the difference in the azimuth DoA angles between two UTs is small. Actually, the difference in the azimuth angles has a little effect on the tracking accuracy when the elevation angles of UTs are distinguishable.Figure 6
The impact of changing rate on the tracking accuracy.Figure 7
The impact of initial DoA angle on the tracking accuracy.


## 7. Conclusion

In this paper, we have proposed a CSI estimation and tracking method for millimeter wave massive MIMO system. We have illustrated the CSI estimation from two aspects: computational complexity and pilot contamination. First, to reduce the computational complexity, a PCA-based algorithm is applied to recover the eigenvectors of the received signal at the BS. Then, with estimated eigenvectors, DoA angles can be estimated by using the MUSIC algorithm and inner products operation. To reduce the pilot overhead, we have utilized the grouping method to divide UTs into different groups. Then, UTs in the same group can share the same pilot sequence and transmit simultaneously without causing interference. In order to express the convergence speed of the proposed algorithm, a BLS algorithm is used for comparison. Compared to this algorithm, the proposed algorithm has faster convergence. Meanwhile, most existing methods are one-time estimation schemes, which need to be rerun after the channel coherence time, but the proposed method possesses good tracking behavior. The numerical simulations are established to verify the superior performance of the proposed method.



---

*Source: 2922359-2021-07-13.xml*