Romuald Rocher

Daniel Menard

Pascal Scalart

Olivier Sentieys Members

Analytical Approach for Numerical Accuracy Estimation of Fixed-Point Systems Based on Smooth Operations

Keywords: Fixed-point arithmetic, quantization noises, adaptive filters, accuracy evaluation

In embedded systems using fixed-point arithmetic, converting applications into fixed-point representations requires a fast and efficient accuracy evaluation. This paper presents a new analytical approach to determine an estimation of the numerical accuracy of a fixed-point system, which is accurate and valid for all systems formulated with smooth operations (e.g. additions, subtractions, multiplications and divisions). The mathematical expression of the system output noise power is determined using matrices to obtain more compact expressions. The proposed approach is based on the determination of the timevarying impulse-response of the system. To speedup computation of the expressions, the impulse response is modelled using a linear prediction approach. The approach is illustrated in the general case of time-varying recursive systems by the Least Mean Square (LMS) algorithm example. Experiments on various and representative applications show the fixedpoint accuracy estimation quality of the proposed approach. Moreover, the approach using the linear-prediction approximation is very fast even for recursive systems. A significant speed-up compared to the best known accuracy evaluation approaches is measured even for the most complex benchmarks.

I. INTRODUCTION

Digital image and signal processing applications are mostly implemented in embedded systems with fixed-point arithmetic to satisfy energy consumption and cost constraints [START_REF] Evans | Modem Design, Implementation, and Testing Using NI's LabVIEW[END_REF], [START_REF] Eyre | The evolution of DSP processors[END_REF], [START_REF] Strauss | DSP chips take on many forms[END_REF], [START_REF] Eyre | DSPs court the consumer[END_REF]. Fixed-point architectures manipulate data with relatively small wordlengths, which could offer the advantages of a significantly lower area, latency, and power consumption. Moreover, floating-point operators are more complex in their processing of exponent and mantissa and, therefore, chip area and latency are more important than for fixed-point operators.

Nevertheless, applications in their earliest design cycle are usually designed, described and simulated with floating-point data types. Consequently, application must be converted into a fixed-point specification where implementation is considered. To reduce application time-to-market, tools to automate fixed-point conversion are needed. Indeed, some experiments [START_REF] Clark | Accelerating Fixed-Point Design for MB-OFDM UWB Systems[END_REF] have shown that the manual fixedpoint conversion process can represent from 25% to 50% of the total development time. Likewise, in a survey conducted by a design tool provider [START_REF] Hill | Acceldsp synthesis tool floating-point to fixed-point conversion of matlab algorithms targeting fpgas[END_REF], the fixed-point conversion was identified by a majority of respondents as one of the most difficult aspects during FPGA implementation of a DSP application. The aim of fixed-point conversion tools is to determine the optimized fixed-point specification that minimizes the implementation cost and provides a sufficient computation accuracy to maintain the application performance.

The fixed-point conversion process is divided into two steps to determine the number of bits for both integer and fractional parts. The first step is the determination of the binary-point position from the data dynamic range [START_REF] Wu | Dynamic-Range Estimation[END_REF]. The integer part word-length of all data inside the application is determined in order to avoid overflow. The second step is the determination of the fractional part. The limitation of the number of bits for the fractional part leads to an unavoidable error resulting from the difference between the expected value and the real finite precision value. Obviously, numerical accuracy effects must be limited as much as possible to ensure algorithm integrity and application performance. Therefore, this stage of the fixed-point conversion flow corresponds to an optimization process in which the implementation cost is minimized under an accuracy constraint. The fractional part word-length is optimized with an iterative process requiring to evaluated the numerical accuracy numerous times. Therefore precise and fast numerical accuracy estimator is a key element of any efficient fixed-point conversion method.

The numerical accuracy of an application can be obtained using two approaches: fixed-point simulations and analytical approaches. Accuracy evaluation based on fixed-point simulations [START_REF] Belanovic | Automated Floating-point to Fixed-point Conversion with the fixify Environment[END_REF], [START_REF] Kim | Fixed-Point Error Analysis and Word Length Optimization of 8x8 IDCT Architectures[END_REF] is very time consuming because a new simulation is needed as soon as a fixed-point format changes in the system resulting in prohibitive optimization time for the fixed-point conversion process. Analytical approaches attempt to determine a mathematical expression of an estimation of a numerical accuracy metric, which leads to short evaluation times compared to fixed-point simulation based approaches. Consequently, the rest of the paper focuses on improving the quality and the application range of analytical approaches for numerical accuracy evaluation.

Numerical accuracy can be evaluated by determining one of the following metrics: quantization error bounds [START_REF] De Figueiredo | Affine Arithmetic : Concepts and Applications[END_REF], [START_REF] Alefeld | Introduction to Interval Computations[END_REF], number of significant bits [START_REF] Chesneaux | The fixed cadna library[END_REF] or quantization noise power [START_REF] Menard | Automatic Evaluation of the Accuracy of Fixed-point Algorithms[END_REF], [START_REF] Shi | A perturbation theory on statistical quantization effects in fixed-point DSP with non-stationary inputs[END_REF], [START_REF] Caraiscos | A Roundoff Error Analysis of the LMS Adaptive Algorithm[END_REF]. The error bound metric (or the number of significant bits) is used for critical systems where the quantization effects could have catastrophic consequences [START_REF] Lions | ARIANE 5, Flight 501 Failure[END_REF], [START_REF] Blair | Patriot Missile Defense: Software Problem Led to System Failure at Dhahran, Saudi Arabia[END_REF]. The error has to be kept small enough to ensure the system behavior within acceptable bounds. The quantization error power is preferred when a slight degradation of application performance due to quantization effects is acceptable. This design approach concerns a great number of digital signal processing applications. Moreover, the quantization error can be linked to the degradation of the application performance using for example the technique presented in [START_REF] Menard | Accuracy Constraint Determination in Fixed-Point System Design[END_REF]. In this paper, the quantization noise power is used as the numerical accuracy metric. This paper proposes an efficient and accurate analytical approach to estimate the numerical accuracy of all types of fixed-point systems based on smooth operations (an operation is considered to be smooth if the output is a continuous and differentiable function of its inputs, as in the case of arithmetic operations). From the signal flow graph of the application, the associated tool automatically generates the source code implementing the analytical expression of the output quantization noise power. In this expression, the gains between the different noise sources and the output are computed from the impulse response (IR), which can be time-varying. This expression is unbiased but introduces infinite sums. To reduce the computation time of the gains, an approach based on linear prediction is also proposed. The approach is based on a matrix model which leads to a compact expression in the case of algorithms expressed with vectors and matrices, such as the Fast Fourier Transform (FFT) and the Discrete Cosine Transform (DCT). Compared to approaches based on fixed-point simulations, the time required to evaluate the fixed-point accuracy is also dramatically reduced. Existing analytical approaches are only valid for linear and time invariant (LTI) systems [START_REF] Lopez | Fast and accurate computation of the roundoff noise of linear time-invariant systems[END_REF], [START_REF] Lopez | Improved Interval-Based Characterisation of Fixed-Point LTI Systems With Feedbacks Loops[END_REF], or for non-LTI and non recursive systems [START_REF] Menard | SQNR determination in non-linear and non-recursive fixed-point systems[END_REF]. In [START_REF] Caffarena | SQNR Estimation of Fixed-Point DSP Algorithms[END_REF], non-LTI systems with feedbacks are handled. But, as mentioned in [START_REF] Caffarena | SQNR Estimation of Fixed-Point DSP Algorithms[END_REF], the execution time can be quite long for systems containing feedbacks. Indeed, to obtain an accurate estimation, this execution time depends on the length of the system impulse response. Compared to the last-mentioned approach, a significant speed-up is obtained with the proposed approach thanks to the proposed linear prediction approach.

This paper is organized as follows. In Section II, existing methods to evaluate the numerical accuracy of fixed-point applications are detailed. The comparison of simulation-based and analytical approaches highlights the benefits of analytical approaches. In Section III, the proposed approach for numerical accuracy evaluation is defined. The quantization noise propagation model is presented and the output noise power expression is expressed. A technique based on linear prediction is also presented to reduce the time to obtain the power expression. In Section IV, the approach is tested on several LTI and non LTI systems. The quality of the estimation and the time to obtain the power expression are evaluated on various applications and benchmarks. Finally, Section V concludes the paper.

II. EXISTING APPROACHES FOR NUMERICAL ACCURACY EVALUATION

In the fixed-point conversion process, the numerical accuracy evaluation is the most critical part, mainly because this evaluation is successively performed in the iterative word-length optimization process. The numerical accuracy can be evaluated using analytical or simulation-based approaches.

A. Accuracy Evaluation Based on Fixed-Point Simulations

In fixed-point simulation-based approaches, the application performance and the numerical accuracy are directly determined from fixedpoint (bit-true) simulations. To obtain accurate results, the simulations must use a high number of samples Nsamp. Let Nops denote the number of operations in the algorithm description, NI -the number of times the accuracy is evaluated and τ -the simulation time for one sample. Therefore, the total system simulation time is

Tsim = τ NsampNopsNI . (1)
To reduce word-length optimization time, approaches based on fixed-point simulations aim at reducing the values of some parameters of (1). First, fixed-point simulators have to minimize the computing time τ by using efficient data types to simulate fixed-point arithmetic. Most of them use object oriented language and operator overloading [START_REF] Kim | Fixed-Point Error Analysis and Word Length Optimization of 8x8 IDCT Architectures[END_REF] to implement fixed-point arithmetic. However, the fixed-point execution time of an algorithm is significantly longer than the one of a floating-point simulation [START_REF] Coors | Integer Code Generation For the TI TMS320C62x[END_REF]. To reduce this overhead, new data types have been introduced in [START_REF] Kim | Fixed-Point Error Analysis and Word Length Optimization of 8x8 IDCT Architectures[END_REF]. They use more efficiently the floating-point units of general-purpose processors. With this new data type, the fixed-point simulation time is still 7.5 times longer than the one obtained with floating-point data types for a simple application like a fourth-order Infinite Impulse Response (IIR) filter [START_REF] Kim | Fixed-Point Error Analysis and Word Length Optimization of 8x8 IDCT Architectures[END_REF]. However, the optimization time to accelerate the fixed-point simulations is not given in [START_REF] Kim | Fixed-Point Error Analysis and Word Length Optimization of 8x8 IDCT Architectures[END_REF], which could actually annihilate the obtained gain for the fixed-point simulation. The approach proposed in [START_REF] Roy | An algorithm for trading off quantization error with hardware resources for Matlab-based FPGA design[END_REF] was developed for MATLAB code but its simulation time is still significantly higher.

Finally, some methods [START_REF] Belanovic | Automated Floating-point to Fixed-point Conversion with the fixify Environment[END_REF] attempt to reduce Ni, the number of iteration of the optimization process by limiting the search space, and thus leading to sub-optimal solutions.

B. Analytical Accuracy Evaluation

The general idea of analytical approaches is to determine a mathematical expression that provides estimations of the output quantization noise power. The computation of this analytical expression could be time-consuming, but it is executed only once. Then, the numerical accuracy is determined very quickly by evaluating the mathematical expression for a given word-length of each data. Existing analytical approaches for numerical accuracy evaluation are based on perturbation theory [START_REF] Constantinides | Perturbation Analysis for Word-length Optimization[END_REF], [START_REF] Wadekar | Accuracy Sensitive Word-Length Selection for Algorithm Optimization[END_REF]. The quantization of a signal x generates a noise bx considered as a signal perturbation. Let xi and y be respectively the n inputs and the output of an operation f , and bx i and by their perturbations. The perturbation of the output value of a differentiable operation is given by by

= bx 1 df dx1 + .... + bx n df dxn . (2)
Then, the expression of the noise power P b corresponding to the second-order moment of the output quantization noise can be expressed as

P b = Ne i=1 σ 2 b i Ki + Ne i,j m b i m b j Lij, (3)
where Ne is the number of quantization noises, m b i and σ 2 b i are respectively the mean value and the variance of the noise bi and the terms Ki and Lij depend on system characteristics. The expressions of Ki and Lij differ according to the considered approach and will be presented in the following subsection. Determining the terms Ki and Lij is crucial for the analytical accuracy evaluation process.

1) Simulation-Based Computation: Approaches based on hybrid techniques [START_REF] Shi | A perturbation theory on statistical quantization effects in fixed-point DSP with non-stationary inputs[END_REF], [START_REF] Constantinides | Word-length optimization for differentiable nonlinear systems[END_REF], [START_REF] Fiore | Efficient approximate wordlength optimization[END_REF] have been proposed to compute the coefficients Ki and Lij of (3) from a set of simulations. In [START_REF] Shi | A perturbation theory on statistical quantization effects in fixed-point DSP with non-stationary inputs[END_REF], these Ne × (Ne + 1) coefficients are obtained by solving a linear system in which Ki and Lij are the variables. The other elements of (3) are determined by carrying out fixed-point simulations. The wordlengths of the data are set to control each quantizer and to analyze its influence on the output. For each simulation, the output noise power is evaluated. The statistical parameters (mean and variance) of each noise source bi are extracted from the data word-length. At least Ne × (Ne + 1) fixed-point simulations are required to solve the system of linear equations. A similar approach to [START_REF] Shi | A perturbation theory on statistical quantization effects in fixed-point DSP with non-stationary inputs[END_REF] is used in [START_REF] Fiore | Efficient approximate wordlength optimization[END_REF] to obtain the coefficients by simulation.

2) Affine Arithmetic Based Simulations: In approach based on affine arithmetic [START_REF] Caffarena | SQNR Estimation of Fixed-Point DSP Algorithms[END_REF], the coefficients Ki and Lij are computed with an affine arithmetic based simulation. This approach was proposed for LTI in [START_REF] Lopez | Fast and accurate computation of the roundoff noise of linear time-invariant systems[END_REF], [START_REF] Lopez | Improved Interval-Based Characterisation of Fixed-Point LTI Systems With Feedbacks Loops[END_REF] and recently, for non-LTI systems in [START_REF] Caffarena | SQNR Estimation of Fixed-Point DSP Algorithms[END_REF].

In affine arithmetic, a data x is defined by a linear relation

x = x0 + ǫ1x1 + ... + ǫnxn, (4)
where ǫi is an uncertainty term included in the interval [-1, 1] and such that ǫi is independent of ǫj ∀i = j. Each noise source is defined by an affine form whose behavior is controlled by the mean and the variance of the noise source. Then, noise source expressions are propagated throughout the system thanks to an affine arithmetic based simulation. The values of Ki and Lii are extracted from the affine form of the output noise. The approach proposed in [START_REF] Caffarena | SQNR Estimation of Fixed-Point DSP Algorithms[END_REF] is general and can handle systems based on smooth operations. This estimation is accurate and the execution time is low for non-recursive systems. Nevertheless, in the case of recursive systems, several iterations for the affine arithmetic based simulation are required to converge to stable values. Thus, the execution time depends on the length of the impulse response and on the number of noise sources. As mentioned in [START_REF] Caffarena | SQNR Estimation of Fixed-Point DSP Algorithms[END_REF], the execution time can be quite long for systems containing feedback loops.

3) Approaches Based on System Characteristics: The approaches presented in [START_REF] Constantinides | Wordlength optimization for linear digital signal processing[END_REF], [START_REF] Menard | Analytical fixed-point accuracy evaluation in linear time-invariant systems[END_REF] provide an approximation of the mathematical expression of the output quantization noise power in the case of Linear and Time Invariant (LTI) systems. The principle of [START_REF] Menard | Analytical fixed-point accuracy evaluation in linear time-invariant systems[END_REF] is based on the system transfer functions and the expressions of the coefficients Ki and Li are

Ki = 1 2π π -π |Hi(e jΩ )| 2 dΩ, (5)
Lij = Hi(1)H * j (1), (6)
where Hi(z) is the transfer function between the white noise source bi and the system output. In [START_REF] Menard | Analytical fixed-point accuracy evaluation in linear time-invariant systems[END_REF], the authors proposed a technique to automatically obtain the transfer function from the application description.

For non-LTI and non-recursive systems, an extended approach was presented in [START_REF] Menard | SQNR determination in non-linear and non-recursive fixed-point systems[END_REF]. Each noise source bi(n) is propagated through Mi operations υM i , which leads to the calculation of the output noise b ′ i (n). The noise b ′ i (n) is the product of each noise source bi(n) and of different signals α k associated to operations υM i included in noise source bi(n) propagation. So, Ki and Lij terms are equal to

Ki = E[ M i k=1 α 2 k ] (7)
Lij = E   M i k=1 M j g=1 α k αg   . (8)
It is noteworthy that this approach is adequate for both non-LTI and non-recursive systems, but not to non-LTI and recursive systems such as the clan of adaptive filters.

C. Summary

While approaches based on fixed-point simulations lead to very long computing time and thus require to limit the search space during word-length optimization, analytical approaches are significantly less time-consuming. Unfortunately, they are limited to a reduced class of systems and are not always efficient in terms of execution time. To cope with the lack of efficient and general approaches, an analytical approach for all types of systems based on smooth operations is presented in the next section.

III. PROPOSED ACCURACY EVALUATION APPROACH

In this section, a new approach for evaluating the quantization noise at the output of any system based on smooth operations (addition, subtraction, multiplication and division) is detailed. First, the models associated with the quantization process and the quantization noise propagation model are introduced. Then, the modelling of the system is presented and the output noise power expression is detailed.

A. Quantization Noise Models 1) Quantization Noise Source: The quantization process can be modelled by the sum of the original signal and of a uniformly distributed white noise [START_REF] Widrow | Statistical Theory of Quantization[END_REF], [START_REF] Sripad | A Necessary and Sufficient Condition for Quantization Error to be Uniform and White[END_REF]. This quantization noise is uncorrelated with the signal and the other noise sources. Such a model is valid provided that the signal has a sufficiently large dynamic range compared to the quantization step. According to the quantization mode, the noise Probability Density Function (PDF) will differ. Three quantization modes are usually considered: truncation, conventional rounding, and convergent rounding. In the truncation mode, the Least Significant Bits (LSB) are directly eliminated. The resulting number is always smaller than or equal to the number obtained before quantization, and, therefore, the quantization noise is always positive. Consequently, the mean of the quantization noise is not equal to zero. To reduce the bias due to truncation, the rounding quantization mode is often used. In conventional rounding, the data are rounded to the nearest value representable in the reducedaccuracy format. For numbers located at the midpoint between two consecutive representable values, the data are rounded-up always to the larger output value. This technique leads to a (small) bias for the quantization noise. To eliminate this quantization noise bias, the convergent rounding can be used as well. In this case, the numbers located at the midpoint between two consecutive representable values are, with equal probability, rounded to the higher or lower output value.

Let wF P denote the number of bits for the fractional part after the quantization process and k-the number of bits eliminated during the quantization. The quantization step q after the quantization is equal to q = 2 -w F P . The quantization noise mean and variance are given in Table I for the three considered quantization modes [START_REF] Constantinides | Truncation Noise in Fixed-Point SFGs[END_REF], [START_REF] Menard | Quantization Mode Opportunities in Fixed-Point System Design[END_REF]. 2) Quantization Noise Propagation: The aim of this part is to define the model for the propagation of the quantization noise for smooth operations. We specifically concentrate on an operation f with two inputs x and y and on the calculation of its output z = f (x, y). The finite precision values x, y and z are respectively the sum of infinite precision values x, y and z with their associated quantization noises bx, by and bz. The propagation model is then defined such as bz is a linear combination of bx and by bz = α1bx + α2by.

Quantization Truncation Conventional Convergent mode rounding rounding Mean q 2 (1 -2 -k ) q 2 (2 -k ) 0 Variance q 2 12 (1 -2 -2k ) q 2 12 (1 -2 -2k ) q 2 12 (1 + 2 -2k+1 )
(

) 9
α1 and α2 are obtained from a first order Taylor development at the output of the differentiable operator f defined as

z = f ( x, y) = f (x, y) + ( x -x) ∂f ∂x (x, y) + ( y -y) ∂f ∂y (x, y). (10)
Thus, the expressions of α1 and α2 are

α1 = ∂f ∂x (x, y) and α2 = ∂f ∂y (x, y). ( 11
) Operation α1 α2 z = x ± y 1 ± 1 z = x × y y x z = x y 1 y -x y 2 TABLE II VALUES α 1 AND α 2 OF EQUATION (9) FOR ARITHMETIC OPERATIONS Operation Ax Dx Ay Dy Z = X ± Y 1 1 ±1 1 Z = X × Y 1 Y X 1 TABLE III VALUES A AND D OF EQUATION (12) FOR OPERATIONS {+, -, ×}
The values of α1 and α2 are summarized in Table II for the arithmetic operations.

For non-scalar noise sources (vectors or matrix), the model of ( 9) is not valid anymore. For example, in the case of the multiplication of two matrices X and Y associated with noise matrices Bx and By , the result Z is equal to Z = XY. Then, the output noise matrix Bz is given by Bz = Bx.Y + X.By.

In the rest of the paper, upper-case letters in bold will represent matrices, lower-case letters in bold will represent vectors, and italic will represent scalars.

Since commutativity is not valid for non-scalar multiplications and since we seek to obtain a general model, each noise source on the operation input is multiplied by a signal term on the left or on the right. Therefore, in the general case, the noise source propagation is expressed as the multiplication of each noise input by two matrices (A on the left and D on the right) as

Bz = AxBxDx + AyByDy. ( 12
)
Table III defines values of A and D according to the type of operation. For the mathematical expressions of the general approach, A and D will be considered as matrices. However, when applied to different examples, depending on the application context, these terms will be noted A and D for matrices, a and d for vectors and a and d for scalars.

To illustrate the proposed approach, the example of a FIR (Finite Impulse Response) filter is firstly considered. In this case, x(n) is an N -size input vector, bx(n) is its associated noise, and h is the filter coefficient vector. So, the output noise by(n) due to input noise propagation is equal to

by(n) = h t bx(n). (13)
The terms ax and dx are therefore given by

ax = h t , (14)
dx = 1. ( 15
)
To illustrate the need of the two matrices A and D, the adaptive filter based on the Affine Projection Algorithm (APA) is then considered. APA has been proposed in [START_REF] Ozeki | An adaptative filtering algorithm using an orthogonal projection to an affine subspace and its properties[END_REF] to have a faster convergence compared to the Least Mean Square (LMS) algorithm and to reduce complexity compared to Recursive Least Square (RLS) algorithm. Let X(n) denote an N × K matrix made-up of the K last observation vectors of the input, w(n)-the coefficient vector and e(n) the error vector. The updated coefficient equation is given by

w(n + 1) = w(n) + µX(n)(X t (n)X(n)) -1 e(n). (16)
The noise matrix Bx(n) associated to the inversion propagation matrix (X t (n)X(n)) -1 is then multiplied by two terms

A = µX(n), (17)
d = e(n). (18)

B. System Model

Each noise source propagates throughout the system and contributes to the global output noise. In this section, the system characterizing each noise propagation is defined to compute afterward the expression of the system output noise power. Let Ne denote the number of noise sources inside the system. As shown in Table II, the noise expression does not contain crossed noise terms. So, each noise source bi(n) at time n contributes to the output noise by a noise b ′ i (n) that is independent of all other noise terms. The system output noise by(n) is the sum of all contributions as presented in Figure 1 and expressed by Each noise contribution b ′ i (n) is obtained from the propagation of noise sources bi(n) through the system Si. Thus, to determine the complete expression of the output noise by(n), each subsystem Si must be characterized analytically.

by(n) = Ne i=1 b ′ i (n). (19)
+ b 1 (n) b i (n) b N e (n) b 1 (n) b i (n) b N e (n) b y (n) S 1 S i S N e
′ i (n) is b ′ i (n) = Q i k=0 gi(k)bi(n -k) + P i m=1 fi(m)b ′ i (n -m), (20)
where gi(k) represents the contribution of noise source bi at time (nk) to the output noise and fi(m)-the contribution of noise b ′ i at time (nm). The terms fi and gi may be time-varying and depend on the system implementation. For LTI systems, the terms fi and gi are constant and correspond to filter coefficients. Nevertheless, to compute the output noise power expression from the statistical parameters of the noise bi(n), (20) must be developed to express contribution b ′ i (n) with only input noise terms bi. This development introduces the time-varying impulse response associated with the system Si.

2) Time-Varying Impulse Response: In this part, the time-varying impulse response hi of the system Si is determined. By developing the recurrence in [START_REF] Kim | Fixed-Point Error Analysis and Word Length Optimization of 8x8 IDCT Architectures[END_REF], the noise term becomes

b ′ i (n) = n k=0 h (n) i (k)bi(n -k) (21)
where h

(n) i (k)
is the time-varying impulse response of the system Si which formulates the noise b ′ i (n) only from the samples of the noise source bi(k). The term h

(n) i (k) represents the contribution of the noise source bi(n -k) at time n -k to generate b ′ i (n) at time n. For the sake of clarity, h (n) i (k) is noted by hi(k) since output time
is always considered at time n. This impulse response is obtained from the terms fi and gi with hi(k)

= P i j=1 fi(j)hi(k -j) + gi(k) (22)
(22) was obtained in the case of scalar noise sources. The model can be extended to vector (or matrix) noise sources by using vector (or matrix) terms for the time-varying impulse response. In this case, the time-varying impulse response is equivalent to the multiplication of the noise source B i (k) (to be more general, noises are considered as matrices in the following) by two terms A i (k) and D i (k). In case of non-scalar noise source, (21) can be rewritten as

B ′ i (n) = n k=0 A i (k)B i (n -k)D i (k). (23)
The output noise By(n) is finally the sum of all noise source contributions and can be expressed as

By(n) = Ne i=1 n k=0 A i (k)B i (n -k)D i (k). (24)
The output noise By(n) and input noises B i (n) have formal dimensions of MyxNy and MixNi. So, A i (k) and D i (k) are MyxMi and NixNy matrices. Their numerical values depend on the system characteristics and can be computed from a analysis of the considered application.

C. Output Noise Power Expression 1) Noise Power:

The output noise power P b is defined as the second order moment of [START_REF] Menard | SQNR determination in non-linear and non-recursive fixed-point systems[END_REF] and can be expressed as

P b = E[T r(By(n)By t (n))], (25)
which is equivalent to

P b = N e i=1 N e j=1 n k=0 n m=0 E T r A i (k)B i (n -k)D i (k) D j t (m)B j t (n -m)A j t (m) . (26)
Equation ( 26) can be split into three parts by separating the noise B i from the other noises B j when j = i and the terms

B i (n -k) from B i (n -m) when m = k. Theorem 1: Each term of D i (k)D i t (k)
is less or equal than its trace value.

Proof: The matrix D i (k)D i t (k)
is real, symmetric and positive and therefore diagonalizable. Thus, eigenvalues of D i (k)D i t (k) are positive, and each of them is less than the trace value. Indeed, this matrix D i (k)D i t (k) can be approximated by its trace T r(D i (k)D i t (k)). In the case where the noises are vectors, then D i (k)D i t (k) is a scalar and is equal to its trace. Moreover, by applying Theorem 1, (26) can be reformulated and leads to [START_REF] Menard | Automatic Evaluation of the Accuracy of Fixed-point Algorithms[END_REF] given at page 6.

Following the method described in [START_REF] Caraiscos | A Roundoff Error Analysis of the LMS Adaptive Algorithm[END_REF], the expected value of the noise and signal terms are computed separately. For the two matrices A and D made-up of signal terms and for a noise B non-correlated with these two signal terms, the following relation is used:

E[ABD] = E[AE[B]D].
Moreover, the different quantization noises B i are assumed to be white and independent so the intercorrelation between two noises B i and B j is equal to

φ B i B j (k -m) = E[B i (n -k)B j (n -m)] (28) = σ 2 b i δ(i -j)δ(k -m)I N + m b i m b j 1 N
where m b i and σ 2 b i represent mean and variance of each element of input noises B i (n), I N is the identity matrix and 1 N is the Nsize matrix composed of 1. Thus, the output noise power P b can be summarized as [START_REF] Ozeki | An adaptative filtering algorithm using an orthogonal projection to an affine subspace and its properties[END_REF] given at page 6.

By introducing the terms Ki and Lij, (29) is reformulated as

P b = N e i=1 σ 2 b i Ki + N e i=1 N e j=1 m b i m bj Lij. ( 30
)
Ki represents the gain on the variance for the noise source bi and is equal to

Ki = n k=0 E T r D i (k)D i t (k) T r A i (k)A i t (k) . (31)
Lij represents the gain on the mean for the couple of noise sources bi and bj and is equal to

Lij = n k=0 n m=0 E T r A i (k)1 N D i (k)D j t (m)1 N t A j t (m) .
(32) The expressions of Ki and Lij are obtained by a single floatingpoint simulation since they only involve signal terms. All the terms included in these expressions are stored from this single floatingpoint simulation and the terms Ki and Lij are then computed. These terms are independent from noise sources and lead to constants in the output noise power expression. Noise statistics m b i and σ 2 b i depend on fixed-point formats and are the variables describing the output noise power expression. They are computed from the data word-length of each data and operation in the application graph.

It should be noted that ( 29) is unbiased since no restrictive assumption was made about the system. Ki and Lij are defined by infinite sums. However, different cases must be taken into account. If the system is LTI, the matrices A and D are constant, [START_REF] Shi | A perturbation theory on statistical quantization effects in fixed-point DSP with non-stationary inputs[END_REF] and [START_REF] Sripad | A Necessary and Sufficient Condition for Quantization Error to be Uniform and White[END_REF] can be simplified, and it leads to the results obtained in equation ( 5). If the system is not recursive, the length of the time-varying impulse responses is finite and the sums in [START_REF] Ozeki | An adaptative filtering algorithm using an orthogonal projection to an affine subspace and its properties[END_REF] are finite and can be computed directly.

For the recursive case, these infinite sums are approximated by computing the first p elements of the sum, and therefore [START_REF] Shi | A perturbation theory on statistical quantization effects in fixed-point DSP with non-stationary inputs[END_REF] and [START_REF] Sripad | A Necessary and Sufficient Condition for Quantization Error to be Uniform and White[END_REF] are computed from k = 0 to p with p depending on the signal correlation inside the terms K and L. Nevertheless, according to the different experimentations that have been carried out, a number p equal to or less than 500 leads to very realistic and accurate results. Moreover, to compute [START_REF] Ozeki | An adaptative filtering algorithm using an orthogonal projection to an affine subspace and its properties[END_REF], the ergodic assumption is used so that

P b = N e i=1 n k=0 E T r(D i (k)D i t (k))T r A i (k)B i (n -k)B i t (n -k)A i t (k) + N e i=1 n k = 0 m = k E T r A i (k)B i (n -k)D i (k)D i t (m)B i t (n -m)A i t (m) + N e i = 1 j = i n k = 0 m = 0 E T r A i (k)B i (n -k)D i (k)D j t (m)B j t (n -m)A j t (m) (27)
P b = N e i=1 n k=0 σ 2 b i E T r(D i (k)D i t (k))T r A i (k)A i t (k) + N e i=1 N e j=1 m b i m b j n k=0 n m=0 E T r A i (k)1 N D i (k)D j t (m)1 N t A j t (m) (29)
the statistical mean is replaced by a time average over Nt samples to get a realistic result. In practice, a number Nt equal to 100 leads to accurate modelling properties, as it will be shown in Section IV.

The complexity of the proposed method has been analyzed and leads to a number of multiplications N mult and additions N add equal to

N mult = Nt.p i (2MyMiNi + MyNy(Mi + Ni)) (33)
N add = Nt.p i (2MyMiNi + MyNy(Mi + Ni) + MyNy) (34)
where MyxNy corresponds to the output noise matrix By(n) size and MixNi those of input noise B i (n). In the next paragraph, an approach to reduce the complexity of the infinite sum computation without a significant loss of accuracy is presented.

2) Linear Prediction Approach: To reduce complexity, a linear prediction approach is proposed to replace the infinite sums. [START_REF] Lopez | Fast and accurate computation of the roundoff noise of linear time-invariant systems[END_REF] has shown that impulse response terms hi(k) are defined by a recurrent equation where the terms Pi corresponds to the maximal delay for the recursive part. For LTI systems, terms f and g are directly equal to transfer function coefficients. Therefore, the rest of the section will focus on non-LTI recursive systems.

As the impulse response is time-varying, the relation ( 22) is non linear. The aim is then to linearize this expression with prediction coefficients λi to obtain an estimate hi(k) of the impulse response hi(k) equal to

hi(k) = P i m=1 λi m hi(k -m). ( 35
)
Computing the coefficient vector λi = [λi 1 ...λi P i ] t gives a model of the impulse response term hi(Pi+Qi-1) with the aim to minimize the mean square error ei between hi(Pi + Qi -1) and its estimate value hi(Pi + Qi -1) defined as

E[e 2 i ] = E[(hi(Pi + Qi -1) -hi(Pi + Qi -1)) 2 ]. (36)
For the sake of clarity on ei, a term k = Pi + Qi -1 is introduced and the index i is removed. Thus, the mean square error is defined as

E[e 2 ] = E[(h(k) -h(k)) 2 ]. ( 37
)
The derivative of the mean square error, with respect to the coefficients λj, ∀j ∈ [1, P ], is

∂E[e 2 ] ∂λj = -2E h(k)h(k -j) - m λmh(k -m)h(k -j) . ( 38
)
These P prediction coefficients λj can be determined by setting (38) to zero, leading to the solution λopt defined as

λopt = U -1 J, (39)
with

U =        E[h 2 (k -1)] . . . E[h(k -1)h(k -P )] E[h(k -1)h(k -2)] . . . E[h(k -2)h(k -P )] . . . E[h(k -1)h(k -P + 1)] . . . E[h(k -P + 1)h(k -P )] E[h(k -1)h(k -P )] . . . E[h 2 (k -P )]        (40)
and

J =       E[h(k)h(k -1)] E[h(k)h(k -2)] . . . E[h(k)h(k -P + 1)] E[h(k)h(k -P )]       (41)
To determine J and U, the first values of the impulse response must be known. Therefore, the approach to determine the prediction coefficients is as follows:

1) The first P +Q-1 terms of the time-varying impulse response are determined using [START_REF] Lopez | Fast and accurate computation of the roundoff noise of linear time-invariant systems[END_REF], with Q and P representing respectively the maximal number of delays for the non-recursive and the recursive part of the system S.

2) The two matrices U et J are determined.

3) The prediction coefficients λopt = [λ1....λP ] are determined by λopt = U -1 J. The coefficients λopt are determined according to Appendices A and B and the infinite sums can therefore be written as

∞ m=0 h(m) = 1 P P j=1 Q-2 m=0 ωm + (I P -Λ) -1 ωQ-1 + θ (j) (42) ∞ m=0 h(m) 2 = 1 P T r Q-2 m=0 ωmω t m + S + Φ (43)
where

ωm =      h(m) h(m + 1)
. . . (46)

h(m + P -1)      (44) Λ =        0
Φr,r = r-2 k=0 h(k) 2 (47)
and with matrix S being the solution of

S -ΛSΛ t = ωQ-1ω t Q-1 . (48)
In practice, this model is applied to matrices A and D for each noise source. The complexity of this method is computed and leads to a number of multiplications N mult and additions N add equal to

N mult = i Nt (Pi + Qi -1) 2MyMiNi + MyNy(Mi + Ni) + Pi(My + Ny)(Mi + Ni) + P 2 i (49)
N add = i Nt (Pi + Qi -1) 2MyMiNi + MyNy(Mi + Ni) + 4Pi (50)

D. Accuracy Evaluation Tool

The proposed approach has been implemented in a software tool to automate this process, its synoptic is given in Figure 2. Our numerical accuracy evaluation tool generates the analytical expression of the output quantization noise from the signal flow graph (SFG) of the application. This analytical expression is implemented through a C function having the fractional part w FP of all data word-lengths as input parameters. This C code can be compiled and dynamically linked to the fixed-point conversion tool for the optimization process. The accuracy evaluation tool uses a technique similar to the one proposed in [START_REF] Menard | Analytical fixed-point accuracy evaluation in linear time-invariant systems[END_REF]. The first step transforms the application SFG Gs into an SFG Gns representing the application at the quantization noise level. All the potential noise sources are included in the graph. The expression of the variance σ 2 b i and the mean m b i of each noise source bi according to w FP are generated in order to be included in the output C function. The expression of the number of bits eliminated between two operations can be found in [START_REF] Menard | Analytical fixed-point accuracy evaluation in linear time-invariant systems[END_REF]. Then, each operation is replaced by its propagation noise model, as presented in Section III-A2. The aim of the second step is to determine the expression of the system impulse response between the output and the different noise sources. To handle recursive systems, cycles in the graph are detected and then dismantled. The third step computes the value of the gains Ki and Lij given in [START_REF] Shi | A perturbation theory on statistical quantization effects in fixed-point DSP with non-stationary inputs[END_REF] and [START_REF] Sripad | A Necessary and Sufficient Condition for Quantization Error to be Uniform and White[END_REF]. The linear prediction approach presented in Section III-C2 is used to approximate the infinite sums. The values of the different terms of the matrices A and D are collected during a single floating-point simulation. The reference simulation and the third step of the tool are performed with the Matlab tool.

IV. RESULTS

In this section, experimentations are carried out to validate the proposed approach on LTI and non-LTI systems. Results are given on various applications: vector normalization illustrates non-LTI non-recursive systems, adaptive filters such as LMS, NLMS and APA illustrate non-LTI recursive systems, and various benchmarks illustrate all kinds of systems. For these applications the Signal-to-Quantization Noise Ratio (SQNR) values are chosen from 30 dB to 90 dB to ensure that the assumptions of the Widrow model are satisfied (the noise variance is significantly smaller than the signal power). Moreover, the optimization time of the proposed approach is compared with approaches based on fixed-point simulations.

A. A non-LTI non-recursive system: Vector Normalization

The first example considers an N -size vector x(n) normalized according to the relation

y(n) = x(n) x t (n)x(n) . (51)
Here to simplify the final expression of the power, only two noises are considered: the noise bx(n) due to the quantization of x(n) and the noise bg(n) generated by the normalization computation. By using the noise model propagation defined in Section III-A, the output noise vector by(n) is equal to

by(n) = bx(n) x t (n)x(n) + 2 x(n)x t (n)bx(n) x t (n)x(n) 2 + bg(n), (52)
which can be simplified as

by(n) = I N x t (n)x(n) + 2 x(n)x t (n) x t (n)x(n) 2 Ax(n) bx(n) + bg(n), (53)
with I N the N -size identity matrix. This example has no delay in the graph so only one impulse term Ax(n) has to be determined. The other terms Ax(k) are then equal to zero for k = n and dx(n) is equal to one. By applying [START_REF] Ozeki | An adaptative filtering algorithm using an orthogonal projection to an affine subspace and its properties[END_REF] to the output noise expression, its power P b is where mx and mg are the means of bx(n) and bg(n) and σ 2 their variances.

P b = E[T r(Ax t (n)Ax(n))]σ 2 x + N σ 2 g + m 2 x E[T r(Ax t (n)Ax(n))] + N m 2 g + 2mgmxE[T r(Ax(n))] (54)
For this example, N is set to 4, and the truncation mode is considered. The fractional part word-length of the input x and for the output y are respectively win and wout. So, using [START_REF] Roy | An algorithm for trading off quantization error with hardware resources for Matlab-based FPGA design[END_REF] and by introducing values of mean and variance of the two considered noises, (54) can be written as :

P b = K1 2 -2w in 12 + K2 2 -2w out 12 (1 -2 -2(k out ) ) + L11 2 -2w in 4 + L22 2 -2w out 4 (1 -2 -(k out ) ) 2 + L12 2 -w out 2 (1 -2 -(k out ) ) 2 -w in 2 + L21 2 -w out 2 (1 -2 -(k out ) ) 2 -w in 2 (55)
where kout represents the number of bits eliminated at the output of the normalization. The expression of the gain K and L are given in Table IV.

The chosen synthetic input signal x(n) is a first-order autoregressive signal defined by

x(n + 1) = c.x(n) + u(n), (56)
with c ∈ [0, 1] and u(n) a zero-mean white noise of variance 0.01. The numerical values of the terms K and L are given in Table IV. As shown by the results, the terms K and L depend on signal characteristics and correlation between them.

B. A non-LTI recursive system: the Least Mean Square (LMS) adaptive filter 1) Fixed-Point LMS Algorithm:

To illustrate the proposed model and to verify its accuracy on a non-LTI and recursive system, the Least Mean Square (LMS) example is detailed. The LMS adaptive algorithm [START_REF] Haykin | Adaptive Filter Theory[END_REF] addresses the problem of estimating a sequence of scalars y(n) from an N -length vector x(n) = [x(n), x(n -1) . . . x(n -N + 1)] t . The linear estimate of y(n) is w t (n)x(n) where w(n) is an N -length vector which converges to the optimal vector wopt in the mean-square error (MSE) sense. This optimal vector is equal to wopt = R -1 v where R corresponds to the correlation matrix of the input signal equal to E[x(n)x t (n)] and where v represents the intercorrelation vector between x(n) and y(n) equal to E[x(n)y(n)]. The vector w(n) is updated according to

w(n + 1) = w(n) + µLMS.x(n)(y(n) -w t (n)x(n)) (57)
where µLMS is a positive constant representing the adaptation step.

In a fixed-point implementation of the LMS, four equivalent noise sources have to be considered as presented in Figure 3. bx(n) and br(n) are noises associated respectively with the input data x(n) quantization and the desired signal y(n) quantization. The vector bw(n) is the noise vector due to the quantization of the products between µ, x(n) and the error e(n) = y(n)w t (n)x(n). Finally, b f (n) groups together all noises generated inside the filtering part during the computation of the inner product w t (n)x(n). The terms µi and σ 2 i represent the mean and variance of each noise source bi. 2) Noise Model for Adaptive Coefficient Computation: In this section, the output noise power expression is determined from the four noise sources: bx(n), b f (n), bw(n) and br(n). For each one, the impulse response between the source and the output is determined. Then, the output noise power can be directly expressed using [START_REF] Roy | An algorithm for trading off quantization error with hardware resources for Matlab-based FPGA design[END_REF] where Ne = 4 in that case. Also, the output noise power can be determined using the linear prediction approach (as presented in Section III-C2). Each impulse response is modelled using the prediction coefficients, leading to a simpler expression of the output noise power. These two approaches are presented on the LMS example in the following. Only the propagation of the noise bw(n) is detailed. The contribution of the other noises can be obtained in the same way. The propagation of noise term bw(n) is shown on Figure 4 and is considered for validating the accuracy of the proposed approach and detailed hereafter. As shown on Figure 3, the noise bw(n) is first delayed with a first recursion. Then, the noise goes through the FIR filter where it is multiplied by the signal x t (n). After, the noise propagates through the second recursion to the system output and comes back to the multiplication by µLMS.x(n).

FIR Q + + - Z -1 Adaptation b x (n) Q Q + Q b r (n) b w (n) b f (n)
Noise propagation can be modelled using the following set of recurrent equations

b ′ w (n) = x t (n)b1w(n -1) (58) b1w(n) = b1w(n -1) + bw(n) -µx(n)b ′ w (n) (59)
where b1w(n) is the noise after the delay. By developing these two previous expressions, the scalar output noise b ′ w (n) is written using its time-varying impulse response hw as The time-varying impulse response hw is defined through two equivalent matrices as defined in Section III. For this noise, the matrix Aw corresponds to an N -size vector aw. This vector is introduced in (61) and is defined in (62).

b ′ w (n) = n-1 k=0 hw(k)bw(k). ( 60
) + Z -1 x x b w '(n) b w (n) b1 w (n) x t (n) µ LMS .x(n)
b ′ w (n) = n-1 k=0 aw(k)bw(k) (61) aw(k) = x t (n)F(n, k) (62)
with

F(n, k) = n-1 m=k+1 IN -µLMSx(m)x t (m) . (63)
The expression of aw can be modelled using linear prediction coefficient as presented in Section III-C2. The first impulse response terms of aw(k) are

aw(n) = 0, aw(n -1) = x t (n), aw(n -2) = x t (n) I N -µ LM S x(n -1)x t (n -1) . (64)
The LMS algorithm is a first-order recursive system. Thus, only one prediction coefficient has to be determined using the approach defined before. The matrix Uw, defined by (40), is a scalar in that case and is equal to

Uw = E(aw(n -1)aw t (n -1)) = E x t (n)x(n) . ( 65
)
The term Jw, defined by (41), is equal to

Jw = E(aw(n -2)aw t (n -1)) = E x t (n) I N -µ LM S x(n -1)x t (n -1) x(n) . ( 66
)
So, from (39), the prediction coefficient λw is expressed as

λw = Jw Uw = E[aw(n -2)aw t (n -1)] E[aw(n -1)aw t (n -1)] = E[x t (n)(IN -x(n -1)x t (n -1))x(n)] E[x t (n)x(n)] . (67)
Thus, with this coefficient, the impulse response term aw can be expressed as

aw(k) = λwaw(k + 1), (68)
for k < n -1 since aw(n) = 0. The infinite sum is obtained using (42) as

n→∞ k=0 aw(k) = 1 1 -λw aw(n -1). ( 69
)
The contributions of the three other terms can be obtained with the same method. The time-varying impulse response associated with these three noises are given in the following expressions.

ax(k) = µLMSx t (n)F(n, k) e(k) -x(k)w t (k) + w t (n)δ(n -k) ar(k) = µLMSx t (n)F(n, k)x(k) a f (k) = -µLMSx t (n)F(n, k)x(k) + δ(n -k) (70)
where δ(k) is the Kronecker operator. It is noteworthy that, as a consequence, the noises b f (n) and br(n) are scalar. Then, the terms A modelling their propagation through the system are also scalars, which let us write T r(AA t ) = A 2 for input noises b f (n) and br(n).

The output noise power is finally computed using [START_REF] Roy | An algorithm for trading off quantization error with hardware resources for Matlab-based FPGA design[END_REF] and leads to (71) given at page 10. The noise mean term M is equal to

M (k) = ax(k)mx + ar(k)mr + a f (k)m f + aw(k)mw
where the mean values mx and mw are N -size vectors and mr and mr are scalars. The complexity of the computation of the gains Ki and Lij of (71) are obtained from ( 33) and [START_REF] Wadekar | Accuracy Sensitive Word-Length Selection for Algorithm Optimization[END_REF]. In this example, input noises are N size vectors for bw and bx (so Mi = N and Ni = 1) and scalar for b f and br ( Mi = Ni = 1), the output noise by is a scalar (My = Ny = 1). Moreover, the number p is set to 500 and Nt = 100. The number of multiplications N mult is equal to

N mult = 50000(6N + 8) (72)
So the complexity is in O(N ). The expression of the output noise power can be modelled using prediction coefficients that can be computed for noises bx(n), b f (n) and br(n) using the same method as the one presented before for bw(n). Using the prediction coefficients, the expression (71) is simplified and leads to (73) given at page 10 in which σ 2 refers to scalar variance value of each noise.

By applying (49), the number of multiplications is then equal to

N mult = 100(8N + 16) (75)
The complexity is divided by a factor 375. Thus, the output noise power can be computed by two different ways: using the impulse response terms of ((71) derived from (30)) or using prediction coefficients using (73) presented in Section III-C2. The quality of these two approaches is compared in the next section.

3) Quality of the Estimation: To analyze the quality of the proposed approach, experiments have been performed. This quality has been evaluated through the measurement of the relative error Er between our analytical estimation P b and the estimation P (sim) b based on fixed-point simulations which is considered to be the reference. The relative error Er is defined as

Er = P b -P (sim) b P (sim) b . (76)
The reference value is obtained by simulations from the difference between a floating-point and a fixed-point simulation. The floating-point simulation is used as the high-precision reference, since floating-point arithmetic generates negligible noises compared to fixed-point arithmetic.

For these experiments, the chosen synthetic input signal x(n) is a first-order autoregressive signal defined by (56). Figure 5 shows the relative error of the proposed approach on a 32-tap LMS adaptive filter. The results are presented according to the number of elements p chosen to compute the infinite sums and the correlation coefficient

E[b 2 y (n)] = n k=0 σ 2 x E T r ax(k)ax t (k) + n k=0 σ 2 f E[a 2 f (k)] + n k=0 σ 2 r E[a 2 r (k)] + n k=0 σ 2 w E T r aw(k)aw t (k) + n k=0 n l=0 E T r M (k)M t (l) (71) E[b 2 y (n)] = σ 2 x E T r(ax(n)ax t (n)) 1 -λ 2 x + σ 2 f E[a 2 f (n)] 1 -λ 2 f + σ 2 br E[a 2 br (n)] 1 -λr σ 2 w E T r(aw(n)aw t (n)) 1 -λ 2 w + T r T T t (73)
with

T = mx E[ax(n)] 1 -λx + mr E[ar(n)] 1 -λr + m f E[a f (n)] 1 -λ f + mw E[aw(n)] 1 -λw . ( 74
)
c of the input signal x(n). The input signal can be a white noise (c = 0), a fairly correlated noise (c = 0.5) or a very correlated noise (c = 0.95). As p increases, the relative error decreases. Indeed, when p increases, more terms are included in the computation of the sums, which leads to a better estimation quality. Moreover, the relative error convergence speed depends on input data correlation.

For non-correlated input data, the relative error convergence is slower than the one for very-correlated input data. In fact, the relative error is less than 20% after 300 samples for very correlated input data, after 350 samples for fairly correlated data and after 550 samples for uncorrelated input data. To explain this result, we recall that the terms A, representing the noise propagation through the system, include the terms F(n, k) defined as

F(n, k) = n-1 m=k+1 IN -µLMSx(m)x t (m) . (77)
For very correlated input data, the term IN -µLMSx(m)x t (m) decreases faster than for uncorrelated data. Thus, the number of elements p required to compute the infinite sums depends on the input data correlation. Nevertheless, with experimentation presented for p > 500, the relative error is less than 25% in all cases, which represents a difference equivalent to less than 1 bit between the noise power obtained with the proposed approach and the real noise power. For the linear prediction approach, the obtained relative error is equal to 21%. For other sizes N of the LMS algorithm, same kinds of results are obtained, thus indicating that the proposed approach is accurate for the LMS algorithm.

Experiments have also been conducted on other adaptive algorithms. The results are presented in Table V for the previous LMS, a 32-tap Normalized LMS (NLMS) and an Affine Projection Algorithm (APA) of size 32 × 12. The expression of the coefficient update for the APA algorithm is given in [START_REF] Eyre | The evolution of DSP processors[END_REF]. The APA application is suitable to show the case where the general model presented in Section III include the two matrices A and D. Both NLMS and APA algorithms moreover include division operations for the normalization module. The floating-point simulation is executed on Nt + p samples where Nt is the number of samples to compute the statistical mean and p is the number of terms to compute the sums. In these experiments, Nt is set to 100 and different values of p are tested.

The approach based on linear prediction leads to a relative error of around 20% (equal to 1 dB). However, the method using the computation of the sums can be really accurate since the relative error is less than 17.65% in mean with p = 100 and less than 6.24% with p = 1000. The accuracy of the estimation can be controlled by adjusting the number of elements p used to compute these sums. We can therefore see the term p as a way to trade-off between the accuracy of the estimation and the time required to compute this estimation.

C. Experiments on LTI and Non-LTI Non-Recursive Systems

In this section, the proposed approach is first evaluated on LTI systems such as FIR and IIR filters, an MP3 coder/decoder and an FFT. Then, the approach is applied to non-LTI and non-recursive systems (square, Volterra filter, correlator). The average and maximum relative errors between the noise power estimated with the proposed approach and the noise power obtained by fixed-point simulations are provided. VI. For the 32-tap FIR filter, rounding and truncation have been tested. The average relative error is 2.49% for rounding and 1.14% for truncating. Moreover, the maximum relative errors are very low (about 8%).

For the MP3 coder/decoder, the average error is equal to 5.24% for the polyphase filter and to 8.42% for the DCT. The maximum errors are 20.57% for the polyphase filter and 24.8% for the DCT. These values represent a difference of 2dB (less than 1 bit) between the real value of the noise power and the one estimated by the proposed approach. For the FFT the maximal relative error is less than 7.89%, which also illustrates the accuracy of the proposed approach on relevant applications such as the FFT. To test the method on a recursive LTI system, we have considered an 8-order IIR filter divided into four 2-order cells. The relative error is less than 3.3% with an average value of 1.68%.

2) Non-LTI and Non-Recursive systems: For the case of non-LTI systems we have tested three applications: a 2 nd -order polynomial evaluation, a 2 nd -order Volterra filter and a correlator. The results are presented in Table VII. For the 2 nd -order polynomial, the estimation is highly accurate, since the error is lower than 0.8%. For the Volterra filter and the correlator, the average relative error is equal to 1.79% and 1.35% with the maximum values of 3.22% and 5.78% respectively. For these three applications, the error is less than 5.78% , thus confirming the accuracy of the proposed approach.

In summary, for all types of systems based on arithmetic operations, the presented results tend to indicate that the method has an overall relative error that is less than 25% even in the worst cases, which can be translated to less than 1 bit of precision.

D. Fixed-Point Specification Optimization Time

The evaluation of the numerical accuracy is used in the wordlength optimization based on an iterative process. The optimization time obtained with the proposed approach is compared with those obtained with fixed-point simulations. For the two approaches, the optimization times are given in Figure 6 for the 32-tap LMS algorithm. A computer with Intel Core 2 Duo T8300 at 2.4 GHz and 2GB of RAM was used for the experiments. For the approach based on fixed-point simulation, a new simulation is required as long as a fixed-point format is modified. So, the optimization time depends linearly on the number of optimization iterations. For approaches based on fixed-point simulations, experimentations were conducted on Matlab to evaluate their computing time. For the proposed approach, the analytical expression of P b needs to be computed first. The determination of this analytical expression is the most time-consuming part and its computation takes 46 seconds for the approach based on infinite sum approximation and 4 seconds for the linear prediction approach. Then, each iteration of this optimization process corresponds only to the noise power expression evaluation, whose computing time is negligible. For the LMS algorithm, the proposed approach leads to a gain after only 100 iterations, which represents an execution time equal to 46 seconds. As an example, for an optimization process with about 30 variables, between 10 3 and 10 5 iterations are required. In that case the proposed approach would require between 10 and 10 3 less time to perform the optimization process. The approach based on linear prediction leads to time gain after only 10 iterations compared to approaches based on fixed-point simulations, thus leading to an execution time equal to 4 seconds. These results clearly show that the proposed approach is faster than the best simulation-based techniques for the process of fixed-point optimization.

A comparison with methods from [START_REF] Caffarena | SQNR Estimation of Fixed-Point DSP Algorithms[END_REF] and [START_REF] Shi | A perturbation theory on statistical quantization effects in fixed-point DSP with non-stationary inputs[END_REF] has also been carried-out. As they are also analytical, only the time to determine the coefficients Ki and Lij is given. For the approach developed in [START_REF] Shi | A perturbation theory on statistical quantization effects in fixed-point DSP with non-stationary inputs[END_REF], no result on the quality and the execution time were given. However, for an application with Ne noise sources, N 2 e fixed-point simulations are required. Thus, the execution time grows quadratically with the complexity of the application. For the 32-tap LMS algorithm presented above, 66 noise sources are considered and thus the approaches proposed in [START_REF] Shi | A perturbation theory on statistical quantization effects in fixed-point DSP with non-stationary inputs[END_REF] would require 4356 fixed-point simulations. Therefore, the proposed method leads to significantly lower optimization time, especially when the application becomes complex.

In the approach based on affine arithmetic simulations proposed in [START_REF] Caffarena | SQNR Estimation of Fixed-Point DSP Algorithms[END_REF], the relative error Er for the estimation and the time Texe required to determine the coefficients Ki and Lij are given in Table VIII for a second-order IIR filter (IIR2), an 8-order IDCT, and for an LMS filter with 2 and 5 taps. For the case of the IIR2 application, the time to execute the third step corresponding to the determination of the noise power expression (computation of Ki and Lij) corresponds to 70 % of the global execution time for the approach based on infinite sums and to 50 % for the approach based on is dedicated to the other steps (front-end, determination of system noise model and expression generation). For the proposed method based on infinite sums with p = 1000, the quality of the estimation is close to the one obtained with the technique of [START_REF] Caffarena | SQNR Estimation of Fixed-Point DSP Algorithms[END_REF]. For our linear prediction based-approach, the quality is slightly lower, but the execution time is reduced. For our two approaches, the execution time is significantly reduced compared to the approach proposed in [START_REF] Caffarena | SQNR Estimation of Fixed-Point DSP Algorithms[END_REF] and a reduction of several order of magnitude is observed.

The results show the benefit of the proposed approaches to reduce the development time of fixed-point systems. We propose two different approaches with different characteristics, as shown in Table VIII. It is seen that the infinite sums approach is really accurate and fast, whereas, the approach based on linear prediction is faster but a little less accurate. So, the proposed framework allows the user to choose between two different approaches with different trade-offs between speed and accuracy. In this paper, a fast general and accurate approach to determine the accuracy of a fixed-point system is presented and validated. The method is analytic and provides a closed-form expression of the output noise power for systems composed of smooth operations. The quantization noises due to fixed-point computation are considered through their statistical properties (mean, variance, etc.) and their propagation through arithmetic operations (addition, subtraction, multiplication, and division) is detailed. The system is characterized by its time-varying impulse response and the expression of the output noise power is obtained. The approach is valid for all types of system made-up of smooth operations and for most common rounding modes (truncation, conventional rounding, convergent rounding). The estimation is unbiased and leads to infinite sums inside the expression of the output noise power which have to be approximated by a finite number of terms. To reduce the complexity due to infinite sums, an approach based on linear prediction is also proposed.

Number of iterations

The proposed approach was applied to the LMS algorithm, a recursive and non-LTI system, to verify its validity and its quality. Moreover, several LTI and non-LTI systems have been tested to show the accuracy of the proposed approach on various types of applications. The maximal relative errors are less than 25% which corresponds to less than 1 bit. Finally, the approach has been compared to approaches based on fixed-point simulations in terms of execution time. After only 10 iterations, the proposed approach is less time-consuming than other approaches based on fixed-point simulations. Thus, it was shown that the proposed approach leads to accurate estimation and allows a significant reduction of the optimization time in the fixed-point conversion process. . . .

h(m + P -1)      . ( 78
)
The first P + Q -1 terms of impulse response are directly determined. Thus, terms h(0) to h(P + Q -2) are computed with [START_REF] Lopez | Fast and accurate computation of the roundoff noise of linear time-invariant systems[END_REF]. The term P + Q -1 is present in the matrix ΩQ and in the following P -1 matrices. This term can also be computed with prediction coefficient λi. The P -size matrix Λ is defined by (81)

Λ =       
To compute the determinant, we develop according to the first column, which leads to the recurrence expression ψP (X) = XψP -1(X ) -λP (-1 P -1 )(-1 P -1 ) = XψP -1(X ) -λP As for P = 1, Λ = λ1, ψ1(X) = X -λ1 and by using the previous recursion, ψP (X) is equal to ψP (X) = X P -P k=1 λ k X P -k .

(82)

The roots of the previous equation are the eigenvalues of Λ. Moreover, these prediction coefficients define a transfer function. Indeed, they linearize the time-varying transfer function. The poles ρm of the transfer function satisfy the equality

where T is an invertible matrix and ρm can be real or complex. However, since the system is stable, the modulus of ρm is less than one so ρ n m converges to 0 when n → ∞. So, Λ n is equal to

Λ n = T       
(85) proves that Λ n tends to zero when n → ∞. Therefore:

∞ m=0 ωm = Q-2 m=0 ωm + ∞ k=0 Λ k ω Q-1 = Q-2 m=0 ωm + (I P -Λ) -1 ω Q-1 (86)
with I P the P -size identity matrix.

The vector

∞ m=0
ωm is also equal to

∞ m=0 ωm =            ∞ m=0 h(m) ∞ m=1 h(m)
. . . This expression is developed for impulse terms h. In practice, the impulse terms are defined by two matrices A and D as shown in (Section III). The previous method is then applied to A1D where 1 is a matrix containing ones (and replacing the noise B).

APPENDIX B

In this appendix, the method modelling the second order infinite sums is presented. Using the same approach as in Appendix A, the following recurrence relation is obtained:

ωm+1ω t m+1 = Λωmω t m Λ t = Λ m+2-Q ωQ-1ω t Q-1 Λ t m+2-Q . ( 90
) Matrix Λ m has a limit equal to zero when m → ∞ as demonstrated previously, which gives

∞ m=0 ωmω t m = Q-2 m=0 ωmω t m + ∞ k=0 Λ k ωQ-1ω t Q-1 Λ t k S . (91)
Yet, the next equality can be determined:

S -ΛSΛ t = ωQ-1ω t Q-1 . (92)
This equality is a Lyapunov equation. So, S is the solution of the Lyapunov equation and can be obtained using a classical solver. As in Appendix A, we introduce the diagonal matrix Φ where each diagonal element is equal to Φr,r =

Fig. 1 .

 1 Fig. 1. System noise model for Ne input noises b i (n) and output noise by(n)

1 )

 1 System Characterization: The contribution b ′ i (n) of the noise source bi(n) depends on the previous samples bi(nk) of the noise source with k ∈ {1, 2, ..Qi}. If the system Si is recursive, the noise b ′ i (n) depends on its previous samples b ′ i (nm) with m ∈ {1, 2, ..Pi}. Thus, the general expression of the quantization noise b

Fig. 2 .

 2 Fig. 2. Synoptic of the Accuracy Evaluation Tool

Fig. 3 .

 3 Fig. 3. Fixed-point LMS algorithm and noise model for input data x ′ (n), estimate data y ′ (n) and coefficients w ′ (n)

Fig. 4 .

 4 Fig. 4. Propagation through the LMS algorithm of the noise source bw(n) due to filter coefficient computing

Fig. 5 .

 5 Fig. 5. Relative error for the 32-tap LMS according to the number of samples p in the computation of the infinite sums

Fig. 6 .

 6 Fig.6. Optimization time (s) for the LMS32 example of our two approaches (linear prediction, infinite sums with p = 500) and the approach based on fixed-point simulations

  APPENDIX AIn this section, infinite sums are modelled for impulse response terms h(k). The vector ωm is defined as ωm =

Theorem 2 :

 2 recurrence relation is obtained:ωm+1 = Λωm = Λ m+2-Q ωQ-1.(80) The matrix Λ n has a limit equal to 0 when n → ∞.Proof:The characteristic polynomial ψP (X) = det(XI P -Λ) of the matrix Λ is equal to ψP (X) = λP -1 . . . . . . X -λ1 .

  solutions of (82) with (83), eigenvalues of Λ are the poles of the transfer function, and Λ can be written as

  vector θ is introduced where each element θm is defined byθm = m-2 k=0 h(k)(88)θ contains the P -1 first impulse response terms which are computed manually. Thus, each line of I P -Λ) -1 ωQ-1 + θ (j) . (89)

TABLE I FIRST

 I

-AND SECOND-ORDER MOMENTS FOR THE THREE CONSIDERED QUANTIZATION MODES.

TABLE IV EXPRESSION

 IV AND VALUES OF THE GAINS K AND L IN THE CASE OF THE

	VECTOR NORMALIZATION EXAMPLE

TABLE VI RELATIVE

 VI ERROR FOR LTI SYSTEMS 1) LTI Systems: First, non-recursive systems, such as FIR filter, MP3 coder/decoder and FFT, are evaluated for different fixed-point specifications. The results are presented in Table

TABLE VIII COMPARISON

 VIII OF THE PROPOSED APPROACHES WITH THE ONE PROPOSED IN [4] IN TERMS OF EXECUTION TIME Texe (SECONDS) AND RELATIVE ERROR Er FOR THE ESTIMATION linear-prediction. The time to execute the second step corresponding to the determination of system expression (computation of Ai and Di) represents to 26 % of the global execution time for the approach based on infinite sums. For this approach, the rest of the time (4 %)

Copyright (c)

IEEE. Personal use of this material is permitted. However, permission to use this material for any other purposes must be obtained from the IEEE by sending an email to pubs-permissions@ieee.org.


Member, IEEE Romuald Rocher

Member, IEEE Pascal Scalart

Noise Probability Density Function in Fixed-Point Systems based on smooth operators

To satisfy cost constraints, application implementation in embedded systems requires fixed point arithmetic. Thus, the application defined in floating point arithmetic must be converted into a fixedpoint specification. This conversion requires accuracy evaluation to ensure algorithm integrity. Indeed, fixed-point arithmetic generates quantization noises due to the elimination of some bits during a cast operation. These noises propagate through the system and degrade computing accuracy. In this paper, a method based on Generalized Gaussian PDF is presented to model and generate the output noise of the system. The accuracy of the proposed model is evaluated through different experiments.

I. INTRODUCTION

Digital signal processing applications are specified in floating point to prevent problems due to computing accuracy. However, to satisfy cost constraints, application implementation in embedded systems requires fixed point arithmetic. Thus, the application defined in floating point arithmetic must be converted into a fixed-point specification. To reduce application time-to-market, tools to automate floating-point to fixed-point conversion are needed. In these tools, an important stage corresponds to accuracy evaluation of fixed-point specification. Indeed, fixed-point arithmetic generates quantization noises due to the elimination of some bits during a cast operation. These noises propagate through the system and modify computing accuracy. Computing accuracy deviations must be limited to ensure algorithm integrity and application performance.

Application accuracy can be evaluated by different manners. On one hand, accuracy evaluation can be obtained with fixed-point simulations [START_REF] Belanovic | Automated Floating-point to Fixed-point Conversion with the fixify Environment[END_REF], [START_REF] Kim | Fixed-Point Optimization Utility for C and C++ Based Digital Signal Processing Programs[END_REF]. However, these methods require high computing time since a new simulation is required as soon as a format changes in the system. So, these approaches lead to very important optimisation time inside a fixed-point conversion process. On the other hand, a fixed-point specification accuracy can be evaluated with analytical methods [START_REF] Rocher | Analytical approach for numerical accuracy estimation of fixed-point systems based on smooth operations[END_REF], [START_REF] Chesneaux | The fixed cadna library[END_REF], [START_REF] Caffarena | SQNR Estimation of Fixed-Point DSP Algorithms[END_REF]. These approaches determine a mathematical expression for the accuracy metric. In this paper, we present a method to study the fixed-point behaviour of a system with a floating-point simulation. Indeed, the fixed-point system can be represented by the floating-point system where a global noise is added at the output. This global noise must represent the quantization noise contributions due to fixed-point arithmetic whose characteristics are analytically determined. So, we need to be able to generate samples following the characteristics (PDF, power) of this global noise. In this way, we consider Generalized Gaussian law which allows us to generate samples following different PDFs (uniform, dirac, gaussian...). The parameters of the Generalized Gaussian are analytically determined from statistics (mean, variance and kurtosis) of the output noise. This paper is organized as follows. We first introduce quantization noise models. The model of a noise source generation is presented and the propagation of these noises through the system is developed. The considered system is supposed to be composed by smooth operations such as additions, subtractions, multiplications and divisions. A general PDF based on Generalized Gaussian is presented and the approach to compute its parameters is defined. We also propose a method to generate samples following this PDF.

Then, the High Order Moments (HOM) are presented. They are necessary to determine the parameters of the Generalized Gaussian. First, the 3 rd (skewness) and 4 th (kurtosis) order moments of the quantization noises are computed for the different quantization types (truncation, rounding and convergent rounding). The 1 st and 2 nd order moments are already known. Then, the computation of the HOM of the output noise are analytically computed from the expression of the output noise.

Finally, the method is applied to different systems such as the FIR and IIR filters and Volterra filters. These results show the accuracy of our model for different noise PDFs.

II. QUANTIZATION NOISE MODELS

A. Quantization Noise Sources

The quantization process can be modelled by the sum of the original signal and of a uniformly distributed white noise [START_REF] Widrow | Statistical Theory of Quantization[END_REF], [START_REF] Sripad | A Necessary and Sufficient Condition for Quantization Error to be Uniform and White[END_REF]. This quantization noise is uncorrelated with the signal and the other noise sources. Such a model is valid provided that the signal has a sufficiently large dynamic range compared to the quantization step. According to the quantization mode, the noise Probability Density Function (PDF) will differ. Three quantization modes are usually considered: truncation, conventional rounding, and convergent rounding. In the truncation mode, the Least Significant Bits (LSB) are directly eliminated. The resulting number is always smaller than or equal to the value available before quantization, and, therefore, the quantization noise is always positive. Consequently, the mean of the quantization noise is not equal to zero. To reduce the bias due to truncation, the rounding quantization mode is often used. In conventional rounding, the data are rounded to the nearest value representable in the reduced-accuracy format. For numbers located at the midpoint between two consecutive representable values, the data are rounded-up always to the larger output value. This technique leads to a (small) bias for the quantization noise. To eliminate this quantization noise bias, the convergent rounding can be used as well. In this case, the numbers located at the midpoint between two consecutive representable values are, with equal probability, rounded to the higher or lower output value.

Let n denotes the number of bits for the fractional part after the quantization process and k the number of bits eliminated during the quantization. The quantization step q after the quantization is equal to q = 2 -n . The quantization noise mean and variance are given in Table I for the three considered quantization modes [START_REF] Constantinides | Truncation Noise in Fixed-Point SFGs[END_REF], [START_REF] Menard | Quantization Mode Opportunities in Fixed-Point System Design[END_REF].

Quantization Truncation Conventional Convergent mode rounding rounding Mean q 2 (1 -2 -k ) q 2 (2 -k ) 0 Variance q 2 12 (1 -2 -2k ) q 2 12 (1 -2 -2k ) q 2 12 (1 + 2 -2k+1 )

B. Output Noise Model

The system output noise b y (n) is the sum of all contributions as presented in Figure 1 and expressed by [START_REF] Rocher | Analytical approach for numerical accuracy estimation of fixed-point systems based on smooth operations[END_REF] b where h i (k) is the time varying impulse response of system S i between the noise source b i (n) and the output noise b y (n). So, as shown in Figure 2, a fixed point system can be represented by a floating-point system whose output is perturbed (in an additive way) by a noise term by(n) modelling all quantization noise contributions.

y (n) = Ne i=1 n k=0 h i (k)b i (n -k) (1)
+ b 1 (n) b i (n) b N e (n) b 1 (n) b i (n) b N e (n) b y (n) S 1 S i S N e
Since fixed-point simulations are very time consuming, it is really interesting to be able to model the fixed-point behaviour using a global noise that is added to the output of a floatingpoint simulation.

C. Generalized Gaussian PDF

This noise must have the same characteristics as the global noise due to all quantization sources (mean, power, PDF). So, to generate this noise, we use the Generalized Gaussian PDF [START_REF] Tesei | The asymmetric generalized gaussian function: A new hos-based model for generic noise pdfs[END_REF] defined by

f G (x) = β 2αΓ(1/β) e -| x-µ α | β (2)
where β and α are parameters defining the shape and the scale respectively, and where µ is the mean value of the random variable x.

The term Γ is the Gamma function defined by

Γ(k) = ∞ 0 e -x x k-1 dx (3)
The variance of the Generalized Gaussian random variable can be expressed as a function of parameters α and β, or equivalently parameter α is defined by:

α = Γ(3/β) σ 2 Γ(1/β) (4)
with σ 2 the variance of the variable x. Its 4 th order moment (normalized kurtosis κ σ 4 ) can be computed [START_REF] Tesei | The asymmetric generalized gaussian function: A new hos-based model for generic noise pdfs[END_REF] and is equal to

κ σ 4 = E(x -µ) 4 σ 4 = Γ(5/β)Γ(1/β) Γ(3/β) 2 (5)
From [START_REF] Tesei | Use of fourth-order statistics for nongaussian noise modelling: The generalized gaussian pdf in terms of kurtosis[END_REF] , β can be approximated as

β = 5 κ σ 4 -1.865 -0.12 (6)
for 1.865 < κ σ 4 < 30. For values less than 1.865, the signal x can be modelled with a uniform PDF, whereas for values higher than 30, it can be represented as a dirac at zero. So, with the determination of the mean, the variance and the kurtosis of the output noise, all parameters defining the Generalized Gaussian PDF can be computed. First, β is computed from the knowledge of the kurtosis value and (6). Then, α can be determined from the values of β and σ, and thus, the PDF is completely expressed. So, in the section III-B, we will present a method to compute analytically these parameters.

As we are interested in the generation of samples x(n) following a Generalized Gaussian PDF, we demonstrate in the Appendix that this operation can easily realized by the multiplication of two random variables u(n) and y(n) as :

x = u.y 1/β (7)
where u(n) are uniform data over [-α, α] and y(n) are samples following the PDF Γ(1/β + 1, 1) at the power 1/β.

III. HIGH ORDER MOMENTS

In this section, we present a model to compute the HOM of the output noise. Indeed, to determine the Generalized Gaussian corresponding to the output noise b y (n) the moments 1 -4 must be evaluated. This evaluation must be analytical to prevent a fixed-point simulation. First, HOM of quantization noises are computed. Then, we can use these values to compute the HOM for the output noise.

A. Input noise High Order Moments

In this subsection the HOM of the quantization are determined for the 3 types. We consider a quantization of k eliminated bits and we note ∆ = 2 -n-k the quantization step before the quantization and q = 2 -n is the quantization step after the quantization. The term n is the bit number of the fractional part after quantization.

1) Quantization by truncation: The noise PDF is represented Figure 3 and is equal to

x .2 0 1 2 3 2 p(x )
2 -1 Fig. 3. PDF of a noise generated by quantization by truncation

p b (x) = 1 2 k 2 k -1 j=0 δ(x -j.∆) (8)
So its 3 rd order moment E(x 3 ) is equal to

E(x 3 ) = R p b (x)x 3 dx = ∆ 3 2 k 2 k -1 j=0 j 3 = ∆ 3 2 k 2 2k (2 k -1) 2 4 = ∆ 3 4 (2 3k -2 (2k+1) + 2 k ) = q 3 4 (1 -2.2 -k + 2 -2k ) (9)
The skewness γ is defined by

γ = E(x -µ) 3 = E(x 3 ) -3µσ 2 -µ 3 (10)
Applying values in Table I and ( 9), the skewness is equal to zero.

The 4 th order moment E(x 4 ) is equal to

E(x 4 ) = R p b (x)x 4 dx = ∆ 4 2 k 2 k -1 j=0 j 4 = ∆ 4 2 k ( 2 k j=0 j 4 -2 4k ) = ∆ 4 2 k 2 k (2 k + 1)(6.2 3k + 9.2 2k + 2 k -1) 30 -2 4k = ∆ 4 (2 k + 1)(6.2 3k + 9.2 2k + 2 k -1) 30 -2 3k = ∆ 4 6.2 4k -15.2 3k + 10.2 2k -1 30 = q 4 5 - q 4 .2 -k 2 + q 4 .2 -2k 3 - q 4 .2 -4k 30 = q 4 5 (1 - 5 2 2 -k + 5 3 2 -2k - 1 6 2 -4k ) (11)
The kurtosis κ is defined by

κ = E(x -µ) 4 = E(x 4 ) -4µγ -6µ 2 σ 2 -µ 4 (12)
Applying values in Table I and [START_REF] Widrow | Statistical Theory of Quantization[END_REF], the kurtosis is equal to

κ = q 4 80 (1 - 10 3 2 -2k + 7 3 2 -4k ) (13)
More generally, the kurtosis is normalized by σ 4 to be able to compare the PDF of different signal leading to

κ σ 4 = 1.8 1 -10 3 2 -2k + 7 3 2 -4k 1 -2.2 -2k + 2 -4k (14)
2) Quantization by rounding: In this subsection, we compute the HOM of noises generated by quantization by rounding. Its PDF is presented in Figure ( 4) and is equal to

p b (x) = 1 2 k 2 k-1 -1 j=-2 k-1 δ(x -j.∆) (15) p(x ) x .2 -1 0 1 -2 2 2 -1
E(x 3 ) = R p b (x)x 3 dx = ∆ 3 2 k 2 k-1 -1 j=-2 k-1 j 3 = ∆ 3 2 k 2 (3k-3) = ∆ 3 8 2 2k = q 3 8 2 -k (16) (17)
It skewness is equal to

γ = 0 (18)
The 4 th order moment is equal to

E(x 4 ) = R p b (x)x 4 dx = ∆ 4 2 k 2 k-1 -1 j=-2 k-1 j 4 = ∆ 4 2 k 2 2 k-1 j=0 j 4 -2 4(k-1) (19)
Using the fact that q = ∆.2 k , this expression is developed and leads to

E(x 4 ) = q 4 80 (1 + 20 3 2 -2k - 8 3 2 -4k ) (20)
So, its kurtosis is computed

κ = q 4 80 (1 - 10 3 2 -2k + 7 3 2 -4k ) (21)
It is equal to the kurtosis of a noise generated by truncation.

3) Quantization by convergent rounding: For quantization by convergent rounding, the PDF is given in Figure [START_REF] Kim | Fixed-Point Optimization Utility for C and C++ Based Digital Signal Processing Programs[END_REF] and is equal to

p b (x) = 1 2 K 2 k-1 -1 j=-2 k-1 +1 δ(x -j.∆) + δ(x -2 k-1 .∆) 2 + δ(x + 2 k-1 .∆) 2 (22) p(x )
x .2 -1 0 1 -2 2 2 Fig. 5. PDF of a noise generated by quantization by convergent rounding Its 3 rd order moment is given by the following relation

E(x 3 ) = R p b (x)x 3 dx = ∆ 3 2 k ( 2 k-1 -1 j=-2 k-1 -1 j 3 - 2 k-1 2 + 2 k-1 2 ) = 0 ( 23
)
It skewness is equal to

γ = 0 (24)
The term E(x 4 ) is also developed

E(x 4 ) = R p b (x)x 4 dx = ∆ 4 2 k 2 k-1 -1 j=-2 k-1 -1 j 4 + 2 2 4(k-1) 2 = ∆ 4 2 k 2 2 k-1 j=0 j 4 -2 4(k-1) (25)
Using the same appraoch as before, it leads to

E(x 4 ) = q 4 80 (1 + 20 3 2 -2k - 8 3 2 -4k ) (26)
Its kurtosis κ is equal to

κ = q 4 80 (1 + 20 3 2 -2k - 8 3 2 -4k ) (27)
The normalized kurtosis is given by

κ σ 4 = 1.8 1 + 20 3 2 -2k -8 3 2 -4k 1 + 4.2 -2k + 4.2 -4k (28)
All previous computation of HOM of quantization noises are summarized in the Table II. It can be noted that, for a only one eliminated bit (k = 1) the normalized kurtosis is equal to 1 for quantization by truncation and rounding and equal to 2 for convergent rounding. For a high number of eliminated bits k → ∞, the normalized kurtosis converges to 1.8, which is the classical value for uniform signals with continuous amplitude.

B. Output noise High Order Moments

In this subsection, the computation of the skewness and kurtosis of the output noise b y

(n) is presented. expression of b y (n) (1) is repeated here b y (n) = Ne i=1 n k=0 h i (k)b i (n -k)
(29) [START_REF] Rocher | Analytical approach for numerical accuracy estimation of fixed-point systems based on smooth operations[END_REF] gives mean µ y and variance σ 2 y of output noise b y (n) equal to

µ y = N i=1 µ i n k=0 E[h i (k)] σ 2 y = N 2 i,j=1 µ i µ j n 2 k,m=0 E[h i (k)h j (m)] - N 2 i,j=1 µ i µ j n 2 k,m=0 E[h i (k)]E[h j (m)] + N i=1 σ 2 i n k=0 E[h 2 i (k)] (30)
Using statistical characteristics of quantization noises, the 3 rd order moment E[b 3 y ] of the output noise is equal to

E[b 3 y ] = 3 N 2 i,j µ i σ 2 j n 2 k,m=0 E[h i (k)h 2 j (m)] + i,j,p µ i µ j µ p n 3 k,m,l=0 E[h i (k)h j (m)h p (l)] (31)
Its skewness value γ y is obtained applying ( 10) to ( 31) and (30) leading to (32) page 6

It can be noted that for a LTI system, γ y is equal to zero. So, in that case, the output noise PDF is always symmetric.

With the same approach, we can compute the 4 th order moment E[b 4 y ] of the output noise given by (33) page 6. Given that the skewness of the quantization noises b i (n) is equal to zero, it does not appear in the expression. Its kurtosis value κ y is obtained applying (12) to (33) and (30) leading to (34) page 6

The HOM of the output noise b y (n) are determined analytically. These values allow us to determine completely Generalized Gaussian corresponding to the output noise. In the next section we present some experiments of our model.

IV. RESULTS

In this section, some experiments are carried out to validate our model on FIR and IIR filters and on a 2 nd order Volterra filter.

For the FIR filter, we first experiment a filter with 2 coefficients 1.5 and 0.5 where the inputs are quantized on 8 bits and the output on 8 bits with quantization by truncation. Figure 6 shows that our model is accurate. In that case, the PDF of the output noise is not uniform and not gaussian. So our model can estimate this case. On Figure 7, we apply a FIR filter with 16 coefficients where inputs are quantized on 16 bits and the output on 8 bits. Here, the output noise predominates leading to a uniform PDF as we can see in the figure.

γ y = 3 N 2 i,j µ i σ 2 j n 2 k,m=0 E[h i (k)h 2 j (m)] -E[h i (k)]E[h 2 j (m)] + i,j,p µ i µ j µ p n 3 k,m,l=0 E[h i (k)h j (m)h p (l)] -E[h i (k)]E[h j (m)]E[h p (l)] -3 i,j,p µ i µ j µ p n 3 k,l,m=0 E[h i (k)h j (m)]E[h p (l)] -E[h i (k)]E[h j (m)]E[h p (l)] (32) E[b 4 y ] = i κ i n k=0 E[h 4 i (k)] + 3 N 2 i,j σ 2 i σ 2 j n 2 k,m=0 E[h 2 i (k)h 2 j (m)] + 3 N i σ 4 i n 2 k =m=0 E[h 2 i (k)h 2 j (m)] + 6 N 3 i,j,p µ i µ p 2 j n 3 k,m,l=0 E[h i (k)h 2 j (m)h p (l)] + i,j,p,q µ i µ j µ p µ q n 4 k,m,l,r=0 n l=0 E[h i (k)h j (m)h p (l)h q (r)] (33)
κ y = i κ i n k=0 E[h 4 i (k)] + 3 N 2 i,j σ 2 i σ 2 j n 2 k,m=0 E[h 2 i (k)h 2 j (m)] + 3 N i σ 4 i n 2 k =m=0 E[h 2 i (k)h 2 j (m)] + 6 N 2 i,j,p µ i µ p σ 2 j n 3 k,m,l=0 E[h i (k)h 2 j (m)h p (l)] -E[h i (k)]E[h 2 j (m)]E[h p (l)] - 6
i,j,p,q µ i µ j µ p µ q n 4 k,m,l,r=0

E[h i (k)]E[h j (m)]E[h p (l)h q (r)] -E[h i (k)]E[h j (m)]E[h p (l)]E[h q (r)] -4 i,j,p,q µ i µ j µ p µ q n 4 k,m,l,r=0 E[h i (k)]E[h j (m)h (l)h q (r)] -E[h i (k)]E[h j (m)]E[h p (l)]E[h q (r)] + i,j,p,q µ i µ j µ p µ q n 4 k,m,l,r=0 E[h i (k)h j (m)h p (l)h q (r)] -E[h i (k)]E[h j (m)]E[h p (l)]E[h q (r)] (34)
bits with quantization by truncation. The recursion leads to an accumulation of quantization noise which tends to a gaussian as it is shown on Figure 8. Our model is also accurate in that case.

Fig. 8. Output noise of a IIR2 with input on 8 bits and output on 8 bits Finally, we have applied our model on a 2 nd order Volterra filter whose equation is given by

y(n) = a 11 x 2 (n) + a 22 x 2 (n -1) + a 1 x(n) + a 2 x(n -1) + a 21 x(n)x(n -1)(35)
For this experiment, inputs are quantized on 16 bits and the output on 16 bits with quantization by rounding. Figure 9 illustrates our results. The PDF is not really a gaussian (due to non-linearities in the filter) but our model is still accurate. These experiments demonstrates the validity of our model for different output noise PDF and for different quantization types.

V. CONCLUSION

In this paper, a model to generate the noise on the output of a system due to fixed-point arithmetic is proposed. This model is based on Generalized Gaussian PDF whose parameters are Fig. 9. Output noise of a volterra filter with inputs on 16 bits and output on 16 bits determined by HOM of the output noise. These HOM are analytically determined which allows to evaluate the fixedpoint behaviour of a system without a fixed-point simulation that are really time consuming. Our model has been tested on various applications (FIR filter, IIR filter, Volterra filter) to demonstrate its accuracy. The model will be tried to be extended to asymmetric PDFs.

APPENDIX

In this appendix, the demonstration of the Generalized Gaussian generator is presented. Let X be a random variable with Generalized Gaussian PDF defined as:

f G (x) = β 2αΓ(1/β) e -| x-µ α | β (36)
Theorem 1: The generation of samples x(n) following this PDF can be obtained using the following equality

x = µ + u.y (1/β) (37)
where u(n) are uniform distributed samples on [-α, α], and y(n) are samples following a gamma PDF Γ(1 + 1/β, 1)

Proof: To demonstrate the theorem, we have to use the three following properties

• The PDF f Z (z) of Z = X + c where c is a constant and X a random variable is equal to

f Z (z) = f X (x -c) (38)
• The PDF of a random variable Z = XY , where X and Y are independent is defined by the following relation.

f Z (z) = ∞ -∞ 1 |x| f X (x)f Y (z/x)dx (39)
• The PDF of Z = h(X) where h is an invertible function is equal to

f Z (z) = 1 h ′ oh -1 (z) f X (h -1 (z)) (40)
So if Z = X a where a is a constant Then applying (41), the PDF f V (v) of samples V = Y 1/β is equal to

f V (v) = β 2 Γ(1/β)|v 1-β | ve -v β (43)
Using (39), the PDF of a random variable Z = U.V = U.Y (1/β) , where U is uniformly distributed on [-α, α] is equal to

f Z (z) = ∞ -∞ 1 |x| β 2 Γ(1/β)|x 1-β | xe -x β 1 2α 1 ]|z/α|,∞] dx (44)
A Gamma PDF is defined on R+. We consider the case where z > 0 so f Z (z) = Now, we consider the case where z < 0. The previous expression is developed as follows

f Z (z) = ∞ 0 1 x β 2 Γ(1/β)x 1-β xe -x β 1 2α 1 [-z/α,∞] dx = ∞ -z/α β 2 Γ(1/β) x β-1 e -x β 1 2α dx = β 2 2αΓ(1/β) ∞ -z/α x β-1 e -x β dx = β 2αΓ(1/β) e -(-z/α) β (46)
So, in the general case, adding a mean value µ and using (38), it can be written as

f Z (z) = β 2αΓ(1/β)
e -(|z-µ|/α) β dx (47)

Fig. 1 .

 1 Fig. 1. System noise model for Ne input noises b i (n) and output noise by(n)

Fig. 2 .

 2 Fig. 2. Equivalence between a fixed-point simulation and a floating-point simulation with an additive noise

Fig. 4 .

 4 Fig. 4. PDF of a noise generated by quantization by rounding Its 3 rd order moment is computed as follows

Fig. 6 .

 6 Fig. 6. Output noise of a FIR2 with input on 8 bits

f

  Z (z) = 1 az 1-1/a f X (z 1/a ) (41)Samples Y follow a Gamma PDF with coefficients Γ(1 + 1/β, 1). Their PDF f Y is equal tof Y (y) = 1 (1+1/β) Γ(1/β + 1)y (1/β) e -y = β Γ(1/β) y (1/β) e -y (42)

  x β-1 e -x β dx = β 2αΓ(1/β) e -(z/α) β(45)

TABLE I FIRST

 I

-AND SECOND-ORDER MOMENTS FOR THE THREE CONSIDERED QUANTIZATION MODES.

TABLE II 3

 II rd AND 4 th -ORDER MOMENTS FOR THE THREE CONSIDERED QUANTIZATION MODES.