1

Random Pilot and Data Access in Massive MIMO
for Machine-type Communications

arXiv:1709.01347v1 [cs.IT] 5 Sep 2017

Elisabeth de Carvalho, Senior Member, IEEE, Emil Björnson, Member, IEEE, Jesper H. Sørensen, Member, IEEE,
Erik G. Larsson, Fellow, IEEE, Petar Popovski, Fellow, IEEE

Abstract—A massive MIMO system, represented by a base
station with hundreds of antennas, is capable of spatially multiplexing many devices and thus naturally suited to serve dense
crowds of wireless devices in emerging applications, such as
machine-type communications. Crowd scenarios pose new challenges in the pilot-based acquisition of channel state information
and call for pilot access protocols that match the intermittent
pattern of device activity. A joint pilot assignment and data
transmission protocol based on random access is proposed in
this paper for the uplink of a massive MIMO system. The
protocol relies on the averaging across multiple transmission slots
of the pilot collision events that result from the random access
process. We derive new uplink sum rate expressions that take
pilot collisions, intermittent device activity, and interference into
account. Simplified bounds are obtained and used to optimize
the device activation probability and pilot length. A performance
analysis indicates how performance scales as a function of the
number of antennas and the transmission slot duration.

I. I NTRODUCTION
In a massive multiple-input multiple-output (MIMO) system [1], the base station (BS) is equipped with a very large
number of antennas that create a very large number of spatial
degrees of freedom (DoFs) under the asymptotic favorable
propagation conditions that most channels seem to offer [2]–
[4]. Those spatial DoFs can be used in different ways. If the
number of antennas at the BS is much larger than the number
of wireless devices [1], then the extra DoFs are exploited
to generate strong spatial beams to the devices, which are
in effect hardened communication channels with negligible
small-scale fading and a minor multi-user interference. When
the number of devices is large, the DoFs can be used for
aggressive spatial multiplexing of devices, which does bring
the largest cell spectral efficiency [5]. However, the devices’
channels are not asymptotically decorrelated in this operating
regime and the spectral efficiency per device is thus lowered.
The enhanced spatial resolution of massive MIMO is essential
E. de Carvalho, J. H. Sørensen, and P. Popovski are with the Department of Electronic Systems, Aalborg University, Aalborg, Denmark
(email: edc@es.aau.dk, jhs@es.aau.dk, petarp@es.aau.dk). E. Björnson and
E. G. Larsson are with the Department of Electrical Engineering (ISY),
Linköping University, Linköping, Sweden (email: emil.bjornson@liu.se,
erik.g.larsson@liu.se).
This work was performed partly in the framework of the Danish Council
for Independent Research (DFF133500273), the Danish National Advanced
Technology Foundation via the VIRTUOSO project, the Horizon 2020 project
FANTASTIC-5G (ICT-671660), the EU FP7 project MAMMOET (ICT619086), ELLIIT, and CENIIT. The authors would like to acknowledge the
contributions of the colleagues in FANTASTIC-5G and MAMMOET.
This work was presented in part at the IEEE International Conference on
Acoustics, Speech and Signal Processing (ICASSP), Shanghai, China, 2016.

in crowd scenarios, such as shopping mall or a stadium [6],
serving high density of devices that are closely spaced. The
scenario that sets the motivation for this article is the one
of massive Machine-Type Communication (mMTC), where
around 100 000 devices could be connected to a single BS.
Acquisition of channel state information (CSI) in massive
MIMO is critical, as a huge number of channel coefficients
have to be estimated per device. Time division duplexing
(TDD) offers a simplified CSI acquisition in a massive MIMO
system with a relatively small number of devices. Assuming
channel reciprocity, the channels are estimated at the BS
relying on uplink training based on orthogonal pilot sequences
and are then used for both uplink and downlink beamforming.
Hence, the length/number of the pilot sequences scales with
the number of devices and not the number of BS antennas.
In crowd and mMTC scenarios, the channel estimation faces
fundamental limits because of two specific features. First, the
number of pilot sequences is limited by the dimensionality of
the coherence interval of the channel, during which the channel
can be considered constant and flat. When the crowd of device
becomes very large, the orthogonal pilots are in severe shortage and the allocation policy of the pilot sequences becomes
a central question. Second, the intermittence of the data traffic
also plays a significant role. In mMTC, each device sends
data to the BS sporadically, at random time instants. Since
not all devices are active simultaneously, the limited set of
orthogonal pilot sequences can highly likely accommodate the
subset of devices that are active at a particular instant. Hence,
pilot allocation has rather to adapt and scale with the traffic
activity pattern and not equal the actual number of devices
present in the system. A natural choice is to decentralize pilot
access to the devices and make it random, which leads to pilot
collision, also known as pilot contamination.
While massive MIMO is a fairly mature research topic [2],
[5], [7]–[10], the existing results on uplink capacity analysis
in the literature [9] assumes a pre-defined pilot allocation
and full data buffers at the devices. Those assumptions are
not applicable to the case we study here. A joint pilot and
data protocol based on random uplink access in a massive
MIMO system has been recently considered in [11], [12],
with the use of coded access and successive interference
cancellation (SIC). In [13], [14], a set of dedicated pilots
is temporarily allocated to the active devices, which in turn
apply random access to the pilot sequences in the set. This
leads to collisions, which are resolved by giving priority to
devices with good channels, while the data transmission is
carried out without pilot collisions. Both schemes [11]–[14]

2

rely on channel hardening. In [15], timing variations between
the devices are used to devise a random access scheme for
initial access, synchronization, and channel estimation. Some
preliminary results on the effect of intermittent device activity
can be found in [16].
In the approach proposed in this paper, transmission is
organized into transmission frames consisting of multiple
transmission slots. When an mMTC device has a data codeword to transmit, it divides the codeword into multiple parts
and, within a transmission frame, transmits the codeword parts
in multiple transmission slots. The channels are estimated from
uplink pilots every time the device transmits. In each time
slot, each active device selects pseudo-randomly a pilot from
a predetermined pilot codebook and, during the rest of the slot,
it sends a part of the data codeword. Hence, a device performs
pilot-hopping over multiple slots and the hopping pattern can
be used to identify the device and appropriately merge and
decode the parts of its codeword at the BS. This approach
is suitable to harvest the large array gain for delay-tolerant
mMTC and low-power devices. It is important to note that our
access scheme relies on non-orthogonal pilot-hopping patterns,
each sequence consisting of orthogonal pilots. Another option
would be to build hopping patterns based on non-orthogonal
pilots and our scheme can be readily extended to that case.
When the number of transmission slots is sufficiently large,
the pilot collision events get averaged, so that we can define
a maximal achievable uplink sum rate. We determine a lower
bound on the sum rate that is tight thanks to channel hardening
and when the total number of terminals is large. This bound
is our main metric for performance assessment. Furthermore,
in order to optimize the system performance, the devices
are allowed to transmit with a certain activation probability
pa which is subject to optimization, along with the number
of training sequences. As the numerical evaluation of our
main performance bound can be cumbersome, optimization
is based on alternative bounds, tailored to account for rather
large variations in channel energy among the devices. We
provide a performance analysis for different asymptotic cases.
In particular, when the number of antennas M and the
duration of transmission
slot τu are of the same order, the
√
sum rate scales of M τu . Heuristic solutions giving robust
performance results indicate that one third of the transmission
slot should be devoted to training while the
√ average number
of active devices should be at the order of M τu . Compared
to [17], we provide performance analysis and solutions for the
case where the channel energy varies across the devices.
The paper is organized as follows. Section II describes the
random access model, the communication protocol, and the
channel model. Section IV presents the main performance
bound, Section V the lower bounds used for optimization
and Section VI the heuristic solutions. In Section VII, scaling
laws are presented as a function of the number of antennas
at the base station and the duration of a transmission slot.
Section VIII contains numerical evaluations and is followed
by the conclusions.

II. S YSTEM M ODEL
A. Random Access Model
We consider the uplink (UL) of a single-cell multi-device
massive MIMO system with random access from a large set of
intermittently active devices. The justification for a single-cell
model is found in Section II-D. The BS is equipped with M
antennas and serves a maximal number of K devices.
A total number of τp orthogonal sequences are available,
denoted as {Φ1 , Φ2 , . . . , Φτp }, where each sequence is τp
symbols long. The duration of a pilot sequence is smaller than
the channel coherence interval. Moreover, K  τp such that
the BS cannot allocate dedicated pilots to particular devices.
The system is assumed to be synchronized to maintain pilot
sequence orthogonality. It is important to note that an OFDM
system can tolerate timing mismatches provided that the cyclic
prefix is larger than the channel delay spread and the maximal
timing mismatch. Such conditions could be satisfied in small
enough cells in LTE as the normal cyclic prefix is designed for
a 750 m cell radius and the extended cyclic prefix corresponds
to a 2.5 km cell radius.
A UL transmission frame, see Figure 1, is divided into
transmission slots of duration τu , where τu is less or equal
to the channel coherence interval. Due to the random access,
the BS does not know a priori for a given slot which devices
transmit or which pilots are activated. Instead, each device
is associated to a unique, predefined pseudo-random pilothopping pattern such that an active device selects the pilot
for a given transmission slot according to this pattern. The BS
knows in advance the pilot-hopping patterns of all potential
transmitters, such that it can buffer the information from
different slots, run a correlation decoder across the slots and
find out which pilot-hopping patterns have been activated.
When a device has data to send, it encodes it into codewords,
so that one codeword is sent within a given transmission frame.
Each codeword is divided into multiple parts and sent over
multiple time slots within the transmission frame. The number
of codeword parts is equal to the number of pilot sequences
in the pilot-hopping identifying patterns. Transmission from
the K devices is sporadic and a device is active within a
given transmission frame with activation probability pa , independently from the other devices. An active device transmits
in all transmission slots of a transmission frame. Within a
transmission frame, a device is active only when it has data
to transmit.
In each active UL transmission slot, the pilot phase is
followed by a data phase, i.e., transmission of a part of a
codeword. Collisions can thus happen in the pilot domain, i.e.,
among contending devices that send the same pilot sequence to
the BS. Pilot collision is in fact similar to pilot contamination
in conventional massive MIMO and induces interference in
the data domain from the set of contaminators.
For an asymptotically large number of time slots, the whole
codeword is affected by an asymptotically large number of
channel fading realizations, pilot collision-induced interference and other interference events. Relying on the ergodicity
of such a process, the system performance can be characterized
as in Section IV. Both the activation probability pa and the

3

'HYLFH
'HYLFH
'HYLFH
'HYLFH

઴૚ ‫ܹܥ‬஽భ ሺͳሻ

઴ଶ ‫ܹܥ‬஽మ ሺͳሻ
઴૚ ‫ܹܥ‬஽య ሺͳሻ

઴ଶ ‫ܹܥ‬஽ర ሺͳሻ

7UDQVPLVVLRQ
VORW

઴૚ ‫ܹܥ‬஽భ ሺʹሻ ઴ଶ ‫ܹܥ‬஽భ ሺ͵ሻ

‫ڮ‬

઴ଶ ‫ܹܥ‬஽ర ሺʹሻ ઴ଶ ‫ܹܥ‬஽ర ሺ͵ሻ

‫ڮ‬

઴ଶ ‫ܹܥ‬஽మ ሺʹሻ ઴૚ ‫ܹܥ‬஽మ ሺ͵ሻ
઴ଶ ‫ܹܥ‬஽య ሺʹሻ ઴૚ ‫ܹܥ‬஽య ሺ͵ሻ
࣎࢖

࢛࣎

‫ڮ‬

઴૚ ‫ܹܥ‬஽భ ሺ‫ܮ‬ሻ

઴ଶ ‫ܹܥ‬஽మ ሺ‫ܮ‬ሻ

઴૚ ‫ܹܥ‬஽య ሺ‫ܮ‬ሻ
઴૚ ‫ܹܥ‬஽ర ሺ‫ܮ‬ሻ

7UDQVPLVVLRQ)UDPH

Fig. 1: Illustration of the transmission frame. In this example,
four devices {D1 , D2 , D3 , D4 } and two mutually orthogonal
pilot sequences {Φ1 , Φ2 } are considered. Transmission of a
codeword is done over multiple channel fades, which enables
averaging over noise, channel fades, and pilot collision events.

number τp of pilots per slot are optimized as elaborated in
Section V.
We model the process of constructing of a pseudo-random
hopping sequence by having each active device in each slot
select randomly one of the τp sequences. Hence, the probability of having Ka active devices within a total of K devices
that become independently active with probability pa is


K
Ka −K
a
p(Ka ) =
pK
.
(1)
a (1 − pa )
Ka
The average number of active devices is pa K and its variance
is pa K(1 − pa ). Given a number of Ka active devices and
considering one given device 0 among the Ka active devices,
the probability of having c colliders to device 0 can be
expressed as the probability that c devices select the same
sequence as device 0 in a population of Ka − 1, i.e.,
Ka −c−1

  c 
1
1
Ka −1
1−
. (2)
p(c|Ka ) =
c
τp
τp

B. Channel Model
A block fading model is adopted where a channel realization
is constant across a transmission slot duration and changes
independently from slot to slot. A transmission slot is defined
as a time-frequency interval that is smaller than the coherence
time and coherence bandwidth of the channel. In our model,
the received signal is divided by the standard deviation of the
receive additive noise, so that the normalized M × 1 receive
noise vector n at the BS is modelled as n ∼ CN (0, IM ). The
channel response between the BS and device j is described
by an M × 1 channel vector gj . The channel realizations
are modeled as circularly symmetric complex Gaussian distributed: gj ∼ CN (0, βj IM ). The variance βj reflects the path
loss, shadowing, transmit power (relative to the received noise
power). One given parameter βj is assumed invariant during
the whole transmission frame.
The devices transmit with different powers: device j transmits with power pj during both the training and data phase. We
recall that the parameter βj incorporates the transmit power
pj . Three models are considered for the distribution of βj .
In the first two models, power control (performed during the
initial access phase [14], [18]) is assumed, meaning that the

transmit power from each device is adjusted so that all the
βj ’s are equal.
Model 1: βj = δ 1 (1 + v), where v is uniformly distributed
between −α and α, with 0 ≤ α ≤ 1. It is assumed that
an imperfect power control is performed so that there
is an error around a nominal value δ 1 . This model was
chosen as large variations of α do not change significantly
the optimization in Section VII, while the following two
models do.
Model 2: βj = δ 2 10v/10 or 10 log10 (β j ) = 10 log10 (δ 2 ) + v,
where v is modelled as a Gaussian random variable with
variance σv2 . In this model, the parameter βj follows a
log-normal distribution, i.e. the conventional distribution
for shadowing. Power control is performed so that the
mean of the logarithm of βj is the same for all devices.
Model 3: βj = δ 3 (d/d0 )−αp , d = d0 (1 + v), where v is
uniformly distributed between −α and α, with 0 ≤ α ≤
1. Parameter d0 is set to 500 meters and αp = 3.76.
This model accounts for the path loss when devices are
distributed uniformly around a nominal distance from the
cell center.
In the simulations, δ 1 = δ 2 = δ 3 and δ j is set to 10 dB. As will
be made clear in Section V, the final results do not significantly
depend on the values selected for δ 1 , δ 2 , δ 3 , provided we are
in interference-limited scenarios. Furthermore, we introduce
the notations E[βj ] = β, E[βj2 ] = β 2 , E[βj4 ] = β 4 . Note that
β is equal to δ 1 for Model 1 but not to δ 2 or δ 3 for Model 2
and Model 3.
We use (·)∗ , (·)T , (·)H , E[·] to denote complex conjugation,
transpose, Hermitian transpose, and the expected value of a
random variable, respectively.
C. Training Phase
During the training phase, the active devices send the
normalized pilot training sequences Φk , ΦH
k Φk = 1, of length
0
τp . These sequences are orthogonal: ΦH
Φ
k0 = 0, for k 6= k .
k
The M × τp matrix of the received pilot signals at the BS is
Yp =

√

τp

Ka
X

gk ΦTk + Np ,

(3)

k=0

where Ka is the number of active devices in a given slot and
Np is an M × τp matrix whose columns are the noise vectors
during the training phase.
One given device, device 0, is affected by a set of contaminators denoted as C0 , which use the same pilot sequence
that we generically denote as Φ[0] , i.e. Φ[0] = Φj , j ∈
{1, · · · , τp }.
D. Multi-Cell Operation
The main line of work on massive MIMO, starting from [1],
has assumed that all devices in a given cell use orthogonal
pilots and analyzed the system performance under inter-cell
pilot collisions due to pilot reuse across cells. Our main
assumption is that intra-cell pilot collision happens because
the crowd of devices calls for random pilot access. This
assumption is by itself a departure from the conventional

4

multi-cell interference models. Note that [19] considers intercell pilot contamination in the sense that it provides a pilot
reuse plan within a cell so that users with very different
spatial characteristics use the same pilot, hence suppressing
pilot contamination. In contrast, our system is uncoordinated
with no specific assumptions of the channel properties, so the
intra-cell pilot contamination can be strong. Concerning the
inter-cell pilot collision, there are two possibilities.
• One possibility is to assign orthogonal pilot sets to each
cell where the sets are also orthogonal across neighboring
cells, so that one can neglect the inter-cell pilot collision.
This protocol has the disadvantage of decreasing the
number of available orthogonal pilot sequences per cell
for a given pilot sequence length, as the orthogonal
sequences are distributed across cells. The advantage is
that pilot contamination can be better controlled when it is
confined to a cell as opposed to multiple cells. This model
is equivalent to the single-cell model in Section II-A as
far as inter-cell pilot contamination is concerned. In data
transmission mode, inter-cell interference is present.
• The other possibility is to have the whole set of orthogonal sequences available in each cell, resulting in
inter-cell pilot collision. Assuming the same activation
probability in all the cells and a uniform spatial coverage,
the difference with a single-cell model is:
– The distribution of channel gains {βj }. Model 3 can
still be applied to multiple cells as it is based on
distance to the BS of interest. However, it is not the
case for Model 1 and Model 2 as the underlying
parameters take as reference the cell the device belongs
to. The distribution would be different for the devices
in the cell and the devices in the interfering cells.
– The sum rate involves only the devices in the cell of
interest.
In this paper, we focus on a single-cell model. As the random
access protocol presented in this paper is new, an initial analysis based on this simplified model is necessary to understand
the main mechanism behind the protocol. The study of a multicell system as presented above is left for future work.

is available: Yp Φ∗j , j = 1, . . . , τp . Pilot j is detected if
kYp Φ∗j k, j = 1, . . . , τp is larger than a pre-determined
threshold and nondetected otherwise. The pilot detection
outcomes are buffered in order to be utilized for device
activity detection, see Step 2.
2) Channel estimation: For each pilot sequence detected,
a scaled version of the MMSE channel estimate is determined. This estimate is the same for all members of a
contamination set.
3) Multi-antenna processing: For each pilot sequence detected, Maximum Ratio Combining (MRC) is applied
to the data symbols in the slot and the MRC output is
buffered along with its associated pilot index.
4) End of transmission: The process is repeated until all
active devices have stopped transmitting their codeword.
End of transmission is jointly detected with the active
device identification.
Step 2: Active device identification. From the pilot-hopping
patterns detected across the multiple slots at Step 1, the
transmitting devices are detected.
Step 3: Decoding. Based on the identified pilot patterns, the
BS identifies which MRC outputs to combine to decode the
data of each transmitting device. For each identified device 0,
data decoding assumes the knowledge of β0 .
B. MMSE Channel Estimation
The MMSE estimation of g[0] proceeds first by correlating
Yp with the pilot sequence Φ0 to obtain
√ X
gj + Np Φ∗[0] .
(4)
yp = Yp Φ∗[0] = τp
j∈C0

The MMSE estimate of g0 is then [20]
σh2 y
−1
ĝ0 = E[yp ypH ]
E[h0 ypH ] yp = 2 0 p yp ,
σyp yp

where E[yp ypH ] = σy2 p yp IM and E[h0 ypH ] = σh2 0 yp IM , with
σy2 p yp

III. R ANDOM ACCESS P ROTOCOL
The pilot and data random access protocol at the BS is
outlined in Section III-A. It relies on the channel estimation,
device identification and decoding procedure presented in
Sections III-B, III-C and III-D.

σh2 0 yp

The following procedure is considered for the transmission
of one codeword per device within a transmission frame. The
notations are defined in Section III-B.
Step 1: Multiple slot processing.
For one transmission slot, the following processing steps are
performed:
1) Pilot sequence detection: The BS detects which pilot
sequences are in use. This is done by correlating the
received signal Yp in (3) with each sequence Φj that

X

= τp β0 + τp
=

√

βj + 1,

(6)

j∈C0

τ p β0 .

(7)

The following results will be used in the computation of the
performance bound in Section IV-B. For j ∈ C0
ĝj =

A. Protocol Description at Receiver

(5)

βj
ĝ0 ,
β0

E[ĝ0 ĝ0H ] = σĝ20 IM ,

σĝ20 =

σh4 0 yp
σy2 p yp

2
E[ε0 εH
0 ] = σĝ0 IM ,
P
σh4 0 yp
τp j∈C0 βj + 1
2
σĝ0 = β0 − 2
=
.
σyp yp
σy2 0 y0

where ε0 = ĝ0 − g0 .
We define the quantity ĝ00 as:
ĝ0 =

β0 ĝ00 ,

ĝ00

=

τp

√
P

τp

j∈{0,C0 }

βj + 1

(8)

(9)

(10)

5

We use ĝ00 as a basis for the MRC operation at each transmission slot (see Section
P III-A).
The sum power j∈{0,C0 } βj can be estimated accurately
as detailed in [14]: indeed, using channel hardening, as the
T
H
∗
number of antenna grows large, the
P term Φj Yp Yp Φj is
2
asymptotically equal to σyp yp = τp j∈{0,C0 } βj + 1. Then,
for each pilot sequence detected, the corresponding channel
estimate ĝ00 is computed using (10).
C. Active Device Identification
Each pseudorandom pilot-hopping pattern serves as a
unique signature/identification for the device that applies that
pattern. The device patterns are known a priori to the BS, such
that the BS combines the pilot sequence detection outcomes
from the slots that follow the pattern in order to detect the
transmitting devices. The pseudorandom pilot-hopping pattern
detection follows the same principle as pseudorandom CDMA
sequence detection. We assume that the number pilot hopping
sequences is sufficiently large so that every user in the cell
has a unique identifying sequence.
D. Decoding Procedure
The MMSE channel estimate of device 0 requires knowledge of parameter β0 which cannot be assigned to device
0 until it is identified. To overcome this limitation, data
decoding is performed as follows. At first, the BS obtains,
for each transmission slot, the output of the MRC receiver
based on the channel estimates ĝ00 defined in (10) which
are estimated at each transmission slot. A pilot sequence
is linked to each of the MRC output at each transmission
slot. MRC is attractive in massive MIMO due to its low
computational complexity and asymptotic rejection of noncoherent multi-device interference and noise [7]. Next, the BS
identifies a transmitting device through its identifying pseudorandom sequence from which it can associate the relevant
MRC outputs over the multiple transmission slots. We assume
that the BS knows the associated β0 parameter (based
on the
P
initial access phase [14], [18]) and can estimate j∈{0,C0 } βj
based on equation (4) (using channel hardening as in [14]).
The BS performs soft combining of the MRC outputs scaled
by β0 .
IV. U PLINK SUM RATE
A. Maximal Achievable Sum Rate and Main Assumptions
Within a transmission frame, the codeword of a given device
experiences all possible contamination events from the Ka
active devices, provided that the number of transmission slots
is sufficiently large. Likewise, for an asymptotically large
number of transmission slots, a device experiences an asymptotically large number of additive noise and fading realizations.
Under those asymptotic conditions, a maximal achievable rate
per device can be defined within each transmission frame.
Achieving this rate assumes the following features:
1) The number of active devices within a transmission frame
can be estimated at the BS.

2) The average channel energy βj per device j can be
estimated at the BS and at the device.
3) Under conditions 1) and 2), the BS can compute a
maximal achievable rate per device. The BS broadcasts
both the channel energy βj and the associated rate for
each active device. As the device itself knows its channel
energy, it can associate the corresponding rate.
The different performance metrics presented in this paper
correspond to the maximal achievable sum rate averaged over
the activation probability of the devices.
We present two types of bounds on the average maximal
achievable sum rate. First, the main bound R1 in (14) tightly
approximates the maximal achievable sum rate and is used
for performance assessment. It necessitates Monte-Carlo simulations of the channel large scale fading realizations (i.e.,
β). Its numerical evaluation can be heavy as it involves the
computation of expected values over binomial distributions
with large number of devices K. A secondary bound, R2 , does
not require Monte-Carlo simulation but its tightness depends
on the distribution of the parameters {βj }. Other bounds R3 in
(20) and Ra in (21) are developed and serve as optimization
tools. They are relatively loose, but analytically simple and
follow the variations of the ergodic sum rate well. Those
bounds are the basis for the heuristic solutions presented in
Section VI.
B. Main Performance Bounds
Here we present the main bound R1 that is used for performance assessment. A secondary bound R2 is also provided
and mainly helps in establishing bound R3 in Section V.
Theorem 1 (Main bound R1 ). An approximation of the
average maximal achievable sum rate (in bits per symbol) is
R1 =

K
X

p(Ka ) Ka

Ka =1

KX
a −1

p(c|Ka ) Eβ [R1 (C0 , Ka , {β})] (14)

c=0

where

τu − τp
log2 1 + SINR1 (C0 , Ka , {β}) (15)
τu
is a lower bound on the maximal achievable rate of device 0
conditioned on a collider set with indices C0 and Ka active devices. The expectation is taken with respect to βj , j ∈ {0, C0 },
denoted as {β}. SINR1 (C0 , Ka , {β}) is given by (11) at the
top of the next page.
Recall that p(c|Ka ) is the probability of having c colliders
to device 0 when Ka devices are active and is given in
Equation (2). p(Ka ) is the probability of having Ka active
devices out of K and is given in Equation (1).

R1 (C0 , Ka , {β}) =

Proof. See Appendix A.
The computation of the sum rate involves averaging w.r.t.
the specific realization of the β’s during the whole transmission frame and computing quantities related to all possible
sets of different size among the K devices (See Appendix A).
As it is prohibitively complex, we instead compute an approximation (14) where the expected value w.r.t. the discrete
probability functions is replaced by an expected value w.r.t.

6

SINR1 (C0 , Ka , {β}) =

τp (M −1)

SINR2 (c, Ka , β0 ) =
SINR3 (β0 ) =

β 2 (M−1)(p

P

j∈C0

P
βj2 + i∈{0,C

0

τp (M −1)β02
P
P
β
(1+τ
p
} i
j∈C βj )+(1+
i∈{0,C
/
i

0}

βi )(1+τp

P

i∈{0,C0 }

τp (M −1)β02
2
2
τp (M−1)β c+β0 (1+τp cβ)−cβ τp +(1+(Ka −1)β)(1+β0 τp +τp cβ)

a K−1)+β0 (1+β(pa K−1))−β

2

τp (M −1)β02

βi )

(11)

(12)

2

(pa K−1)+(1+(pa K−1)β)(1+β0 τp )+(pa K−1)β+β (p2a K(K−1)−(pa K−1))

continuous probability functions. As K becomes large, the
approximation becomes very good. The results show that R1
approaches a limit when K → ∞.
Furthermore, after a certain threshold on K, R1 does not
increase significantly. This can be understood as the tail
of the binomial distribution becomes very small for large
values
K and brings the associated expected value of term
PKofa −1
Ka c=0
p(c|Ka ) Eβ [R1 ({β}, C0 |Ka )] in (14) to zero.
This result also means that after this threshold, the individual
rate per device diminishes.
Under the assumptions described in Section IV-A, device 0
computes its transmission rate as

V. B OUNDS FOR O PTIMIZATION AND P ERFORMANCE
A NALYSIS

We now determine the values of the activation probability
pa and the pilot sequence length τp that maximize the performance bound R1 . As the numerical evaluation of R1 is heavy
for large values of K, we use two lower bounds on R1 that are
the main bounds used for optimization as well as performance
analysis. Such an optimization assumes the knowledge of the
distribution of the β’s. As pa K is in fact the optimization
parameter, the knowledge of the total number of devices K is
not necessary for the optimization itself, but it is necessary to
decide the reliable rate for each device.
KX
K
We mostly treat the case where τp is not linked to an actual
a −1
X
1
channel,
p(Ka )Ka
p(c|Ka ) E{β}\β0 [R1 (C0 , Ka , {β})]
(16) meaning that we do not integrate in the optimization
R0 =
M
c=0
problem
that the constraint τp should be smaller than the
Ka =1
coherence
time or is limited by the transmit power of the
where E{β}\β0 designates the expectation w.r.t. {βk , k 6= 0}.
devices. The main results from the latter case are mentioned
in Appendix B. Furthermore, as mentioned in Section II-B, we
C. Secondary Performance Bounds
assume that a device always has data to transmit at activation
The next bound R2 is a secondary bound that is used as a time. To summarize, the optimization problem is:
basis for the bounds in Section V. Using Jensen’s inequality
max
Ro (τp , pa ).
(19)
on log2 (1 + 1/x), a lower bound is obtained by taking the
0≤τp ≤τu ,0≤pa ≤1
expected value of the denominator in (12) w.r.t. all sets
Ro denotes generically the different cost functions considered
of contaminators to device 0. The main advantage of R2
in Sections IV, V and VI.
compared with R1 is that the equivalent SINR, SINR2 (β0 ), in
In the next bound R3 , the expectation is taken in the
(12) can now be explicitly written as a function of the number
denominator of (12) w.r.t. the distribution of the number of
of contaminators to device 0 and number of active devices.
contaminators and the number of active devices. The bound
Theorem 2 (Secondary bound: R2 ). A lower bound on R1 R3 is relatively loose as compared to R1 and R2 , since it
averages over the number of colliders and active devices in
(in bits per symbols) is
the interference variances. However, R3 follows very well the
KX
K
a −1
X
variations
of R1 and provides very good optimization results
R2 =
p(Ka ) Ka
p(c|Ka )Eβ0 [R2 (c, Ka , β0 )] (17)
provided
that
the variance of the parameter β is not too large.
c=0
Ka =1
This aspect is highlighted in the numerical results section.
where
Theorem 3 (Bound for optimization: R3 ). A lower bound on
τu − τp
R2 (c, Ka , β0 ) =
log2 (1 + SINR2 (c, Ka , β0 )) (18) R (in bits per symbol) is
1
τu
τu − τp
and SINR2 (c, Ka , β0 ) is shown in (12) at top of the next page.
(20)
R3 =
pa K Eβ0 [log2 (1 + SINR3 (β0 ))]
τu
p(c|Ka ) and p(Ka ) are defined in (1) and (2).
In (17), we maintain the expected value w.r.t. β0 outside
of the log term, as its impact is not easily captured through
simplifying tools like Jensen’s inequality, if the β’s have
a large variance across the device population. Note that
SINR2 (c, Ka , β0 )) is a fraction depending on β0 and β02 and
its integral w.r.t. β0 can be written as a closed form expression,
so that Eβ0 [log2 (1 + SINR2 (c, Ka , β0 ))] can be written in
closed form for Model 1, i.e., a uniform distribution of β.

where SINR3 (β0 ) is given in (13) at the top of the next page.
Proof. In the denominator of SINR2 , we take the expected
value w.r.t. the probability mass of the binomial distribution
p(c) = p(Ka )p(c|Ka ). More specifically, we take first the
expected value of c conditioned on a number of active devices
Ka . It is the average number of contaminators to one given
device and is equal to (Ka −1)/τp . Then, we take the expected
value w.r.t. Ka , i.e., the average number of active devices out

(13)

7

of K devices which is equal to pa K. Hence, E[c] = (pa K −
1)/τp . We obtain a lower bound on the ergodic rate of a single
device. Next, we multiply by the average number of devices
pa K.

Asymptotic bound: Ra . If K  1, τu  1, M  1, τp  1
and pa K  1, a lower bound on the ergodic sum rate is
τu − τp
(21)
Ra = pa K
Eβ0 [log2 (1 + SINRa (β0 ))] ,
τu

60

50
50

Number of pilot sequences τp

From (13), a simpler bound can be deduced given asymptotic conditions on the parameters for which we keep the
equivalent expression of SINR3 (β0 ).

Asymptotic rate Ra (bits/symbols)

100

40
150

200

30

250
20
300
10
350

SINRa (β0 ) =

M τp β02
2

β 2 M pa K + β p2a K 2 + ββ0 pa Kτp

.

(22)

100

200

300

400

500

600

700

800

900

1000

0

Average number of users pa K

The expression of SINRa (β0 ) allows for an easy identification of the different terms and their effect. For that purpose,
we rewrite SINRa (β0 ) as
2

1
β 2 pa K
β p2a K 2
ββ0 pa K
=
+
+
.
2
2
SINRa (β0 )
τp β0
M τp β0
M β02

400

(23)

The first term in the sum comes from interference induced
by pilot collision. The second and third terms come from the
residual multi-device interference left by MRC.
We illustrate the variations of Ra in Figure 2 as a function
of τp and pa K, for M = τu = 100 and no variations across the
β’s (or perfect power control). For fixed pa K, the sum rate
first increases as τp increases since the probability of pilot
collision decreases. Hence, the first and second interference
τ −τ
terms decrease. Then, the rate penalty uτu p coming from
longer pilot sequences results in a sum rate decrease.
For a fixed number of pilot sequences τp , the sum rate
first increases through the pre-log term: the number of devices
increases faster than the SINR decreases from increased interference. When the decrease of the SINR becomes too large,
the sum rate decreases. This
decrease comes from the multi2
β p2 K 2
device interference term M τap β 2 . This last observation will be
0
the basis of heuristic solution (26) in Section VI.
VI. H EURISTIC S OLUTIONS
In this section, we present heuristic solutions to our optimization problem. To simplify the notations, we employ the
superscript ()o to denote the solutions of the optimization
problem regardless of the optimization cost function.
A. Small variations of {β}
Consider Model 1 and Model 3 in Section II-B with α  1.
A Taylor approximation of the log term in (21) is of the form
log2 (1 + SINRa (β0 )) = log2 (1 + SINRa (β)) + c1 α + O(α2 ),
where c1 is a scalar involving all the parameters except α.
Therefore, Eβ0 [log2 (1+SINRa (β0 )] = log2 (1+SINRa (β))+
O(α3 ). The same kind of derivation also holds for Model 3.
We note that the approximation error is of order α3 , which
explains why the approximation is quite robust to values of α
that are not strictly speaking “small”, as seen in the numerical

Fig. 2: Ra as a function of τp and pa K for M = τu = 400.
results. For Model 1, we found that the performance and
optimization are quite insensitive to large values of α. For
Model 3, a value of α = 0.25 results in optimized performance
that is close to the ones with perfect power control (α = 0),
while, for Model 2, this value is α = 0.1.
For small variance in β’s, the approximate cost function is
τu − τp
Ra ≈ pa K
(24)
log2 (1 + SINRa (β)).
τu
In Section V, we noted that the variations of Ra depend on
the term in p2a K 2 in the denominator of SINRa . The heuristic
solution presented next is obtained based on (24) where only
the term in p2a K 2 is kept.
Heuristic Solution 1. A heuristic solution to the optimization
of parameters τp and pa maximizes the modified cost function


τu − τp
M τp
[0]
Rh = pa K
(25)
log2 1 + 2 2 .
τu
pa K
The solution is
τpo

τu
=
3

poa K

√
τu M
= √
3s0

(26)

x
where s0 ≈ 3.92 is the solution of log(1 + x) = 2 1+x
.

Proof. We look for the expression of τp and pa K maximizing
the following rate function:
M τp
Rh = pa K(τu − τp ) log2 (1 + X), X = 2
. (27)
β β −2 p2a K 2
The
derivative of Rh are
( partial
∂Rh
∂X 1
= − log2 (1 + X) + (τu − τp ) ∂τ
,
∂τp
p 1+X
(28)
h
∂R
∂X 1
= log2 (1 + X) + pa ∂pa 1+X .
∂pa
∂X
∂X
Noting that pa ∂p
= −2X and ∂τ
= τ1p X, we obtain
a
p

τu
1 + 1+X
X log(1 + X) = τp ,
(29)
1+X
X log(1 + X) = 2.

From those equations, we obtain (26).
An advantage of (24) is that it does not depend on the
statistics of the β’s, such that they need not to be known.

8

B. Large variations of {β}
For large values of α, the performance of the solution
previously described degrades. We observe that the optimal
solution involves a larger average number of active devices
pa K while τu /3 remains a good approximation for τp . In the
following heuristic solutions, we set τpo = τu /3, while pa K is
[0]
optimized based on Ra or heuristic Rh . The superscript 1D
below indicates that the optimization is now one dimensional.
Asymptotic 1D. A heuristic solution to the optimization of
parameters τp and pa is
p
τu
poa K = bo τu M
(30)
τpo =
3
where bo maximizes the sum rate Ra where τp is set to τu /3:
!
√
2
M
τ
β
/3
u
0
R[1D]
= b Eβ0 log2 1 +
a
2√
bβ 2 M + b2 β M τu + bββ0 τu /3
(31)
bo depends on M , τu and the statistics of the β’s.
Heuristic Solution 2. A heuristic solution to the optimization
of parameters τp and pa is
p
τu
τpo =
poa K = bo τu M
(32)
3
where bo maximizes the heuristic cost function (26) where τp
is set to τu /3:



β02
[1D]
Rh = b Eβ0 log2 1 +
(33)
3βb2
bo does not depend on M or τu but solely on the statistics of
the β’s.
VII. S UM R ATE S CALING L AWS
In this section, we give a summary of the scaling laws of
the optimized sum rate and optimized parameters based on
Ra and (19). As seen through numerical results, Ra -based
optimization gives results that are close to direct optimization
based on R1 , provided that the variance of the β’s is not
very large. Hence such an analysis makes sense under such
conditions. The details can be found in Appendix B.
Scaling Laws. Assuming τu  1, M  1, the scaling laws
of the optimized sum rate Ra and optimized parameters τpo
and poa K are summarized as follows. Parameter τp is not
constrained by the channel coherence time.
Case 1: M  τu :
 r 
τu
τu
o
τp ≈
+ O τu
(34)
2
M
s
β4 1 p
poa K ≈
M τu + O(τu )
(35)
2
β β2 2
 r 
τu
τu
Roa ≈
+ O τu
(36)
4 log(2)
M
r
τ 
τu
u
SINRa ≈
+O
(37)
M
M

Case 2: M  τu :
2/3
M
τu1/3 + o(τu ),
=
2
s

5/6
β4
M
o
pa K =
+ O(M )
2
β β 2 2τu
 2/3 !
M
Roa = M + O M
τu
 1/6
 1/3 !
M
M
+O
SINRa = 21/3
τu
τu
τpo



(38)
(39)
(40)
(41)

Case 3: M ∼ τu :
τpo = aτu ,

p
a = O(1), b = O(1) (42)
poa K = b M τu ,
p
o
Ra = O( M τu ), SINRa = O(1) (43)

In Case 1, the average number of active devices and the
sum rate is limited by τu . In Case 2, it is is limited by M .
In Case 3, the optimal number of pilot sequences and average
number of active devices become comparable. In Case 1 and
Case 2, the rate of each device becomes asymptotically small
but the average number of active devices that the system can
accommodate grows faster. In Case 3, the rate of each device
becomes constant while the average number of active devices
increases. The quality of service requirement should dictate
the choice of M and τu .
In Case 1 and Case 2, we notice that the optimization results
depend on the distribution of the parameters β through the
4
factor β2 2 . In Model 1, this factor has a small value range
β β
even for a wide range of values of α. Hence, for this model, the
optimization will not significantly depend on the distribution
of the β’s. On the other hand, for Model 2 and Model 3,
4
the factor β2 2 varies significantly with varying values of σv2
β β
and α and so do the optimized values. For Case 3, we notice
the same kind of behavior via simulations. Note that all those
points are detailed in the simulation section.
VIII. N UMERICAL R ESULTS
In this section we illustrate the performance of the random
pilot and data access protocol at optimized value of τp and
pa K. Due to the computational limitations in evaluating R1 ,
the value of K is not as large as could be expected in an
mMTC scenario. However, the value of K is chosen so that R1
is close to its limit when K grows to infinity (see Section IV)
so that our results are still representative. Furthermore, τu and
M are chosen to be of the same order.
The curves show three different metrics as a function of
the transmission slot duration τu and for a fixed number of
antennas M . The metrics are evaluated using the different
optimization methods described in the previous sections. The
metrics are: (1) Rate R1 in bits per symbols; (2) The optimal
number of pilot sequences; (3) The optimal activation probability through the quantity poa K, i.e., the average number of
active devices. For easy reference, the different optimization
methods are summarized here:
1) R1 -opt: direct optimization based on the main bound R1 .

9

Average sum rate - Model 1

50

40

K=800, α=0
K=400, α=0.25
K=400, α=0.5
K=600, α=0.25
K=600, α=0.5
K=800, α=0.25
K=800, α=0.5
K=1000, α=0.25
K=1000, α=0.5

30

20

10

100

200

300

400

500

600

700

50

40

30

20
K=800, α=0
K=400, α=0.5
K=800, α=0.5
K=1000, α=0.5

10

0
0

100

200

300

400

500

600

700

800

Transmission slot duration τ u

0
0

Optimal sum rate - Model 2

60

Optimal sum rate (bits/symbol)

60

800

Transmission slot duration τ u
Fig. 3: Model 1. Optimal sum rate R1 , for M = 100, K =
{400, 600, 800, 1000}, α = {0, 0.25, 0.5}.

R3 -opt: optimization based on bound R3 .
Ra -opt: optimization based on asymptotic bound Ra .
[1D]
Ra -opt: τpo is fixed to τu /3, poa K determined from Ra .
√
√
[0]
Rh -opt: τpo = τu /3, poa K = τu M / 3s0 in (26).
[1D]
Rh -opt: τpo = τu /3, poa K is based on the heuristic rate
[1D]
function Rh in (33).
In all the figures, except Figure 8, the number of antennas is
M = 100. We show performance for perfect power control
and different values of the variance of the β’s. For Model 1,
we show only one curve as performance is quite insensitive
to the variance of the β’s. Model 2 and Model 3 exhibit the
same kind of behavior and we only show performance for one
of the models to illustrate one particular performance aspect.
2)
3)
4)
5)
6)

A. Optimization of Main Bound R1
We first show performance results with direct optimization
of R1 . Our main purpose is to show the variations of the
optimized points w.r.t. the distribution of the β’s and w.r.t. to
total number of devices K. Those points are highlighted in
Section IV-B and VII. The heavy computations involved in
the numerical evaluations have imposed a limit on the number
of Monte-Carlo runs on parameters {βj } in (14). Hence,
some curves fluctuate but still indicate significant performance
results.
In Figure 3, for Model 1, the optimal sum rate is shown
for the case of perfect power control (α = 0) and α = 0.25,
α = 0.5 for different values of the total number of devices
K. We see that performance is insensitive to the values of
parameters α and K. The optimized values of τp and pa K
are also approximately the same.
In Figure 4, for Model 2, the optimal sum rate is shown for
σv2 = 0 and σv2 = 0.5, i.e., a large variation of the variance,
for different values of the total number of devices K. As

Fig. 4: Model 2. Optimal sum rate R1 , for M = 100, K =
{400, 800, 1000}, σn2 = {0, 0.5}.
mentioned in Section IV, there is K above which the rate
does not change significantly. In numerical evaluations that
are not shown in this paper, we also observed that the larger
the variance of the β’s, the larger the number of devices to
reach this steady state.
In Figure 5, for Model 2, the sum rate curves, optimized τp
and pa K are shown for different values of σv2 and a large K =
1000. Performance degrades as σv2 increases. However, the
degradation is small, even with large σv2 . The optimal number
of pilot sequences appears quite insensitive to the variance of
β, while the optimal average number of devices increases with
the variance. The same was observed for Model 3.
B. Optimization Methods
We now show the performance metrics using the six different methods previously described. In Figures 6, the sum
rate and optimized τp is shown for Model 3, α = 0.25. In
Figures 7, for the same model, α = 0.5, the sum rate and
optimized pa K is shown. The total number of devices K is
equal to 800. For all the parameter settings that we tested, the
performance using R3 or Ra are almost the same, so that we
refer only to Ra from now on.
For the medium value α = 0.25, Ra -opt follows closely
R1 -opt. There is a slight degradation for the heuristic solutions
for larger values of τu while heuristic solutions perform
slightly better for smaller values of τu . This behavior is
confirmed and magnified for α = 0.5. The optimal number
of pilot sequences tends to be overestimated compared to R1 opt when the transmission slot duration becomes large. The
gap in the estimation of pa K is quite large, but this does not
show significantly in the sum rate R1 . This can be explained
by the fact that R1 is quite flat around the optimal point.
At last, we show a more extreme case, with M = 400,
Model 2 and σv2 = 0.25, where in the selected range of τu , Ra
based methods perform significantly worse than the heuristic
methods. In such a case, power control is the preferred solution
as it leads to a robust and simple optimization. If power control

10

Optimal sum rate - Model 2

50

40

30

20
K=800, σ2v =0
K=1000, σ2v =0.25

10

Optimal sum rate - Model 3, α = 0.25

60

Optimal sum rate (bits/symbol)

Optimal sum rate (bits/symbol)

60

R1-opt
R 3 -opt

50

R a -opt
R [1D]
-opt
a

40

R [0]
-opt
h
R [1D]
-opt
h

30

20

10

K=1000, σ2v =0.5
0

0
0

100

200

300

400

500

600

700

800

0

100

Transmission slot duration τ u

200

300

(a)

Number of pilot sequences

Average number of users

140
120
100
80
60

K=800, σ2v =0
K=1000, σ2v =0.25
K=1000, σ2v =0.5

40
20

R1-opt
R 3 -opt

250

R a -opt
R [1D]
-opt
a

200

R [0]
-opt
h
R [1D]
-opt
h

150

100

50

0

300

400

500

600

700

800

Transmission slot duration τ u

Optimal number of pilot

sequences τ op

100

200

300

400

500

600

700

800

Transmission slot duration τ u

(b)

(b)

200

800

0

0
200

700

300

160

100

600

Optimal number of pilot sequences τ op - Model 3, α = 0.25

180

0

500

(a)

Average number of users poa K - Model 2

200

400

Transmission slot duration τ u

- Model 2

Fig. 6: Model 3. (a) Optimal sum rate R1 and (b) τp , for
M = 100, K = 800, α = 0.25.

Number of pilot sequences

180

is not possible, it is clear that resorting to direct optimization
of R1 is preferable, while heuristic solutions could be of
significant value.

160
140
120

IX. C ONCLUSION

100
80
60

K=800, σ2v =0

40

K=1000, σ2v =0.25

20

K=1000, σ2v =0.5

0
0

100

200

300

400

500

600

700

800

Transmission slot duration τ u

(c)

Fig. 5: Model 2. (a) Optimal sum rate R1 , (b) pa K and (c)
τp , for M = 100, K = {800, 1000}, σn2 = {0, 0.25, 0.5}.

Massive MIMO is recognized as a primary technology
in 5G for its ability to provide very high data rates. It is
also a major enabling tool for machine-type communications
where the sheer number of devices can be accommodated
relying on the high multiplexing gain of massive MIMO and
low-power devices can be served relying on its large array
gain. In this paper, we have proposed a joint pilot and data
transmission protocol based on random access that is adapted
to the intermittent activity pattern of the devices. This protocol
is organized in transmission slots and relies on the averaging
of the pilot collision events across the transmission slots.
It is suited for delay-tolerant and low-rate applications. We
have provided performance expressions as well as optimization

11

140

50

R1-opt
R 3 -opt

45

R a -opt

Optimal sum rate (bits/symbol)

Optimal sum rate (bits/symbol)

Optimal sum rate - Model 2, σ 2v = 0.25 - M=400

Optimal sum rate - Model 3, α = 0.5

55

R [1D]
-opt
a
R [0]
-opt
h
R [1D]
-opt
h

40
35
30
25
20
15

R 3 -opt
120

R a -opt
R [1D]
-opt
a

100

R [0]
-opt
h
R [1D]
-opt
h

80

60

40

20

10
0

5
0

100

200

300

400

500

600

700

0

800

100

200

300

400

500

600

700

800

Transmission slot duration τ u

Transmission slot duration τ u

Fig. 8: Model 2. Optimal average sum rate R1 for M = 400,
K = 800, σv2 = 0.25.

(a)
Average number of active users poa K - Model 3, α = 0.5

Average number of active users

300

where xk is the transmitted symbol from device k and n is the
additive noise vector. We isolate the contributions with index
set {0, C0 } and its complement and write gk as ĝk −εk , where
εk is the error vector in the MMSE estimation. The received
signal can be rewritten as
X
X
yd =
(ĝj − εj )xj +
g j xj + n
(45)

R1-opt
R 3 -opt

250

R a -opt
R [1D]
-opt
a

200

R [0]
-opt
h
R [1D]
-opt
h

150

j∈{0,C0 }

100

j ∈{0,C
/
0}

The BS applies MRC based on channel estimate ĝ0 in (5):
X
ĝ0H yd = ĝ0H ĝ0 x0 +
ĝ0H ĝj xj +
(46)

50

j∈C0

0
0

100

200

300

400

500

600

700

800

X

Transmission slot duration τ u

Fig. 7: Model 3. (a) Optimal sum rate R1 and (b) pa K, for
M = 100, K = 800, α = 0.5.

tools that are particularly important for a system where the
activity of the devices and the number of pilots have to obey
certain statistical rules.

We first compute the achievable sum rate for a fixed number
of active devices Ka . The M × 1 vectorial received signal in
the data phase is
(44)

(47)

j ∈{0,C
/
0}

R = E{hk } log2 (1 + SINR)

(48)

where SINR is equal to

j∈C0

A. Derivation of SINR1

gk xk + n

ĝ0H gj xj + ĝ0H n

where ĝ0H ĝ0 x0 is the signal of interest. The other additive
terms are treated as a Gaussian noise. Using the independence
of εj and ĝ0 in MMSE estimation, those terms are independent
from the signal of interest. From (8), the achievable rate is

X

The derivation of SINR1 ({β}) follows [5], [9]. The pre-log
term in (14) accounts for the activity probability of devices.
The main derivation steps are reported below.

i=1

X

|ĝ0H ĝ0 |2

A PPENDIX A
P ROOF OF T HEOREM 1

yd =

+

j∈{j,C0 }

(b)

Ka
X

ĝ0H εj xj

|ĝ0H ĝi |2 +

X

|ĝ0H εi |2 +

j∈{k,C0 }

X

|ĝ0H gi |2 + kĝ0 k2

.

(49)

j ∈{k,C
/
0}

The notation E{hk } indicates the expected
√ value w.r.t. small
scale fading terms hk , defined as gk = βk hk .
Following [9], using Jensen’s inequality on function
log2 (1 + 1/x) allows averaging over the channel fast fading
of the interferers. We obtain a lower bound as:
R1 ({β}) = log2 (1 + SINR1 ({β}))
where 1/SINR1 ({β}) is equal to
P
P
P
2
2
βj + 1
j∈{0,C0 }σεj
j ∈{0,C
/
j∈C0βj
0}
+
+
.
2
2
β0
(M − 1)σĝ0
(M − 1)σĝ20

(50)

(51)

12

Replacing the expression of σε2j and σĝ20 in (8) and (9), we
obtain (11).
The first term in (51) represents the interference induced by
pilot contamination, the second comes from channel estimation
errors and the third term comes is due to noise and residual
interference after MRC.

•

Finally, the sum rate is

PK

K
X

k=1

R01 (βk ) equal to

p(Ka |K)Ka

Ka =1

Ka−1
X

p(c|Ka0 ) (56)

c=0

K


1 X
EKak ,Ck R1 βk , Ck (c, Kak ({β}))|Kak ({β}) (57)
K
k=1

•

B. Derivation of the pre-log term
(50) gives a lower bound on the maximal achievable rate
of one given device 0 with a population of Ka active devices,
with exactly c contaminators and a given realization of the
β’s. This bound is tight compared to the maximal achievable
rate thanks to channel hardening. Next, we give the different
steps leading to (14). Note that when the notations become
too heavy and when there is no ambiguity, we give up the
dependency w.r.t. the parameters. Furthermore, the notations
in equations (14) and (15) have been simplified compared to
the notations adopted in this appendix.
•

We compute the achievable rate for device 0 averaged
over the possible sets of contaminators within the set of
Ka − 1 active devices excluding device 0. This last set
is denoted as Ka0 ({β}) and the set of contaminators is
denoted as C0 (c, Ka0 ({β})). The rate is


EC0 R1 β0 , C0 (c, Ka0 ({β}))|Ka0 ({β})

•

where EC0 represents the expected value w.r.t. all possible
sets of c contaminators within Ka0 .
Next, we compute the expected value of (52) accounting
for the random selection of the pilot sequences within the
set Ka0 and obtain
KX
a −1

•

K
1 X
EKak ,Ck [R1 ] ≈ E{β} [R1 ]
K

(58)

k=1

where E{β} is the expected value of all βk s with k 6= 0
w.r.t. their probability density function. Given that the
βk ’s are all independent, this approximation becomes
better as K grows large. Expression (58) gives us (14).
A PPENDIX B
S CALING L AWS
An asymptotic analysis is performed with the three cases:
• Case 1: M  τu , i.e., τu /M close to 0.
• Case 2: M  τu , i.e., M/τu close to 0.
• Case 3: M ∼ τu , i.e., M and τu are of the same order.
The analysis is directly based on the asymptotic expression
of the sum rate Ra in (21), instead of R3 , which is justified
by examining the asymptotic conditions τu  1 and M  1
and the ones considered below. Most of the results are based
on first or second order Taylor approximations. To ease the
presentation, we write the order of the approximations only
for the final optimized quantities. Furthermore, τp is not
constrained by the channel coherence time in Appendix B-A,
Appendix B-B, Appendix B-C. We treat this case in Appendix B-D.



A. Case 1: M  τu (or M  τp )
p(c|Ka0 )EC0 R1 β0 , C0 (c, Ka0 ({β}))|Ka0 ({β})
c=0
We consider different variation domains for pa K and τp in
(53)
order to study the optimization of the rate function Ra .
0
where p(c|Ka ) is the probability of having c contamina1) pa K  τp : In the expression Ra , SINRa can be
tors to device 0 within the set Ka0 .
approximated as:
Next, we compute the expected value of (53) accounting
M τp β02
for the activation probability, i.e., the probability to have
SINRa ≈
(59)
2
Ka − 1 active devices within a total number of K − 1
β 2 M Kpa + β p2a K 2
devices (as device 0 is excluded). We get
provided that β 2 M  ββ0 τp . Furthermore, if τp β02 
2
K
X
β pa K, then SINRa  1.
c
p(Ka − 1|K − 1)EKa0 [R1 ]
(54)
The two conditions on β0 are verified for Model 1 and
Ka =1
Model 3 (see Section II-B) provided α is such that β0 keeps
values of same order as β 2 /β and β. For Model 2, we assume
where EKa0 represents the expected value w.r.t. all possithat the largest value of β0 , denoted as βT , verifying the two
ble sets of Ka − 1 active devices excluding device 0.
conditions is sufficiently large so that the expected value of the
Now, accounting for the activation probability pa of delog-term in (21) restricted to domain [0, βT ] is approximately
vice 0 and noting that pa p(Ka −1|K −1) = KKa p(Ka |K),
equal to the expected value over the whole domain. This result
the achievable rate for device 0 averaged over the conis used below where the expected value w.r.t. β0 over the
tamination events and activation probability is
restricted domain is approximated as β. The same type of
conditions
on β0 are met in some of the derivations that will
K
Ka−1
X
1 X
0
0
follow,
and
can be treated in the same way so that we directly
R1 (β0 )=
p(Ka |K)Ka
p(c|Ka )EKa0 ,C0 [R1 ] (55)
K
consider
the
case of perfect power control to simplify.
c=0
Ka =1
Rc1=

•

(52)

We make the approximation

13

Next, we proceed to Taylor series approximations in order
to find the optimal values of parameter pa K first and then τp .
2
We assume β 2 M  β pa K, equivalent to M  pa K for the
models that we consider. A second order Taylor approximation
of the term in log gives:
"
#
2

1
1 τp β02
τu − τp
M τp β02
−
Ra log(2) ≈
Eβ0
2
τu
Kpa
β 2 M + β pa K 2 β 2
(60)
Taking the derivative w.r.t. the variable pa K, we get:
#
"
2
∂Ra
M τp β02 β
1 τp2 β04
+
≈0
≈ 0 ⇔ Eβ0 −
2
∂pa K
2 β 2 2 p2 K 2
[β 2 M + β pa K]2
a
(61)
2
We assume first β 2 M  β pa K. The other cases are treated
later. We obtain the expression of the optimal value of pa K:
s
1 β4 p
o
pa K =
M τp + O(τp ).
(62)
2 β2β2
Next, keeping only the first order terms in τp in (60), the
optimal value of τp is the one maximizing (τu − τp )τp , i.e.,
 r 
τu
τu
+ O τu
.
(63)
τp =
2
M
poa K

τpo

It can easily be verified that solutions
and are within
2
2
the definition intervals pa K  τp , β M  β pa K for M 
τu . The maximum value of the sum rate is
 r 
τu
τu
o
+ O τu
Ra =
.
(64)
4 log(2)
M
The last step is to verify that Ra does not get larger values
for the subcase where pa K  M or pa K ∼ M . For pa K 
M:


τu − τp
M τp
Ra ≈
pa K log2 1 + 2 2 .
(65)
τu
pa K
o
For
p a fixed
√ τp , Ra has its maximal value foro pa K =
M τp / 3s0 , s0 ≈ 3.92 (see Section VI). As pa K  M ,
then pa K > poa K and Ra decreases with pa K and is maximal
τ −τ
when pa K ∼ M , in which case Ra = uτu p τp /M . This
value is smaller than the maximal value found in the interval
pa K  M . For pa K ∼ M , the value of the sum rate is also
smaller.

2) pa K  τp : The sum rate is approximated as


τu − τp
τp
Ra ≈
pa K log2 1 +
.
τu
Kpa

(66)

For a fixed value of τp , Ra increases with pa K and takes its
maximal value when pa K ∼ τp , which case we consider next.
3) pa K = ατp , α = O(1):
Ra ≈

τu − τp
τp α log2 (1 + 1/α).
τu

(67)

The maximal value is smaller than the one obtained in the
case pa K  τp .

B. Case 2: M  τu
If M  τp , then the results for Case 1 tells that the optimal
value is of the form (62) and (63). This solution can be
verified to give a smaller sum rate than the solution that is
given when M  τp below. The case M ∼ τp is treated in
Appendix B-C. We now assume M  τp . Next we consider
different definition intervals for the parameters.
1) pa K  τp :
Ra ≈



M τp
τu − τp
pa K log2 1 + 2 2 .
τu
pa K

(68)

As already mentioned in Appendixp
B-A1, for
√ a fixed τp , Ra
has its maximal value for poa K = M τp / 3s0 , s0 ≈ 3.92.
As poa K  M , then pa K < poa K and Ra increases with
pa K and is maximal when pa K ∼ τp , in which case Ra =
τu −τp
τu M . This value is smaller than the maximal value found
in the interval pa K  τp .
2) pa K  τp : The SINR can be approximated as:
SINRa ≈

M τp β02
2

.

β p2a K 2 + β0 βτp pa K

(69)

With the same steps as in Case 1, pa K  τp , we obtain
s
1 β4 p
M τp + O(M ).
(70)
poa K =
2 β2β2
For the optimization w.r.t. τp , a first order Taylor approximation of Ra leads to
!
τu − τp
M τp β02
Ra log(2) ≈
(71)
2
τu
β pa K + β0 βτp
√ !
τu − τp
M
β0
≈
(72)
M Eβ0
− √
τu
τp
β
√ !
τu − τp
M
≈
M 1− √
(73)
τu
τp
Taking the derivative, we find that the optimal point τpo is the
positive root of 3rd order polynomial equation:
√
√
−2x3/2 + α M x + τu α M = 0.
(74)
Writing τp as τp = aτu , we get: −2 τMu a3/2 + a + 1 = 0.
Assuming a  1 and a  τMu a3/2 , we obtain a =
 2/3
M
. One can further prove that values of a that are
2τu
not close to 0 lead to invalid solutions. Finally, the solution
is:
2/3

M
τu + o(τu ).
(75)
τpo =
2τu
From (73), the sum rate can be written as:
 2/3 !
M
o
Ra log(2) = M + O M
.
τu

(76)

14

3) pa K = ατp , α = O(1):


τu − τp
M τp
Ra ≈
(77)
ατp log2 1 + 2 2
τu
α τp + ατp2


τu − τp
M
1
M2
≈
α
−
/ log(2). (78)
τu
α + 1 2 (α + 1)2 τp

0.9
0.8
0.7
0.6

α=0.25
0.5

One can prove that the associated maximal value is smaller
than the one obtained in the case pa K  τp .

0.4

α=0

0.3

C. Case 3: M and τu are of same order and M/τu = δ,
δ = O(1).
We show that the solution is of the following form:
p
τpo = aτu ; poa Ka = b M τu

α=0, α=0.25
0.1
10-2

(79)

a
b

0.2

10-1

100

101

102

103

delta = M/τu

Fig. 9: Optimization of parameters a and b in Equation (82)
where a and b are scalars of order 1 and are solution of the as a function of δ. Model 3 with α = 0 and α = 0.25.
maximization of the function:
"
!#
aβ02 δ
. (80)
Ra = (1 − a)b Eβ0 log2 1 +
√
√
[7] E. G. Larsson, F. Tufvesson, O. Edfors, and T. L. Marzetta, “Massive
2
bβ 2 δ δ + b2 β δ + abββ0 δ
MIMO for next generation wireless systems,” IEEE Commun. Mag.,
vol. 52, no. 2, pp. 186–195, 2014.
[8] H. Huh, G. Caire, H.C. Papadopoulos, and S.A. Ramprashad, “Achieving
“massive MIMO” spectral efficiency with a not-so-large number of
antennas,” IEEE Trans. Wireless Commun., vol. 11, no. 9, pp. 3226–
3239, 2012.
[9] H. Ngo, E. Larsson, and T. Marzetta, “Energy and spectral efficiency of
very large multiuser MIMO systems,” IEEE Trans. Commun., vol. 61,
no. 4, pp. 1436–1449, 2013.
[10] J. Hoydis, S. ten Brink, and M. Debbah, “Massive MIMO in the UL/DL
of cellular networks: How many antennas do we need?,” IEEE J. Sel.
Areas Commun., vol. 31, no. 2, pp. 160–171, 2013.
[11] J.H. Sørensen, E. de Carvalho, and P. Popovski, “Massive MIMO for
crowd scenarios: A solution based on random access,” in Globecom
Workshops (GC Wkshps), Dec 2014, pp. 352–357.
D. τp constrained by the channel coherence time: τp ≤ τpmax . [12] J. H. Sørensen, E. de Carvalho, C. Stefanović, and P. Popovski,
“Coded pilot access: A random access solution for massive MIMO
For cases M  τpmax and M  τpmax , using equations (62)
systems,” IEEE Trans. Wireless Commun., 2016, Submitted, Available:
http://arxiv.org/pdf/1605.05862.
and (70)
[13] E. Björnson, E. de Carvalho, E. G. Larsson, and P. Popovski, “Random
s
access protocol for massive MIMO: Strongest-user collision resolution
1 β4 q
max + O(M ).
(SUCR),” in Proc. IEEE ICC, 2016.
M
τ
(81)
poa K =
p
[14] E. Björnson, E. de Carvalho, J. H. Sørensen, E. G. Larsson, and
2 β2β2
P. Popovski, “A random access protocol for pilot allocation in crowded
p
massive mimo systems,” IEEE Transactions on Wireless CommunicaFor case M ∼ τpmax , poa Ka = b M τpmax where b optimizes
tions, vol. 16, no. 4, pp. 2220–2234, April 2017.
"
!#
[15] L. Sanguinetti, A. A. D’Amico, M. Morelli, and M. Debbah, “Random
2 0
β0 δ
access in uplink massive MIMO systems: How to exploit asynchronicity
Ra = b Eβ0 log2 1 +
. (82)
√
√
2
and excess antennas,” in IEEE Globecom, 2016.
2
0
2
bβ δ δ + b β δ 0 + bββ0 δ 0
[16] E. Björnson and E. G. Larsson, “Three practical aspects of massive
MIMO: Intermittent user activity, pilot synchronism, and asymmetric
with δ 0 = M/τpmax .
deployment,” in Globecom Workshops (GC Wkshps), 2015.
[17] E. de Carvalho, E. Björnson, E. G. Larsson, and P. Popovski, “Random
access for massive MIMO systems with intra-cell pilot contamination,”
R EFERENCES
in Proc. IEEE ICASSP, 2016.
[1] T. L. Marzetta, “Noncooperative cellular wireless with unlimited
[18] E. Dahlman, S. Parkvall, and J. Skold, 4G: LTE/LTE-Advanced for
numbers of base station antennas,” IEEE Trans. Wireless Commun.,
Mobile Broadband, Academic Press, 1st edition, 2011.
vol. 9, no. 11, pp. 3590–3600, 2010.
[19] A. Adhikary, J. Nam, J. Y. Ahn, and G. Caire, “Joint spatial division
[2] E. Björnson, E. G. Larsson, and T. L. Marzetta, “Massive MIMO: 10
and multiplexing - the large-scale array regime,” IEEE Transactions on
myths and one critical question,” IEEE Commun. Mag., vol. 54, no. 2,
Information Theory, vol. 59, no. 10, pp. 6441–6463, Oct 2013.
2016.
[20] S.M. Kay, Fundamentals of Statistical Signal Processing: Estimation
[3] A.O. Martinez, E. De Carvalho, and J. Nielsen, “Towards very large
Theory, Prentice Hall, 1993.
aperture massive mimo: A measurement based study,” in 2014 IEEE
Globecom Workshops (GC Wkshps), Dec 2014, pp. 281–286.
[4] H. Q. Ngo, E. G. Larsson, and T. L. Marzetta, “Aspects of favorable
propagation in massive MIMO,” in Proc. IEEE EUSIPCO, 2014.
[5] E. Björnson, E. G. Larsson, and M. Debbah, “Massive MIMO for
maximal spectral efficiency: How many users and pilots should be
allocated?,” IEEE Trans. Wireless Commun., vol. 15, no. 2, pp. 1293–
1308, 2016.
[6] 5G Mobile and Wireless Communications Technology, Cambridge
University Press, 2016.

The expression (80) is obtained from the expression of Ra in
(21). One key point is that is only depends on M/τu = δ.
We show that a and b are of order 1 provided that δ is
of order 1 through numerical evaluations. In Figure 9, the
optimized values of a and b are shown as a function of δ. For
values of δ around 1, the values of a and b are indeed of order
1. When δ tends to ∞, the value saturates to a = b = 1/2 as
predicted in Appendix B-A, while the values obtained when
δ tends to 0 correspond to Appendix B-B2.

