Parametric channel estimation for massive MIMO
Abstract
Channel state information is crucial to achieving the capacity of multiantenna (MIMO) wireless communication systems. It requires estimating the channel matrix. This estimation task is studied, considering a sparse physical channel model, as well as a general measurement model taking into account hybrid architectures. The contribution is twofold. First, the CramérRao bound in this context is derived. Second, interpretation of the Fisher Information Matrix structure allows to assess the role of system parameters, as well as to propose asymptotically optimal and computationally efficient estimation algorithms.
Parametric channel estimation for massive MIMO
Luc Le Magoarou, Stéphane Paquelet ^{†}^{†}thanks: 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. 
b com, Rennes, France 
Index Terms— CramérRao bound, Channel estimation, MIMO.
1 Introduction
MultipleInput MultipleOutput (MIMO) wireless communication systems allow for a dramatic increase in channel capacity, by adding the spatial dimension to the classical time and frequency ones [1, 2]. This is done by sampling space with several antenna elements, forming antenna arrays both at the transmitter (with antennas) and receiver (with antennas). Capacity gains over single antenna systems are at most proportional to .
Millimeter wavelengths have recently appeared as a viable solution for the fifth generation (5G) wireless communication systems [3, 4]. Indeed, smaller wavelengths allow to densify halfwavelength separated antennas, resulting in higher angular resolution and capacity for a given array size. This observation has given rise to the massive MIMO field, i.e. the study of systems with up to hundreds or even thousands of antennas.
Massive MIMO systems are very promising in terms of capacity. However, they pose several challenges to the research community [5, 6], in particular for channel estimation. Indeed, maximal capacity gains are obtained in the case of perfect knowledge of the channel state by both the transmitter and the receiver. The estimation task amounts to determine a complex gain between each transmit/receive antenna pair, the narrowband (single carrier) MIMO channel as a whole being usually represented as a complex matrix of such complex gains. Without a parametric model, the number of real parameters to estimate is thus , which is very large for massive MIMO systems.
Contributions and organization. In this work, massive MIMO channel estimation is studied, and its performance limits are sought, as well as their dependency on key system parameters. In order to answer this question, the framework of parametric estimation [7] is used. A physical channel model is first presented, with the general considered observation model, and the objective is precisely stated. The CramérRao bound for is then derived, which bounds the variance of any unbiased estimator. Then, the interpretation of the bound allows to precisely assess the role of system design on estimation performance, as well as to propose new computationally efficient channel estimation algorithms showing asymptotic performance equivalent to classical ones based on sparse recovery.
2 Problem formulation
Notations. Matrices and vectors are denoted by bold uppercase and lowercase letters: and (except 3D “spatial” vectors that are denoted ); the th column of a matrix by: ; its entry at the th line and th column by: or . A matrix transpose, conjugate and transconjugate is denoted by: , and respectively. The image, rank and trace of a linear transformation represented by are denoted: , and respectively. For matrices and , means that is positive semidefinite. The linear span of a set of vectors is denoted: . The Kronecker product, standard vectorization and diagonalization operators are denoted by , , and respectively. The identity matrix, the matrix of zeros and ones are denoted by , and respectively. denotes the standard complex gaussian distribution with mean and covariance . denotes expectation and the covariance of its argument.
2.1 Parametric physical channel model
Consider a narrowband block fading channel between a transmitter and a receiver with respectively and antennas. It is represented by the matrix , in which corresponds to the channel between the th transmit and th receive antennas.
Classically, for MIMO systems with few antennas, i.e. when the quantity is small (up to a few dozens), estimators such as the Least Squares (LS) or the Linear Minimum Mean Squared Error (LMMSE) are used [8].
However, for massive MIMO systems, the quantity is large (typically several hundreds), and resorting to classical estimators may become computationally intractable. In that case, a parametric model may be used. Establishing it consists in defining a set of parameters that describe the channel as for a given function , where the approximation is inherent to the model structure and neglected in the sequel (considering ). Channel estimation then amounts to estimate the parameters instead of the channel matrix directly. The parametrization is particularly useful if , without harming accuracy of the channel description. Inspired by the physics of wave propagation under the plane waves assumption, it has been proposed to express the channel matrix as a sum of rank1 matrices, each corresponding to a single physical path between transmitter and receiver [9]. Adopting this kind of modeling and generalizing it to take into account any threedimensional antenna array geometry, channel matrices take the form
(1) 
where is the total number of considered paths (no more than a few dozens), is the complex gain of the th path, is the unit vector corresponding to its Direction of Departure (DoD) and the unit vector corresponding to its Direction of Arrival (DoA). Any unit vector is described in spherical coordinates by an azimuth angle and an elevation angle . The complex response and steering vectors and are defined as for . The set gathers the positions of the antennas with respect to the centroid of the considered array (transmit if , receive if ). In order to lighten notations, the matrix is introduced. It simplifies the steering/response vector expression to , where the exponential function is applied componentwise. In order to further lighten notations, the th atomic channel is defined as , and its vectorized version . Therefore, defining the vectorized channel , yields . Note that the channel description used here is very general, as it handles any threedimensional antenna array geometry, and not only Uniform Linear Arrays (ULA) or Uniform Planar Arrays (UPA) as is sometimes proposed.
In short, the physical channel model can be seen as a parametric model with . There are thus real parameters in this model (the complex gain, DoD and DoA of every path are described with two parameters each). Of course, the model is most useful for estimation in the case where , since the number of parameters is thus greatly reduced.
Note that most classical massive MIMO channel estimation methods assume a similar physical model, but discretize a priori the DoDs and DoAs, so that the problem fits the framework of sparse recovery [10, 11, 12]. The approach used here is different, in the sense that no discretization is assumed for the analysis.
2.2 Observation model
In order to carry out channel estimation, known pilot symbols are sent through the channel by each transmit antenna. The corresponding training matrix is denoted . The signal at the receive antennas is thus expressed as , where is a noise matrix with . Due to the high cost and power consumption of millimeter wave Radio Frequency (RF) chains, it has been proposed to have less RF chains than antennas in both the transmitter and receiver [13, 14, 15, 16]. Such systems are often referred to as hybrid architectures. Mathematically speaking, this translates into specific constraints on the training matrix (which has to “sense” the channel through analog precoders , , being the number of RF chains on the transmit side), as well as observing the signal at the receiver through analog combiners. Let us denote , the used analog combiners, the observed data is thus expressed in all generality as
(2) 
where and the training matrix is constrained to be of the form , where is the digital training matrix.
2.3 Objective: bounding the variance of unbiased estimators
In order to assess the fundamental performance limits of channel estimation, the considered performance measure is the relative Mean Squared Error (rMSE). Denoting indifferently or the true channel ( or in vectorized form) and or its estimate ( or in vectorized form) in order to lighten notations, rMSE is expressed
(3) 
where the bias/variance decomposition can be done independently of the considered model [7]. The goal here is to lowerbound the variance term, considering the physical model introduced in the previous subsection. The bias term is not studied in details here, but its role is evoked in section 3.3.
3 CramérRao lower bound
In this section, the variance term of eq. (3) is bounded using the CramérRao Bound (CRB) [17, 18], which is valid for any unbiased estimator of the true parameter . The complex CRB [19] states,
with the Fisher Information Matrix (FIM), where denotes the model likelihood, and is any complex differentiable vector function. In particular, regarding the variance term of eq. (3),
(4) 
with A model independent expression for the FIM is provided in section 3.1, and particularized in section 3.2 to the model of section: 2.1. Finally, the bound is derived from eq. (4) in section 3.3.
3.1 General derivation
First, notice that vectorizing eq. (2), the observation matrix follows a complex gaussian distribution,
In that particular case, the SlepianBangs formula [20, 21] yields:
(5) 
with where is the average transmit power per time step. Note that the expression can be simplified to using elementary properties of the Kronecker product. The matrix is a projection matrix onto the range of . In order to ease further interpretation, assume that . This assumption means that the transmit power is constant during training time (, ) and that pilots sent at different time instants are mutually orthogonal (, ). This way, is a projection matrix onto the range of , and can itself be interpreted as a projection, being the Kronecker product of two projection matrices [22, p.112] (it is an orthogonal projection since ).
3.2 Fisher information matrix for a sparse channel model
Consider now the parametric channel model of section 2.1, where , with .
Intrapath couplings. The derivatives of with respect to parameters of the th path can be determined using matrix differentiation rules [23]:

Regarding the complex gain , the model yields the expressions and .

Regarding the DoA, and , where and are the unit vectors in the azimuth and elevation directions at , respectively.

Regarding the DoD, and , where and are the unit vectors in the azimuth and elevation directions at , respectively.
Denoting , the part of the FIM corresponding to couplings between the parameters (intrapath couplings) is expressed as
(6) 
Let us now particularize this expression. First of all, in order to ease interpretations carried out in section 4, consider the case of optimal observation conditions (when the range of contains the range of ). This allows indeed to interpret separately the role of the observation matrices and the antenna arrays geometries. Second, consider for example the entry corresponding to the coupling between the departure azimuth angle and the arrival azimuth angle of the th path. It is expressed under the optimal observation assumption as . Moreover,
since and by construction (because the antennas positions are taken with respect to the array centroid). This means that the parameters and are statistically uncoupled, i.e. orthogonal parameters [24]. Computing all couplings for yields
(7) 
where
(8) 
with . These expressions are thoroughly interpreted in section 4.
Global FIM. Taking into account couplings between all paths, The global FIM is easily deduced from the previous calculations and block structured,
where contains the couplings between parameters of the th and th paths and is expressed . The offdiagonal blocks of , corresponding to couplings between parameters of distinct paths, or interpath couplings, can be expressed explicitly (as in eq. (7) for intrapath couplings). However, the obtained expressions are less prone to interesting interpretations, and interpaths couplings have been observed to be negligible in most cases. They are thus not displayed in the present paper, for brevity reasons. Note that a similar FIM computation was recently carried out in the particular case of linear arrays [25]. However, the form of the FIM (in particular parameter orthogonality) was not exploited in [25], as is done here in sections 4 and 5.
3.3 Bound on the variance
The variance of channel estimators remains to be bounded, using eq. (4). From eq. (5), the FIM can be expressed more conveniently only with real matrices as with
where is also a projection matrix. Finally, injecting eq. (5) into eq. (4) assuming the FIM is invertible, gives for the relative variance
(9) 
where the second inequality comes from the fact that being an orthogonal projection matrix, (using elementary properties of the ordering of semidefinite positive matrices, in particular [26, Theorem 4.3]). The first equality comes from the fact that . Finally, the second equality is justified by considering the sparse channel model, and by taking (this is actually an optimal SNR, only attained with perfect precoding and combining).
Optimal bound. The first inequality in eq. (9) becomes an equality if an efficient estimator is used [7]. Moreover, the second inequality is an equality if the condition is fulfilled (this corresponds to optimal observations, further discussed in section 4). Remarkably, under optimal observations, the lower bound on the relative variance is directly proportional to the considered number of paths and inversely proportional to the SNR, and does not depend on the specific model structure, since the influence of the derivative matrix cancels out in the derivation.
Sparse recovery CRB. It is interesting to notice that the bound obtained here is similar to the CRB for sparse recovery [27] (corresponding to an intrinsically discrete model), that is proportional to the sparsity of the estimated vector, analogous here to the number of paths.
4 Interpretations
The main results of sections 3.2 and 3.3 are interpreted in this section, ultimately guiding the design of efficient estimation algorithms.
Parameterization choice. The particular expression of the FIM allows to assess precisely the chosen parameterization. First of all, has to be invertible and wellconditioned, for the model to be theoretically and practically identifiable [28, 29], respectively. As a counterexample, imagine two paths indexed by and share the same DoD and DoA, then proportional columns appear in , which implies noninvertibility of the FIM. However, it is possible to summarize the effect of these two paths with a single virtual path of complex gain without any accuracy loss in channel description, yielding an invertible FIM. Similarly, two paths with very close DoD and DoA yield an illconditioned FIM (since the corresponding steering vectors are close to colinear), but can be merged into a single virtual path with a limited accuracy loss, improving the conditioning. Interestingly, in most channel models, paths are assumed to be grouped into clusters, in which all DoDs and DoAs are close to a principal direction [30, 31, 32]. Considering the MSE, merging close paths indeed decreases the variance term (lowering the total number of parameters), without increasing significantly the bias term (because their effects on the channel matrix are very correlated). These considerations suggest dissociating the number of paths considered in the model from the number of physical paths, denoted , taking by merging paths. This is one motivation behind the famous virtual channel representation [9], where the resolution at which paths are merged is fixed and given by the number of antennas. The theoretical framework of this paper suggests to set (and thus the merging resolution) so as to minimize the MSE. A theoretical study of the bias term of the MSE (which should decrease when increases) could thus allow to calibrate models, choosing an optimal number of paths for estimation. Such a quest for is carried out empirically in section 5.
Optimal observations. The matrices and (pilot symbols and analog combiners) determine the quality of channel observation. Indeed, it was shown in section 3.3 that the lowest CRB is obtained when , with . In case of sparse channel model, using the expressions for derived above, this is equivalent to two distinct conditions for the training matrix:
and for the analog combiners:
where with and . These conditions are fairly intuitive: to estimate accurately parameters corresponding to a given DoD (respectively DoA), the sent pilot sequence (respectively analog combiners) should span the corresponding steering vector and its derivatives (to “sense” small changes). To accurately estimate all the channel parameters, it should be met for each atomic channel.
Array geometry. Under optimal observation conditions, performance limits on DoD/DoA estimation are given by eq. (8). The lower the diagonal entries , the better the bound. This implies the bound is better if the diagonal entries of are large and the offdiagonal entries are small (in absolute value). Since the unit vectors and are by definition orthogonal, having with maximal is optimal, and yields uniform performance limits for any DoD/DoA. Moreover, in this situation, is proportional to , the mean squared norm of antenna positions with respect to the array centroid. Having a larger antenna array is thus beneficial (as expected), because the furthest antennas are from the array centroid, the larger is.
Orthogonality of DoA and DoD. Section 3.2 shows that the matrix corresponding to intrapath couplings (eq. (7)) is block diagonal, meaning that for a given path, parameters corresponding to gain, phase, DoD and DoA are mutually orthogonal. Maximum Likelihood (ML) estimators of orthogonal parameters are asymptotically independent [24] (when the number of observations, or equivalently the SNR goes to infinity). Classically, channel estimation in massive MIMO systems is done using greedy sparse recovery algorithms [10, 11, 12]. Such algorithms can be cast into ML estimation with discretized directions, in which the DoD and DoA (coefficient support) are estimated jointly first (which is costly), and then the gain and phase are deduced (coefficient value), iteratively for each path. Orthogonality between the DoD and DoA parameters is thus not exploited by classical channel estimation methods. We propose here to exploit it via a sequential decoupled DoD/DoA estimation, that can be inserted in any sparse recovery algorithm in place of the support estimation step, without loss of optimality in the ML sense. In the proposed method, one direction (DoD or DoA) is estimated first using an ML criterion considering the other direction as a nuisance parameter, and the other one is deduced using the joint ML criterion. Such a strategy is presented in algorithm 1. It can be verified that lines and of the algorithm actually correspond to ML estimation of the DoA and joint ML estimation, respectively. The overall complexity of the sequential directions estimation is thus , compared to for the joint estimation with the same test directions. Note that a similar approach, in which DoAs for all paths are estimated at once first, was recently proposed [33] (without theoretical justification).
5 Preliminary experiment
Let us compare the proposed sequential direction estimation to the classical joint estimation. This experiment must be seen as an example illustrating the potential of the approach, and not as an extensive experimental validation.
Experimental settings. Consider synthetic channels generated using the NYUSIM channel simulator [34] (setting GHz, the distance between transmitter and receiver to m) to obtain the DoDs, DoAs, gains and phases of each path. The channel matrix is then obtained from eq. (1), considering square Uniform Planar Arrays (UPAs) with halfwavelength separated antennas, with and . Optimal observations are considered, taking both and as the identity. Moreover, the noise variance is set so as to get an SNR of dB. Finally, the two aforementioned direction estimation strategies are inserted in the Matching Pursuit (MP) algorithm [10], discretizing the directions taking , and varying the total number of estimated paths.
Joint estimation  Sequential estimation  

rMSE  Time  rMSE  Time  
0.077  1.24  0.092  0.11  
0.031  2.40  0.039  0.16  
0.017  4.66  0.021  0.24  
0.025  9.50  0.023  0.42 
Results. Table 1 shows the obtained relative MSE and estimation times (Python implementation on a laptop with an Intel(R) Core(TM) i73740QM CPU @ 2.70 GHz). First of all, for , the estimation error decreases and the estimation time increases with , exhibiting a tradeoff between accuracy and time. However, increasing beyond a certain point seems useless, since the error reincreases, as shown by the MSE for , echoing the tradeoff evoked in section 3.3, and indicating that is certainly between and for both methods in this setting. Finally, for any value of , while the relative errors of the sequential and joint estimation methods are very similar, the estimation time is much lower (between ten and twenty times) for sequential estimation. This observation validates experimentally the theoretical claims made in the previous section.
6 Conclusions and perspectives
In this paper, the performance limits of massive MIMO channel estimation were studied. To this end, training based estimation with a physical channel model and an hybrid architecture was considered. The Fisher Information Matrix and the CramérRao bound were derived, yielding several results. The CRB ended up being proportional to the number of parameters in the model and independent from the precise model structure. The FIM allowed to draw several conclusions regarding the observation matrices and the arrays geometries. Moreover, it suggested computationally efficient algorithm which are asymptotically as accurate as classical ones.
This paper is obviously only a first step toward a deep theoretical understanding of massive MIMO channel estimation. Apart from more extensive experimental evaluations and optimized algorithms, a theoretical study of the bias term of the MSE would be needed to calibrate models, and the interpretations of section 4 could be leveraged to guide system design.
Acknowledgments. The authors wish to thank Matthieu Crussière for the fruitful discussions that greatly helped improving this work.
References
 [1] Emre Telatar, “Capacity of multiantenna gaussian channels,” European transactions on telecommunications, vol. 10, no. 6, pp. 585–595, 1999.
 [2] David Tse and Pramod Viswanath, Fundamentals of wireless communication, Cambridge university press, 2005.
 [3] Theodore S Rappaport, Shu Sun, Rimma Mayzus, Hang Zhao, Yaniv Azar, Kevin Wang, George N Wong, Jocelyn K Schulz, Mathew Samimi, and Felix Gutierrez, “Millimeter wave mobile communications for 5g cellular: It will work!,” IEEE access, vol. 1, pp. 335–349, 2013.
 [4] A Lee Swindlehurst, Ender Ayanoglu, Payam Heydari, and Filippo Capolino, “Millimeterwave massive mimo: the next wireless revolution?,” IEEE Communications Magazine, vol. 52, no. 9, pp. 56–62, 2014.
 [5] Fredrik Rusek, Daniel Persson, Buon Kiong Lau, Erik G Larsson, Thomas L Marzetta, Ove Edfors, and Fredrik Tufvesson, “Scaling up mimo: Opportunities and challenges with very large arrays,” IEEE Signal Processing Magazine, vol. 30, no. 1, pp. 40–60, 2013.
 [6] Erik G Larsson, Ove Edfors, Fredrik Tufvesson, and Thomas L Marzetta, “Massive mimo for next generation wireless systems,” IEEE Communications Magazine, vol. 52, no. 2, pp. 186–195, 2014.
 [7] Steven M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory, PrenticeHall, Inc., Upper Saddle River, NJ, USA, 1993.
 [8] Mehrzad Biguesh and Alex B Gershman, “Trainingbased mimo channel estimation: a study of estimator tradeoffs and optimal training signals,” IEEE transactions on signal processing, vol. 54, no. 3, pp. 884–893, 2006.
 [9] Akbar M Sayeed, “Deconstructing multiantenna fading channels,” IEEE Transactions on Signal Processing, vol. 50, no. 10, pp. 2563–2579, 2002.
 [10] S.G. Mallat and Z. Zhang, “Matching pursuits with timefrequency dictionaries,” Signal Processing, IEEE Transactions on, vol. 41, no. 12, pp. 3397–3415, Dec 1993.
 [11] J.A. Tropp and A.C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” Information Theory, IEEE Transactions on, vol. 53, no. 12, pp. 4655–4666, Dec 2007.
 [12] W. U. Bajwa, J. Haupt, A. M. Sayeed, and R. Nowak, “Compressed channel sensing: A new approach to estimating sparse multipath channels,” Proceedings of the IEEE, vol. 98, no. 6, pp. 1058–1076, June 2010.
 [13] Omar El Ayach, Sridhar Rajagopal, Shadi AbuSurra, Zhouyue Pi, and Robert W Heath, “Spatially sparse precoding in millimeter wave mimo systems,” IEEE Transactions on Wireless Communications, vol. 13, no. 3, pp. 1499–1513, 2014.
 [14] Ahmed Alkhateeb, Omar El Ayach, Geert Leus, and Robert W Heath, “Channel estimation and hybrid precoding for millimeter wave cellular systems,” IEEE Journal of Selected Topics in Signal Processing, vol. 8, no. 5, pp. 831–846, 2014.
 [15] Robert W Heath, Nuria GonzalezPrelcic, Sundeep Rangan, Wonil Roh, and Akbar M Sayeed, “An overview of signal processing techniques for millimeter wave mimo systems,” IEEE journal of selected topics in signal processing, vol. 10, no. 3, pp. 436–453, 2016.
 [16] Akbar M. Sayeed and John H. Brady, MillimeterWave MIMO Transceivers: Theory, Design and Implementation, pp. 231–253, John Wiley & Sons, Ltd, 2016.
 [17] Calyampudi Radakrishna Rao, “Information and the accuracy attainable in the estimation of statistical parameters,” Bulletin of the Calcutta Mathematical Society, vol. 37, pp. 81–89, 1945.
 [18] Harald Cramér, Mathematical Methods of Statistics, vol. 9, Princeton university press, 1946.
 [19] Adriaan Van den Bos, “A cramérrao lower bound for complex parameters,” IEEE Transactions on Signal Processing [see also Acoustics, Speech, and Signal Processing, IEEE Transactions on], 42 (10), 1994.
 [20] David Slepian, “Estimation of signal parameters in the presence of noise,” Transactions of the IRE Professional Group on Information Theory, vol. 3, no. 3, pp. 68–89, 1954.
 [21] G. W. Bangs, Array Processing With Generalized Beamformers, Ph.D. thesis, Yale university, CT, USA, 1971.
 [22] WilliHans Steeb and Yorick Hardy, Matrix calculus and Kronecker product: a practical approach to linear and multilinear algebra, World Scientific, 2011.
 [23] Kaare Brandt Petersen, Michael Syskind Pedersen, et al., “The matrix cookbook,” Technical University of Denmark, vol. 7, pp. 15, 2008.
 [24] David Roxbee Cox and Nancy Reid, “Parameter orthogonality and approximate conditional inference,” Journal of the Royal Statistical Society. Series B (Methodological), pp. 1–39, 1987.
 [25] Nil Garcia, Henk Wymeersch, and Dirk Slock, “Optimal robust precoders for tracking the aod and aoa of a mmwave path,” arXiv preprint arXiv:1703.10978, 2017.
 [26] Jerzy K Baksalary, Friedrich Pukelsheim, and George PH Styan, “Some properties of matrix partial orderings,” Linear Algebra and its Applications, vol. 119, pp. 57–85, 1989.
 [27] Zvika BenHaim and Yonina C Eldar, “The cramérrao bound for estimating a sparse parameter vector,” IEEE Transactions on Signal Processing, vol. 58, no. 6, pp. 3384–3389, 2010.
 [28] Thomas J Rothenberg, “Identification in parametric models,” Econometrica: Journal of the Econometric Society, pp. 577–591, 1971.
 [29] Costas Kravaris, Juergen Hahn, and Yunfei Chu, “Advances and selected recent developments in state and parameter estimation,” Computers & chemical engineering, vol. 51, pp. 111–123, 2013.
 [30] Adel AM Saleh and Reinaldo Valenzuela, “A statistical model for indoor multipath propagation,” IEEE Journal on selected areas in communications, vol. 5, no. 2, pp. 128–137, 1987.
 [31] Jon W Wallace and Michael A Jensen, “Modeling the indoor mimo wireless channel,” IEEE Transactions on Antennas and Propagation, vol. 50, no. 5, pp. 591–599, 2002.
 [32] Michael A Jensen and Jon W Wallace, “A review of antennas and propagation for mimo wireless communications,” IEEE Transactions on Antennas and Propagation, vol. 52, no. 11, pp. 2810–2824, 2004.
 [33] Hadi Noureddine and Qianrui Li, “A twostep compressed sensing based channel estimation solution for millimeter wave mimo systems,” in Colloque GRETSI, 2017.
 [34] Mathew K Samimi and Theodore S Rappaport, “3d millimeterwave statistical channel model for 5g wireless system design,” IEEE Transactions on Microwave Theory and Techniques, vol. 64, no. 7, pp. 2207–2225, 2016.