I Introduction
As a promising technique for the fifth generation (5G) cellular network, massive multipleinput multipleoutput (MIMO), a MIMO system with tens or hundreds of antennas, has attracted attentions from both the academia and industry [1][2]. Although massive MIMO is able to maximize its usage of spatial resources due to its high spatial resolution, previous studies suggested that pilot contamination would be a key limiting factor to massive MIMO networks [3], i.e., the spectral efficiency in a massive MIMO network will be bounded even if the number of antenna grows to infinity.
However, recent research progresses in massive MIMO [4] have shown that massive MIMO has unlimited capacity as long as the channels satisfy a certain condition, i.e., the covariance matrices of user equipments (UEs) using the same pilot are asymptotically linearly independent. In [4], the authors claimed that this condition can be satisfied in ordinary systems. Then, a multicell minimum mean squared error (MMMSE) receiver was developed and it was shown that MMMSE can achieve unbounded spectral efficiency.
Additionally, the impact of imperfect covariance matrix in massive MIMO was discussed in [5], where a widely used covariance matrix estimation method known as the viaQ method was proposed. The viaQ method estimated the covariance matrix by calculating the weighted average of the sample covariance matrix and its diagonal. However, the viaQ method requires the true knowledge of the covariance matrix during estimation. Authors in [6] proposed a twostep procedure to reconstruct the covariance matrix. However, this procedure requires proper pilot allocation to UEs, which increases its complexity. In the latest 5G systems, in addition to uniform linear arrays (ULAs), uniform planar arrays (UPAs) are used as well [7]. Performance of these covariance matrix methods for UPAs is unanswered.
Therefore, to fill these research gaps, a novel lowcomplexity and practical covariance matrix estimation method named the antennalayoutaware (ALA) method is proposed in the paper. This method requires the BS to have access to the layout of its antenna array only, which can be conveniently set when a BS is deployed. The contributions of this paper are summarized as follows.

This proposed method has been shown effective for both ULA and UPA layouts, which is compatible to most 5G systems.

Most importantly, the proposed method does not require the true knowledge of the covariance matrix.
Ii System Model
This paper considers a mobile network with base stations (BSs) in a timedivision duplexing (TDD) mode, i.e., channel reciprocity is assumed. Each BS is equipped with antennas and corresponds to a cell. Also, there is one singleantenna UE in each cell. All UEs are assumed to share the same pilot symbols and timefrequency resources. UEs are arranged in a way such that all UEs are at the cell edge of the center cell. Moreover, the target UE lies in the line connecting an interfering UE and a neighbor BS. These settings are as shown in Fig. 1. As a result, the worst scenario, i.e., BSs experiencing the largest interference in uplink and UEs experiencing the largest interference in downlink, is considered in this paper.
Let
be the channel vector between the
BS and the th UE () with , where is the corresponding spatial covariance matrix. Considering the Kronecker channel model and UPAs, can be presented as(1) 
where denotes the Kronecker product, is the elevation covariance matrix, and is the azimuth covariance matrix. The exponential spatial correlation model [4][8] is used in this paper. Therefore, the element in the th row and the th column of is shown as
(2) 
where is the correlation factor and is the angle of arrival (AoA) in the azimuth (elevation). In the special case of a ULA, equals .
Assuming that the signaltonoise ratio (SNR) between the th BS and the th UE is , the decorrelated received signal vector is presented as
(3) 
where is the Gaussian noise in the uplink following . Then, assuming that both channel estimation and data transmission phases are within the coherence time of the channel, the minimum mean squared error (MMSE) estimate of the channel vector with ideal covariance matrix can be expressed as [9]
(4) 
where
is the sum of covariance matrices of users using the same pilot and a scaled identity matrix, i.e.,
(5) 
Using the orthogonality property of MMSE, the estimated channel vector is distributed as , with covariance matrix
(6) 
In practice, the assumption that BSs have knowledge of ideal covariance matrices is not realistic. The estimated covariance matrices and can only be obtained via processing the sample covariance matrices and . The sample covariance matrices are calculated by
(7) 
and
(8) 
where is the number of pilot symbols and is the received signal excluding the associated UE’s contribution. This can be done by using two slots for channel estimation. In the first slot, all UEs in the network will transmit pilot symbols simultaneously. In the second slot, only UEs in the neighboring cells will transmit pilot symbols. The processing function is an algorithm operating on the sample covariance matrix, which will be the viaQ method [5] or the ALA method in this paper.
After estimating the channel in the uplink, a BS will transmit information to its UE in downlink. Meanwhile, a UE will receive both the useful signal from its associated BS and interference from neighboring BSs. Let denote the downlink SNR between the th BS and the th UE, the received signal of the th UE can be expressed as
(9) 
where
is the zero mean unit variance Gaussian noise and
is the precoding vector of the th BS.Here, we follow the downlink precoder design in [10], where the sum of the signal detection error and signal leakage is minimized, i.e.,
(10) 
This optimization can be solved in closed form as in [10] and the optimal precoding vector can be computed as
(11) 
where is a normalization factor for the precoding vector. The downlink signaltonoiseplusinterference ratio (SINR) of the th UE with the optimal precoding vector can be expressed as
(12) 
where is the variance operator. It should be mentioned that the effects of nonideal channel estimation and spatial covariance matrix estimation have been factored in (12).
Iii ALA Covariance Matrix Estimation
After obtaining and , further processing can be performed to calculate the final estimated covariance matrices and to improve the channel estimation performance. The viaQ method proposed in [5][11] computes as a weighted average of and its diagonal, i.e.,
(13) 
It was shown in [5] and [11] that the optimal weight (regularization factor) can be computed in closed form. However, during the calculation of the optimal weight, true knowledge of and is required by the BS. This requirement is not practical.
On the contrary, the ALA method proposed in this paper does not require the true knowledge of any covariance matrices. The only additional information needed in the estimation process is the antenna layout. This can conveniently be accessed when a BS is deployed. Moveover, unlike the viaQ method, which assumes full degrees of freedom in the estimation process, the proposed ALA method takes the advantage that a covariance matrix does not have full degrees of freedom. This means that if certain antenna pairs whose relative positions are the same, then their covariance values are statistically equivalent.
As both ULA and UPA are the most widely used antenna layouts, this paper focuses on discussing the ALA covariance matrix estimation for these two layouts.
Iiia Ula
For a antenna ULA, its antennas are placed in a horizontal line. Moreover, since the antenna elements are equally spaced, the covariance matrix of a ULA is a Toeplitz matrix. In this case, elements in an offdiagonal are statistically equivalent. Let us consider two pairs of antennas, PAIR() and PAIR() for instance, if the two interelement spacings are the same (), the two covariance values will be equal, i.e., . Therefore, all elements in the covariance matrix satisfying the above condition can be replaced by the average of them. Using the Toeplitz structure of a ULA covariance matrix, this is equivalent to computing the average value along offdiagonal lines, i.e.,
(14) 
Pseudo codes of the ALA covariance matrix estimation for ULA are shown in Fig. 2. The same procedure can be used to estimate .
IiiB Upa
For a antenna UPA, its antennas are arranged in a plane in a columnmajor manner. The location of each antenna () can be represented as in an by grid, where is the number of antennas in each column and is the number of antennas in each row, with and . An example of a antenna array arranged as a by panel is depicted in Fig. 3. Given the antenna index , the coordinates can be calculated as
where is the maximum integer less than and is the modulo operation with respect to . On the contrary, when the coordinates of an antenna are given, the antenna index can be computed as .
To describe the principles of ALA covariance matrix estimation for UPA, let us consider an antenna pair PAIR. It can be noticed that, when an angular spectrum is provided, the covariance value of PAIR depends only on the relative positions of antennas and . As a result, the covariance value of an antenna pair is translationinvariant in the 2D grid. An example is illustrated in Fig. 3. Let us consider PAIR() and PAIR(). The antenna relative positions in these two pairs are the same and their covariance values are statistically equivalent, i.e., . However, covariance values are not rotationinvariant to antenna pairs, because the angular spread in horizontal may not be the same as that in elevation. Rotating the antenna pair will change its covariance value. Then, the estimated covariance value can be calculated as the average of all covariance values of antenna pairs with the same relative positions.
Next, assuming a given PAIR(), the remaining question is how to identify the set consisting of all antenna pairs with the same relative antenna positions as PAIR(). PAIR() is a translation of PAIR() if the coordinates in PAIR() can be expressed as shifted versions of coordinates in PAIR(). These can be presented as
(15) 
(16) 
Considering that the ranges of the horizontal direction and the vertical direction are bounded by and , the set can be expressed as
(17) 
Moreover, the cardinality can be calculated as
(18) 
As a result, elements in the estimated covariance matrix using ALA can be computed as
(19) 
The pseudo codes of the ALA covariance matrix estimation for UPA are shown in Fig. 4. It can be seen in Fig. 4 that a by boolean assistance matrix is used to mark which elements have been calculated. If an element has been already calculated, the loop will be skipped. This assistance matrix can be used to avoid repeated computations, minimizing the complexity.
IiiC General antenna layouts
Although the most two typical antenna layouts, i.e., ULA and UPA, are discussed in this paper, the proposed ALA covariance matrix estimation algorithm can be applied to more general antenna layouts as long as the following two conditions are satisfied. First, the antenna layout information needs to be provided to the estimator as prior information. Second, only antenna pairs with translation operations can contribute when calculating the average covariance value.
IiiD Complexity analysis of the ALA estimation method
Iv Results and Analysis
Simulations are performed in a onering sevencell network. Uplink power settings are assumed to be the same as [5], i.e., if and if . In downlink, the SNR of the attached UE is dB and the SNRs of UEs in neighbor cells are calculated proportionally to their distances to the the BS assuming the path loss exponent is . For example, in Fig. 1, let us assume that the distance between the BS and its attached UE (center UE) is and the downlink SNR is dB. Then, the distance between BS2 and the center UE is . With path loss exponent , it can be obtained that the SNR between BS2 and the center UE is dB. All other SNRs are computed in the same way. Moreover, correlation factors are set as and for azimuth and elevation. Azimuth AoAs
are drawn from a uniform distribution within
and elevation AoAs are drawn from a uniform distribution within .The normalized channel estimation errors with different covariance matrix estimation methods, number of pilot samples, and antenna layouts are depicted in Fig. 5. The normalized channel estimation error is measured by the mean squared error (MSE) and is calculated as [5]
(20) 
The MSE curves with ideal covariance matrix serve as lower bounds and remain flat in terms of as expected. When is relatively small (), the viaQ method outperforms the proposed ALA method. The reason for this is that multiple entries in the estimated covariance matrix of the ALA method share a common value. When this common value is not accurate enough due to insufficient number of pilots, the estimated error propagates, causing the slightly poorer performance in the ALA method when is relatively small. However, when continues to grow, the proposed ALA method dominates and is able to close its gap to the ideal covariance matrix curve more quickly than the viaQ method. The viaQ method suffers from its nature in (13) that its offdiagonal entries are always scaled and times smaller than those in the sample covariance matrix. The impact of these scaled entries is shown to be significant when is sufficiently large.
Downlink spectral efficiency with ULA is illustrated in Fig. 6. The spectral efficiency with ideal covariance matrix serves as the upper bound and grows with the antenna number. When the antenna number is small, the ALA and viaQ achieve similar performance. However, when the antenna number is moderate and large , the ALA method achieves better spectral efficiency than the viaQ method. It can also be observed that the spectral efficiency of the viaQ method starts to saturate, while the spectral efficiency of the ALA method continues to increase with antenna number.
Downlink spectral efficiency with UPA is shown in Fig. 7. Spectral efficiencies of both the ideal and ALA method show a similar trend as those using ULA. What is more, the spectral efficiency of the viaQ method is able to grow with the antenna number as well. This shows that UPA is more resilient to covariance matrix estimation errors than ULA. The reason for this can be explained in Table I. It can be seen that the viaQ for ULA converges to after 128 antennas, meaning that the viaQ estimated covariance matrix is fully determined by the diagonal matrix of the sample covariance matrix. On the other hand, the viaQ for UPA converges much slower. As a result, the offdiagonal entries of the sample covariance matrix still contribute to the estimated covariance matrix.
4  8  16  32  64  128  256  

ULA  0.29  0.38  0.54  0.70  0.84  0.96  1.00  1.00 
UPA  0.29  0.26  0.37  0.45  0.60  0.71  0.84  0.93 
V Conclusions
A lowcomplexity ALA covariance matrix estimation method has been presented in this paper. Since covariance matrices do not have full degrees of freedom, the proposed ALA method can maximize the benefit of this property by allowing a BS to have knowledge of its antenna layout. It has been shown that the ALA method is effective to both ULA and UPA layouts, which are the two most widely used in practice. The proposed ALA method has been applied to a multicell network. Simulations have shown that the proposed ALA method has lower MSE of the estimated channel than the viaQ method when the number of pilot symbols is moderate. Additionally, the proposed ALA method significantly outperforms the viaQ method in terms of spectral efficiency. For future work, the application of the ALA method to multicell multiuser network in both uplink and downlink can be investigated. Moreover, impact of effects such as antenna rotations, misalignment, and imperfect knowledge of antenna geometry can be studied.
Acknowledgment
Part of this work has been performed in the framework of the Horizon 2020 project ONE5G (ICT760809) receiving funds from the European Union. The authors would like to acknowledge the contributions of their colleagues in the project, although the views expressed in this contribution are those of the authors and do not necessarily represent the project.
References
 [1] J. G. Andrews, S. Buzzi, W. Choi, S. V. Hanly, A. Lozano, A. C. K. Soong, and J. C. Zhang, “What will 5G be?,” IEEE J. Sel. Areas Commun., vol. 32, no. 6, pp. 1065–1082, June 2014.
 [2] A. Zaidi, R. Baldemair, M. Andersson, S. Faxer, V. Molescases, and Z. Wang, “5G new radio: designing for the future,” Ericsson Technol. Rev., 1–14, June 2017.
 [3] T. L. Marzetta, “Noncooperative cellular wireless with unlimited numbers of BS antennas,” IEEE Trans. Wireless Commun., vol. 9, no. 11, pp. 3590–3600, Nov. 2010.
 [4] E. Bjrnson, J. Hoydis, and L. Sanguinetti, “Massive MIMO has unlimited capacity,” IEEE Trans. Wireless Commun., vol. 17, no. 1, pp. 574–590, Nov. 2018.
 [5] E. Bjrnson, L. Sanguinetti and M. Debbah, “Massive MIMO with imperfect channel covariance information,” in Proc. ACSSC’16., Pacific Grove, USA, Nov. 2016, pp. 974–978.
 [6] D. Neumann, M. Joham, and W. Utschick, “Covariance matrix estimation in massive MIMO,” IEEE Signal Process. Lett., vol. 25, no. 6, pp. 863–867, June 2018.
 [7] 3GPP T. R. 38.901, “Study on channel model for frequencies from 0.5 to 100 GHz ,” V14.0.0, Mar. 2017.
 [8] S. L. Loyka, “Channel capacity of MIMO architecture using the exponential correlation matrix,” IEEE Commun. Lett., vol. 5, no. 9, pp. 369–371, Sep. 2001.
 [9] T. Kailath, A. H. Sayed, and B. Hassibi, Linear Estimation., Prentice Hall, New Jersey, 2000.
 [10] J. Jose, A. Ashikhmin, T. L. Marzetta, and S. Vishwanath, “Pilot contamination and precoding in multicell TDD systems,” IEEE Trans. Wireless Commun., vol. 10, no. 8, pp. 2640–2651, Aug. 2011.
 [11] N. Shariati, E. Bjrnson, M. Bengtsson, and M. Debbah, “Lowcomplexity polynomial channel estimation in largescale MIMO with arbitrary statistics,” IEEE J. Sel. Topics Signal Proc., vol. 8, no. 5, pp. 815–830, Oct. 2014.