Rémi Garcia

Anastasia Volkova

Martin Kumm

Alexandre Goldsztejn

Jonas Kühle

Hardware-aware Design of Multiplierless Second-Order IIR Filters with Minimum Adders

Keywords: Digital filters, IIR, optimal design, multiplierless hardware, ILP

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

D IGITAL filters are essential components of modern tech- nology, from medical equipment to scientific instruments. Filter design is a core topic in digital signal processing and control, and efficient filter implementation in software and hardware has received a significant research interest for half a century. Infinite Impulse Response (IIR) filters are a class of widely-used recursive linear time-invariant filters (see Fig. 1). IIR filters can be relatively easily designed in software, but hardware implementation is essential for embedded systems, where performance/power constraints are critical. While classic mobile or medical applications often rely on Application Specific Integrated Circuits (ASICs), other application domains, such as 5G/6G backbones and autonomous vehicles, rely on reconfigurable hardware using Field Programmable Gate Arrays (FPGA).

Classic Fixed-Point (FxP) filter design and implementation flow follows three separate steps:

1) Filter Design (FD) consists in finding real (in practice, double precision) filter coefficients, adhering to the given frequency specification. IIR filters are defined by coefficients of a rational transfer function, for which a stability criteria must be also satisfied. In general, a large R. Garcia, A. Volkova and A. Goldsztejn were with Nantes Université, Centrale Nantes, CNRS, LS2N, 44000 Nantes. Email: firstname.lastname@ls2n.fr M. Kumm and J. Kühle were with Faculty of Applied Computer Science, Fulda University of Applied Sciences, Germany. Email: firstname.lastname@informatik.hs-fulda.de

x n b 0 y n z -1 b 1 -a 1 z -1 b 2 -a 2
Fig. 1: Transposed Direct Form II for a second-order IIR filter.

amount of different filter coefficient sets can realize a given frequency specification; 2) Quantization (Q) converts the coefficients to a FxP format such that they still respect the given frequency response and they still lead to a stable filter; 3) Implementation (I) consists in generating, using quantized coefficients, a valid hardware description. This step exposes a high number of parameters, e. g., the type of multipliers used. The combination of FD and Q steps has been studied extensively since 1960s [START_REF] Smékal | Optimized models of IIR digital filters for fixed-point digital signal processor[END_REF]- [START_REF] Wang | Fixed-point digital IIR filter design using multi-objective optimization evolutionary algorithm[END_REF]. For certain structures of IIR filters it can even be considered solved [START_REF] Gevers | Parametrizations in Control, Estimation and Filtering Problems: Accuracy Aspects[END_REF], but with respect to a signal to quantization noise ratio, which is a probabilistic measure, not guaranteeing numerical safety. The general approach for quantization of a transfer function is nevertheless quite straightforward, passing through iterative increase of coefficient word length.

A large body of work exists for the I step. Hardware filters involve multiplications with constants, for which optimization techniques have been extensively explored. In the multiplierless shift-and-add-based methods, constant multiplications are replaced by additions, subtractions and bit shifts. First, simple digit-based methods, such as Canonical Signed Digit, were used, but later it was shown that graph-based approaches are more general and can achieve better, and even optimal, results [START_REF] Aksoy | Search algorithms for the multiple constant multiplications problem: Exact and approximate[END_REF], [START_REF] Dempster | Constant Integer Multiplication Using Minimum Adders[END_REF]. The associated optimization problem is known as the multiple constant multiplication (MCM) problem, for which heuristics [START_REF] Dempster | Use of Minimum-Adder Multiplier Blocks in FIR Digital Filters[END_REF]- [START_REF] Kumm | Pipelined adder graph optimization for high speed multiple constant multiplication[END_REF] as well as optimal approaches exist [START_REF] Aksoy | Search algorithms for the multiple constant multiplications problem: Exact and approximate[END_REF]- [START_REF]Optimal Constant Multiplication Using Integer Linear Programming[END_REF]. Various cost functions are possible: the highlevel ones count the number of adders required to perform all multiplications, and low-level ones counting the number of full adders. This problem has been successfully modeled as (Mixed) Integer Linear Programming (ILP) and solved using efficient solvers such as CPLEX or Gurobi. Another constant multiplication method, especially relevant for FPGAs, is based on precomputed tables and is called Ken Chapman multiplier (KCM), after its inventor [START_REF] Chapman | Fast Integer Multipliers Fit in FPGAs[END_REF]. It has also been successfully applied to digital filtering [START_REF] Kumm | FIR Filter Optimization for Video Processing on FPGAs[END_REF], [START_REF] Faust | Bit-parallel Multiple Constant Multiplication using Look-Up Tables on FPGA[END_REF], [START_REF] Dinechin | Table-Based versus Shift-And-Add Constant Multipliers for FPGAs[END_REF] and using this method, an approach for optimization of combined Q & I steps has been proposed [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF]. Both KCM and shift-and-add reduce arithmetic resources, by sharing intermediate results, and improve the delay of the circuit. As recursive filters are very problematic to pipeline [START_REF] Parhi | VLSI Digital Signal Processing Systems: Design and Implementation[END_REF], their multiplierless implementation is an attractive alternative to achieve a higher throughput without increasing the implementation cost. However, the resource reduction strongly depends on the coefficient values, which are typically fixed in the previous FD & Q steps. Hence, the obtained implementations are optimized only for one filter instance, or a small sub-set of the overall design space, not permitting overall optimal solution.

The most common model in design of IIR filters is based on the decomposition of a higher-order filter into a cascade of second-order sections implemented with the Direct Form structure [START_REF] Wang | Fixed-point digital IIR filter design using multi-objective optimization evolutionary algorithm[END_REF], [START_REF] Oppenheim | Discrete-Time Signal Processing[END_REF]- [START_REF] Dempster | Comparison of IIR filter structure complexities using multiplier blocks[END_REF]. We focus here on the optimal design of individual second-order filters (see Fig. 1) as a building brick for methods using cascaded forms, e. g., [START_REF] Ansari | Design of multiplierless digital iir filter using modified cuckoo search algorithm[END_REF]. Our ambition is to solve the combined FD & Q & I steps for a second-order IIR within one global optimization. Even for one second-order filter this is a difficult non-linear combinatorial problem. Recent work [START_REF] Ansari | Design of multiplierless digital iir filter using modified cuckoo search algorithm[END_REF] attempts to solve the IIR filter design and cascading problem, while targeting coefficients with the minimal binary writing, but does not lead to optimal solutions due to restricted design-space and a heuristic solver.

In contrast, in this paper we model the overall process as an ILP problem, which is consequently optimally solved using generic solvers. With this work, the whole design space for the FD & Q & I problem is covered: we search for FxP coefficients (i. e., the quantization stage is implicit) that satisfy filter specifications while minimizing the cost of both structural and multiplier blocks in a multiplierless Transposed Direct Form II (tDFII) realization of a second-order IIR filter. Fig. 2 exhibits an example of filter specifications, for which the standard FD & Q design for 5-bit coefficients fails to satisfy the specification after truncation, whilst our approach, directly searching among 5-bit coefficients, finds a solution and it needs only 8 adders.

This paper presents a novel approach, implemented as an automated tool flow, for the simultaneous design of secondorder IIR filters and their hardware implementation. The main contributions of this paper are:

• linearization of frequency constraints upon the rational transfer function. Based on ideas for linearization of FIR filters in [START_REF] Kumm | Design of Optimal Multiplierless FIR Filters with Minimal Number of Adders[END_REF] we formally rewrite the frequency constraints and to deal with inevitable nonlinearities we modify the integer product linearization technique [START_REF] Billionnet | Linear Reformulations of Integer Quadratic Programs[END_REF]; • detection of ranges for all coefficients is indispensable to build an ILP model; we show that bounds on the transfer function's numerator depend on the target frequency specification and the corresponding optimization problem is quadratic convex, which is efficiently solved; • as the IIR filter has two multiplier blocks, we propose a solution that ensures a simultaneous design of two independent multiplier blocks within one MCM model; • we propose a novel search-space reduction method based on detection of symmetries in the solutions. We present and prove a symmetry breaking constraint that leaves the solved with only a quarter of the original design space; • a new operator was developed within FloPoCo1 [START_REF] De Dinechin | Designing Custom Arithmetic Data Paths with FloPoCo[END_REF],

which is a state of the art code generator for hardware arithmetic units. Given the optimal filter coefficients computed with our ILP, FloPoCo now generates faithfullyrounded second-order IIR filters, i. e., the internal data formats are automatically adapted such that the output is accurate up to its last bit. The complete tool flow is available open-source 2 . We validate our approach on real-life and artificial examples by automatically generating and synthesizing VHDL code for FPGA and ASIC targets.

II. ILP FORMULATION FOR SECOND-ORDER IIR FILTER DESIGN A. Problem Definition

The tDFII structure of second order digital filters is represented in Fig. 1. Its output is computed as

y n = 2 k=0 b k u n-k - 2 k=1 a k y n-k , (1)
where a k , b k ∈ R are the filter coefficients. The corresponding transfer function is

H(z) = B(z) A(z) = b 0 + b 1 z -1 + b 2 z -2 1 + a 1 z -1 + a 2 z -2 , (2)
with z ∈ C. Designing a filter means finding the coefficients a k and b k such that the filter is stable and given frequency constraints hold. These constraints can be expressed as bounds on the filter frequency response

β(ω) ≤ H e jω ≤ β(ω) , ∀ω ∈ [0; π] , (3)
where β and β are the given frequency dependent lower and upper bounds. They typically encode constant bounds valid inside frequency intervals, but our definition (3) allows to model them as functions of ω.

There exists a large number of real-valued coefficient sets satisfying the frequency specifications (3) but only a subset of those is representable in FxP arithmetic with a given word length. In our approach, we directly search for a k and b k in a FxP representation and introduce their integer counterparts, a k , b k ∈ Z, that are linked with the real-valued coefficients via

a k = 2 -w+1 2 ga a k , a k ∈ -2 w-1 ; 2 w-1 -1 , (4) b k = 2 -w+1 2 g b b k , b k ∈ -2 w-1 ; 2 w-1 -1 , ( 5
)
where w is the word length (including sign bit), g a and g b are the largest Most Significant Bit (MSB) positions of the FxP representations of a k and b k , respectively. It is typical to use the same word length for a k and b k but we allow different FxP formats in the numerator and denominator of the transfer function. The optimal design process proposed here is cast into the integer programming framework by searching for the integer counterparts a k and b k within the set of all representable numbers defined by w, g b and g a . While the word length w is a user-given parameter in our setting and is typically desired to be as small as possible in order to save resources, the values of the MSB positions g b and g a are not known. Hence, in order to provide an ILP model one must first predict the precision and the range of the FxP coefficients. More precisely, the goal is to predict the MSB positions that, on the one hand, encompass all representable coefficients, i. e., do not lose solutions, and that maximize the precision, i. e., do not overestimate the MSB, on the other hand. We find upper bounds for g a and g b via a pre-processing routine described in Section II-C. The FxP coefficients must satisfy frequency specifications [START_REF] Wang | Fixed-point digital IIR filter design using multi-objective optimization evolutionary algorithm[END_REF], which are non-linear and involve a complex modulus and variable products. Section II-B presents linearized frequency specification constraints, for which we detect a symmetry in the obtained solution. We present and prove a symmetry-breaking constraint in Section II-E, significantly reducing the search-space.

Finally, we search for FxP coefficients realizing the filter with the optimal number of adders in its hardware implementation. As can be inferred from Fig. 1, two MCM blocks, for coefficients a k and b k , need to be optimized. In Section II-D we demonstrate how to combine or improve the state-of-theart ILP formulations [START_REF] Kumm | Multiple Constant Multiplication Optimizations for Field Programmable Gate Arrays[END_REF], [START_REF]Optimal Constant Multiplication Using Integer Linear Programming[END_REF] in order to efficiently model both multiplier blocks in the IIR.

Finally, a short summary of the workflow is presented in Section II-F.

B. Formulating Linear Frequency Constraints

In order to obtain an ILP model, we have to solve two issues related to the frequency constraints [START_REF] Wang | Fixed-point digital IIR filter design using multi-objective optimization evolutionary algorithm[END_REF]. First of all, (3) actually enforces one constraint for each ω ∈ [0; π]. There are infinitely many such constraints to be satisfied. The constraint (3) is called semi-infinite (see [START_REF] Hettich | Semi-infinite programming: Theory, methods, and applications[END_REF] for a survey on this class of constraints). The standard discretization approach to handling semi-infinite constraints consists in discretizing [0; π] into a finite set of frequencies Ω d , leading to a finite number of constraints. This discretization can be dynamically updated, but we consider here a fixed discretization and a a posteriori verification of the semi-infinite frequency constraint [START_REF] Volkova | Reliable Verification of Digital Implemented Filters Against Frequency Specifications[END_REF]. In the rare cases where the verification fails, the verification procedure proposed in [START_REF] Volkova | Reliable Verification of Digital Implemented Filters Against Frequency Specifications[END_REF] outputs a faulty frequency that can be added to the discretized set of frequencies for a new trial. Continuous variables generically lead to active constraints, which require some theory and accurate algorithms with some adaptive discretization to allow discovering active frequencies for optimal coefficients (see [START_REF] Hettich | Semi-infinite programming: Theory, methods, and applications[END_REF], [START_REF] Marendet | A standard branch-and-bound approach for nonlinear semi-infinite problems[END_REF] for details). Integer variables generically do not lead to active constraints. This means that only a finite number of faulty frequencies can happen after the ILP design, before we obtain coefficients verifying rigorously the frequency specifications. This reduction to a finite problem of the semi-infinite constraint is a positive side effect of searching directly FxP representations of the coefficients.

After this discretization of the semi-infinite constraint, a finite number of nonlinear constraints needs to be handled. For a fixed ω ∈ Ω d , the constraint (3) includes a complex absolute value, which involves a square-root of the sum of squared terms, and a fraction. By incorporating (2) into (3), then multiplying with the denominator A e jω , and finally squaring the result we obtain a constraint equivalent to (3)

A e jω 2 β(ω) 2 ≤ B e jω 2 ≤ A e jω 2 β(ω) 2 , (6)
where,

B e jω 2 = 2 k=0 2 l=0 b k b l cos((k -l) ω), (7)
A e jω 2 = 2 k=0 2 l=0 a k a l cos((k -l) ω), with a 0 = 1. (8)
Yet, in [START_REF] Kumm | Pipelined adder graph optimization for high speed multiple constant multiplication[END_REF] and [START_REF] Aksoy | Search algorithms for the multiple constant multiplications problem: Exact and approximate[END_REF], the filter coefficients are still involved in bilinear terms b k b l and a k a l . Billionnet et al. [START_REF] Billionnet | Linear Reformulations of Integer Quadratic Programs[END_REF] proposed a method that allows linearizing products of positive integers assuming that bounds on these numbers are known: consider the product z = xy where x, y, z ∈ N, with x ≤ x and y ≤ y. The linearization [START_REF] Billionnet | Linear Reformulations of Integer Quadratic Programs[END_REF] basically consists in first rewriting one of the positive integers into its binary representation

x = log 2 x +1 i=0 2 i t x,i , (9)
where t x,i are binary auxiliary variables. This constraint ensures that the bits t x,i encode the value of x. Then, the product z = xy becomes a sum of products between the binary variables t x,i and the positive integer y. Finally, such a binary by integer product is common and its well-known linearization involves indicator or big M constraints [START_REF] Glover | Improved Linear Integer Programming Formulations of Nonlinear Integer Problems[END_REF], [START_REF] Oral | A Linearization Procedure for Quadratic and Cubic Mixed-Integer Problems[END_REF].

Here however, x and y correspond to filter coefficients that have no sign restriction. We extend the linearization exposed above to signed integers by adding the auxiliary variables x + , y + ∈ N and x sg , y sg ∈ {0, 1}, and link them by following constraints

x + = |x| , y + = |y| , (10)
x sg = sign(x) , y sg = sign(y) , (11)
z + = x + y + , (12)
where the linearization of the absolute values [START_REF] Kumm | Multiple Constant Multiplication Optimizations for Field Programmable Gate Arrays[END_REF] and the sign constraints [START_REF]Optimal Constant Multiplication Using Integer Linear Programming[END_REF] are well-known and involve indicator or big M constraints [START_REF] Bertsimas | Introduction to linear optimization[END_REF], [START_REF] Mangasarian | Absolute value programming[END_REF], and where z + = x + y + is the positive case we already presented. Finally, z = ±z + and the sign is determined by the values of x sg and y sg directly in the model. This whole linearization relies on the fact that bounds on x and y are known. This is addressed in the next subsection.

C. Stability and Bounds on Filter Coefficients

Necessary and sufficient stability conditions for secondorder filters are well-known [START_REF] Antoniou | Digital Filters: Analysis, Design, and Signal Processing Applications[END_REF]Section 16.8] to be

-2 < a 1 < 2, ( 13
) |a 1 | -1 < a 2 < 1. ( 14
)
As explained before, the absolute value is standardly linearized using indicator or big M constraints, thus ( 13)-( 14) actually fit in an ILP model. From these constraints it is straightforward to derive bounds on a k : a 1 ∈ ]-2; 2[ and a 2 ∈ ]-1; 1[. These bounds are independent of the frequency specification of the filter and yield an upper bound, g a = 1, for the MSB of the coefficients a k . Bounds on b k , however, cannot be obtained independently of the filter specifications. Using the bounds ( 13)-( 14) and the fact the cosine in [START_REF] Aksoy | Search algorithms for the multiple constant multiplications problem: Exact and approximate[END_REF] belongs to [-1; 1], we deduce that A e jω 2 ≤ 16. This bound together with the frequency specification constraints (6) lead to the constraint

B e jω 2 ≤ 16β(ω) 2 , (15)
that needs to be satisfied by the coefficients b k As can be seen from ( 7), B e jω 2 is a quadratic form b T Qb with respect to the variables b k . Its characteristic matrix Q, whose entries are Q kl = cos((k-l)ω), is symmetric and its spectrum is {0, 1-cos(2ω), 2+cos(2ω)}. Its eigenvalues being non-negative, the inequality constraint ( 15) is convex. As a consequence, some lower and upper bounds on the coefficients b k can be computed by solving the convex quadratic problems consisting in minimizing or maximizing b k subject to the convex quadratic constraints [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF] for all frequencies ω ∈ Ω d where Ω d is a discretization of Ω. The global minimum and maximum of these problems are required to be used as valid lower and upper bounds, such global extrema being easily computed by local solvers since the quadratic constraints are all convex and the cost is linear. In particular, common mixed ILP solvers can solve this kind of nonlinear problem. In addition to allowing the linearization of frequency constraints, these bounds permit to decide a first upper bound, g b , for the MSB of the coefficients b k .

These bounds do not fully take into account the specificity of the filter we are designing and use the worst case bounds in [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF]. In any case, it is possible to compute tighter bounds on the coefficients in order to reduce the search space as explained in Section II-E.

D. MCM for Direct-Form IIR Filters

The FxP coefficients satisfying the above frequency and stability constraints are also subject to hardware constraints. We use the ILP-based hardware models for optimal multiplierless constant multiplication presented in [START_REF] Kumm | Multiple Constant Multiplication Optimizations for Field Programmable Gate Arrays[END_REF], [START_REF]Optimal Constant Multiplication Using Integer Linear Programming[END_REF]. These models, given a set of constant coefficients, solve either a minimization w. r. t. number of adders [START_REF] Kumm | Multiple Constant Multiplication Optimizations for Field Programmable Gate Arrays[END_REF] (with bounded adder graph depth) or a satisfiability problem [START_REF]Optimal Constant Multiplication Using Integer Linear Programming[END_REF] (whether an implementation with a given number of adders is feasible).

For the global IIR filter design problem, the coefficients are, however, the unknowns. Similarly to [START_REF] Kumm | Design of Optimal Multiplierless FIR Filters with Minimal Number of Adders[END_REF], we first design a number of linking constraints that bind the coefficient variables for filter design with the inputs of an MCM problem. We refer to them as glue constraints in Fig. 3.

The challenge is that for the IIR filters two multiplier blocks should be optimized simultaneously, one for coefficients a k and one for b k , with the goal of minimization of the total number of adders, i. e., in both multiplier blocks and in the filter structure. For the minimization model [START_REF] Kumm | Multiple Constant Multiplication Optimizations for Field Programmable Gate Arrays[END_REF] this can be achieved by incorporating twice the MCM constraints, one for each multiplier block, and combining the objective functions

min A Ma + A M b - 2 k=0 ζ b k - 2 k=1 ζ a k , (16)
where If a general ILP model [START_REF]Optimal Constant Multiplication Using Integer Linear Programming[END_REF] solving a satisfiability problem is used, a simple duplication of constraints is not possible. When a target number of total adders is given to the satisfiability problem, a trade-off problem arises: how are they distributed between the two multiplier blocks? In order to explore the trade-off space, a simple but costly solution would be to verify all possible distributions via an external loop, e. g., distributing four adders between two multiplier blocks would require solving different models five times.

A
We propose a more efficient solution, where we enhance the MCM model such that the distibution of adders among multiplier blocks is done by the solver, instead. The idea is to produce one adder graph, which can be separated into two independent multiplier blocks, one for a k and one for b k . In other words, we propose an MCM model that, given a certain number of input coefficient variables, adds the constraint that one part of those is separable from the other. Then, the satisfability problem will be solved for both MCM block simultaneously, by the solver. This is achieved in the following way. In the orginal model, a variable c p is associated to every adder p and these variables are linked together with binary variables c p,i,q where c p,i,q = 1 if, and only if, the adder q is the i-th input of adder p. In the new model, we need to ensure that the adders related to the coefficients a k are not computed using the adders related to the coefficients b k , and vice versa.

We propose to add binary variables mcm a p and mcm b p , for each adder p, in order to ensure that the adder can be used for one adder graph or the other, but not both. We add a set of binary variables, s p,q , which encodes the information about two adders p and q being used in the same adder graph, and a set of the following constraints c p,i,q ≤ s p,q ∀p, q, i,

which ensure that an adder cannot be used as the input of another if they are not in the same adder graph. In a similar way, using mcm a p and mcm b p , we ensure that the glue between the filter coefficients and the adders is only possible for adders in the corresponding adder graph.

Both solutions, for the satisfibality or minimzation problems, prioritize coefficient sparsity and even permit an FIR filter design when possible, as a k can be equal to 0.

E. Search Space Reduction

1) Linearized Specifications Projections: The bounds for a k and b k computed in the previous section allow implementing an ILP model by fixing an FxP format and linearizing the bilinear terms involved in the frequency specifications. This first incomplete ILP model does not include the geometry of the adder graph that represents the multiplierless solution as defined in Section II-D. Yet, in order to speed up the solving of the complete ILP model, we tighten the bounds on the coefficients by solving these incomplete ILPs that are simpler than the complete ILP model in the sense that it does not include the geometry of the adder graph defined in Section II-D: they consist in minimizing or maximizing b k subject to the stability constraints and the linearized frequency constraints. These ILPs are obviously more difficult to solve than the continuous convex quadratic problems used in the previous section to obtain crude bounds, but still much easier to solve than the final complete model. As it has been shown in [START_REF] Kumm | Design of Optimal Multiplierless FIR Filters with Minimal Number of Adders[END_REF] in the context of the design of FIR filters, solving these simpler ILPs to obtain tighter bounds for the solving of the final ILP is worthwhile.

2) Symmetry Breaking: The model contains some symmetries that can be broken in order to reduce the search space: the first symmetry consists in simultaneously changing the sign of the values taken by b 0 , b 1 and b 2 , the second symmetry consists in exchanging the values taken by b 0 and b 2 . These two symmetries leave the constraints on the coefficients b k unchanged, as can be easily seen on the explicit expression

b 2 0 + b 2 1 + b 2 2 + 2b 0 b 1 cos(ω) + 2b 1 b 2 cos(ω) + 2b 0 b 2 cos(2ω), (18)
of ( 6) involved in the constraint (15): indeed, ( 18) is insensitive to changing all coefficients sign simultaneously and to exchanging variables b 0 and b 2 . Both symmetries have no incidence on the MCM problem defined in Section II-D, the existence of symmetric adder graphs being obvious for opposed or exchanged coefficients.

As a consequence, from an arbitrary solution with

b 0 = b * 0 , b 1 = b * 1 and b 2 = b * 2 , (19)
we can build three new solutions by simply applying these symmetries to obtain

b 0 = -b * 0 , b 1 = -b * 1 and b 2 = -b * 2 , (20) b 0 = b * 2 , b 1 = b * 1 and b 2 = b * 0 , (21) b 0 = -b * 2 , b 1 = -b * 1 and b 2 = -b * 0 . (22)
The fourth solution ( 22) is obtained by applying the two symmetries consecutively, in any order. Note that some of these four symmetric solutions ( 19)-( 22) may be equal in some special cases, e. g., when b 0 = b 1 = b 2 = 0. Breaking these symmetries means finding additional constraints, called symmetry breaking constraints (SBCs), that remove some symmetric solutions but keep at least one of them [START_REF] Walsh | General Symmetry Breaking Constraints[END_REF]. If necessary, symmetric solutions that have not been calculated due to the SBCs can be built afterward. In the best case, SBCs keep only one solution among all symmetric solutions, in that case the SBCs are called total. Therefore, we expect here to reduce the size of the search space by a ratio of four since solutions come within symmetry classes containing four symmetric solutions. Finding SBCs for general symmetry groups is difficult, e. g., SBC generation for symmetries consisting only of variable permutations rely on modern group theoretic algorithms [START_REF] Goldsztejn | Variable symmetry breaking in numerical constraint problems[END_REF]. In our case, the symmetry group is generated by one variable permutation and one central symmetry. Up to our knowledge, SBCs involving both variable permutations and central symmetries have not yet been investigated and no SBC defined for them.

In order to derive these SBCs, the search space for b 0 , b 1 and b 2 , which is R 3 , is divided into four areas:

Σ 1 = (b 0 , b 1 , b 2 ) ∈ R 3 | b 0 ≥ |b 2 | , ( 23
) Σ 2 = (b 0 , b 1 , b 2 ) ∈ R 3 | -b 0 ≥ |b 2 | , ( 24
) Σ 3 = (b 0 , b 1 , b 2 ) ∈ R 3 | b 2 ≥ |b 0 | , ( 25
) Σ 4 = (b 0 , b 1 , b 2 ) ∈ R 3 | -b 2 ≥ |b 0 | . (26)
Then, one can verify that: a solution lying inside Σ 2 moves to Σ 1 applying the sign symmetry; a solution lying inside Σ 3 moves to Σ 1 applying the exchange symmetry; a solution lying inside Σ 4 moves to Σ 1 applying both symmetries consecutively. As a consequence, one can restrict the search to Σ 1 and reconstruct all solutions using symmetries. This restriction to Σ 1 is achieved by adding the SBC

b 0 ≥ |b 2 | (27)
to the model. Restricting to another Σ i would be lead to another SBC, with an equivalent improvement of the resolution process.

For completeness, it should be noted that on the frontier Σ k ∩ Σ l between areas, two equivalent solutions might still be kept despite the SBC. However, resolving this issue is counterproductive, because the great majority of the search space lies in the interior of the sets Σ i , and removing the symmetries on the boundaries would introduce many additional constraints.

F. Wrapping Up

To wrap it up, our approach is to model both the filter design and the design of a constant multiplication scheme through shift-and-add using a global ILP problem. On top of that, we search for coefficients directly as integers, i. e., in a FxP format with user-given word length. Since filter coefficients, depending on the filter specification and word length, can have different MSB positions, we first perform a pre-processing. This pre-processing consists in solving a quadratic convex optimization problem yielding tight and rigorous bounds on filter coefficients and permitting to define coherent MSB positions and construct the design space for the main model. As Fig. 3 shows, the high-level IIR model consists of the following constraints:

• linearized frequency-specification constraints ( 6)-( 8) ;

• stability constraints ( 13)-( 14);

• (optional) symmetry breaking constraints [START_REF] Oral | A Linearization Procedure for Quadratic and Cubic Mixed-Integer Problems[END_REF] enabling a significant reduction of the design space; • constraints responsible for the design of optimal constantmultiplication blocks for a k and b k using [START_REF] Kumm | Multiple Constant Multiplication Optimizations for Field Programmable Gate Arrays[END_REF] or our modification of [START_REF]Optimal Constant Multiplication Using Integer Linear Programming[END_REF]; • so-called glue constraints, connecting the unknown filter coefficients with multiplier blocks analogously to [START_REF] Kumm | Design of Optimal Multiplierless FIR Filters with Minimal Number of Adders[END_REF]. The above ILP model can now be solved using any available ILP solver. Moreover, despite its non-linearity, the preprocessing problem can be solved using most generic ILP solvers as well, thanks to its convexity.

III. EXPERIMENTAL RESULTS AND DISCUSSION

A. Implemented Toolflow

We implemented our approach into a tool, whose minimalist interface is shown in Fig. 4. The input specification includes the frequency-domain specification of the filter, and also the information for hardware implementation, i. e., the coefficient word length, the input/output formats and performance parameters, such as the required target frequency. All the other parameters (filter coefficients and their FxP format, the multiplierless operator architecture, the representation of intermediate data, etc.) are determined automatically as a part of a global optimization process. One of the goals of our tool is to bring the attention of the filter designer to the higherlevel parameters, e. g., filter specifications, while relying on our optimal implementations. Our new ILP model and the front-end of the proposed tool are implemented in the tool IIRoptim using the julia language, which offers a unified access to the major ILP solvers through the JuMP library. Given the user input, we first construct an ILP model and solve it with one of many open-access or commercial solvers available through JuMP, such as Gurobi, CPLEX, GLPK, etc. The result of the global optimization problem, i. e., the list of filter coefficients and the adder graphs defining the optimal shift-and-add architectures, is then passed on to the FloPoCo tool [START_REF] De Dinechin | Designing Custom Arithmetic Data Paths with FloPoCo[END_REF], which was extended for that. FloPoCo is the state-of-the-art tool for the design and automatic generation of fixed/floating-point arithmetic cores. We implemented a new operator FixIIRShiftAdd within FloPoCo, generating faithfully-rounded multiplierless IIR filters, i. e., only the last output bit might be erroneous and all other bits are guaranteed to be correct. The new operator alleviates the filter designer from all internal architectural decisions and presents a final VHDL code. FloPoCo also requires a target frequency to perform an automatic pipelining [START_REF] Dinechin | Automating the pipeline of arithmetic datapaths[END_REF]. However, this parameter has only effect on the non-recursive parts of the filter as no additional registers can be added within the recursive part, so this parameter has limited effect.

B. Set of Benchmarks and Comparison Approaches

Benchmarks: Although the design of second-order IIR filters is an important part of the design of larger order filters, benchmarks are rarely targeting frequency specifications of individual second-order sections. Hence, we use three sets of filter specifications with increasing filter design difficulty that could without doubt be used in real-life applications. In addition to that, we add another artificial low-pass filter, and a real-life example from [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF].

The normalized low-pass filter specifications are here defined as

1 -δ ≤ H(e iω ) ≤ 1 + δ, ∀ω ∈ [0; ω p ] , (passband) 0 ≤ H(e iω ) ≤ δ, ∀ω ∈ [ω s ; 1] . (stopband)
We fix the initial passband to [0; 0.3], stopband to [0.7; 1] and δ = 0.1. Then, for each of the families of filter specifications, we vary one of the parameters in dependence of a variable k to increase the filter design difficulty. The detailed frequency specifications for each family of filters are given in Table I, and their graphical representation is sketched in Fig. 5. For example, in family lp1, the δ varies from 0.1 to 0.04 with step

0.01k

where k = 0, . . . , 6. However, the designs were possible only up to k = 5, reaching the maximum design possibilities for second-order IIR filters. Analogously, the families lp2 and lp3 increase/reduce the pass/stopband, respectively. The filter specification lp4 is a lowpass with a short stopband.

Finally, our last benchmark hp0 is a real-life application, it is a highpass filter (Fig. 6), which is a compensator used in a magnetic-bearing control system. It was initially derived by discretizing the analog controller [START_REF] Krach | FPGA-based implementation of digital control for a magnetic bearing[END_REF] and further used as an example in word length optimization literature [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF], [START_REF] Sarbishei | Analytical Optimization of Bit-Widths in Fixed-Point LTI Systems[END_REF]. It is most recently used in [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF] to demonstrate a KCM-based faithfully-rounded implementation of IIR filters. This filter is not defined in terms of frequency specifications but by its frequency response (sole poles and zeros are given in [START_REF] Krach | FPGA-based implementation of digital control for a magnetic bearing[END_REF]), hence to implement the three-step approaches we approximate the frequency specification in Matlab using the coefficients from [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF], [START_REF] Sarbishei | Analytical Optimization of Bit-Widths in Fixed-Point LTI Systems[END_REF]. For our implementation, the versatility of an ILP modeling permits to easily integrate frequency response bounds as functions of ω and similarly approximate the filter hp0.

Comparison approaches: We aim at comparing with the classical and state-of-the-art approaches for the IIR design comparing high level metrics like adder counts and implementation results for FPGAs and ASICs. The classical approach passes through three steps 1) FD: in our case double-precision filter coefficients are obtained, if possible, with Matlab's elliptic method; 2) Q: as in Matlab's Fixed-Point design Toolbox, we convert the double coefficients to FxP with the user-specified word length and round-to-nearest mode. Post-design quantization often leads to errors in frequency response or instabilities, hence we increase the word length s. t. the frequency-domain error is below a threshold (< 10 -7 ) and the filter is stable;

3) I: generic multipliers provided by hardware manufacturers are often used. On FPGAs, the digital signal processing (DSP) blocks (whose availability and number on each FPGA vary), are allowed but it might be interesting to disable them to provide comparison with purely LUT-based implementations. We also compare with the state-of-the-art constant multiplications based on KCM [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF] and MCM [START_REF]Optimal Constant Multiplication Using Integer Linear Programming[END_REF] methods.

Different combinations of the above steps are possible. Each of the benchmark filter specifications will be implemented using the following approaches:

• 3-step Generic: Matlab's FD + Q + implementation using generic multipliers (using the VHDL ' * ' operator, potentially using DSP blocks on FPGA); • 3-step Generic NoDSP: same as above with disabled DSP blocks for synthesis and implementation; • 3-step MCM: Matlab's FD + Q + optimal shift-and-add implementation for FxP coefficients; • 2-step KCM: Matlab's FD + approach from [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF], which directly obtains a KCM-based implementation for real coefficients; • Ours: 1-step combined approach, performing directly FD & Q & I using optimal shift-and-add multipliers.

The 3-step Generic (NoDSP) can be seen as a baseline while 3-step MCM and 2-step KCM is the state-of-the-art. With such a setting, our goal is to analyze the benefits of the global approach compared to 2-or 3-step approaches that first fix filter coefficients to some values and only then optimize for implementation. It should be noted that beside the 3-step MCM approach, the existing 2-and 3-step approaches take less than a second for second order IIR filter design, and even the 3-step MCM approach should produce a solution in a few minutes in the worst case.

Bit widths: The coefficient word length is an input parameter for our tool, hence, for each specification we will explore a range of different coefficient word lengths. For the input/output data in hardware implementation, we used 8-, 12and 16-bit configurations but, due to the lack of space, only the results for 16-bit experiments are reported in the paper. The same observations hold for the 8-and 12-bit experiments.

All benchmarks are available and reproducible, see the tool's git repository for full benchmark information.

C. Evaluation of the ILP Model and Design Results

In the following we evaluate the performance of our ILP model, explore the design space and compare with the 3step MCM-based approach to see the benefits w. r. t. number of adders in a shift-and-add implementation. All experiments were performed on a Linux laptop with i7-10810U processor and 32 GB RAM. Any generic solver with an [START_REF] Kumm | Pipelined adder graph optimization for high speed multiple constant multiplication[END_REF] interface for JuMP library [START_REF] Dunning | Jump: A modeling language for mathematical optimization[END_REF] can be used, here we used CPLEX 12.10 [38].

The first remark concerning our tool is that the running times are quite reasonable, varying from 10 seconds for small word lengths (4-5 bits) and going up to a few minutes in the worst case if we push the word length to unnecessarily large values. These solving time results are similar to the ones for the 3-step MCM approach. By default, we use the symmetry breaking constraints as in general we observed a significant improvement (around 2×-20×) in running times, depending on problem complexity. Obviously, the complexity of the ILP model is increasing with increasing the word length, since the ranges of integer variables are doubled with each new coefficient bit, and a few additional variables and constraints are added as well. However, it is not the model complexity but numerical instabilities that represent the main bottleneck in pushing the coefficient word lengths further than 10-11 bits. Indeed, our ILP model for the MCM design makes intensive use of the so-called big-M constraints, that are limited by certain floating-point tolerances internal to the solver, beyond which the solver cannot use efficient floatingpoint arithmetic for integer programming. An alternative to the big-M are the indicator constraints but the drawback is the increased computational time, leading to a similar bound on the maximum coefficient word lengths.

The second remark is that with our tool, an infeasibility of the design problem can be quickly proven. For example, for the specification lp1 4 in just a few seconds we prove that no implementation that perfectly fits the specification with word length 4 is possible. This is an important feature, since when trying to lower the coefficient word length as much as possible, the filter designer can quickly stop the exploration. Inversely, when searching the smallest feasible word length, the design iteration will quickly move on from infeasible ones.

The goal of our tool is to provide optimal architectures w. r. t. the number of adders in the multiplierless implementation. The number of adders is not a fine-grained metric but it enables the design-space exploration a priori, before any hardware synthesis and experiments. It is a good indicator of the performance of implemented systems, as the number of adders is correlated with the number of LUTs.

We claim that optimizing the filter coefficients simultaneously for the filter design and MCM problem yields smaller total number of adders than the 3-step MCM approach, as it covers the whole design space. Table II presents the optimization results for our method, and for the 3-step MCM method. We report results for the smallest word lengths possible, in which our method has a feasible result (and, by construction, no frequency domain error) and the Q step in 3-step MCM yields stable filters with frequency-domain error smaller than 10 -7 .

It can be noticed that our method finds solutions for smaller coefficient word lengths than the 3-step MCM, and with significantly smaller total number of adders for most of the cases (≈46% on average in Table II). In some cases, marked with asterisk, target word lengths for 3-step MCM were too big s. t. the optimal MCM timed-out and the RPAG heuristic [START_REF] Kumm | Pipelined adder graph optimization for high speed multiple constant multiplication[END_REF] was used to obtain adder graphs instead. For the specification lp4 both approaches find designs with 4-bit coefficients, but thanks to efficiently covering the whole design space of all possible FxP coefficients, our method determines coefficients that require fewer adders in the multiplier block and the total of 5 adders, instead of 7. Finally, the benchmark lp3 1 is the only case when the coefficients in 3-step MCM coincide with the ones found by our tool.

Another advantage of the proposed ILP formulation is privileging sparse implementations, which even for second-order filters largely improves performance due to smaller number of structural adders. For instance, the benchmark hp0 was traditionally implemented with all non-zero coefficients [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF], [START_REF] Sarbishei | Analytical Optimization of Bit-Widths in Fixed-Point LTI Systems[END_REF]. We have approximated the frequency response of this compensator using the poles and zeros from the original paper [START_REF] Krach | FPGA-based implementation of digital control for a magnetic bearing[END_REF] and used it as reference in our ILP. As a result, a sparse implementation with 6-bit coefficients was possible, having as coefficients b 0 = 1, b 1 = -1, b 2 = 0, a 1 = -31/32, a 2 = 0 and leading to the total of 3 adders for the whole filter. To compare, the 3-step MCM can provide at minimum 11-bit coefficients implemented with 10 adders. Moreover, our design has poles further apart from the unit circle, which improves the numerical quality of the time-domain implementation.

Our tool provides results for difficult specifications even when FD with Matlab fails. For benchmarks lp1 5 , lp2 3 and lp3 3 no IIR filter design method in Matlab could find doubleprecision coefficients for a frequency response of sufficient quality. Our method, however, successfully completes the task.

It is important to note that our tool provides optimal implementations w. r. t. the total number of adders for a given coefficient word length. In other words, finding the best coefficient word length is still the filter designer task. While typical flow is to stop at the smallest possible coefficient word length, as in Table II, it does not necessarily lead to the optimal implementation w. r. t. all possible word lengths, and increasing the coefficient word length can actually lead to fewer adders.

To illustrate this, Fig. 7 shows synthesis results for the lp1 4 benchmark with coefficient word lengths varying from 4 to 10 bits (the hardware implementation is described in details in Section III-D). First, one can note that our results start at 5-bit coefficients. Indeed, our tool proved that no stable solution satisfying the given frequency specifications is possible for coefficients with word length 4. It can be seen that increasing the word length from 5 bits to 7 bits actually reduces the total number of adders from 8 to 7 and improves the number of LUTs. Moreover, our approach provides a stable behavior: once a 7-adder solution is found, increasing coefficient word length will never lead to a worse implementation. In practice, the ILP either finds different coefficients with the same number of adders, or simply multiplies values by two for each additional bit (which is the case in Fig. 7 starting 7 bits). In other words, even if the user specifies a larger word length than required, our tool finds the best coefficients that might fit in a smaller format and guarantees that trying smaller word length will not give smaller number of adders.

For the classical 3-step approaches, design-space exploration is more difficult and irregular. For these methods, coefficient quantization introduces a frequency-domain error meaning that the quantized filter does not satisfy the frequency specifications any more (see the red line in Fig. 7). This error is highly non-linear, and a typical intuition that increasing coefficient word length improves the quality of filter is simply not true (see the frequency-domain error for 7-and 8-bit coefficients in Fig. 7). Hence, the search for the best coefficient size must be exhaustive for 3-step methods.

D. Hardware Implementation and Discussion

In the following, we first describe in details the faithfullyrounded architectures that we implement in FloPoCo and then discuss the synthesis results obtained for our benchmarks and each approach.

Implemented architecture: The addition in the center of the high-level description of the tDFII filter in Fig. 1 requires the addition of three inputs, hence it must be first separated into two consecutive two-input add operators. The critical path then contains one multiplication and two additions. To cut the critical path, pipeline registers are usually added, leading to three additional registers for tDFII.

In the perspective of extending this work to higher-order IIR filters decomposed into second-order sections, we propose to slightly modify the classic tDFII structure and separate the recursion into its own branch as in Fig. 8. The number of twoinput adders stays the same, and the two additional registers are compensated by the fact that only one pipeline register is required to cut the critical path, instead of three for the tDFII. Hence, at the same hardware cost we get a better modularity of the feedforward/feedback loops.

We aim at providing faithfully-rounded implementations, i. e., the precision of the output y (the Least Significant Bit (LSB) position out ) serves as the accuracy constraint. There are different ways to assign the sizes of data paths for a filter implementation, and it is important to not underestimate the sizes (otherwise overflow occurs) but, on the other hand, assigning larger data sizes will waste resources to compute unnecessary bits. Hence, our goal is to provide a code generator that, given the input/output format and filter coefficients, automatically computes the word sizes of all internal data paths to guarantee the time-domain error smaller than 2 out but not more.

Fig. 8 presents our approach for multiplierless hardware IIR on the example of lp1 4 benchmark, fulfilling frequency specifications in Fig. 2 with 7 adders. Its transfer function was obtained as

H lp14 (z) = 25 × 2 -7 + 40 × 2 -7 z -1 + 25 × 2 -7 z -2 1 -40 × 2 -6 z -1 + 20 × 2 -6 z -2 .
(28) The inputs to the architecture generator are the MSB and LSB positions of the input x and output y, the adder graphs for multiplier blocks a k and b k , and their corresponding LSBs. For instance, here b = -7, a = -6.

Obviously, one cannot compute exactly (or with some fixed precision) on each iteration, truncate to out and simply feed truncated values back into the loop, as this will degrade tremendously numerical quality and accumulated errors will explode. The Worst-Case Peak Gain (WCPG) measure for IIR filters [START_REF] Volkova | Reliable Evaluation of the Worst-Case Peak Gain Matrix in Multiple Precision[END_REF], which has been applied for hardware IIR filters implemented with KCM multipliers [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF], permits to determine the necessary extended internal precision ext s.t. the propagated error never reaches the LSB of the output. For example, in Fig. 8 ext = out + G, where the number of guard bits G for the filter lp1 4 determined with its WCPG is G = 3. Then, the output of the multiplier blocks needs to guarantee its result with accuracy ext . In our architecture we perform all

x (min, in ) < < b < < a ext y (mout, out) out ← 2 5x z -1 z -1 ← 2 ← 3 5y ← 2 ← 3 25x z -1 z -1 ← 2
Fig. 8: The lp1 4 benchmark can be implemented with mere 7 adders for 7-bit coefficients. All additions are exact, the truncation to internal extended format ext is determined s.t. the output is faithfully rounded to out .

additions and shifts exactly, increasing the size of data paths until their truncation to ext . For the generic approach, based on plain VHDL multipliers (using the * operator), we adopt a similar approach.

Synthesis results: FPGA synthesis was performed using Vivado v2020.2 for a Kintex 7 device (xc7k70tfbv484-3). The delay and power results are obtained after place and route. The power was evaluated for a 100 MHz clock for all dynamic parts incl. clocks, logic, signals and DSPs but no static or I/O power to ease the comparison of the cores. ASIC synthesis was performed using Cadence Genus v18.13 targeting a commercial standard-cell library in 14 nm FinFET technology from Global Foundry with 0.8 V supply voltage and regular VTH. Power dissipation is estimated by simulating the final netlist with 1000 random input vectors using Cadence NCSIM to obtain the switching activity. The timing constraints were set to 10 ns to obtain rather compact circuits.

To leave the process of (non-optimal) search for the best coefficient size out of scope, we further compare the designs with the smallest coefficient word lengths, for which our ILP has a feasible solution, and the 3-step approach has a frequency-domain error smaller than 10 -7 . Fig. 9 summarizes the obtained resources, the critical path delay and the power for FPGA and ASIC targets. For each benchmark, on x-axis we see the coefficient word length for our result (left) and for the 3-step quantization (right). For instance, we recognize the values 5 and 9 for the lp1 4 . While for FPGA targets comparison w. r. t. all approaches is reasonable, for ASICs comparison with implementations using DSPs or the KCM table-based method is not representative due to technology differences.

We can observe that on both FPGAs and ASICs our method is always superior to any of the 3-step or KCM-based methods w. r. t. all metrics. On our FPGA benchmarks we observed the average LUT improvement of 48% compared to the best results of other methods, excluding the one using DSPs. In addition to that, the proposed approach offers the lowest delay with an improvement of 27%, which is not as drastic as LUT improvement but expected due to the similar number of delays. Finally, the resource reductions translate to the significant power reduction of 57%, on average, and up to 95% in the best case.

On ASICs benchmarks, our designs yield similar results, reducing the area by 48% and the delay by 27% on average, with the peak reduction by 76% and 51%, respectively. Our new designs demonstrate an impressive average power improvement of 65%, and in the best cases of 90%. This means that our designs divide the power consumption by 4 on average, and by almost 10 in the best observed cases.

Of course, our method is more efficient partly due to the smallest possible coefficient word lengths, but even when using the same word sizes (see benchmark lp4) our design requires less LUTs and has smaller delay.

The KCM-based approach, representing the state-of-the-art for faithfully-rounded IIR filters, had worse general performance than our method in all benchmarks. While the LUT consumption of KCM-based IIR filters is significantly bigger than for our approach, it is comparable to the 3-step Generic NoDSP approach. However, in terms of delay, KCM-based multipliers prove to be generally faster or comparable to 3step methods. It should be noted that the superiority of our approach over KCM-based IIR filters is not a surprise for the small word lengths, according to the recent result [START_REF] Dinechin | Table-Based versus Shift-And-Add Constant Multipliers for FPGAs[END_REF], and while it can be expected that our approach is possible only for word lengths less than 16 bits due to optimization timeout, the KCM-based multipliers will have no issue dealing with large word lengths.

It can be noted that with the increase in frequency specification difficulty, the 3-step methods generally degrade in performance, while our design methods have a more regular behavior, providing small and fast implementations even for complicated filters.

One of the best performance improvements was achieved for the hp0 compensator. This is due to the much smaller coefficient word length than in all the previous literature, the sparsity and the fact that we succeeded in finding a filter with poles further from the unit circle. As a consequence, the number of internal guard bits was smaller too, resulting in only 105 LUTs compared to 286 in the 3-step MCM and 760 in the KCM-based state of the art implementation. Remark that in [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF], hp0 is implemented for 10-bit I/O word lengths but for the uniformity of experiments and without loss of generality, Fig. 9 features the 16-bit I/O implementation, generated with the FloPoCo operator presented in [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF].

IV. CONCLUSION AND PERSPECTIVES

We proposed a new method for the optimal design of multiplierless second-order IIR filters w. r. t. the number of adders. Our approach is based on a combined global optimization problem, which searches for stable filter coefficients directly in FxP format such that the number of adders in a shiftand-add implementation is minimized. Furthermore, we proposed an automated tool which combines our approach with FloPoCo and provides automatic code generation for FPGAs and ASICs. With our tool, the instabilities and coefficient quantization effects on second-order IIR filters become a thing of the past, and efficient faithfully-rounded implementations automatically generated. We proposed a linearized formulation of the combined filter design and MCM problem as one ILP model, which provides a convenient way for extensions. Several design space reduction techniques were proposed, including a novel symmetry breaking constraint, which we formally proved. As a result, the filter design and optimization takes mere seconds.

Our approach has several useful and important consequences. First, with our ILP it is easy and quick to prove the impossibility of a design, with a given coefficient word length and large enough range of MSBs, such that the filter specifications are fully satisfied. Moreover, this would be true for any smaller coefficient word length, providing a filter designer the unprecedented assurance that the design space exploration in that direction can be stopped. Second, if the filter design is possible with a certain given word length, increasing the word length will never yield a larger number of adders, and often the coefficients will remain the same or shifted by one bit. With this property and the fact that all our designs are guaranteed to have zero frequency-domain error, the design-space exploration becomes a more regular process.

For hardware experiments, we provided a faithfully-rounded multiplierless operator for IIR filters within FloPoCo. The synthesis results confirmed that a global optimization approach is superior to the multi-step FD & Q & I classical methods, and even to the state-of-the-art KCM-based IIR filters. With the ILP formulation we search directly the FxP coefficient design space, privileging sparse implementations and sometimes finding the filters that even double-precision Matlab fails to find. After testing the tool on numerous benchmarks, we observed a 48% reduction in number of LUTs, a 27% reduction in delay and a 57% reduction in power, on average, for FPGA targets. Implementing the same benchmarks on 14 nm ASIC confirmed the efficiency of our tool, demonstrating a 48% reduction in area, 27% reduction in delay and 4× smaller power consumption, on average.

The superiority over the KCM-based IIR [START_REF] Volkova | Towards Hardware IIR Filters Computing Just Right: Direct Form I Case Study[END_REF], which introduced faithfully-rounded filters using analysis of the worstcase rounding errors, advances the progress towards reliable IIR filters and demonstrates again that numerical guarantees do not necessarily come at a higher cost. Some efforts are still required to extend our method, in particular to higher order filters. We see two possible directions for that extension: first, a single ILP model which would permit the design of cascaded second order sections and, second, an external loop for the decomposition of specifications into simpler specifications that are reachable by a second order filter. Extension to other structures than the Direct Forms is also a promising direction that should be tackled in the future.

Although the number of adders is a reliable high level metric, optimizing the number full adders instead would be better. We are fairly optimistic on the fact that our method can be refined to minimize that criteria, first as a postdesign optimization of truncations, and then as one global optimization problem. Furthermore, we plan to introduce the truncations of data paths and model the rounding-error in the ILP model. This will permit further performance gains in the implemented filters.

Fig. 2 :

 2 Fig.2: Specification lp1 4 and its reference Matlab, quantized to 5-bit and our 5-bit implementations.

  Ma and A M b denote the number of adders in multiplier blocks for a k and b k , respectively; and ζ a k and ζ b k are binary variables validating whether the respective filter coefficients are zero. In other words, the number of zerovalued coefficients is maximized, privileging sparsity in the implemented filters.

Fig. 4 :

 4 Fig. 4: Interface of the proposed tool.

Fig. 5 :Fig. 6 :

 56 Fig. 5: Proposed families of benchmarks.

8 Fig. 7 :

 87 Fig. 7: Implementations for lp1 4 benchmark with coefficient word lengths varying from 4 to 10 bits. The bars correspond to the number of LUTs (left axis) and labels on bars denote number of adders. The frequency-domain error of 3-step methods is the red line (right axis).

lp1 1 lp1 2 lp1 3 lp1 4 lp1 5 lp2 0 lp2 1 lp2 2 lp2 3 lp3 0 lp3 1 lp3 2 lp3 3 Fig. 9 :

 39 Fig.9: Benchmark results of resources, power and critical path delay using I/O word length 16 bits for FPGA (top) and ASIC (bottom) targets. The x-axis features the benchmark name and the smallest possible coefficient word length for our designs (left) and 3-step approaches (right).

TABLE I :

 I Sets of lowpass filters used for the IIR experiments. First the set with decreasing δ, next the sets with increasing passband and stopband. Finally, a lowpass filter.

	Benchmarks	lp1 k	lp2 k	lp3 k	lp4
	k	{0, 1, . . . , 6}	{0, 1, . . . , 4}	{0, 1, . . . , 4}	-
	passband/π	[0; 0.3]	[0; 0.3 + 0.05k]	[0; 0.3]	[0; 0.5]
	stopband/π	[0.7; 1]	[0.7; 1]	[0.7 -0.05k; 1]	0.91
	δ	0.1 -0.01k	0.1	0.1	0.1

TABLE II :

 II Results for our global optimization method vs. applying optimal MCM upon quantized coefficients. The total number of adders A (A = A M +A S ) consists in the multiplier block A M , and structural A S adders. Results are reported for the smallest coefficient word length W possible.

	Benchmark		Our method			3-step MCM
		W AM AS A W AM AS A
	lp10	4	1	4	5	10 8	4	12
	lp11	4	1	4	5	16 11	4	15 *
	lp12	4	2	4	6	6	5	4	9
	lp13	4	3	4	7	10 7	4	11
	lp14	5	4	4	8	9	7	4	11
	lp15	5	4	4	8	-	-	-	-
	lp20	4	1	4	5	10 8	4	12
	lp21	4	1	4	5	10 8	4	12
	lp22	5	3	4	7	10 8	4	12
	lp23	6	4	4	8	-	-	-	-
	lp30	4	1	4	5	10 8	4	12
	lp31	4	2	4	6	4	2	4	6
	lp32	4	3	4	7	23 18	4	22 *
	lp33	5	3	4	7	-	-	-	-
	lp4	4	1	4	5	4	3	4	7
	hp0	6	1	2	3	11 6	4	10

* heuristic solution using

http://flopoco.org/

https://gitlab.univ-nantes.fr/volkova-a/jiir2hw