Juan F P J Abascal

Manuel Desco

Keywords: Compressed sensing, lung diffusion MRI, hyperpolarized gas MRI, sparsity, split Bregman method

Diffusion MRI data are generally acquired using hyperpolarized gases during patient breath-hold, which yields a compromise between achievable image resolution, lung coverage and number of b-values. In this study, we propose a novel method that accelerates the acquisition of diffusion MRI data by undersampling in both the spatial and b-value dimensions and incorporating knowledge about signal decay into the reconstruction (SIDER). SIDER is compared with total variation (TV) reconstruction by assessing its effect on both the recovery of ventilation images and the estimated mean alveolar dimensions (MAD). Both methods are assessed by retrospectively undersampling diffusion datasets of healthy volunteers and COPD patients (n=8) for acceleration factors between x2 and x10. TV led to large errors and artifacts for acceleration factors equal to or larger than x5. SIDER improved TV, with a lower solution error and MAD histograms closer to those obtained from fully sampled data for acceleration factors up to x10. SIDER preserved image quality at all acceleration factors, although images were slightly smoothed and some details were lost at x10. In conclusion, we developed and validated a novel compressed sensing method for lung MRI imaging and achieved high acceleration factors, which can be used to increase the amount of data acquired during breath-hold. This methodology is expected to improve the accuracy of estimated lung microstructure dimensions and provide more options in the study of lung diseases with MRI.

I. INTRODUCTION

HE increase in MRI sensitivity using hyperpolarized contrast media has enabled the development of imaging techniques to assess anatomical features and functional processes beyond the limits of conventional MRI [START_REF] Salerno | Hyperpolarized noble gas MR imaging of the lung: potential clinical applications[END_REF], [START_REF] Roos | Hyperpolarized Gas MR Imaging[END_REF]. In particular, hyperpolarized noble gas MRI can provide quantitative maps of clinically relevant anatomical and physiological parameters (e.g. ventilation distribution [START_REF] Deninger | Quantitative measurement of regional lung ventilation using 3 He MRI[END_REF], [START_REF] Horn | Quantification of regional fractional ventilation in human subjects by measurement of hyperpolarized 3He washout with 2D and 3D MRI[END_REF], acinar airway dimensions [START_REF] Wu | Accelerated MR diffusion tensor imaging using distributed t sensing[END_REF], alveolar oxygen partial pressure (pO2) [START_REF] Deninger | Quantification of regional intrapulmonary oxygen partial pressure evolution during apnea by (3)He MRI[END_REF], [START_REF] Marshall | 3)He pO2 mapping is limited by delayed-ventilation and diffusion in chronic obstructive pulmonary disease[END_REF], and gas washout [START_REF] Deppe | Combined measurement of pulmonary inert gas washout and regional ventilation heterogeneity by MR of a single dose of hyperpolarized 3He[END_REF]). An increase in the MRI signal is achieved using laser polarization techniques (e.g. optical pumping-spin exchange [START_REF] Walker | Spin-exchange optical pumping of noble-gas nuclei[END_REF], [START_REF] Parnell | Enhancement of X129e polarization by off-resonant spin exchange optical pumping[END_REF]), which result in a nonequilibrium net magnetization that is up to 5 orders of magnitude higher than in conventional MRI (thermally polarized).

However, the non-renewable nature of the magnetization in hyperpolarized gases imposes limits on the duration of MRI acquisition. In the presence of oxygen, the T1 of hyperpolarized gases (15-20 s) is of the order of duration of the breath-hold that can be achieved by patients within the scanner. As a consequence of this limitation, together with the high cost of the gas (e.g. 4 He and isotopically enriched 129 Xe), most hyperpolarized gas methods aim to perform a complete acquisition during a single breath-hold using a single hyperpolarized gas dose. Furthermore, there are advantages in acquisition with more than one sequence or even different nuclei during a single breath-hold [11]. This need for rapid acquisition has been addressed using various accelerated J.F.P.J. Abascal was with the Departamento de Bioingeniería e Ingeniería Aeroespacial, Universidad Carlos III de Madrid (UC3M), 28903, Leganes, Spain and with Instituto de Investigación Sanitaria Gregorio Marañón (IiSGM), 28007 Madrid, Spain, and is now with Université de Lyon, INSA-Lyon, Université Claude Bernard Lyon 1, UJM-Saint Etienne, CNRS, Inserm, CREATIS UMR 5220, U1206, F-69100, Lyon, France (email: juan.abascal@creatis.insa-lyon.fr).

M. Desco and J. Parra-Robles are with Departamento de Bioingeniería e Ingeniería Aeroespacial, Universidad Carlos III de Madrid (UC3M), 28903 Leganes, Spain and with Instituto de Investigación Sanitaria Gregorio Marañón (IiSGM), 28007 Madrid, Spain (email: desco@hggm.es), (email: juparrar@ing.uc3m.es).

M. Desco is also with Centro de Investigaciones Cardiovasculares Carlos III (CNIC) and with Centro de Investigación en Red de Salud Mental (CIBERSAM), 28029 Madrid, Spain.

Incorporation of prior knowledge of signal behavior into the reconstruction to accelerate the acquisition of diffusion MRI data acquisition approaches, including parallel imaging [START_REF] Lee | Advantages of parallel imaging in conjunction with hyperpolarized helium--a new approach to MRI of the lung[END_REF], [START_REF] Chang | In vivo lung morphometry with accelerated hyperpolarized 3 He diffusion MRI: A preliminary study[END_REF] and compressed sensing (CS) [START_REF] Ajraoui | Compressed sensing in hyperpolarized 3He lung MRI[END_REF], [START_REF] Ajraoui | Incorporation of prior knowledge in compressed sensing for faster acquisition of hyperpolarized gas images[END_REF]. In CS, the acquisition was accelerated in the spatial encoding direction and the images were reconstructed using spatial total variation (TV). More recently, Chan et al [START_REF] Chan | Whole lung morphometry with 3D multiple b-value hyperpolarized gas MRI and compressed sensing[END_REF] used CS to acquire diffusion images of hyperpolarized gases in the lungs. Diffusion MRI images are sensitive to changes in lung microstructure [START_REF] Swift | Emphysematous changes and normal variation in smokers and COPD patients using diffusion 3He MRI[END_REF] due to disease and are used to estimate the dimensions of acinar airways based on theoretical models obtained from numerical simulations. These theoretical models require the acquisition of images for several diffusion sensitization values (b-values), which, together with the longer duration of the diffusion scan (due to the presence of diffusion gradients), result in long acquisition times. Owing to the limitation imposed by breath-hold duration, there is a compromise between achievable image resolution, number of slices and number of b-values, which limits the accuracy and number of parameters of the theoretical models [START_REF] Parra-Robles | The influence of lung airways branching structure and diffusion time on measurements and models of short-range 3He gas MR diffusion[END_REF]. Typically, most protocols acquire 5 slices (10 mm thick, spacing 10 mm) with a resolution of 64x64 pixels and 4 to 6 bvalues, thus leading lung coverage to be sacrificed [START_REF] Parra-Robles | The influence of field strength on the apparent diffusion coefficient of 3He gas in human lungs[END_REF]. Chan et al. achieved full lung coverage by using CS with a 3D diffusion acquisition and undersampling along the two phaseencoded spatial directions [START_REF] Chan | Whole lung morphometry with 3D multiple b-value hyperpolarized gas MRI and compressed sensing[END_REF]. CS has also been proposed for accelerating diffusion-weighted imaging [START_REF] Wu | Accelerated MR diffusion tensor imaging using distributed t sensing[END_REF], [START_REF] Michailovich | Spatially regularized compressed sensing for high angular resolution diffusion imaging[END_REF]- [START_REF] Ning | A joint compressed-sensing and super-resolution approach for very high-resolution diffusion imaging[END_REF]. In diffusion spectrum imaging, CS has been applied using adaptive dictionaries [START_REF] Bilgic | Accelerated diffusion spectrum imaging with compressed sensing using adaptive dictionaries[END_REF] and different types of undersampling patterns and sparsifying transforms [START_REF] Menzel | Accelerated diffusion spectrum imaging in the human brain using compressed sensing[END_REF], [START_REF] Paquette | Comparison of sampling strategies and sparsifying transforms to improve compressed sensing diffusion spectrum imaging[END_REF].

Although 3 He lung diffusion images are sparser in the bdirection than in the spatial domain, the feasibility of exploiting sparsity along both the spatial encoding and the b-value directions has not been studied in hyperpolarized gas MRI. In other MRI applications, such as in cardiac cine MRI [START_REF] Lingala | Accelerated dynamic MRI exploiting sparsity and low-rank structure: k-t SLR[END_REF]- [START_REF] Abascal | Comparison of Total Variation with a Motion Estimation Based Compressed Sensing Approach for Self-Gated Cardiac Cine MRI in Small Animal Studies[END_REF] and fMRI [START_REF] Chavarrías | Exploitation of temporal redundancy in compressed sensing reconstruction of fMRI studies with a prior-based algorithm (PICCS)[END_REF], [START_REF] Chiew | k-t FASTER: Acceleration of functional MRI data acquisition using low rank constraints[END_REF], CS leads to large acceleration factors by exploiting high data dimensionality.

In this study, we propose a novel CS method that incorporates a model of the signal decay as prior information into the reconstruction (SIDER) to accelerate the acquisition of diffusion-weighted MRI data by undersampling in both the spatial and the b-value dimensions. We incorporate knowledge of diffusion signal behavior into the reconstruction to accelerate the acquisition of diffusion MRI data. The method we propose is compared with TV and zero filling reconstructions by assessing its effect on the estimated parameters of a stretched exponential model, which has been used to estimate mean alveolar dimensions (MAD) [START_REF] Parra-Robles | Characterization of 3 He Diffusion in Lungs using a Stretched Exponential Model[END_REF], [START_REF] Parra-Robles | Quantification of lung microstructure in asthma using a 3He fractional diffusion approach[END_REF]. The methods were applied in control and COPD patient data sets (n=8 in total) using retrospective undersampling simulations, with fully sampled data as the gold standard.

II. METHODS

A. Image reconstruction methods 1) Total variation

Previous CS studies for MRI using hyperpolarized gases assumed that each ventilation image ui is sparse under a transformation , which accounts for spatial sparsity [START_REF] Ajraoui | Compressed sensing in hyperpolarized 3He lung MRI[END_REF]- [START_REF] Chan | Whole lung morphometry with 3D multiple b-value hyperpolarized gas MRI and compressed sensing[END_REF]. The most common choice for  is the gradient that leads to the total variation functional. If F represents the undersampled Fourier transform and fi represents the undersampled k-space corresponding to the i-th b-value, then the total variation problem is given by ,

where =(x,y), i=1,…,B, and B is the total number of bvalues.

2) Signal decay-based reconstruction

We propose a novel CS method based on undersampling in both the spatial and the b-value dimensions that incorporates a model of the signal decay into the reconstruction (SIDER). It combines TV with a penalty function that promotes sparsity across the b-direction as follows:

(

where u and f are the ventilation images and undersampled data corresponding to all values of b, u=[u1

H ,…,uB H ] H , f =[f1 H ,…, f B H ] H
, u H denotes the Hermitian transpose of u, F is the undersampled multislice Fourier transform, and M is an operator that encodes the relationship between ventilation images for consecutive values of b. This relationship can be approximated using a stretched exponential model [START_REF] Chan | Whole lung morphometry with 3D multiple b-value hyperpolarized gas MRI and compressed sensing[END_REF], [START_REF] Parra-Robles | Characterization of 3 He Diffusion in Lungs using a Stretched Exponential Model[END_REF] as [START_REF] Deninger | Quantitative measurement of regional lung ventilation using 3 He MRI[END_REF] where uj is the ventilation image for the j-th b-value and 𝐷 and are the estimated average value of diffusivity and the heterogeneity index, respectively, across the whole volume, which can be obtained from a previously reconstructed image (here we used the image provided by the TV method). 1 This global estimation for the model parameter provides sufficient a priori information to guide the image reconstruction. The transform M is a generalization of the gradient across the bdimension, which takes the signal decay into account. Ventilation images are actually sparser on the stretched exponential domain defined by M than under the spatial gradient domain.

3) Split Bregman formulation

The problems in (1) and (2) were solved using the Split Bregman formulation, which efficiently handles L1-based constrained problems [START_REF] Montesinos | Application of the compressed sensing technique to selfgated cardiac cine sequences in small animals[END_REF], [START_REF] Abascal | Comparison of Total Variation with a Motion Estimation Based Compressed Sensing Approach for Self-Gated Cardiac Cine MRI in Small Animal Studies[END_REF], [START_REF] Osher | An Iterative Regularization Method for Total Variation-Based Image Restoration[END_REF], [START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF]. Using this formulation, constrained problems are converted to equivalent unconstrained problems, where constraints are imposed

i i i i u u Fu f     2 2 1 1 min such that , u u Mu Fu f              1 1 exp , j j j j j Mu u u Db Db              
 iteratively using the Bregman iteration. L2-and L1-norm functionals are separated into several subproblems, which are solved analytically in alternating steps. The subproblem including the L2-norm functionals results in a linear system that can be efficiently solved using iterative Krylov solvers; subproblems including L1-norm functionals are solved using shrinkage formulas. As TV can be obtained from SIDER by making β=0, we develop the formulation for the general case of SIDER.

To perform the split, we include the new variables dx, dy, and w and formulate a new problem that is equivalent to (2):

(4) Equation ( 4) is now easily managed using an equivalent unconstrained optimization approach, in which constraints are imposed by adding Bregman iterations bi, [START_REF] Wu | Accelerated MR diffusion tensor imaging using distributed t sensing[END_REF] where µ is a parameter that weights the data fidelity term,  is another parameter that weights the terms imposing the constraints for the dummy variables, k is the iteration number and the Bregman iterations are updated as [START_REF] Deninger | Quantification of regional intrapulmonary oxygen partial pressure evolution during apnea by (3)He MRI[END_REF] Since u and auxiliary variables w, dx, and dy are independent of each other, equation ( 5) can now be split into several equations (one for each variable) that are solved sequentially, as follows: [START_REF] Marshall | 3)He pO2 mapping is limited by delayed-ventilation and diffusion in chronic obstructive pulmonary disease[END_REF] Since the solution of u only involves L2-norm functionals, it can be obtained exactly as the solution of the linear system [START_REF] Deppe | Combined measurement of pulmonary inert gas washout and regional ventilation heterogeneity by MR of a single dose of hyperpolarized 3He[END_REF] where w H denotes the Hermitian transpose of w. It is common practice to modify [START_REF] Marshall | 3)He pO2 mapping is limited by delayed-ventilation and diffusion in chronic obstructive pulmonary disease[END_REF] by adding a quadratic term to improve the stability of the algorithm [START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF]. This replaces the update for u k+1 (8) as follows: [START_REF] Walker | Spin-exchange optical pumping of noble-gas nuclei[END_REF] Note that (9) constitutes a very large-scale problem, where K=NxN and N is the number of pixels, yet it can be solved efficiently using a Krylov solver, such as the biconjugate gradient stabilized method, which involves only matrix-vector multiplications: [START_REF] Parnell | Enhancement of X129e polarization by off-resonant spin exchange optical pumping[END_REF] The auxiliary variables dx, dy, and w are solved analytically using shrinkage formulas, which are thresholding operations (1,2) (11)

B. Datasets and retrospective undersampling

Eight fully sampled diffusion datasets were available from a previous study, 3 from healthy volunteers [START_REF] Parra-Robles | The influence of field strength on the apparent diffusion coefficient of 3He gas in human lungs[END_REF] and 5 from three patients with COPD (2 patients had 2 acquisitions at different sessions) [START_REF] Parra-Robles | Characterization of 3 He Diffusion in Lungs using a Stretched Exponential Model[END_REF]. Data consisted of 5 slices (10 mm thick with 10 mm gap between slices), 64x64 resolution and 5 b-values (0, 1.6, 3.2, 4.8 and 6.4 s/cm 2 ), obtained with a diffusion time of 1.6 ms. These data were acquired in a GE HDx 1.5 T scanner (GE Healthcare, USA) using a 2D spoiled gradient echo sequence with bipolar diffusion sensitization gradients [START_REF] Parra-Robles | The influence of field strength on the apparent diffusion coefficient of 3He gas in human lungs[END_REF],

  2 2 1 , , , 1 min , such that , , , x y
x y u d d w [START_REF] Parra-Robles | Characterization of 3 He Diffusion in Lungs using a Stretched Exponential Model[END_REF]. 3 He gas was polarized using a SEOP commercial polarizer (Helispin, GE Healthcare, USA), which achieved polarizations of 30-40%. These datasets were retrospectively undersampled to simulate CS acquisition and reconstruction. Quasi-random undersampling patterns were created using a variable density pattern as proposed in [START_REF] Lustig | Sparse MRI: The application of compressed sensing for rapid MR imaging[END_REF], [START_REF] Montesinos | Application of the compressed sensing technique to selfgated cardiac cine sequences in small animals[END_REF]. The patterns were based on the polynomial probability density function of the form (1-r) p , which assigned the highest probability to the center of the kspace and then assigned probabilities with a decay rate of p as a function to the center. Randomization was performed in the phase-encoding direction and through the b-direction (Fig. 1), thus enabling us to exploit data redundancy in 2 dimensions. We analyzed the results for acceleration factors of x2, x4, x5, x7, and x10. The decay rate was set to p=4 (for x2), p=4 (for x4), p=4 (for x5), p=6 (for x7) and p=9 (for x10).

x x y y d d w Fu f d u d u w Mu              2 1 2 , , , 1
2 2 2 2 2 2 min , 2 2 2 2 x y k x y u d d w k k k x x x y y y w d d w Fu f d u b d u b w Mu b                     1 1 1 1 1 1 1 1 1 1 1 k k k k x x x x k k k k y y y y k k k k w w k k k b b u d b b u d b b Mu w f f f Fu                            2 2 1 2 2 2 2 2 2 2 1 1 1 1 2 2 1 2 2 1 1 1 2 arg min 2 2 2 2 , arg min ( , ) 2
2 arg min 2 w x y w k k k k x x x u k k k k y y y k k k k x y x y x x x d ,d k k y y y k k k w u Fu f d D u b d D u b w Mu b d d d d d D u b d D u b w w w Mu b                                               1 H H H H H H H H k k x x y y k k k k k k x x x y y y k k w Ku r K F F D D D D M M r F f D d b D d b M w b                            1 H H H H H H H H k k x x y y k k k k k k x x x y y y k k k w Ku r K F F D D D D M M I r F f D d b D d b M w b u                                  k H H H x x y y H r F Fu D D u D D u M Mu u                    1 1 2 2 1 1 1 1 1 1 max , 0 , , , shrink , max / , 0 sign k k j j k k j k k k k k k x x y y k k k w k k k k w w D u b d s j x y s s D u b D u b w Mu b Mu b Mu b                             

C. Image analysis and evaluation

To evaluate the results, we first fitted the reconstructed signal u(b) on a pixel-by-pixel basis to the stretched exponential model, which estimates maps of the distributed diffusion coefficient D and heterogeneity index α [START_REF] Parra-Robles | Characterization of 3 He Diffusion in Lungs using a Stretched Exponential Model[END_REF]:

. ( 12
)
As reconstructed images u(b) are noisy, especially for patients and larger values of b, images u(b) were smoothed using a Gaussian filter (window of 3 neighboring pixels and SD of 1 pixel) before fitting the model in [START_REF] Lee | Advantages of parallel imaging in conjunction with hyperpolarized helium--a new approach to MRI of the lung[END_REF]. D and α were estimated only within a mask that had been created by segmenting the lungs in ventilation images from fully sampled data. Then, we estimated the MAD, Lm, from D and α as described in [START_REF] Chan | Whole lung morphometry with 3D multiple b-value hyperpolarized gas MRI and compressed sensing[END_REF], [START_REF] Parra-Robles | Quantification of lung microstructure in asthma using a 3He fractional diffusion approach[END_REF]. Lm was estimated only in the ranges 0<D<0.9 and 0.3<α<1, which were assumed to be physically meaningful. (The upper limits for these parameters correspond to diffusion within structures much larger than acinar airspaces, i.e. free diffusion with D0=0.87 cm 2 /s and α=1). Areas of the image that were outside this range appear black on the estimated maps. The process for estimating Lm from ventilation images, u(b), as well as the differences between the patient and control data sets are shown in Fig. 2.

All retrospectively undersampled data sets (n=8 for acceleration factors x2, x4, x5, x7, and x10) were reconstructed with zero filling (ZF), TV and SIDER. These methods were evaluated in terms of the following metrics: 1) relative meansquare error (MSE) of the recovered ventilation images, 2) relative MSE of MAD, 3) histograms and 4) mean values of the estimated maps of MAD. MSE was computed by adopting the images and maps obtained from the fully sampled data as the gold standard. The results are given as mean and SD across all data sets. To assess the statistical significance of the difference between the methods, we used a Mann-Whitney test, which is robust and avoids the assumption of normality in the data.

Histograms and images are shown for 1 control and 1 patient data set for all methods. The mean and SE of Lm in a region of interest (1 slice) are shown for the 3 patients to verify that errors due to undersampling were smaller than patient variability.

To determine the maximum acceleration factor achieved by the SIDER method, we studied its ability to discriminate between patients and controls and compared it with the reproducibility of the experiment. for the SIDER method and for the maximum acceleration factor.

To validate the feasibility of SIDER for differentiation between control and COPD patients, we run a Mann-Whitney test to assess statistical significance between control and patients and between SIDER and fully sampled data.

III. RESULTS

1) Selection of the reconstruction algorithm parameters

The parameters necessary for the TV and SIDER methods are µ, , γ, and the number of iterations. The parameter µ that weights the data fidelity term was selected following suggestions from previous studies [START_REF] Abascal | Comparison of Total Variation with a Motion Estimation Based Compressed Sensing Approach for Self-Gated Cardiac Cine MRI in Small Animal Studies[END_REF], [START_REF] Osher | An Iterative Regularization Method for Total Variation-Based Image Restoration[END_REF], [START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF], [START_REF] Montesinos | Application of the compressed sensing technique to selfgated cardiac cine sequences in small animals[END_REF]- [START_REF] Abascal | A Novel Prior-and Motion-Based Compressed Sensing Method for Small-Animal Respiratory Gated CT[END_REF], which showed that the algorithm converges to the same solution independently of its value, as long as it is sufficiently small. In [START_REF]Fluorescence diffuse optical tomography using the split Bregman method[END_REF], a two-fold decrease in µ doubled the number of iterations needed to converge to the same solution. Hence, the larger the value of µ, the faster the convergence, although if µ is excessively large, the algorithm converges to a noisy solution on the first iterations. In contrast to unconstrained optimization, where optimal selection of the regularization parameters is needed based on the L-curve or a similar method [START_REF] Hansen | The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems[END_REF], [START_REF] Abascal | Comparison of methods for optimal choice of the regularization parameter for linear electrical impedance tomography of brain function[END_REF], we only need to select the number of iterations (stopping criterion).

The stopping criterion used in previous studies was the discrepancy principle [START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF] or a fixed number of iterations [START_REF] Abascal | Comparison of Total Variation with a Motion Estimation Based Compressed Sensing Approach for Self-Gated Cardiac Cine MRI in Small Animal Studies[END_REF]. In order to compare TV and SIDER in this study, the number of iterations was chosen to minimize the MSE, considering the fully sampled image as the correct solution. For µ≤2, the method converged to the same solution at different iteration numbers. For µ>2 overfitting occurred at the first iterations and an optimal solution was not found. Increasing  (higher weight to the TV constraint) lowers the threshold in the shrinkage formulas (11), resulting in slightly faster convergence, but for >20, the convergence became unstable. Values of <1 resulted in low smoothing at the first iterations and an excessively high threshold in the shrinkage formulas, which also led to unstable convergence and large errors. Values in the range 1≤≤20 led to similar results. The parameter γ makes it possible to keep the Hessian matrix K (9) well-conditioned when µ is relatively larger than . A small value of γ is usually chosen (10 -4 ≤γ≤10 - 2 ).

The weighting parameters that control the relative degree of sparsity between TV (α) and model decay sparsity (β) were heuristically determined as follows. Increasing α>1 (higher sparsity to TV) led to images largely affected by cartoon-like artifacts and large solution errors. Selecting α<0.1 led to excessively fast convergence, compromising the robustness of the algorithm. Values of α in the range 0.1<α<1 suppressed most noise in the image and yielded similar results. Increasing β imposed higher sparsity on the model of signal decay than on TV, leading to lower solution errors, images with fewer cartoon-like artifacts and images less affected by noise (especially images corresponding to large b-values). However, the larger the value of β, the slower the convergence. We found that values in the range 0.2≤β≤1 were a good compromise. Table 1 shows a summary of the parameter values used for both TV and SIDER. In order to use the same range of parameters for all data sets, the algorithm normalizes the data item f as f/ǁf/nǁ, where n is the square root of the number of pixels in the image, following the suggestions from [Tom Goldstein. Split Bregman.

Retrieved in 2009 from http://www.ece.rice.edu/~tag7/Tom_Goldstein/Split_Bregman .html].

B. Comparison of methods

Fig. 3 (left) shows the MSE of the reconstruction of ventilation images (for b=0) with the various methods for all the acceleration factors tested. SIDER led to a significantly lower MSE than ZF and TV in all cases (for b=0); SIDER also led to a lower MSE for large values of b, although the differences were significant only at high accelerations. SIDER still preserved image quality with small errors in the estimated maps at x10. Fig. 5 shows the histograms of estimated MAD from images reconstructed with the algorithms for a control and patient dataset. ZF and TV led to larger errors as the acceleration increased, becoming very noticeable above x7. In contrast, SIDER presented histograms close to the target for accelerations up to x10. Fig. 6 shows the mean and SE of MAD across one slice for three patients. For accelerations larger than five, errors by ZF and TV were larger than differences between datasets. In contrast, SIDER led to small errors for all acceleration factors.

C. Analysis of the acceleration factor in retrospective data

Fig. 7 shows a slice of the map of estimated MAD and image differences obtained by the SIDER method for the remaining datasets (two controls and four patient data sets) and for the highest acceleration factors. Image quality was preserved for all accelerations, although there were small deviations from the fully sampled map for an acceleration factor of x10.

Errors for the estimation of the mean and the histogram of Lm are shown in Fig. 8. SIDER provides a very accurate estimation of the mean value of Lm, with errors below 1% for all accelerations. These errors were comparable to the experimental error (dotted line). An estimation of the histogram led to a value that was larger than the mean value. The crosscorrelation error for x7 was 2%, which was similar to the experimental error. Higher accelerations (x10) led to larger errors.

Fig 9 . shows histograms of the estimated MAD obtained from fully sampled data and using the SIDER method for acceleration factor x7 in all datasets. Histograms were very close to those of fully sampled data, and errors derived from the undersampling were much smaller than the intragroup differences between the control and patient datasets. Fig. 10 shows that the proposed method can be used to differentiate between controls and COPD patients based on the mean and FWHM of Lm (p<0.01). There were no significant differences between the proposed method and the complete data for these metrics (p<0.01).

D. Computation time

The code was implemented in MATLAB on a Windows computer with a 64-bit operating system, i7-3770, 3.40 GHz CPU and 16 GB RAM. SIDER took 28 s to reconstruct all ventilation images (for 5 values of b) for one slice. We used a straightforward parallelization of SIDER over the five slices to reduce computation time.

IV. DISCUSSION

We proposed and validated a novel CS method that incorporates a model of the signal decay into the reconstruction method as prior information. The proposed method, SIDER, was validated on both control and COPD patient data (n=8). Using retrospective undersampling, we found that accelerations of at least x7 are achievable with a negligible effect on the estimates of ventilation images and estimated MAD maps. This acceleration factor (x7) is very relevant, as it could be used to increase the resolution by two-fold and simultaneously in the xy plane, in the number of slices and in the number of b-values. Such an approach would not only increase the resolution and volume coverage, but also improve the accuracy of estimated microstructural dimensions and may enable the use of models with a larger number of parameters.

The acceleration factor depends on the criterion chosen to decide whether a result is acceptable. Adopting the experimental repeatability error as an upper bound, SIDER accurately calculated the mean of Lm up to x10 and of the entire histogram up to x7. Nevertheless, this criterion may be very restrictive, as the experimental repeatability is very high. In fact, we showed that MSE remained small and histograms were similar to those obtained from fully sampled data for acceleration factors up to x10. Image quality was preserved with small errors in the estimated maps at a factor x10. A more practical criterion would be the ability to differentiate between controls and patients. In this case, SIDER provided the same results for x7 as those provided by fully sampled data, in terms of mean and FWHM of Lm, which are used in clinical practice. For accelerations higher than x10, errors were larger and the number of encoding lines acquired was very small for the present resolution (images of 64x64), which was not considered acceptable. An acceleration factor of x7, which would enable doubling of the spatial resolution, number of slices and the number of b-values, presented negligible errors.

Previous CS implementations achieved lower acceleration factors: x2 using spatial TV in 2D [START_REF] Ajraoui | Compressed sensing in hyperpolarized 3He lung MRI[END_REF], x3 using spatial TV in 2D with prior knowledge of a proton image acquired during the same breath-hold [START_REF] Ajraoui | Incorporation of prior knowledge in compressed sensing for faster acquisition of hyperpolarized gas images[END_REF], and x3 using TV in 3D [START_REF] Chan | Whole lung morphometry with 3D multiple b-value hyperpolarized gas MRI and compressed sensing[END_REF]. The acceleration factors achieved in this study, from x7 to x10, are superior because of the exploitation of undersampling along the b-dimension and the use of a reconstruction model that incorporates prior knowledge of signal decay. SIDER could also be extended to the 3D case as in [START_REF] Chan | Whole lung morphometry with 3D multiple b-value hyperpolarized gas MRI and compressed sensing[END_REF], potentially achieving even higher accelerations.

In this study, we aimed to exploit sparsity not only in the spatial domain, but also across the b-dimension, by both modeling the signal decay and randomizing the acquisition along the b-direction. The proposed SIDER method is a combination of spatial TV and a functional that models signal decay. Given that SIDER and TV are both solved using the split Bregman formulation, the only difference between these two methodologies lies in the exploitation of the b-dimension. TV has been chosen as spatial regularization, which had been previously proposed for hyperpolarized gas MRI [START_REF] Ajraoui | Compressed sensing in hyperpolarized 3He lung MRI[END_REF]- [START_REF] Chan | Whole lung morphometry with 3D multiple b-value hyperpolarized gas MRI and compressed sensing[END_REF]. Other functionals for spatial regularization, such as total generalized variation (TGV) [START_REF] Bredies | Total Generalized Variation[END_REF], [START_REF] Knoll | Second order total generalized variation (TGV) for MRI[END_REF], could lead to results that are superior to those of TV. Nevertheless, TGV or different spatial regularization could be included in the proposed methodology by replacing the TV functional.

Our study is subject to several limitations. Various approximations have been adopted for modeling signal decay. First, we assumed a stretched exponential model ( 3) with global parameters D and averaged across the whole volume. Second, model parameters were estimated from a previously reconstructed image, which we chose as that provided by the TV method. We performed tests to assess these limitations. We compared the use of the global parameters obtained from TV and from SIDER (by running SIDER twice), although no significant improvements were observed. This could be explained by the fact that average parameters provided by TV were sufficient to guide the reconstruction algorithm. We also found that using zero filling instead of TV led to large errors for high acceleration factors. Thus, zero filling could not provide a good estimate of average parameters for the stretched exponential model. A more simplified assumption could be to model the signal decay for consecutive values of b with an exponential model, .

In this case, using zero filling to estimate D across the whole volume led to reasonably good results. This approximation is faster and simpler, as it requires a less accurate estimation of D , although it provided slightly worse results for some values of b than the proposed model. A plausible improvement in the method proposed could be to model the parameters of the stretched model (D and α) on a pixel-by-pixel basis (instead of average parameters) and to include them as unknowns in the reconstruction procedure. This method could use the reconstructed image to update the signal decay model in an iterative fashion and iterate until convergence. Another improvement could be to rewrite the reconstruction problem in terms of D and α or in terms of Lm. However, such an approach lies beyond the scope of the present study.

Another limitation is that we validated the method by retrospectively undersampling control and patient data. However, implementation of the proposed CS sequence in an MR scanner will be validated in future studies. In addition, the potential gain in microstructure information is promising but requires a dedicated study.

In future studies, this method will be extended to the acquisition of 3D diffusion data, where higher accelerations can be achieved through undersampling over 3 directions (e.g. 2 spatial directions plus b-value), and to other hyperpolarized gas MR applications, where a model of the signal behavior has been assessed (e.g., pO2 mapping [START_REF] Marshall | 3)He pO2 mapping is limited by delayed-ventilation and diffusion in chronic obstructive pulmonary disease[END_REF]). Our method is not restricted to hyperpolarized gas imaging, but could also be used in other diffusion MRI applications and in metabolic imaging based on hyperpolarized 13 C.

In conclusion, we validated a novel CS method for lung MRI imaging and achieved high acceleration factors that can be used to increase the amount of data acquired during breath-hold. This methodology is expected to improve the accuracy of estimated lung microstructure data and widen the possibilities of studying lung diseases with MRI.

Fig. 1 .

 1 Fig. 1. Undersampling pattern for an acceleration factor of x7. Randomization is performed in the phase-encoding direction and through the b-direction, for a 64x64 resolution and 5 b-values.

Fig. 2 .

 2 Fig. 2. Process for estimating mean alveolar dimensions, Lm, from ventilation images, u(b), for datasets from a control and COPD patient. Top left: slice of ventilation images for different values of b. Top right: mean value of u(b) across the slice to illustrate the typical difference in signal decay between controls and patients. Bottom: Estimated maps of distributed diffusion coefficient D, heterogeneity index α (12) and Lm for one control and one patient.

Fig. 3 .Fig. 5 .

 35 Fig. 3. MSE of ventilation images for b=0 (left) and estimated mean alveolar dimension (right) versus the acceleration factor for the different reconstruction algorithms. Results show mean and SD across all data sets (n=8).

Fig 4 .

 4 Fig 4. One slice of an estimated mean alveolar length map (left) and image differences (right) for each of the methods and acceleration factors, for a control subject (top) and for a COPD patient (bottom). The arrows point to areas where errors are more visible: (1) Black regions where estimated values were outside what is physically reliable, (2) shift in the mean values (3) and other errors in the estimated maps.

Fig. 6 .

 6 Fig. 6. Mean and SE of mean alveolar length across one slice vs. the acceleration factor for three different patients. Solid and dashed lines represent mean and SE, respectively, for fully sampled data.

Fig. 7 .

 7 Fig. 7. One slice of estimated mean alveolar dimensions (left) and image differences (right) with SIDER for four patient and two control datasets.

Fig 8 .

 8 Fig 8. Left: Percentage error of the estimated average of the mean alveolar dimensions (Lm) versus the acceleration factor. The mean and SD values shown were computed across all datasets. The dotted line represents the experimental error, which was estimated from two patients who had two acquisitions acquired at different sessions. Right: Percentage error of the estimation of the histogram of Lm. The estimation of the histogram is measured by the cross-correlation between the histogram for each acceleration factor and the histogram for the complete data. The error is computed with respect to the autocorrelation for the complete data.

Fig. 9 .

 9 Fig. 9. Histograms of estimated mean alveolar dimensions for fully sampled data (left) and SIDER method for an acceleration factor x7 (right). Controls are depicted using dashed lines and COPD patients with solid lines.

Fig 10 .

 10 Fig 10. Comparison between COPD patients and controls for mean and FWHM of histograms of Lm. These results correspond to the histograms shown in Fig. 9.

  This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2017.2765281, IEEE Transactions on Medical Imaging

4

© 20xx IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. DOI 10.1109/TMI.2017.2765281

Results and MATLAB code for TV and SIDER methods are available from the GitHub repository https://github.com/HGGM-LIM/compressed-sensingdiffusion-lung-MRI

© 20xx IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. DOI 10.1109/TMI.2017.2765281

© 20xx IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works. DOI 10.1109/TMI.2017.2765281

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2017.2765281, IEEE Transactions on Medical Imaging

ACKNOWLEDGMENT

The MRI diffusion data used in this work were kindly made available by Prof. Jim M. Wild (POLARIS, Academic Unit of Radiology, University of Sheffield, UK).

This project has received funding from the Universidad Carlos III de Madrid, the European Union's Seventh Framework Programme for research, technological development and demostration under CONEX programme, Banco Santander an other external funds: Marie Curie-COFUND Programme (PCOFUND-GA-2012-60371) and Ministerio de Economía y Competitividad (COFUND2013-40258, TEC2014-56600-R), Ministerio de Educación, cultura y Deporte (CEI-15-17), Comunidad de Madrid (BRADE-CM S2013/ICE-2958) and Instituto de Salud Carlos III, cofinanciado con fondos FEDER, Una manera de hacer Europa (PI16/02037). The research leading to these results has received funding from the Innovative Medicines Initiative (www.imi.europa.eu) Joint Undertaking under grant agreement n°115337, resources of which are composed of financial contribution from the European Union's Seventh Framework Programme (FP7/2007-2013) and EFPIA companies' in kind contribution. This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement N° 701915.


Veronica Quintuna Rodriguez
email: veronicakarina.quintunarodriguez@orange.com

Fabrice Guillemin
email: fabrice.guillemin@orange.com

Towards the deployment of a fully centralized Cloud-RAN architecture

Keywords: Cloud-RAN, NFV, VNF, BBU, channel decoding, multi-core, global scheduling

In the framework of Network Function Virtualization (NFV), we address in this work the design and sizing of Cloud-RAN architectures. We concretely investigate the execution time of software-based Base Band Units (BBUs) on multicore systems. Since Cloud-RAN requires a real-time behavior, we use parallel programming techniques in order to minimize the runtime of BBU functions. For an efficient utilization of computing resources, we investigate the relevance of resource pooling where a global scheduling algorithm allocates processing units to runnable BBU-jobs. We specifically examine the gain that can be obtained when applying data parallelism on the channel decoding BBU-function which is the most expensive one in terms of processing time. Performance results show a significant reduction in the runtime of PHY functions which enables the deployment of a fully centralized Cloud-RAN architecture.

I. INTRODUCTION

Network Function Virtualization (NFV) [START_REF] Han | Network function virtualization: Challenges and opportunities for innovations[END_REF] is an emerging paradigm which advocates the separation between network functions and hosting hardware. This approach is intended to make the network architecture and operation more flexible since network functions can be instantiated on the fly in cloud infrastructures and no more on dedicated servers via complex procedures. Evolved Packet Core (EPC), IP Multimedia Subsystem (IMS), etc. are typical candidates for virtualization. In this paper, we pay special attention to Radio Access Network (RAN) functions. The cloudification of RAN functions is of utmost interest for several reasons. In the context of network slicing, RAN functions could be instantiated on the fly for a specific use case (e.g., a temporary cellular network for the needs of a company). Cloud-RAN enables the deployment of radio networks with a tailored behavior for a particular service or client, i.e., RAN as a Service (RANaaS).

Cloud-RAN is based on network functions which are defined in software. They run on general purpose hardware, e.g., GPU/CPU based servers. The cloudification of RAN can be seen as a key architectural evolution of mobile networks (especially in areas with a dense mesh of small cells). It allows hundreds of remote Radio Remote Heads (RRHs) to be connected to a centralized BBU-pool hosted in a Central Office (CO).

For a network operator, the main advantage of Cloud-RAN is the native support of collaborative radio technologies such as Coordinated Multi-point (CoMP), coordinated beam-forming, joint transmission, among others. These advanced algorithms are difficult to implement in current mobile architectures due to the latency between cell sites. Resource saving is another important benefit of this evolution. Energy consumption reduction is possible via BBU pooling and even thanks to coordinated ON/OFF switching strategies of eNodeBs (eNBs). Computing resources in data centers can also be saved via resource pooling. The key point is to apply statistical multiplexing when dimensioning future COs.

Centralizing BBUs brings some challenges, a key one being the front-haul sizing. The front-haul dimension represents the distance between the CO (BBUs) and antennas (RRHs). This size is constrained by the Hybrid Automatic Repeat reQuest (HARQ) process defined in LTE, i.e., the Round Trip Time (RTT). As a matter of fact, Cloud-RAN requires a real-time behavior. The BBU should complete the base-band processing within 3 ms where 2 ms are available for the reception process and 1 ms for the transmission process. The front-haul capacity should also be studied. It depends on the configuration of each eNB hosted in the BBU-pool (e.g., MIMO, SISO, LTE bandwidth of 5 MHz, 20 MHz, etc).

The time budget (1 ms and 2 ms for down-link and uplink, respectively) should actually be shared between the propagation time and the processing time. If the processing time is reduced, the front-haul size can be increased and BBU functions can be located higher in the network, typically in Points of Presence at the edge of an IP network.

In [START_REF] Quintuna | VNF modeling towards the Cloud-RAN implementation[END_REF] a threading model has been proposed for reducing the execution time of BBU functions in the down-link (DL) direction, notably the Channel Coding function. In the present work, we continue along the same line of investigation by introducing a solution to minimize the execution time of BBU functions in the up-link (UL) direction. We concretely investigate the execution of BBU functions on multi-core systems when applying parallel programming. First, we identify the driving factors which determine the runtime of BBU functions. Then, we implement a threading model with data parallelism in order to decrease the execution time of BBU functions in the up-link direction. Finally, we evaluate by simulation the performance of virtualized BBU functions in terms of latency (i.e., processing time). The knowledge of runtime of BBU functions is needed for dimensioning the Cloud-RAN architecture, i.e., the computing resources of data centers, the front-haul size and capacity. This document is organized as follows: The Cloud-RAN architecture is discussed in Section II. The runtime of BBU functions and parallelization are described in Section III. The performance analysis in terms of latency is presented in Section IV. Concluding remarks are presented in Section V.

II. CLOUD-RAN ARCHITECTURE

Current mobile networks have a distributed architecture where BBUs are located near to antennas. BBUs are implemented on proprietary hardware and are provided by a single vendor. On the contrary, Cloud-RAN architectures are based on open-platforms where the base-band functions can be instantiated on demand. In the same way the computing resources can be dynamically allocated.

The Cloud-RAN architecture can support selective centralization of BBU functions. Several functional splits of the BBU have been widely considered in the literature [START_REF] Checko | Cloud RAN for mobile networks technology overview[END_REF], [START_REF] Nikaein | Processing radio access network functions in the cloud: Critical issues and modeling[END_REF], [START_REF] Sabella | RAN as a service: Challenges of designing a flexible RAN architecture in a cloud-based heterogeneous mobile network[END_REF]. We can roughly classify Cloud-RAN architectures as fully and partially centralized; in this paper, we focus our study in the case of full centralization which moves all base-band functions (BBUs) higher in the network. See Figure 1 for an illustration.

A. BBU-pool

The BBU-pool is deployed on commodity hardware, i.e., multi-core GPU/CPU-based servers. Base-band functions are defined in software and run as applications. Cloud-RAN employs virtualization technologies which can be based on Virtual Machines (VMs) and/or containers. In this work, we take advantage of the performance provided by containers which, unlike VMs run on a common single kernel. This gives them the benefit of being faster and more resource-efficient. This point has been studied in [START_REF] Islam Alyafawi | Critical issues of centralized and cloudified lte-fdd radio access networks[END_REF]. In a container-based Cloud-RAN platform, BBU functions are represented as runnable-tasks (processes or jobs) which are placed in the highest layer of the Cloud-RAN system as shown in Figure 1. These tasks correspond to all physical, MAC and network functionalities, i.e., IFFT/FFT, modulation/demodulation, encoding/decoding, radio scheduling, HARQ management, radio link control, data convergence procedures, and User Equipments (UEs) measurement reporting and paging maintained by the Radio Resource Control (RRC) protocol.

Since strict real-time constraints are present in the execution of BBU functions, the behavior of the Operating System (OS) plays a crucial role in the performance of Cloud-RAN systems. In fact, each BBU-function represents a chunk of jobs or tasks to be scheduled by a global scheduling algorithm embedded into the kernel of the OS. This algorithm determines the order of execution of BBU-jobs and allocates them computingresources.

B. Front-haul size

The front-haul size, i.e., the distance between the BBUpool and antennas, is limited by the time-budget of the RTT defined by LTE which includes the acknowledgment of each sub-frame. In LTE, the acknowledgment messages and retransmission procedures in case of errors are handled by the HARQ process. As shown in Figure 2, the BBU-pool has less than 3 ms for the whole base-band processing (namely, decoding, checking the Cyclic Redundancy Check (CRC), and encoding the ACK/NACK). The reception process (Rx) has a budget of 2 ms and the transmission process (Tx) 1 ms, denoted by T Rx and T T x , respectively; the turnaround time is 8 ms.

In fact, to dimension the front-haul size, it is first necessary to know the response time T r of the BBU-pool, i.e., the required time to execute all BBU functions. Thus, the timebudget for the propagation of IQ signals, i.e., the front-haul delay, is the remaining time after the base-band processing in the BBU-pool. Since the BBU-pool (CO) is linked with antennas by optic-fiber, the front-haul size d can easily be obtained from the front-haul time-budget, so-called front-haul delay T F h , and the speed of light c, as d = c * T F h .

The HARQ mechanism considers an advancing time T A in order to align signals in time due to propagation delay between the UE and the eNB. In LTE there are 8 HARQ processes executed at the same time with an offset of 1 ms each which corresponds to the acquisition time of a sub-frame, T Aq . Thus, T HARQ = RT T + T A , where

RT T = T T x + T Rx + 2 * T Aq + BBU proc -2 * T F h .

C. Front-haul capacity

The front-haul capacity is defined by the number of eNBs hosted in the CO. The current widely used protocol for data transmission between antennas and BBUs is Common Public Radio Interface (CPRI) which transmits IQ signals. The transmission rate is constant since CPRI is a serial Constant Bit Rate (CBR) interface. It is then independent of the mobile network load [START_REF] Duan | Performance analysis of several functional splits in C-RAN[END_REF]. The data throughput depends on the cell configuration and can be obtained as

R IQ = 2 * M * f s * F coding * F control * N ant * N sec ,
where the factor 2 corresponds to I and Q radio signals; M is the number of bits per sample used in the quantisation process (current LTE implementations use M = 15 [START_REF] Fabio | 4G mobile networks: An analysis of spectrum allocation, software radio architectures and interfacing technology[END_REF]); f s is the sampling frequency, which is a multiple of the nominal chip rate of LTE, f c = 3.84 MHz, e.g., the sampling frequency for 10 MHz and 20 MHz is 15.36 MHz and 30.72 MHz, respectively; F coding is the line coding factor which can be either 10/8 or 66/64 [START_REF][END_REF]; F control is the control factor which defines the number of data words and control words (F control = 16/15 when using CPRI); N ant and N sec are the number of antennas and sectors, respectively.

For example, an eNB with 3 sectors and 8x8 MIMO antennas needs a data rate of 29.49 Gbps for 20 MHz. Research studies are focusing on bandwidth compression techniques and packetisation of CPRI data via Ethernet [START_REF] Guo | LTE/LTE-A signal compression on the CPRI interface[END_REF]. Others solutions propose different functional split architectures in order to reduce the required bandwidth [START_REF] Wubben | Benefits and impact of cloud computing on 5G signal processing: Flexible centralization through cloud-ran[END_REF], [START_REF] Duan | Performance analysis of several functional splits in C-RAN[END_REF]. For example, when including the demapping process in the RRH, it is possible to adapt the bandwidth as a function of the traffic load in the cell, then the required front-haul capacity is directly given by the fraction of utilized radio resources [START_REF] Wubben | Benefits and impact of cloud computing on 5G signal processing: Flexible centralization through cloud-ran[END_REF].

In general, the required data rate significantly decreases when the functional split is shifted after the PHY layer or even after the MAC layer [START_REF] Wubben | Benefits and impact of cloud computing on 5G signal processing: Flexible centralization through cloud-ran[END_REF].

III. PROCESSING BBU FUNCTIONS IN THE CLOUD

This work is focused on the implementation of BBU functions carried out during the Rx process (up-link). It includes the channel decoding which is the most expensive sub-function in terms of latency [START_REF] Nikaein | Processing radio access network functions in the cloud: Critical issues and modeling[END_REF], [START_REF] Islam Alyafawi | Critical issues of centralized and cloudified LTE-FDD radio access networks[END_REF], [START_REF] Nikaein | Openairinterface 4G: an open LTE network in a PC[END_REF]. The performance of BBU functions for the Tx process (down-link) is studied in [START_REF] Quintuna | VNF modeling towards the Cloud-RAN implementation[END_REF].

In order to minimize the runtime of BBU functions, we investigate the relevance of parallel programing and resource pooling in multi-core systems. The main goal of the runtime reduction is to enlarge the distance between the BBU-pool and antennas while saving computing resources.

In [START_REF] Bhaumik | Cloudiq: A framework for processing base stations in a data center[END_REF], namely CloudIQ, is shown that at least 22% of computing resources can be saved only when using statistical multiplexing and resource pooling in Cloud-RAN systems. It exploits the variations in the processing load of individuals BBUs to use fewer computing resources.

While in CloudIQ, an LTE sub-frame is completely processed on a single computing resource, we are interested on splitting the sub-frame processing in parallel runnable tasks in order to decrease the runtime of the entire sub-frame, namely data parallelism. Note that functional parallelism cannot be applied since BBU functions (e.g, decoding, demodulation) need to be executed sequentially.

We propose the parallel processing of sub-frames as follows: The workload of a sub-frame is divided in slices where each of them corresponds to the data of a single UE. These slices, so-called sub-tasks, are executed simultaneously on different cores. For further performance improvement, we propose splitting the workload of a single UE in parallel runnable Code Blocks (CBs).

Hence, the resulting threading model executes one thread per UE, and/or, one thread per CB for the channel decoding function. The threading model is illustrated in Figure 3. A global scheduler allocates a dedicated single core to each subtask (either per UE or per CB) in order to avoid context switching overhead. The Rx base-band processing begins with the demapping process which selects the signal corresponding to a UE in order to demodulate and decode the received data. The MAC layer which handles the radio scheduling and HARQ mechanisms, performs demultiplexing of the Up-Link Shared Channel (UL-SCH) to restore the various logical channels containing the control messages and the IP packets. The size of data per UE in a sub-frame is determined by the radio scheduler which builds the resource grid of a cell. The radio scheduler allocates resources (time and frequency), namely Physical Resource Blocks (PRB), to a UE in function of the data transmission needs and the radio channel quality. The number of allocated PRB together with the Modulation and Coding Scheme (MCS) determine the Transport Block Size (TBS).

The execution time of each PHY function (i.e., FFT, demodulation, and decoding) depends of the TBS attributed to each UE. In other words, when applying data parallelism, the runtime of BBU functions (i.e., the Cloud-RAN performance in terms of latency) depends on the structure of the resource grid which is determined by the radio scheduler. Hence, the radio resource grid represents the workload of a Cloud-RAN system. It is important to note that scheduling strategies are not defined by LTE and as a consequence they are vendor specific. Most common radio schedulers takes into account the Channel Quality Indicator (CQI), the QoS Class Identifier (QCI), the number of spatial layers, Inter-Cell Interference Coordination (ICIC), among others.

IV. PERFORMANCE ANALYSIS

In this section we investigate the runtime of the PHY BBU sub-functions on a multi-core platform. In order to evaluate the gain in terms of latency, we apply data parallelism for the execution of PHY UL functions, i.e., FFT, demodulation and channel decoding. We use Open Air Interface (OAI) data [START_REF] Nikaein | Openairinterface 4G: an open LTE network in a PC[END_REF] in order to determine the runtime of BBU sub-functions.

OAI is an open-source solution that implements the RAN functionality in software and can be used as a benchmark.

A. Simulation settings

The simulation is performed for a single cell with SISO 20 MHz of bandwidth and Frequency Division Duplex (FDD) transmission mode. The simulator builds the resource grid considering the behavior of current commercial radio scheduler deployed in large cities. It selects both the MCS and the number of PRBs based on emulated radio conditions, the traffic in the cell, and the data load per UE. The simulator determines the TBS using tables 7.1.7.1-1 and 7.1.7.2.1.1-1 defined in the LTE Physical Layer Specification 3GPP TS 36.213 version 12.4.0 Release 12 [START_REF]LTE, evolved universal terrestrial radio access, physical layer procedures[END_REF]. In order to study the worst case, we consider non empty scheduling, i.e., all UL resources designated to a UE are used and need to be processed. We evaluate two scenarios.

1) Scenario A: Real traffic emulation during busy hours: This scenario is carried out under real traffic conditions. We emulate the behavior of current deployed eNBs. We use the MCS and Channel Quality Indicator (CQI) patterns performed by an eNB during busy-hours as shown in Figure 4(a).

In general, radio schedulers use the CQI to determine the MCS. When the MCS takes values from 0 to 9, data is modulated in QPSK, while for values between 10 -16 the eNB uses 16-QAM modulation. In order to build the resource grid we have also obtained the number of UEs scheduled by Transmission Time Interval (TTI) which is illustrated in Figure 4(b).

Note that only one or two UEs are scheduled each millisecond. A fragment of the obtained resource grid is shown in Figure 5(a).

2) Scenario B: No traffic patterns: In order to investigate the worst case, the simulator does not limit the behavior of eNBs to a specific traffic pattern. The resource grid is built under variable radio conditions and elastic number of connected UEs. The number of PRBs allocated to a UE ranges from 6 to 110 [START_REF]LTE, evolved universal terrestrial radio access, physical layer procedures[END_REF], [START_REF]LTE, evolved universal terrestrial radio access (e-utra) and evolved universal terrestrial radio access network(e-utran)[END_REF]. The MCS varies between 0 and 27 which enables QPSK, 16-QAM, and 64-QAM modulation orders [START_REF]LTE, evolved universal terrestrial radio access, physical layer procedures[END_REF]. The obtained resource grid is shown in Figure 5

B. PHY UL runtime

We evaluate the runtime of the PHY layer during the reception process. In a first step, we do not consider the data parallel model presented in Section III. Figure 6(a) and Figure 6(b) show the execution time for the demodulation, channel decoding and FFT sub-functions for Scenarios A and B, respectively. Besides the fact that the FFT runtime is constant and that the demodulation processing time is less fluctuating than the decoding process, it is worth noting that the runtime of the whole PHY layer is essentially determined by the decoding processing time. The FFT runtime only depends on the number of PRBs. The demodulation sub-function is directly influenced by the modulation order and the number of PRBs. The decoding process depends on the TBS, i.e., the CQI, the data load by UE and the traffic in the cell.

Simulation results show that the execution time of the PHY UL sub-functions can reach values of up to 1.8 ms (Scenario B). It does not leave enough margin to deploy a fully centralized RAN architecture which should additionally consider the front-haul delay and the processing time of L2 and L3 functions within the LTE Rx budget, i.e., 2 ms. In order to enable the deployment of a fully centralized RAN architecture, we evaluate in the following the performance gain obtained when using data parallelism.

C. Data parallelism: Performance gain

We focus our analysis on the most expensive function in terms of latency, i.e., the channel decoding process. We evaluate the performance gain when using data parallelism on a multi-core platform with c cores. We decompose the channel decoding function in a set of jobs where each of them processes a portion of data. In our proposed Cloud-RAN system, each runnable job accesses to computing resources according to the scheduling strategy. The scheduler selects them in first-in-first-out order. We use a non-preemptive scheduling algorithm which allocates a single core to each job. After allocation, a core is busy until that the job is completed, i.e., a job cannot be interrupted.

We evaluate the performance of the decoding sub-function under two degrees of data parallelism (1) by UEs and (2) by Code Blocks (CBs). In general, a sub-frame contains more than one UE, hence the PHY processing can be executed in parallel threads where each job corresponds to a UE. In the same way, when a Transport Block allocated to a UE exceeds 6120 bits, it is split in CBs of 6120 bits which can be executed in parallel. Hence, BBU-jobs arrive to the computing platform in batches. The batch size when applying data parallelism by UEs corresponds to the number of UEs per TTI. The batch size when processing CBs in parallel is the product of active UEs and the number of CBs per UE. The number of parallel threads for both (1) and ( 2) varies every millisecond in function of the number of UEs allocated by TTI and the data load per UE.

We study the runtime of the decoding sub-function when using a multi-core platform with 6 cores for both above described scenarios A and B. We evaluate the gain of the parallel execution by UEs and by CBs with respect to noparallelism, i.e., the serial execution of jobs. Let us consider the Scenario A which reflects the behavior of current deployed eNBs in busy hours. Here, the mean runtime of the decoding process can be reduced by a factor of five when using data parallelism by CBs. It is shown in Figure 7(a). When more than one UE is scheduled by TTI, data parallelism by UEs can also considerably decrease the execution time. Parallelism by UEs requires notably less developing effort than parallelism by CBs. Parallelism by UEs can then be considered as advantageous when taking into account the complexity in the implementation.

Let us now consider the Scenario B, performance results are shown in Figure 7(b). The runtime of the decoding function is divided by 3 when executing CBs in parallel and halved when applying threads per UEs.

The study of the probability density function of the channel decoding runtime shows that when executing one thread per CB, runtime values are more concentrated around the mean. Even if the reduction experimented when running UEs in parallel is already very interesting in terms of latency, the parallelism by CBs offers better performance, especially when comparing the tail of both probability density functions. See Figure 8 for an illustration. Note that when applying parallelism by UEs, the probability density function exhibits a long tail that is similar to that obtained when there is no parallelism.

In other words, parallelism by CBs presents less statistical dispersion, i.e., the execution of the channel decoding function is near to a specific value (e.g., 300 ms for a multi-core platform with 6 cores). This fact is crucial for the cloudification of RAN functions.

V. CONCLUSION

In this work, we have studied a fully centralized Cloud-RAN architecture which meets the need for improving the performance of software-based BBUs in terms of latency. To reduce the runtime of the most expensive BBU function, i.e., the channel decoding process, we have applied data parallelism and resource pooling in a multi-core system.

We have evaluated two methods of parallelism. The first one uses one thread per UE, and the second one goes further, employing one thread per CB. We have implemented a global scheduler to determine which core executes a particular job and in what order the BBU-jobs are processed.

Performance results show that a finer granularity in the parallelization model of the decoding function can significantly reduce the runtime and the variability. This performance gain enables longer distances between the BBU-pool and antennas, as well as, the efficient utilization of computing resources.

The main conclusion is that the cloudification of RAN functions for the up-link is feasible in dense areas, which is in accordance with [START_REF] Quintuna | VNF modeling towards the Cloud-RAN implementation[END_REF].

Fig. 1 .

 1 Fig. 1. Fully-centralized Cloud-RAN architecture.

Fig. 2 .

 2 Fig. 2. HARQ process in Cloud-RAN architectures.

Fig. 3 .

 3 Fig. 3. Threading model of BBU-UL functions using data parallelism.

Fig. 4 .

 4 Fig. 4. Radio scheduling indicators.

  (b).

Fig. 5 .

 5 Fig. 5. LTE resource grid (Cloud-RAN workload).

Fig. 6 .

 6 Fig. 6. Performance of virtual-BBUs (PHY layer, up-link).

Fig. 7 .

 7 Fig. 7. Performance of channel decoding when using data parallelism.

Fig. 8 .

 8 Fig. 8. Probability density of channel decoding runtime, Scenario B.