Achille Mouaffo

Leonardo Cardoso

Hervé Boeglen

Guillaume Villemaud

Rodolphe Vauzelle

Radio link characterization of the CorteXlab testbed with a large number of software defined radio nodes

Keywords: propagation, measurement, software defined radio. I

This paper presents the first implementation on software defined radio nodes in the large scale testbed called CorteXlab of a radio link estimation technique based on OFDM transmissions. The purpose of this large scale testbed is to offer to the whole scientific community an open tool to test new techniques for multi-user, cooperative and cognitive radio networks in a controlled environment. As the experimentation room was defined in order to offer reproducible measurements, it is important to be able to characterize each radio link between all transceivers. Therefore, we present here the development of a channel sounder directly implemented on the software radio nodes.

INTRODUCTION

In the future internet of things, wireless technologies will represent a large part of the market in weak competition with power line or infrared communications. Up to now, standards have been developed with a bottom-up approach and the radio spectrum sharing policies are mostly static. Since billions of objects are expected to use wireless links, the present way the wireless medium is shared has to be revisited. For such, radio systems, algorithms and protocols have to be deeply transformed.

The most promising approach stems from the cognitive radio paradigm, which relies on three complementary mandatory properties: all radio systems real-time reconfigurability, wireless environment awareness behavior, and self-organization capability. The first item is referred to as software defined radio (SDR), the second as cognitive radio (CR) and the third as self-optimization networks (SON).

In this paper, we will present the CorteXlab testbed [1], a large scale testbed open to the scientific community to test, evaluate and optimize current or future wireless systems, particularly focusing on multi-user, cooperative and cognitive networks. We will first present the globality of this testbed, then detail the problem and implementation of a channel sounder on this kind of SDR nodes, and finally present some preliminary results and draw some conclusions.

II. PRESENTATION OF THE CORTEXLAB TESTBED

A. Motivation for this large scale testbed

A valuable platform able to test realistic future scenarios, should offer the possibility to evaluate simultaneously these three items to support the actual theoretical developments from an experimental point of view. Unfortunately, a platform that contemplates all of these topics at the same time does not currently exist. CorteXlab aims to fill in this gap. The CorteXlab testbed will be hosted at INSA-Lyon in France, benefiting from the Senslab (now FIT/IoT-lab) experience and from the 5 years experience on developing a MIMO (Multiple Input Multiple Output) high data rate reconfigurable platform. CorteXlab will use the network architecture developed in IoTlab and will integrate SDR nodes to offer a remotely accessible development platform for distributed Cognitive Radio (CR). Reconfigurability, compatibility, coexistence and even cooperation between SDR nodes will be evaluable. A large set of heterogeneous SDR nodes (MIMO nodes, SISO nodes and Wireless Sensor Network (WSN) nodes) will permit a full experimental evaluation.

CorteXlab will allow remote users to test their own algorithms on the existing nodes, but the architecture will be also opened to industry third party to deploy their own frontend (RF or UWB) or baseband systems to test and validate their developments. A clear expected result is offering a remote access to all equipments in a comprehensive way, such that many scenarios can be evaluated by remote users. A second target is to create a "network of SDR" development community including people from digital communications, networking and embedded systems, to provide a complete set of functionalities and also enroll interesting industrial partners to include in the platform new SDR or front end components.

B. Numbers and facts

To address all of these issues we introduce CorteXlab [1], a testbed composed of 82 heterogeneous and high performance radio nodes deployed in an electromagnetically isolated room, for cutting edge radio experimentation. With CorteXlab, our main objective is to provide a unified and comprehensive access to a large set of heterogeneous nodes in a reproducible environment, in the aim to foster experimental development of future radio techniques. CorteXlab offers 22 USRPs 2932 from National Instrument, 18 PicoSDRs from Nutaq and 42 IoT-Lab wireless sensor nodes from Hikob. It makes available the full potential of SDR, through the widely accepted GNU Radio toolkit as well as high performance real-time field programmable gate array (FPGA) development. The SDR hardware available in CorteXlab ranges from simplistic wireless sensor network (WSN) nodes to full blown 4x4 MIMO SDR nodes with agile radio capabilities. Through a carefully designed backbone network, those radio nodes are capable of cooperating to emulate complex radio technologies such as network multiple input -multiple output (MIMO) [START_REF] Venkatesan | Network mimo: Overcoming intercell interference in indoor wireless systems[END_REF], interference alignment (IA) [START_REF] Cadambe | Interference alignment and the degrees of freedom for the k user interference channel[END_REF] and physical layer based relay networks.

CorteXlab is developed, along with 8 other testbeds, under the framework of a nationwide French program Future Internet of Things (FIT) [4]. FIT aims to develop an experimental facility, a federated and competitive infrastructure with international visibility and a broad panel of users.

The inaugural meeting of this CorteXlab testbed will be held the 28 th of October 2014, and then CorteXlab will be officially open to the community for experimentations. In order to characterize all possible radio links between the different wireless nodes in this testbed, we have developed a channel sounding technique directly based on the USRP nodes. This sounder is obviously limited in quality by the intrinsic radio characteristics of these USRP, but this will enable an easy-to-use and easy-to-deploy solution for all future users in order to have an evaluation of the linking each pairs of transceivers. This channel sounder has been adapted from a previous work presented in [START_REF] Boeglen | A poor man's VANET channel sounder[END_REF] and redesign for this special purpose.

A. Experimentation room

CorteXlab counts with a completely shielded room to ensure EM isolation with the outside world. This allows the use of any (restricted or open) frequency between 300 MHz and 5 GHz, which might appeal to researchers looking to test specific frequency or standard related techniques. It also improves reproducibility since all sources of interference and channel impairments are confined to the interior of the room. To avoid excessive reflections, EM wave absorbing foams cover all walls and roof, leaving room from some reflections to occur. A glimpse of the experimentation room can be seen in Figure 1.

A measurement campaign has shown an attenuation of more than 80 dB is ensured in the band [800MHz-5GHz], with at least 60 dB throughout the whole target band. We are well aware that by shielding the room, a risk is taken on producing an unrealistic propagation environment. To properly assess the radio environment with the EM shielding, two approaches are currently under study. In the first approach, an extensive raytracing simulation campaign could provide an initial understanding of the expected kind of propagation environment, and then allow (if needed) the subsequent placement of absorbers or reflection surfaces to improve the radio characteristics of the room. The second approach, subject of the presented study, consists on running a campaign of channel measurements to properly assess the effectiveness of the changes made in the first step, as well as to provide real channel measurements to our future users. Finally, cameras will also be installed in the room to diminish human intervention to a minimum.

B. USRP software defined radio nodes

Represented by the National Instruments USRP 2932, the general-purpose SDR nodes will use (but are not limited to) the GNU Radio toolkit for rapid prototyping of transmission techniques mostly reliant on the general purpose processor (GPP) of the host PC. The USRP 2932 [6] is a high end radio platform, counting with a 400 MHz -4.4 GHz RF board, data rates of up to 40 MHz (with reduced dynamic range, nominal band of 20 MHz), a precise OCXO clock source and a 1 gigabit ethernet (GigE) connection to the host PC. The host PC is based on a Linux environment and will allow users to test not only PHY layer techniques, but also MAC, NET, TRA and APP, by deploying custom kernel modules that can coexist with the stock implementations easily.

Both the use of Linux and GNU Radio will facilitate the development and test at the user's own computer before bringing the experiment over to CorteXlab. Furthermore, the GNU Radio and Linux communities have progressed a lot in the last decade, offering good and free support to users as well as a multitude of open and accessible examples, which can be downloaded off the internet and used without additional charges.

C. Channel sounder

The GNU radio API is based on two languages. Python is in charge of setting up the signal processing flow graph forming the transceiver chain and for the synchronization between blocks. The signal processing blocks which require a high processing power are written in C++. The relation between Python and C++ is made by SWIG [7].

There already exist a lot of applications running on USRPs in several domains ranging from ZigBee to GSM. There is also an 802.11p transmitter application designed by Fuxjäger et al. [START_REF] Fuxjäger | IEEE 802.11p Transmission Using GNURadio[END_REF]. We have used this implementation as a starting point for our design. The IEEE 802.11p physical layer has similar specifications as IEEE 802.11a with some changes. In IEEE 802.11p, a 10MHz frequency bandwidth is used, instead of 20MHz bandwidth in IEEE 802.11a, thus all parameters in the time domain for IEEE 802.11p are doubled compared with IEEE 802.11a. Intrinsically, 802.11p uses OFDM.

When implementing an 802.11p transceiver, the best method we found was to use Annex G of the standard [START_REF]IEEE Standard for Information technology --Part 11: Wireless LAN Medium Access Control (MAC) and Physical Laye (PHY) Specification[END_REF]. It is a complete example of a packet encoding which gives the results of the encoding at the end of each block ending up with the baseband IQ vectors. We did the same for the corresponding blocks of the decoder. At this point, we had everything needed to simulate an 802.11p transmission [START_REF] Boeglen | A poor man's VANET channel sounder[END_REF]. In order to perform a real transmission, Automatic Gain Control (AGC) and synchronization at the receiver is needed.

But for a channel sounding application, we only need to capture the signal after it has passed through a 10MHz bandpass filter and the AGC. The synchronization and the decoding of the packets can therefore be performed off line using Matlab. Of course, one important limitation of using USRPs for this channel sounding is that the resolution is relative to the 10 MHz of bandwidth used. As a matter of example of this limitation, we can show in Fig. 3 that we can only resolve two paths if we plot the impulse response of the system in the largest dimensions of the room. This test was using a directive antenna at the emitter, creating a direct path to the receiver and one bouncing path on a reflector. Then we can achieve to separate these two paths only if we use all the length of the room.

D. Some first results

In order to promote the automated control of all nodes during an experiment, an experimentation control plane was created, called Minus. It allows for the scheduling, deployment, start, finish and results collection of user experiments.

Each experiment is organized as a Minus task, and contains:

• Scenario description: a textual file with the list of nodes, their roles and their parameters; • Scenario roles: the actual firmware, GNU Radio Python scripts, pre-compiled libraries, pre-compiled C code and FPGA bitstreams that define the node behavior during the experiment; • Role parameters: the configuration values that refine the behavior of each role.

Therefore, the GNU radio code developed to realize the channel sounding was used as a source for Minus in order to attribute roles of one emitter and 21 receivers in each experimentation. Up to now, a first run of experimentations was performed, with 6 nodes as emitters, so 126 links in total. After analysis, we could already say that, on those 126, 110 could be considered as mainly LOS, only 10 have a strong NLOS behavior and finally 6 are too weak to allow demodulation of the signal. Another interesting and simple result is the ability with this scheme to evaluate the attenuation of the signal on each link, resulting in a global value of the mean pathloss. Moreover, if we want to characterize the effect of shadowing and fading for each link, this could be more difficult with the limitation of this platform, but not necessarily meaningful for network-level users. The mean power received on each subcarrier could be sufficient information (as shown in Fig. 4) as the controlled environment is static.

But in order to give an evaluation of the diversity order offered depending on the chosen link, an interesting analysis could be based on the phase of each subcarrier. Fig. 5 shows a typical example of a NLOS scenario where the phase distribution is non-linear, resulting in higher diversity.

CONCLUSION

This paper has presented the development of a channel sounder in order to characterized all radio links in the new CorteXlab testbed. Some first results are exposed, now an extensive campaign has to be launched to fully characterized all possible links, and potentially to increase the number of metrics that are transmitted to the users.

Fig. 1 .

 1 Fig. 1. Pictures of the CorteXlab testbed, part of the FIT equipment of excellence

Fig. 3 .

 3 Fig. 3. Impulse responses of a transmission between two USRPs in the CorteXlab testbed, with only one direct path (blue) or with two paths (red) using a reflector

Fig. 4 .Fig. 5 .

 45 Fig. 4. Radio link evaluation in frequency domain between nodes 37 and 18


Hanna Becker

Ahmad Karfoul

Laurent Albera

Rémi Gribonval

Julien Fleureau

Philippe Guillotel

Amar Kachenoura

Lotfi Senhadji

Isabelle Merlet

Tensor Decomposition Exploiting Structural Constraints for Brain Source Imaging

The separation of Electroencephalography (EEG) sources is a typical application of tensor decompositions in biomedical engineering. The objective of most approaches studied in the literature consists in providing separate spatial maps and time signatures for the identified sources. However, for some applications, a precise localization of each source is required. To achieve this, a two-step approach has been proposed. The idea of this approach is to separate the sources using the canonical polyadic decomposition in the first step and to employ the results of the tensor decomposition to estimate distributed sources in the second step, using the so-called disk algorithm. In this paper, we propose to combine the tensor decomposition and the source localization in a single step. To this end, we directly impose structural constraints, which are based on a priori information on the possible source locations, on the factor matrix of spatial characteristics. The resulting optimization problem is solved using the alternating direction method of multipliers, which is incorporated in the alternating least squares tensor decomposition algorithm. Realistic simulations with epileptic EEG data confirm that the proposed single-step source localization approach outperforms the previously developed two-step approach.

I. INTRODUCTION

An important application domain of tensor decompositions in biomedical engineering is the separation of brain sources in the analysis of Electroencephalography (EEG) recordings. Over the past two decades several tensor-based approaches have been proposed to this end. These methods differ in the dimensions which are exploited in addition to space and time, e.g., frequency, wave vector, subject or realization, and in the employed decomposition model such as the Canonical Polyadic (CP) decomposition [START_REF] Möcks | Decomposing event-related potentials: a new topographic components model[END_REF], [START_REF] Miwakeichi | Decomposing EEG data into spacetime-frequency components using parallel factor analysis[END_REF], [START_REF] Morup | Parallel factor analysis as an exploratory tool for wavelet transformed event-related EEG[END_REF], [START_REF] Vos | Canonical decomposition of ictal scalp EEG and accurate source localisation: Principles and simulations study[END_REF], [START_REF] Deburchgraeve | Neonatal seizure localization using Parafac decomposition[END_REF], the PARAFAC2 model [START_REF] Weis | Multi-dimensional PARAFAC2 component analysis of multi-channel EEG data including temporal tracking[END_REF], the Shift-invariant CP (SCP) approach [START_REF] Morup | Shift-invariant multilinear decomposition of neuroimaging data[END_REF], or the block term decomposition [START_REF] Hunyadi | Block term decomposition for modelling epileptic seizures[END_REF]. While most authors contend themselves with the identification of the spatial maps and time signals of the different sources, the interpretation of the EEG data can be further improved by determining the positions and spatial extents of the sources within the brain. This is, for example, of particular importance in the pre-surgical evaluation of drug-resistant epileptic patients, where a precise delineation of the epileptogenic zone is sought. To identify the source positions, a specially developed brain source localization algorithm, referred to as the Disk Algorithm (DA) has been applied in a second step, after the tensor decomposition, to the estimated spatial maps of the sources [START_REF] Becker | Multi-way space-time-wave-vector analysis for EEG source separation[END_REF], [START_REF] Becker | EEG extended source localization: tensor-based vs. conventional methods[END_REF]. However, the results of this two-step approach depend crucially on the correct separation of all source regions in the first step [START_REF] Becker | EEG extended source localization: tensor-based vs. conventional methods[END_REF], without which the sources cannot be accurately localized. This problem arises especially in the context of close sources and sources with highly correlated time signals as involved, for example, in the propagation of epileptic phenomena.

In this paper, we propose to combine the tensor decomposition step and the source localization step in a single approach. Our idea is to incorporate a priori knowledge on the possible source positions and configurations in the tensor decomposition. The proposed method could generally be applied to any tensor decomposition approach which assumes a separate, constant vector of spatial characteristics for each source. For simplicity, in this paper, we will focus on the CP decomposition and a third order tensor with dimensions space, time, and realization. Our objective consists in imposing structural constraints based on the possible source locations within the brain on the factor matrix of spatial characteristics. Furthermore, we regularize the source localization problem following the fused LASSO strategy, which has previously been employed in the Source Imaging based on Structured Sparsity (SISSY) algorithm1 [START_REF] Becker | Fast, variation-based methods for the analysis of extended brain sources[END_REF], leading to a good performance for the localization of extended sources [START_REF] Becker | Fast, variation-based methods for the analysis of extended brain sources[END_REF]. The tensor decomposition is then computed using the Alternating Least Squares (ALS) algorithm [START_REF] Bro | Multi-way analysis in the food industry: Models, algorithms and applications[END_REF] where the classical update equation for the matrix of spatial characteristics is replaced by an Alternating Direction Method of Multipliers (ADMM) [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite elements approximations[END_REF], [START_REF] Glowinski | Sur l'approximation, par éléments finis d'ordre un, et la résolution, par pénalisation-dualité, d'une classe de problèmes de dirichlet non linéaires[END_REF] update to solve the constrained optimization problem for the matrix of spatial characteristics. As will be shown in this paper, the proposed method leads to a more accurate localization of correlated sources than the DA-based two-step algorithm. Note that a similar tensor decomposition approach exploiting (nonnegativity) constraints on the factor matrices has previously been treated in [START_REF] Liavas | Parallel algorithms for constrained tensor factorization via the alternating direction method of multipliers[END_REF], [START_REF] Sidiropoulos | Two takes on constrained CP decomposition using the alternating direction method of multipliers[END_REF], but the application context and the structural constraint considered in this paper are different.

II. DATA MODEL

EEG is a non-invasive technique that records the electrical activity of the brain with an array of N sensors placed on the scalp. Assuming that T time samples of data have been collected, the EEG recordings can be stored into the data matrix X ∈ R N ×T . In the following, we suppose that the cerebral activity can be modeled by a grid of D dipoles with positions spread all over the cortical surface and orientations perpendicular to this surface. These dipoles form the source space. The signals emitted by all dipoles of the source space are described by the matrix S ∈ R D×T . Then the EEG data constitute a superposition of these dipole signals:

X = G S. (1)
The lead field matrix G ∈ R N ×D characterizes the attenuation that is inflicted on the signal of each dipole before it can be measured by the sensors. For a given head model and source space, the lead field matrix can be computed numerically (see, e.g., [START_REF] Gramfort | Mapping, timing and tracking cortical activations with MEG and EEG: Methods and application to human vision[END_REF]) and is therefore assumed to be known in the following.

In practice, for a source signal that has sufficiently high amplitude to be measurable at the surface, a certain area of cortex needs to emit synchronized activity. A contiguous cortical area with synchronized activity, also called patch, can be modeled by a certain number of adjacent grid dipoles with highly correlated signals. Furthermore, if several patches emit highly correlated activity, they form a distributed source. Assuming that the EEG data (for one spike realization) are generated by P distributed sources and attributing the average signal s p to the p-th source, we obtain the following data model:

X = P p=1 h p s T p + X b (2)
where

h p = Gψ p (3)
corresponds to the distributed source lead field vector (in the following called spatial mixing vector). The latter is characterized by the sparse coefficient vector ψ p , whose nonzero elements describe the contributions of the associated grid dipoles to the distributed source. The signals of all dipoles that do not belong to a distributed source contribute to the measurements through the matrix of background activity X b .

The objective of brain source localization consists in identifying the coefficient vectors ψ p of all distributed sources. In this paper, this will be accomplished using tensor decomposition approaches.

III. TENSOR-BASED BRAIN SOURCE LOCALIZATION

In order to construct a tensor, we either need to apply a suitable transformation, such as the wavelet transform or the short term Fourier transform, to the data or collect an additional dimension from the measurements. Subsequently, we consider EEG data containing interictal epileptic spike signals, which repeatedly occur in irregular time intervals, involving the same source regions at each occurrence. Following [START_REF] Deburchgraeve | Neonatal seizure localization using Parafac decomposition[END_REF], [START_REF] Albera | Localization of spatially distributed brain sources after a tensor-based preprocessing of interictal epileptic eeg data[END_REF], it is then possible to build a Space-Time-Spike (STS) tensor by stacking the spike-like signals observed at R different time instants along the third dimension of the tensor X . Assuming that the signals s p of the P sources are the same for each spike realization and that the spike amplitudes of the different realizations are described by the vector c p ∈ R R , we obtain the following trilinear tensor model:

X ≈ P p=1 h p • s p • c p . ( 4
)
The spatial maps of the P distributed sources, the associated source signals, and the amplitudes can be summarized in the factor matrices H = [h 1 , . . . , h P ], S = [s 1 , . . . , s P ], and C = [c 1 , . . . , c P ].

A. Classical two-step approach

The classical approach for tensor-based distributed source localization consists of two steps. In the first step, the data tensor X is decomposed using the CP decomposition. This permits to obtain estimates of the spatial mixing vectors ĥp , p = 1, . . . , P . In the second step, these estimates can be used to localize each distributed source separately. In [START_REF] Becker | EEG extended source localization: tensor-based vs. conventional methods[END_REF], this is achieved by DA, which is based on a dictionary of M circularshaped potential distributed sources of varying sizes, the socalled disks, characterized by the coefficient vectors ψ m , m = 1, . . . , M . The spatial mixing vectors Gψ m of these disks are then compared to the estimated spatial mixing vector ĥp of the p-th source based on a metric. Finally, an estimate of the p-th distributed source is obtained by merging all disks for which the metric exceeds a given threshold value.

B. The proposed single-step approch

To perform tensor decomposition and source localization in a single step, we propose to decompose the tensor X exploiting the fact that the spatial mixing vectors h p are linear combinations of the lead field vectors as described by equation (3). To this end, we impose the structural constraint H = GΨ, where Ψ = [ψ 1 , . . . , ψ P ], on the spatial mixing matrix H. Furthermore, to obtain a sparse, piece-wise constant source distribution, which enables us to easily delineate the active patches, we adopt the fused LASSO regularization strategy of the SISSY algorithm [START_REF] Becker | Fast, variation-based methods for the analysis of extended brain sources[END_REF]. Hence we include the regularization term λ(||TΨ|| 1 + α||Ψ|| 1 ), where T is the operator implementing the variational map (see [START_REF] Ding | Reconstructing cortical current density by exploring sparseness in the transform domain[END_REF]), which describes the difference in amplitude for all pairs of adjacent grid dipoles, and λ and α are regularization parameters. In order to incorporate the structural constraint on the spatial mixing matrix and the regularization term in the tensor decomposition, we employ the ALS algorithm. More particularly, to estimate the spatial mixing matrix, we aim at solving the following optimization problem:

min H,Ψ ||X (1) -H(C S) T || 2 F + λ(||TΨ|| 1 + α||Ψ|| 1 ) s. t. H = GΨ (5)
where X (1) ∈ R N ×T R denotes the mode-1 unfolding matrix of the tensor X . This can be done using the ADMM algorithm, which introduces the latent variables Y and Z and replaces problem [START_REF] Deburchgraeve | Neonatal seizure localization using Parafac decomposition[END_REF] by

min H,Y,Z ||X (1) -H(C S) T || 2 F + λ(||Y|| 1 + α||Z|| 1 ) s. t. H = GΨ, Y = TΨ, Z = Ψ. (6)
Based on the augmented Lagrangian associated with (6), the following update equations can be derived for the spatial mixing matrix H and the coefficient matrix Ψ, the latent variables Y and Z and the Lagrangian multipliers U, V, and W:

H = (X (1) (C S) + ρGΨ + V)Θ -1 (7) Ψ = (ρT T T + ρI D + ρG T G) -1 Φ (8) Y = prox f, λ 1 ρ (TΨ + U/ρ) (9)
Z = prox f, λ 2 ρ (Ψ + W/ρ) (10) ∆U = ρ(TΨ -Y) (11) ∆V = ρ(GΨ -H) (12) ∆W = ρ(Ψ -Z). (13)
Here Θ = (C S) T (C S) + ρI P , Φ = ρ(T T Y + Z + G T (H -V)) -T T U -W, and ρ > 0 denotes a penalty parameter. Also, the prox stands for the proximity operator dealing with nonsmooth functions f (here f (M) = ||M|| 1 ) and initially introduced in [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF]. Regarding the Lagrangian multipliers U, V and W, their expressions above are given in terms of their updates using the dual ascent method where ∆M denotes the difference in the value of M between two successive iterations. The distributed source signal matrix S and the amplitude matrix C are updated using the classical ALS equations:

S = X (2) ((C H) T ) + (14) C = X (3) ((S H) T ) + . ( 15
)
Here, X (2) ∈ R T ×N R and X (3) ∈ R R×T N are the mode-2 and mode-3 unfolding matrices of X , respectively, and + denotes the Moore-Penrose pseudo-inverse.

After (random) initialization, the updates for H, S, and C using equations ( 7) to ( 15) are repeated alternatingly until convergence or a maximal number of iterations are reached. At the end of the ALS algorithm, in addition to the matrices H, S, and C, estimates of all distributed sources are directly available and correspond to the columns of the matrix Ψ.

IV. SIMULATIONS

To demonstrate the good performance of the proposed single-step tensor-based source localization approach in comparison to the classical two-step approach, we conduct a simulation study. To this end, we simulate interictal epileptic EEG data for N = 91 sensors, T = 200 time samples at a sampling rate of 256 Hz, and R = 50 epileptiform spikes. We employ a source space composed of D = 19626 grid dipoles that are located on the cortical surface. The lead field matrix G ∈ R 91×19626 is computed numerically using the ASA software (ASA, ANT, Enschede, Netherlands) and a realistic head model composed of three compartments representing the brain, the skull, and the scalp. We consider two different scenarios with three patches. The first scenario consists of a patch on the superior occipital gyrus (SupOcc) and two close patches on the superior frontal gyrus (SupFr) and the inferior frontal gyrus (InfFr). The second scenario consists of the patch InfFr and two close patches on the mid-temporal gyrus (MidTe) and the occipital temporal gyrus (OccTe). Because of the lack of space, only the second scenario is depicted in Figure 1. Each patch is composed of 100 adjacent grid dipoles, corresponding to a cortical area of about 5 cm 2 . The first two patches of each scenario are attributed epileptic spike signals of slightly different morphology that were segmented from stereotactic EEG (SEEG) recordings of a patient suffering from epilepsy. We then generate 100 different realizations of both signals for the 100 dipoles of each patch by introducing small variations in amplitude and delay. Assuming that the third patch is activated due to a propagation of the epileptic activity of the second patch, we use the same signals for the dipoles of the second and third patch, but introduce a delay of 4 to 8 ms due to the small distance between the two patches. To simulate different spike realizations, we multiply the amplitudes of the spike signals by a factor that is randomly drawn from a Gaussian distribution with mean 1 and variance 0.1. All source dipoles that do not belong to a patch are attributed Gaussian background activity with an amplitude that is adjusted to the amplitude of the SEEG signals between epileptic spikes, thus leading to a realistic Signal to Noise Ratio (SNR) of about 1.

Due to the high correlation of the signals of both the second and the third patch and to the close distance between these latter, we consider these patches to belong to the same distributed source. Therefore, we decide to separate only P = 2 CP components using the ALS algorithm for both the two-step and the single-step approaches. It is noteworthy that the same random point is used as initial guess for both approcaches. For the DA source localization in the two-step approach, we consider disks of varying sizes ranging from 1 to 100 grid dipoles. The threshold value for the metric is adjusted for each distributed source based on the goodness-of-fit criterion (cf. [START_REF] Becker | EEG extended source localization: tensor-based vs. conventional methods[END_REF]). In the single-step approach, the regularization parameter α is set to 0.67, the parameter λ is automatically selected from a range of tested values based on a heuristic criterion, and the penalty parameter ρ is set to 1.

To evaluate the source localization results obtained by the two different methods, we employ the Dipole Localization Error (DLE), which provides a measure of similarity between the original and the estimated source configuration. The DLE is defined as

DLE = 1 2Q k∈I min ∈ Î ||r k -r || 2 + 1 2 Q ∈ Î min k∈I ||r k -r || 2
where I and Î denote the original and the estimated sets of indices of all dipoles belonging to an active patch, Q and Q are the numbers of original and estimated active dipoles, and r k denotes the position of the k-th source dipole. To determine the estimated sets of active dipoles for the single-step approach, we threshold the amplitudes of the estimated source distributions.

Table I (left) lists the DLE values, averaged over 30 realizations with different signals and background activity, of the proposed single-step approach in comparison to the classical two-step method with source localization using DA. In addition, Table I (right) lists the execution time elapsed for each method considered in both scenarios using a PC of 2 GHz Quad-Core Intel Xeon with 32 GB of memory. Due to its higher localization accuracy and relatively small execution time, it is clear that the proposed single-step method gives the best localization accuracy and execution time compared to the two-step strategy.

Figure 1 shows the estimated source distributions associated with the two CP components for the tested single-step and two-step source localization algorithms for the second scenario (InfFr, MidTe, OccTe). For all algorithms, it can be seen that the patch InfFr may be associated with the first CP component whereas the patches MidTe and OccTe, that emit highly correlated signals, are combined in the second CP component. The patch InfFr is well delineated by both the single-step approach and the two-step STS-DA method. The patches MidTe and OccTe, on the other hand, are only well localized by the single-step algorithm which clearly shows two distinct source regions associated with the second CP component. STS-DA identifies a single patch in this caselocated between patches MidTe and OccTe. Note that similar behavior was observed for these methods in the case of the second scenario. V. CONCLUSION

We have proposed a new source imaging algorithm that permits to separate and localize extended EEG sources in a single step. To do so, we have employed a tensor decomposition approach which takes into account structural constraints by making use of the ALS and the ADMM optimization strategies. Simulations conducted using realistic EEG data, in the context of drug-resistent epilepsy, have shown that the proposed single-step approach provides the best localization accuracy and execution time compared to the classical twostep tensor-based source localization approach.

Fig. 1 .

 1 Fig. 1. Scenario 2. (Top) from left to right: patches InfFr and patches MidTe and OccTe with highly correlated signals. (Bottom) source localization results for both CP components obtained with the different tested algorithms. C1 stands for the first CP component while C2 denotes the second one.

TABLE I .

 I

	DLE AND EXECUTION TIME OF SOURCE IMAGING
	ALGORITHMS FOR SCENARIO 1 (OCCSUP, SUPFR, INFFR) AND SCENARIO
	2 (INFFR, MIDTE, OCCTE)
		DLE (in cm)	CPU time (in sec)
	scenario	1	2	1	2
	Single step (STS-SISSY)	1.32	1.37	75.5	78.5
	Two-step (STS-DA)	18.23 7.29	238.0	240.0

Please note that this algorithm was originally introduced under the name Sparse Variation-Based Sparse Cortical Current Distribution (SVB-SCCD) in[START_REF] Becker | Fast, variation-based methods for the analysis of extended brain sources[END_REF] and is here renamed to SISSY.

ACKNOWLEDGMENT

The work of R. Gribonval was funded by the FP7 European Research Council Programme, PLEASE project, under grant ERC-StG-2011-277906.


Thanh Phuong Nguyen
email: thanh-phuong.nguyen@ensta-paristech.fr

Antoine Manzanera
email: antoine.manzanera@ensta-paristech.fr

Noise Tolerant Descriptor for Texture Classification

Keywords: Local binary pattern, texture classification, multi-resolution, noise robustness

Among many texture descriptors, the LBP-based representation emerged as an attractive approach thanks to its low complexity and effectiveness. Many variants have been proposed to deal with several limitations of the basic approach like the small spatial support or the noise sensitivity. This paper presents a new method to construct an effective texture descriptor addressing those limitations by combining three features:

(1) a circular average filter is applied before calculating the Complemented Local Binary Pattern (CLBP), (2) the histogram of CLBPs is calculated by weighting the contribution of every local pattern according to the gradient magnitude, and (3) the image features are calculated at different scales using a pyramidal framework. An efficient calculation of the pyramid using integral images, together with a simple construction of the multi-scale histogram based on concatenation, make the proposed approach both fast and memory efficient. Experimental results on different texture classification databases show the good results of the method, and its excellent noise robustness, compared to recent LBP-based methods.

I. INTRODUCTION

Texture analysis is a very active research topic in computer vision and pattern recognition. A large number of models have been proposed [START_REF] Mirmehdi | Handbook of Texture Analysis[END_REF] for texture classification. One of the most popular approaches is based on feature distribution with Local Binary Pattern (LBP), proposed by Ojala et al. [START_REF] Ojala | A comparative study of texture measures with classification based on featured distributions[END_REF], [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF]. Thanks to its low computational cost and invariance to illumination, LBP is widely considered as an efficient descriptor for capturing local properties of images. Besides, thanks to the flexibility of the LBP framework, many variants have been proposed to fit a specific problem or to improve the robustness and discriminative power.

The two main limitations of the classic LBP codes are: (1) they capture only local features of images, and (2) they are sensitive to noise and/or uniform regions of image. Many aspects have been considered to overcome these limitations. Pre-processing steps, like Gabor filters [START_REF] Riaz | Texture classification using rotation-and scale-invariant gabor texture features[END_REF] are widely used for both capturing more spatial information and decreasing the noise. Different neighborhood topologies have been used, such as elliptical neighborhood [START_REF] Liao | Face recognition by using elongated local binary patterns with average maximum distance gradient magnitude[END_REF], three-patch or four-patch approach [START_REF] Wolf | Descriptor based methods in the wild[END_REF] to exploit anisotropic information. Multi-resolution or multi-structure approaches [START_REF] Mäenpää | Multi-scale binary patterns for texture analysis[END_REF]- [START_REF] Raja | Sparse multiscale local binary patterns[END_REF] were considered to gather information at larger scales. In encoding step, three values (-1, 0, 1) were used in Local Ternary Pattern [START_REF] Tan | Enhanced local texture feature sets for face recognition under difficult lighting conditions[END_REF] to address the sensitivity to near constant image areas. Liao [START_REF] Liao | Dominant local binary patterns for texture classification[END_REF] coded only the most frequent patterns to increase the robustness. Guo et al. used complementary component about magnitude information [START_REF] Guo | A completed modeling of local binary pattern operator for texture classification[END_REF] and also variance [START_REF] Guo | Rotation invariant texture classification using LBP variance (LBPV) with global matching[END_REF] to enhance the descriptive power and then improve the texture classification.

We propose a new approach to improve discrimination power of the classic LBP while enhacing its robustness to noise and near uniform regions. It combines spatial filtering and adaptive construction of the LBP histogram, which is weighted according to the gradient magnitude, with a multi-resolution analysis. Thanks to efficient calculation of the pyramid and construction of the multiscale histogram, it provides a fast and memory efficient texture descriptor. Its performance and noise robustness are evaluated on three texture data sets and compared with other state-of-the-art LBP-based descriptors. The next section recalls the most relevant related works. Section III introduces the proposed descriptor, and Section IV evaluates our method and compares it to state-of-the-art descriptors for the texture classification task.

II. RELATED WORKS

A. LBP and CLBP

The LBP operator [START_REF] Ojala | A comparative study of texture measures with classification based on featured distributions[END_REF], [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF] was introduced to represent local relations between a pixel and its neighbors by a binary code. The basic LBP is a P-bit code described as follows:

LBP p P,R = s(g p -g c ), s(x) = 1, x ≥ 0 0, x < 0 (1)
where g c denotes the gray value of the center pixel and g p (p = 0, ..., P -1) the gray value of the P neighbor pixels on a circle of radius R. If (x c , y c ) are the coordinate of the center pixel, then the coordinates of its neighbors are (x c + R cos(2πp/P ), y c + R sin(2πp/P )). The gray values g p are estimated by interpolation. The uniformity measure of a LBP is defined as follows:

U (LBP P,R ) = P p=1 |LBP p P,R -LBP p-1 P,R |,
Ojala et al. [START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF] then proposed the rotation invariant uniform (riu2) pattern as follows:

LBP riu2 P,R =        P -1 p=0 LBP p P,R , if U (LBP P,R ) ≤ 2 P + 1, otherwise.
In texture representation, LBP riu2 P,R are often used to reduce the feature vector length (only P + 2 distinct values instead of 2 P ), and to make the descriptor more robust to unsignificant variations.

Guo et al. [START_REF] Guo | A completed modeling of local binary pattern operator for texture classification[END_REF] significantly improved the performance of LBP riu2 P,R by considering a completed representation (CLBP). The local differences are decomposed into two complementary components: the signs CLBP S and the magnitudes CLBP M. In addition, the gray value of the center pixel is also converted into a binary code: CLBP C, using a global thresholding. CLBP codes then have 2(P + 2) 2 distinct values.

B. Noise robustness

Different methods have been proposed to address the sensivity of LBP to noise and near uniform regions. Tan and Triggs [START_REF] Tan | Enhanced local texture feature sets for face recognition under difficult lighting conditions[END_REF] proposed local ternary patterns (LTP), using ternary {-1, 0, 1} instead of binary encoding, 0 being used when the pixel difference is less than a given threshold. Ren et al. [START_REF] Ren | Noise-resistant local binary pattern with an embedded error-correction mechanism[END_REF] proposed a mechanism to recover the corrupted patterns by encoding small difference as an uncertain bit determined later using the other bits. Mean filters have been widely used in different steps like filtering [START_REF] Liao | Learning multiscale block local binary patterns for face recognition[END_REF], thresholding [START_REF] Zhao | Completed robust local binary pattern for texture classification[END_REF], or sampling [START_REF] Liu | BRINT: Binary Rotation Invariant and Noise Tolerant Texture Classification[END_REF]. Song et al. [START_REF] Song | Noise-robust texture description using local contrast patterns via global measures[END_REF] explored a set of local contrast patterns via global measures to handle noise effect.

III. MULTIRESOLUTION ADAPTIVE CLBP

This section details the proposed descriptor, whose principles are illustrated on Figures 1 and2.

A. Construction of the descriptor

The proposed Multiresolution Adaptive Complemented LBP (MACLBP) is presented in three steps: (1) Spatial filtering, (2) Weighted histogram of CLBP codes, and (3) Multiresolution histograms.

1) Spatial filtering: Mean filters are widely used to reduce the effect of additive noise. Inspired by [START_REF] Liao | Learning multiscale block local binary patterns for face recognition[END_REF], we calculate the mean filter on a spatial support made of a union of discrete circles. The size of the support should be proportionate to the neighborhood used to calculate the LBP. Aside from noise reduction, the interest of the spatial filtering step is to capture more information in a rotation invariant manner thanks to the isotropic support.

2) Contrast weighted CLBP histogram: In the descriptors based on distribution of LBP codes, every pattern classically have the same weight, whatever its significance. However many patterns due to noise or to near uniform regions should be discarded. Considering only uniform patterns is not fully satisfying, since uniform patterns may also be due to noise and conversely, a significant part of the non uniform patterns may be highly relevant. Considering only the most frequent patterns do not avoid the contribution of spurious patterns in the most frequent bins. Our proposal to address the sensitivity to noise and near uniform regions is to weight the histogram of every pattern by a significance index corresponding to the contrast measured on the filtered image, using a Sobel filter as shown on Figure 1. This can obviously be applied with any LBP variant. In practice, CLBP [START_REF] Guo | A completed modeling of local binary pattern operator for texture classification[END_REF] is chosen for its known performance in texture classification. The resulting descriptor is named Adaptive CLBP (ACLBP).

3) Pyramidal approach: To overcome the limitation of LBP to local description, a pyramidal approach is applied. The initial image is recursively reduced to half its size at each step (see Figure 2). At each resolution, ACLBP histogram is calculated. The Multiresolution Adaptive CLBP (MACLBP) descriptor is finally constructed by concatenating the ACLBP descriptors at each resolution. An important question is: should the histograms be normalized or not? Our experiments proved that the MACLBP descriptor performed better by keeping the histograms un-normalized, so that the contribution of each scale is proportional to the number of patterns, and so drops exponentially. In practice the number of resolutions is limited to 4.

B. Comments and links to other methods

1) Use of the contrast: Gradient information has been much used to design robust local descriptors. The orientation of the gradient is the main component of descriptors like HoG [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] or POEM [START_REF] Vu | Face recognition using the POEM descriptor[END_REF]. Gradient was also used in LBP-based variants to highlight edge structures. Yao and Chen [START_REF] Yao | Retrieval of translated, rotated and scaled color textures[END_REF] proposed an LBP-like computation on binary edges to describe texture. In Sobel-LBP [START_REF] Zhao | Sobel-LBP[END_REF], LBP is applied on gradient magnitude images. Weighting the LBP histogram by a contrast measure was previously done by [START_REF] Guo | Rotation invariant texture classification using LBP variance (LBPV) with global matching[END_REF] who used the gray level variance on the LBP support.

2) Multi-scale approaches: Different extensions to multiscale were proposed for LBP in the literature. Mäenpää et al. [START_REF] Mäenpää | Multi-scale binary patterns for texture analysis[END_REF] combined Gaussian low-pass filters with exponentially growing circular neighborhoods. Liao et al. [START_REF] Liao | Learning multiscale block local binary patterns for face recognition[END_REF] extended the pixel-wise classic LBP to block-wise LBP by comparing average values of blocks instead of pixel gray levels. He et al. [START_REF] He | Pyramid-based multi-structure local binary pattern for texture classification[END_REF] used a pyramidal multi-structure based on five convolution templates including anisotropic filters. The closest approach to ours is Turtinen et al. [START_REF] Turtinen | Contextual analysis of textured scene images[END_REF], who also calculated a descriptor by concatenating histograms of LBP codes at different scales, but in their case, the image of all scales are interpolated and rescaled to the initial resolution, so that the statistical contribution of the different scales are the same. 3) Computational aspects: The pyramid is efficiently computed using the integral image technique similarly as [START_REF] Viola | Robust real-time face detection[END_REF]. Once the integral image calculated, the average value of any 2 p × 2 p block is obtained by three additions and one division by a power of 2 (right bit shifting). The MACLBP descriptor is limited to a reasonable size by using the concatenation of ACLBP at the different resolution. If P is the number of neighbors considered in the basic LBP, then the number of bins of the ACLBP histogram is 2(P + 2) 2 , and so the size of the MACLBP descriptor using S different scales is 2S(P + 2) 2 .

IV. TEXTURE CLASSIFICATION EXPERIMENTS A. Datasets and protocols

The effectiveness of the proposed method is assessed by a series of experiments on three large and representative databases: Outex [START_REF] Ojala | Outex -new framework for empirical evaluation of texture analysis algorithms[END_REF], CUReT [START_REF] Dana | Reflectance and texture of real-world surfaces[END_REF] and UIUC [START_REF] Lazebnik | A sparse texture representation using local affine regions[END_REF].

The Outex dataset contains texture images from a wide variety of real materials. We consider two common test suites, TC10 and TC12, containing 24 texture classes collected under three different illuminations and nine different rotation angles. The images at angle 0 0 are chosen for training the classifier, the remaining images for testing. The CUReT (resp. UIUC) texture dataset includes 61 (resp. 25) classes with 92 (resp. 40) images in each class. The resolution of each image is 200 × 200 (resp. 640 × 480). These datasets contain images taken under significant viewpoint and illumination variations. In the training phase, 46 (resp. 23) images per class are randomly selected while the remainder is used for testing. The reported classification rate correspond to the average result over 100 random training/testing splits.

B. Classifier

For comparison purposes with other LBP-based descriptors, the 1-nearest neighbor classifier is used, using the χ 2 distance between histograms to measure the similarity between two texture images.

C. Parameter setting and component evaluation

Different neighborhood configurations have been evaluated. In the reported results, we retained the union of two discrete circles {(5, 1), (8, 2)} (as seen in Figure 1) for the spatial support of the average filter and (P, R) = (24, 3) for the support of the LBP. Four levels of resolution were used in the pyramid. We also evaluated the contribution of the three steps composing the MACLBP descriptor as described in Section III, by testing different partial versions of the descriptor on Outex TC10 without or with additive Gaussian noise. The results can be seen in Table I, showing the complementarity of the three components (Filtering, Weighting, Multiresolution, respectively denoted F, W and M in the table).

D. Comparative evaluation

Table II and Table III report the experimental results of our descriptor compared with different methods1 , on the different datasets. It can be seen that MACLBP outperforms the stateof-the-art results on the two test suites of Outex, with an improvement up to 1-3 % with respect to recent methods, and also provides very good results on CUReT and UIUC datasets, with the best average classification rate on the two datasets. [START_REF] Khellah | Texture classification using dominant neighborhood structure[END_REF] 94.52 Xu et al. [START_REF] Xu | Viewpoint invariant texture description using fractal analysis[END_REF] 92.74 BRINT2 S M (MS9) [START_REF] Liu | BRINT: Binary Rotation Invariant and Noise Tolerant Texture Classification[END_REF] 97.86 Lazebnik et al. [START_REF] Lazebnik | A sparse texture representation using local affine regions[END_REF] 72.5 96.03 [START_REF] Guo | A completed modeling of local binary pattern operator for texture classification[END_REF] 95.88 95.58 95.53 94.80 LTP [START_REF] Tan | Enhanced local texture feature sets for face recognition under difficult lighting conditions[END_REF] 88.55 87.97 87.86 86.37 LCP R=5 [START_REF] Song | Noise-robust texture description using local contrast patterns via global measures[END_REF] 94.30 92.22 90.28 85.00 ENRLBP [START_REF] Ren | Noise-resistant local binary pattern with an embedded error-correction mechanism[END_REF] 83.55 77.67 73.76 64.90 DNS [START_REF] Khellah | Texture classification using dominant neighborhood structure[END_REF] 83.03 80.00 77.80 71.60 LSEP [START_REF] He | Rotation invariant texture descriptor using local shearlet-based energy histograms[END_REF] 87.58 83.45 -72.08

Table IV evaluates the noise robustness of different methods on Outex TC10 and CUReT datasets by comparing the classification rates for different noise levels (measured using SNR i.e. Signal to Noise Ratio). It can be seen that MACLBP is much more robust than the other methods: the classification rate is nearly constant when the SNR decreases, whereas it significantly drops for the other methods.

V. CONCLUSIONS

We have proposed an efficient variant of LBP for texture representation called MACLBP. The combined use of spatial filtering, histogram weighting based on the gradient magnitude, and multi-scale approach based on dyadic pyramid allows to increase the expressiveness of the descriptor while reducing its sensitivity to noise or near uniform regions. Thanks to efficient calculation of the pyramid and concatenation of a limited number of scale specific descriptors, MACLBP has low computational and memory footprint. The experimental results show very promising results for texture classification.

Fig. 1 .

 1 Fig. 1. Spatial filtering and calculation of the weighted histogram of CLBP codes.

Fig. 2 .

 2 Fig. 2. From local to global features using the pyramidal approach.

TABLE I CONTRIBUTION

 I OF DIFFERENT COMPONENTS ON OUTEX TC10.

	Noise	F	W	M	F+W	F+M	W+M F+W+M
	Without	99.14 98.75 99.37 99.45 99.69	99.38	99.79
	SNR=5	98.96 98.07 98.78 99.35 99.43	99.30	99.74

TABLE II CLASSIFICATION

 II RATE (%) ON OUTEX DATASET.

	Method	TC10 TC12 h TC12 t Average
	MACLBP	99.79	98.98	99.35	99.37
	LNIRP/LBP/DCI PTP [30]	99.77	98.87	98.70	99.11
	LTP [12]	98.20	93.59	89.42	93.74
	DLBP+NGF [13]	98.20	91.60	87.40	92.37
	CLBP [14]	98.93	95.32	94.53	96.26
	CRLBP [17]	99.48	97.57	97.34	98.14
	LBPV [15]	97.76	95.39	95.57	96.24
	BF+CLBP [29]	99.32	96.83	96.50	97.55

We reimplemented or used existing codes for the following methods:[START_REF] Rahtu | Local phase quantization for blur-insensitive image analysis[END_REF],[START_REF] Guo | A completed modeling of local binary pattern operator for texture classification[END_REF],[START_REF] Ojala | Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[END_REF],[START_REF] Tan | Enhanced local texture feature sets for face recognition under difficult lighting conditions[END_REF],[START_REF] Vu | Improving texture categorization with biologically-inspired filtering[END_REF]. The other results are taken from the literature.


Nivine Abbas
email: nivine.abbas@telecom-paristech.fr

Thomas Bonald
email: thomas.bonald@telecom-paristech.fr

Berna Sayrac
email: berna.sayrac@orange.com

Opportunistic Gains of Mobility in Cellular Data Networks

Keywords: Cellular data networks, mobility, opportunistic scheduling, flow-level modeling, queuing theory

In this paper, we assess the performance gains of mobility on the downlink of cellular data networks. These gains are only due to the elastic nature of traffic and thus observed even under a blind, fair scheduling scheme: data are more likely transmitted when users are close to the base stations, in good radio conditions. This phenomenon is further amplified by opportunistic scheduling schemes that exploit multi-user diversity. The results are based on the analysis of flow-level traffic models and validated by system-level simulations.

I. INTRODUCTION

The efficiency of wireless data transmissions strongly depends on the variations of the radio channel quality, especially the fast variations caused by the constructive and destructive combinations of the multiple signal paths between the transmitter and the receiver. While traditional wireless systems use diversity techniques to mitigate these channel variations and thus provide constant bit rates, advanced wireless systems exploit these variations through so-called opportunistic schedulers, taking advantage of the inherent "elasticity" of data transfers to increase the overall system capacity [START_REF] Tse | Multiuser Diversity in Wireless Networks[END_REF], [START_REF] Ayesta | A modeling framework for optimizing the flow-level scheduling with time-varying channels[END_REF], [START_REF] Bonald | A score-based opportunistic scheduler for fading radio channels[END_REF].

In this paper, we show that slow channel variations, due to the mobility of the device, can be exploited as well. In fact, they are already exploited by current packet schedulers, even blind strategies like the round-robin policy. The reason is that, due to mobility, data are more likely transmitted when users are close to the base stations, in good radio conditions. This phenomenon is further amplified by opportunistic schedulers.

The fact that mobility may improve performance has already been observed in the literature, see for instance [START_REF] Bonald | How mobility impacts the flow-level performance of wireless data systems[END_REF], [START_REF] Bonald | Flow-level performance and capacity of wireless networks with user mobility[END_REF], [START_REF] Borst | Capacity of wireless data networks with intra-and inter-cell mobility[END_REF], [START_REF] Borst | Mobility-driven scheduling in wireless networks[END_REF], [START_REF] Karray | User's mobility effect on the performance of wireless cellular networks serving elastic traffic[END_REF], [START_REF] Kim | Stability of flow-level scheduling with Markovian time-varying channels[END_REF]. However, these papers mainly give qualitative results on the impact of mobility. In the present paper, we assess the intra-cell and inter-cell mobility-induced gain under a fair scheduling scheme (the round-robin policy) and an opportunistic scheduling scheme (the max C/I policy, see [START_REF] Kolding | High speed downlink packet access: WCDMA evolution[END_REF]) for both the cell center and the cell edge. We show in particular that under both scheduling policies, mobility improves throughput performance at the cell edge. This is a rather surprising result for the max C/I policy, which gives priority to cell-center users. The reason is that these users quickly complete their file transfers, which saves scheduling resources for cell-edge users. The results are based on the analysis of flow-level traffic models and validated by systemlevel simulations.

The rest of the paper is organized as follows. We present the reference model in the absence of mobility in the next section. We analyse the impacts of intra-cell and inter-cell mobility in Sections III, IV and V. Section VI investigates the impact of fast fading. Section VII presents the system-level simulations and Section VIII concludes the paper.

II. REFERENCE MODEL

We first present the reference model in the absence of mobility and the corresponding performance results under round-robin and max C/I scheduling schemes.

A. Cellular network

We model each cell by a set of N regions. In each region, radio conditions are supposed to be homogeneous and thus users are served at the same physical data rate on the downlink. In the simple case of two regions illustrated by Figure 1, users may be close to the cell center and experience good radio conditions (light gray) or close to the cell edge and suffer from bad radio conditions (dark gray). We model each region by a queue with a specific service rate corresponding to the physical data rate in this region; since all users in the cell share the same radio resources, each cell can be viewed as a set of N parallel queues with coupled processors. The precise coupling depends on the scheduling policy, as explained below.

B. Traffic model

We consider elastic traffic only, corresponding to data transfers. Specifically, we assume that new data flows are generated in region i at the random times of a Poisson process of intensity λ i . We denote by λ = i λ i the total flow arrival rate in the cell and by p i = λ i /λ the probability that a new data flow is generated by a user in region i.

For the sake of tractability, we use a flow-level model where each data flow is viewed as a fluid of random volume to be transmitted. The volumes have an exponential distribution with mean σ (in bits). When region i is served, flows are completed at rate µ i in the absence of fast fading, corresponding to the physical rate µ i σ (in bit/s). We assume regions are numbered in decreasing order of physical rates, that is

µ 1 > µ 2 > . . . > µ N .
Queue i represents the number of active flows in region i. The load of queue i is the ratio ρ i = λ i /µ i of arrival rate to service rate. Since the radio resources are shared in time, the actual service rate of queue i is modulated by φ i , the fraction of time spent by the scheduler on users in region i. This depends on the system state and the scheduling policy. For work-conserving policies, we have i φ i = 1. The overall stability condition is then ρ < 1 in the absence of mobility, with ρ = i ρ i . We shall see in Sections III and IV that mobility may increase the stability region. Observe that the cell load may be written ρ = λ/µ where µ is the harmonic mean service rate weighted by the distribution of arrivals in the N regions:

µ = 1 i p i /µ i .

C. Throughput metrics

We measure performance in terms of mean throughputs in the different regions. We denote by X(t) the system state at time t. This is an N -dimensional vector whose component i gives the length of queue i at time t. It is an irreducible Markov process with stationary distribution π. In any state x such that x i > 0, each user in region i has throughput µ i σφ i (x)/x i . Now the distribution seen by users in region i is the sizebiased distribution [START_REF] Bonald | Network performance analysis[END_REF]:

π i (x) ∝ x i π(x).
Observe in particular that π i (x) = 0 in all states x such that x i = 0. We denote by E i the corresponding expectation. The mean throughput of users in region i is then given by

γ i = E i µ i σφ i (X) X i = E(µ i σφ i (X)) E(X i ) . (1)
By the traffic conservation equation

λ i = E(µ i φ i (X)), (2)
we deduce

γ i = λ i σ E(X i ) . (3)
Observe that, by Little's law, this is the ratio of mean flow size (in bits) to mean flow duration in region i.

By a similar argument, the mean throughput in the cell is

γ = E( i µ i σφ i (X)) E( i X i ) . (4)
By the traffic conservation equation

λ = i E(µ i φ i (X)), (5)
we deduce

γ = λσ i E(X i ) . (6)
We shall see that, in the absence of mobility, the mean throughput in region i decreases from the physical data rate µ i σ to 0 when the cell load ρ grows from 0 to 1, while the mean throughput in the cell decreases from the harmonic mean physical data rate µσ to 0, the precise impact of the cell load depending on the scheduling policy.

D. Round-robin policy

Under the round-robin policy, users share the radio resources equally, independently of their radio conditions. Thus users in region i are allocated a fraction

φ i (x) =
x i j x j of radio resources in state x. The transition graph of the Markov process X(t) is shown in Figure 4 (omitting the blue arrows, corresponding to mobility). This Markov process is reversible and corresponds to the state of a multi-class processor-sharing queue. The stationary distribution of X(t) is:

π(x) = (1 -ρ) i x i x 1 , . . . , x N ρ x1 1 . . . ρ x N N .
We deduce the mean number of flows in region i:

E(X i ) = ρ i 1 -ρ ,
and, in view of (3), the mean throughput in region i:

γ i = µ i σ(1 -ρ). (7)
Observe that the mean throughputs decrease linearly with the cell load in all regions.

In view of ( 6), the mean throughput in the cell is

γ = µσ(1 -ρ). (8)
It also decreases linearly with the cell load.

E. Max C/I policy

The max C/I policy is an opportunistic scheduling strategy that prioritizes those users with the best radio conditions. For two regions for instance, cell-center users are scheduled first and are allocated all the resources whenever active; cell-edge users are served only when there are no active cell-center users. The transition graph of the Markov process X(t) is shown in Figure 5 (omitting the blue arrows). The Markov process is no longer reversible and corresponds to a preemptive priority queue. Applying known results of queuing theory [14], we get the mean number of flows in region i:

E(X i ) = ρ i 1 -j<i ρ j   1 + ρ 1 -j<i ρ j µ i µ j<i p j µ 2 j  
and, in view of (3), the mean throughput in region i:

γ i = µ i σ 1 -j<i ρ j 1 + ρ 1- j≤i ρj µ i µ j<i pj µ 2 j .
Observe that the mean throughput in region i is positive whenever j≤i ρ j < 1: it decreases from µ i σ to 0 when the load generated by users in regions 1, . . . , i grows from 0 to 1.

Similarly, we get the mean throughput in the cell

γ = µσ   i ρ i /ρ 1 -j<i ρ j   1 + ρ 1 -j≤i ρ j µ i µ j<i p j µ 2 j     -1
which decreases from µσ to 0 when the cell load ρ grows from 0 to 1. This is the harmonic mean of the flow throughputs in the N regions.

F. Numerical results

Figure 2 compares the throughput performance of both policies for two regions, with p 1 = p 2 = 1/2, µ 1 = 8, µ 2 = 2 and σ = 1. Observe that, under max C/I scheduling, the cell-center throughput is positive whenever ρ 1 < 1, that is ρ < 5, and decreases approximately linearly. The throughput gain compared to the round-robin policy is huge for cell-center users, at the expense of a slight throughput degradation for cell-edge users. The cell-average throughput remains approximately the same since it corresponds to the harmonic mean throughput over the cell and thus is mainly affected by cell-edge users.

III. IMPACT OF INTRA-CELL MOBILITY

In this section, we add intra-cell mobility to the previous model, as shown in Figure 3 for two regions. We suppose that users move from the center of the cell to the edge and vice versa. We still assume that there is no fast fading.

A. Mobility model

We assume that each user in region i moves to region i -1 (for i > 1) and to region i + 1 (for i < N ) after exponential durations, at respective rates ν i,i+1 and ν i,i-1 . The probability that a user is in region i then satisfies:

q i ∝ i-1 j=1 ν j,j+1 ν j+1,j
.

Note that this is not the probability that an active user is in region i, which is given by

p i = E(X i )
j E(X j ) and depends both on the cell load and on the scheduling policy. In state x, the total mobility rate from region i to region i + 1 and from region i + 1 to region i (for i < N ) are x i ν i,i+1 and x i+1 ν i+1,i , respectively. For notational convenience, we write ν 1 ≡ ν 1,2 and ν 2 ≡ ν 2,1 for two regions. The corresponding state diagrams are given by Figures 4 and5, where the blue arrows correspond to the transitions due to mobility.

It is worth noting that, for mobile users, the traffic conservation equation (2) no longer applies: the traffic arriving in region i is not equal in general to the traffic served in region i. Thus the mean throughput in region i is given by (1) but not by [START_REF] Bonald | A score-based opportunistic scheduler for fading radio channels[END_REF]. Now the overall traffic conservation equation [START_REF] Bonald | Flow-level performance and capacity of wireless networks with user mobility[END_REF] still applies (in the absence of inter-cell mobility) so that the mean throughput in the cell is given by [START_REF] Bonald | Network performance analysis[END_REF]. It follows from ( 1) and ( 4) that γ = i p i γ i .

B. Round-robin policy

Under the round-robin policy, there is no longer an explicit expression for the stationary distribution of the Markov process X(t), except in the limiting regime of infinite mobility where ν i,i+1 , ν i+1,i → ∞ for all i < N . In this case, the N queues are equivalent to one processor-sharing queue with service rate:

μ = i q i µ i . 0, 0 1, 0 0, 1 1, 1 2, 0 0, 2 2, 1 1, 2 2, 2 λ 1 λ 1 µ 1 µ 1 λ 2 λ 2 µ 2 µ 2 2 3 µ 2 1 2 µ 2 1 2 µ 2 1 3 µ 2 λ 1 λ 1 2 3 µ 1 1 2 µ 1 λ 1 λ 1 1 2 µ 1 1 3 µ 1 λ 2 λ 2 λ 2 λ 2 ν1 2ν1 ν1 2ν1 ν2 ν2 2ν2 2ν2
Fig. 4. State diagram for two regions under the round-robin policy.

We denote by ρ = λ/μ the corresponding load. The stationary distribution of X(t) becomes:

π(x) = (1 -ρ)ρ n n x 1 , . . . , x N q x1 1 . . . q x N N ,
with n = i x i , under the stability condition ρ < 1. The flow throughput is independent of the region where the flow starts and is given by: γ = μσ (1 -ρ) .

The stability condition, which can be written

ρ < μ µ , (9)
is actually independent of the mobility rates: at maximum load, flows tend to be infinitely long so that users move in the whole cell during their file transfers and have the mean physical rate μσ. In the particular case where p i = q i for all i (that is, the flow arrivals are distributed according to the distribution of the user positions), then μ ≥ µ (the weighted arithmetic mean is larger than the corresponding weighted harmonic mean) and the maximum load increases, as observed in [START_REF] Bonald | How mobility impacts the flow-level performance of wireless data systems[END_REF]. Under the max C/I policy, it is also possible to derive the stationary distribution of the Markov process X(t) in the limiting regime of infinite mobility. Given n active users, the probability that the user(s) with the best radio conditions is (are) in region i is given by

C. Max C/I policy 0, 0 1, 0 0, 1 1, 1 2, 0 0, 2 2, 1 1, 2 2, 2 λ 1 λ 1 µ 1 µ 1 λ 2 λ 2 µ 2 µ 2 λ 1 λ 1 µ 1 µ 1 λ 1 λ 1 µ 1 µ 1 λ 2 λ 2 λ 2 λ 2 ν1 2ν1 ν1 2ν1 ν2 ν2 2ν2
q i (n) =   j≥i q i   n -   j>i q i   n .
Since the max C/I policy only serves this (these) user(s), the mean service rate becomes:

μ(n) = i q i (n)µ i .
Observe that μ(n) is increasing, with μ(1) = μ and lim

n→∞ μ(n) = µ 1 .
The stationary distribution of X(t) is:

π(x) = π(0) n k=1 λ μ(k) × n x 1 , . . . , x N q x1 1 . . . q x N N ,
with n = i x i , under the stability condition λ < µ 1 .

This stability condition is in fact independent of the mobility rates. The maximum load is:

ρ < µ 1 µ , (10)
which is always larger than in the absence of mobility.

D. Throughput in light traffic

The performance in light traffic (that is, when ρ → 0) is independent of the scheduling policy. Indeed, a user when alone in the system is always allocated all radio resources. The mean throughput in region i is then equal to µ i σ. For the mean throughput in the cell, we consider the simple case of two regions. The probability that a user in region 1 moves to region 2 before leaving the system is:

α 1 = ν 1 ν 1 + µ 1 .
Similarly, the probability that a user in region 2 moves to region 1 before leaving the system is:

α 2 = ν 2 ν 2 + µ 2 .
We deduce the mean duration of a flow initiated in region 1

1 ν 1 + µ 1 + α 1 1 ν 2 + µ 2 + α 1 α 2 1 ν 1 + µ 1 + . . . = n≥0 (α 1 α 2 ) n 1 ν 1 + µ 1 + α 1 n≥0 (α 1 α 2 ) n 1 ν 2 + µ 2 = 1 1 -α 1 α 2 1 ν 1 + µ 1 + α 1 ν 2 + µ 2 .
Similarly, the mean duration of a flow initiated in region 2 is

1 1 -α 1 α 2 1 ν 2 + µ 2 + α 2 ν 1 + µ 1 .
By ( 6), the cell-average throughput in light traffic is

γ = ν 2 µ 1 + ν 1 µ 2 + µ 1 µ 2 ν 1 + ν 2 + p 1 µ 2 + p 2 µ 1 σ. (11)
Observe that this throughput depends on user's mobility: γ = µσ when ν 1 , ν 2 → 0 (no mobility) while γ → μσ when ν 1 , ν 2 → +∞ (infinite mobility).

E. Numerical results

Consider two regions with p 1 = p 2 = 1/2, µ 1 = 8, µ 2 = 2 and σ = 1. The mobility rates are symmetric, that is ν 1 = ν 2 = ν. We consider the static case where ν = 0, a case with moderate mobility ν = 0.1, and the case of very high (nearly infinite) mobility, ν = 100. The results are obtained by the numerical evaluation of the stationary distribution of the Markov process X(t) and shown in Figure 6. As expected, the maximum load when users move is independent of the mobility rate ν; it is approximately equal to 1.6 under the round-robin policy, in accordance with [START_REF] Elayoubi | Mobility-aware admission control scheme in the downlink of third-generation wireless systems[END_REF], and equal to 2.5 under the max C/I policy, in accordance with [START_REF] Grossglauser | Mobility increases the capacity of ad-hoc wireless networks[END_REF], both corresponding to huge gains compared to the static case. The cell-average throughput in light traffic is independent of the scheduling policy but very sensitive to the mobility rate, as predicted by [START_REF] Karray | User's mobility effect on the performance of wireless cellular networks serving elastic traffic[END_REF].

We observe that mobility improves the cell-average throughput for both scheduling policies, with a higher gain under the max C/I policy. However, mobility may decrease the cell-center throughput (depending on the load) since users may then suffer from bad radio conditions. Surprisingly, the max C/I policy outperforms the round-robin policy for cell-edge users in case of mobility, even in case of moderate mobility ν = 0.1: the fact that cell-edge users are not scheduled (in the presence of center-cell users) is compensated by the fact that cell-center users complete more rapidly their file transfers.

IV. MULTI-CLASS EXTENSION

We have assumed so far that all users have the same mobility behavior. In particular, they are all static or all mobile. We here extend the results to multiple classes of mobility, where each user may be static or mobile, as illustrated by Figure 7.

A. Mobility model

We consider K classes of mobility. For instance, each class represents a range of speeds (pedestrian, train, etc). Class-1 users are static. For any k > 1, a class-k user moves from region i to region i + 1 and from region i + 1 to region i (for i < N ) at respective rates ν k i,i+1 and ν k i+1,i . We denote by w k the probability that a new flow is initiated by a class-k user. Let X k i (t) be the number of active class-k flows in region i at time t. The Markov process X(t) now has dimension N × K. The probability that a class-k user is in region i satisfies

q k i ∝ i-1 j=1 ν k j,j+1 ν k j+1,j
.

B. Round-robin policy

In the limiting regime where mobile users have infinite speeds, class-k users are served at rate

μk = i q k i µ i ,
for any k > 1, while class-1 users are static and thus served at the mean rate µ. Let ρ1 = w 1 λ/µ and ρk = w k λ/μ k for k > 1. The queuing system corresponds to a multi-class processor-sharing queue of load ρ = k ρk . We deduce the stationary distribution of the Markov process X(t):

π(x) = (1 -ρ)ρ n × i x 1 i x 1 1 ,...,x 1 N p1µ1 µ x 1 1 . . . p N µ N µ x 1 N × k>1 i x k i x k 1 , . . . , x k N q k 1 x k 1 . . . q k N x k N ,
under the stability condition ρ < 1, that is

ρ < 1 w 1 + k>1 w k µ/μ k . (12)
The mean throughput of static users in region i is

γ 1 i = µ i σ(1 -ρ),
while the mean throughput of mobile users depends on their class k > 1:

γ k = μk σ(1 -ρ).
The cell-average throughput is

γ = μσ(1 -ρ),
where

μ = 1 w 1 /µ + k>1 w k /μ k .

C. Max C/I policy

These results can hardly be extended to the max C/I policy. Only the stability condition is explicit. At maximum load, mobile users are served only in region 1; we deduce the stability condition

λ(w 1 p 1 + 1 -w 1 ) < µ 1 , that is ρ < µ 1 µ(w 1 p 1 + 1 -w 1 ) . (13)
For static users in region i, the stability condition is

ρ < 1 w 1 j≤i pj µ µj + (1 -w 1 ) µ µ1 , (14)
which is more restrictive.

D. Throughput in light traffic

When ρ → 0, the mean throughput of static users in region i is µ i σ while the mean throughput of mobile users of class k > 1 is, for two regions,

γ k = ν k 2 µ 1 + ν k 1 µ 2 + µ 1 µ 2 ν k 1 + ν k 2 + p 1 µ 2 + p 2 µ 1
σ.

E. Numerical Results

Figure 8 gives the results for the same example as in §III-E but with two classes of users: a static class and a mobile class, with w 1 = w 2 = 50%. The maximum load is approximately 1.2 for the round-robin policy, as given by ( 12), and 1.4 under the max C/I policy, as given by ( 13), while static cell-center users have positive throughput up to load 3.2, in accordance with (14). We observe that, even for static users at the cell edge, the max C/I policy outperforms the round-robin in the presence of mobile users.

V. IMPACT OF INTER-CELL MOBILITY

We now take hand-overs into account: users may leave the cell during the file transfer. Similarly, some users may arrive from neighboring cells. We assume that handovers occur only in region N (the cell edge) and that the incoming and outgoing handover rates are equal.

Let ν h be the outgoing handover rate: each user in region N tends to leave the cell at rate ν h . To balance incoming and outgoing handovers, we assume that incoming handovers are generated according to a Poisson process of intensity E(X N )ν h . Note that the net flow arrival rate in region N is still equal to λ N . In particular, the traffic conservation equation [START_REF] Bonald | Flow-level performance and capacity of wireless networks with user mobility[END_REF] applies and the stability condition remains the same as without handovers for both policies.

These properties suggest that inter-cell mobility has a limited impact on throughput performance, in the presence of intra-cell mobility. This is confirmed by the results of Figure 9, obtained for the example of §III-E with a handover rate ν h = ν = 0.1, using fixed-point iterations to estimate the incoming handover rate.

VI. IMPACT OF FAST FADING

We have so far ignored the impact of fast fading, which is the essence opportunistic schedulers [START_REF] Tse | Multiuser Diversity in Wireless Networks[END_REF]. Fast fading occurs on a time scale of the order of ms, much shorter than the typical time scale of slow fading, of the order of second. This is why these two phenomena can be considered separately.

Consider n active users in region i, without any active users in regions 1, . . . , i -1 so that only these n users are scheduled. Assuming Rayleigh fading, users' SINR (signalto-interference-plus-noise ratios) S 1 , S 2 ...S n are exponentially distributed. The resulting SINR gain is given by:

F (n) = E [max(S 1 , S 2 , . . . , S n )] = 1 + 1 2 + ... + 1 n .
The throughput gain is typically lower due to the concavity of the data rate with respect to SINR. Assuming that the physical data rate is proportional to log(1 + S) for SINR S, we deduce the rate gain of opportunistic scheduling:

G(n) = E[ log 2 (1 + max(S 1 , S 2 , . . . , S n )) ] E[ log 2 (1 + S 1 ) ] .
The corresponding gains are shown in Figure 10.

In any state x such that x i > 0, the total throughput of users in region i becomes:

µ i σφ i (x)G(x i ).
The results obtained for the same scenario as in §III-E (without inter-cell mobility) are shown in Figure 11. Compared to Figure 6, the max C/I policy brings a significant throughput gain due to fast fading, in addition to that due to slow fading.

VII. SYSTEM SIMULATIONS

A. Simulation setting

We now validate the results of previous sections by systemlevel simulations based on the LTE technology. We consider the 21 hexagonal cells formed by 7 tri-sector sites (a reference site surrounded by 6 interferring sites). Both eNodeB and mobile devices are equipped with 2 antennas. Spatially correlated shadowing maps for different sites are generated each 3s, with time correlation. Each TTI (of 1 ms), 4 spatially correlated fading complex values are generated for each eNodeB to mobile device pair. The main parameters are summarized in Table I Traffic consists of file transfers only. Flows arrive according to a Poisson process with uniform spatial distribution. File sizes are generated from an exponential distribution with mean σ = 0.5 MB. Scheduling decisions are taken in each eNodeB at each TTI. In order to evaluate the proposed model, the two previous schedulers are implemented. The round-robin scheduler serves users in a cycle way, independently of their radio conditions; the max C/I scheduler chooses at each TTI the user with the highest instantaneous SINR.

We study two scenarios, with and without mobility. With mobility, each user chooses a direction at random and moves along this direction at a constant speed of 30, 60, 100 and 300 km/h until the file transfer is completed. The network is a torus in this case to avoid border effects (that is, a user leaving cell 20 from the south arrives in cell 6, 7 or 13 to the north according to the direction). We simulate 10 minutes and estimate for each scenario the cell-average throughput as the ratio of mean flow size to mean flow duration.

B. Results

Figure 12 shows the normalized cell-average throughput (in bit/s/Hz) as a function of the arrival rate per eNodeB in each scenario. Results are very close to those obtained by analysis.

VIII. CONCLUSION

We have quantified the gains induced by mobility in cellular data networks. Both the round-robin and max C/I policies improve throughput performance, with a proportionally higher gain for the latter. By prioritizing cell-center users, the max C/I gives the opportunity to mobile cell-edge users to be served in better radio conditions. Although the focus of this paper is on the downlink, we think that similar performance gains can be obtained on the uplink. Other interesting problems for future research include opportunistic scheduling policies that prioritize short flows, as considered in [START_REF] Aalto | On the optimal trade-off between SRPT and opportunistic scheduling[END_REF] for instance, as well as heterogeneous networks consisting of both picocells and macro cells [START_REF] Kavitha | Spatial queueing for analysis, design and dimensioning of Picocell networks with mobile users[END_REF].

Fig. 1 .

 1 Fig. 1. A simple model with two cell regions.

Fig. 2 .

 2 Fig. 2. Throughput performance for two regions without mobility.

Fig. 3 .

 3 Fig. 3. Queuing model for two regions with intra-cell mobility.

2ν2 Fig. 5 .

 2ν25 Fig. 5. State diagram for two regions under the max C/I policy.

Fig. 6 .

 6 Fig. 6. Throughput performance for two regions with and without mobility, under round-robin (black) and max C/I (blue).

Fig. 7 .

 7 Fig. 7. Queuing model for two regions and two classes of mobility.

Fig. 8 .

 8 Fig. 8. Throughput performance for two regions with two classes of mobility, under round-robin (black) and max C/I (blue).

Fig. 9 .

 9 Fig.9. Throughput performance for two regions with both intra-cell and inter-cell mobility, under round-robin (black) and max C/I (blue).

Fig. 10 .Fig. 11 .

 1011 Fig. 10. SINR gain F (n) and rate gain G(n) for n scheduled users.

Fig. 12 .

 12 Fig. 12. Throughput performance obtained by system-level simulation.

  .

	Network topology	Macro cells only
	Environment	Urban
	Context		outdoor
	Inter-site distance	500 m
	Resolution		5 m
	PathLoss		ITU Model
	shadowing std	4 dB
	shadowing cross correlation (between sites)	0.5
	shadowing correlation distance	30 m
	Radio access technology	LTE
	Number of tx/rx antennas	MIMO (2x2)
	number of streams	1
	Receiver		MRC
	Codebook		3GPP
	carrier frequency	2 GHz
	Bandwidth		10 MHz
	TABLE I.	SIMULATION SETTING.