Approximations of the Aggregated Interference Statistics for Outage Analysis in Massive MTC †

This paper presents several analytic closed-form approximations of the aggregated interference statistics within the framework of uplink massive machine-type-communications (mMTC), taking into account the random activity of the sensors. Given its discrete nature and the large number of devices involved, a continuous approximation based on the Gram–Charlier series expansion of a truncated Gaussian kernel is proposed. We use this approximation to derive an analytic closed-form expression for the outage probability, corresponding to the event of the signal-to-interference-and-noise ratio being below a detection threshold. This metric is useful since it can be used for evaluating the performance of mMTC systems. We analyze, as an illustrative application of the previous approximation, a scenario with several multi-antenna collector nodes, each equipped with a set of predefined spatial beams. We consider two setups, namely single- and multiple-resource, in reference to the number of resources that are allocated to each beam. A graph-based approach that minimizes the average outage probability, and that is based on the statistics approximation, is used as allocation strategy. Finally, we describe an access protocol where the resource identifiers are broadcast (distributed) through the beams. Numerical simulations prove the accuracy of the approximations and the benefits of the allocation strategy.


Introduction
Machine-type-communications (MTC) have drawn a lot of attention in the past years among academic and industrial communities. They can be defined as a set of transmissions between connected terminals with no human interaction [1], which will enable the creation of a myriad of applications such as the Internet-of-Things (IoT) [2,3]. This is the reason they have become an essential part of the evolution towards future mobile communications. In fact, they are one of the objects of study in the current development of fifth-generation (5G) systems [4][5][6][7]. 3GPP standards such as long-term evolution for machines (LTE-M), also known as enhanced MTC, and narrow-band IoT, are only two examples of the impact that MTC are having on cellular communications [8][9][10]. Other standards proposed by different entities can be found in [11]. Coexistence with current systems will then play an important role in the entire progress of development of next mobile generations [12,13].
In this framework, we can distinguish a class of MTC where a large number of devices try to access the network simultaneously, the so-called massive MTC (mMTC) [14][15][16]. From now on, we consider them to be sensors that collect information from the environment and send it to a central unit (CU). Unlike typical human-to-human (H2H) communications, low complexity and high energy efficiency are preferred in mMTC instead of high data rates [15]. Thus, it is crucial to look for magnitudes that measure these figures of merit reliably and, at the same time, strategies that try to optimize them feasibly. For instance, the authors of [14] argued that non-orthogonal medium access schemes, such as sparse code multiple access, together with grant-free protocols represent a good candidate to meet the previous requirements of mMTC networks. In this work, although the motivation has the same origin, we focus on the communication aspects, which have an indirect impact on the complexity and energy efficiency.
Given that sensors transmit in a sporadic way [17], the exact transmission time of each one of them may be difficult to know in some applications (for example, in event-drive transmissions). To facilitate the analysis, particularly in the case of a massive number of sensors, we model the sporadic transmissions of different sensors as Bernoulli random variables (RVs) with known probability. As a result, we model the state of these devices as active or asleep (on/off). This intermittent behavior conditions the communication between the sensors and the CU. In particular, the received signal at the CU from each sensor is affected by a random aggregated interference coming from the other active sensors. As a result, a transmission for a given sensor can be sometimes in outage, which means the interference level is high enough to make the correct detection of the signal unfeasible. In this framework, it is then desirable to characterize the statistics of the aggregated interference to properly analyze the system performance.
In that sense, the outage probability, defined as the probability that the receiver is unable to decode the transmit message properly [18], represents an adequate metric of the system performance. It captures the random nature of the activity of the sensors, which is an intrinsic property of MTC networks, and is completely defined by the statistical distribution of the aggregated interference. In addition, note that lower values of the outage probability will lead to less retransmissions and, thus, a smaller power consumption. In addition, the energy used in this kind of systems during the idle state is very small compared to that when the device is active [10]. Hence, the outage probability is related to the network energy efficiency and its optimization can help to improve this figure of merit. Nevertheless, as mentioned above, an explicit study of the energy efficiency is beyond the scope of this paper.
On the other hand, another issue in mMTC is the medium access coordination. Given that common H2H solutions are no longer valid, a lot of effort is put into the task of coordinating the interactions in these networks. In particular, dedicated random access channels (RACH) are no longer feasible given the amount of signalling they need and the large number of these devices [19,20]. The resulting overhead, when compared to the short packet length, yields to a reduction of the overall efficiency [21]. That is why strategies that control the access in a grant-free manner seem more interesting for these systems [14]. However, these approaches entail a large amount of collisions and, thus, might lead to network failure. In this work, we focus on access methods that use orthogonal resources to reduce those collisions.
Finally, when considering that a certain set of resources can be distributed within the network, their specific allocation is still an open, yet popular, optimization problem [22]. In this work, we focus on strategies that minimize the average outage probability of the sensors. Given the large number of transmitters in mMTC systems, the solution to this problem is not trivial.

Prior Work
The statistical modeling of the aggregated interference has been studied in several cellular wireless networks considered in the literature (see, e.g., [23][24][25][26][27][28]). For example, in [23], the authors described this distribution in the context of cognitive radio with the help of the cumulants and a truncated-stable

Contributions
With the above considerations, this paper can be seen as an extension of the work presented by the same authors in the conference paper [43]. In particular, the main contributions of this work are the following: (1) Several analytic closed-form approximations of the statistics of the aggregated interference that include the random activity of the sensors. We present their advantages over the ones introduced in [43] by the same authors.
Our approach is firstly formulated for a generic scenario and is later particularized to an uplink multi-antenna mMTC setup in order to present the second novelty of this work: (2) An analytic closed-form expression for the outage probability of the sensors. It is used to evaluate the performance of this type of systems.
Finally, in this work, we also address the following issues: (3) A medium access scheme with spatial beamforming.
(4) A graph-based resource allocation strategy that minimizes the average outage probability.

Organization
The paper is organized as follows. In Section 2, the system model is described. Next, the approximations of the aggregated interference statistics are presented in Section 3. The outage probability is derived in Section 4 for an illustrative reference scenario. Next, the resource allocation problem is formulated and solved in Section 5. Numerical results are shown in Section 6 to illustrate the accuracy of the approximations and evaluate the performance of the proposed resource allocation strategy. Section 7 is devoted to conclusions.

Notation
In this paper, scalars are denoted by italic letters. Boldface lowercase and uppercase letters denote vectors and matrices, respectively. For a given vector or matrix, the operations (·) T and (·) H denote their transpose and Hermitian, respectively. Matrix I M denotes the identity matrix of size M × M. For given sets A and B, the union and intersection are denoted by A ∪ B and A ∩ B, respectively. The cardinality of A is denoted by |A|. C m×n and N m×n denote the m by n dimensional complex space and natural space, respectively. The Gaussian and the Gaussian circularly symmetric complex distributions are denoted by N (·, ·) and N C (·, ·), respectively. The Bernoulli distribution with parameter p is denoted by Ber(p). Besides, for the sake of clarity in the explanation, in Table 1, we include a list with most of the relevant variables used in further sections of this work.

System Model
In this paper, we start by considering a generic scenario in which a set of sensors transmit towards a CU through possibly different orthogonal resources. Then, only the sensors sharing the same resources will interfere each other. It is important to highlight that, under these first few assumptions, the upcoming analysis applies to many network topologies.
To model the activity of the sensors, we consider them to be in active or sleep mode (on/off) with the help of a Bernoulli RV β j ∼ Ber(p j ), where j refers to the sensor index, and p j represents the probability that sensor j is active and transmitting. Besides, we assume that the RVs corresponding to different sensors are independent. Note that this practice is commonly used in the literature to describe the sporadic nature of transmissions (cf. [23,24,28]).
At the detection stage, the received signal at the CU from each sensor is affected by an aggregated interference coming from the other active sensors, which follows a certain statistical distribution. The corresponding probability mass function (pmf) of the aggregated interference is the result of the sum of scaled independent Bernoulli RVs.
Let us consider i as the index of the sensor under analysis, the transmit signal of which is to be detected. Besides, to represent actual communication from the ith sensor perspective, we consider that this sensor is active. Hence, its activity variable β i is set to one, i.e., β i = 1.
With the above considerations, the SINR corresponding to the received signal at the CU from sensor i when it is transmitting is where a i,i is used to denote the received power of the signal from sensor i, i.e., the gain of the channel multiplied by the transmit power. The second sub-index indicates the detector in charge of detecting the signal from sensor i. In fact, σ 2 n,i corresponds to the power of the noise at the detector where the signal from sensor i is detected. Finally, the term Γ i represents the aggregated interference: where J i is the set of sensors that can potentially interfere with sensor i, i.e., those that use the same resources as the ones used by sensor i. Likewise, a j,i represents the received power of the signal from sensor j at the detector where the signal coming from sensor i is detected. The distribution of the aggregated interference Γ i can be obtained through a set of discrete convolutions given that all the individual addends β j a j,i are binary and independent [43]. However, such operation can be tedious since the complexity depends on the number of interfering devices |J i | and grows exponentially with this magnitude. Thus, it becomes rapidly unfeasible, even for a small |J i |. Even if a Monte-Carlo based numerical approach were always available, it would still carry a large computational complexity. This is why, in the following, we propose two different and reasonable alternatives to express the previous statistics using closed-form approximations. In addition, thanks to having these closed-form approximations, and not only numerical ones, an optimization of the resource allocation can be carried out, as shown below.

Approximations of the Statistics of the Aggregated Interference
The purpose of this work is to provide a closed-form expression approximating the pmf of the aggregated interference Γ i , i.e., p Γ i (γ i ). To this end, we propose two alternatives. The first one is based on the characteristic function of Γ i , whereas the second one relies on the Gram-Charlier series expansion of a Gaussian kernel. The latter is used in Section 4 to derive an analytic closed-form expression for the outage probability (i.e., the probability of not being able to decode the sensors transmit signal because the corresponding SINR at the detection stage is lower than a certain predefined threshold). This is the main novelty of the paper and represents the core of our work. To the best of our knowledge, no similar work has been done in this direction.
The approximations for the statistics of the aggregated interference that we develop in this section could be used in many different applications. In the following, we list some illustrative examples:

•
Outage Probability: Thanks to the approximation of the interference statistics, we are able to find an analytic closed-form expression for the outage probability (defined as the probability that the SINR is below a certain detection threshold). Based on that, and considering an orthogonal multiple access with limited number of available resources, an allocation scheme could be designed to minimize this metric improving, therefore, the system performance. • Throughput: The approximated statistics of the aggregated interference could also be useful to obtain an analytic closed-form expression for the throughput of the sensors (e.g., through Shannon's capacity). Note that this is indeed related to the outage probability, yet it captures different aspects of the communication.

•
Harvested Energy: Given the large number of active devices in mMTC, the energy coming from the transmitted signals can be recycled (i.e., harvested). Then, we could employ the derived statistics to characterize the amount of harvested energy. • Power Consumption: Taking into account the number of retransmissions (e.g., due to a high outage probability) and the different energy supplies (e.g., harvested energy), a power consumption model of the sensors could be derived to study and optimize the management of available energy within the network.
As already mentioned, these are only few examples of possible uses of the approximations of the interferences statistics that are developed in this work. In particular, in this paper, we concentrate on the first application, which is first discussed in Section 4 (outage probability) and then in Section 5 (resource allocation). The rest might be object of study of future works.

Characteristic Function
To obtain the pmf of Γ i , we can use the characteristic function. For X ∼ Ber(p), its characteristic function reads as where i denotes the square root of −1 (not to be confused with the sensor index i). Then, when introducing the weights a j,i with the RVs β j , the total characteristic function of Γ i results assuming independence among individuals. Since this can be interpreted as the Fourier transform of the pmf (with opposite sign), we can just invert this transformation to obtain the pmf: This way, we go from a set of convolutions to simple products and the inverse Fourier transform, which can be calculated numerically with the inverse fast Fourier transform (IFFT). The number of operations are then significantly reduced and, given that the number of points used in the IFFT for the discretization of the continuous inverse Fourier transform is actually limited, an approximation of the pmf can now be obtained. In fact, note that the necessary points for good precision in the IFFT increase with the number of interfering devices |J i |. This can be translated into a high computational cost, yet bearable in finite time. As a result, assuming that enough points are used in the IFFT routine, this method is only used to validate the accuracy of the next alternative. Finally, we can relate this approach to the moment generating function presented in [24] in the sense that it describes a similar transformation of the statistical distribution of the aggregated interference, i.e., it can be interpreted as the Laplace transform of the pmf.

Gram-Charlier Approximation
Another way to find a suitable and more computationally efficient expression for the pmf of the aggregated interference Γ i is by means of a continuous approximation. To this end, in this work, we make use of the Gram-Charlier series expansion of a Gaussian kernel [44]: where φ(γ i ; µ i , σ i ) is the probability density function (pdf) of the Gaussian distribution with mean µ i and standard deviation σ i , and κ i n , B n , and H n are the nth cumulant of Γ i , Bell, and Hermite polynomials, respectively [45].
This approach allows the approximation of a distribution through its statistical moments. Given that these magnitudes can be determined for Γ i , this inference method represents a good option for approximating its pmf. Since all β j are independent, the first two moments of Γ i are indeed The expression in Equation (6) leads to a perfect approximation for infinite addends as it converges to the actual distribution [44]. Since adding more terms reduces the error in the approximation, we truncate the series up to a finite number of addends to work with a tractable expression. Thereby, similar to the work described in [23], we obtain an expression based on the cumulants of the aggregated interference, which can be found recursively in terms of the first nth non-centralized order moments µ i,n = ∑ j∈J i p j a n j,i [46]: It is important to highlight that other kernels can be employed for the expansion. The reason for choosing the Gaussian kernel follows from the reasoning in [43], where the authors justified, through the proof of Lyapunov's central limit theorem (CLT), the use of a Gaussian distribution to properly approximate the interference statistics. In fact, note that Liesegang et al. [43] presented a particular case of the Gram-Charlier series expansion defined in Equation (6) when the order is set to 0. However, the approach in [43] (i.e., order 0) might not be adequate for some scenarios as the necessary condition might not be fulfilled (see Expression (15) in [43]). This is the case when the number of interfering devices is not sufficiently large or when sensors have different transmission probabilities, in which case Expression (15) According to Liesegang et al. [43], Equation (9) tends to zero for equal p j not close to zero and a large number of interfering devices |J i |. In that case, Condition (15) in [43] is satisfied and the Gaussian kernel with order set to 0 suffices. However, since here the probabilities are different and can be small, the simple Gaussian distribution may not be enough for good accuracy. That is why we propose an extension of that approach, i.e., an expansion of a Gaussian kernel with order higher than 0.
Until now, we have assumed that the pmf of Γ i has an infinite support. However, it is actually lower and upper bounded by 0 and J i = ∑ j∈J i a j,i , respectively. These are the extreme possible values for the aggregated interference. For analytic consistency, we must use a truncated Gaussian kernel [47]. Let X ∼ N (µ, σ 2 ) be a Gaussian RV defined between a and b with pdf in the interval a ≤ x ≤ b and 0 otherwise. The term F represents the normalization factor used to achieve unit area, i.e., . Note that φ s (·) ≡ φ(·, 0, 1) refers to the pdf of the standard Gaussian RV and Φ s (·) is the corresponding cumulative distribution function (cdf).
In our case, by defining for 0 ≤ γ i ≤ J i . Thereby, it can be shown that a similar expansion can be found for a finite support. We only need to introduce the cumulants of the truncated Gaussian distribution since, differently from before, they are not zero for orders higher than 2 and the first two do not equal κ i 1 and κ i 2 , i.e., where we define η i n = κ i n − λ i n with λ i n being the cumulants of the new truncated kernel. As mentioned, they are related to the non-centralized moments, which can be obtained in an analytic closed-form manner [47]. In particular, for a general X ∼ N (µ, σ 2 ) with support [a, b], we have where Finally, when truncating the series up to order N, we obtain the continuous approximation for the pmf of the aggregated interference, denoted by f N Γ i (γ i ). For instance, for N = 5, we have: where In the next section, this approximation is used to express the outage probability of the sensors (i.e., the probability of the SINR being below a certain predefined threshold). Then, in Section 5, this expression is used to design a resource allocation strategy.

Outage Probability
This section is devoted to presenting an analytic closed-form expression for the outage probability taking into account the massive uplink (UL) communication in mMTC systems. This figure of merit represents the probability that the receiver is unable to decode the transmit message and, thus, it can also be interpreted as the portion of time for which the communication fails. Thereby, this magnitude can be used as a valuable performance indicator in this kind of systems (cf. [28]). Furthermore, recall that lower values of the outage probability result in less retransmissions and, thus, in a lower power consumption and also in lower delays. That is why, in Section 5, for a better efficiency, we aim to find an allocation that minimizes this metric.
The outage probability is defined as [32]: where ρ i is the SINR of the signal from sensor i defined in Equation (1), and δ is the detection threshold, i.e., the minimum required SINR for the message to be successfully decoded without retransmission. Hence, we have an outage whenever the quality of the communication is not enough for correct decoding. This may happen in cases of large interference, either due to high channel gains and/or high activities of non-intended sources (represented by the set of interfering sensors J i ).
As mentioned above, we use the Gram-Charlier approximation derived in Section 3.2 to find an analytic closed-form expression for the outage probability defined in Equation (16). Note that this magnitude may be reduced when mitigating the interference. Therefore, we consider a scenario with a set of multi-antenna collector nodes (CNs), each equipped with a set of predefined spatial beams. This way, we can help the CU overcome the problem of massive access and obtain a reasonable outage probability P i out . For the sake of clarity in the explanation, we start our analysis from a simple setup, and we then sophisticate it towards the more generalized multiple CN multi-antenna scenario. Note that each scenario is a particular case of the following one, but we have decided to present the setups in this constructive way to avoid any possible confusion in the description of the notation and the developments. Thus, all previous setups are particular cases of the last scenario. In that sense, although at each step we characterize the system model, the expression of P i out is only presented for the general case in order to avoid redundancy.

Single-Antenna CU
Let us consider a single-antenna CU collecting the information from a set of M transmitting single-antenna sensors using the same access resources. In that case, the received signal is where h j ∈ C is the channel of sensor j, x j is the transmit signal with zero mean and power P, independent for each sensor, and w ∼ N C (0, σ 2 w ) is the additive noise, independent of x j . Thereby, since no further processing is performed at the CU, the received power a i,i of the signal from sensor i defined in Section 2 and the addends a j,i in the aggregated interference Γ i from Equation (2) read as and, given that all sensors share the same resources, the interfering set results J i = {j = i}. In addition, in this case, the noise at the detection stage is the same for all sensors, i.e., σ 2 n,i = σ 2 w .

Multiple-Antenna CU
To reduce the interference coming from the rest of sensors, we now consider that the CU is equipped with N antennas. The received signal is given by Given the degrees of freedom provided by the multi-antenna technology, linear processing is employed in this setup. More specifically, the CU now has L predefined spatial beams implemented through a spatial filter represented by the matrix G = [g 1 , . . . , g L ] ∈ C N×L . A possible design option is to construct the different beams so that their pointing directions are equispaced. However, for the sake of generality, we keep the filtering scheme generic and represented by G.
The signal coming from each sensor is then detected at the output of a given spatial beam. The L outputs at the L beams can be collected in a single vector given by: where y = [y 1 , . . . , y L ] T ∈ C L is the processed signal and n = G H w ∈ C L is the filtered noise. In this work, and for the sake of simplicity, we assume that the signal from a given sensor is decoded using the output signal of a single beam. Then, to determine which is the detecting beam for each sensor, a possible criterion is to choose the one leading to the largest signal-to-noise ratio (SNR) at the output of the spatial beam. Intuitively, it represents the beam where the quality of the signal might be better. Let l(i) represent the index of the beam used to detect the signal coming from the ith sensor: where g l 2 denotes the L2-norm of the filter g l . Therefore, we denote g l(i) as the spatial filter used for the detection of sensor i. Note that, if a different beam selection strategy is employed, we have a different expression for l(i). For the sake of generality, we just use l(i); hence, the upcoming analysis holds regardless of the beam selection criterion. Accordingly, we can express the term a i,i for the received power of sensor i and the interfering power terms a j,i as and, assuming that there is still a complete reuse of resources, i.e., all sensors share the same ones, the interfering set is again J i = {j = i}. In addition, the power of the noise affecting the ith sensor is given by σ 2 n,i = σ 2 w g l(i)

Multiple-Antenna CNs and CU
To further enhance the performance of the system, we now consider an extension of the previous setup with K data CNs, each one also equipped with N antennas. Each of them is responsible for collecting the information from a subset of M k single-antenna sensors for later retransmitting it to a CU, with M = ∑ K k=1 M k . Thereby, we have a two-hop communication system, as illustrated in Figure 1. Such multi-hop scheme is largely exploited in the literature [39,41] and in standards such as LTE-M [20,48]. In this work, we focus on the communication in the first hop (solid line) and leave the second stage (dashed line) for further studies. Each CN has L predefined spatial beams, represented with the spatial filter matrix G k = [g k,1 , . . . , g k,L ] ∈ C N×L , where k represents the CN index, i.e., k = {1, . . . , K}. Now, the signal coming from each sensor is detected at the output of a given spatial beam of a given CN, i.e., at a given CN-beam pair or tuple. Besides, we assume that there is no cooperation among beams and no cooperation among CNs in the signal detection stage.
The signal received at the kth CN can be written as where w k ∼ N C (0, σ 2 w I N ) is the noise vector and H k = [h k,1 , . . . , h k,M ] ∈ C N×M is the matrix containing the channels of each sensor w.r.t. each CN, i.e., h k,i is the channel between sensor i and the kth CN. Again, x ∈ C M is the vector of independent transmit signals with zero mean and power P, independent of w k , and β = diag(β 1 , . . . , β M ) is the matrix of RVs β j .
The L outputs at the L beams of the kth CN can be collected in a single vector given by: with y k = [y k,1 , . . . , y k,L ] T ∈ C L the processed signal at the kth CN and n k = G H k w k ∈ C L the filtered noise. The channels and filters are specified in Section 6 for an exemplifying scenario.
Similar to above, for each sensor, we focus on the detection at the CN-beam pair leading to the largest SNR after the spatial filter. Recall that other criteria to choose the CN-beam pair could also be used. Let now k(i) and l(i) represent the indexes of the CN and the beam used to detect the signal coming from the ith sensor using the previous criterion: where g k,l 2 denotes the L2-norm of the filter g k,l . Hence, g k(i),l(i) represents the spatial filter used for the detection of sensor i. In addition, the noise power is given by σ 2 n,i = σ 2 w g k(i),l(i) On the other hand, note that, until now, we have considered all sensors to interfere among them. In other words, we have assumed that they are all using the same orthogonal resource. Nevertheless, in real systems, more than one resource is usually available. Therefore, from now on, we consider that we dispose of R orthogonal resources with indexes {1, . . . , R}, and that they are allocated to the different CN-beam tuples.
In that sense, the sensor that is detected at a certain CN-beam pair uses the resources allocated to that tuple. For the moment, we consider that sensors know which resources they can employ, and we ignore the way this information is acquired. An example of a mechanism that provides this knowledge at the sensors side is described in Section 6.
Furthermore, we distinguish between the case where only one resource is allocated to each CN-beam pair and the situation where these tuples can use more than one resource. They are referred to as singleand multiple-resource scenario, respectively. In the latter, each sensor chooses one of the resources available for their pair, defined in Equation (25), at random. In both setups, resources are allocated to reduce the impact of the interference, as discussed in Section 5.
With the above considerations, and considering β i = 1, the received signal from the ith sensor at the tuple (k(i), l(i)) can be written as (cf. Equation (24)) where y (i) ≡ y k(i),l(i) ∈ C is the received signal; g k(i),l(i) and h k(i),i are the filter and channel of sensor i, respectively; and n (i) ≡ g H k(i),l(i) n k(i) ∈ C is the noise at that tuple, with power σ 2 n,i . Each of the y (i) signals experiences an interference I i coming from the sensors that transmit through the same resources. This interference can be decomposed in what follows: (i) the interference coming from the sensors to be detected at the same beam and CN, namely intra-beam I intra i ; and (ii) the interference coming from the sensors to be detected at the rest of beams and CNs that share the same resources, namely inter-beam I inter i . As a result, we can write where The set J k,l represents all the sensors detected at the CN-beam pair (k, l), i.e., J k,l = {j : (k(j), l(j)) = (k, l)}.
Thereby, to be consistent, we need to extract the signal from the ith sensor from the set J k(i),l(i) , as shown in Equation (28). The term t i ∈ {1, . . . , R} in Equation (28) denotes the identifier of the resource that the ith sensor is using, and the set T k,l contains the resources allocated to the CN-beam pair (k, l).
Overall, the aggregated interference Γ i from Equation (2) is decomposed into the following: where the terms Γ intra i , Γ inter i follow directly from Equation (28): Similar to above, the received power terms a i,i and a j,i read as follows Now, the interfering set J i is: At this point, we can particularize P i out from Equation (16) for the single-and multiple-resource scenarios.

Single-Resource Scenario
In this case, only a single resource is allowed per CN-beam tuple. Accordingly, the resource sets T k,l have one element and sensor i only uses a certain resource, denoted by t i , which constitutes the set T k(i),l(i) , i.e., T k(i),l(i) = {t i }. Hence, all sensors detected at the tuples (k, l) with the same resource t i , i.e., T k,j = {t i }, create interference when detecting the signal coming from the ith sensor (cf. Equation (34)).
As a result, P i out from Equation (16) is completely defined by the pmf of the aggregated interference Γ i , regardless of the resource t i , which only determines the interfering set: where we use ξ i = a i,i /δ for the sake of brevity in the notation. Note that, using the characteristic function method, a numerical approximation of P i out can be found. However, we use the continuous approximation f N Γ i (γ i ) of the pmf p Γ i (γ i ) to express the outage probability from Equation (16) in closed-form: where J i refers to the upper bound of the finite support of p Γ i (γ i ). This way, using Equation (36) together with the expression of f N Γ i (γ i ) derived in Section 3.2, we can obtain an analytic closed-form approximation for the expression of the outage probability P i out defined in Equation (16). For instance, using f 5 Γ i (γ i ) from Equation (14), i.e., set order to N = 5, the approximation in Equation (36) yields where we gather together the coefficients of the polynomials of equal order to get a more compact expression. The integral terms G i n can be found via Owen's T function [49] and the terms A i n are listed in the Table 2 below, where B i n = B n (η i 1 , . . . , η i n )/(n!σ n i ) include the set of Bell polynomials: Table 2. A i n terms used in Equation (37), listed from n = 0 to n = 5.
It is important to highlight that this approximation allows us to work directly with the statistical moments of the aggregated interference instead of the instantaneous power values a j,i . In fact, the cumulants needed for the Gram-Charlier series expansion are also obtained with these statistical moments. Therefore, the outage probability defined in Equation (36) is completely characterized by those parameters together with the term a i,i (i.e., the received power of the signal from sensor i).

Multiple-Resource Scenario
Equation (36) is valid only when a single resource is allowed per CN-beam tuple. However, when considering that multiple resources can be allocated to each tuple, we need to generalize. Up to now, β j has been a Bernoulli RV modeling the activity of sensor j. In the multiple-resource case, β j is also a RV modeling the activity of the jth sensor when it is actually creating interference (cf. Equation (34)). Accordingly, for the cases t i ∈ T k(j),l(j) , β j can be decomposed as follows: where α j ∼ Ber(p act j ) represents the event of being active and transmitting. This RV depends only on the sensor itself and is equivalent to the RV in the single-resource setup. On the other hand, τ j is a Bernoulli RV, independent of α j , that is equal to 1 whenever sensor j selects randomly the same resource that is using the ith sensor, i.e., t i . Then, in the cases that t i ∈ T k(j),l(j) , we have τ j ∼ Ber(p res j ), where p res j = 1/|T k(j),l(j) | assuming that sensors choose one of the resources within T k(j),l(j) with equal probability. Hence, τ j depends on the number of possible resources that sensor j can equally choose, i.e., 1 ≤ |T k(j),l(j) | ≤ R and, thus, on the tuple (k(j), l(j)) and the resource allocation. The case |T k(j),l(j) | = 0 is not considered since it corresponds to the situation when no resources are allocated to the tuple (k(j), l(j)). The sensors detected at that pair are not included in the interfering set since t i ∈ T k(j),l(j) in that case. Besides, all RVs τ j are assumed to be independent, including those from sensors detected at the same CN-beam tuple and that share the same parameter 1/|T k(j),l(j) |.
As a result, given that both Bernoulli RVs are independent, β j is still a Bernoulli RV with parameter p j = p act j p res j . In this case, p j represents the probability of being active and also of transmitting through the resource t i within the set T k(j),l(j) . Thus, the use of multiple resources entails a reduction of sensors activity which, in turn, reduces interference. In addition, note that, for the single-resource case, we only need to set τ j = 1 ∀j or, equivalently, |T k(j),l(j) | = 1 ∀j. That is why the multiple-resource can be seen as a generalization of the single-resource setup.
Furthermore, now interference may take place whenever there is non-null intersection between the resource sets, i.e., T k(i),l(i) ∩ T k(j),l(j) = ∅. In fact, since t i ∈ T k(i),l(i) might no longer be unique, the interfering set J i changes accordingly (cf. Equation (34)). Let us denote this sets by J i (t i ) to include the dependence with the resource t i . Therefore, although the expression of the SINR in Equation (1) is still valid when the resource t i is used, we have |T k(i),l(i) | different SINRs, one for each resource available for the ith sensor that is detected at the CN-beam tuple (k(i), l(i)).
Overall, since the reference sensor i decides equally among |T k(i),l(i) | resources, we need to include this random selection in the outage probability. To this end, we average over the different possibilities, where the resource t i changes and so does the interfering set J i (t i ): where now the RV Γ i (t i ) also depends on the resource chosen by the ith sensor from the set of resources available at the CN-beam tuple (k(i), l(i)). Note that P i out (t i ) is used to denote the outage probability when the resource t i ∈ T k(i),l(i) is used and can be approximated using Equation (36). Thus, here we also end up with an analytic closed-form approximation for P i out formulated in Equation (39). The previous analysis is of special interest in the next section, where we formulate and describe a possible solution for the resource allocation problem. Note that this task consists in deciding which resources from the set {1, . . . , R} are allocated to each CN-beam pair (k, l). Accordingly, it ultimately defines the sets T k,l . To do so, we use a graph-based approach that minimizes the average outage probability of all sensors within the network.

Resource Allocation Problem
Once the outage probability has been defined, we can design a strategy based on this magnitude to allocate the different resources in order to enhance the performance of the whole network. For this task, in the following, we present a graph-based approach that relies on the minimization of the previous approximation of the outage probability P i out . It is noteworthy to mention that the implementation of the proposed algorithm for the resource allocation is possible thanks to the analytic closed-form expression of the outage probability found in the previous section (which relies on the statistical moments of the aggregated interference).
Note that a good performance of the resulting resource allocation highlights the accuracy of the proposed Gram-Charlier approximation. Besides, given its relation to the sensors power consumption, lower values of P i out result in a better energy efficiency, which is crucial in mMTC. As mentioned, we distinguish two scenarios: singleand multiple-resource scenarios. In both setups, we seek for resource allocations that aim to minimize P i out . In particular, we adopt a strategy in which the average network performance is optimized. However, any other approach could be used, e.g., a fair strategy where a certain minimum QoS is satisfied for all sensors.
Given that the positions of sensors are assumed to be fixed, the CN-beam tuples where the sensors are detected are known. Recall that they refer to the pair with the highest received SNR according to Equation (25). To illustrate these tuples, let us consider a CN equipped with a uniform circular array (UCA). Note that this specific configuration is used here as an example, but any other structure could be used. In addition, to better understand the graphical representation, we assume the channels to be the steering vectors computed with the angle-of-arrival (AoA) of the sensors signals. A simple set of spatial beam filters is then constructed with equispaced pointing directions. It is important to highlight that, for a different array configuration and beamforming scheme, the shape of the beams (i.e., the radiation pattern) changes and so do the portions in which the whole space is divided. This is ultimately represented by the set of tuples and, thus, once they are defined, the following formulation is valid and the resource allocation mechanism remains the same. In Figure 2, the received SNR at the beams of the CN is depicted for N = L = 8, unitary transmit and noise powers, and a free-space path loss model. Sensors are deployed uniformly in a circle of radius 100 m centered at the location of the CN.

Formulation in the Single and Multiple-Resource Setups
Let us start by defining an allocation matrix C ∈ N K×L containing the resources allocated to each CN-beam tuple. Thereby, the rows and columns represent the CNs and the beams, respectively. The element [C] k,l corresponds to the resources of the pair (k, l).
In the single-resource scenario, this matrix takes values within the set {0, . . . , R}, where the zero refers to the case where no resource is allocated. Recall that here the sets T k,l have a unique element.
Hence, given that T k(i),l(i) = {t i }, the element [C] k(i),l(i) is t i . As before, t i represents the identifier of the resource that the ith sensor is employing.
We can relate the allocation matrix C to the interfering sets J i in the following way: On the other hand, a different notation must be used in the multiple-resource scenario. Now, the sets T k,l can have more than one element. To represent all the possible combinations of resources, the elements in the matrix C take values between 0 and 2 R − 1. Each element [C] k,l corresponds to the decimal value of the binary vector c k,l ∈ {0, 1} R , i.e., Thereby, the specific resources allocated to the tuple (k, l) are indicated by the positions of the nonzero elements of the vector c k,l . For instance, for R = 6 and c k,l = [101001], the element [C] k,l would be 37. In that case, we would use the first, third, and sixth resource, i.e., T k,l = {1, 3, 6}. Note that the all zero vector is also allowed as the solution might "switch off" completely some CN-beam tuples by no allocating resources to them.
Regarding the interfering sets J i , we have that where R t i contains all the elements in the set {0, . . . , 2 R − 1} that include the usage of the resource t i , i.e., Finally, since the purpose of this paper is to optimize the overall performance of the network, we look for allocation strategies that minimize the average outage probability of all the sensors. Consequently, the resource allocation can be formulated as the following optimization problem: where P i out follows the definition in Equations (36) and (39) for the single-and multiple-resource cases, respectively. Given the previous definitions, the outage probability is completely defined by the resource allocation matrix C. Note that, even though our approach has been formulated based on the objective function defined in Equation (44), any other objective function could have been considered. For instance, we could have used the maximum P i out , i.e., It is straightforward to see that the problem formulated in Equation (44) is combinatorial and that the optimal solution has an exponential complexity: O((R + 1) KN ) and O(2 RKN ) for the singleand multiple-resource cases, respectively. Hence, a brute force approach is not affordable as trying all possibilities becomes quickly unfeasible. That is why we need to seek suboptimal strategies that provide a proper solution with far lower complexity. As mentioned above, in the forthcoming subsection, we present a graph-based approach for this task, which uses coloring techniques to achieve a feasible allocation given that we have a limited number of resources.

Solution Based on Graph Coloring
One way to solve the problem stated in Equation (44) is by means of graph coloring methods. In this work, we present an approach that relies on graph structures that capture the previous setup reliably. We denote the CN-beam tuples as the nodes or vertices, and an edge or connection is established whenever two pairs can potentially interfere (i.e., whenever the sensors signals detected at a pair can potentially interfere the other pair). Accordingly, we use geometrical measures to determine whether the interference is large enough to create a connection.
The whole graph can be represented with an adjacency matrix A ∈ {0, 1} D×D [50], where D = LK is the number of nodes. Each entry [A] n,m is 1 for connecting nodes n and m.
Given that CNs might be deployed at the center of areas where the sensor concentration is high, we define a circle of radius R dep around each CN to represent these regions. Then, with the help of an interference radius R int = cR dep , we decide which of the interference coming from the sensors UL signals might be significant. The factor c ≥ 1 essentially determines how far the interfering sensors should be to be considered negligible (e.g., when the power of their received signal is orders of magnitude lower than that of sensor i).
Thereby, whenever the distance between two CNs is smaller than R int + R dep , we label them as potentially interfering CNs [38]. Note that this is the condition for intersection between the circle describing the deployment area of sensors around one CN (R dep ) and the circle representing the range of the interference coming from the sensors deployed around another CN (R int ). In addition, we denote their beams as interfering if one of them is pointing towards the other. Besides, all beams in the same CN are connected to each other. The reason is that sensors detected by a CN-beam pair are probably close to that CN. Thus, it is likely that they create a large interference to the other beams in that CN (due to the secondary lobes of the beam radiation pattern). This procedure determines A.
An example of the previous procedure is shown in Figures 3 and 4 for L = N = 4, K = 2 and c = 2. In addition, we consider a simple beamforming where the spatial filters are constructed with the steering vectors computed at equispaced pointing directions. However, the following approach does not depend on the filtering scheme. From now on, we consider that CNs use a generic beamforming. For a different scheme, a similar procedure can be used to find A.   Once the graph is created, we color it with R resources. In the single-resource case, the colors are directly the resource identifier. On the contrary, in the multiple-resource, they refer to the decimal value from Equation (41) that represents the set of resources of each pair. As a result, we have R + 1 and 2 R colors for each scenario, respectively. In both cases, we try to find an allocation such that two neighbors, i.e., connected nodes, do not share resources.
In the single-resource scenario, having different resources at neighboring nodes is equivalent to having different colors. However, in the multiple-resource, not only we need different colors, but also we need to minimize the number of resources in common that they represent. For example, for R = 6 and colors 37 and 53, the resource sets are {1, 3, 6} and {1, 3, 5, 6}, respectively. As a result, even though the colors are different, the resource sets have a non-null intersection.
That is why we need to seek strategies that avoid any reuse of resources between neighbors (and not only colors). In the single-resource case, this is known as proper graph coloring [50]. However, in the multiple-resource, standard proper coloring does not guarantee that there is no reuse (i.e., null intersections between resource sets of connected nodes). Given the limited number of resources (e.g., in LTE-M, only R = 6 resources are destined to MTC [9]), we allow two neighbors to share resources [51]. Consequently, we must introduce some criterion to choose which nodes (i.e., CN-beam tuples) can reuse resources.
The proposed solution for the optimization problem in Equation (44) is the result of an iterative algorithm similar to the well-known first-fit (or greedy) approach used for graph coloring [52]. To ease of notation, we use the vectorized form of the allocation matrix C, i.e., c = vec(C).
The first step is to order the nodes according to the number of neighbors n = [n 1 , . . . , n D ] in a descent way, where d ∈ {1, . . . , D} is the node index. This magnitude is usually referred to as degree and captures roughly the amount of interference that the nodes can suffer. It can be computed as n = Adiag(I D ). This ordering is represented with the vector o = [o 1 , . . . , o D ].
Next, we allocate random colors to each node. Afterwards, following the order given by o, we iterate over the nodes by assigning to each of them the color that leads to the minimum average outage probability. Note that this metric takes into account the non-null intersection of resources between different colors. This way, we follow the criterion in Equation (44), and the o d th element of c, i.e., c(o d ), is updated as: where W is R and 2 R − 1 in the single-and multiple-resource case, respectively. At each iteration, only one of the elements of vector c (that denoted by c(o d )) is allowed to change. The procedure is summarized in Algorithm 1, whereP (u) out is the average outage probability at the uth iteration. Note that we always follow the direction ofP out decrease and that the routine terminates when the decrease becomes smaller than a threshold χ or when a number of iterations U is exceeded.

Numerical Results
This section is devoted to present several numerical results to validate the approximation of the aggregated interference statistics introduced in Section 3.2 and, thus, to sustain the adequacy of this tool for calculating the outage probability derived in Section 4.3. Later, simulations to evaluate the performance of the allocation strategy described in Section 5.2 are shown.
In particular, we compare our approximation with the experimental results obtained with the characteristic function method from Section 3.1. Regarding the resource allocation, we compare it w.r.t. a random allocation in order to highlight the performance of our approach. For both studies, we consider the system model from Section 4.3.
On the other hand, as mentioned in Section 4.3, here we present a practical implementation of the mechanism that informs the sensors about the resources they can use. In addition, to faithfully represent a realistic scenario, we use parameters and guidelines specified by 3GPP and ITU standards. That is why we dedicate an initial subsection to discuss all these issues.

Practical Issues
In LTE/LTE-A [19,53], the smallest resource unit is the physical resource block (PRB). It corresponds to a time-frequency orthogonal resource that occupies a 0.5 ms slot and a 180 kHz bandwidth [17]. To include the coexistence of MTC systems in cellular communications in our study, we adopt the frame structure specified by that standard.
As mentioned, we assume we dispose of R available PRBs, which are allocated to the different CN-beam tuples through the graph-based approach described in Section 5.2. However, the process through which sensors know the resources they can use has not been specified yet.
Given that typical RACH based approaches are not suited for mMTC systems, we propose a methodology similar to that described in [34], where the resource identifiers are distributed among the spatial beams. Once resources are allocated to the beams, they are broadcast resource grants. Sensors detecting a certain grant, which means they are located in the beams pointing directions, use the associated PRBs to communicate. Recall that, in the multiple-resource scenario, sensors choose one of the available PRBs at random. In the event of receiving PRBs grants from more than one beam, sensors may choose that with the highest SNR (assuming normalized power per beam), i.e., that from the CN-beam tuple defined in Equation (25).
On the other hand, we consider that sensors are deployed uniformly in a circle of radius R dep centered at the CN location, and that M k is approximately equal for all CNs. The CNs are also uniformly distributed in a square area of side R t , which is set to 1 km to represent the typical dimensions of a LTE macrocell. In turn, R dep = 100 m to match those of a microcell.
To represent the multi-antenna technology used at the CNs, we use an UCA configuration [54]. In addition to that, given the low mobility of sensors [40,41], we consider channels to change very slowly and, thus, constant and known during the data transmission. For simplicity, we assume that we have a Line-of-Sight communication. A possible extension to this work could be to analyze the case where the channel varies, incorporating fading in the communication, and/or where the channel is not perfectly known.
Considering free-space propagation, the sensors UL channels are expressed using the steering vector and the path-loss coefficient. For the sake of simplicity, we consider the previous simple beamforming with L = N, where the spatial filters are constructed with the steering vectors computed at equispaced pointing directions. The factor c used to generate all the graphs is set to 2 since the interfering signals coming from distances d j ≥ R int = 2R dep are received with a sufficiently high attenuation to be considered negligible.
Finally, the probability p act j is assumed the same for all sensors and equal to an activity factor p (not to be confused with the transmit power P). Nonetheless, in the multiple-resource case, the probabilities p res j can be different as they depend on the number of resources allocated to the corresponding CN-beam tuple (k(j), l(j)). The set of simulation parameters are listed in Table 3. Table 3. Simulation Parameters. The transmit power P and the number of PRBs R are selected following the LTE-M standard [9]. According to the 3GPP indications in [17], low order constellations are used. For a QPSK modulation, a SINR of −6.7 dB is needed to achieve standard block error probabilities less than 10 % [53]. The rest of parameters (e.g., carrier frequency) can be found in [55].

Aggregated Interference Statistics Approximations
To show the accuracy of the Gram-Charlier approximation, we start by plotting the resulting cdf together with that obtained with the characteristic function. Recall that, for a sufficient number of points in the IFFT, the method of the characteristic function represents a more precise approximation. That is why, since the real pmf of Γ i is not available, this approach is used as reference. All plots are done for a certain sensor i at a random location. Note that only the single-resource case is shown as it is a sufficient representation of the setup. In comparison with the single-resource case, in the multiple-resource case only the probabilities p j change. However, since the proposed Gram-Charlier approximation already takes into account that situation (i.e., is valid for any value of the probabilities p j ), the shape of the resulting cdf in the multiple-resource does not differ from that of the single-resource. Thus, the multiple-resource case is omitted to avoid redundancy.
Results are shown in Figure 5, where we can observe the accuracy of our approximation. Orders up to 5 are presented to illustrate that the Gram-Charlier series expansion converges towards the actual pmf when using more addends. Note that order 0 corresponds to the approach described in [43], where the Gaussian kernel is employed without any expansion. Besides, the cdf of the Chernoff-based approximation from [43] is included to make a broader comparison. As expected, and following the discussion in [43], it yields a large error. In turn, our approach reveals an accurate approximation and a promising performance, specially for high orders.

Resource Allocation
In this section, we present different results to assess the performance of the allocation strategy described in Section 5.2. These simulations are used to illustrate the enhancement w.r.t. the trivial allocation where the entries of C are selected randomly between 0 and W.
We first plot the average outage probabilityP out from Equation (44) obtained with both allocation strategies when changing M and p. This is depicted in Figures 8 and 9, respectively. It can be seen that a substantial improvement is obtained with our proposal. In turn, the multiple-resource yields lowerP out as the activity is reduced and more degrees of freedom (i.e., colors) are available. Note that the resource allocation is always done with the Gram-Charlier approximation of order 5 and that the resulting probability values are computed using the characteristic function method.
On the other hand, to get richer insights, in Figure 10, we show the cdf of the outage probability: As we can see, our strategy helps to decrease considerably the outage probabilities within the network. Therefore, given the relation between the outage probability and the consumed power, our approach could also be useful to improve the energy efficiency of mMTC systems and extend the battery lifetime of sensors.

Conclusions
In this paper, we address the problem of how to model the aggregated interference statistics, which captures the sporadic activity of sensors in the context of UL mMTC. Given its discrete nature and the large number of devices, the expression of this distribution can be difficult to find. That is why we propose a Gram-Charlier series expansion of a truncated Gaussian kernel to approximate the aggregated interference statistics. Thanks to that, we derive an analytic closed-form expression for the outage probability, which is a valuable figure of merit.
In particular, we consider a scenario with several multi-antenna CNs that receive information from a group of sensors. Each of the CNs is equipped with a set of predefined spatial beams. We distinguish two scenarios, single-or multiple-resource, depending on the number of resources allocated to each beam. Since the total number of resources is limited, we present a graph coloring technique that tries to minimize the average outage probability as allocation strategy. Finally, we describe a practical mechanism where resource grants are sent through the beams in a broadcast way. Sensors located at their pointing directions will receive those permissions and use those resources to communicate. This non-dedicated scheduling approach can serve as an alternative to typical RACH schemes, which fail in the presence of massive requests.
Simulation results show that our proposal yields an accurate approximation and that our allocation method can improve the overall system performance.
Based on the justifications detailed in the paper, the random nature of the scenario comes from the sporadic activity and the high number of devices, while the channel is assumed to be known and fading is not included. Future work will focus on the extension to the cases where these assumptions do not hold. Note, however, that the performance predicted by the results in this paper can be taken as a valid benchmark for comparison purposes in those cases.