Amelia Carolina Sparavigna

Bernoulli Numbers: from Ada Lovelace to the Debye Functions

Keywords: History of Physics, Ada Lovelace, Debye Functions, Debye specific heat

Jacob Bernoulli owes his fame for the numerous contributions to calculus and for his discoveries in the field of probability. Here we will discuss one of his contributions to the theory of numbers, the Bernoulli numbers. They were proposed as a case study by Ada Lovelace in her analysis of Menabrea's report on Babbage Analytical Engine. It is probable that it was this Lovelace's work, that inspired Hans Thirring in using the Bernoulli numbers in the calculus of the Debye functions.

Introduction

Jacob Bernoulli (1655-1705) was one of the many prominent mathematicians in the Bernoulli family, a family of merchants and scholars, originally from Antwerp that resettled in Basel, Switzerland. Another member and prominent scholar of this family was Daniel Bernoulli that worked in fluid dynamics, publishing his results in the book Hydrodynamica in 1738.

Jacob Bernoulli owes his fame for the numerous contributions to calculus; along with his brother Johann, was one of the founders of the calculus of variations. Jacob's most important contribution is generally considered in the field of probability. In fact, he derived the first version of the law of large numbers in his work Ars Conjectandi [START_REF] Jacob | The MacTutor History of Mathematics archive[END_REF].

Here we will discuss one of his contributions to the theory of numbers, the Bernoulli numbers. They were proposed as a case study by Ada Lovelace in her analysis of Menabrea's report on Babbage Analytical Engine. Among their several uses, we have the calculus of the Debye functions, in theory of the lattice specific heat. Probably, their use had been inspired by the Lovelace's work.

Bernoulli and the integers

One of the numbers that Jacob Bernoulli investigated was e. In 1683, Bernoulli discovered the constant e by studying a question about compound interests, which required him to find the value of the following limit, having as a value e:

(1)

e n n =       + ∞ → 1 1 lim
Number e is an important mathematical constant that is the base of the natural logarithm. It is approximately equal to 2.71828. Sometimes, this number is called Euler's number after Leonhard Euler, but it is not to be confused with γ, the Euler-Mascheroni constant, sometimes called simply Euler's constant. In fact, number e is also known as the Napier's constant (see Appendix for e). However, here our aim is discussing the Bernoulli's numbers.

The manner in which Bernoulli discovered then is discussed in detail in [START_REF] Ward | s Mathematics Pages, Series Sums of the Powers of the Natural Numbers and Bernoulli Numbers[END_REF] by K. Ward. Bernoulli was working on the sums of the natural number power series. He noted certain facts about them and discovered a close formula by means of which the sum of a given power of natural number is calculated. This formula requires some special numbers, that is, the Bernoulli numbers. In his original formulae, Bernoulli used for then the letters A,B,C, etc.. Ward gives us the following formula:

(2) ....

3 4 1 2 1 1 3 4 1 1 2 1 1 0 +         +         + - + = - = - + ∑ m n k m m m m n m B n m B n B n m B k
In (2), we have (with p and q integers):

. )! ( !

! q p q p q p - =        

Ada Lovelace and the Bernoulli numbers

Augusta Ada King-Noel, Countess of Lovelace (née Byron; 1815-1852) was an English mathematician known for her work on Charles Babbage's mechanical computer, the Analytical Engine. Babbage (1791-1871) is considered a computer pioneer, because of the design of this first automatic computing engine.

In 1840, Giovanni Plana of University of Turin invited Babbage for a talk about his Analytical Engine, in occasion of the Second Congress of Italian Scientists invited by the King Carlo Alberto of Savoy, congress that was held at the Academy of Science. "La presentazione appassionò gli scienziati italiani e proseguì in seminari ristretti. Particolarmente interessati a questi seminari, nei quali per la prima volta si discusse di concatenamento delle operazioni, potremmo dire di programmazione, furono il fisico Mossotti e Luigi Menabrea" [3].

Luigi Menabrea, an engineer who became a Prime Minister of Italy, wrote an article upon the Babbage's lecture in French, and this transcript was published in the Bibliothèque Universelle de Genève in October 1842. Babbage's friend Charles Wheatstone commissioned Ada Lovelace to translate Menabrea's paper into English. She added some notes to the translation, with a discussion on the calculation of Bernoulli numbers [4,[START_REF] Menabrea | Sketch of The Analytical Engine Invented by Charles Babbage[END_REF]. So we have that the Bernoulli numbers had been recognized as important for calculations that we can obtain by concatenated operations.

In reading the Ref. [START_REF] Menabrea | Sketch of The Analytical Engine Invented by Charles Babbage[END_REF], we can see that Ada Lovelace proposed one of this operations, that, as we will see, is used for the Debye functions (see Figure 1). . The older definition of the Bernoulli numbers is given using equation:

(

3) ... ! 6 ! 4 ! 2 )! 2 ( ) 1 ( 1 2 1 6 * 3 4 * 2 2 * 1 1 2 * 1 + + - = - = - + - ∑ ∞ = - x B x B x B n x B x e x n n n n x or (4) ... ! 6 ! 4 ! 2 )! 2 ( 2 cot 2 1 6 * 3 4 * 2 2 * 1 1 2 * + + + = =       - ∑ ∞ = x B x B x B n x B x x n n n
In the modern use, the Bernoulli number s are given by the general of the generating function: The formula for the case n=0 is the case used by Ada Lovelace (see Figure 1).

As told in [START_REF]CRC Concise Encyclopedia of Mathematics, Second Edition[END_REF], the modern Bernoulli numbers appears in the calculation of the Debye functions, which are:

(8)         + + + - = - ∑ ∫ ∞ =1 2 2 0 ) ! 2 )( 2 ( ) 1 ( 2 1 1 k k k n x t n k n k x B n x n x e dt t
where |x|<2π and B n are the Bernoulli numbers.

Let us see when and where this modern use is mentioned for the first time.

Debye theory:

The first use of the Bernoulli numbers in the Debye theory of the specific heat was given by Hans Thirring in 1913, as told in [START_REF] Thirring | Zur Theorie der Raumgitterschwingungen und der spezifischen Wärme fester Körper[END_REF][START_REF] Blackman | On the relation of Debye theory and the lattice theory of specific heats[END_REF]. In solid state physics, the Debye model is a method developed by Peter Debye in 1912.

From the Debye model we can estimate the phonon contribution to the specific heat in a solid [START_REF] Debye | Zur Theorie der spezifischen Waerme[END_REF]. The model correctly predicts the low temperature dependence of the heat capacity, and, like the Einstein model, it is giving the Dulong-Petit law at high temperatures. As told in [START_REF] Blackman | On the relation of Debye theory and the lattice theory of specific heats[END_REF], for intermediate temperatures, Thirring proposed a specific heat given by:

(9)                 - - + = ∑ ∞ =1 2 ! ) 2 ( ) 1 2 ( ) 1 ( 1 n n n n kT hv n n B k c
where hν/kT<2π, h is the Planck constant and k the Boltzmann constant. B n are the Bernoulli numbers. In [START_REF] Blackman | On the relation of Debye theory and the lattice theory of specific heats[END_REF], ν is the frequency of a linear oscillator.

In the Ref.9, B 1 is equal to 1/6. In fact, it is giving the following for the heat capacity:

(10)           + - =           + - = ∑ ∑ K K ν ν ν ν 2 2 2 2 2 2 2 2 1 3 1 12 1 1 3 2 3 N T k h Nk T k h B N k C
It means that in [START_REF] Blackman | On the relation of Debye theory and the lattice theory of specific heats[END_REF], it was an old notation used for the Bernoulli numbers. Let us write Eq.10 with the notation of Ada Lovelace of these numbers:

(11)                 - - + = ∑ ∞ = - 1 2 1 2 ! 2 ) 1 2 ( ) 1 ( 1 n n A n n kT hv n n B k c
Probably, their use had been inspired to Thirring by the Lovelace's work.

In general, the expression of the heat capacity is: [START_REF] Fumi | On the Mie-Grüneisen and Hildebrand approximations to the equation of state of cubic solids[END_REF].

(12)                 - - + = ∑ ∞ = - n n n A n n kT h n n B k c 2 1 2 1 2 ! 2 ) 1 2 ( ) 1 ( 1 µ In (12), we have ∑ = ν ν µ n n
Appendix: Number e has a very interesting history [12]. In [12], it is told that the number e first comes into mathematics in a very minor way. This was in 1618 when, in an appendix to Napier's work on logarithms, a table appeared giving the natural logarithms of various numbers. A few years later, in 1624, Briggs gave a numerical approximation to the base 10 logarithm of e but did not mention this number in his work. The next possible occurrence of e is dubious, when, in 1647, Saint-Vincent computed the area under a rectangular hyperbola. Certainly by 1661, Huygens understood the relation between the rectangular hyperbola and the logarithm. Huygens made another advance in 1661. He defined a curve that he calls "logarithmic", that in today terminology is an exponential curve, having the form y = k a x . In 1668, Nicolaus Mercator published Logarithmotechnia which contains the series expansion of log(1+x). In this work Mercator uses the term "natural logarithm" for the first time for logarithms to base e. As told in [12], it is curious that, since the works on logarithms had come so close to recognizing the number e, that e was first "discovered" not through the notion of logarithm but rather through a study of compound interest, as did in 1683 Jacob Bernoulli looked at the problem of compound interest [12].

Figure 1 :

 1 Figure 1: Snapshot from Ref.5.


Yassine Hacha
email: yassine.hachaichi@ipeiem.rnu.tn

Younes Lahbib
email: younes.lahbib@enicarthage.rnu.tn

In order to approximate transandental functions, several algorithms were proposed. Historically, polynomial interpolation, infinite series, • • • and other +, ×, -and / based algorithms were studied for this purpose.

The CORDIC (COordinate Rotation DIgital Computer) introduced by Jack E. Volder in 1959, and generalized by J. S. Walther a few years later, is a hardware based algorithm for the approximation of trigonometric, hyperbolic and logarithmic functions.

As a consequence, CORDIC is used for applications in diverse areas such as signal and image processing. For these reasons, several modified versions were proposed.

In this article, we present an overview of the CORDIC algorithm for the computation of the circular functions, essentially the scaling free version, and we will give a substential improvement to the commonly used one.

Introduction

In 1959, Volder [START_REF] Volder | The CORDIC trigonometric computing technique[END_REF], introduced the CORDIC algorithm in order to compute approximations of trigonometric functions. This method is still used because of its adequacy to hardware design. It is a recursive method using only shiftand-add operations.

A decade later, J. S. Walther in [START_REF] Walther | A unified algorithm for elementary functions[END_REF], generalized this method to other transendantal functions used in engineering fields.

The development of the CORDIC algorithm and architectures [START_REF] Lakshmi | CORDIC Architectures: A Survey[END_REF] has taken place for achieving the highest throughput rate and reduction of hardwarecomplexity as well as the computational latency of implementation. Some of the typical approaches for reducing complexity implementation are targeted on minimization of using the scaling-operation and complexity of barrel-shifters and adders in the CORDIC engine. However, one of the problems associated with the classical CORDIC formulation is that the scale factor depends of the angle, and is not constant. The complexity of the computation of the scale factor is in principle comparable to that of the basic CORDIC process itself. In a recent work, [START_REF] Garrido | CORDIC II: A New Improved CORDIC Algorithm[END_REF], a new algorithm, CORDIC II, is proposed that substitutes the CORDIC micro-rotation by a new angle set.

Aiming to eliminate scale multiplication in conventional CORDIC, scale free CORDIC was used to eliminate the scale factor, see the piooneering papres [START_REF] Aggarwal | Area-time efficient scaling-free cordic using generalized micro-rotation selection. Very Large Scale Integration[END_REF][START_REF] Aggarwal | Hardware efficient architecture for generating sine/cosine waves[END_REF] and also [START_REF] Mokhtar | Scaling free CORDIC algorithm implementation of sine and cosine function[END_REF][START_REF] Maharatna | Virtually scalingfree adaptive CORDIC rotator[END_REF][START_REF] Maharatna | Modified virtually scaling-free adaptive CORDIC rotator algorithm and architecture[END_REF]. The scale free CORDIC algorithm for cosine and sine functions is proved to be faster and efficient in terms of area and accuracy compared to conventional CORDIC.

We give in this paper a method in order to minimize the number of iterations ine the CORDIC method. This is given by computing the closest elementary angle to the residual one at each iteration. Our second contribution is the correction of the Taylor series used for the composed functions. We will prove that with our polynomial approximation, we will get faster computation for the same acuity. The CORDIC algorithm operates either in, rotation mode or vectoring mode, following linear, circular or hyperbolic coordinate trajectories. In this paper, we focus on rotation mode CORDIC using circular trajectories.

2 The CORDIC algorithm.

The idea behind conventional CORDIC algorithm is the rotation of a vector [x in y in ] T in cartesian coordinate which can be expressed in [START_REF] Aggarwal | Hardware efficient architecture for generating sine/cosine waves[END_REF], where [x out y out ] T is the output vector produced after rotation and θ is the angle of rotation.

x out y out = cos(θ) -sin(θ) sin(θ) cos(θ)

x in y in (1)
This can also be written as

x out y out = cos(θ) 1 -tan(θ) tan(θ) 1 x in y in (2)
We split the rotation angle in a sum of angles, and carries out the rotation by a series of the so called micro-rotation by these angles. The idea is to decompose any angle θ into a sum of some "elemntary" angles

θ = α 1 + • • • + α n (3)
where

α k = ± arctan(2 -l ).
If we use the fact that, if R θ denotes the matrix of the 2D rotation of angle θ:

R θ × R θ ′ = R θ+θ ′
We can translate the equation (3) into the matrix product :

R θ = R α1 × • • • × R αn 2.1 The conventional CORDIC
The conventional CORDIC method performs a sequence of rotations by elementary angles. Any rotation θ on the plan can be decomposed into a composition (matrix product) of n elementary rotations. When taking θ k = arctan(2 -k ), the equation (2) becomes:

R θ k = 1 √ 1 + 2 2k 1 -2 -k 2 -k 1
Using also the identity

R -θ = cos(θ) 1 tan(θ) -tan(θ) 1 we obtain, for ε k = ±1, R ε k θ k = 1 √ 1 + 2 2k 1 -ε k 2 -k ε k 2 -k 1
The idea is that the angles used are constant, so we have a constant scale

K = n k=1 1 √ 1+2 2k
, which approximately equals, according to the litterature [START_REF] Surapong | Design and Analysis of Extension-Rotation CORDIC Algorithms based on Non-Redundant Method[END_REF], 0.60725. For this aim, we construct a sequence of vectors [x k y k z k ] T by the recurrence schema:

  x 0 y 0 z 0   =   x in y in θ   (4)
and   

x k+1 = x k -ε k 2 -k y k y k+1 = ε k 2 -k x k + y k z k+1 = z k -ε k θ k (5)
After the fixed number of iterations, we mutiply the resulting vector by the constant K, this means [

x out y out ] = K[x n y n ].
The essence of the CORDIC algorithm is that he is multiplication free (only shift-and-add operations). The scale multiplication, also called compensation, to get an output vector isometric to the input one, causes a problem.

The introduction of the scale free CORDIC is then legitimated.

The correct scale free CORDIC for sine and cosine

The scale free CORDIC for circular functions is based on the Taylor series

sin(x) = ∞ n=0 (-1) n • x 2n+1 (2n + 1)! = x - x 3 6 + x 5 120 + • • • cos(x) = ∞ n=0 (-1) n • x 2n (2n)! = 1 - x 2 2 + x 4 24 + • • • arctan(x) = ∞ n=0 (-1) n • x 2n+1 2n + 1 = x - x 3 3 + x 5 5 + • • •
The approximation Taylor polynomial of the composed functions to order 5:

sin(arctan(x)) ≈ x - x 3 3 + x 5 5 - 1 6 x - x 3 3 + x 5 5 3 + + 1 120 x - x 3 3 + x 5 5 5 cos(arctan(x)) ≈ 1 - 1 2 x - x 3 3 + x 5 5 2 + 1 24 x - x 3 3 + x 5 5
4

When we truncate the polynomials to the order 5, we obtain the right equations:

cos(arctan(x)) ≈ 1 - 1 2 x 2 + 3 8 x 4 (6) sin(arctan
(x)) ≈ x - 1 2 x 3 + 3 8 x 5 (7)
We can observe that sin(arctan(x)) = x cos(arctan(x)), this is simply due to the fact that sin(arctan(x)) cos(arctan

(x)) = tan(arctan(x)) = x (8)
We the obtain, for the elementary angles θ k = arctan(2 -k ), and remarking

that 3 8 = 1 4 + 1 8 sin(arctan(2 -k )) ≈ 2 -k -2 -3k-1 + 2 -5k-2 + 2 -5k-3 cos(arctan(2 -k )) ≈ 1 -2 -2k-1 + 2 -4k-2 + 2 -4k-3 The rotation matrix M θ k becomes: 1 -2 -2k-1 + 3 • 2 -4k-3 -2 -k + 2 -3k-1 -3 • 2 -5k-3 2 -k -2 -3k-1 + 3 • 2 -5k-3 1 -2 -2k-1 + 3 • 2 -4k-3 (9)
As we know, all the works we have seen uses the Taylor series for sine and cosine functions and replace θ k = arctan(2 -k ) by 2 -k , see [START_REF] Aggarwal | Area-time efficient scaling-free cordic using generalized micro-rotation selection. Very Large Scale Integration[END_REF][START_REF] Aggarwal | Hardware efficient architecture for generating sine/cosine waves[END_REF] and also [START_REF] Mokhtar | Scaling free CORDIC algorithm implementation of sine and cosine function[END_REF][START_REF] Maharatna | Virtually scalingfree adaptive CORDIC rotator[END_REF][START_REF] Maharatna | Modified virtually scaling-free adaptive CORDIC rotator algorithm and architecture[END_REF]. The error is that when using a Taylor polynomial of a composite function f • g, we have to use the same degree and truncate the resulting polynomial at the demanded degree, you can see [START_REF] Monier | Analyse MPSI: cours et 1000 exercices corrigés[END_REF].

In order to give an empirical proof, we will compare the orders 3, 4 and 5 of ou method to the recent works [START_REF] Mokhtar | Scaling free CORDIC algorithm implementation of sine and cosine function[END_REF][START_REF] Aggarwal | Hardware efficient architecture for generating sine/cosine waves[END_REF].

Benchmark of scale free CORDIC for circular functions

In order to minimize the number of iterations of the CORDIC algorithm, we choose the microrotations to be the closest arctan(2 -k ) to the residual angle. This can be done by choosing the closest power of 2 to tan(θ), where θ represent the risidual angle.

Due of the continuity of the function arctan, if tan(θ) is close to 2 -k , then so is θ to θ k . This leads us to choose k such that, the closest θ k = arctan(2 -k ) to θ the following way:

k = Round log 2 1 |θ| (10)
we replace arctan(θ) by θ without any loss of acuity because θ is very close to arctan(θ) for θ in [0, π 4 ]. As an example:

π 16 = θ 1 -θ 4 -θ 7 -θ 10 + θ 12 (±10 -5 )
For a hardware design, the translation of our method for the binary represen-

tation |θ| = 0.ε 1 • • • ε i • • • where ε i ∈ {0, 1}
if for i ≥ 1 we have ∀j < i; ε j = 0 and ε i = 1 and ε i+1 = 0 then k = i if for i ≥ 1 we have ∀j < i; ε j = 0 and

ε i = ε i+1 = 1 then k = i -1
In this section, we will compare our approximation to the one given in [START_REF] Aggarwal | Hardware efficient architecture for generating sine/cosine waves[END_REF] and [START_REF] Mokhtar | Scaling free CORDIC algorithm implementation of sine and cosine function[END_REF] for order 3 Taylor approximations. The range of angles used is [0, π 4 ]. This range is enough, using the trigonometric identities, to can calculate any sine or cosine of any angle.

In [START_REF] Aggarwal | Hardware efficient architecture for generating sine/cosine waves[END_REF], cos(arctan(2 -k )) is approximated by 1-2 -2k-1 , and sin(arctan(2 -k )) is approximated by 2 -k -2 -3k-3 . The authors in [START_REF] Mokhtar | Scaling free CORDIC algorithm implementation of sine and cosine function[END_REF], use the approximation cos(arctan(2 -k )) ≈ 1 -2 -2k-1 , and sin(arctan

(2 -k )) ≈ 2 -k -2 -3k-2 .
The proposed method use the approximation cos(arctan(2 -k )) ≈ 1 -2 -2k-1 , and sin(arctan(2 -k )) ≈ 2 -k -2 -3k-1 , which is the mathematically correct developpement.

The quadratic errors for cosine and sine function for the three method are summarized in the tables (1 and 2).

Method in [START_REF] Aggarwal | Hardware efficient architecture for generating sine/cosine waves[END_REF] Method in [START_REF] Mokhtar | Scaling free CORDIC algorithm implementation of sine and cosine function[END_REF] In the table (3) below, we compare the quadratic errors of our method in different order Taylor approximation.

In figures (2 and 1), a MatLab simulation of our method is given. In blue, the graph of our method, in red, the graph of the matlab function and in green the difference between them. Method in [START_REF] Aggarwal | Hardware efficient architecture for generating sine/cosine waves[END_REF] Method in [START_REF] Mokhtar | Scaling free CORDIC algorithm implementation of sine and cosine function[END_REF] A new design of the main computation unit is proposed in this paper and compared to the conventional CORDIC one. This is mainly a study to check if the theoretical results are feasible and simple to embed. Optimizations, complete CORDIC computation schemes, advanced CORDIC architectures and comparisons, which are based on the underlined computing unit, are future works.

The proposed scale-free CORDIC algorithm is based on Taylor polynomials. Three orders are evaluated for benchmarking the theoretical study (see Table 3). As the complexity of hardware architecture grows in function of the development order, only the order 3 is implemented in the hardware side. However, the impact of Taylor's order on hardware performances is also ongoing. It is not the focus of this paper.

Hence, the implemented hardware architecture is restricted to the main computation unit and the order 3 of Taylor series. It is composed of 4 blocks: dynamic index predictor, shifting processor, storing angles ROM and FSM controller. Figure 3 gives general description of the architecture. A detailed description is presented in the following sub-sections. The section ends with a summary of the main hardware results.

FSM controller

The controller is a finite state machine with three states. In the initial state, the signal load is set to initiate the initial values of the CORDIC core, namely X 0 = 1, Y 0 = 0 and the angle θ. The second state is processed 2 N cycles according to an N -bits counter which fixes the number of CORDIC iterations.

For a given iteration, new intermediate values X n and Y n are obtained by shifting previous X n-1 and Y n-1 according to the closest micro-rotation as given in equation 10. The third and final state sets the done signal. Cosine and sine of the angle θ are computed and stored in the output registers.

Dynamic Index predictor

The main theoretical result proved in section 2. (see theoretical result) is implemented in the Index Predictor. The computation of the next index is the main improvement of the proposed hardware architecture. It estimates the optimal index with which we address the ROM and read the closest CORDIC microrotation for a given iteration. The determined angle is compared with the ongoing error Z in order to compute the new direction of the micro-rotation. X n and Y n are then shifted index positions to right. The next listing describes the behavior of this block. The hardware implementation is based on a 32 bits comparator. When the first sequence of '11' bits are detected the block returns the corresponding K-index. Otherwise the first sequence of '10' is cheeked, but the (K -1)-index is returned in this case. In this way, the most significant power of 2 micro-rotation is obtained.

IndexEstim

Synthesis of the iterative CORDIC design

The proposed iterative design based on the described computation unit is implemented on the Nexys3 spartan6 FPGA. Compared to the conventional architecture we reach the same cosine and sine results for only 4 iterations versus 32 iterations in the conventional case. The precision error is about 10 -3 in both cases. Hardware performances such as power, frequency and area are summarized in table 4. Implementation results show a degradation of the In this paper two improvements were made. First, we have minimized the number of iterations for some fixed error by calculating the closest elementary angle to the residual one. Second, we gave the correct polynomial approximation for the scale free CORDIC. The comparison between our method and two other famous methods is given to confirm empirically our theoretical proof. In our simulation, we remark that the order of approximation of Taylor series used meets the accuracy requirements.

In section 3, we showed that these methods have a simple hardware implementation, in order to meet the objectives of the CORDIC's introduction. The iterative and pipelined architectures were implemented, and significant improvements of hardware performance were denoted in the pipelined case. The future works will focus on the improvement of the hardware architecture.

Fig. 1 Fig. 2

 12 Fig.1sin approximation (in green the difference)

Fig. 3

 3 Fig. 3 Design of the CORDIC architecture.

Fig. 4

 4 Fig. 4 Design of the pipelined CORDIC architecture.

Table 1

 1 Comparison of methods for cosine function

	Our method

Table 2

 2 Comparison of methods for sine function4 Hardware implementationCommon Hardware implementations of CORDIC algorithms are either iterative or pipelined[START_REF] Valls | Evaluation of CORDIC Algorithms for FPGA Design[END_REF][START_REF] Antonius | FPGA implementation of CORDIC algorithms for sine and cosine generator[END_REF]. The main computation CORDIC unit is iterated in both cases. It is unrolled in the first class and rolled in the former using pipelined registers to store intermediate computations[START_REF] Dian-Marie | Design rules for implementing CORDIC on FPGAs[END_REF][START_REF] Ray | A survey of CORDIC algorithms for FPGA based computers[END_REF].

	Our method

Table 3

 3 Comparison of different orders of Taylor polynomials

  -k are stored in a ROM as constants, for k within the range 0 to 31. The ROM values are not addressed inclemently as done in the conventional implementation of CORDIC. The closest 4 microrotations are read instead of 32 in the case of conventional CORDIC; The 4 indexes are estimated dynamically by the Index Predictor Block.Listing below shows the addressing behavior the stored arctangent values.X and Y intermediate signals are computed as presented in listing below. This corresponds to processing equations in section 2. They are 32-bit coded signals in 2-complement format. Control signals sgn and Index are dynamically estimated and set by the Index predictor block. Note that listing below gives the behavior computing of the signal Y , it is similar in the case of X.The VHDL code behavior of the corresponding processing is similar to the listing below.

	: process(Zerror) for index in 31 downto 1 loop begin if(Zerror (index)='1' and Zerror (index-1) ='1') then Inxt <= std_logic_vector(to_unsigned(31-index,32)); exit; elsif(Zerror(index) ='1' and Zerror(index-1) ='0') then Inxt<= std_logic_vector(to_unsigned(31-(index-1),4)); exit; else Inxt <= (others => '0'); end if; end loop; end process IndexEstim; 4.3 Storing angles ROM begin case Index is when "00000" => MicroRotI <= X"3243f6a8"; _____ when "11000" => MicroRotI <= X"0000003f"; _____ when "11111" => MicroRotI <= X"00000000"; end Case; end process IndexAccess; 4.4 Shifting processor sinShift:process(Xsgn) begin Yshift_tmp((31-I) downto 0) <= Xsgn(31 downto I); if Xsgn(31) = '0' then Yshift_tmp(31 downto (32-I)) <= (others => '0'); else Yshift_tmp(31 downto (32-I)) <= (others => '1'); end if; Radix values of arctangent of 2 IndexAccess: process(INDEX) end process sinShift;

hardware performances of the proposed architecture. By an in-depth Analysis of sub-blocks, we find that the Index predictor block consumes alone 20 mw and uses 240 slices. The shifting processor block on the other side consumes also 20 mw and uses 540 slices.

The index predictor is an extra block in our case which explains the extra values against the conventional architecture. However, we think that the main reason is our coding style of the VHDL design which was behavioral. The behavioral synthesis infers the use of LUTs rather than basic logics. Hence, more optimized implementation should lead to less logic slices.

Synthesis of the pipelined CORDIC design

A pipelined CORDIC architecture consists of rolling the main computation unit by storing intermediate computation into registers. In the case of the conventional architecture the main unit is rolled 32 times when processed data is 32-bit coded. The main improvement of our proposed architecture is rolling the same unit 3 or 4 times whatever the size of the processed data. More rolled units can be implemented if more precision is needed. With 3 units a precision of 10 -2 is reached and 10 -3 when 4 units are used. Figure 4 shows the proposed pipelined architecture. The index predictor which is resources consuming is instantiated only one time. The FSM controller enables the communication with only one pipeline stage. As shown in table 5 significant results compared to the iterative architecture are obtained. We save almost 50% of area and power with a speedup of 10 Mhz, against only 16-bit conventional architecture.

Conclusion

CORDIC algorithm has several applications in several domains, for an the reader can read [START_REF] Meher | 50 years of CORDIC: Algorithms, architectures, and applications[END_REF]. The popularity of this method is due to the simplicity of its hardware implementation, see [START_REF] Lakshmi | CORDIC Architectures: A Survey[END_REF] for example.


Imen Ben Saad
email: bensaad_imen1@yahoo.fr

Younes Lahbib
email: younes.lahbib@enicarthage.rnu.tn

Yassine Hachaïchi
email: yassine.hachaichi@ipeiem.rnu.tn

Sonia Mami
email: soniammkd@gmail.com

Abdelkader Mami
email: abdelkader.mami@fst.rnu.tn

Generic-Precision algorithm for DCT-Cordic architectures

In this paper we propose a generic algorithm to calculate the rotation parameters of CORDIC angles required for the Discrete Cosine Transform algorithm (DCT). This leads us to increase the precision of calculation meeting any accuracy.Our contribution is to use this decomposition in CORDIC based DCT which is appropriate for domains which require high quality and top precision.

We then propose a hardware implementation of the novel transformation, and as expected, a substantial improvement in PSNR quality is found.

Introduction

The Discrete Cosine Transform DCT was developed by Ahmed et.al in 1974 [START_REF] Ahmed | Discrete cosine transform[END_REF]. It is a robust approximation of the optimal Karhunen-Loeve Transform (KLT) [START_REF] Rao | Discrete Cosine Transform Algorithms, Advantages, Applications[END_REF]. It has become one of the most widely used techniques of transforms in digital signal processing. The DCT is one of the computationally intensive transforms. It requires many multiplication and additions. Many researches had been done on low-power DCT designs [START_REF] August | Low power design of DCT and IDCT for low bit rate video codecs[END_REF][START_REF] Jeong | Low-power multiplierless DCT architecture using image correlation[END_REF]. In consideration of VLSI-implementation, Flow-Graph Algorithm (FGA) is the most popular way to realize the fast DCT (FDCT) [START_REF] Chen | A Fast Computational Algorithm for the Discrete Cosine Transform[END_REF][START_REF] Wang | Fast algorithms for the discrete W transform and for the discrete Fourier transform[END_REF].

Most existing works handled the reduction of energy consumption of the DCT. As the multiplications are energy expensive operations, several algorithms are based on additions and shifts instead of multiplications.

In 2004, Jeong et .al [START_REF] Jeong | Low-power multiplierless DCT architecture using image correlation[END_REF] suggested improving a Cordic based implementation of the DCT. COordinate Rotation Digital Computer (CORDIC) is an algorithm which can be used to evaluate various functions in signal processing [START_REF] Volder | The CORDIC trigonometric computing technique[END_REF][START_REF] Walther | A unified algorithm for elementary functions[END_REF]. In [START_REF] Jeong | Low-power multiplierless DCT architecture using image correlation[END_REF], authors proposed a low-complexity CORDIC based DCT algorithm based on FGA.It requires only 38 add and 16 shift operations and consumes about26.1 % less power compared to [START_REF] Mariatos | A fast DCT processor, based on special purpose CORDIC rotators[END_REF],with a minor image quality degradation of 0.04 dB.

In the same direction, Sun et .al [START_REF] Sun | Low-power and high quality Cordic-based Loeffler DCT for signal processing[END_REF][START_REF] Sun | VLSI implementation of a configurable IP Core for quantized discrete cosine and integer transforms[END_REF] proposed a new flow graph for Cordic based Loeffler DCT implementation. A new table of parameters is obtained with new choice of the elementary rotations. However, Authors did not give any demonstrations or a method to compute the underlying elementary rotations.Their experimental result shows that the Cordic-based Loeffler DCT consumes 16%of energy compared to [START_REF] Loeffler | Practical fast 1-D DCT algorithms with 11-multiplications[END_REF] with a minor image quality degradation of 0.03 dB.

After this analysis of state of the art, we remark that previous works chose to reduce precision of calculation to decrease the energy consumption.Themain disadvantage of these methods is the lack of precision, which restricted their applications in Medical, Biometrics and High frequency for example.

Despite of 0.01w power loss, our proposed DCT architecture reaches an improved image quality of 6.55 dB PSNR compared to [START_REF] Sun | Low-power and high quality Cordic-based Loeffler DCT for signal processing[END_REF][START_REF] Sun | VLSI implementation of a configurable IP Core for quantized discrete cosine and integer transforms[END_REF].This gain is very significant compared to PSNR evaluation of previous works, (see TABLE III in [START_REF] Lee | Reconfigurable CORDIC-Based Low-Power DCT Architecture Based on Data Priority[END_REF]). For thisreason, we proposeincreasing the precision degree ofDCT-Cordic micro-rotations with an exact and generic method.According to a chosen precision degree, exact DCT parameters are computed using the proposed algorithm.

This paper is organized as follows. Section 2 briefly introduces the algorithms of conventional Cordic-Based DCT Architecture. In Section 3, we present the solution approach of the high precision Cordic-Rotation. The experimental results are shown in Section 4 while Section 5 concludes this paper.

Conventional Cordic-Based DCT Architecture

Cordic Algorithm

The conventional Cordic algorithm [START_REF] Volder | The CORDIC trigonometric computing technique[END_REF] is hardware-efficient used for the approximation computation of the transcendental functions. It only uses shift and addition operations. The Cordic algorithm can operate in two modes, namely vectoring and rotation. In this paper we focus on the first mode.

In the conventional Cordic algorithm, a rotation angle is decomposedinto combination of micro-rotation angles of arctangent radix. When the vector is rotated by an angle θ, the coordinate changed from (𝑋 𝑖 ,𝑌 𝑖 ) to (𝑋 𝑖+1 ,𝑌 𝑖+1 ).

The value of vector after this micro rotation can be represented as:

𝑥 𝑖+1 𝑦 𝑖+1 = 𝐾 𝑖 1 -𝜎 𝑖 tan 𝜃 𝑖 𝜎 𝑖 tan 𝜃 𝑖 1 𝑥 𝑖 𝑦 𝑖 ( 1
)
where𝐾 𝑖 =𝑐𝑜𝑠 𝜃 𝑖 and 𝜃 𝑖 =𝑡𝑎𝑛 -1 (2 -𝑖 )
The circular rotation angle is depicted as:

= 𝜎 𝑖 . 𝑖 𝜃 𝑖 , where 𝜎 𝑖 =±1 (2)
Then, the vector rotation can be performed iteratively as follows:

𝑥 𝑖+1 𝑦 𝑖+1 = 𝐾 𝑖 1 -𝜎 𝑖 2 -𝑖 𝜎 𝑖 2 -𝑖 1 𝑥 𝑖 𝑦 𝑖 (3)
Fig. 1.The direct implementation of equation 3.

In the equation 3, only shift and add operations are required to perform the rotation angle described in Fig. 1. But, the results of the rotation iterations need to be scaled by a compensation factor K. This can be done by using the following iterative method.

K= 𝑘 𝑖 𝑖 = 1/ 1 + 2 -2𝑖 𝑖 (4)
The scale factor K which can be interpreted as a constant gain (hence not data dependent) can be tolerated in many digital signal processing applications. Hence, it should be carefully investigated whether it is necessary to compensate for the scaling at all. If scale factor correction cannot be avoided, two possibilities are known. The first approach consists on performing a constant factor multiplication with 1/𝐾 𝑖 . The second method is based on extending the Cordic iteration in a way that the resulting inverse of the scale factor takes a value. In other words, we need to write the scaling factor as a sum of 2 -i where i must be determined so that the error is minimized. In the rotation mode, the angle accumulator is initialized with the desired rotation angle. The rotation decision at each iteration is made to diminish the magnitude of the resi-dual angle in the accumulator one. The decision at each iteration is therefore based on the sign of the residual angle after each step [START_REF] Meher | 50 years of CORDIC: Algorithms, architectures and applications[END_REF].

Cordic-Based DCT Architecture

The One-dimensional DCT for 8x8 sub-images is defined as

F(k)= 1 2 C(k) 𝑓 𝑥 cos [ 2𝑥+1 𝑘𝜋 16 ] 7 𝑥=0 C(k)= 1 2 𝑖𝑓 𝑘 = 0 1 𝑜𝑡𝑒𝑟𝑤𝑖𝑠𝑒 (5)
Where x (i) is the input data and X (t) is 1-D DCT transformed output data. The two-dimensional DCT is a separable transform. It can be executed by onedimensional DCT in a serial manner as shown in the Fig. 2. .

The rearranged 1-D DCT equation is now represented as a vector rotation matrix together with the consecutive CORDIC iterations as shown in Fig. 3. The Cordic array performs the fixed-angle rotation in the DCT algorithm. Therefore, the general signal flow graph of Cordic-based DCT is presented by Fig. 4.

Proposed High precision Cordic-Based Loeffler DCT Architecture

In this section we will present the proposed algorithm which calculates the Cordic Rotations. The main result of this algorithm is enhancing the degree of precision by improving the selected parameters in order to find the exact values of the rotations.

Computation of Micro-Rotation decomposition

The proposed algorithm takes as input the rotation angle and the precision degree. The choice of this last parameter will depend on the application field. For instance, medical or biometrics applications require high precisions whereas compression process or video surveillance can be treated with low a greater tolerance error. This algorithm provides the Cordic parameters (iterations and direction) corresponding to the angle and the selected precision. The iterations, in other words, the micro-rotations are identified by the algorithm with their orientation, clockwise or anticlockwise.

Algorithm obtains the optimal parameter of micro-rotations. 𝐼𝑛𝑝𝑢𝑡 ∶ 𝜃 = 𝑎𝑛𝑔𝑙𝑒 𝑎𝑛𝑑 𝜀 = 𝑝𝑟𝑒𝑐𝑖𝑠𝑖𝑜𝑛

When applying our algorithm we obtain the new parameters table which depends on the required precision. (See Table 1)

> > 2 > > 2 - + > > 4 > > 4 + - > > 6 > > 6 - + > > 9 > > 9 + - > > 1 3 > > 1 3 - + x y X Y
In order to check the correctness of the computed rotations matrices according to the proposed algorithm, obtained values are compared to ideal ones. Like seen in the equation 3 and since we have 5 iterations, the Cordic of the angle  = 𝜋 16 is obtained by multiplying 5 rotation matrices: 13 1

𝑋 𝑌 = 𝐾 1 -2 -2 2 -2 1 1 2 -4 -2 -4 1 1 -2 -6 2 -6 1 1 2 -9 -2 -9 1 1 -2 -13 2 -

𝑥 𝑦

Then we get the following matrix coefficients according to the proposed algorithm (equation 8

As it can be denoted, the ideal rotation matrix (equation 8) has exactly the same coefficients of the computed rotation matrix (equation 9), for an approximation of 10 -3

Experimental Results

Our algorithm has been implemented on Virtex5 xc5vlx30-3ff676 using Xilinx System Generator (XSG). The power consumption of our DCT architecture with different precision degrees is shown in the In order to demonstrate the high-quality feature of the proposed DCT algorithm, it has been evaluated considering a JPEG2000 compression chain [14].

Table 3 demonstrates the comparison of the average PSNR of the proposed DCT algorithm in tow precision degrees, with the CORDIC based Loeffler DCT [START_REF] Sun | Low-power and high quality Cordic-based Loeffler DCT for signal processing[END_REF][START_REF] Sun | VLSI implementation of a configurable IP Core for quantized discrete cosine and integer transforms[END_REF]. Checked results considerhigh-to-low quality compression (i.e. quantization factors from 95 to 75) using some well-known test images. For instance, Fig. 8 gives the experimental results based on Lena image.

It can be easily noticed from Table 2 that the performance of the proposed DCT algorithm has better quality about 6.55 for Q=95 dB than the Cordic-based Loeffler for the JPEG2000 standard, without substantial loss of power.

Conclusion

In this paper, we proposed a high quality algorithm which calculates the rotation parameters of CORDIC angles required for the Cordic-based Loeffler DCT algorithm.

The obtained results show a significant improvement in the PSNR (6, 55 dB for Q=95) without a substantial loss of Power. These enhancements make our algorithm appropriate for high accuracy applications such as medical or biometrics domains. The proposed algorithm calculates the rotation parameters of CORDIC for different precision degrees from 10 -3 up to 10 -7 . The precision degree results higher than 10 -4 they don't give a great improvement in terms of PSNR with respect to the precision 10 -3 . In fact, for Q=95, the PSNR obtained for a precision of 10 -6 is equal to 43.53831 dB, which means that the variation is in the order of 10 -4 dB. This improvement is small especially that it requires the addition of some layers of rotations in the CORDIC architectures. Hence, it is useless to go higher than a precision of 10 -3 .

Fig. 2 .

 2 Fig.2. 8 x 8 2-D DCT processor with separable 1-D DCT The 1-D DCT transform is represented as follows.

Fig. 3 .

 3 Fig. 3. Hardware architecture of CORDIC-based 1-D DCT[START_REF] Lee | Reconfigurable CORDIC-Based Low-Power DCT Architecture Based on Data Priority[END_REF]

Fig. 4 .

 4 Fig. 4. The general signal flow graph CORDIC-based DCT

1 :

 1 while 𝜃 > 𝜀 do 2: 𝑥 = 𝑡𝑎𝑛 |𝜃| 3: 𝑖 = -𝑙𝑜𝑔(𝑥) + 1 4: 𝑖𝑓 𝜃 > 𝑡𝑎𝑛 -1 (2 -𝑖 ) then 5:  = + 6: 𝜃 = 𝜃 -𝑡𝑎𝑛 -1 (2 -𝑖 ) 7: else 8:𝜎 = -9: 𝜃 = 𝜃 + 𝑡𝑎𝑛 -1 (2 -𝑖 ) 10: 𝑤𝑟𝑖𝑡𝑒 𝑖 𝑎𝑛𝑑  11: 𝑒𝑛𝑑 𝑤𝑖𝑙𝑒

=

  𝜃 2 -𝜃 4 + 𝜃 6 -𝜃 9 + 𝜃 13 = 0.196349 So for the precision P = 10 -4 and since the real value of  = 0.1963 we can conclude that we have a zero-error approximation in the angle. Based on the previous computed micro-rotations of the π 16 angle, the Cordic architecture computing π 16 angle is given in Fig. 5.

Fig. 5 .

 5 Fig. 5. Unfolded flow graph of the π 16 angle

Fig. 8 .

 8 Fig. 8.Lena images obtained using the proposed Cordic-based Loeffler DCT for P =10 -3 . (a) Q=75. (b) Q=85. (c) Q=95.(d) Lena image obtained using[START_REF] Sun | Low-power and high quality Cordic-based Loeffler DCT for signal processing[END_REF][START_REF] Sun | VLSI implementation of a configurable IP Core for quantized discrete cosine and integer transforms[END_REF]

Table 1 .

 1 Required Iterations and Directions for Vector Rotation

	Angle	𝜽 =	𝝅 𝟒	𝜽 =	𝟑𝝅 𝟖	𝜽 =	𝝅 𝟏𝟔	𝜽 =	𝟑𝝅 𝟏𝟔
	Precision
		i=0		i=0/1/4/7	i=2/4/6/9	i= 1/3/10
	P=10 -3	𝜎 𝑖 = -		𝜎 𝑖 =-+ --	𝜎 𝑖 =+ -+ -	𝜎 𝑖 =+ + +
		i=0		i=0/1/4/7/10/12	i=2/4/6/9/13	i= 1/3/10
	P=10 -4	𝜎 𝑖 = -		𝜎 𝑖 =+ + ---+	𝜎 𝑖 =+ -+ -+	𝜎 𝑖 =+ + +
	Let's take an example and show its implementation. (See Table 2)

Table 2 .

 2 Rotator parameter of the

				π 16	angle for precisionP = 10 -4
		Input			Output
		𝑷	i	2	4	6	9	13
	𝜋 16	10 -4	 i	+	-	+	-	+
	So we can write the rotation angle as the weighted sum of micro-rotations i as
	shown in the equation 2 and 7, we obtain:

  ).

	𝑋 𝑌	= 𝐾	1.013067933963612 -0.2015148886130191 0.2015148886130191 1.013067933963612	𝑥 𝑦 (8)
		The ideal computed coefficients of	𝜋 16	rotation matrix are depicted in equation 9.
	𝑋 𝑌	= 𝐾	1.014359843600569 -0.201768717865449 0.201768717865449 1.014359843600569	𝑥 𝑦

Table 3 .

 3 The power consumption is measured with Xpower Analyzer with 100 Mhz clock cycles, 1V supply power.

	(7)

Table 3 .

 3 Hardware Implementation and PSNR from high-to-low compression quality in JPEG2000 for Lena image

		Power (W)	Quality factor	PSNR
			Q=95	36.982
	[10, 11]	0.642	Q=90 Q=85	36.019 35.113
			Q=80	34.295
			Q=75	33.708
			Q=95	43.535
	Proposed algorithm (P=𝟏𝟎 -𝟑 )	0.656	Q=90 Q=85	39.458 36.992
			Q=80	35.355
			Q=75	34.280
			Q=95	43.538
	Proposed algorithm	0.659	Q=90	39.459
	(P=𝟏𝟎 -𝟒 )		Q=85	36.994
			Q=80	35.355
			Q=75	34.281


Ezzahra Fatma

Salem

Abdoulaye Tall
email: abdoulaye.tall@orange.com

Zwi Altman
email: zwi.altman@orange.com

Azeddine Gati
email: azeddine.gati@orange.com

Energy consumption optimization in 5G networks using multilevel beamforming and large scale antenna systems

Keywords: Energy consumption, Green cellular networks, base station power model, multilevel codebook, Green policy management Throughput limited region Power limited region

Cellular networks are witnessing an exponential traffic growth leading to an increase in Energy Consumption (EC), and having both environmental and economic impact. Recently, different approaches have been studied to build Green cellular networks focusing mainly on the Base Stations (BSs) as the access network represents 80% of the total wireless network consumption [1] [2]. One of the promising solutions for increasing throughput and reducing EC is the deployment of large antenna arrays, known as Large Scale Antenna Systems (LSAS) [3] [4] [5] [6], that can transmit highly focused beams. In the present work, we focus on the use of an advanced BS power model developed within the GreenTouch project [7] [8]. This model allows to quantify the power consumption of a reference scenario comprising multiple sites with standard BS antennas, and then compare it to a LSAS solution implementing multilevel beamforming. We then exploit the LSAS merits to get greener deployment with less BSs, some of which can be turned off as a function of the traffic demand. The coverage areas of each cell can be modified, by updating the codebook of beams. We finally investigate different network configurations which represent distinct trade-offs between EC and capacity. We propose a methodology to represent and design green policies for managing the network which select the desired operating points. Detailed simulation results illustrate the proposed methodology.

I. INTRODUCTION

Nowadays, cellular networks are facing a tremendous increase in traffic due to the multiplicity of mobile devices and services. Obviously, this growth is accompanied by a high increase in the EC which has a major impact on the overall operational expenditures (OPEX) of the network operators and leads to higher CO 2 emissions [START_REF] Litjens | Assessment of the energy efficiency enhancement of future mobile networks[END_REF].

The increasing energy costs, the carbon footprint and the concerns for environment led to a new research direction, called Green Radio, aiming at finding solutions towards high Energy-Efficient Networks. Within this direction, several workshops and projects have been organized, for instance Celtic-OperaNet, FP7-Earth and more recently GreenTouch. Moreover, many activities have been triggered in standardization and regulatory bodies such as 3GPP, ITU, ETSI and ATIS [START_REF] Correia | Challenges and enabling technologies for energy aware mobile radio networks[END_REF] and were dedicated to greening the mobile communications.

The purpose of this study is to contribute to this emerging issue by designing, and, evaluating some advanced solutions focusing mainly on the access network as it consumes most of the energy of the network. In this work, we make use of an advanced power model developed by Interuniversity MicroElectronics Centre (IMEC) [11] [START_REF] Desset | Towards a flexible and futureproof power model for cellular base stations[END_REF] [13] [START_REF] Debaillie | A flexible and futureproof power model for cellular base stations[END_REF] within the GreenTouch project that allows to evaluate the BS power consumption. The IMEC model will serve to evaluate the performance of a reference scenario comprising BSs with standard antennas, and compare it to that of a recently proposed technology, namely Multilevel Beamforming which implements a multilevel codebook of beams [START_REF] Tall | Multilevel beamforming for high data rate communication in 5G networks[END_REF]. The multilevel codebook is constructed according to coverage requirements that enables to select the adequate focused beam for each user. The merits of the proposed beamforming approach are the following: The received signal is amplified by the high gain of the beams, specially the beams of the highest levels; The focused beams create very low interference on neighboring cells, resulting in much higher Spectral Efficiency (SE) and Energy Efficiency (EE). The proposed approach also benefits from low complexity and hence potentially low cost.

The SE is defined as the system throughput per unit of bandwidth. It has always been an efficient criterion enabling to evaluate the performance of wireless networks. The EE has become of major interest for communication systems. However, the two metrics do not progress in the same direction (as shown in Figure 1). As a result a clear trade-off should be reached in order to satisfy the network operators needs. The problem of representing and devising green policies to manage the network is of particular importance and is addressed in the paper.

The contributions of the work can be summarized as follows:

• Exploit the technology evolution in Large Scale Antenna Systems (LSAS) to build more energy-efficient networks.

• Design a LSAS solution with a multilevel codebook according to coverage requirements. • Analyze the EC of a network implementing the LSAS.

• Design green policies for selecting desired network operation trade-offs between capacity and EC. The rest of this paper is organized as follows: Section II presents the multilevel beamforming technology and the power model used in this work. Section III describes the two scenarios considered in this paper as well as a policy design framework for making the network more energy efficient. The performance results obtained from an event-based Long Term Evolution (LTE) network simulator are presented in Section IV and Section V concludes the paper.

II. SYSTEM MODEL

This Section presents the system model including the beamforming model using the multilevel codebook and the corresponding beam selection algorithm; the radio modeling and the EC calculation.

A. Multilevel beamforming

In order to reduce the complexity of beamforming implementation with LSAS, multilevel beamforming has been proposed in [START_REF] Tall | Multilevel beamforming for high data rate communication in 5G networks[END_REF]. It consists in iteratively focusing the beam towards a given location by navigating through a beam codebook designed in a multilevel manner. Figure 2 represents the coverage maps of the beams at each level of the codebook for one sector. The LSAS comprises a two dimensional array of Nx x Nz equi-spaced radiating dipoles in the x (horizontal) and z (vertical) directions of the array. The antenna model is described in [START_REF] Tall | Multilevel beamforming for high data rate communication in 5G networks[END_REF].

The best beam is selected in an iterative manner: for each user u, we navigate in the multilevel codebook starting from the lowest level (the full cell) to the narrowest beam. At each step, the beam that maximizes the Signal to Interference plus Noise Ratio (SINR) of the user u among the best beams obtained so far is selected. The algorithm stops when the best beam does not change or the narrowest beams' level is reached. This can be represented as the following optimization problem:

B * = argmax B∈Codebook SIN R u (B) (1)
where SIN R u (B) is the SINR of user u if he is served by the beam B.

B. Data rate model

We consider the downlink in this paper. The data rate for a given user i located in a sector s * at time slot t is computed according to the Shannon formula:

R i,t = α i,t W log(1 + SINR i,t ) (2)
where W represents the frequency bandwidth, α i,t -the fraction of the bandwidth allocated to user i in time slot t and SINR i,t -the SINR of a user i during a time slot t. α i,t is determined by the scheduling algorithm used and the configuration of users served by sector s * . We focus on Round-Robin scheduling thus

α i,t = 1 n s * ,t (3)
where n s,t is the number of users served by sector s * in time slot t.

The SINR is computed as follows:

SINR i,t = P s * h s * ,i (i) I i,t + N 0 (4)
where :

• P s * is the transmit power of sector s * .

• h s * ,l (k) is the channel plus antenna gain from sector s * to user l with the antenna of sector s * configured for user k. The antenna configuration mainly concerns the beam used. • N 0 is the Noise variance.

• I i,t is the interference computed as follows

I i,t = s =s * 1 n s,t j∈s P s h s,i (j) (5)
Remark: The expression for the interference corresponds to an approximation that we consider in order to simplify the computation. Indeed, the exact computation in the simulator would require to know for each BS which users are scheduled in all the other BSs (which beams are used) at each scheduling time of a millisecond. Numerical simulations have shown that the proposed approximation is tight.

C. Energy Consumption model

The EC for a single BS s is computed as follows:

E k = Ne k=1 τ k P s,k (ρ s,τ k , T s , N ant ) (6)
where N e is the number of events, τ k -the duration between two events k and k -1 and P s,k -the power consumption of the BS s during τ k . An event is defined by an arrival or a departure of a user in the network. P s,k is computed using the IMEC power model tool and it depends on ρ s,τ k which is the load of the BS s during the period τ k (we assume data traffic where a scheduled user utilizes all the BS resources. So if BS s serves at least one user during this period, then ρ s,τ k = 1, otherwise ρ s,τ k = 0). P s,k depends also on the type of the BS which we denote by T s , and the number of radiating dipoles in the LSAS, N ant . For a traditional BS with 2x2 MIMO, we have N ant = 2 and for a LSAS BS, N ant depends on the beamforming level used (see Table I). For a given BS s, there can be more than one user scheduled for transmission during the period τ k and each user may use a different beamforming level. As a result, P s,k must combine the power consumption of the different scheduled users as follows

P s,k = n s,k j=1 ω j,s,k p j,s,k (7)
where:

• n s,k is the number of users served by the BS s during τ k . • ω j,s,k is the fraction of time during which the user j is served by the BS s. As we are considering a Round Robin scheduling, ω j = 1 n s,k .

• p j,s,k is the power consumption of the user j during period k in which he was served by the BS s. It takes into account which beam is used by user j so that we can deduce exactly his power consumption value according to the IMEC tool. Figure 3 presents the values of P s,k for different configurations: traditional BS denoted as macro, LSAS BS with increasing beamforming levels denoted as Level 0 through Level 3 and with n s,k = 0 (OFF) and n s,k ≥ 1 (ON). The total EC of the whole network is the sum of the ECs of all the powered BSs

E T otal = NBS k=1 E k (8)
where N BS is the number of BSs.

III. GREEN NETWORK MODEL In order to reduce EC, we study two scenarios making use of the multilevel beamforming technique. We start with a dense hexagonal network of traditional BSs which serves as the reference scenario. FTP type elastic traffic is assumed. The coverage map of the reference scenario is presented in Figure 4.

The first scenario denoted as green scenario consists in switching off one tier of BSs (i.e. of sites) and then using multilevel beamforming to cover the areas previously occupied by the switched off BSs. The coverage map for the level 0 of beams in the multilevel codebook for this scenario is given in Figure 5. The objective of this scenario is to preserve energy by switching off some BSs but preserve the Quality of Service (QoS) using multilevel beamforming.

In order to use multilevel beamforming in the green scenario, a new beam codebook must be designed tailored specifically to the cells in Figure 5 since each cell has now a specific coverage.

The second scenario denoted as capacitive scenario consists in keeping all BSs switched on but replacing the traditional macro BSs with LSAS BSs and using the multilevel beamforming in every cell. The idea here is to significantly improve the users' SINR and thus reduce the resource utilization due to shorter sojourn time. A reduced resource utilization implies a reduced EC since the system transmits less time at full transmission power. The beam codebook presented in Figure 2 is used for the capacitive scenario. Other scenarios can be devised in order to exploit the multilevel beamforming technology and each has its advantages and drawbacks. They will perform differently in terms energy efficiency or spectral efficiency as a function of the traffic demand. We propose a framework for comparing and selecting the best scenario according to a green policy that defines a trade-off between energy efficiency and spectral efficiency. Let S be a given scenario, we define a utility function as :

µ S,a,b : χ E × χ T -→ R * + (x, y) -→ ax + by
where :

• χ E and χ T are, respectively, the sets of energy reduction factors (defined as the ratio between the EC of the reference scenario and that of the considered scenario, capacitive or green), and mean user throughput improvement factors (defined as the ratio between the throughput of the considered scenario and the reference one). • a and b are, respectively, the weights attributed to each factor that defines the green policy. The best scenario is chosen as the maximizer of the utility function. The two scenarios we consider represent two policies for operating the network, S 1 and S 2 . If µ S1,a,b (x 1 , y 1 ) ≥ µ S2,a,b (x 2 , y 2 ) then S 1 is to be deployed and vice versa. The network operator applies its policy by adjusting the weights a and b.

IV. NUMERICAL RESULTS

A. Simulation setup

We use a LTE event-based simulator coded in Matlab. An event corresponds to an arrival or a departure in the whole network. We assume that the users arrive according to a Poisson process, download a file of exponential size with mean of 4 Mbits, and leave the network as soon as their download is complete. No mobility is considered in the simulations but the users arrive at random locations in the whole network. We focus on a dense urban case with an inter-site distance of 500 m.

In the reference scenario, we consider a hexagonal network of 19 tri-sectorized conventional macro BSs (see Figure 4) with 2x2 MIMO antennas, antenna heights of 30 m and maximum antenna gains of 18 dBi. In the capacitive and green scenarios, we replace the macro BSs of the reference scenario with LSAS BSs in order to perform multilevel beamforming.

The simulation parameters are summarized in Table II.

B. Performance results

We consider four performance indicators to compare the different scenarios: Mean User Throughput (MUT), Cell-Edge Throughput (CET), load and EC. Table III presents the performance results for each scenario with varying arrival rates.

The performance results show that multilevel beamforming is a promising solution for increasing the throughput and decreasing the EC of the network. The user throughput improvement comes from the increase of the antenna gain obtained by using a large antenna array. The sojourn time of the users in the network is reduced resulting in lower load values in the capacitive scenario as shown in Table III. It is noted that the load value is higher in the green scenario compared to the reference scenario. This can be explained by the fact that we shut down some cells so each of the remaining ones will serve more users.

The EC reduction of the capacitive scenario with respect to the reference one is mainly due to the shorter mean sojourn time of the former. In the green scenario, the reduction of the EC is the result of the decrease in the number of the active BSs.

The EC reduction in both capacitive and green scenarios is also linked to the transmission schemes in terms of the beam level selection probabilities.

Multilevel beamforming is clearly an efficient solution, both for improving energy and throughput. However, the capacitive scenario can become greener (with lower energy consumption) than the green scenario in high load scenarios. The best scenario must then be selected and deployed according to the traffic demand. The energy and spectral efficiency of the capacitive and the green scenarios depend on the traffic demand. For very high loads (λ ≥ 15 in Table III), the capacitive scenario is both more spectrally efficient and more energy efficient than the green scenario. This is due to the fact that in multilevel beamforming, the most used beams to serve users are in the highest level (see Figure 6) which consumes 8 times more energy than the level 0 beams (see Figure 3) which are used during idle times. Thus, the EC will highly depend on the load of the system. In high traffic demand scenarios, the load of the capacitive scenario is much lower than in the green scenario, thus the capacitive scenario uses less often the highest levels of beams leading to a lower EC.

C. Energy-Efficient Network management

Figure 7 presents the performance of the green and capacitive scenarios in the EC-MUT plane. Zone one (on the left of the dashed red line) corresponds to the zone where it can be beneficial to shift from one scenario to the other. This shift corresponds to the trade-off between the EC reduction and the user throughput improvement. The choice of the scenario depends on the operator policy, namely on the relative importance given to each criterion.

Zone 2 corresponds to the values of λ where the capacitive scenario is better both for EC and MUT. The red dashed line sets the limit between Zones 1 and 2. It corresponds to the arrival rate λ 4 at which at least one BS in the Green Scenario has a load equal to 80%, and starting from which the green scenario is no longer interesting.

We apply the green network policy framework discussed in Section III. Three different policies are considered with corresponding parameters specified in Table IV. The best policy performance is highlighted in bold. The results for applying these policies are also shown in Figure 7. For a given policy, the utility can be monitored, and the shift between the two scenarios occurs when one utility gets better than the other. Alternatively, one can compute offline a table which maps arrival rates onto utility values and identify the switching points from one scenario to the other.

V. CONCLUSION

We have presented in this paper two scenarios in which multilevel beamforming has been used to decrease EC either by switching off some BSs or by increasing the system capacity thus reducing its load. We also presented a framework for designing a green policy that allows to switch between different scenarios according to the network operator objectives. This framework implements a trade-off between energy efficiency and spectral efficiency. More scenarios can be devised for leveraging the potential of multilevel beamforming in making mobile networks more energy efficient such as reducing the transmission power when using beamforming or load balancing using adequately designed codebooks. It is noted that the performance results obtained for the energy consumption are highly dependent on the power model used, i.e. here the IMEC power model. Other power models can be considered and the performance of the scenarios proposed can be evaluated with these models. This paper has the ambition to foster research on green wireless networks management using self-organized features. The paper has shown that a trade-off has to be found between SE and EE. Moreover, different scenarios have to be studied e.g. rural, in order to evaluate their operability.

Fig. 1 :

 1 Fig. 1: Trade-off between Energy Efficiency and Spectral Efficiency.

Fig. 2 :

 2 Fig. 2: Multilevel beam coverage maps.

Fig. 3 :

 3 Fig. 3: Power consumption values (Watt). The terms ON and OFF stand respectively for the cases with and without traffic in the cell.

Fig. 4 :

 4 Fig. 4: Hexagonal network deployment.

Fig. 5 :

 5 Fig. 5: Coverage map after switching off the middle ring.

Fig. 6 :

 6 Fig. 6: Probabilities of selection of the different beam levels in the Capacitive and the Green scenarios (for λ 6 as an example).

Fig. 7 :

 7 Fig. 7: Illustration of the different deployment policies.

TABLE I :

 I Number of antennas in the different beamforming levels

	Beamforming level	Nx x Nz = Nant
	Level 0	28
	Level 1	5 x 14 = 70
	Level 2	10 x 14 = 140
	Level 3	20 x 14 = 280

TABLE II :

 II Simulator characteristics

	Network parameters
	Macro Cell layout	hexagonal trisector
	Antenna height	30 m
	Bandwidth	10MHz
	Scheduling type	Round Robin
	Channel characteristics
	Thermal noise	-174 dBm/Hz
	Path loss (d in km) 128.1 + 37.6 log 10 (d) dB
	Shadowing	Log-normal (6dB)
	Traffic characteristics
	Service type	FTP
	Average file size	4 Mbits

TABLE III :

 III Simulation results

	Arrival rate	Scenario	EC	MUT	CET	Load
	[users/s/Km 2 ]		[Wh]	[Mbps]	[Mbps]	[%]
		Reference	135.67	26.76	6.7	4
	λ1 = 3	Capacitive	42.92	40.74	11.89	3
		Green	32.06	30.67	7.84	5.5
		Reference	139.73	25	5.98	9.5
	λ2 = 6	Capacitive	48.43	39.7	11.6	6
		Green	39.48	29.6	5.8	11.7
		Reference	146.17	21.9	4.6	18
	λ3 = 10	Capacitive	53.94	38.86	9.98	9
		Green	48.49	26.23	3.8	20
		Reference	151.66	19.7	4.08	25
	λ4 = 13	Capacitive	59.5	36.7	8.1	13
		Green	57.11	24.14	2.7	30
		Reference	154.4	18.68	3.79	28.7
	λ5 = 15	Capacitive	62.17	35.9	7.2	15
		Green	65.78	21.16	1.3	35
		Reference	166.65	14.18	2.15	45
	λ6 = 20	Capacitive	70.14	33.9	5.7	22
		Green	81.64	14.8	0.24	49

TABLE IV :

 IV Utility values for the two scenarios

	Policies	Arrival rate	µ capacitive,a,b	µ green,a,b
	Policy 1	λ1	4.68	5.38
	(a,b) = (1,1)	λ2	4.47	4.72
		λ3	4.48	4.2
	Policy 2	λ1	2.67	3.31
	(a,b) = (0.3 , 0.7)	λ2	2.5	2.83
		λ3	2.43	2.47
	Policy 3	λ1	2.01	2.07
	(a,b) = (0.7 , 0.3)	λ2	1.98	1.89
		λ3	2.06	1.74
	Policy 3
	Policy 1	Policy 2