An Analytical Model for the Intercell Interference Power in the Downlink of Wireless Cellular Networks

In this paper, we propose a methodology for estimating the statistics of the intercell interference power in the downlink of a multicellular network. We first establish an analytical expression for the probability law of the interference power when only Rayleigh multipath fading is considered. Next, focusing on a propagation environment where small-scale Rayleigh fading as well as large-scale effects, including attenuation with distance and lognormal shadowing, are taken into consideration, we elaborate a semi-analytical method to build up the histogram of the interference power distribution. From the results obtained for this combined small- and large-scale fading context, we then develop a statistical model for the interference power distribution. The interest of this model lies in the fact that it can be applied to a large range of values of the shadowing parameter. The proposed methods can also be easily extended to other types of networks.


I. INTRODUCTION
In the emerging wireless communication standards LTE-Advanced and Mobile WiMAX, aggressive spectrum reuse is mandatory in order to achieve the increased spectral efficiency required by IMT-Advanced for the 4th generation of standard telephony. However, since spectrum reuse comes at the expense of increased intercell interference, these standards explicitely require interference management as a basic system functionality [1]- [3]. The research area related to the development and analysis of interference management techniques, mostly in relation with the more general subject of radio ressource management, is very dynamic, as witnessed by the high number of relevant recent contributions in this area [4]- [10]. All these new standards use OFDMA as the modulation and the multiple access scheme. In an OFDMA system, there is no intracell interference as the users remain orthogonal, even through multipath channels. However, when users from different cells are present at the same time on the same subchannel, which is the case under aggressive frequency reuse, signals superpose, leading to some form of intercell interference.
Providing statistical models of the interference power is essential to allow for an accurate evaluation of networks performances without the need for lenghtly and costly Monte Carlo simulations. The statistical characterization of the interferences has been investigated for a long time, under lots of different scenarios, and following several approaches. The distribution of cumulated instantaneous interference power in a Rayleigh fading channel was investigated in [11], where an infinite number of interfering stations was considered. In [12], the interference power statistics is obtained analytically for the uplink and downlink of a cellular system, but in the presence of large-scale fading only. Interference modeling when considering only largescale fading effects has also been investigated in [13]- [15], where the emphasis is on finding a good approximation of the lognormal sum distribution. In [16], an analytical derivation of the probability density function (pdf) of the adjacent channel interference is derived for the uplink. More recently, in [17] the pdf of the downlink SINR was derived in the context of randomly located femtocells via a semi-analytical method. Other contributions have focused directly on the analysis of a particular performance measure that is influenced by intercell interference, like the probability of outage and the radio spectrum efficiency [18]- [20]. The analysis of interference in dense asynchronous networks, such as ad-hoc networks, is also an active research area, for which a deep review of the recent developments can be found in [21], [22].
In this paper, we derive a semi-analytical methodology to estimate the statistics of the intercell interference power in a wireless cellular network, when the combined effects of large-scale and small-scale multipath fading are taken into consideration. Large-scale effects include attenuation with distance (path-loss) as well as lognormal shadowing, and the small-scale fading is Rayleigh distributed. We consider a distributed wireless multicellular network, in both cases where power control and no power control is applied. The proposed methodology is semi-analytical, in that the statistical estimate of the interference power resulting from N > 1 interferers is obtained by numerical techniques from an analytically-derived interference model for one interferer. The methodology is valid in a quite general framework; we have chosen to present it using a hexagonal network layout, although it can handle any other topology. We validate the proposed methods by comparing the moments of the estimates to the exact moments of the distribution which can be derived analytically. Using this methodology, we are able to provide a very good estimate of the pdf of the interference power, for different values of the shadowing standard deviation, σ dB . Based on these estimates, we then propose an analytical statistical model of the interference power, based on a modified Burr distribution, which includes 5 parameters. This analytical, parameterized by σ dB , model will hopefully serve as a practical tool for the assesment and simulation of wireless cellular networks when the effect of shadowing is to be considered.
The main contributions of this paper are as follows.
• In the special situation where only path-loss and Rayleigh fading are considered (no shadowing), we derive a very accurate approximated analytical expression for the pdf and the cumulative distribution function (cdf) of the intercell interference power; • We propose a semi-analytical method for the estimation of the pdf of the intercell interference power in a multicellular network when the combined propagation effects of path-loss, Rayleigh fading and lognormal shadowing are considered; • Based on this method, we derive an analytical model for the pdf of the intercell interference power by slightly modifying a Burr probability distribution. This model is parameterized by the lognormal standard deviation σ dB and its interest resides in the fact that it is valid on the whole [0, 12]-dB range of values.
The remainder of this paper is organized as follows. In Section II, we describe the multicell downlink transmission environment, and we provide the expression of the interference power for which we want to find a statistical model. In Section III, the original methodology for estimating the statistics of the interference power is presented. For this purpose, we examine in Section III-A the particular case where path loss and Rayleigh fast-fading are the only fading phenomena considered. In Section III-B, we include the shadowing effect and we consider in the first instance the contribution of one interfering cell. We then generalize to N > 1 interferers. In Section IV, we apply the proposed method to estimate the pdf of the interference power in a typical multicellular network, under two frequency reuse scenarios. Section V is dedicated to the parametric analytical modeling of the interference power. Section VI concludes the paper by summarizing the proposed methods and by presenting some perspectives.
We will use the following notation for the rest of the paper. Non-bold letters such as x are used to denote scalar variables, and |x| is the magnitude of x. Bold letters like x denote vectors. We use E {X} to denote the expectation of X. The pdf and cdf of the random variable (r.v.) X will be denoted p X (x) and F X (x) respectively.

II. MULTICELL DOWNLINK TRANSMISSION MODEL
We consider the downlink of an OFDMA-based 19-cell cellular network having the 2D hexagonal layout depicted on Fig. 1. We assume a unit-gain omnidirectional SISO (single input, single output) antenna pattern, both for the fixed access points (APs) and the mobile user terminals (UTs) which are supposed to be uniformly distributed over the service area. As OFDM is used for intracell communication, we assume an orthogonal transmission scheme within a cell. We consider a synchronous discrete-time communication model in which active APs at any given time slot send information symbols to their respective UTs over a shared spectral resource, which gives rise to an interference-limited environment. In this framework, we will focus on the statistics of the so-called intercell interference power undergone by a typical UT. In this regard, we will consider UT in cell 0 (denoted UT 0 , see Fig. 1), for it is surrounded by 18 potential interferers. For UT 0 , the received signal on OFDMA subchannel at time slot m can be modeled as Here x 0 (m, ) represents the information symbol intended to UT 0 and x n (m, ), n = 0, the nth interfering symbol (this symbol is sent from AP n to its respective user). The coefficient h n (m, ) denotes the instantaneous gain of the th (interfering) subchannel from AP n to UT 0 . Each subchannel is subject to additive white Gaussian noise w (m, ). In the following, we will focus without loss of generality on a single OFDMA subchannel, thereby omitting subchannel index in all subsequent notations.
Two frequency reuse scenarios will be considered (see Fig. 1): • the full frequency reuse pattern, denoted FR1, where all APs in the network transmit at the same time using the same frequency range (N = 18 intercell interferers); • a partial frequency reuse pattern, denoted FR3, with reuse factor 3 (N = 6 interferers). Each channel is assumed to be flat-fading, possibly experiencing small-scale multipath fading and/or large-scale effects. For the rest of the paper, we concentrate on the instantaneous channel power gain 1 G n (r n ) , which is proportional to |h n (m)| 2 and can be expressed as a three-factor product: G n (r n ) = G pl,n (r n ) G f,n G s,n , n = 1, 2, . . . , N.
In the above equation, r n denotes the distance between UT 0 and AP n (distances r n are functions of UT 0 's position within its cell). G pl,n (r n ) = K (1/r n ) γ is the (deterministic) path loss (normalized with distance, see Appendix A), where K is a constant and γ represents the path loss exponent. The Rayleigh fading gain G f,n is modeled by an exponential distribution with rate parameter equal to 1, i.e., E {G f,n } = 1; we denote the corresponding pdf by p G f,n (x). The shadowing gain G s,n is modeled by a lognormal distribution whose pdf can be written where ξ = 10/ ln(10) [23]. Note that the importance of the shadowing phenomenon is directly related to the standard deviation σ dB . For a given σ dB , the parameter µ dB is determined to ensure a unit mean shadowing gain: E {G s,n } = 1, which leads to µ dB = −σ 2 dB / (2ξ). As r.v.'s G f,n and G s,n are independent from each other, and as E {G f,n } = E {G s,n } = 1, we have, from (1), E {G n (r n )} = G pl,n (r n ), which reflects the fact that the nth interfering channel's Rayleigh fading and shadowing components cause the actual gain G n (r n ) to fluctuate about its mean value G pl,n (r n ).
The total interference power undergone by UT 0 can then be written I = N n=1 P n G n (r n ), where P n = E |x n | 2 is the power emitted by AP n. In what follows, we consider that all APs transmit at the same power, i.e., P n = P for all n. This corresponds to, e.g., a fast-fading environment where no channel state information feeds back from mobile users to APs, which results in a no power control scheme where all APs transmit at the maximum power; although crude, this scheme can be seen as a lower bound on performance for real systems. Considering that each AP transmits at the same power P also applies to a more practical scenario where APs have access to channel state information and power control is associated with the opportunistic scheduling policy proposed (and proved to be sum-rate optimal) in [10], when the number of users per cell is high (since in this case, it can be expected that the channels between users scheduled at the same time and their serving APs have about the same power gains). Thus, the interference simplifies to I = P N n=1 G n (r n ). We now define the interference gain -which will be denoted G -as being the sum of the channel power gains between the interested user and the N interferers, i.e., (Note that G is a function of UT 0 's location through the distances r n .) So, as I = P G, characterizing the interference power I is equivalent to studying the interference gain G. We will concentrate on the latter in the subsequent sections.

III. METHODOLOGY
We are now interested in finding an estimate of the pdf of the random interference gain (2). Since direct calculation of the pdf does not seem possible, we aim at producing an accurate histogram for the interference gain G that will then be modeled using a specified statistical distribution. Such a histogram is constructed from a set of samples called a typical set, i.e., a discrete ensemble of values that accurately represents a random phenomenon. Traditionally (and especially in the telecommunications area), this typical set is issued from Monte Carlo simulations, which might, at first sight, produce satisfying results. However, in a propagation environment that is subject to intense shadowing (i.e., for large values of the [0, 12]-dB range under consideration), the classical Monte Carlo method fails at producing a representative set of sampled gains [24], [25]. This can be explained by examining the particular distribution involved, for one single as well as for multiple interfering cells. A typical cdf of the interference gain (single or multiple interferers) for a high value of σ dB belongs to the class of heavy-tailed distributions [26], for which the least-frequently occuring values -also called rare events -are the most important ones, as a proportion of the total population, in terms of moments. A finitetime random drawing process performed on this cdf never produces these rare events because of their very low probabilities, which causes the resulting set to be not typical. Hence the need for a new approach.
As will be seen in Subsection III-B2, the pdf and the cdf of the interference gain for one single interferer may be expressed in its integral form. From this expression, we propose the following two-step approach: 1) Produce a typical set of gains for one interferer using the generalized inverse method.
This method consists in generating a typical set of samples corresponding to an arbitrary continuous cdf F , and is based upon the following property: if U is a uniform [0, 1] r.v., then F −1 (U ) has cdf F ; 2) Produce a typical set for multiple interferers by adequately combining typical sets from single interferers and the Monte Carlo computational technique.

A. Special case: No shadowing
We start this section by considering a propagation environment in which the only fading phenomenon is due to Rayleigh multipath fading. In this particular case, (1) simplifies to ( We first note that, because of the symmetry of the network geometry, we need only study the interference power distribution for UT 0 located within one of the twelve triangular sectors depicted on Fig. 2; in the following, we will consider the grey-shaded region for illustration purposes. We now introduce an original approximation that will help simplify further computations. We can see that in (3), it is UT 0 's random position that makes the path loss G pl,n (r n ) fluctuate, when the randomness of G f,n is due to Rayleigh fading. But it is worth noting that, although both phenomena are random, path loss fluctuations differ from multipath fading in an important way : the path loss takes values in a finite set (related to UT 0 's location within its cell) whereas the variations due to fading have an (theoretically) infinite dynamic range. Since pathloss fluctuations' dynamics are very small compared to fading's, we propose to approximate (3) by replacing each gain G pl,n (r n ) by its average value, which leads to using the notation r n = f n (r 0 , θ), n = 1, 2, ..., N , where (r 0 , θ) are UT 0 's polar coordinates, as depicted in Fig. 2. By examining (4), we see that, under this approximation, G n does not depend on UT 0 's varying position anymore.
We further note that G n , as expressed in (4), is an exponentially distributed r.v. with rate parameter 1/λ n [27], λ n -which we call the average path loss -being defined as follows: Using (5), (4) can also be written and the intercell interference gain (2) can be reduced to a sum of independant (but not identically distributed) exponential r.v.'s: G, as expressed in (7), is a r.v. whose cdf, denoted F G (g), has a closed form expression available in the literature [28]; it can be expressed as where The pdf, denoted p G (g), can be easily calculated by deriving (8): In Section IV-A, it is first shown that approximation (4) is valid in the case of one single interfering cell. This consequently validates the proposed model (7) in the case of multiple interfering cells, which we show for both frequency reuse patterns FR1 and FR3.

B. General case: Attenuation with distance, shadowing and multipath fading
Let us now focus on characterizing the distribution of the intercell interference gain G in a propagation environment where Rayleigh fading as well as shadowing (due to obstacles between the transmitter and receiver that attenuate signal power) are taken into account. To the best of our knowledge, no closed form expression for the interference gain G exists in the literature. But, as will be seen in Section III-B2, we determine an analytical formula (under integral form) of the distribution of the interference gain for one interferer. Using this result, we are able to obtain a histogram for G's distribution in the presence of multiple interferers.
For this purpose, we proceed in two steps: first, we compute a typical set for the interference gain produced by one single interferer. As described in Section III-B2, this is done by numerical computation (from the integral-form cdf), followed by non-uniform partitioning, and then inversion, of the cdf. Then we generate a typical set for N interferers using an appropriate combination of the (weighted by λ n ) typical sets of each single interferer (Section III-B3). The accuracy of the proposed method will be evaluated in both single-and multiple-interferer cases by comparing the actual moments computed from the typical sets with the exact moments of the interference gain distribution (which can be formulated analytically, as will be seen in Section III-B1). 1) Preliminaries: We begin this section by examining two important points. When taking into account multipath fading as well as shadowing as the fading effects in the propagation environment, a question arises about the validity of the original approximation (6). Fortunately, our approximation is being strengthened by this additional contribution due to shadowing, since this phenomenon is just another source of infinite-dynamics randomness. Taking shadowing into consideration amounts to introducing an additional term in (6) that can now be written A second point pertains to the moments of both statistical distributions of G n (single interferer) and G (multiple interferers). Using approximation (10), it is shown in Appendix B that the kthorder moment of G n 's distribution has the following expression: Computation of the kth-order moment of G's distribution is done in Appendix C and leads to the following formula: where a = (α 1 , α 2 , . . . , α N ), α n ∈ N, n = 1, 2, . . . , N , is an N -dimensional vector whose sum of components is written |a| = N n=1 α n , and λ α = λ α 1 1 λ α 2 2 . . . λ α N N . So the summation in Eq (12) is taken over all sequences of non-negative integer indices α 1 through α N such that the sum of all α n is k. Note that the 1st-order moment, is a quantity of particular interest because it is proportional to the average power of the interference signal.
As closed form expressions of moments have been determined, they may be used in evaluating the accuracy of typical sets for both single-and multiple-interferer statistical laws.
2) Single interferer: We now turn on to computing a typical set for the interference gain produced by one interferer. For convenience, the average path loss (5) for this single interferer is normalized to 1, i.e., λ n = 1, so (10) reduces to As G n is the product of two independent r.v.'s, its cdf can be written where F Gs,n (x/u) denotes the shadowing gain's cdf. Recalling that G s,n is modeled as a lognormal r.v., we have, using the same notations as in Section II, Replacing p G f,n (u) and F Gs,n (x/u) by their respective expression in (15), we obtain an integral-form expression for the cdf of the intercell interference gain produced by one single interfer: We are now interested in generating a typical set of the interference gain G n ; we denote this typical set by S n , where is the number of elements in the set. It was mentioned in Section III that, though widely used in telecommunications, the Monte Carlo computational technique proves inefficient for large values of σ dB . An interesting alternative method is the generalized inverse method, for which an -element typical set for a given distribution is obtained by an -level uniform partitioning, followed by inversion, of the cdf. Now we know that, for large values of σ dB , the distribution of G n exhibits the heavy-tailed property, which means, as described before, that the least-frequently occuring values (i.e., the highest gains) are the most important ones in terms of moments. Therefore, taking these highest amplitudes into consideration using the 'classical' generalized inverse method would require a finer partitioning of the cdf, which would produce a typical set made up of a huge amount of elements.
In order to construct a typical set with a reasonable value for , we propose to accomodate the above-mentioned method by performing a non-uniform partitioning of G n 's cdf, and, as high amplitudes are important in terms of moments, we proceed with a finer partitioning of the [0, 1] segment for values close to 1. The implementation details of the method are described on Fig. 3; they result from a good compromise between accuracy and simplicity. We first divide the interval [0 1] of the cdf into J intervals, numbered j = 1, . . . , J, of different lengths: the jth interval has a length d j = 9 × 10 −j , j = 1, . . . , J − 1; the last interval has a length d J = 10 −J to ensure J j=1 δ j = 1. We next perform a P -level uniform partitioning on each interval, i.e., each interval is now partitioned by P equally-spaced points. Finally, we invert the partitioned cdf to obtain a typical set S n of cardinality = J × P . Also, as the proposed partitioning is nonuniform, S n needs to be associated a probability set: the probability of an element computed from the jth interval is δ j = d j /P . It can be shown (see Section IV-B) that using J = 25 intervals containing P = 900 points each -which results in a typical set that contains only = 25 × 900 = 22, 500 elements 2 -guarantees that up to third-order moments derived from the typical set are within 1% of the exact values for all σ dB .
3) Multiple interferers: We now focus on finding an L-element typical set -denoted S Lfor the interference gain G that must be computed from N typical sets S n , n = 1, 2, . . . , N .
We first note that interferer n's typical set can be directly obtained by weighting each element of S n by its average path loss λ n ; we will denote interferer n's typical set by λ n S n . Let us now find a way to produce the ensemble S L from the typical sets λ n S n .
Ideally, S L should be constructed by considering all combinations of the elements of the typical sets λ n S n , but the cardinality of the resulting set, L = N = (JP ) N , would rapidly become prohibitive as the number N of interferers increases.
To get rid of this complexity, we point out that the above-mentioned ideal (exhaustive) solution can also be viewed as an exhaustive combination of intervals (J N combinations) associated with an exhaustive combination of elements within each interval combination (P N combinations). And we observe that the most important part of this exhaustive solution pertains to the combination of intervals, i.e., the combination of elements belonging to interval j of typical set λ n S n with elements belonging to interval k, k = j, of typical set λ m S m , m = n. So a way to construct a (near optimal) typical set for G could be to perform exhaustive combinations of the intervals (as in the exhaustive solution), and to approximate the exhaustive combination of the elements within each interval combination by the following procedure: for each of the J N combinations of N P -point intervals, • Perform a random permutation of the P elements within each of the N P -point intervals 3 ; • Add up these N permuted P -point intervals to obtain one resulting P-element interval. This last P -element interval approximates the P N -element interval that would have resulted from an exhaustive combination of elements within the considered interval combination. Now, as there are J N interval combinations, the resulting typical set would contain J N P elements, which can still be prohibitive, so this second solution -which we will refer to as the near-optimal solution -can not be applied as such.
We eventually propose a novel approach which makes use of this near-optimal solution and is based on the following two-step algorithm: Step 1 Apply exhaustive combinations of intervals to a subset of M interfering links; Step 2 Perform Monte Carlo simulations for the N − M remaining links.
We now detail the principle of the proposed method. In Step 1, we apply the near-optimal solution described above, but to a subset of M < N interfering links which we will call compelled links. The compelled links are chosen to have the highest average path losses (λ 1 ≥ · · · λ M ≥ · · · λ N ) so as to minimize errors in other (non-compelled) interfering links. The exhaustive combination of the J intervals for M compelled links obtained from the near-optimal solution thus results in one set of J M P elements. In Step 2, we build up a J M P -element set for each of the N − M remaining, non-compelled, links by performing J M random drawings of intervals according to the probability set {δ j }, j = 1, 2, . . . , J. As in the near-optimal solution, a random permutation of the elements is applied at each drawing. The ensemble of amplitudes of the intercell interference gain G -the so-called typical set S L -is then constructed by adding up these N − M + 1 sets; it is of cardinality L = J M P . Associated to S L is a probability set determined as follows: to each interval is associated a weight which is the product of probabilities δ k of intervals issued from compelled links (for non-compelled links, probabilities are accounted for by means of the random selection process); these weights are then normalized to obtain probabilities. Finally, the histogram of the interference gain G can be constructed from these resulting amplitude and probabiliy sets. It is important to note, however, that, as a random drawing process is involved, a number of iterations might be needed in order for this process to converge (elements of S L and associated probabilities are averaged at each iteration). We will call this semi-analytical technique the Monte Carlo-panel method (MCP, in short) 4 .
The MCP method is illustrated on Fig. 4 for N = 4 interfering cells, M = 2 compelled links, and J = 2 intervals per typical set (these intervals -denoted A and B -have probabilities δ 1 = 0.9 and δ 2 = 0.1 respectively, and each one of them contains P elements).
Step 1 of the alogrithm is summarized in the light-grey shaded box: intervals from typical sets S 1 and S 2 (corresponding to compelled interfering links 1 and 2, and weighted by their respective average path losses λ 1 and λ 2 ) are combined together, as described in the near-optimal solution, to obtain a set of amplitudes of cardinality 4P representative of the two compelled links; associated to this set of amplitudes is a set of weights {0.81, 0.09, 0.09, 0.01}. The dark-grey shaded box summarizes Step 2: for each non-compelled interfering link, a 4P -element set of amplitudes is made up by 4 intervals (A or B) drawn according to the probability set {0.9, 0.1} and applied random permutations. The typical set S L (with L = 4P in our example) is then obtained by summing up together all these sets. The histogram of the interference gain G is constructed from S 4P and the associated probability set 5 . Note that one random permutation of the interval (permuted intervals have been assigned the prime symbol) is performed at each (compelled or random) manipulation of an interval.
Implementing the MCP method however requires cautiousness. In non-compelled links, random drawings of intervals are performed based on the probability set {δ j }, j = 1, 2, . . . , J. In this process, lowest-probability intervals, which contain the highest interference gains, are totally ignored for two reasons. The first reason pertains to the fact that obtaining a significant frequency of appearance of such rare events would require a prohibitive number of simulation runs. The second reason is due to limitations inherent to software simulation tools which use pseudo-random number generators to generate sequences of 'random' numbers belonging to a fixed set of values. In order to take into account the ignorance of the contribution of the highest interference gains of the N − M non-compelled interfering links in the probability set {δ j }, we suggest the following workaround: in these links, we intentionally make exclusive use of the J , 1 ≤ J < J, first intervals, and we associate them a loaded probability set δ j defined as follows: is a normalizing constant such that J j=1 δ j = 1 (using the particular non-uniform partitioning described previously, we have: α = 1/ 1 − 0.1 J 1). Now, as was mentioned before, high amplitudes play an important role in terms of moments. Although the impact of neglecting them in non-compelled links is globally limited because these links are weighted by smaller average path losses λ n (n = M +1, . . . , N ), it has to be compensated in order to satisfy the 1st-order moment constraint (i.e., the sampled mean has to converge to the exact value 6 ). For this purpose, small (resp. large) amplitudes need to be underweighted (resp. overweighted). Thus, an underweighting multiplicative factor, denoted f − , is applied to amplitudes of the J first intervals of compelled links; similarly, an overweighting multiplicative factor f + is applied to amplitudes of the last N − J intervals. (Computation details of factors f − and f + are given in Appendix D.) Let us last notice that the choice for values of M and J is a trade-off between differents aspects: cardinality of the resulting typical set (i.e., tractable number of points), number of simulation runs and accuracy of the histogram. We have determined that M = 2 and J = 3 meet all these 5 The probability set is obtained by normalizing the set of weigths. 6 We recall that the mean E {G} is of particular importance because it is proportional to the average interference power. requirements.

IV. NUMERICAL RESULTS
In this section, we present numerical results related to the different methods introduced in the preceding section. In Section IV-A, we first examine the validity of the original approximation introduced in Section III, stating that the interference gain G n (and, consequently, G) does not depend on the user's position within its cell. For this purpose, we compare the approximation of G given by (6) with the 'exact' formula (3). Then, in Section IV-B, we obtain the histogram of the interference gain G n (one single interferer) by applying the non-uniform partitioning generalized inverse method described in III-B2. Finally, the MCP method (see III-B3) is used to build up the histogram of the interference gain for multiple interferers in Section IV-C.
We use the following simulation parameters. We consider a system functioning at 1 GHz. We fix the cell radius to R = 700 m, d 0 = 10 m, and the pathloss exponent to 3.2, which corresponds to a typical urban environment, as described in the COST-231 reference model [29]. The reference distance is chosen to be equal to 2R. Average path losses λ n , n = 1, 2, . . . , N , are determined numerically using (5) and are summarized in Table I.

A. No shadowing
In this section, we evaluate the proposed approximation (6) against Monte Carlo simulations performed on (3). We first consider the contribution of one interfering cell and, in this regard, we examine two opposite scenarios: one for which the investigated interferer (i.e., AP 1) produces the largest dynamic range for the intercell interference power undergone by a user in the greyshaded triangular area of Fig. 2; the other one for which the investigated cell (i.e., AP 13) has the smallest dynamics. Obviously, both dynamics differently impact the accuracy of our model. Note that, in both cases, the sum of interference gains (7) reduces to one exponential r.v. Modeled and simulated pdf's for above-mentioned cases (a) and (b) are plotted in Fig. 5 and Fig. 6 respectively, and the good match of the curves shows that the proposed method is a good approximation.
We then consider the whole set of interfering cells (N interferers) under frequency reuse patterns FR1 and then FR3, for which results are shown in Fig. 7 and Fig. 8 respectively. We see that simulated and modeled probability laws (2) and (7) respectively closely match for both frequency reuse patterns. We also note that simulated and approximated curves are closer to one another for FR3 than they are for FR1. As explained before for the single-interferer scenario, fluctuations of actual pathlosses G pl,n (r n ), n = 7, ..., 18, can be assumed to have about the same dynamic range, but these dynamics are smaller than those of gains G pl,n (r n ), n = 1, ..., 6.

B. Shadowing, one interferer
In this section, we make use of the non-uniform partitioning generalized inversion method introduced in Section III-B2 to obtain a typical set for the interference gain of one interferer. Table II presents the three first moments computed from typical set S n , as compared with the exact moments of the distribution of the interference gain G n . We see that moments issued from the typical set are far beyond the 1% accuracy requirement. The proposed method also outperforms the Monte Carlo simulation technique, which cannot be guaranteed to converge for such a small number of points.
Histograms of the interference gain G n computed from typical set S n is illustrated on Fig. 9 for different values of σ dB .

C. Shadowing, multiple interferers
We now evaluate the MCP method developed in Section III-B3. We have determined that 20, 000 iterations of the base MCP algorithm guarantee that the 1st-order moment computed from any typical set (whatever σ dB value is considered) converges to its exact value (13). Table II presents the values of the 1st-order moment of G, both exact (analytical) and approximated (computed from the typical set). We can see that the proposed method performs very well for the whole range of σ dB values.
Histograms of the interference gain G computed from typical sets obtained by the MCP method are illustrated on Fig. 10 (FR1 scenario) and Fig. 11 (FR3 scenario) for different values of σ dB .

V. STATISTICAL MODEL
In Section III, we developed analytical and numerical methods to build up a good approximation of the histogram of the interference gain G. In this section, we aim at using this result to elaborate a statistical model for G, i.e., a closed form expression of the probability law, characterized by the shadowing parameter σ dB . This task is challenging in that one single parametric law is required, that is valid for propagation environments which considerably vary depending upon the shadowing phenomenon (parameter σ dB ), and that is applicable to various frequency reuse scenarios (FR1 and FR3).
We initialize the modeling process by extracting usefull information from a carefull analysis of the histograms of the interference gain G (see Fig.'s 10 and 11). We first note that G is a positive continuous r.v. We then observe that all curves are asymmetric, and this property is even more pronounced for large values of σ dB . In this case, G's pdf's also have a sharper peak and a longer, fatter tail, the last of which being a characteristic of heavy-tailed distributions (a.k.a. power distributions), as already mentioned.
Due to the strongly skewed nature of the interference gain distribution for large σ dB 's, a power-type statistical model turns out to be suitable here. In this regard, a Pareto-like distribution seems to be a good candidate, so we focus, in first approximation, on a 3-parameter Burr-type XII distribution [27]. The Burr distribution has a flexible shape and controllable location and scale, which makes it appealing to fit any given set of unimodal data that exhibits a heavy-tail behavior (e.g., it is an appropriate model for characterizing insurance claim sizes). However, as 3 parameters seem to not be sufficient to correctly characterize the interference gain distribution under those particularly tight constraints, another law is required, which offers greater flexibility to match the whole range of σ dB values. Such a flexibility is provided by introducing an additional shape parameter into the Burr distribution, based on the following property [30]: if F (x) is a cdf, so is (F (x)) η , ∀η > 0. Thus, we have established a new Burr-based probability law, whose cdf -denoted F G (x) -is given by where η, α and k are the shape parameters, and β is the scale parameter of the distribution. G's pdf -denoted p G (x) -can be easily obtained by deriving (19): We next establish a parametric family of functions (parameterized by σ dB ) for the interference gain G by determining empirical formulas for parameters η, α, k, and β. For this purpose, we propose that all parameters (whatever frequency scenario is considered) be modeled by the same 6-parameter function f that has the following expression: where coefficients a i , i = 1, 2, . . . , 6 have been determined empirically and are summarized in Table III. Corresponding empirical laws f , as functions of σ dB , are plotted on Fig. 12 (FR1 scenario) and Fig. 13 (FR3 scenario). The pdf's of the proposed statistical model are superimposed on histograms obtained by the MCP method for different values of the shadowing parameter σ dB on Fig. 14 (resp. Fig. 15) for the FR1 (resp. FR3) scenario. We now come to the last step of our modeling process. As seen earlier, MCP-obtained histograms and the proposed Burr-based distributions closely match for the whole range of σ dB . However, care must be taken in defining the range of gains for which our model is valid. And indeed, the Burr-based statistical law needs to be truncated at a maximum value -denoted x t -defined in such a way that the 1st-order moment constraint holds, which we can write where E {G} is the exact mean (13). As a consequence of this truncation process, a normalizing factor, has to be incorporated in both the cdf and pdf of the elaborated model, which are then written AF G (x) and Ap G (x) respectively. Regarding the empirical law x t as a function of σ dB , we also propose the same 5-parameter function for both FR1 and FR3 scenarios: where coefficients a i , i = 1, 2, . . . , 6 have been determined empirically and are summarized in Table IV. Empirical laws x t , as functions of σ dB , are plotted on Fig. 16 (FR1 scenario) and Fig. 17 (FR3 scenario). The normalizing factor may be easily computed by replacing x t by its actual value in (22).

VI. CONCLUSION AND FUTURE WORK
In this paper, we have proposed a methodology to estimate the statistics of the intercell interference power in the downlink of a multicellular network. In a propagation environment subject only to path loss and multipath Rayleigh fading, we have established an accurate approximated analytical expression for the interference power distribution. Then, considering the combined effects of path loss, lognormal shadowing and Rayleigh fading, we have proposed a semi-analytical method for the estimation of the pdf of the interference power. Finally, we have developed a statistical model parameterized by the shadowing parameter σ dB and valid on a large range of values ([0, 12] dB). It is our hope that the methods described in this paper are sufficiently detailed to enable the reader to apply them to other types of environments.
A future work will pertain to improving the statistical interference power model by more closely linking the proposed model developed for a combined Rayleigh fading-lognormal shadowing environment to the 'exact' analytical formula obtained in the case where only Rayleigh fading was considered. Another perspective is to apply the proposed methods to other wireless network topologies (e.g., ad hoc networks,...).

A. Normalized channel power gain
In this paper, we concentrate on the channel power gain H n (r n ) = |h n (m)| 2 , where h n (m) is the instantaneous gain of the channel between AP n and UT 0 . H n (r n ) can be expressed as a three-factor product:

B. Computation of moments for one interferer
We find the closed form expression of the kth-order moment E (G n ) k of the statistical distribution of the interference gain G n (one interfering cell). We have: where (30) follows from the independance property of the r.v.'s G f,n and G s,n . As G f,n is exponentially distributed with unit mean, its kth-order moment is given by: As for G s,n , it has a lognormal distribution with parameters −σ dB /2 and σ dB ; its raw moment can be written: Replacing (31) and (32) in (30) leads to (11).

C. Computation of moments for multiple interferers
We establish the analytical formula of the kth-order moment E G k of the statistical distribution of the interference gain G (multiple interferers). Using approximation (10), we can write: where the following notation is used: • a = (α 1 , α 2 , . . . , α N ), α n ∈ N, n = 1, 2, . . . , N , is an N -dimensional vector whose sum of components is |a| = N n=1 α n ; • the multifactorial a! is such that (α n !) ; • the variable Z a is defined as follows: Using (30), we can further develop (33), which gives (12).

D. Computation of correction factors
We determine the correction factors used in the MCP method described in Section III-B3. Recall that the technique consists, for non-compelled links, in randomly selecting intervals from a subset containing only the J highest-probability (i.e., smallest-amplitude) intervals. But, as high-amplitude intervals never appear in this random process, small amplitudes get overweighted in non-compelled links, which must be compensated in compelled links, where small (resp. large) amplitudes need to be underweighted (resp. overweighted), in such a way that the 1st-order sampled moment converges to its exact value. Thus, in order to satisfy the mean constraint, an underweighting multiplicative factor, denoted f − , is applied to amplitudes of the J first intervals of compelled links; similarly, an overweighting multiplicative factor f + is applied to amplitudes of the last N − J intervals. We now compute these two correction factors.
Let us first see how each interfering link contributes to the 1st-order moment of the intercell interference gain G. For each compelled link n, n = 1, . . . , M , we can write 7 : where G n = G f,n G s,n (approximation (14), with λ n = 1), and, by construction of the typical set S n , A + B = 1, ∀σ dB . For each non-compelled link n, n = M + 1, . . . , N , G n 's mean is where the probability set δ j is given by (17). So, if no correction factors are introduced, the contribution of all (compelled and non-compelled) links to the intercell interference gain G gives the following mean: is the exact mean (13).
Let us now introduce the correction factors f − and f + into compelled links, as described previously. G's 1st-order moment -denoted E cor {G} -then becomes: In order for both exact and actual means to be equivalent (i.e., (34)≡(35)), we need to solve the following system: which leads to Note that we have f + > 1 and, as α 1, f − 1.     Fig. 2. Because of the particular symmetry of the network geometry, we need only study the interference gain distribution for a user located within one the twelve dashed triangular areas. For illustration purposes, we will consider the grey-shaded sector.