Downlink channel spatial covariance estimation in realistic FDD massive MIMO systems

DO WNLINK CHANNEL SP A TIAL CO V ARIANCE E STIMA TION IN REALISTIC FDD MASSIVE MIMO SYSTEMS Lor enzo Mir etti ⋆ Renato L.G. Cavalcante † Slawomir Sta ´ nczak † ⋆ EURECOM † Fraunhofer Heinrich Hertz Instit ute and T echnical University of Be rlin ABSTRA CT The kno wledge of the downlink (DL) chan nel spatial cova riance ma- trix at the BS is of fundamental importance for large-scale array systems operating in frequency division duplex i ng (FDD) mode. In particular , this kno wledge plays a key role in the DL channel state information (CS I) acquisition. In the massiv e MIMO regime, tradi- tional schemes based on DL pilots are sev erely limited by the co- v ari ance feedback and the DL training overhead . T o ov ercome this problem, many authors have proposed to obtain an estimate of the DL spatial cov ariance based on uplink (UL) measurements. Ho w- e ver , many of these approaches rely on simple channel models, and they are difficult to extend to more complex models that take into account impo rtant ef fects of propaga tion in 3 D en vironments and of dual-polarized antenna arrays. In this study we propose a no vel technique that takes into account the aforemention ed effects, in com- pliance with the requirements of modern 4G and 5G system design s. Numerical simulations sho w t he ef fectiv eness of our approach. Index T erms — Massi ve MIMO, FD D, cov ari ance matrix, 3 D propagation, dual-p olarized arrays 1. INT RO DUCTION In this study we propose a technique to estimate the do w nli nk (DL) channel spatial co variance matrix R d in realistic massi ve MIMO systems operating in frequency div ision duple xing (FDD) mode. Al- though t ypical massiv e MIMO systems operate in time div i sion du- plexing (TDD) mode, the e xtension of this technology to FDD mode is of great practical interest [1] [2, Chapter 8.4]. The capability of the base station (BS) to access accurate and efficient R d estimates has emerged as an enabling technology to address practical implemen- tation issues of FDD l arge-scale array systems [3, 4], as it prov ides long-term information that is essential for beamforming and for CSI acquisition [5–9]. Con ventional FDD systems typically acquire R d by using tr a- ditional DL training and uplink ( U L ) cov ariance feedback schemes. Ho we ver , in massi ve MIMO systems, due to the large si ze of the cov ariance matrices and to the large DL training ov erhead, the tra- ditional schemes become unfeasible. T o ov ercome the drawbacks of con ventional systems, in [10], we propose a scheme to infer R d from t he observed UL cov ariance R u . This approach is based on projection methods and has many benefits. In particular, i t elimi- nates the continuou s DL training and cov ariance f eedback loop re- quired by con ventiona l direct R d estimation t echniques. Moreove r , it is com pletely transparent to the user equip ment (UE), hence it can be implemented in compliance with current standards. Related stat e-of-t he-art techniques for UL to DL cov ariance ma- trix con version i n FDD systems have been considered in [11–14]. The main li mitation of [10] and of the related works [11 –13] is that they are based on simple channel models that do not meet the requirements of modern 4G and 5G system designs. More pre- cisely , the approaches in [11–13] seem hard to generalize to ar- rays with arbitrary geometries and non-isotropic antennas. Further- more, [10–13] do not consider propagation effects of 3D en viron- ments and, most i mportantly , dual-polarized antenna arrays. In con- trast, the techniqu e proposed in [14] is able to cope with these design requirements. Ho wever , it is a machine learning approach that relies on the acquisition of a training set, and hence it is significantly more complex . In t his study , we propose a simple t r ai ning-free approach that takes into accoun t the aforementioned ef fects. T o this end, i n Sect. 2 we present a realistic multipath channel model, and we deriv e ex- pressions f or R d and R u under the assumption of bo th narro w-band and wide-band OFDM systems. Then , in Sect. 3, we describe the proposed scheme which infers R d from R u by exp l oiting the pro- posed cova riance model. The resulti ng algorithm is based on the joint estimation of the two angular po wer spectra for the v ertical (V - APS) and for the horizontal (H-AP S) polarization (defined in Sect. 2). T he ke y idea behind our approach is the definition of a suitable Hilbert space that allows us to f ormalize the joint V -APS and H-APS estimation problem as a con vex feasibility problem. This enables us to adopt standard projection-based solutions inspired by [10]. F ur- thermore, i n S ect. 4 we provid e implementation details for t he case of a cross-polarized uniform planar array (UP A) at the BS. Finally , in S ect. 5 we ev aluate the prop osed approach by means of numerical simulations. Notation: W e use boldface to denote vectors and matrices. ( · ) T and ( · ) H denote respecti vely the transpose an d Hermitian transpo se. By defining the set I ⊂ R 2 , C 0 [ I ] and L 2 [ I ] denote, respectiv ely , the set of all continuous functions and the set of all squa re Lebesgue integrable functions ov er I . ℜ [ · ] and ℑ [ · ] denote respecti vely the real and the imaginary parts. W e denote the imaginary unit by j . Throughout the paper , superscripts ( · ) u and ( · ) d indicate respec- tiv ely UL and DL matrices, vectors, or functions when we need to emphasize the dependen cy on the carrier frequency . 2. S YSTEM MODEL W e consider a MU-MIMO channel between a BS w i th N ≫ 1 an- tennas and single-antenna UE s. W e denote by h the channel vector of an arbitrary user . In the remainder of this section we describe the underlying models that we use for designing our scheme. Fi rst, in S ect. 2.1 we re vie w a widely-conside red narro w-band multipath model that takes into account 3D propagation and polarization ef- fects. Then, in Sect. 2.2 we present analytical expressions for R d and R u based on the considered channel model. Finally , in Sect. 2.3 we obtain analogou s expression s also for wide-band OFDM sys- tems. Due to the space limitation, all proo fs of this section are omit- ted. 2.1. Realistic Directional Multi-path Model In this section we consider an extension of classical 2D directional multi-path channel models ( e. g., the ones adopted in [6] and [15]) that take into acco unt 3D propagation and polarization effects. With this ex tension, we sho w later in Sect. 3 that we are able to address the UL-DL covarian ce con version problem by using projection meth- ods in a Hilbert space different from that in [10]. In more detail, by dropping the frequenc y dependent superscript f or simplicity , we model the channel vector h at an arbitr ary time t = t 0 according to the 3GPP narrow-band clustered directional multi-path model [16, Eq. (7.3-22)]. In this model, we hav e h = P N c c =1 h c , where h c := r α c N p N p X i =1 A ( θ ic )     e j ϕ V V ,ic e j ϕ V H,ic √ K ic e j ϕ H V ,ic √ K ic e j ϕ H H,ic     B ( φ ic ) H . The notation used here is defined as follows: • N c ∈ N denotes the number of clusters of scatterers, and N p ∈ N denotes the associated number of subpaths. T his terminology de- riv es from the classical geometry-based stochastic channel model (GSCM) [17, Chapter 7]. • θ ic ∈ R and φ ic ∈ R are, respective ly , either the direction of de- parture (DoD) and of arri val (DoA) of subpath i of cluster c for the DL case, or the DoA and DoD of sub path i of cluster c for the UL case. T he directions θ ic and φ ic are defined as tuples taking val- ues in the set Ω := [ − π , π ] × [0 , π ] , which represents the azimuth and the zenith of a spherical coordinate system. They are drawn independe ntly fr om a continuou s j oint distri bution f c ( θ , φ ) , and they are assumed to be equal for UL and DL. This DoD/DoA statistical modeling approach, which is very popular in the litera- ture [6, 15, 17], generalizes the model gi ven by 3GPP [16], where only the main cluster angles are random and the subpaths angles are obtained from tables. • A : Ω → C N × 2 is the dual polarized antenna array response of the BS. In FDD systems, A d is different f rom A u . The columns of A are denoted by [ a V , a H ] := A , and they rep- resent the array responses for, r espectively , the vertical and the horizontal polarization. Giv en an element a ij of A , we assume ℜ{ a ij } , ℑ{ a ij } ∈ C 0 [Ω] . • α c > 0 is the av erage power of all t he subpaths of cluster c , and it is assumed to be equal for UL and DL, which i s a reasonable assumption for current FDD systems [18]. • B : Ω → R 1 × 2 is the fr equency independent dual polarized an- tenna radiation pattern of the UE . The columns of B are denoted by [ b V , b H ] := B , and they represent, respecti vely , t he radiation patterns for the vertical and for the horizontal polarization. W e assume ℜ{ b V } , ℑ{ b H } ∈ C 0 [Ω] . • The random matrix M ic :=    e j ϕ V V ,ic 1 √ K i e j ϕ V H,ic 1 √ K i e j ϕ H V ,ic e j ϕ H H,ic    , models the fading of the vertical and horizontal polarization, and also of the cross-polarization terms caused by the polar- ization changes t hat the electromagnetic wav es undergo dur- ing the propagation. The components of the tuple ϕ ic := { ϕ V V ,ic , ϕ V H ,ic , ϕ H V ,ic , ϕ H H ,ic } are i.i.d. random variables, uniformly distributed in [ − π , π ] . The UL and DL phases are assumed independent. The parameters K ic ∈ R , usually termed as cross polarization po wer ratios ( X P Rs), are assumed to be i.i.d. r andom v ariables and to be equal for UL and DL. This polarization propagation model is i dentical to the one suggested by [17, Chapter 7], where the two polarizations are assumed to experien ce independent fading. W e point out that, in contrast to [16, Eq. (7.3-22)], this model does not take into account the time dependent phase t erm e j 2 πν ic t , where t is the time and ν ic is the Doppler shift of subpath i of clus- ter c , which models deterministically the short-term time ev olution of the channel. Howe ver , as the focus of this work is on the long- term channel statist i cs, we consider only a l ong-term time ev olution model, give n in a statistical sense. More precisely , we model the time ev olution of the channel as follo ws. T he fast time-varyin g pa- rameters θ ic , φ ic , K ic and ϕ ic are drawn independently and kept fixed at intervals corresponding to the coherence time T c (“block- fading” assumption). The slow ti me-varying parameters α c and f c are assumed constant ov er a windo w T W S S , with T W S S ≫ T c . This model reflects the classical “windowed WSS” assumption, which ap- proximates the channel as wide-sense stationary (W SS) for a giv en time windo w T W S S , which is usually several order of magnitude larger t han T c [5, 6]. 2.2. E xpression f or the Spatial Cov ariance Matrix In the next proposition we present an expression for the spatial co- v ari ance matrices R d and R u on which the DL covarian ce estima- tion scheme proposed in Sect. 3 i s based. Proposition 1. By assuming t he m odel intro duced in Sect. 2.1, the spatial covariance matrices R d := E  h d ( h d ) H  and R u := E  h u ( h u ) H  take the following forms: R d = Z Ω ρ V ( θ ) a d V ( θ ) a d V ( θ ) H d θ + Z Ω ρ H ( θ ) a d H ( θ ) a d H ( θ ) H d θ , (1) R u = Z Ω ρ V ( θ ) a u V ( θ ) a u V ( θ ) H d θ + Z Ω ρ H ( θ ) a u H ( θ ) a u H ( θ ) H d θ , (2) wher e the f unctions ρ V , ρ H : Ω → R + , referr ed to, r espectively , as “vertical polarization angular power spectrum” (V -A PS) and “hor- izontal polarization angular power spectrum” (H-APS), are define d to be ρ V ( θ ) := N c X c =1 α c Z Ω f c ( θ , φ )  b 2 V ( φ ) + 1 K b 2 H ( φ )  d φ , ρ H ( θ ) := N c X c =1 α c Z Ω f c ( θ , φ )  b 2 H ( φ ) + 1 K b 2 V ( φ )  d φ . Her e 1 /K := E [1 /K ic ] is the avera ge effect of the XPRs K ic . By recalling the notation defined in Sect. 2.1, we highlight that the V -APS and the H-AP S do not depend on the carrier frequenc y . Furthermore, we hav e that ρ V , ρ H ∈ L 2 [Ω] . 2.3. OF DM Systems W e no w show that expressions (1) and (2) (and hence the algorithms in Sect. 3) carry ove r to wide-band OFDM systems by extending the model in S ect. 2 with t he approach in [16] and [17, Chapter 6] for the “tapped delay line” model. More precisely , we consider a wide-band channel in an under-spread en vironment; i.e., with delay spread T s ≪ T c . By denoting wit h l ∈ N the discrete time index of the l th tap of the sampled impulse response, the channel vector ˜ h [ k ] in the sub-carrier domain is giv en by [19, Chapter 3.4]: ˜ h [ k ] = L − 1 X l =0 h [ l ] e − j 2 π k l N s , h [ l ] = N c X c =1 h c δ [ l − l c ] , (3) where { h c } c =1 ,...,N c are defined in Sect. 2.1, l c ∈ N denotes the discrete time de lay of all the subpaths belonging to cluster c , L is the impulse response length, N s is the chosen OFDM block length, and k = 0 , . . . , ( N s − 1) is the sub-carrier index. With t his model in hand, we can deriv e expression s f or the spatial cov ariance matrices in the sub-carrier domain. T hey are equi valent to the ones giv en by (1) and (2), and they do not depend on the sub-carrier index. More precisely , we hav e: Proposition 2. By assuming the wide-band OFDM channel model in (3 ) , the spatial covariance matrices R d k := E h ˜ h d [ k ]( ˜ h d [ k ]) H i and R u k := E h ˜ h u [ k ]( ˜ h u [ k ]) H i for a given sub-carrier k satisfy R d k = R d , R u k = R u , wher e R d and R u ar e given by (1) and (2) , and they do not depend on the sub-carrier index. 3. CHANN EL SP A TIAL CO V ARIANCE CONVERSION W e no w propose a practical FDD DL cov ariance esti mation scheme based on the channel model described in Sect. 2. The estimates of the DL channel cov ariance matrix R d are obtained from the UL channel cov ariance matrix R u by performing the following two-step scheme: 1. Giv en R u , we obtain an estimate ( ˆ ρ V , ˆ ρ H ) of ( ρ V , ρ H ) from (2) and kno wn properties of ( ρ V , ρ H ) . 2. W e compute the estimated R d by using (1) with ( ρ V , ρ H ) replaced by their estimates ( ˆ ρ V , ˆ ρ H ) . In this section, we assume perfect kno wledge of A u , A d , and R u , while later i n Sect. 4 and 5 we assume that the BS have access only to noisy estimates of R u . The core idea of the proposed scheme is that it is possible to address the joint V - APS and H-APS estimation problem of the first step as a con vex feasibility pr oblem , which enables us to apply so- lutions based on projection methods. W e point out that the related approaches in [10] cannot address properly the problem considered in t his paper because they are based on a Hilbert space that i s not ap- propriate to represent the esti mandum ( ρ V , ρ H ) resulting from the channel model we consider here. T o derive the proposed approaches, we first re write (2) as a sys- tem of equations of the form r u m = Z Ω ρ V ( θ ) g u V ,m ( θ ) d 2 θ + Z Ω ρ H ( θ ) g u H,m ( θ ) d 2 θ , (4) where r u m ∈ R is the m th element of r u := vec (  ℜ{ R u } ℑ{ R u }  ) , g u ( · ) ,m : Ω → R is the corresponding m t h coordinate function of vec (  ℜ{ a u ( · ) ( θ )a u ( · ) ( θ ) H } ℑ{ a u ( · ) ( θ )a u ( · ) ( θ ) H }  ) , and m = 1 , . . . , M , with M = 2 N 2 . No w l et H := L 2 [Ω] × L 2 [Ω] be the Hilbert space of tuples of biv ariate square-integrable real func- tions equipped with the following i nner product h ( f V , f H ) , ( g V , g H ) i := Z Ω f V ( θ ) g V ( θ ) d θ + Z Ω f H ( θ ) g H ( θ ) d θ . (5) Based on the model in Sect. 2, ( ρ V , ρ H ) and ( g u V ,m , g u H,m ) are mem- bers of H , thus (4) can can be rewritten as r u m = h ( ρ V , ρ H ) , ( g u V ,m , g u H,m ) i m = 1 , . . . , M . By using the set-theoretic paradigm [20–23], we obtain an estimate ( ˆ ρ V , ˆ ρ H ) of ( ρ V , ρ H ) by solving on e of the two f oll o wi ng feasibility problems: find ( ˆ ρ V , ˆ ρ H ) ∈ V := ∩ M m =1 V m 6 = ∅ , (6) find ( ˆ ρ V , ˆ ρ H ) ∈ C := V ∩ Z 6 = ∅ , (7) where V m := { ( h V , h H ) ∈ H : h ( h V , h H ) , ( g u V ,m , g u H,m ) i = r u m } are hyperplanes and Z := { ( h V , h H ) ∈ H : ( ∀ θ ∈ Ω) h V ( θ ) ≥ 0 , h H ( θ ) ≥ 0 } is the cone of tuples of non-nega tiv e functions. W e solve problem (6 ) by computing the projection onto the linear variety V , while problem (7), which takes i nto account also the positivity of ρ V and ρ H , is solved via an iterative projection method called extr apolated alternating pr ojection method (EAP M) . For the details about the solutions of the considered feasibility problems, we refer to [10] and to the references herein. The choice of solving either (6) or (7) leads to two variants of the proposed scheme with different comple xit y and accuracy , and t hey are referred here as Algorithm 1 and Algorithm 2 . More precisely , Algorithm 1 can be implemented as a simple matrix multiplication of the form r d = F r u , where F depends just on the array geometry and can be computed once for the entire system lifetime. In con- trast, A lgorithm 2 requires iterativ ely the ev aluation of i ntegrals of the form R Ω x ( θ ) d 2 θ (see [10] for details). 4. IM PLEMENT A TION F OR UNIFORM PLANAR ARRA Y WITH P AIRS OF CR OSS-POLARIZED ANTENNAS In this section we describe implementation aspects for a cross- polarized uniform planar array (UP A) , defined here as a rectangular grid of i dentical and equispaced antenna elements, each of them composed of a pair of two vertically polarized antennas with a po- larization slant of ± 45 ◦ . W e denote by N V and N H the number of vertical and horizontal elements, respectiv el y , and by d the inter- antenna spacing. W e further denote by x ( u, v , 1) the antenna in position ( u, v ) , u = 1 , . . . , N V and v = 1 , . . . , N H , with +45 ◦ po- larization slant, and by u ( u, v , 2) the co-located antenna with − 45 ◦ polarization slant. For this antenna array , the cov ariance matrix has the followin g structure: Proposition 3 (St ructure of the U P A Cov ari ance Matrix) . By let- ting h :=  h T 1 h T 2  T , wher e the channe l coefficien t for an- tenna x ( u, v , k ) corr esponds to the n th element of the vector h k ∈ C N V N H × 1 , with n = ( u − 1) N H + v , and by assuming without l oss of generality that N V ≥ N H , the covariance matrix takes on the following block struc tur e: R =  B 1 B H 2 B 2 B 3  ∈ C 2 N V N H × 2 N V N H , wher e every m acro-b lock B l ∈ C N V N H × N V N H , l = 1 , 2 , 3 , is Hermitian and it has the following block structur e: B l =        B l, 1 B l, 2 B l, 1 B l, 3 B l, 2 B l, 1 . . . . . . . . . . . . B l,N V . . . B l, 3 B l, 2 B l, 1        , wher e every block B l,i ∈ C N H × N H , i = 1 , . . . N V has identical diago nal entries b li , and every bloc k B l, 1 is Hermitian T oeplitz. The proof is omitted here but we point out that it follows by direct inspection of the matrices a V ( θ )a V ( θ ) H and a H ( θ )a H ( θ ) H of (1 ) and (2), where the elements of the array responses a V and a H are arranged with the same scheme adopted for h . T able 1 . General simulation parameters Carrier frequenc y ( f c ) 1.8 GHz for UL, 1.9 GHz for DL System type Narro w-band or wide-band OFDM BS 8x4 cross-polarized UP A d = λ u / 2 UE Single antenna, vertically polarized Antennas radiation pattern 3GPP [16, Section 7.1], 3D-UMa The structure of the UP A cov ariance matrix R described in Prop. 3 has the following consequences in practical implementations of the algorithms presented in Section 3: • R can be bijectiv ely vectorized by using only M = 6( N H + ( N V − 1)( N 2 H − N H + 1)) real numbers, comp ared to the M = 2( N V N H ) 2 elements giv en by the vectorization operation de- fined in Sect. 3. • Any estimate ˆ R of t he cov ariance matrix R u (for example, ob- tained from the sample cova riance matrix similarly to [10, Sect. 4.2.]) can be further improved by substituting each element ( i, j ) with the arithmetic av erage of all the elements that are assumed to be identical. 5. S IMULA TION In this secti on we ev aluate the proposed algorithms by simulating a communication scenario wi th system parameters giv en in T able 1. T he channel coefficients are gi ven by the narrow-band multipath model described in Sect. 2.1, with parameters randomly drawn as follo ws: • Cluster po wers α c are drawn uniformly from [0 , 1] and further normalized such that P N c c =1 α c = 1 . • The X P Rs values K ic are drawn f rom a log-Normal distr i bution with parameters ( µ XPR , σ XPR ) = (7 , 3)[ dB ] . This is identical to the 3GPP model [16, Sect. 7.3, Step 9], wi t h parameters for 3D- UMa, NLOS propaga tion. • The angles θ ic , φ ic are generated from t he jointly Gaussian distribution f c ( θ , φ ) = f BS ,c ( θ ) f UE ,c ( φ ) , where f BS ,c ∼ N ( µ BS , σ 2 BS I ) and f UE ,c ∼ N ( µ UE , σ 2 UE I ) , and where the clus- ters means and angular spreads µ BS := [ µ BS ,a µ BS ,z ] , σ 2 BS := [ σ 2 BS ,a σ 2 BS ,z ] , µ UE := [ µ UE ,a µ UE ,z ] , σ 2 UE := [ σ 2 UE ,a σ 2 UE ,z ] , are dra wn as follows: µ BTS ,a , µ UE ,a are uniformly drawn from  − 2 π 3 , 2 π 3  , µ BTS ,z , µ UE ,z from  π 4 , 3 π 4  , σ BTS ,a from [3 ◦ , 5 ◦ ] , σ UE ,a from [5 ◦ , 10 ◦ ] , σ BTS ,z from [1 ◦ , 3 ◦ ] , and σ UE ,z from [3 ◦ , 5 ◦ ] . This choice of parameters is inspired by experimen tal properties of ρ V and ρ H gi ven by [16], e. g. the elev ation angular spread is usually narrower than the azimuth one. • T o simulate different UE antenna orientation, the UE antenna ar- ray response is giv en by applying a 3D rotation t o the antenna radiation pattern as described in 3GPP [16, Sect. 5.1.3], w i th pa- rameters α, β , γ ∼ U  0 , π 6  . The BS is assumed to ha ve access t o t he estimated UL cova ri- ance matri x ˆ R u obtained from N s = 1000 noisy channel esti- mates ˆ h u = h u + z , z ∼ C N ( 0 , σ 2 z I ) i.i.d., w i th noise po wer defined by setting an averag e per-antenna SNR est to S NR est := tr { R u } / ( N σ 2 z ) = 10 [dB], where N = 2 N V N H denotes the num- ber of BS antennas. The estimate ˆ R u is computed by projecting the sample cov ariance matrix onto the space of positiv e semi-definite 10 -4 10 -3 10 -2 10 -1 0 0.2 0.4 0.6 0.8 1 (a) Normal ized Eucli dean distance 10 -4 10 -2 10 0 10 2 0 0.2 0.4 0.6 0.8 1 (b) Principal subspaces distance Fig. 1 . Empirical CDF of the squared error (SE ) matrix C as described in [10 , Sect. 4.2], and by applying the cor- rection procedure described i n Sect. 4. Furthermore, the proposed algorithms are implemented by explo iting the efficient v ectorization for UP A described in Sect. 4. The accuracy of an estimate ˆ R d of R d is ev aluated in terms of the square error SE := e 2 ( R d , ˆ R d ) , where e ( · , · ) is a giv en error metric. In particular, we consider as error metrics the normalized Frobenius norm and the 90% Grassmanian principal subspace dis- tance defined in [10, Sect. 5]. T o ev aluate the proposed algorithms we use as a baseline an es- timate of the DL cov ariance matri x obtained from ˆ h d samples with the same technique f or the estimati on of ˆ R u . Furthermore, we also compare t he proposed approaches with a solution that relies on a pre-stored dictionary of cov ariance matrices ( ˆ R u , ˆ R d ) and based on the Wiener filter, similar to the approach proposed in [14] that was already analyzed with the preliminary results in [10, Sect. 5]. The results are shown in Figure 1, which sho ws the empirical cu mulativ e distribution function (CDF ) of the SE for t he two chosen metrics, obtained by drawing independent realizations of the quantities that are assumed to stay fixed for T WSS (i.e. by dra wing a ne w V -AP S and H-APS). The simulation confirms that the proposed algorithms are able t o provide an accurate DL estimate by using only UL training, thus it can be used as an ef f ectiv e solution to the DL channel cova ri- ance acquisition problem. W it h respect to Algorithm 1 , Al gorithm 2 is slightly more accurate, as it considers also the non-negati vity property of the V -APS and H-APS, but it pays a price in terms of increased complexity . W e point out that, as opposed to t he Wiener filter approach and to all other techniques based on supervised ma- chine learning t ools, the two proposed algorithms do no t require an y training phase. 6. RE FERENCES [1] E. Bj ¨ ornson, E.G. Larsson, and T .L. Marzetta, “Massi ve MIMO: T en myths and one critical question, ” IEEE Commu- nications Ma gazine , vol. 54, no. 2, pp. 114–123 , 2016. [2] T .L. Marzetta, E.G. Larsson, H. Y ang, and H.Q. Ngo, Funda- mentals of Massive MIMO , Cambridge Unive rsity Press, 2016. [3] F . R usek, D. Persson, B. K. Lau, E. G. Larsson, T . L. Marzetta, O. Edfors, and F . Tufv esson, “S caling up MIMO: Opportuni- ties and challenges wi th very large arrays, ” IE EE Signal Pro - cessing Ma gazine , vol. 30, no. 1, pp. 40–60, 2013. [4] H. Huh, A. M. Tu l ino, and G. Cai r e, “Network MIMO with linear zero-forcing beamforming: Large system analysis, im- pact of channel estimation, and reduced-comple xity schedul- ing, ” IEEE T ransactions on Information Theory , vol. 58, no. 5, pp. 2911–2 934, 2012. [5] A. Adhikary , J. Nam, J.-Y . Ahn, and G. C ai r e, “Joint spatial di vision and multiplexing –The large-scale array regime, ” IEE E transac tions on information theory , vol. 59, no. 10, pp. 6441– 6463, 2013 . [6] H. Y in, D. Gesbert, and L. Cottatellucci, “Dealing with inter- ference in distributed large-scale MIMO systems: A statistical approach, ” IE EE J ournal of Selected T opics in Signal Pr ocess- ing , v ol. 8, no. 5, pp. 942–953, 2014. [7] M. R. Akdeniz, Y . Liu, M. K. Samimi, S. Sun, S. Rang an, T . S. Rappaport, and E. Erkip, “Millimeter wav e channel modeling and cellular capacity ev aluation, ” IEEE journal on selected ar eas in communications , vol. 32 , no. 6, pp. 1164–1179, 2014. [8] J. Fang, X. Li, H. Li, and F . Gao, “Lo w-rank cov ariance- assisted do wnlink training and channel estimation for F D D massi ve MIMO systems, ” IEE E T ransactions on W ireless Communications , vo l . 16, no. 3, pp. 1935–1947, 2017. [9] Jayesh H Kotecha and Akbar M S ayeed, “T ransmit signal design for optimal estimation of correlated MIMO channels, ” IEEE T ransa ctions on Signal P r ocessing , vol. 52, no. 2, pp. 546–55 7, 2004. [10] L. Miretti, R. L. G. Cav alcante, and S . Stanczak, “FD D Mas- siv e MIMO channel spatial cov ariance con version using pro- jection methods, ” IEEE ICASSP , April 2018. [11] Y ing-Chang Li ang and Francois P . S. Chin, “Do wnlink channel cov ariance matr i x (DCCM) estimation and its applications in wireless DS-CDMA systems, ” IEE E Jou rnal on Selected Ar eas in Communications , vo l. 19, no. 2, pp. 222–232, 2001. [12] M. Jordan, A. Dimofte, X. Gong , and G. Ascheid, “Con version from uplink to downlink spatio-temporal correlation with cu- bic splines, ” in IEEE 69th V ehicular T echno logy Confer ence , 2009, pp. 1–5. [13] S. Haghighatshoar , M.B. Khalilsarai, and G. Caire, “Multi- band cov ariance interpolation with applications in massive MIMO, ” CoRR , vol. abs/1801 . 03714, 2018. [14] A. Decurninge, M. Guillaud, and D.T .M. Slock, “Channel co- v ari ance estimation in massi ve MIMO frequenc y di vision du- plex systems, ” in IEEE Globecom , 2015, pp. 1–6. [15] S. Haghighatshoar and G. Caire, “Massiv e MIMO channel subspace estimation from low-dimensional projections, ” IEE E T ransactions on Signal Pr ocessing , vol. 65, no. 2, pp. 303–3 18, 2017. [16] 3GPP , “Study on 3D channel model for L T E (release 12), ” T ech. Rep. TR 36.873 V12.6.0, 3rd Generation Partnership Project (3GPP), 2017. [17] A.F . Molish, W ire less Commun i cations , W iley , second edition, 2010. [18] S. Imtiaz, G. S . Dahman, F . Rusek, and F . T ufvesson, “On the directional reciprocity of uplink and do w nli nk channels in fre- quenc y division duplex systems, ” in 2014 IEEE 25th A nnual International Symposium on P ersonal, Indoor , and Mobile Ra- dio Communication (PIMRC) , Sept 2014, pp. 172–176. [19] D. Tse and P . V iswanath, Fundamentals of wireless commun i- cation , Cambridge univ ersity press, 2005. [20] P . L. Combettes, “The foundations of set theoretic estimation, ” Pr oceedings of the IEEE , vol. 81, no. 2, pp. 182–208 , 1993. [21] H. S tark, Y . Y ang, and Y . Y ang, V ector space pr ojections: a nu- merical appr oach to signal and image pr ocessing , neur al nets, and optics , John Wile y & Sons, Inc., 1998. [22] I. Y amada and N. Ogura, “ Adapti ve projected subgrad i ent method for asymptotic minimization of sequence of nonneg - ativ e con vex functions, ” Numerical functional analysis and optimization , vol. 25 , pp. 593–617 , 2005. [23] S. Theodoridis, K. Slav akis, an d I. Y amada, “ Adapti ve learning in a world of projections, ” IEEE Signal Pr ocessing Maga zine , vol. 28, no. 1, pp. 97–123 , 2011.