Random Access in Massive MIMO by Exploiting Timing Offsets and Excess Antennas

1 Random Access in Massi v e MIMO by Exploiting T i ming Of fsets and Ex cess Antennas Luca Sanguinetti, Senior Member , IEEE , Antonio A. D’Amico, Michele Morelli Member , IEEE , and Merouane Debbah, F ellow , IEEE Abstract Massi ve MIMO systems, wher e ba se stations are eq uipped with hun dreds of antennas, are an attractive way to handle th e rapid growth o f data traffic. As th e num ber of u ser equip ments (UEs) increases, the initial access and handover in contempor ary networks will be flooded by user co llisions. In this paper, a r andom a ccess protoco l is p r oposed th at resolves collision s and perform s timing estimation by simp ly utilizing th e large numbe r of antennas envisioned in Massi ve MIMO n e tworks. UEs enter ing the ne twork per form spreading in both time and freq uency domains, and their timing offsets are estimated at the base station in closed-f orm using a subspace decompo sition app roach. T h is in formation is u sed to compute channel estimates that are subsequently e m ployed by the base station to commu nicate with the detected UE s. T he fa vorable propag ation c o nditions of Massive MIMO suppr ess interfer ence amon g UEs whereas the inherent timing misalignments improve the de te c tio n capabilities of th e protocol. Numerical resu lts are used to validate the perfo r mance of the pr oposed proc e dure in cellular network s under u ncorre lated an d correlate d fading channels. W ith 2 . 5 × 10 3 UEs th at may simultaneously be come activ e with probability 1% and a total of 16 frequency-time codes (in a giv en rand om access block), it turns out that, with 1 00 antennas, th e proposed procedu re successfully detects a gi ven UE with probab ility 75% while providing reliable timing e stimates. I . I N T R O D U C T I O N Massive M IMO is consi d ered as one of th e mo s t promising soluti on to handl e th e dramati c increase of mobile data traffic in t he years to come [1], [2]. The basic premise behind Massiv e L. Sanguinetti, A. A. D’Amico, M. Morelli are with the Univ ersity of Pi sa, Dipartimento di Ingegn eria dell’Informazione, Italy (luca.sanguinetti@unipi.it). L S anguinetti is also with the Large Systems and Networks Group (LAN E AS), Centrale- Sup ´ elec, Univ ersit ´ e Paris-Saclay , 3 r ue Joliot-Curie, 91192 Gi f-sur-Yv ette, France. M. Debbah is with t he Large Systems and Networks Grou p (LANEAS), CentraleSup ´ elec, Universit ´ e Paris-Saclay , 3 rue Joliot-Curie, 91192 Gif-sur -Yvette, France (m.debbah@centralesu pelec.fr) and also with the Mathematical and Algorithmic Sciences Lab, Huawei T echnologies Co. L td., France (merouane.debbah @huawei.com). This research has been supported by the ERC S tarting Grant 305123 MORE. A preliminary version of this paper was presented at the IEEE G l obal Communication Conference, W ashin gton DC, USA, Dec. 2016. 2 MIMO is to reap all t h e benefits of con ventional M IM O, but on a much greater s cale: a fe w hundred antennas are used at the base stati on (BS) to simultaneously serve many tens o f user equipment terminals (UEs) in t h e same frequency-time resource using a time d ivision duplexing protocol. The b enefits of Mass iv e MIMO in t erms of area throughput, power consumption and ener gy effi ciency ha ve been extensively s tudied in recent years and are now adays well understood [2]–[7]. On the other hand, the potential benefits of Massive MIMO in the network access functionalities ha ve receiv ed little attention so far [8]–[10]. These network access functionalit i es refer t o all the functions that a UE needs to go through in order to establish a commu nication link with the BS for data transmis s ion and reception. Next, we first re vise the n et work entry procedure specified by the L TE standards, and then briefly describe how this procedure has been recently extended to M ass iv e M IMO systems. A. Random access in L TE The L TE standards specify a n et work entry p rocedure called random access (RA) by which uplink (UL) signals can arrive at the BS aligned in t ime and with approximately th e same power lev el [11 ]. In i t s basic form, t h e RA function is a contention-based procedure, which essent ially dev elo p s through the four s teps specified in Fig. 1(a). In Step 1, each UE tryin g to establish a com m unication link first acquires b as i c synchronizatio n from eNodeB (e.g. determining L T E parameters, frequenc y synchronization, and frame ti m ing), and then accesses th e network using the so-called RA block (or RA channel ), which is composed of a specified set of con s ecutiv e symbols and adjacent s ubcarriers. Each UE m akes us e of th e RA blo ck to transmit a pilot sequence random ly chosen from a predefined set. As a consequence of the di fferent UEs’ positions within the cell, RA signals are su b ject to UEs’ specific propagation delays and arrive at t he eNod eB at dif ferent time instants. In Step 2, the eNodeB detects each pil ot sequence and extracts the associated p hysical parameters (e.g., ti m ing advance and received power). T h en, it broadcasts a RA response message informing t he UEs associat ed to the detected sequences that the procedure has been successfully completed and giving instructions for subsequent data transmissio n . In Step 3, all the UEs that have selected one of the detected sequences, adjust their physical parameters and send a conn ection request. If multip l e UEs a ccess the network with the same pilo t sequence, then collision s occur at the eNodeB. The centralized contention resol u tion in Step 4 is a dem anding procedure m eant for identifying th e UEs that ha ve been detected in 3 !" #$%&' () *+, ' -)./01/%2/ 3.- /4 )5 6 #$ %&'() 7 22/88)#/84'%8/ 3.- /4 )9 6 :'% %/2-+' % )#/01/8- 3.- /4 ); 6 :/%-<$,+=/&):'%-/%-+'%) # / 8 ',1- + ' %)$%&)#/8'1 <2 /) 7, , ' 2 $ -+'%) 3.- /4 )> 6 (a) R A protocol i n L TE systems !" #$%%&'()#*#+),- .$/01 2) 3&4 1 5)-(67(/8( 9-5 (: ); < 3=(8 1 0(0).$/012)>88(%%) . ( % :1/%( 9-5 (: )? < @&%5 =&A75(0)B1/5(/5 &1 /) .(%1475&1/)$/0)3&415) .(:(5&5&1/ 9-5 (: )C < .(%1 7=8()>4418$5&1/ 9-5 (: )D < (b) T he SUCR RA protocol for Massive MIMO Fig. 1 : RA pro tocols for L TE systems and Ma ssi ve MIMO. Step 2 and for allocating to them resources for data transmission. The undetected UEs repeat the RA procedure after a random wa iting time. The RA procedure d escrib ed above is used in L TE for various functionalities: ini t ial access, handover , maintaining UL synchronization, and scheduling request. For each of them, there exists a variety of d i ff erent solu t ions in the literature. In the context of in i tial access and handover , an example of possible soluti ons is giv en by [11]–[18]. The methods illustrated in [12] perform code detection and ti ming re covery by correlating t he received samp l es with time-shifted versions of a training sequence. A simp le energy detector is employed i n [13] whereas a timing recovery scheme specifically devised for th e L TE UL is discuss ed in [11], and further enhanced in [14]. Schemes for ini tial access, based on subspace m et h ods, are proposed in [15] and [16]. A solution based on the generalized like lihood ratio test is proposed i n [17], whereas [18] illustrates a RA algorithm that exploits a u nique ranging symbol with a re petitive structure in the time-domain . B. Random access in Massive MIMO All the aforementio ned solutions can be applied to cellular networks wherein the number of UEs, that may potentially enter the network, is relatively small compared to the number of a vailable p i lot sequences. On the other hand, they will be flooded by collision s for a mu ch larger number of UEs , e.g., in the order of h undreds or thousands as en visioned in future networks. In 4 the conte xt of Massiv e MIMO, rece nt attempts in the abo ve direction can be found in [8]–[10]. The papers [8], [9] consider a crowded network in which UEs intermittently enter the network, whenev er they want to, by selecting a pilot sequence from a common pool. In particular , a coded RA proto col is presented in [8] le veraging the channel hardening p roperties of Massiv e MIMO, which allow t o view a set of contaminated RA signals as a graph code on which iterati ve belief propagation can be performed. The proposed solution outp erforms the con ventional ALOHA method at the price of an increased error rate, due to accumul ation of estimation errors in the belief propagation algo ri t hm. In [9], sum UL rate expressions are derived that take intra-cell pilot coll isions, intermit tent UE activity , and interference into account . These expressions are used to op t imize the UE activation p rob ability and pilot length. In [10], the authors consider a con ventional cellular network in which a UE needs to be assigned to a dedicated pilot s equence before transmitting data. In this context, the channel hardening and spatial re solutio n properties of M assiv e MIMO are used to derive a new protocol, called strongest -us er collisio n resoluti on (SUCR), which enables d istributed collis i on detection and resoluti on at the UEs. The four s teps of t he SUCR protocol are ill ustrated in Fig. 1(b). Unlike the L TE RA proto col in Fig. 1(a), in Step 2 the BS sends prec oded signals to all the RA pil o ts that are detected b y the BS in Step 1. Only the UE with the strongest recei ved s ignal (among thos e using the same code) retransmits in Step 3. If correctly detected by the BS, the UE will be admitted t o the payload transmission phase in Step 4. If not, it will repeat th e RA procedure after a random waiting tim e. T h e SUCR protocol can be u s ed as an add-on to con ventional L TE RA mechanisms . T wo extensions of the SUCR protocol are presented in [19], [20]. Both solutions ai m at im proving the d etection probability of the weaker UEs. This is achiev ed b y allowing these UEs to randomly select pilots from those that are not selected by an y UE in the initial step. All the above works consider perfectly frequency- and time-synchronized networks such that the orth o gonality of RA pi lot sequences is preserved at th e BS. Frequency errors during RA are mainly due to Doppler shifts and /or estimation errors occurring in the initi al downlink synchronization process. As such, they are normally small and result o nly in negligible phase rotations over one symb o l [15], [21]–[23]. Howe ver , phase rotations become significant over a RA blo ck spanning sev eral consecutive s ymbols. Timing errors are due to the different position s of UEs within the cell. In these circumstances, the receiv ed RA pilots are af fected by differe nt linear phase shifts in frequency-domain [21]. Therefore, in the presence of frequency and timin g errors the receiv ed RA pil ots transmitt ed over adj acent sub carriers and consecutive symb ols 5 are no longer orthogonal at the BS side. As a consequence, the performance of the proposed solutions may be substantially deteriorated. C. Contr i butions and outl i ne In thi s work, we propose a novel RA protocol which operates through the following three steps. In Step 1, each UE that wants to enter the network randomly selects a pair of predefined RA codes and perform spreading over t he RA block in both the frequency and time domains. At the BS, the spatial degrees of freedom provided by Massive MIMO system s are used toget h er with the inherent different ti m e instants of reception of UEs ’ signals (before the data transmiss i on begins) to resolve collis i ons. In particular , the large number of antennas at the BS i s first u sed to compute a reasonable approxi m ation of a sample covariance matrix, whi ch i s then employed by the minimum descripti on length (MDL) algorithm [30] to determi n e the number of frequency- domain codes for each give n tim e-domain code. This information is used for timing recovery through the estimation of signal parameters via rotati o n al in variance technique (ESPRIT) [24] that allows to comput e estimate of the tim ing offsets in clos ed-form . These estimates are exploit ed to comput e the least-sq uare (LS) estim ate of the channels of al l detected codes. Step 2 of the proposed procedure operates according to Step 2 of the SUCR prot ocol proposed in [10] and illustrated in Fig. 1(b); t hat is, the BS respo n ds by sending DL pilots th at are precoded using the channel estim ates. This allows to detect UEs us i ng the same RA codes in a distrib uted way; that is, only t h e UE with the strongest signal sho uld repeat the RA codes. Compared to [10], howe ver , a colli s ion occurs when two UEs select the same pair of cod es (in time and f requency domains) and are characterized by (nearly) the s ame ti ming offset. If this latter case does not apply , no collision occurs among the t wo UEs and each one i s allowed to retransmit t h e selected pair of codes (followed b y an UL mess age containing the uni que i dentity number of the UE). The two UEs will be discrim inated in Step 3 by u s ing the LS channel estim ates obtained in Step 1. This improves the detection capabilities of the proposed protocol. All this is achiev ed at the price of an increased com putation complexity compared t o [10]. Numerical results show that the proposed RA procedure largely outperforms a “baseline” approach in which colli sion-av oid ance entirely rests on the choice of different code sequences, while providing reliable timing estimates. Compared t o its prelimin ary version presented in [25], this work is substantiall y different because of the foll owing reasons: (i) i t contains more technical details and applies to a multicell 6 network; (ii) it is developed and e valuated ov er a general correlated Rayleigh fading channel model; (iii) the full procedure is described until success n o tification is broadcasted by the BS . The remainder of this paper is organized as follows. Next section introduces basic notation and describes the Massiv e MIMO network with the underlying assu m ptions. Section III develops the proposed RA protocol by exploiting the large number of antennas at the BS and by assum ing that the resolution of the ESPRIT algorithm is sufficiently high such that all the t i ming of fset s of t he recei ved RA signals are accurately estimated. In Section IV, we consider a sim ple case study , in which two UEs choose the same time- and frequency-domain codes and are characterized by the same timing off set, and show what are t h e practical consequences of the finite reso l ution of the ESPRIT algorithm. Numerical results are giv en in Se ction V to validate the perf ormance of the proposed RA procedure in a Massiv e MIMO network with a finite number of BS ant enn as un der uncorrelated and correlated fading c hannels. Finally , the major conclusi ons and impli cations are drawn i n Section VI. Notation: Matrices and vectors are denoted by boldface letters, with I N being the identity matrix of order N . The transpose, conjugate-transpo se, and conj ugate of a matrix X are denoted by X T , X H , and X ∗ , respectively . W e use k·k and | · | to indicate the Eucli dean norm and the cardinality of a set, respectively . A random vector x ∼ C N ( ¯ x , R ) is com plex Gaussian distributed wi th mean ¯ x and cov ariance m atrix R . W e use x ⊙ y and x ⊗ y to denote t he Hadamard and Krono cker products between vec tors x and y , respecti vely . W e use a n ≍ b n to denote a n − b n → n →∞ 0 almost surely (a.s.) for two random sequences a n , b n . I I . N E T W O R K M O D E L Consider a Massiv e MIMO network based on OFDM with L cells, each comprising a BS with M antennas. W e deno t e by N FFT the numb er of subcarriers with frequency spacing ∆ f and call U j the total s et of UE s (activ e and inactive) t h at are in cell j . The network operates according to a tim e-division duplexing (TDD) protocol and th e ti me-frequency resources are divided i n to payload and RA blocks. The payload bl ocks are used for data t ransmission and consist s of τ C samples. 1 W e assume that, at any given time, only a subs et A j ⊆ U j of UEs is activ e for data transmissio n , with | A j | < τ C . The RA blocks are reserved for the inactive UEs, i.e. those i n the 1 The number of samples per block depends on the coherence bandwidth and coherence time of al l U E s. Since it is hard to dynamically adapt the netwo rk to these values becaus e the same protocol should apply to all UEs, a practical solution i s t o design the coheren ce block for the worst-case propagation scenario that the network should support. 7 set I j = U j \ A j that may become active , and consi s t of τ < τ C samples. W e further assume that the RA blocks of diff erent cells are allocated ov er different time and frequency resources such that no inter -cell interference arises among inactive UEs while accessing the n et work. Nev erth el ess , the inactive UEs in each cell j will be af fected by the int er -cell interference generated by the activ e UEs in A j ′ , wi th j ′ 6 = j . W i t hout lo s s o f generality , in the sequel w e concentrate on a generic cell j and omit the cell inde x for simplicity . A. Random access bloc k W e ass u me that each UE in I may become acti ve in a g iven RA block wi t h probabili ty p A and that the τ samples of each bl ock consis t s of Q con s ecutiv e OFDM symbols and N adjacent su bcarriers such that τ = QN . 2 After do wnlink synchroni zati on, a gi ven UE k in I , that would l ike to access the network, selects randomly a pair of cod es from the ort hogonal sets C N = { f 0 , . . . , f N − 1 } and C Q = { t 0 , . . . , t Q − 1 } , with { f i ∈ C N : f H i f i = N ∀ i } and { t i ∈ C Q : t H i t i = Q ∀ i } . W e denote by l k ∈ { 0 , . . . , N − 1 } and i k ∈ { 0 , . . . , Q − 1 } the code indices selected by UE k , and assume th at f l k and t i k are used in t h e frequency and time domain, respectiv ely , ov er the RA block. W e further assume that f l k belongs to the F o u rier basis with [ f l k ] n = e  2 π N nl k n = 0 , 1 , . . . , N − 1 (1) while no particul ar structure is assum ed for t i k ∈ C Q . Th e frequenc y codes { f 0 , . . . , f N − 1 } are selected from the Fourier basis because they allow us to us e an effi cient frequency d omain algorithm for the estimation of the tim ing m isalignment s between BS and UEs. The estimati on algorithm is based on the ESP RIT method and is described in Section III.B.2. An access att empt from UE k consist s in transmit ting t h e code m atrix f l k t T i k with a certain power level ρ k > 0 where √ ρ k  f l k t T i k  n,q = √ ρ k t i k ( q ) e  2 π N nl k (2) is transmitted over s ubcarrier n d uring OFDM symbol q . The value of ρ k depends on the number of RA attempts already made by UE k . Indeed, we assume t h at UE k ent ers the n et work with a relatively low power lev el ρ k = ρ min . If not admitted imm ediately , it retransm its in the next a vailable RA block by e xpo nentially increasing ρ k . If the maximu m power lev el ρ max is reached 2 Notice that an L T E resource block, over which the channel is assumed to be constant over t ime and f requenc y , spans Q = 14 OFDM symbols and N = 12 subcarriers, for a t otal of τ = 168 samples. 8 10 3 10 4 10 5 0 0.2 0.4 0.6 0.8 1 UE density , µ [UE/km 2 ] Pr  collision | f j t T i  p A = 0 . 5% p A = 1% p A = 2% Fig. 2: Probability of collision fo r a gi ven pa ir ( f j , t i ) vs. UE de n sity µ with activation pro bability p A = 0 . 5% , 1% and 2% for Q = 2 and N = 8 . A squ are cell of side length 500 m is conside r ed. and s till U E k has not succeeded, then i t starts the process again from the mi nimum power lev el ρ min . Clearly , ρ k = 0 if UE k does not want to ent er the network. W e call K ij = { k : i k = i, l k = j, ρ k > 0 } (3) the s et t h at contains the i ndices of all UEs that ut ilize code f j t T i , wi th K = ∪ i,j K ij being the index set of UEs transmi tting in the considered RA blo ck. Accordingly , the cardinality o f K ij is a binomial random variable dis tributed as |K ij | ∼ B  |I | , p A / ( QN )  where |I | is the number of inactive UEs in the considered cell and p A / ( QN ) is the probability th at e ach of t hem selects code f j t T i . Based on this model, a collision f or f j t T i occurs with probability [10] Pr ( collis ion | f j t T i ) = 1 −  1 − p A QN  |I | − |I | p A QN  1 − p A QN  |I |− 1 (4) and the a verage number of UEs selecting the same code f j t T i is E {|K ij |} = |I | p A / ( QN ) . T o provide realisti c values for these quant ities, l et u s consider a squ are cell of si de l ength 500 m wherein codes of length Q = 2 and N = 8 are used. Fig. 2 il lustrates the p robability of collision with p A = 0 . 5% , 1% and 2% for differe nt va lues of UE density 3 , µ [measured in UE/km 2 ]. W ith 3 Note that UE densities from µ = 10 2 UE/km 2 (in rural areas) to µ = 10 5 UE/km 2 (in shopping malls) ha ve been predicted in the METIS project [26]. 9 µ = 10 4 (which corresponds to | I | = 2500 ), the ave rage n u mber of i nactiv e UEs selecting t h e same code i s 0 . 78 , 1 . 56 and 3 . 12 for p A = 0 . 5% , 1% and 2% , respecti vely , leading to a collision with prob abi lity 0 . 18 , 0 . 4 7 and 0 . 82 . All coll isions mus t be detected and resolved before any UE can establish a data communication link with the BS . B. Channel model W e assume that the channel response can be approximated as constant and fl at-fading wi t hin a RA block and denot e by h k = [ h k 1 , . . . , h k M ] T ∈ C M the channel frequency response o f UE k at the BS antenna array over the considered RA block. W e assume correlated Rayleigh fading such that h k ∼ C N ( 0 , R k ) where R k ∈ C M × M is a positi ve semi-definite matrix with bounded spectral no rm [3]. The Gaussian distribution models the small -s cale f ading whereas R k is the spatial channel covariance matrix, which describes the macroscopic propagation effects (path loss and s hadowing), including the antenna gains and radiation patterns at the BS and UE. The normalized trace β k = 1 M tr ( R k ) (5) determines t he a verage channel gain from the BS to UE k . W e further assume t hat channel vectors { h k } satisfy the two foll owing conditi o ns [27]: 1 M h H k h k ≍ β k ∀ k (6) 1 M h H k h i ≍ 0 ∀ k , i, k 6 = i. (7) The first on e is known as channel hardening 4 and should be interpreted in the sense that the relativ e d eviation o f k h k k 2 from E {k h k k 2 } = tr( R k ) vanishes asym ptotically . The second condition is kno wn as fa vorable propagation 5 and makes the channels of two UEs orthogonal when the number of antennas gro ws unboundedly . This property makes interference between UEs vanish asymptotically . Note that channel hardening and fa vorable propagati on are two related but diffe rent properties. Generally speaking, a channel model can have both properties, o ne of them, 4 W ith correlated Rayleigh fading, a sufficient condition for asymptotic chann el hardening is that the spectral norm k R k k 2 of the channel co v ariance matrix remains bounded and β k = 1 M tr( R k ) remains strictly positiv e as M → ∞ . 5 For correlated Rayleigh fading channe ls, a sufficient condition for (7) is that the cov ariance matrices R i , R k hav e spectral norms tha t remain boun ded and the a verage channel gains β i = 1 M tr( R i ) an d β k = 1 M tr( R k ) rema in strictly positi ve as M → ∞ . 10 or no ne of them. The ke yhole channel that is st u died in [28] provides fa vorable prop agation, but not channel hardening. Uncorrelated R ayleigh fading with R k = β k I M satisfies both conditions and is often con s idered in the literature. In add i tion to uncorrelated Rayleigh fading, the two conditions are s at i sfied by a variety o f other channel models [7, Sec. 2.5] [27] such as correlated Rayleigh fading and line-of-sight (Lo S) with uniform ly random angles-of-arri val. C. Si gnal model The RA si g nal transmitted by UE k arrives at the BS wit h a specific carrier frequenc y of fset (CFO) ω k and a n ormalized (wit h respect t o the samp l ing period) timing misali g nment θ k . Follo wing [17], we assum e t hat ω k is within 2% of ∆ f such that its impact can reasonably be neglected if the RA block spans o n ly a fe w consecutive OFDM s ymbols [21]. On the other hand, tim ing errors { θ k } depend on the dis tances of UEs from the BS, and their maximu m value can reasonably be approximated as θ max = 2 D / ( cT s ) , where D is the boundary distance of the considered cell, T s = 1 / (∆ f N FFT ) is the sampli ng period and c = 3 × 10 8 m/s is the speed of light. A sim ple way to counteract the ef fects of { θ k } relies on the us e of a sufficiently long cyclic prefix comprising N G ≥ θ max + ∆ max sampling in t erv als, wi th ∆ max being the maximum expected delay spread within t he considered cell. 6 In d o ing so, t iming errors { θ k } only appear as ph ase shifts at the outp ut of t he receiv e dis crete Fourier transform (DFT) unit [21]. Notice that the presence of { θ k } d est roys the orthogonality amo n g t he frequency-domain cod es { f l k } and giv es rise to i n terference. Under the abov e ass umptions , i n the UL the DFT out put y ul m ( n ) ∈ C Q at antenn a m of the BS over subcarrier n during the Q OFDM symbols takes the form: y ul T m ( n ) = X k ∈K √ ρ k h k m Phase shift due to the timing error θ k z }| { e −  2 π N F F T nθ k e  2 π N nl k t T i k | {z } Intra-ce ll RA s ignals + i ul T m ( n ) | {z } Inter -cell inte rference due to acti ve UEs in the UL of all other cells + w T m ( n ) | {z } Noise (8) = X k ∈K √ ρ k h k m e  2 π nǫ k t T i k + i ul T m ( n ) + w T m ( n ) (9) 6 Note that such a solution is possible only for RA blocks. The CP of payload blocks must be made just greater than the channel length to minimize unnecessary overh ead. T his is why accurate timing estimates must be obtained during RA in order to avoid inter-block i nterference in the subsequent data transmission phase. 11 where we recall that K denotes the set of all UEs transmitting in the RA block, ǫ k = l k N − θ k N FFT | {z } Effe cti ve timing of fset of UE k (10) is the effective timing offset of UE k , w m ( n ) ∼ C N ( 0 Q , σ 2 I Q ) is the th ermal noise, and the vector i ul m ( n ) ∈ C Q accounts for the inter-ce ll interference generated in the UL by the activ e UEs i n all other cells. In writing the intra-cell RA signals in (8), we have assumed, wit hout any l oss of generality , that the first subcarrier of the consi d ered RA block has index 0 . For later con venience, let us denote Y ul m = [ y ul m (0) , . . . , y ul m ( N − 1)] T ∈ C N × Q the matrix collecting the DFT outputs at antenna m ov er the RA block, i .e. Y ul m = X k ∈K √ ρ k h k m c ( ǫ k ) t T i k + I ul m + W m (11) where c ( ǫ k ) = [1 , . . . , e  2 π ( N − 1) ǫ k ] T ∈ C N is the effecti ve frequency-domain code of UE k . As it is seen, the recei ved signal Y ul m depends on l k and θ k through the ef fective tim ing offset ǫ k . From (10) it follows that, in general, l k and θ k cannot be uni vocally determined from ǫ k . Ho wev er , under the assum ption that the maximum t i ming error θ max ≤ N F F T / N , the following result holds. Lemma 1. If θ max ≤ N F F T / N , then ǫ k in (10) can be univocally mapped i nto a si n gle pair ( l k , θ k ) as follows: l k = ceil  N ǫ k  (12) θ k = N F F T  l k N − ǫ k  . (13) Pr oof. T aking (10) into account and assuming θ max ≤ N F F T / N , one gets l k − 1 ≤ l k − N N F F T θ k ≤ N ǫ k ≤ l k (14) from which (12) and (13) are easily deri ved. Observe that the condi tion θ max ≤ N F F T / N in Lemma 1 is satisfied in practical scenarios. Consider , for example, a typical L TE system in which the subcarrier spacing i s ∆ f = 15 kHz and the DFT size i s N F F T = 1024 such that T s = 1 / (∆ f N FFT ) ≈ 65 . 1 ns. Recall also that θ max = 2 D max / ( cT s ) with D max being the cell boundary distance. Therefore, θ max = 2 D max / ( cT s ) ≤ N F F T / N is satisfied by D max ≤ 10 4 / N m . W i th N = 8 or 1 2 subcarriers per RA block, one gets 12 D max ≤ 555 m or D max ≤ 833 m. Both conditions are surely m et by future cellular networks for which cells of radius between 1 00 and 250 m ar e e xp ected in urban deplo yments. Therefore, i n the remainder we ass ume that t he result of Lemm a 1 holds t rue and thus th at a unique mappi ng exists b et w een ǫ k and ( l k , θ k ) . I I I . R A N D O M A C C E S S P RO C E D U R E In wh at fol lows, we show how the recei ved matrices { Y ul m } in (11) can be u s ed t o de velop an RA algo ri t hm that allo ws to detect the active codes { ( f l k , t i k ) : k ∈ K} , estim ate the tim ing offsets { θ k ; k ∈ K} of UEs’ signals, and resolve pos sible collisions by exploiting the l ar ge numb er M of antennas. In doing so, we exploit the fact that t he orthogonality o f the time-domain codes i s not destroyed 7 by the propagation channels and thus UEs that ha ve selected differe nt codes t i k do not in terfere with each other . W it hout loss of generality , we can t hus only focus on the subset of UEs that ha ve selected the same time-dom ain RA code and neglect the presence of the other UEs. This amounts to ass uming that t here is a sin g le t i me-domain RA cod e, and accordingly we can drop the i ndex i k to simplify the notat i on and e x position. In parti cular , we may re write (11) as Y ul m = K X k =1 √ ρ k h k m c ( ǫ k ) ! t T + I ul m + W m (15) where the sum is over the UEs sharing the same time-domain code t , whose number has been denoted by K . The RA procedure is designed by consi dering that any giv en UE k is identified by the triplet ( i k , l k , θ k ) . In particular , it de velops through t he following three steps. A. Step 1 - At the BS 1) Determinat ion of the nu mber o f UEs that ar e using the code t : The first problem is to determine the number of UEs that are transmitting in the RA blocks using the code t . F or this purpose, we start by c orrelating the receiv ed signal Y ul m in the time-dom ain w i th code t 8 , which amounts to computing vector z m = Y ul m t ∗ / || t || . B y taking (15) into account yields z m = Y ul m t ∗ || t || = K X k =1 h ′ k m c ( ǫ k ) + n m (16) 7 This is true only if the CFOs are relatively small (within 2% of subcarrier spacing ∆ f ) and the time-domain codes span only a few OF DM symbols. 8 W e stress that the BS performs such a correlation for all of the possible time-domain codes, but only the generic code t is considered here for simplicity . 13 where n m = ( I ul m + W m ) t ∗ / || t || ∈ C N , and h ′ k m = p ρ k Qh k m (17) denotes the ef fective channel of U E k at antenna m after time-domain d esp reading. From (16), it follows that z m has the same structure as the measurem ent model for a uniform linear array of passive sens ors in the p resence of m ultiple uncorrelated sources. W e can t h us identify the activ ated UEs and estim ate thei r corresponding ef fective timin g offsets by applying subspace- based methods [29]. T o see how this comes about, let u s compute the sampl e correlation matrix ˆ R z = 1 M P M m =1 z m z H m . By taking t h e limit M → ∞ and exploit ing the channel hardening and fa vorable propagation properties gi ven in (6) and (7), respectiv ely , yields ˆ R z ≍ R z = A ǫ + σ 2 I M (18) with A ǫ = K P k =1 ρ k Qβ k c ( ǫ k ) c H ( ǫ k ) . Let λ 1 ≥ λ 2 ≥ · · · ≥ λ N be the eigen values of R z arranged in non-increasing order . Then, from (18) it follows that λ j = µ j + σ 2 j = 1 , . . . , ra nk ( A ǫ ) (19) λ j = σ 2 j = rank ( A ǫ ) + 1 , . . . , N (20) where rank ( A ǫ ) ≤ K and µ j > 0 are, respectively , t he rank and t he non-zero eig env alues of the matrix A ǫ . Such a matrix is of rank rank ( A ǫ ) = K , iff ǫ k 6 = ǫ ℓ for ℓ 6 = k . Since the timing of fsets { θ k } are continuous random variables, from (10) it follows that the probabil ity that ǫ k = ǫ ℓ for ℓ 6 = k is equal t o zero, and hence rank ( A ǫ ) = K with probability 1. This means that, if R z were a vailable, all collisions could in principle be resolved provided that K ≤ N − 1 . In practice, howe ver , R z is not av ail able at t he BS and mus t be replaced with ˆ R z . The latter , howe ver , provides a good approximation o f R z when M i s suf ficiently large (as it is the case in Massive MIMO). Performing the EVD of ˆ R z and arranging the correspondin g eigen values ˆ λ 1 ≥ . . . ≥ ˆ λ N in non-i n creasing order , we can find an estimate of K through information-theoretic criteria. T wo prominent solutions in this sense are based on the Akaike and MDL criteria. Here, we adopt the MDL approach, which looks for the minimum of the follo wi ng objecti ve function [30 ]: ˆ K = arg min ℓ =0 ,...,N − 1 MDL ( ℓ ) = arg min ℓ =0 ,...,N − 1  1 2 ℓ (2 N − ℓ ) ln M − M ( N − ℓ ) ln ˆ g ( ℓ )  (21) 14 where ˆ g ( ℓ ) =  N Q n = ℓ +1 ˆ λ n  1 N − ℓ 1 N − ℓ N P n = ℓ +1 ˆ λ n (22) is the ratio between t he geometric and arithmetic means of { ˆ λ n ; n = ℓ + 1 , . . . , N } . In the remainder , we assume that M is suffi ciently lar ge such that ˆ K = K . Remark 1 (Asymptotic analysis of MDL est imator) . V arious works ana l yzed t he p erf ormance of the MDL est i mator (see for example [31]–[34]), which was pr oven to be strongly consis- tent [3 5], namely that lim M →∞ Pr ( ˆ K = K ) = 1 . F or fin i te M , it was ob s erved empiricall y that the main sour ce o f err or in the MDL estimator is und ere stimati o n of t he K sign als by e xactly one. F ollowing this observation , the authors in [31], [ 3 3], [34] studied the p roperties of ∆ MDL = MDL ( K − 1) − MDL ( K ) to show that , a s ymptotically as M → ∞ , ∆ MDL follows a Gaussian distribution with mean η MDL and standar d d eviation σ 2 MDL such that Pr  ˆ K 6 = K  can be appr oximat ed with Pr  ∆ MDL < 0  = Q  η MDL σ MDL  . Both η MDL and σ 2 MDL ar e given in explicit form in [34] as a function of N , K , an d the smallest eig en value of A ǫ . 2) Identificati on of the fr equency-domain codes { f l k } a nd estimation of the timing of fsets { θ k } : From (16), it foll ows that the observation space of z m can be decompo s ed into a signal subspace S , which is spanned by vectors { c ( ǫ k ) } , plus a noise subspace S such that any vec tor in S is orthog o nal to an y other one in S . Subspace-based methods like th e MUSIC (Multiple Signal Classification) [36] or ESPRIT [29] algorithms can be applied to the model (16 ) to find an estimate of { ǫ k : k ∈ K} . Com p ared to the MUSIC estimator [36], ESPRIT exhibits similar accurac y while dis pensing w i th any peak search procedure. A fundamental assu m ption behind both m eth ods estimator is that the dimension of the noise subspace S is at least unit ary . Thi s implies K < N , which means that the n umber of UEs selecting the same code t cannot e x ceed N − 1 . W e begin by arranging the eigen vectors of ˆ R z associated to the ˆ K largest eigen values into a matrix V = [ v 1 v 2 · · · v ˆ K ] ∈ C N × ˆ K . Th en, we app l y t he ESPRIT method t o (16) and retrieve the ef fectiv e t iming offsets in a decoupled fashion as ˆ ǫ ( j ) = arg { ψ j } 2 π j = 1 , . . . , ˆ K (23) 15 where { ψ 1 , . . . , ψ ˆ K } are the eigen values o f V =  V H 1 V 1  − 1 V H 1 V 2 , and the matrices V 1 and V 2 are obtained by collecting t he first and the l ast N − 1 rows of V , respectiv ely . Notice that the BS does not know which acti vated UEs the estimates { ˆ ǫ ( j ) } are associated to. This task will only be accomplished in step t ree of the RA process. In fact, we ha ve used the not ation ˆ ǫ ( j ) on purpose to emphasize that ˆ ǫ ( k ) is not i n general the estim ate of ǫ k . This is evident when ˆ K 6 = K , b ut it holds true even wh en ˆ K = K simply because the estimates provided by the ESPRIT algorithm are arranged arbitrarily . Howev er , there exists a bijective mappi ng (unknown at the BS) between the sets { ˆ ǫ ( j ) } and { ǫ k } . For simpl icity , we denote by j k the value of the index j corresponding to UE k such that ˆ ǫ ( j k ) → ǫ k . If θ max ≤ N F F T / N , from Lemma 1 we obtain ˆ l ( j k ) → l k and ˆ θ ( j k ) → θ k with ˆ l ( j k ) = ceil  N ˆ ǫ ( j k )  (24) ˆ θ ( j k ) = N F F T ˆ l ( j k ) N − ˆ ǫ ( j k ) ! . (25) For a giv en N and K < N , the estimation errors ˆ ǫ ( j k ) − ǫ k are asymptoticall y (e.g., M → ∞ ) jointly Gaussian distributed with zero mean and variances [37 ] V AR ( ˆ ǫ ( j k ) ) = 1 4 π 2 1 2 N SNR k 1 γ ( ǫ k ) 1 +  ( C H C ) − 1  k ,k SNR k ! (26) where γ ( ǫ k ) = c H ( ǫ k )  I N + C ( C H C ) − 1 C H  c ( ǫ k ) , C ∈ C N × K collects th e vectors { c ( ǫ k ) : j = 1 , . . . , K } , and SNR k = ρ k β k /σ 2 is the recei ved SNR. Notice that V AR (ˆ ǫ ( j k ) ) decreases monotically as N increases. When N → ∞ and SNR k takes relativ ely large values, (26) reduces to [37, Appendix G] V AR (ˆ ǫ ( j k ) ) → N →∞ 1 4 π 2 6 N 3 M 1 SNR k . (27) According to [37 ], the ESPRI T algorithm is un l ikely to resolve signals for wh i ch 8 p V AR ( ˆ ǫ ( j k ) ) ≥ min j k ,j i   ǫ ( j k ) − ǫ ( j i )   . This m eans that, if th e effe ctive timing of fsets ǫ k and ǫ i of two (or more) UEs are di f ferent (i.e., ǫ k 6 = ǫ i ) but such t hat   ǫ k − ǫ i   is sm aller than the resolution provided by the ESPRIT algorithm , then the two UEs are undistingui shable. Not ice that the ESPRIT resolution increases cubi call y with N and l inearly with M ; howe ver , an infinite resol ution is achiev ed only if bo th M and N grow to infinity . In the s equ el, we assume that N and M are large enough such that th e ESPRIT algorit hm is able to resolve UEs for which ǫ k 6 = ǫ i . Nu merical resul ts will be used in Section V to v al i date t he impact of the finite resolution of the ESP RIT algorithm . 16 3) Channel estimat ion: The estimates { ( ˆ l ( j ) , ˆ θ ( j ) ) : j = 1 , . . . , ˆ K } are used to acquire infor- mation abo ut the correspondi n g chann el vectors. From (16), the LS estimate o f th e channel gain h ′ ( j ) m associated to the pair ( ˆ l ( j ) , ˆ θ ( j ) ) is found to be [38] ˆ h ′ ( j ) m = e T j  ˆ C H ˆ C  − 1 ˆ C H z m j = 1 , . . . , ˆ K (28) where e j denotes the j th component of the c anonical basis, and ˆ C ∈ C N × ˆ K collects the vectors { c (ˆ ǫ ( j ) ) : j = 1 , . . . , ˆ K } . W e observe that ˆ C is a V andermonde matrix, so that the full-rank condition, needed for the computation o f ( ˆ C H ˆ C ) − 1 , is m et i f and only if ˆ ǫ ( j ) 6 = ˆ ǫ ( j ′ ) ∀ j 6 = j ′ . This happens with probability one since the ES PRIT algorithm provides ˆ K dis t inct estimates. Under the a ssumpt ion that ˆ K = K and the ef fective timing of fsets ǫ k are perfectly e stimated (i.e. ˆ ǫ ( j k ) = ǫ k ), by plugging (16) into (21) yields: ˆ h ′ ( j k ) m = h ′ k m + η ( j k ) m (29) where η ( j k ) m = e T j k  ˆ C H ˆ C  − 1 ˆ C H n m . In matrix form, (29) can be re written as ˆ h ′ ( j k ) = h ′ k + η ( j k ) (30) with η ( j k ) = [ η ( j k )1 , . . . , η ( j k ) M ] T . W e emphasi ze that (30) hold s true only wh en ˆ K = K and ˆ ǫ ( j k ) = ǫ k , and such conditions can be met only when both M and N tend t o infinity . B. Step 2 (SUCR) The estim at ed channel vectors ˆ h ′ ( j ) = [ ˆ h ′ ( j )1 , . . . , ˆ h ′ ( j ) M ] T with j = 1 , . . . , ˆ K , are used by the BS in Step 2 to respond to the possibly identified UEs by s end ing a DL precoded version of th e frequency- and tim e-domain codes. The DL precoded matrix V ∈ C M N × Q over the RA block from all transmit antennas is 9 V = √ ρ dl ˆ K X j =1 ˆ h ′ ( j ) || ˆ h ′ ( j ) || ⊗  f ˆ l ( j ) t T  (31) where ρ dl > 0 denotes th e DL transmit power . The DL p recoded matrix V is transmit ted in a multicast fashion by the BS, and i s exploited by Step 3 of the SUCR protocol (a t each activ ated UE) as shown next. 9 The signal transmitted by the BS is actually obtained as the sum of signals like (31); that is, one for each detected time-domain code. A code is detected when the MDL algorithm estimates that at least one UE is using t hat code. 17 The recei ved signal R dl k ∈ C N × Q over the RA block at UE k is R dl k = √ ρ dl ˆ K X j =1 h H k ˆ h ′ ( j ) || ˆ h ′ ( j ) || f ˆ l ( j ) t T + I dl k + W dl k (32) where I dl k ∈ C N × Q accounts for the i n ter -cell interference i n the DL receiv ed from all other cells at UE k and W dl k ∈ C N × Q is the recei ver noise matrix. The recei ved si gnal R dl k is used by UE k to implem ent the SUCR protocol proposed in [10], which allows to resolve possible collisions and to enable retra nsmiss ion of detected UEs. By correlating the recei ved sig nal with its selected (and normalized) random codes f l k and t , UE k gets r dl k = f H l k || f l k || R dl k t ∗ || t || . (33) W e assume that (29) holds true. This requires ˆ K = K and ˆ ǫ ( j k ) = ǫ k , which is achieve d for M → ∞ and N → ∞ . By normalizing r dl k with √ M and taking the limit M → ∞ yields r dl k √ M ≍ √ ρ dl τ √ ρ k Qβ k √ α k (34) where α k ≍ 1 M || ˆ h ′ ( j k ) || 2 (35) and we ha ve used the prop erty 1 M h H k ˆ h ′ ( j k ) ≍ p ρ k Qβ k . (36) In writing (3 4), w e ha ve ass umed (as in [10] ) that th e inter -cell i n terference does n ot scale with M s o that the noise-plus -interference term in (33) (after t he normalizati on by √ M ) v anish es as M → ∞ 1 √ M · f H l k || f l k || ( I dl m + W dl m ) t ∗ || t || ≍ 0 . (37) Based on the approach i n [10], we propose that UE k appli es t he following rule to d ecide whether to reply or not to the DL RA signal transmitted b y the BS: R k : ρ k β k Q > ˆ α k / 2 + ǫ k ( Rep eat ) (38) W k : ρ k β k Q ≤ ˆ α k / 2 + ǫ k ( W ait and sta rt o v er ) (39) where ˆ α k is an estimate of α k giv en by ˆ α k = max M ρ dl τ ρ k β 2 k Q ( ℜ e( r dl k )) 2 − σ 2 , ρ k β k Q ! (40) 18 and ǫ k is a bias parameter that can be used to t une the syst em beha v i or to the final performance criterion [10 ]. Specifica lly , i f W k is true, UE k will pick up ne w RA codes and retransmit after a random w aiting time. On the other hand, i f R k is true, it n otifies the BS by retransmitting the code f l k t T , follo wed by an UL message that contains the uni que identity number of the UE. C. St ep 3 The BS receives the pilot codes from each UE that decided in fa vour of the repeat hypothesis in the pre vious s tep. The recei ved si gnal Y ul m ∈ C N × Q at antenn a m ov er the RA block take s the form Y ul m = X k ∈R √ ρ k , ul h k m c ( ǫ k ) t T + I ul m + W m (41) where the elements of R are the indices of U E s for which R k is true, and ρ k , ul is giv en by ρ k , ul = 1 ρ k ˆ α k τ Qβ 2 k . (42) Notice that the comp u tation of ρ k , ul at UE k requires knowledge of the large-sca le fading coef ficient β k . This information can be acquired by the UE on the basis of the DL control channel [17]. The receiv ed signal Y ul m is first correlated with the detected (ef fectiv e) frequency- and time- domain sequences { c (ˆ ǫ ( j ) ) , t } yielding Z ul ( j ) m = c H (ˆ ǫ ( j ) ) || c ( ˆ ǫ ( j ) ) || Y ul m t ∗ || t || . (43) The correlatio n with the effecti ve code c ( ˆ ǫ ( j ) ) all ows the BS to discriminate UEs on the basis of both t h e selected codes and the tim ing offsets. By correlating Z ul ( j ) = [ Z ul ( j )1 , . . . , Z ul ( j ) M ] T ∈ C M with the corresponding estimated channel vec tor ˆ h ′ ( j ) produces r ul ( j ) = ( ˆ h ′ ( j ) ) H || ˆ h ′ ( j ) || Z ul ( j ) . (44) Under the assu mption that M and N are su fficiently lar ge, we hav e ˆ K = K and ˆ ǫ ( j k ) = ǫ k such that, by taking (41) and (43) – (44) i nto account, we obtain r ul ( j k ) = √ ρ k , ul τ ( ˆ h ′ ( j k ) ) H h k || ˆ h ′ ( j k ) || + X ν ∈R ,ν 6 = k p Qρ ν, ul ( ˆ h ′ ( j k ) ) H h ν || ˆ h ′ ( j k ) || c H (ˆ ǫ ( j k ) ) c ( ǫ ν ) k c (ˆ ǫ ( j k ) ) k + ξ k (45) 19 Fig. 3: Propo sed RA p rotocol for Massi ve MIMO. Unlike the SUCR pro tocol reported in Fig. 1(b), the propo sed proced u re aims at d etecting the n umber o f active codes thr o ugh th e MDL alg o rithm an d , at the sam e time, perfo rming timing estimation b y means of the ESPRIT algo rithm. T iming estimates are exploited to comp ute the LS estimate of the c hannels of all detected codes. where ξ k accounts for the interference and noi se terms in (41). By norm alizing r ul ( j k ) with √ M and taking the limit M → ∞ , yields r ul ( j k ) √ M ( a ) ≍ √ ρ k , ul τ √ ρ k Qβ k √ ˆ α k ( b ) = 1 (46) where ( a ) follows from (34) and from 1 √ M ( ˆ h ′ ( j k ) ) H h ν || ˆ h ′ ( j k ) || ≍ 0 for ν 6 = k (47) whereas ( b ) is due to (42). Based on (46), the BS adopts t h e following rule to decide whether there is or not an acti ve UE in the RA block associated to the inde x j = j k or , equiv alently , t o the pair ( ˆ l ( j k ) , ˆ θ ( j k ) ) : ( ˆ l ( j k ) , ˆ θ ( j k ) ) is declared as    Detected Undetected if δ 1 < r ul ( j k ) √ M < δ 2 otherwise (48) where the thresholds δ 1 < 1 and δ 2 > 1 shoul d be properly designed to tune the system behavior to the final performance crit eria; for example, to maximize the a verage nu m ber of resolved 20 Algorithm 1: Th e proposed RA protocol 1: Com p ute z m in (16) for m = 1 , . . . , M . # Step 1 2: Com p ute the SVD of the sample correlation matrix ˆ R z = 1 M P M m =1 z m z H m . 3: Com p ute ˆ K throug h the MDL alg o rithm in (21). 4: Com p ute { ˆ ǫ ( j ) ; j = 1 , . . . , ˆ K } by applying the ESPRIT algorithm to (16). 5: Us e { ˆ ǫ ( j ) ; j = 1 , . . . , ˆ K } in (28) to obtain LS channel estimates. 6: BS u ses the LS channel estimates to send t he precoded signal (31). # Step 2 7: Each UE correlates the received sign al with its selected codes as in (33). 8: Each UE di stributi vely computes (40) and decides wh eth er to reply or not according to (38) and (39). 9: Com p ute (43) for each detected pair of codes. # Step 3 10: Compu t e (44) by correlating with the corre sponding LS channel estimate. 11: Use (48) to decide whether there is or not an acti ve UE. collisions or to minim i ze the risk of false p ositives (or negativ es). Once a pair ( ˆ l ( j k ) , ˆ θ ( j k ) ) is declared as detected, t he BS proceeds recovering the unique ident ity num ber , cont ai n ed i n t he recei ved signal, and uses i t to perform autho ri zati on and registration of the ass o ciated UE. Then, the BS broadcasts a DL response message indicating wh ich UEs have been detected and giving the correspondi n g instructions for timing adjustment. Those UEs that do not recei ve the notification wi ll pick up new RA codes and retransmi t after a random waiting time. Thi s is don e until success notification. The steps t h rough which th e proposed RA operates are reported in Fig. 3 and also in Algorithm 1. Unlike t he SUCR protocol reported in Fig. 1 (b), t h e additional blocks 10 allow to estimate the t iming of fsets and to inherently exploit them to improve the detection capabilities of the protocol it s elf. This latter poi nt is di scussed furth er in the next section. W e conclud e by recalling (as mentioned in the Introduction) that the SUCR protocol [10] cannot be applied in the presence of timing offsets { θ k } since the orthogo n ality among the frequency-domain codes would be destroyed. Such a loss of orthogonality g ives rise to interference, which highly degra des the detection performance of the protocol itself. 10 Notice that no signaling is exchan ged between the BS and UE s, ex cept for in S t ep 4 where the BS broadcasts a DL response message indicating which UEs have been detected and giv ing the corresponding instructions f or timing adjustment. 21 I V . C A S E S T U D Y - T W O C O L L I D I N G U E S The rati onale behind the proposed RA protocol relies on the assumptio n that M and N are suffi ciently large s uch that ˆ K = K and ˆ ǫ j k = ǫ k . While the asymptotic re gi m e M → ∞ can be virtually achiev ed in Massive MIMO, the conditio n N → ∞ is not granted due to the limitation s imposed by the coherence bandwidth of the propagation channel. In Section V, the performance of the p rop osed RA protocol wi ll be in vestigated by means of nu merical results for p ractical values of M and N . In order to understand the ef fect o f a finite resolut ion o f the ESPRIT algorithm, let us cons ider the following simple case s tudy . Ass ume that th ere are K = 2 UEs, which have s elected the same time- and frequency-domain codes, namely , t and f l and wh i ch are also characterized by the same timing offset θ . The latt er assumption adequately model s a practical si tuation in which the t wo UEs ha ve sli ghtly different ti ming offsets θ 1 and θ 2 , such that the quantity | θ 1 − θ 2 | is m uch smaller than t he resolution provided by the ESPRIT algorithm. Accordingly , the two UEs are approximately seen as a single UE with a sing l e ti ming of fset. The ESPRIT algo rithm provides an estimate ˆ ǫ of ǫ = l / N − θ / N F F T , which is first used for the computation of ˆ l and ˆ θ through (24) and (25), and then by the channel estimation algorithm. By using (28) – (30), one gets ˆ h ′ = κ h ′ + η (49) where κ = e j π ( N − 1)( ǫ − ˆ ǫ ) sin( π N ( ǫ − ˆ ǫ ) ) / [ N sin( π ( ǫ − ˆ ǫ ))] , and h ′ is the ef fective composite channel giv en by h ′ = p Q ( √ ρ 1 h 1 + √ ρ 2 h 2 ) (50) with ( ρ 1 , ρ 2 ), and ( h 1 , h 2 ) being, respectively , t he powers and the channels of UE 1 and UE 2. During Step 2, by using (32) – ( 33) UE k com putes r dl k = √ ρ dl τ h H k ˆ h ′ || ˆ h ′ || + ζ k k = 1 , 2 (51) where ζ k results from interference and nois e. In the asy m ptotic regime ( M → ∞ ), r dl k can be approximated as follows: r dl k ≍ √ M √ ρ dl τ √ ρ k β k √ ρ 1 β 1 + ρ 2 β 2 k = 1 , 2 (52) from which, by plugging (50) into (40), one gets ˆ α 1 = ˆ α 2 ≍ ( ρ 1 β 1 + ρ 2 β 2 ) Q. (53) 22 By using t he asymptotic resul t (53) into (38) and assuming ρ 1 β 1 6 = ρ 2 β 2 , it fol l ows that only the strong est UE would retransmit to the BS after Step 2 (as it is expected from the application of the SUCR algorithm) and thus it would be detected as explained above during Step 3. On the other hand, when ρ 1 β 1 ≈ ρ 2 β 2 and M is not s u fficiently large, i t may happen that both UEs decide to retransmit their RA codes in response to the DL s ignal from t he BS. In s uch a case, (45) reduces to r ul = κ √ τ ( ˆ h ′ ) H || ˆ h ′ ||  √ ρ 1 , ul h 1 + √ ρ 2 , ul h 2  + ξ (54) where ρ 1 , ul and ρ 2 , ul are computed according to (42). By t aki ng the lim it M → ∞ into (54) yields r ul ≍ √ τ √ M β 1 √ ρ 1 ρ 1 , ul + β 2 √ ρ 2 ρ 2 , ul √ ρ 1 β 1 + ρ 2 β 2 ( a ) = 2 (55) where ( a ) follows from (42) and (53). In this situation, there is no way for the BS to distingui sh between UE 1 and UE 2, and hence it must dis card all the signals associated wit h the pair ( ˆ l, ˆ θ ) . This is only possible if we set δ 2 < 2 in (48). Th e above reasoni ng can straightforwardly be extended to the case of more than two UEs sharing t h e same pair of codes and ha ve (nearly) the same timing of fsets . V . N U M E R I C A L R E S U LT S Numerical result s are used t o assess the p erformance of the proposed RA protocol. W e consider a cellul ar network op eratin g over a bandwidth B = 20 MHz and composed of 9 cells distributed on a regular grid with an in ter -site distance of 500 m; each cell cov ers a squ are area centered at the BS w i th side l ength D = 50 0 m. The DFT s ize is N FFT = 1 024 and the no ise power is σ 2 = − 97 . 8 dBm. The UE density is µ = 10 4 UE/km 2 (which corresponds t o |I | = 2500 ) and UEs are uni formly distributed in each cell at locations further than 25 m from the serving BS. W e denote by d k the distance of UE k from its own BS. The RA block is composed of Q = 2 consecutiv e OFDM symbols (such that the impact of th e residual C FO errors is negligible) and N = 8 o r 12 adjacent subcarriers. A W alsh-Hadamard codeboo k is u s ed in t he time-domain whereas frequency-domain codes belong to the Fourier basis. Unless otherwise specified, each UE decides to access the n et work with prob ability p A = 1 % , meaning that 2 5 UEs on average try to ent er t he network and the probability of having a collision is around 0 . 47 with N = 8 and 0 . 27 with N = 12 . The timi ng error θ k of UE k is com puted on the basis of its dis tance d k from the BS as θ k = round (2 d k B /c ) where c = 3 × 10 8 m/s is the speed of li ght. Accordingly , 23 T ABLE I: Network and system par ameters Parameter V alue Network layout Square pattern Number of cells L = 9 Cell area 500 × 500 m 2 Bandwidth B = 20 MHz DFT size N FFT = 1024 UE density µ = 10 4 UE/km 2 Probability of acti v ation p A = 0 . 5% , 1% , 2% W alsh-Hadamard ti me-domain codes Q = 2 Fourier frequenc y-domain codes N = 8 , 12 Minimum RA po wer ρ min = 100 mW Maximum RA power ρ max = 1 W DL transmit po wer ρ dl = 1 W the maximum timing error i s θ max = max k θ k = round ( √ 2 D B /c ) = 47 s amples and is achie ved by a UE positio ned in the cell corner at a dis tance of √ 2 / (2 D ) . Not ice t h at θ max satisfies the condition in Lem ma 1 for both N = 8 and N = 12 . W e assume that the m i nimum and m axi mum power l e vels during the RA procedure are ρ min = 100 mW and ρ max = 1 W , respectively . T o emulate a netw ork with UEs that ha ve made di f ferent attempts in the RA procedure, the power lev el ρ k of UE k in (2) is s elected with uni form probability from the set [∆ 0 , ∆ 1 , . . . , ∆ 10 ] with ∆ i = ρ min e i ∆ and ∆ = 0 . 1 ln( ρ max /ρ min ) . 11 The DL and UL transmit po wers ρ dl and ρ k , ul in (31) and (42) of Step 2 are r espective ly set to ρ dl = 1 W and ρ k , ul = min  1 ρ k ˆ α k τ Qβ 2 k , ρ max  . (56) All the above parameter values are summ arized in T able I. The performance of the proposed RA procedure are measured in t erms of the a verage number of codes declared as Detecte d over a giv en time-domain code on the basi s o f (48) in Step 3. In d o ing thi s, we rest ri ct t o those UEs for which the recei ved SNR SNR k = ρ k β k /σ 2 is larger than 5 dB. The results are ob tained av eraging over 1000 different channel realizations and UE positions . T wo channel models are considered. The first one is uncorrelated Rayleigh fading and is such t h at h k ∼ C N ( 0 , β k I M ) where β k is the path l oss function obtained as β k = Ω d − κ k where κ = 3 . 7 i s the path loss exponent and Ω = − 148 . 1 dB i s the path loss at a reference 11 This choice al l o ws to emulate UEs that retransmit, if not admitted, by ex ponentially i ncreasing their transmit po wers. 24 distance of 1 km. The second one is correlated Rayleigh fading with h k ∼ C N ( 0 , R k ) , where we consider a uniform linear array (ULA) at the BS modeled by the exponential correlation model wit h correlatio n factor r between adjacent antennas, a verage lar ge-scale fading β k , and angle-of-arri val φ k [39]. This leads to [ R k ] m,n = β k r | n − m | e j φ k ( n − m ) . (57) Both cases are considered: i ) the adjacent cells are si l ent during the R A procedure (i.e., without inter-ce ll interference); ii ) they perform regular data transmis sion (i.e., with inter -cell interfer - ence). In the latt er case, we assum e that there are ten active UEs i n each of th e neighboring cells and the propagation channels are modeled as un correlated Rayleigh fading (using the same power lev els and path loss m odels as above). Following [10], the a verage UL i nterference ¯ ω = E n || I ul m t ∗ || t || || 2 / M o is assumed to be k nown at the UE (it is the same for all UEs) and is subtracted from ˆ α k / 2 by setting ǫ k = − ¯ ω / 2 . Comparisons are made with a “baseline” procedure wherein UE s are detected by the BS independently of their power le vels (and hence independently of the SNRs) w h ene ver d iffe rent codes are selected. Therefore, the p rob ability of successfully detecting a given UE with s uch a “baseline” procedure coincides with t he probability that a given code is selected by a UE on l y (under the assumption that the number of ac tiv ated UEs is not zero), which is giv en by |I | p A QN  1 − p A QN  |I |− 1 1 −  1 − p A QN  |I | . (58) A. Impact of timing of fsets W e b egin by in vestigati ng to what extent tim ing of fsets improves t h e detection capabilities of the RA procedure. T o this end, we assume that two UEs (amon g the total number of activ ated UEs) with timing offsets θ 1 and θ 2 hav e selected the same code f j t T i . While the ti m ing of fsets of all other UEs entering the network are computed as described above (as a function o f distances), we assu me for simpli city that θ 1 = 0 whereas θ 2 var ies from 0 to 32 . Fig. 4 plots the av erage number of UEs declared as Detected ov er f j t T i as a function of ∆ θ = θ 2 when Q = 2 , N = 8 and M = 100 . Th ree dif ferent values of p A are considered. The case p A = 2% correspon d s to a high-overloaded scenario in th e sense that th e avera ge number of UEs entering the network, which is equal t o |I | p A = 5 0 , is greater (more than double) than all possi ble codes QN = 16 . 25 0 4 8 12 16 20 24 28 32 0 0.5 1 1.5 2 T iming offset differen ce, ∆ θ A verage number of detected U E s p A = 0 . 5% p A = 1% p A = 2% Fig. 4 : A verage number of detected UEs vs. ∆ θ for a given cod e f j t T i . In particular, we assume that the code f j t T i has been selected by two UEs (among the total nu mber of activ ated UEs) with timing o ffsets θ 1 = 0 and θ 2 = ∆ θ . W e assume that N = 8 , M = 10 0 an d that th e pro bability of acti vation is p A = 0 . 5 , 1% and 2% . Uncorrelated Rayleigh fading is considered. The other two cases p A = 0 . 5% and 1% can b e considered as low- and medium-overloaded scenarios (wi th |I | p A = 12 . 5 and |I | p A = 25 , respecti vely). As anticipated in Section IV, the results of Fig. 4 show th at the detection capabilities of the RA proto col improves as ∆ θ gets lar ger . F or p A = 0 . 5% and 1% , the a verage number of detected U E s is larger than 1 already for ∆ θ > 4 and ∆ θ > 12 samples, respectively . This proves that, when the i n herent timing offsets are suffi ciently dif ferent, the RA procedure is able to resolve the two UEs, though both have selected t he s ame code. For a high-overloaded scenario wit h p A = 2% , the ave rage nu mber of detected UEs is approximately 0 . 75 and increases very slo wly wit h ∆ θ . This is due to the high interference created by the other U E s entering the network. Notice that if two UEs s elect th e same code f j t T i and the SUCR protocol in [10] is used, then at mo st one of them can be declared as Detected . This is beca use the SUCR allows only the strongest between the two t o retransmit its code to the BS. B. P erfor mance e valuation Fig. 5(a) plots the p robability that a UE is declared as Detecte d as a function of M when N = 8 or 12 . Uncorrelated Rayleigh fading is assumed. As expected, addi ng more antennas at the BS improves the RA p erformance in all cases, but at a slow p ace for M > 50 . W ith M = 100 , the probabili ty of success without inter -cell interference is 0 . 75 wit h N = 8 and 0 . 83 with N = 12 . W ith in ter- cell interference, it reduces to 0 . 66 and 0 . 76 , which is sti ll relativ ely hi g h 26 10 50 100 150 200 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Number of BS antennas, M Probability of detecting a giv en UE N = 8 N = 12 w/o inte rcell inte rf. w intercell inte rf. baselin e (a) Probability that a giv en UE is declared as Detected i n (48). 10 50 100 150 200 2 4 6 8 10 12 14 16 Number of BS antennas, M RMSE of timing estimates N = 8 N = 12 w/o inte rcell inte rf. w intercell inte rf. (b) Timing estimation accuracy of the detected UE. Fig. 5: Perfor m ance of the RA proposed pr otocol vs. number of BS anten nas fo r a giv en UE when Q = 2 and N = 8 , 1 6 with p A = 1% . Un correlated Rayleigh fading with and without inter c ell interferenc e is assumed. taking i n to account that the av erage num ber of activ ated UEs is 25 whil e the number of t i me- frequency codes is 16 and 24 . F or the considered set u p, the “baseline” approach yields a l ower detection probability . Specifically , it provides 0 . 53 and 0 . 71 with N = 8 and 1 2 , respectively . Notice that for th e “baselin e” system the resul ts are highly optim istic s ince it has been assumed that, if a RA code is selected by a single UE, this UE is detected by the BS independently o f its power level (and h ence independently o f the SNR). On the other hand, for the proposed RA protocol the results in Fig. 5(a) take into account the power l e vels of t h e diff erent UEs trying to access the network as well as the interference coming from other ce lls. Fig. 5(b) ill ustrates the 27 T ABLE II: Probability to resolve a collision betwe en UEs that h ave selected the same pair of cod es with Q = 2 , N = 8 , p A = 1% , M = 100 and uncorrelated Rayleigh fading. Setup Probability to resolve a collision N = 8 w/o intercell interf. 0 . 81 N = 12 w/o intercell interf. 0 . 91 N = 8 w intercell interf. 0 . 55 N = 12 w intercell interf. 0 . 65 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 M = 100 M = 50 w/o intercell interf . w intercel l interf. baselin e Fig. 6 : Performance of the pr oposed RA protocol when the correlated Ray leigh fading m odel in (57) is considered with Q = 2 , N = 8 , p A = 1% and M = 50 or 10 0 . Both cases with and withou t in tercell interfer ence ar e considered . root mean-square-error (RMSE) of the tim i ng estimates ve rsus M for N = 8 and 12 . In both cases (wit h and wi thout inter-cell int erference), the results show that the RMSE decreases fast as M grows large, and it is s maller than a fe w sampling interva ls for M > 50 with both N = 8 and 1 2 . This provides evidence of the f act that, un like existing s o lutions, the proposed protocol allows to comput e reliable estimates of t he timing of fsets. T o further hig h light the capability of the proposed RA procedure in identifyi ng UEs that hav e selected the sam e pair of tim e- and frequency-domain codes (by exploiting timin g mi s- alignments), T abl e II report s the probability t o resolve collisions with N = 8 and 12 (wit h and without intercell i nterference) under t he condition that t wo or more UEs (sharin g the same codes) reply to the DL RA signal transmitted by the BS (according to (38) and (39)). As seen, with N = 8 collisions are resolved with probability 0 . 55 and 0 . 81 with and wi t hout intercell interference, respecti vely . 28 0.005 0.01 0.015 0.02 0 1 2 3 4 5 6 7 8 9 10 Probability of acti v ation, p A A verage number of attempts, η A N = 8 N = 12 w/o intercell interf . w intercel l interf. baselin e Fig. 7: A verag e n umber of attempts req uired by a given UE to be successfully detected with with Q = 2 , N = 8 or 12 , an d M = 100 . Both cases with and without in tercell interference ar e considered with unco rrelated Rayleigh fading. Fig. 6 eva luates the i mpact of correlation at BS antennas when the exponential correlation model in (57) i s used with M = 50 and M = 100 . Fig. 6 shows t hat with M = 1 00 the detection probability is marginally af fected by v alues o f the correlation factor r up to 0 . 8 for both cases (with and without int ercell interference). If M is reduced to 50 , then the performance deteriorates as soon as r ≥ 0 . 6 . This is because with M = 50 and r = 0 . 6 the number of independent observ ations b ecom es o n the same order of N and thus the estim ation accuracy of the sample correlation matrix ˆ R z decreases. Numerical results (not reported for s pace limit at i ons) show that the R MSE of timin g e stimates ke eps constant for all the considered v al u es of r . This makes the proposed RA protocol well suited for both uncorrelated and correlated propagation channels. The main purpos e of an RA prot ocol is that eve ry UE shou l d be admitt ed t o data transmis sion after as few RA attempts as possible. Fig. 7 shows the av erage number of RA attempts, η A , that each UE makes as a function of p A , wi th N = 8 and 12 , and in bot h cases with and without intercell interference. Uncorrelated Rayleigh fading is consi dered. As expected, η A increases as p A grows. W ith p A = 1% , 2 ≤ η A ≤ 3 attem p ts are required for all the in vestigated scenarios. W ith the baseline procedure, η A rapidly increases with p A (this is particularly evident wi th N = 8 ), and a significantly larger number of retransmissions is required compared to t he proposed RA protocol. 29 T ABLE III: Comp utational complexity of Steps 1 and 3 Number of complex multiplications and divisions Step 1 Q  K 2 + K 2 ( N − 1) + K 2 ( N − 1) + K 3 + K 3 − K 3 + M  N + N 2 + N 2 + K 2 N + K N 2 + K 2 + K 3 − K 3   Step 3 Q ( M K ( N + N Q ) + M K ) C. Complexity analysis As illustrated in Section III, the proposed RA procedure operates through three steps o f which Step 2 i s exactly the same as Step 2 of the SUCR prot ocol [10]. The additional complexity of Steps 1 and 3 is assessed in t erms of complex mul t iplications and divisions as follows. 12 In Step 1 for each given ti me-domain code, e valuating z ( m ) in (16) for m = 1 , . . . , M and requires M N complex multipli cations while the complexity in volved in the computation of ˆ R z is approxi mately M ( N 2 + N ) / 2 . The com p u tation of the eigen vectors of ˆ R z requires N 3 operations w h ereas e valuating V =  V H 1 V 1  − 1 V H 1 V 2 needs 13 approximately K 2 + K 2 ( N − 1) + K 2 ( N − 1) + K 3 + K 3 − K 3 under the assumption that ˆ K = K . T h e computational burden of channel estimati on i s M ( K 2 N + K N 2 + K 2 + K 3 − K 3 ) . In Step 3 for each t ime-domain code, the computation of { Z ul ( j ) m } in (43) for m = 1 , . . . , M requires M K ( N + N Q ) complex multipl i cations whereas M K mu ltiplication s are required for { r ul ( j ) } in (44). The number of com plex operatio n s required by th e two steps i s reported in T able III. As we can see, it scales li nearly with M for both. On the other hand, it increases as N 2 for Step 1 and as N for Step 3. Als o, the functi onal dependence wi th respect to K is linear only for Step 3 while i t is cubi c for Step 1. Fig. 8 illustrates the num ber of complex operations as a function of M with N = 8 and N = 12 when the number of active UEs is alw ays K = |I | p A . As expected, Step 1 has t he hi ghest com p lexity . W ith M = 1 00 , pass i ng from N = 8 and N = 12 increases the complexity of Steps 1 and 3 by a factor 1 . 36 and 1 . 48 , respectiv ely . Note that the additi o nal complexity due t o Steps 1 and 3 is the price to pay for detecting UEs while performing timing estimation wi th high accurac y . 12 Consider the matrices A ∈ C N 1 × N 2 and B ∈ C N 2 × N 3 . The matr i x-matrix multiplication AB requires N 1 N 2 N 3 complex multiplications. The multiplication AA H only requires N 2 1 + N 1 2 N 2 complex multiplicati ons, by utilizing the Hermitian symmetry . 13 Consider the Hermitian positiv e semi-definite matrix A ∈ C N 1 × N 1 and the matri x B ∈ C N 1 × N 2 . The LDL H decomposition of A can be com puted using N 3 1 − N 1 3 complex multiplications. The ma trix A − 1 B can be computed using N 2 1 N 2 complex multiplications and N 1 complex di visions if the LDL H decomposition of A is known. 30 10 50 100 150 200 10 4 10 5 10 6 10 7 Number of antennas, M Number of comple x operations Step 1, N = 8 Step 3, N = 8 Step 1, N = 12 Step 3, N = 12 Fig. 8: Numb e r of com plex ope r ations ( multiplication s an d divisions) per Steps 1 and 3, when N = 8 an d N = 12 . V I . C O N C L U S I O N S W e proposed an RA procedure for initial access and handover in the uplink of Massive MIMO systems. Each UE that wants to enter the network random ly selects a pair of predefined RA codes and perform spreading ov er the RA block in b oth the frequency and time domains. By exploiting t h e fav orable propagation conditions offered by Massiv e MIMO systems as well as the inherent different ti m e offsets in the reception of uplin k s ignals, t he proposed RA procedure aimed at detecting the largest numb er of UEs and, at the same t i me, performing timing est imation. Numerical results showed t h at a few tens of antennas are enough t o succe ssfully detect a given UE, while providing reliable timing estim ates (smaller than a few sampling int erv als). W ith 2 . 5 × 10 3 UEs that may simult aneously become acti ve with probability 1% and a total of 16 frequency-time codes (in a given random access b l ock), it turns o ut that, with 100 antennas, the proposed procedure successfully detects a gi ven UE with probabi lity 7 5 % wh ile providing reliable timing estimates. The p rice to pay for all this is a certain increase of sys tem complexity that scales li n early with th e number of antennas and quadratically with the lengt h of frequency- domain codes. R E F E R E N C E S [1] T . Marzetta, “Noncoop erativ e cellular wireless with unlimited numbers of base station antennas, ” IEEE T ransac tions on W ir eless C ommunications , vol. 9, no. 11, pp. 3590 – 3600, Nov . 2010. [2] E . Larsson, O. Edfors, F . Tu fvesson, and T . Marzetta, “Massiv e MIMO for next generation wireless systems, ” IEEE Communications Magazine , v ol. 52, no. 2, pp. 186 – 195, February 2014. 31 [3] J. Hoydis, S. ten B rink, and M. Debbah, “Massi ve MIMO in the UL/DL of cellular networks: Ho w many antennas do we need?” IEE E J. Sel. Ar eas Commun. , vo l. 31, no. 2, pp. 160–171 , Feb . 2013. [4] T . L. Marzetta, “Massiv e MIMO: An i ntroduction, ” Bell Labs T echnical Jou rnal , vol. 20, pp. 11–22, 2015 . [5] B j ¨ ornson, E. G. Larsson, and T . L. Marzetta, “Massi ve MIMO: ten myths and one crit i cal question, ” IEEE Communications Maga zine , vol. 54, no. 2, pp. 114–123, February 2016. [6] E . Bj ¨ ornson, J. H oydis, and L. Sanguinetti, “Mass iv e MIMO has unlimited capa city , ” IEEE T ransaction s on W ireless Communications , vol. 17, no. 1, pp. 574 –590, Jan 2018. [7] E . Bj ¨ ornson, J. Hoydis, and L. Sanguinetti, “Massiv e MIMO networks: S pectral, energy , and hardware efficiency , ” F oundations and T r ends in Signal Pr ocessing , v ol. 11, no. 3-4, pp. 154–655, 2017 . [Online]. A v ailable: http://dx.doi.org/10 .1561/20000 00093 [8] J. H. Sorensen, E . de Carvalho, and P . Popovski, “Ma ssiv e MIMO for cro wd scenarios: A solution based on random access, ” in Globecom W orkshops (GC Wkshps), 2014 , Dec 2014, pp. 352–357. [9] E . de Carv alho, Bj ¨ orn son, E. G. Larsson, and P . Popo vski, “Random access for Massi ve MIMO systems with intra-cell pilot contamination, ” in 2016 IE E E International Confer ence on Acoustics, Speec h and Signal Pr ocessing (ICASSP) , March 2016, pp. 3361 –3365. [10] E. Bj ¨ ornson, E. de Carvalho, J. H. Srensen, E . G. Larsson, and P . Popovski, “ A random access protocol for pilot allocation in crowde d massi ve MIMO sy stems, ” IEEE T ransactions on W ireless Communications , vol. 16, no. 4 , pp. 2 220–223 4, April 2017. [11] S. Sesia, I. T oufik, and M. Baker , LTE - the UMTS l ong term evolution: from theory to pr actice . Chichester: Wiley , 2009. [Onli ne]. A vailable: http://opac.inria.fr/record=b11309 16 [12] Y . Zhou, Z. Zhang, and X. Zhou, “OFDMA initial ranging for IEEE 802.16e based on time-domain and frequenc y-domain approaches , ” in International Confer ence on Communication T echnology , Nov 2006, pp. 1–5. [13] H. A. Mahmoud, H. Arslan, and M. K. Ozdemir , “Initial ranging for W iMAX (802.16e) OFDMA, ” in MILCOM 2006 - 2006 IEE E Military Communications confer ence , Oct 2006, pp. 1–7. [14] L. Sanguinetti, M. Morelli, and L. Marchetti, “ A random access algorithm for L TE systems, ” T rans. Emer ging T elecommun. T echno l. , vol. 24, no. 1, 2013. [15] M. Morelli, L . Sanguinetti, and H. V . Poor, “ A robust ranging scheme for OFDMA-based networks, ” I E EE Tr ans. Commun. , vol. 57, no. 8, pp. 2441 – 2452, Aug 2009 . [16] L. S anguinetti, M. Morelli, and H. V . Poor, “ An ES PRIT-based appro ach for initial ranging in OFDMA systems, ” IEEE T rans. Commun. , vo l. 57, no. 11, pp. 3225 – 3229, Nov 2009. [17] L. S anguinetti and M. Morelli, “ An initial ranging scheme for the IEEE 802.16 OFDMA uplink, ” IEEE T rans. W ir eless Commun. , vol. 11, no. 9, pp. 320 4–3215, Sept. 2012. [18] R. Miao, L. Gui, J. S un, and J. Xiong, “ A ranging method for OFD MA uplink system, ” IEEE T rans. C onsumer El ectr onics , vol. 56, no. 3, pp. 1223–1228, Aug 2010. [19] H. Han, X. Guo, and Y . Li, “ A high throughput pilot allocation for m2m communication in crowde d massiv e mimo systems, ” IEEE T ransaction s on V ehicular T echn ology , vol. 66, pp. 9572–957 6, 2017. [20] H. Han, Y . Li, and X. Guo, “ A graph-based random access protocol for crowd ed massiv e mimo systems, ” IEEE Tr ansactions on W ireless Communications , vol. 16, pp. 7348–73 61, 2017. [21] M. Morelli, C. C. J. Kuo, and M. O. Pun, “Synchronization techniques for orthogon al frequenc y division multiple access (OFDMA): A tutorial revie w , ” P roceed ings of the IEEE , vol. 95, no. 7, pp. 1394–1427, July 2007. [22] X. F u, Y . L i, and H. Minn, “ A new ranging method for OF DMA systems, ” IEEE T rans. W ireless Commun. , v ol. 6, no. 2, pp. 659 – 669, Feb 2007. 32 [23] M. Ruan, M. C. Reed, and Z. Shi, “Successi ve multiuser detection and interference cancelation for conten tion based OFDMA ranging channel [ transactions letters], ” IEEE T rans. W ir eless Commun. , vo l. 9, no. 2, pp. 481 – 487, February 2010. [24] R. R oy , A. A. Paulraj, and T . Kailath, “ESP RIT – direction-of-arri val estimation by subspace rotation methods, ” IEE E T rans. Acoustic, Speech Signal Pr oc. , vol. 37, no. 7, pp. 984 – 995, July 1989. [25] L. Sanguinetti, A. A. D’Amico, M. Morelli, and M. Debbah, “Random access in uplink Massiv e MIMO systems: How to exploit asynchronicity and excess antennas, ” in 2016 IEE E Global Communications C onferen ce (GLOBECOM) , Dec 2016, pp. 1–5. [26] M. Fallgren, B. Timus et al. , D1.1: Scenarios, requir ements and K PIs for 5G mobile and wir eless system . ICT -317669-METIS, 2013. [Online]. A vailable: https://www .metis2020.com/ [27] H. Q. Ngo, E. G. Larsson, and T . L. Marzetta, “ Aspects of fav orable propagation in massiv e MIMO, ” in Eur opean Signal Pr ocessing Confer ence (EUSIPCO) , Sept 2014, pp. 76–80. [28] H. Q. Ngo and E. G. Larsson, “No downlink pilots are needed in TDD Massi ve MIMO, ” IEEE Tr ans. W ireless Commun. , vol. 16, no. 5, pp. 2921–2935, May 2017. [29] P . Stoica and R. L . Moses, Introd uction to spectral analysis . Upper Saddle Ri ver , N.J. Prentice Hall, 1 997. [Online]. A vailable: http://opac.inria.fr/record=b10924 27 [30] M. W ax and T . Kail ath, “Detection of si gnals by information theoretic criteria, ” IEEE T rans. Acoustic, Speech Signal Pro c. , vol. ASSP - 33, pp. 387 – 392, April 1985. [31] Q. T . Zhang, K. M. W ong, P . C. Y ip, and J. P . R eilly , “Statistical analysis of the performance of information theoretic criteria in t he detection of the number of signals in array processing, ” IEEE T rans. Signa l Pr ocess. , v ol. 37, no. 10, Oct 1989. [32] A. P . Liav as and P . A. Reg alia, “On the behavior of information theoretic crit eria for model order selection, ” IEEE T rans. Signal Pro cess. , vol. 49, no. 8, Aug 2001. [33] E. Fishler, M. Grosmann, and H. Messer , “Detection of signals by information theoretic criteria: general asymptotic performance analysis, ” IEEE T rans. Signal Pr ocess. , vol. 50, no. 5, May 2002. [34] B. Nadler , “Nonparametric detection of signals by information theoretic criteria: Performance analysis and an improved estimator , ” IEEE T rans. Signal Pr ocess. , vol. 58, no. 5, May 2010. [35] L. Zhao, P . Krishnaiah, and Z . Bai, “On detection of the number of signals in presence of white noise, ” J ournal of Multivariate Analysis , vol. 20, no. 1, pp. 1 – 25, 1986. [Online]. A vailable: http://www .sciencedirect.com/science/article/pii/00472 59X869001 75 [36] R. O. Schmidt, “Multiple emitter l ocation and signal parameter estimation, ” IEEE Tr ansactions on Antennas and Pr opagation , v ol. 34, no. 3, pp. 276–280, 1986. [37] P . Stoica and A. Nehorai, “MUSIC, maximum likelihood, and Cramer-Rao bou nd, ” IEE E T ransactions on Acoustics, Speec h, and Signal P r ocessing , vol. 37, no. 5, pp. 720–74 1, May 1989 . [38] S. M. Kay , Fundamentals of Statistical Signal Pr ocessing: Estimation Theory . Prentice Hall, 1993 . [39] S. L. L oyka , “Channel capacity of MIMO architecture using the expon ential co rrelation matrix, ” IEEE Commun. Lett. , vol. 5, no. 9, pp. 369–371, Sept 2001.