Design of Optimal Multiplierless FIR Filters with Minimal Number of Adders

Martin Kumm, Member, IEEE, Anastasia Volkova, and Silviu-Ioan Filip, Member, IEEE Abstract-This work presents two novel methods that simultaneously optimize both the design of a finite impulse response (FIR) filter and its multiplierless hardware implementation. We use integer linear programming (ILP) to minimize the number of adders used to implement a direct/transposed FIR filter adhering to a given frequency specification. The proposed algorithms work by either fixing the number of adders used to implement the products (multiplier block adders) or by bounding the adder depth (AD) used for these products. The latter can be used to design filters with minimal AD for low power applications. In contrast to previous multiplierless FIR filter approaches, the methods introduced here ensure adder count optimality. We perform extensive numerical experiments which demonstrate that our simultaneous filter design approach yields results which are in many cases on par or better than those in the literature.

Index Terms-FIR filters, multiplierless implementation, ILP optimization, MCM problem

I. INTRODUCTION

F Inite impulse response (FIR) filters are fundamental building blocks in digital signal processing (DSP). They provide strict stability and phase linearity, enabling many applications. However, their flexibility typically comes at the expense of a large number of multiplications, making them compute-intensive. Hence, many attempts have been made in the last four decades to avoid costly multiplications and to implement FIR filters in a multiplierless way [START_REF] Lim | FIR filter design over a discrete powers-of-two coefficient space[END_REF]- [START_REF]Exact and Approximate Algorithms for the Filter Design Optimization Problem[END_REF].

One of the most widespread ways to do so is to replace constant multiplications by additions, subtractions and bit shifts. Consider, for instance, the multiplication by a constant coefficient 23. It can be computed without dedicated multipliers as 23x = 8 • (2x + x)x = ((x << 1) + x) << 3)x, [START_REF] Lim | FIR filter design over a discrete powers-of-two coefficient space[END_REF] where (x << b) denotes the arithmetic left shift of x by b bits. This computation uses one addition and one subtraction. Since in this context the add and subtract operations both have similar hardware cost, the total number of add/subtract units is usually referred to as adder cost. Bit shifts can be hard-wired in hardware implementations and do not contribute any cost. In general, the task of finding a minimal adder circuit for a given constant is known as the single constant multiplication (SCM) problem and is already an NP-complete optimization Manuscript received xx xx, 2021; revised xx xx, 2021. M. Kumm is with the Fulda University of Applied Sciences, 36037 Fulda, Germany (e-mail: martin.kumm@cs.hs-fulda.de).

A. Volkova is with Nantes Université, CNRS, LS2N, 44000 Nantes, France. (e-mail: anastasia.volkova@univ-nantes.fr).

S.I. Filip is with University of Rennes, Inria, CNRS, IRISA, Rennes, France (e-mail: silviu.filip@inria.fr). problem [START_REF] Cappello | Some Complexity Issues in Digital Signal Processing[END_REF]. For [START_REF] Lim | FIR filter design over a discrete powers-of-two coefficient space[END_REF], the corresponding SCM adder circuit is illustrated in Fig. 1a. Such a problem extends to multiplication with multiple constants, which is necessary when implementing FIR filters. It is called multiple constant multiplication (MCM). Here, some of the intermediate factors like the adder computing 3x in Fig. 1a can be shared among different outputs. Take for example the coefficients {7, 23}; Fig. 1b shows a solution for multiplying with both coefficients at an adder cost of only two. The corresponding optimization problem is called the MCM problem and has been addressed by numerous heuristic [START_REF] Dempster | Use of Minimum-Adder Multiplier Blocks in FIR Digital Filters[END_REF]- [START_REF] Gustafsson | A Difference Based Adder Graph Heuristic for Multiple Constant Multiplication Problems[END_REF] and optimal [START_REF]Towards Optimal Multiple Constant Multiplication: A Hypergraph Approach[END_REF]- [START_REF] Kumm | Optimal Constant Multiplication using Integer Linear Programming[END_REF] approaches.

Fig. 2 shows the two most popular structures used to implement FIR filters: the direct and transposed forms. The result of an MCM solution can be directly placed in the multiplier block of the transposed form from Fig. 2b. The total adder cost can be modeled as the sum of the number of multiplier block adders and the remaining ones, commonly called structural adders. The transposed form can be obtained from the direct form by transposition [START_REF] Crochiere | Analysis of Linear Digital Networks[END_REF]. As the transposition of a single-input-single-output system does not change the adder count, it leads to the same adder cost. In the end, it does not matter which one of these two filter structures is actually optimized with respect to the adder count.

In the MCM optimization problem, it is assumed that the coefficients are known and already quantized to a fixed-point (or integer) representation. The design of FIR filters with fixed-point coefficients and a minimum frequency response approximation error is itself a well-known optimization problem, going back to at least [START_REF] Kodek | Design of optimal finite wordlength FIR digital filters using integer programming techniques[END_REF] with subsequent extensions and improvements [START_REF] Lim | FIR filter design over a discrete powers-of-two coefficient space[END_REF]- [START_REF] Redmill | Genetic synthesis of reduced complexity filters and filter banks using primitive operator directed graphs[END_REF], [START_REF] Gustafsson | Design of linear-phase FIR filters combining subexpression sharing with MILP[END_REF], [START_REF] Aksoy | Optimization of Area Under a Delay Constraint in Digital Filter Synthesis Using SAT-Based Integer Linear Programming[END_REF]- [START_REF] Xu | Design of Low-Complexity FIR Filters Based on Signed-Powers-of-Two Coefficients With Reusable Common Subexpressions[END_REF], [START_REF] Aksoy | Optimization Algorithms for the Multiple Constant Multiplications Problem[END_REF], [START_REF] Aksoy | Design of Low-Power Multiple Constant Multiplications Using Low-Complexity Minimum Depth Operations[END_REF], [START_REF] Shahein | A Novel Hybrid Monotonic Local Search Algorithm for FIR Filter Coefficients Optimization[END_REF], [START_REF] Ye | Design of Low Power Multiplierless Linear-Phase FIR Filters[END_REF], [START_REF] Aksoy | A Tutorial on Multiplierless Design of FIR Filters: Algorithms and Architectures[END_REF], [START_REF]Exact and Approximate Algorithms for the Filter Design Optimization Problem[END_REF], [START_REF] Lim | Design of discrete-coefficient-value linear phase FIR filters with optimum normalized peak ripple magnitude[END_REF], [START_REF] Kodek | Design of optimal finite wordlength FIR digital filters[END_REF]. It is however often the case in practice that a bounded frequency response is acceptable. In fact, there may be a large number (often hundreds or more) of different fixed-point coefficient sets that meet such a specification. Because of this, the most widely used filter design technique is a 2-step approach in which: (a) a filter with real-valued coefficients is derived using standard approaches (e.g. windowing or Chebyshev) after which (b) the obtained values are quantized and optimized for a minimum number of adders using an MCM approach.

Still, when optimizing for resources (in this particular case, total number of adders), the obtained results are usually far from optimal (see Section VI-C for a comparison). Therefore, a lot of effort in fixed-point FIR filter design has gone into aggregate 1-step methods that put emphasis on resource use. Finding a minimal adder circuit for a given filter specification has thus been approached by several authors [START_REF] Lim | FIR filter design over a discrete powers-of-two coefficient space[END_REF]- [START_REF] Redmill | Genetic synthesis of reduced complexity filters and filter banks using primitive operator directed graphs[END_REF], [START_REF] Gustafsson | Design of linear-phase FIR filters combining subexpression sharing with MILP[END_REF], [START_REF] Yu | Design of Linear Phase FIR Filters in Subexpression Space Using Mixed Integer Linear Programming[END_REF], [START_REF] Xu | Design of Low-Complexity FIR Filters Based on Signed-Powers-of-Two Coefficients With Reusable Common Subexpressions[END_REF], [START_REF] Aksoy | Design of Low-Power Multiple Constant Multiplications Using Low-Complexity Minimum Depth Operations[END_REF], [START_REF] Shahein | A Novel Hybrid Monotonic Local Search Algorithm for FIR Filter Coefficients Optimization[END_REF], [START_REF] Ye | Design of Low Power Multiplierless Linear-Phase FIR Filters[END_REF], [START_REF]Exact and Approximate Algorithms for the Filter Design Optimization Problem[END_REF], [START_REF] Aksoy | SIREN: a depth-first search algorithm for the filter design optimization problem[END_REF].

However, to the best of our knowledge, only one approach, SIREN [START_REF]Exact and Approximate Algorithms for the Filter Design Optimization Problem[END_REF], [START_REF] Aksoy | SIREN: a depth-first search algorithm for the filter design optimization problem[END_REF] , has actually addressed this 1-step design of multiplierless FIR filters in an optimal way, using a custom branch and bound algorithm. Here, by optimal multiplierless filter we mean a direct/transposed form FIR filter requiring a minimum number of addition/subtraction operations to meet a target frequency specification, as well as constraints on the maximum coefficient word size and filter order.

This work presents for the first time a closed form integer linear programming (ILP) formulation for solving this problem. The authors believe that having a mathematical ILP formulation allows for easy re-implementations and offers a convenient framework for further extensions. Examples are minimizing the power consumption [START_REF] Demirsoy | Transition Analysis on FPGA for Multiplier-Block Based FIR Filter Structures[END_REF]- [START_REF] Faust | Minimal Logic Depth Adder Tree Optimization for Multiple Constant Multiplication[END_REF], the inclusion of lookup table-based multipliers [START_REF] Kumm | FIR Filter Optimization for Video Processing on FPGAs[END_REF], targeting 3-input adders for FPGAs [START_REF] Simkins | Structures and Methods for Implementing Ternary Adders/Subtractors in Programmable Logic Devices[END_REF], [START_REF] Baeckler | Logic Cell Supporting Addition of Three Binary Words[END_REF], pipelining [START_REF] Kumm | FIR Filter Optimization for Video Processing on FPGAs[END_REF], optimizing FIR cascaded filters [START_REF] Mehrnia | Optimal Factoring of FIR Filters[END_REF], etc. Besides that, any performance improvements to ILP solvers directly translate into faster run-times for solving the multiplierless FIR filter design problem.

The main contributions of our work can be summarized as:

• We present for the first time a solution for the optimal multiplierless design (in terms of adder count) of FIR filters from a frequency specification using a closed form ILP formulation. • We provide another ILP formulation that is capable of additionally limiting the adder depth inside the FIR filter. • We show that relevant problem sizes can be addressed by current ILP solvers and that the adder complexity of well-known FIR filters can be further reduced compared to the most advanced methods.

Our approach builds on previous ILP formulations [START_REF] Kumm | Optimal Constant Multiplication using Integer Linear Programming[END_REF], [START_REF] Kumm | FIR Filter Optimization for Video Processing on FPGAs[END_REF] for optimally solving the MCM problem, which we extend here to multiplierless FIR design. The major challenge is that such previous work assumes that the constants are known in advance, while in the FIR filter setting they are the unknowns. To address it we formalized the search space for the coefficients, incorporated the frequency domain constraints into the problem and "linked" the unknowns with the MCM problem inputs. We also formalized the notion of structural adders in order to efficiently exploit sparse filter implementations. In addition, the proposed tool takes advantage of a variety of techniques for reducing the search space of the obtained models.

In the following, we will give background information about previous work this paper is based on. In Section III and Section IV we describe the two ILP formulations that are at the core of the paper, whereas in Section V we talk about ideas meant to improve the practical run-time of the proposed algorithms. We then present experimental results accompanied by a comparison with the state-of-the-art (Section VI), followed by concluding remarks (Section VII).

II. BACKGROUND

Multiplierless filter design problems usually start with a functional specification of the frequency domain behavior, together with the number of filter coefficients and their word lengths. An optimization procedure is applied to get a set of bounded integer coefficients together with their associated adder circuits for the constant multiplications needed in the final implementation. Summarized in Fig. 3, this section overviews these parameters and their interactions, together with the state-of-the-art design methods found in the literature.

A. Linear Phase FIR Filters

An N -th order linear phase FIR filter can be described by its zero-phase frequency response

H R (ω) = M -1 m=0 h m c m (ω), ω ∈ [0, π], (2)
I sym. even N 2 + 1 cm(ω) = 1 for m = 0 2 cos(ωm) for m > 0 II sym. odd N +1 2 cm(ω) = 2 cos(ω(m + 1/2)) III asym. even N 2 cm(ω) = 2 sin(ω(m -1)) IV asym. odd N +1 2 cm(ω) = 2 sin(ω(m + 1/2))
which has the property that its magnitude is identical to that of the transfer function, i.e.,

|H(e jω )| = |H R (ω)| . (3)
The c m (ω) terms are trigonometric functions and M denotes the number of independent coefficients after removing identical or negated ones due to symmetry. Both depend on the filter symmetry and on the parity of N as given in Table I. Let D(ω) and D(ω) be the desired lower and upper bounds of the output frequency response H R (ω). The associated frequency specification-based FIR filter design problem consists of finding coefficients h m , m = 0, . . . , M -1 that fulfill the constraints

D(ω) ⩽ H R (ω) ⩽ D(ω), ∀ω ∈ Ω, (4)
where Ω ⊆ [0, π] is a set of target frequency bands (usually pass and stopbands). A standard approach in practice is to work with Ω d ⊆ Ω, a uniform discretization of Ω. One number for the size of Ω d found in the literature is 16M [START_REF] Kodek | Design of optimal finite wordlength FIR digital filters[END_REF].

B. Fixed-point Constraints

Fixed-point (integer coefficient) FIR filter design problems further restrict the search space to integer variables h

′ m ∈ Z with |h ′ m | < 2 B , where the coefficients of H R (ω) are h m = 2 -B h ′ m (5)
and B ∈ N is the maximum effective word length of each coefficient (excluding sign bit). Note that we do not rely on any number representations or other limited number spaces like many previous approaches [START_REF] Lim | FIR filter design over a discrete powers-of-two coefficient space[END_REF], [START_REF] Samueli | An improved search algorithm for the design of multiplierless FIR filters with powers-of-two coefficients[END_REF], [START_REF] Yli-Kaakinen | A systematic algorithm for the design of multiplierless FIR filters[END_REF]- [START_REF] Vinod | On the implementation of efficient channel filters for wideband receivers by optimizing common subexpression elimination methods[END_REF], [START_REF] Yu | Design of Linear Phase FIR Filters in Subexpression Space Using Mixed Integer Linear Programming[END_REF]- [START_REF] Aktan | An algorithm for the design of low-power hardware-efficient FIR filters[END_REF], [START_REF] Yu | Optimization of Linear Phase FIR Filters in Dynamically Expanding Subexpression Space[END_REF], [START_REF] Shi | Design of Linear Phase FIR Filters With High Probability of Achieving Minimum Number of Adders[END_REF], [START_REF] Shahein | A Novel Hybrid Monotonic Local Search Algorithm for FIR Filter Coefficients Optimization[END_REF].

To broaden the feasible set of efficient designs, some applications allow the use of a real-valued scaling factor G > 0 when computing the quantized fractional coefficients h m . Equation (4) thus becomes

GD(ω) ⩽ H R (ω) ⩽ GD(ω), ∀ω ∈ Ω . (6)
When the frequency specification contains a passband, it is called the passband gain [START_REF] Lim | Design of discrete-coefficient-value linear phase FIR filters with optimum normalized peak ripple magnitude[END_REF]. Finding adequate bounds for G is dependent on the set/format of feasible h ′ m coefficient values. If these values are constrained to a power of two space, the ratio between the upper and the lower bound on G does not need to be larger than 2 [33, Lemma 1]. Even when this is not the case, the interval [0.7, 1.4] is frequently used [START_REF] Shi | Design of Linear Phase FIR Filters With High Probability of Achieving Minimum Number of Adders[END_REF], [START_REF] Ye | Design of Low Power Multiplierless Linear-Phase FIR Filters[END_REF], [START_REF] Lim | Design of discrete-coefficient-value linear phase FIR filters with optimum normalized peak ripple magnitude[END_REF]. For our tests, unless otherwise stated, we project the polytope described by [START_REF] Gustafsson | Design of linear-phase FIR filters combining subexpression sharing with MILP[END_REF] onto G in order to obtain a sufficiently large search domain [G, G] (see Section V-A). In case a unity or fixed-gain filter is required we use G = 1.

C. Multiplierless FIR Filters

Formulas ( 5) and ( 6) are easily expressed as constraints in an ILP formulation. However, to ensure an optimal multiplierless design, further constraints are needed.

The way these constraints are constructed and used has varied over the years. Early research in this direction looked at multiplierless designs where each coefficient was represented by a limited number of signed power-of-two terms, optimized using branch-and-bound techniques [START_REF] Lim | FIR filter design over a discrete powers-of-two coefficient space[END_REF]. Later, minimum signed digit (MSD) representations like the canonic signed digit (CSD) representation characterized by a minimum number of non-zero digits were quickly adopted for this purpose [START_REF] Samueli | An improved search algorithm for the design of multiplierless FIR filters with powers-of-two coefficients[END_REF], [START_REF] Yli-Kaakinen | A systematic algorithm for the design of multiplierless FIR filters[END_REF], [START_REF] Aktan | An algorithm for the design of low-power hardware-efficient FIR filters[END_REF], [START_REF] Shahein | A Novel Hybrid Monotonic Local Search Algorithm for FIR Filter Coefficients Optimization[END_REF].

The MSD representations of a number include all possible representations that have a minimum number of non-zero digits. This includes the CSD representation, but also alternatives. It can be used to find sharing opportunities of intermediate computations like the 7x term shown in Fig. 1b. One way is by searching and eliminating redundant bit patterns common to several coefficients, a technique called common subexpression elimination (CSE). Savings are obtained by performing the computation specified by the bit pattern and distributing the result to all coefficients depending on it [START_REF] Gustafsson | Design of linear-phase FIR filters combining subexpression sharing with MILP[END_REF], [START_REF] Vinod | On the implementation of efficient channel filters for wideband receivers by optimizing common subexpression elimination methods[END_REF], [START_REF] Dempster | Use of Minimum-Adder Multiplier Blocks in FIR Digital Filters[END_REF]. However, the CSE search cannot deliver all possible sharing opportunities due to its dependency on the number representation [START_REF] Voronenko | Multiplierless Multiple Constant Multiplication[END_REF] and the effect of hidden non-zeros [START_REF] Faust | Minimal Logic Depth Adder Tree Optimization for Multiple Constant Multiplication[END_REF]. To avoid them, graph-based approaches are commonly used in state-of-the-art MCM methods [START_REF] Dempster | Use of Minimum-Adder Multiplier Blocks in FIR Digital Filters[END_REF]- [START_REF] Kumm | Optimal Constant Multiplication using Integer Linear Programming[END_REF]. Some early work on multiplierless FIR filter design already considered this by incorporating the graph-based MCM algorithm of [START_REF] Dempster | Use of Minimum-Adder Multiplier Blocks in FIR Digital Filters[END_REF] into a genetic algorithm that optimizes the filter coefficients according to the adder cost [START_REF] Redmill | Genetic synthesis of reduced complexity filters and filter banks using primitive operator directed graphs[END_REF]. A different approach is followed by [START_REF] Yu | Design of Linear Phase FIR Filters in Subexpression Space Using Mixed Integer Linear Programming[END_REF], where a branch-and-bound-based ILP optimization is used; here, a pre-specified set of integer terms that can be shared among the different coefficient expansions, called the subexpression space, has to be provided. This work was later extended with a dynamic subexpression space expansion algorithm [START_REF] Yu | Optimization of Linear Phase FIR Filters in Dynamically Expanding Subexpression Space[END_REF], [START_REF] Shi | Design of Linear Phase FIR Filters With High Probability of Achieving Minimum Number of Adders[END_REF], which, at least in the case of [START_REF] Shi | Design of Linear Phase FIR Filters With High Probability of Achieving Minimum Number of Adders[END_REF], claims to usually produce designs with a minimal number of adders. In contrast to these potentially slow branch-and-bound approaches, in [START_REF] Bin Ye | A polynomial-time algorithm for the design of multiplierless linear-phase FIR filters with low hardware cost[END_REF], a fast polynomial-time heuristic for the design of low complexity multiplierless linear-phase FIR filters was proposed.

Although many of the approaches described above use optimal branch-and-bound or ILP methods, they are only used on a limited search space corresponding to the selected number representation. The branch-and-bound method SIREN, first described in [START_REF] Aksoy | SIREN: a depth-first search algorithm for the filter design optimization problem[END_REF] and later refined in [START_REF]Exact and Approximate Algorithms for the Filter Design Optimization Problem[END_REF], is the only other work the authors are aware of, besides this one, which addresses the optimal multiplierless FIR filter design problem regardless of the number representation. It performs a depthfirst search on a search tree where each level consists of the possible values of one of the coefficients. When reaching the bottom of the tree, their optimal MCM algorithm [START_REF] Aksoy | Search Algorithms for the Multiple Constant Multiplications Problem: Exact and Approximate[END_REF] is used to determine the adder cost. Lower and upper bounds of these coefficients are computed using linear programming (LP) and clever ways are proposed to further prune the search tree. The first objective in [START_REF]Exact and Approximate Algorithms for the Filter Design Optimization Problem[END_REF] is to find the minimal effective word length B, with the number of adders being a secondary objective. The method could work in principle with any B, but run-time will become prohibitive, as the search tree increases exponentially with B.

Recent work has also focused on integrating filter coefficient sparsity, which can also have a big impact on the complexity of the final design [START_REF] Ye | Bit-level multiplierless FIR filter optimization incorporating sparse filter technique[END_REF], [START_REF] Chen | Design of sparse fir filters with reduced effective length[END_REF], [START_REF] Chen | Cascaded Form Sparse FIR Filter Design[END_REF] by reducing the number of structural adders. Also, other structures than the direct and transposed forms (see Fig. 2) have been shown to possess good properties. The factoring of FIR filters into a cascade of relatively small subsections can lead to a lower bit-level complexity [START_REF] Mehrnia | Optimal Factoring of FIR Filters[END_REF]. Alternative structures have also been proposed [START_REF] Lou | Novel Structure for Area-Efficient Implementation of FIR Filters[END_REF], [START_REF] Chen | Cascaded Form Sparse FIR Filter Design[END_REF]- [START_REF] Wang | Matrix decomposition based lowcomplexity fir filter: Further results[END_REF]; they provide lower word sizes for the structural adders, reducing resource use.

Besides optimizing the adder count, it was shown early that the power consumption of the resulting filter also strongly depends on the adder depth (AD), which is defined as the number of cascaded adders in the multiplier block [START_REF] Demirsoy | Transition Analysis on FPGA for Multiplier-Block Based FIR Filter Structures[END_REF]. Since then, many works have focused on limiting the AD either in MCM algorithms [START_REF] Dempster | Designing Multiplier Blocks with Low Logic Depth[END_REF]- [START_REF] Kumm | FIR Filter Optimization for Video Processing on FPGAs[END_REF] or directly in multiplierless filter design methods [START_REF] Shi | Design of Linear Phase FIR Filters With High Probability of Achieving Minimum Number of Adders[END_REF]. Again, all of those approaches are heuristics that provide minimal AD, but do not guarantee minimal adder cost.

A low-level metric, e.g. minimizing the number of full adders [START_REF] Johansson | A Detailed Complexity Model for Multiple Constant Multiplication and an Algorithm to Minimize the Complexity[END_REF] or gates [START_REF] Aksoy | Optimization of Area in Digital FIR Filters using Gate-Level Metrics[END_REF], would lead to more hardwareefficient results. However, modeling with respect to these metrics significantly increases the size and complexity of the optimization problem, limiting its practicality. While counting adders is a larger-grain approach, it nevertheless gives a good indication of the hardware resources needed in practice, and is better-suited for efficient ILP modeling.

III. MULTIPLIERLESS FIR FILTERS WITH FIXED NUMBER

OF MULTIPLIER BLOCK ADDERS Our first ILP model targets the design of generic multiplierless FIR filters regardless of their adder depth. It is based on a recently proposed MCM ILP formulation [START_REF] Kumm | Optimal Constant Multiplication using Integer Linear Programming[END_REF], where the goal is to directly compute the parameters of an MCM adder graph, if feasible, for a given number of adders. This idea is extended here for multiplierless FIR filter design by adding constraints on the frequency specification. As a result, we get an ILP model to design a multiplierless filter for a fixed number of adders in the multiplier block. To optimize the total number of adders, this ILP model is applied several times using an overall algorithm discussed in Section III-B. In the following, we first present the ILP formulation for the fixed number of multiplier block adders.

A. ILP Formulation for Fixed Multiplier Block Adder Count

The proposed ILP formulation is given in ILP Formulation 1 and uses the constants and variables listed in Table II. The objective is, given a fixed number of multiplier block adders A M , to minimize the number of structural adders A S (which

ILP Formulation 1 Multiplierless FIR filters with fixed A M minimize A S (h m,0 ) subject to C1a: G2 B D(ω) ⩽ M -1 m=0 h ′ m c m (ω) ⩽ G2 B D(ω), ∀ω ∈ Ω d C1b: h m ⩽ h ′ m ⩽ h m , ∀m = 0, . . . , M - C2: G ⩽ G ⩽ G C3a: h ′ m = (-1) ϕ 2 s c a if o a,m,s,ϕ = ∀a = 0, . . . , A M , m = 0, . . . , M - C3b: h ′ m = 0 if h m,0 = 1, ∀m = 0, . . . , M - C3c: AM a=0 Smax s=Smin 1 ϕ=0 o a,m,s,ϕ + h m,0 = 1, ∀m = 0, . . . , M - C4: c 0 = C5: c a = c sh,sg a,ℓ + c sh,sg a,r , ∀a = 1, . . . , A M C6a: c a,i = c k if c a,i,k = 1, ∀a = 1, . . . , A M , i ∈ {ℓ, r} k = 0, . . . , a - C6b: a-1 k=1 c a,i,k = 1, ∀a = 1, . . . , A M , i ∈ {ℓ, r} C7a: c sh a,i = 2 s c a,i if φ a,i,s = ∀a = 1, . . . , A M , i ∈ {ℓ, r}, s = S min , . . . , S max C7b: Smax-1 s=Smin φ a,i,s = 1, ∀a = 1, . . . , A M , i ∈ {ℓ, r} C7c: φ a,ℓ,s = 0 ∀s > C7d: φ a,ℓ,s = φ a,r,s ∀s < C8a: c sh,sg a,i = -c sh a,i if ϕ a,i = 1, ∀a = 1, . . . , A M , i ∈ {ℓ, r} C8b: c sh,sg a,i = c sh a,i if ϕ a,i = 0, ∀a = 1, . . . , A M , i ∈ {ℓ, r} C8c: ϕ a,ℓ + ϕ a,r ⩽ 1, ∀a = 1, . . . , A M
depend on the number of zero filter coefficients, encoded by the binary decision variables h m,0 ). The resulting constraints can be roughly divided into frequency response conditions (C1, C2), equations linking the filter coefficients with the coefficients of the multiplier block (C3) and formulas describing the multiplierless realization of the multiplier block (C4 -C8).

The integer coefficients h ′ m (m = 0, . . . , M -1) of the FIR filter are directly used as integer variables in the ILP formulation. The resulting frequency response is constrained in C1a by setting (2) and ( 5) into [START_REF] Gustafsson | Design of linear-phase FIR filters combining subexpression sharing with MILP[END_REF]. Constraints C1b are called lifting constraints. These are actually not required to solve the problem, but can significantly reduce the search space and improve run-time performance. Specifically, they limit the range of the coefficients to lower h m and upper h m bounds. Constraint C2 similarly limits the range of the gain G. The computation of these bounds is considered in Section V-A.

Constraints C3a to C3c provide the connection between

A S ∈ N Number of structural adders h ′ m ∈ Z Integer representation of filter coefficient h m,0 ∈ {0, 1} true, if coefficient h ′ m is zero ca ∈ N Constant computed in adder a c a,i ∈ N Constant of input i ∈ {ℓ,
= 0 . . . M -1 G ∈ [G, G]
Gain of a flexible gain filter (G = 1 when the gain is fixed)

the filter coefficient h ′ m and the (potentially shifted and signcorrected) multiples computed in the multiplier block c a . For that, the binary decision variables o a,m,s,ϕ ∈ {0, 1} encode if h ′ m is connected to adder a of the multiplier block, shifted by s, and either added (ϕ = 0) or subtracted (ϕ = 1) in the structural adders (C3a). In case the coefficient is zero, a single binary decision variable h m,0 is used (C3b). This encoding allows the optimization of structural adders by considering the h m,0 variables in the objective function. For every zero coefficient, the corresponding structural adder(s) can be saved depending on the coefficient and filter type. Table III shows the number of structural adders for the different filter types. Overall, constraints C3c ensure that only one of the above cases is valid for each filter coefficient.

The remaining constraints C4 -C8 are identical to the ones used for solving the MCM problem from [START_REF] Kumm | Optimal Constant Multiplication using Integer Linear Programming[END_REF]. We give a brief description here, but refer the reader to [START_REF] Kumm | Optimal Constant Multiplication using Integer Linear Programming[END_REF] for a more detailed presentation. The multiplier block input is viewed as a multiplication by factor one (c 0 = 1) and is defined with constraint C4. Constraints C5 represent the actual add operation of adder a and its corresponding factor c a . It is obtained by adding the shifted and possibly sign corrected factors of its left input c sh,sg a,ℓ and its right input c sh,sg a,r . The source of the adder inputs is encoded by C6a/b. Indicator constraints C6a are used to set the value c a,i of the adder input i ∈ {ℓ, r} to the actual factor when the corresponding decision variable c a,i,k is set. Indicator constraints are special constraints in which a binary variable controls whether or not a specified linear constraint is active. They are in-fact non-linear, but are supported by many modern ILP solvers and are also simple to linearize for other solvers (see [START_REF] Kumm | Optimal Constant Multiplication using Integer Linear Programming[END_REF]). Constraints C6b make sure that only one source is selected. The actual shift is

I N -h 0,0 -2 M -1 m=1 h m,0 II N -2 M -1 m=0 h m,0 III N -2 M -1 m=1 h m,0 IV N -2 M -1 m=0 h m,0
constrained by C7a/b in a similar way: indicator constraints C7a are used to set the shifted factor c sh a,i according to the corresponding decision variable φ a,i,s .

Constraints C7c and C7d are both optional lifting constraints used to reduce the search space. As the filter coefficients can be shifted in constraint C3a, we can limit the constants of the multiplier block to odd numbers. This allows us to use the well-known fact that odd coefficients can be computed from odd numbers using one addition where either one operand is left shifted and the other operand is not shifted, or both operands are right shifted by the same value [START_REF] Dempster | Constant Integer Multiplication Using Minimum Adders[END_REF]Theorem 3].

To support subtractions, indicator constraints C8a/b are used to set the sign according to decision variable ϕ a,i . Finally, constraints C8c ensure that at most one input of the adder can be negative, as subtracting both inputs is typically more hardware demanding.

All of the integer variables from Table II are computed from boolean variables weighted by integer constants. They can thus be relaxed to real variables without changing the problem, while also speeding up the optimization (by reducing the number of integer variables).

B. Minimizing the Total Number of Adders

As the number of adders in the multiplier block A M is fixed in ILP Formulation 1, we need to iterate over various values A M to find the minimum number of total adders

A = A M + A S . (7)
For that, we first search for a solution with minimal number of multiplier block adders by solving ILP Formulation 1 using a lower bound for the multiplier block adders A M = A M,LB and, if infeasible, increase A M by one until we obtain the first feasible solution. A feasible lower bound A M,LB is zero, but a better lower bound is provided later in Section V-B. This solution with minimum multiplier block adders A M,min is not necessarily the global optimum, as there might be a solution with A M > A M,min and a smaller A S . To account for this, we need a lower bound for the structural adders A S,min . This is obtained once at the beginning of the overall algorithm by solving the problem for a maximally sparse FIR filter, which we do by taking ILP Formulation 1 where only the constraints C1 -C3 are considered.

In case the structural adders A S of solutions with A M = A M,min are not identical with A S,min , we increment A M until we find a solution with A S = A S,min . This is a safe stopping point since, by the optimality of ILP Formulation 1, there is no solution with larger A M and smaller A S . The solution with minimum total adders A found so far is hence also globally optimal. Typically, only a few iterations are necessary to reach this point.

IV. MULTIPLIERLESS FIR FILTERS WITH BOUNDED ADDER DEPTH

As discussed above, limiting the AD is important to reduce the power consumption of a filter. Unfortunately, adapting ILP Formulation 1 to limit the AD is not straightforward, as the topology of the adders and thus the AD is left open. We present in this section a novel ILP model for the design of multiplierless FIR filter with limited AD which is based on a formulation that was initially designed for optimizing pipelined MCM (PMCM) circuits [START_REF] Kumm | FIR Filter Optimization for Video Processing on FPGAs[END_REF].

In contrast to ILP Formulation 1, the possible coefficients are precomputed for each adder stage s and selected using binary decision variables. The computation of the corresponding coefficient sets is given next.

A. Definition of Coefficient Sets

We use some notation and definitions originally introduced in [START_REF] Voronenko | Multiplierless Multiple Constant Multiplication[END_REF]. First, we define the generalized add operation called A-operation, which includes shifts. An A-operation has two input coefficients u, v ∈ N and computes

A q (u, v) = |2 lu u + (-1) sv 2 lv v|2 -r , (8)
where q = (l u , l v , r, s v ) is a configuration vector which determines the left shifts l u , l v ∈ N 0 of the inputs, r ∈ N 0 is the output right-shift and s v ∈ {0, 1} is a sign bit which denotes whether an addition or subtraction is performed. Next, we define the set A * (u, v) containing all possible coefficients which can be obtained from u and v by using exactly one A-operation:

A * (u, v) := {A q (u, v) | q is a valid configuration} . (9)
A valid configuration is a combination of l u , l v , r and s v such that the result is a positive odd integer A q (u, v) ⩽ c max . The reason for limiting the integers to odd values is that we can compute every even multiple by shifting the corresponding odd multiple to the left. The c max limit is used to keep A * (u, v) finite. It is chosen as a power-of-two value which is usually set to the maximum coefficient bit width B plus one [START_REF] Dempster | Use of Minimum-Adder Multiplier Blocks in FIR Digital Filters[END_REF], [START_REF] Voronenko | Multiplierless Multiple Constant Multiplication[END_REF]

c max := 2 B+1 . ( 10
)
For convenience, the A * set is also defined for an input set X ⊆ N as

A * (X) := u,v∈X A * (u, v) . (11)
We can now define the coefficients that can be computed at adder stage s, denoted as A s , by recursively computing the A * sets

A 0 := {1} (12) A s := A * (A s-1 ) . (13)
In addition, let T s denote the set of (u, v, w) triplets for which w ∈ A s can be computed using u and v from the previous stage (i.e., u, v ∈ A s-1 ). T s can be computed recursively, starting from the last stage s, which is equal to the maximum allowable AD:

T s := {(u, v, w) | w = A q (u, v), u, v ∈ A s , u ⩽ v, w ∈ A s+1 }. (14)
To give an example, the first elements of T 1 are T 1 = {(1, 1, 1), (1, 1, 3), (1, 1, 5), (1, 1, 7), (1, 1, 9), (1, 1, 15), . . .}. This set contains all the possible rules for computing multiples from the input within one stage of additions, while set T 2 = T 1 ∪ { [START_REF] Lim | FIR filter design over a discrete powers-of-two coefficient space[END_REF][START_REF] Chen | A trellis search algorithm for the design of FIR filters with signed-powers-of-two coefficients[END_REF][START_REF] Aktan | An algorithm for the design of low-power hardware-efficient FIR filters[END_REF], [START_REF] Lim | FIR filter design over a discrete powers-of-two coefficient space[END_REF][START_REF] Yli-Kaakinen | A systematic algorithm for the design of multiplierless FIR filters[END_REF][START_REF] Aktan | An algorithm for the design of low-power hardware-efficient FIR filters[END_REF], . . . , [START_REF] Chen | A trellis search algorithm for the design of FIR filters with signed-powers-of-two coefficients[END_REF][START_REF] Yli-Kaakinen | A systematic algorithm for the design of multiplierless FIR filters[END_REF][START_REF] Aktan | An algorithm for the design of low-power hardware-efficient FIR filters[END_REF], . . .} contains all the combinations of how elements in the next stage can be computed.

B. ILP Formulation for Fixed Adder Depth

The bounded AD model is given in ILP Formulation 2, while the corresponding constants and variables are given in Table IV.

In contrast to ILP Formulation 1, the objective is to directly minimize the total number of adders A, which is separated into adders in the multiplier block (A M ) and structural adders A S . Similar to ILP Formulation 1, the constraints are divided into frequency response conditions (C1, C2), the link between the filter coefficients and the coefficients of the multiplier block (C3, C4) and the equations describing the multiplierless realization of the multiplier block (C5 -C8).

Constraints C1a/b and C2 are identical to the ones in ILP Formulation 1. Now, the connection between the odd multiplier block coefficients of the pre-computed sets and the filter coefficients is performed using binary decision variables. Let h m,w ∈ {0, 1} be a binary decision variable that is true if the magnitude of h ′ m is identical to w, i.e.,

h m,w = 1 when |h ′ m | = w 0 otherwise ( 15
)
for m = 0, . . . , M -1 and w = 0, . . . , 2 B -1. Furthermore, let ϕ m determine the sign of h ′ m as follows

ϕ m = 0 when h ′ m ≥ 0 1 otherwise . ( 16
)
The value of each integer coefficient h ′ m is selected by the indicator constraints C3a. In addition, constraints C3b make sure that only one value per filter coefficient is selected.

Next, we distinguish between coefficients that are computed for the selected stage (by using an addition) and coefficients that are just duplicated from a previous stage. Hence

G2 B D(ω) ⩽ M -1 m=0 h ′ m c m (ω) ⩽ G2 B D(ω), ∀ω ∈ Ω d C1b: h m ⩽ h ′ m ⩽ h m , ∀m = 0, . . . , M -1 C2: G ⩽ G ⩽ G C3a: h ′ m =              2 B -1 w=0 wh m,w if ϕ m = 0 - 2 B -1 w=1 wh m,w if ϕ m = 1 ∀m = 0, . . . , M - 1
C3b:

2 B -1 w=0 h m,w = 1, ∀m = 0, . . . , M -1 C4: d AD odd(w) + a AD odd(w) ⩾ 1 M M -1 m=0 h m,w , ∀w = 0, . . . , 2 B -1 C5: d s w = 0 ∀w ∈ A s \ s-1 s ′ =0
A s ′ and s = 1, . . . , AD -1

C6: d s w -a s-1 w -d s-1 w ⩽ 0, ∀w ∈ A s \{0}, s = 2, . . . , AD C7: a s w - (u,v,w ′ )∈T s | w ′ =w x s-1 (u,v) ⩽ 0, ∀w ∈ A s , s = 2, . . . , AD C8: x s (u,v) -d s u -a s u ⩽ 0 x s (u,v) -d s v -a s v ⩽ 0 ∀(u, v,
w) ∈ T s with s = 1, . . . , AD -1 using an adder or just duplicated, respectively. Duplicates do not contribute to the cost in this work, but note that the decision variable can be used to model the cost in a pipelined implementation of the multiplier block.

The connection to the filter coefficients h m,w is made through C4. As several of the M h ′ m coefficients can have the same w value, the right hand side of C4 is scaled by 1/M to keep it less than one. Whenever the right hand side of C4 is non-zero it forces the realization of coefficient w in the output stage AD, either as an adder or as a register/wire.

Constraints C5 and C6 consider the realization as register/wire: they require that a value w can only be replicated from a previous stage if it was computed or replicated before.

The realization as an adder computing constant w from the inputs u and v requires the presence of both inputs in the previous stage. For that, the binary variables x s (u,v) are introduced to determine if both are available in stage s:

x s (u,v) = 1 if
G ∈ [G, G]
Gain of a flexible gain filter (G = 1 when the gain is fixed)

Now, constraint C7 specifies that if w is computed by w = A(u, v) in stage s, the pair (u, v) has to be available in the previous stage. If a pair (u, v) is required in stage s, constraints C8 make sure that u and v have been realized in the previous stage either as register or adder. Note that instead of using constraint C5, it is more practical to remove all variables d s w = 0 from the cost function and their related constraints. Also note that the binary variables x s (u,v) and the integer variables h ′ m can be relaxed to real numbers to speed up the optimization.

Finally, this model can be extended to optimize an approximate low-level cost. For that, the A S and a s w terms have to be weighted with their corresponding low-level cost. While this is exact for the multiplier blocks (as their coefficients are pre-computed) the costs for the structural adders have to be approximated. This can be done by taking the input word size for computing the cost of each structural adder and considering their word size increase due to the multiplier block in the a s w terms. This works, of course, only when the structural adder is actually realized.

C. Selecting the Adder Depth

The AD is often selected to be as small as possible, typically at the expense of a higher adder cost. It is well known that the minimal AD needed when multiplying with a given coefficient can be realized by using a binary tree [START_REF] Kang | FIR Filter Synthesis Algorithms for Minimizing the Delay and the Number of Adders[END_REF]. Therefore, it cannot be lower than the base two logarithm of the nonzero digit count of its MSD representation. Unfortunately, as the coefficients are not known in advance, the minimum AD cannot be derived from the filter specification. However, the upper bound of the AD can be computed from the coefficient word size B as follows. A B bit binary number can have up to B + 1 digits when represented as an MSD number and up to ⌊(B + 1)/2⌋ + 1 non-zeros in the worst case. This leads to a maximum adder depth of

AD max = log 2 B + 1 2 + 1 . ( 18
)
Using this bound, a search from AD = 0, . . . , AD max can be performed until the first feasible solution is found. For practical FIR filters, early studies have shown that coefficient word sizes between 15 bit to 20 bit are sufficient to achieve approximation errors between -70 and -100 dB. Using [START_REF] Ye | Bit-level multiplierless FIR filter optimization incorporating sparse filter technique[END_REF], this translates to ADs of at most three to four. In our experiments, we found very good solutions with AD = 2 for most of the filters from practice.

V. REDUCING THE PROBLEM COMPLEXITY A. Reducing the Coefficient Range

Following [START_REF]Exact and Approximate Algorithms for the Filter Design Optimization Problem[END_REF]Sec. 3], we bound the search space for the gain and coefficient values, respectively, by projecting the polytope corresponding to the discretized versions of ( 4) or ( 6) onto G and each h ′ m . For instance, in the case of the coefficients, the goal is a tight interval enclosure [h m , h m ] for the feasible values of h ′ m . This corresponds to the LPs:

minimize h ′ m or maximize h ′ m subject to G2 B D(ω) ⩽ M -1 k=0 h ′ k c k (ω) ⩽ G2 B D(ω), ∀ω ∈ Ω d ,
where h ′ k ∈ R for k = 0, . . . , M -1 and G ∈ G, G (or G = 1 when unity gain is used). We get [h m , h m ] by taking

h m = ⌈h ′ m ⌉ from minimize h ′ m , h m = ⌊h ′ m ⌋ from maximize h ′ m .
We do the same for G in computing the enclosure [G, G], which is done before bounding the filter coefficients.

B. Computing a Lower Bound for the Number of Adders

Having a good lower bound for the number of adders helps to reduce the number of ILP runs for ILP Formulation 1 as described in Section III-B. There are well-known bounds for the number of adders that are required to multiply with a pre-determined set of constants [START_REF] Gustafsson | Lower Bounds for Constant Multiplication Problems[END_REF]. Unfortunately, these can not be directly applied as we do not know the coefficients in advance. However, we can adopt them to the ranges we obtained in Section V-A.

The initial lower bound in [START_REF] Gustafsson | Lower Bounds for Constant Multiplication Problems[END_REF] is the following:

A M,LB = min m ⌈AD min (h ′ m )⌉ + M uq -1 (19)
The function AD min (h ′ m ) computes the minimum AD for the coefficient h ′ m and M uq is the number of positive odd unique coefficients excluding zero and one. To transfer this to the coefficient ranges we obtained, we have to compute lower bounds for AD min (h ′ m ) and M uq . A lower bound of AD min (h ′ m ) can be obtained by simply evaluating the minimum AD for all the possible values h

′ m ∈ [h m , h m ], i.e., A M,LB = min m,h ′ m ∈[hm,hm] ⌈AD min (h ′ m )⌉ + M uq -1 (20)
To obtain a lower bound of the number of unique coefficients M uq , we compute the following: 1) We initialize M uq = 0 and O m = {}, m = 0, . . . , M -1.

2) For each h ′ m ∈ [h m , h m ], we compute its positive odd representation by dividing its absolute value by two until it is odd and add it to the O m set.

3) For each O m , we check whether it contains neither 0 nor 1 and does not have any intersection with another set O n with n > m. In this case, we know that we need at least one adder to compute this coefficient and we increase M uq by one.

C. Discretizing the Frequency Domain

Even though Ω is replaced by a finite set Ω d , we perform a rigorous posteriori validation of the result over Ω [START_REF] Volkova | Reliable verification of digital implemented filters against frequency specifications[END_REF]. Still, the typically large size of Ω d (16M is a common value found in the literature) can have a big impact on the run-time of the filter design routine. This is shown for instance in the context of an optimal branch-and-bound algorithm for FIR filter design with fixed-point coefficients [START_REF] Kodek | Design of optimal finite wordlength FIR digital filters[END_REF]Table 2]. A too small number of points can, on the other hand, lead to an invalid solution over Ω and a larger feasible set, potentially incurring a larger run-time as well. It is thus important to consider a discretization of reasonable size that is unlikely to lead to invalid solutions over Ω (i.e., equations ( 4) or ( 6) do not hold) and does not increase the search space by a too large factor. To this effect, we use so-called approximate Fekete points (AFPs), which contain the most critical frequencies for a given filter that needs to fit a target frequency response. They have recently been used to improve the robustness of the classic Parks-McClellan Chebyshev FIR filter design algorithm [START_REF] Filip | A robust and scalable implementation of the Parks-McClellan algorithm for designing FIR filters[END_REF] and for a fast and efficient heuristic for FIR fixedpoint coefficient optimization [START_REF] Brisebarre | A Lattice Basis Reduction Approach for the Design of Finite Wordlength FIR Filters[END_REF]. They are efficient choices when performing polynomial interpolation/approximation on domains such as Ω. This is relevant in our context since H R (ω) in ( 2) is a polynomial in cos(ω). For details on how to compute them we refer the reader to [START_REF] Filip | A robust and scalable implementation of the Parks-McClellan algorithm for designing FIR filters[END_REF], [START_REF] Brisebarre | A Lattice Basis Reduction Approach for the Design of Finite Wordlength FIR Filters[END_REF] and the references therein.

D. An Adaptive Search Strategy

Even if the current Ω d leads to a solution that does not pass a posteriori validation, it might still be possible to rescale the gain factor G such that (6) holds. By taking a point ω max ∈ Ω where GD(ω max ) -H R (ω max ) or H R (ω max ) -GD(ω max ) is largest (i.e., the point of largest deviation from the specification) we first update G to take a value close to H R (ω max )/D(ω max ) or H R (ω max )/D(ω max ), depending on where the deviation occurs. If this new gain leads to a valid solution over Ω, then it is optimal. If not, we update Ω d by adding the points of largest deviation for each frequency subdomain. We rerun the optimization with this new Ω d , repeating until either (a) there are no more invalid frequency points or (b) the problem becomes infeasible, meaning no solution with the imposed constraints over Ω exists.

We should mention that running the result validation code of [START_REF] Volkova | Reliable verification of digital implemented filters against frequency specifications[END_REF] at each iteration of the adaptive routine is computationally expensive. This is why at each iteration we perform a fast, non rigorous test consisting of verifying (6) on a much denser discretization of Ω than Ω d . We found this to usually be sufficient in ensuring that the a posteriori validation [START_REF] Volkova | Reliable verification of digital implemented filters against frequency specifications[END_REF] done at the end of optimization is successful. This is in stark contrast with the rest of the literature, which generally only uses a small discretization of Ω throughout the design process. While good for performance, such an approach will sometimes lead to designs which actually fail to satisfy the specification (see results in Table VI).

VI. EXPERIMENTAL RESULTS

To test the ILP formulations discussed above, we have implemented them in a C++ filter design tool 1 . It features a flexible command-line interface.

A. Experimental Setup and Parameter Choices

The proposed implementation supports several popular (M)ILP solvers, such as Gurobi2 and CPLEX 3 . For convenience, these solvers are accessed through the ScaLP [START_REF] Sittel | ScaLP: A Light-Weighted (MI)LP Library[END_REF] library, which acts as a frontend. Based on our experiments, Gurobi usually proved to be the fastest, which is why, apart from a few exceptions, we use it on all the examples below.

All experiments use the AFP-based frequency grid discretization mentioned in Section V-C. As discussed before, the number of frequency points in Ω d is run-time critical. To determine an appropriate size, we ran an experiment using a typical design scenario with an Ω d size of kM points and k = 1, . . . , 32. They start large for very low k, as in these cases the frequency grid usually has to be extended to address violations, which require re-running the optimization routine on a new grid. If k is large enough (e.g., k ⩾ 4), invalid results become rare, meaning just one optimization pass is sufficient. Further increasing k at this point just leads to more constraints in the model and likely a larger run-time. Based on these results, we start with 4M points. This choice usually delivers a good balance between optimizer run-time and number of iterations needed to obtain a valid solution over Ω.

B. Benchmark Set

Several multiplierless filter designs were computed to evaluate our methods. They are introduced next.

1) A Family of Specifications from [4, Example 1]: We consider a family of low-pass linear-phase filter specifications from Redmill et al. [START_REF] Redmill | Genetic synthesis of reduced complexity filters and filter banks using primitive operator directed graphs[END_REF]. These specifications are defined by: where δ is a parameter regulating error. We set δ = 10 -p 20 , where p > 0 is the error in decibels (dB).

1 -δ ⩽ H R (ω) ⩽ 1 + δ, ω ∈ [0, 0.3] (passband) -δ ⩽ H R (ω) ⩽ δ, ω ∈ [0.5, 0.1] (stopband)
Our goal with this benchmark is to explore the tradeoff between the error (p), the filter order (N ) and the word length (B) in terms of the total number of adders.

2) A Set of State-of-the-art Specifications: We also test our tool on a set of reference specifications from the literature [START_REF] Lim | FIR filter design over a discrete powers-of-two coefficient space[END_REF], [START_REF] Samueli | An improved search algorithm for the design of multiplierless FIR filters with powers-of-two coefficients[END_REF], [START_REF] Redmill | Genetic synthesis of reduced complexity filters and filter banks using primitive operator directed graphs[END_REF], [START_REF] Gustafsson | Design of linear-phase FIR filters combining subexpression sharing with MILP[END_REF], [START_REF] Yu | Design of Linear Phase FIR Filters in Subexpression Space Using Mixed Integer Linear Programming[END_REF], [START_REF] Xu | Design of Low-Complexity FIR Filters Based on Signed-Powers-of-Two Coefficients With Reusable Common Subexpressions[END_REF], [START_REF] Shi | Design of Linear Phase FIR Filters With High Probability of Achieving Minimum Number of Adders[END_REF], [START_REF] Shahein | A Novel Hybrid Monotonic Local Search Algorithm for FIR Filter Coefficients Optimization[END_REF], [START_REF]Exact and Approximate Algorithms for the Filter Design Optimization Problem[END_REF], referred to as S1, S2, L1, L2, L3, X1, G1, Y 1 and Y 2. They are all low-pass filters defined by

1 -δ p ⩽ H R (ω) ⩽ 1 + δ p , ω ∈ Ω p (passband), -δ s ⩽ H R (ω) ⩽ δ s , ω ∈ Ω s (stopband),
where the values of δ p , δ s , Ω p , Ω s for each specification are given in Table V. Over time, some of these reference filter specifications were slightly modified by different publications.

We compare against the most recent version in each case, updating it to account for any violation of the specification in the results reported in the literature. The restriction to low-pass filters comes only from the existing literature. Our tool can also be used to design other types of filters, such as multiband filters or decimators (since we generalize constraints on the frequency response as functions of frequencies).

C. 1-step vs 2-step Design Optimization Approach

We start by demonstrating the advantages of an overall 1step optimization over the classic 2-step filter design process. For each frequency specification, there exist numerous coefficient sets satisfying it. Take, for example, frequency specification S1, when realized as a type I filter with N = 24: there exist 237 sets of coefficients of word length B = 9 that satisfy the constraints. We brute-force explored the design space and applied the optimal MCM solver [START_REF] Kumm | Optimal Constant Multiplication using Integer Linear Programming[END_REF] on each possible filter coefficient set to design the multiplierless implementations. For this particular example, the total adder cost varies from 26 up to 34 adders. Figure 4 presents the histogram illustrating the number of coefficient sets falling into each adder cost category. It can be interpreted in the following way: when selecting coefficient sets out of the pool of feasible results, with high probability one obtains the design of cost 32, 33 or perhaps 31 adders. Indeed, only 15 out of 237 filters have cost less than 30 adders and the lowest cost is achieved by only one filter. Generally speaking, for an arbitrary filter specification, the 2-step approach that first selects a filter coefficient set (without a priori knowledge of the cost) and then designs an optimal MCM architecture, has a high probability to be far from optimal. Brute-force design exploration as presented here is not practical for higher-order and high-wordlength solutions, hence a 1-step optimization approach is necessary to navigate the search space towards an optimal solution.

D. ILP1 versus ILP2

1) Model Complexity and Run-Time Comparison: To estimate the model complexity, we evaluated the model sizes as well as the run-times. This is done using the family of specifications from [4, Example 1] as described in Section VI-B1. The value of p is varied from -2 dB to -46 dB using an effective word length of up to B = 10 bits. A timeout of 2 hours was set for each ILP run. Fig. 5 shows the ILP model sizes in terms of number of variables and constraints, while Fig. 6 shows the run-times obtained on an 8-core Intel Core i9 notebook CPU. The run-times include the ILP runs of the coefficient range reduction (Section V-A) as well as the adaptive search strategy (Section V-D). Coefficient range reduction always took less than a second to compute the lower and upper bounds for each coefficient. For the adaptive search strategy, at least one re-run was necessary in 27 out of the 125 cases (21.6%). At most four re-runs were necessary to meet frequency specification, in two of the cases.

The following can be observed: for ILP1, the model size grows with decreasing error as more adders are required, leading to longer run-times. For errors smaller than -38 dB, the ILP timeout was exceeded. For ILP2, the model size remains approximately constant for a fixed AD. The fluctuations in model size can be explained by the coefficient range reduction of Section V-A.

For AD=2, fast run-times were observed, outperforming ILP1 and scaling to -46 dB errors within a few seconds (for smaller errors, the word size has to be increased). For AD=3, This shows the limitations of the proposed methods: ILP1 execution is dominated by the total number of adders while ILP2 execution is dominated by the adder depth. As filters with low adder depth are desired for low power applications, this later limitation is less impactful.

2) Optimization Results Comparison: Both ILP1 and ILP2 can be used to optimize for the total number of adders given fixed parameters like filter order N , filter type and the effective word length B (see Fig. 3). In case of ILP2, the adder depth is an additional constraint. Therefore, in practice, the two approaches can sometimes lead to different results. This is exemplified in Fig. 7, where we again design a set of filters using the family of specifications from [4, Example 1]. We consider 37 filters with error varying from -2 dB to -38 dB, with a 9-bit effective word length and gain G = 1. In each case, a type I filter with smallest N that leads to a feasible Fig. 7: Total adder count comparison when using ILP1 (adapted to minimize the total number of adders) and ILP2 (with AD limit set to 2) on the Redmill set of filters [START_REF] Redmill | Genetic synthesis of reduced complexity filters and filter banks using primitive operator directed graphs[END_REF]. solution under the given constraints was used. For ILP2, the upper bounds on the AD were set to 2 and 3, respectively. The ILP solver timeout was again set to 2 hours.

For most error targets the total adder count is identical between the three variants. The exceptions are as follows:

• At -22 dB the minimum filter order is N = 10, leading to a higher adder count than for -23 dB. . .-26 dB where a filter order of N = 12 is necessary. • At -27 dB and -32 dB, AD = 2 was simply not sufficient to reach the minimum adders. Similarly, at -31 dB and -36 dB, AD = 3 was not sufficient. • At -37 dB, ILP1 finds an optimal solution with N = 20 and 25 adders while an AD = 2 solution for ILP2 is only possible starting with N = 22 and achieving 23 adders. As we further illustrate in the next section, the combination of the filter order, word length and adder depth for a given specification forms a non-linear discrete design-space, increasing the difficulty of the search for the overall minimal-adder implementation.

E. Optimization Results for Benchmark Sets

For comparison with previous work, we use ILP2 with AD = 2 (unless otherwise stated).

1) The Family of Specifications from [4, Example 1]: The experiment setting from Section VI-D2 is expanded upon. We compare our best results (with effective word lengths B ∈ {8, 9, 10, 11}) with those from [4, Example 1]. We start off by considering only type I filters (just like in [START_REF] Redmill | Genetic synthesis of reduced complexity filters and filter banks using primitive operator directed graphs[END_REF]), flexible gain G ∈ [2/3, 4/3] and minimal order N for each error target. The results are illustrated in Fig. 8. We note that there are certain cases where, for a given B, taking the minimal filter order leading to a feasible solution does not minimize the adder cost. This is most visible for B = 11 and a -30 dB error target, where a minimal order N = 14 filter requires 24 adders. For B = 10, the minimal N is 16, leading to only 17 total adders, a 7 adder improvement. Taking N = 16 for B = 11 also results in a 17 adder solution. A lower implementation cost Total number of adders best results extracted from [START_REF] Redmill | Genetic synthesis of reduced complexity filters and filter banks using primitive operator directed graphs[END_REF] best results with fixed gain best results with flexible gain Fig. 9: Comparison between our best design space exploration results and the best results from [4, Example 1]. Our tool improves designs from [START_REF] Redmill | Genetic synthesis of reduced complexity filters and filter banks using primitive operator directed graphs[END_REF] or proves them optimal. is sometimes possible when increasing the filter order leads to a sparser filter and/or a more economical MCM design. Such solutions better optimize the objective functions in the proposed ILP models. We nevertheless remark that increasing the filter order beyond a certain threshold will not lead to different solutions.

Of course, increasing the word length can also lead to a significant improvement in the results. For instance, the optimal -50 dB attenuation results for B ∈ {9, 10, 11} require 42, 32 and 31 adders, respectively. This nonlinearity of the word length/cost relation means that the user should favor a comprehensive exploration of the design space, varying the design parameters (especially B, filter type and N ) and examine the various trade-offs. This is possible with our tool. Fig. 9 shows the results of such an experiment where, with respect to the setting of Fig. 8, we additionally allow N to vary and also consider type II filters. We also added the fixed and flexible gain as well as the results from the genetic algorithm presented in [START_REF] Redmill | Genetic synthesis of reduced complexity filters and filter banks using primitive operator directed graphs[END_REF]. Compared to [START_REF] Redmill | Genetic synthesis of reduced complexity filters and filter banks using primitive operator directed graphs[END_REF], we could improve all of the results, except two of them where we obtained the same adder cost (-9 dB and -25 dB). It is clearly visible that allowing flexible gain designs can have a major influence on the quality of the results.

2) The Set of State-of-the-art Specifications: Table VI presents the comparison between our designs and the best results from literature [START_REF] Lim | FIR filter design over a discrete powers-of-two coefficient space[END_REF], [START_REF] Samueli | An improved search algorithm for the design of multiplierless FIR filters with powers-of-two coefficients[END_REF], [START_REF] Redmill | Genetic synthesis of reduced complexity filters and filter banks using primitive operator directed graphs[END_REF], [START_REF] Yu | Design of Linear Phase FIR Filters in Subexpression Space Using Mixed Integer Linear Programming[END_REF], [START_REF] Shi | Design of Linear Phase FIR Filters With High Probability of Achieving Minimum Number of Adders[END_REF], [START_REF] Shahein | A Novel Hybrid Monotonic Local Search Algorithm for FIR Filter Coefficients Optimization[END_REF] for the filter specifications from Table V. The following information for each implementation is given: filter order (N ), filter type, number of multiplier adders (A M ), number of structural adders (A S ), total number of adders (A), adder depth of the multiplier block (AD), gain (G), effective coefficient word length (B) and, finally, coefficients of the filter. In the following, we discuss each instance in detail.

S1: for this specification we show that the AD can be reduced from 3 to 2 stages, while keeping the same effective word length. We also reduce the number of adders by two (one structural and one MB).

S2: in [START_REF] Yu | Design of Linear Phase FIR Filters in Subexpression Space Using Mixed Integer Linear Programming[END_REF], implementations with adder depths 3 and 2 are proposed, at the cost of 78 and 80 adders, respectively. These results are improved in [START_REF] Shi | Design of Linear Phase FIR Filters With High Probability of Achieving Minimum Number of Adders[END_REF], with the authors claiming that a 3-stage implementation at the cost of 76 adders has a high probability to be optimal. We demonstrate that a 2-stage design with a cost of only 66 adders is possible. Again, our result has higher sparsity than previous designs.

L1: the tool timed out for this specification before giving a feasible result, hence it is not presented in Table VI. The best known result from literature is a 120-tap filter [START_REF] Yu | Design of Linear Phase FIR Filters in Subexpression Space Using Mixed Integer Linear Programming[END_REF] and the size of an instance of the corresponding ILP formulation goes beyond the current capabilities of the solvers we tried, showing its limitations.

L2: we could not achieve a better solution than the state of the art [START_REF] Shi | Design of Linear Phase FIR Filters With High Probability of Achieving Minimum Number of Adders[END_REF] before the timeout. This example showcases the scalability limits of our tool with a high filter order and AD = 3, as also shown in Section VI-D1.

L3: in [START_REF] Shi | Design of Linear Phase FIR Filters With High Probability of Achieving Minimum Number of Adders[END_REF], an implementation with 35 adders is provided. We show that a 34 adder solution is possible.

X1&G1: we obtain the same results as the optimal ones reported in [START_REF]Exact and Approximate Algorithms for the Filter Design Optimization Problem[END_REF].

Y1: the verification of the frequency response of the result of [START_REF]Exact and Approximate Algorithms for the Filter Design Optimization Problem[END_REF] revealed a slight stopband violation. By adjusting the frequency specification (denoted as Y1 * ), we could obtain the same result. The solution of the initial frequency specification with AD= 2 requires one adder more while setting AD= 3 as in [START_REF] Shi | Design of Linear Phase FIR Filters With High Probability of Achieving Minimum Number of Adders[END_REF] yield to an optimal solution with 29 adders.

Y2: we obtain a result with same adder count as the optimal one reported in [START_REF]Exact and Approximate Algorithms for the Filter Design Optimization Problem[END_REF].

Overall, the proposed tool achieves improvements to several of the considered filter design problems and confirms optimal results reported in literature. Most importantly, the user can explore a large design space by setting different implementation parameters, e.g., adder depth, coefficient word length, filter type, etc. The required run-time however, will depend greatly on the problem. In the case of Table VI, it varied from several seconds for the smallest filters (S1) up to several days for the largest ones (S2 and L2).

VII. CONCLUSION AND PERSPECTIVES

In this paper we have introduced two new algorithms for the design of optimal multiplierless FIR filters. Relying on ILP formulations stemming from the MCM literature, our algorithms minimize either (a) the number of structural adders given a fixed budget of multiplier block adders or (b) the total number of adders (multiplier block + structural adders) given a fixed adder depth. We further show how (a) can be applied iteratively to optimally minimize the total number of adders (without any adder count or adder depth constraints). Extensive numerical tests with example design problems from the stateof-the-art show that our approaches can offer in many cases better results. We also make available an open-source C++ implementation of the proposed methods.

In the future we plan to extend the current solution to other filter structures and use other efficiency metrics, in particular to optimize the number of full adders. Regarding new metrics, this calls for modeling the impact of the precision choice on the numerical quality of the result, which is non-trivial due to its highly non-linear nature. Developing a dedicated branch and bound solver for instance, inspired by [START_REF]Exact and Approximate Algorithms for the Filter Design Optimization Problem[END_REF], has the potential of helping us deal with this heterogenous design problem. Regarding filter structures, we plan to adapt our framework to cascaded forms and recursive filters.

Fig. 1 :

 1 Fig. 1: Different adder circuits for constant multiplications.

Fig. 2 :

 2 Fig. 2: Structure of FIR filters.

Fig. 3 :

 3 Fig. 3: Simplified multiplierless FIR filter design flow.

Fig. 4 :

 4 Fig. 4: Histogram demonstrating the total adder cost of feasible filter designs for S1 specification with word length B = 9.

× 10 6 1 × 10 7 1 × 10 8 Fig. 5 :

 185 Fig. 5: Model size in terms of the number of variables and constraints for different pass/stopband errors for ILP1 and ILP2 with AD=2 and AD=3.

Fig. 6 :

 6 Fig. 6: Run-time for different pass/stopband errors for ILP1 and ILP2 with AD=2 and AD=3.

  pass/stopband (dB) Total number of adders best results for B = 8 best results for B = 9 best results for B = 10 best results for B = 11

Fig. 8 :

 8 Fig. 8: Our designs with effective word lengths varying from 8 to 11, filter type I and smallest feasible filter order.

TABLE I :

 I Relation between filter order N , number of coefficients M and function c m (ω) for different filter types

	Type Sym.	N	M	cm(ω)

TABLE II :

 II Used constants (top) and variables (bottom) in ILP Formulation 1

	Constant/Variable	Meaning
	A M ∈ N	Number of adders in the multiplier block
	M ∈ N	Number of filter coefficients
	S min , Smax ∈ Z	Minimum and maximum shift
	B ∈ Z	max. effective word length of coefficients
	G, G ∈ R	lower and upper bounds on the filter gain

TABLE III :

 III Structural adder count for the different filter types Type no. of structural adders, A S (h m,0 )

TABLE V :

 V Specifications of the reference filters

	Name	Source	Ωp/π	Ωs/π	δp	δs
	S1	[4]	[0, 0.3]	[0.5, 1]	0.00636 0.00636
	S2	[9], [14], [23]	[0, 0.042]	[0.14, 1]	0.026	0.001
	L1	[9], [11], [23]	[0.8, 1]	[0, 0.74]	0.0057	0.0001
	L2	[1], [16], [23]	[0, 0.2]	[0.28, 1]	0.02800 0.001
			[0; 0.15]		0.0165
	L3	[1], [16]	[0.15; 0.1875]	[0.2875; 1]	0.0296	0.0316
			[0.1875; 0.2125]		0.0546
	X1	[10], [14], [23]	[0, 0.2]	[0.8, 1]	0.0001	0.0001
	G1	[6], [14], [23]	[0, 0.2]	[0.5, 1]	0.01	0.01
	Y1	[14], [23]	[0, 0.3]	[0.5, 1]	0.00316 0.00316
	Y2	[14], [23]	[0, 0.3]	[0.5, 1]	0.00115 0.00115

TABLE VI :

 VI Comparison between our method and the state-of-the-art results for the specifications in Table V. -5 -13 0 24 20 -26 -55 0 91 73 -96 -212 0 534 1004 * The solution slightly fails the specification

	Name	Source	N	Type A M	A S	A	AD	G	B	Coefficients
	S1	[4]	24	I	6	20	26	3	2.4570	9	2 8 0 -16 -14 20 43 0 -80 -71 112 377 502
	S1	ours	23	II	5	19	24	2	2.46492	9	6 6 -8 -21 0 36 32 -42 -96 0 248 472
	S2	[14]	59	II	17	59	76	3	10.47032 10 5 5 6 5 2 -2 -10 -20 -32 -48 -64 -78 -92 -98 -87 -65
											-26 26 93 174 267 368 472 575 672 757 826 874 898
	S2	ours	59	II	15	51	66	2	7.5904	10 0 0 0 -2 -5 -10 -16 -23 -32 -40 -50 -58 -64 -64 -61
											-50 -29 0 38 86 143 206 274 344 412 476 532 576 608 624
	L2	[14]	62	I	17	56	73	3	4.1991	10 4 8 12 13 9 0 -10 -16 -13 0 19 35 36 18 -15 -49 -64 -48 0
											60 102 96 32 -72 -170 -203 -124 79 371 678 911 998
	L2	ours	62	I	16	62	78	3	2.6668	11 4 9 13 12 4 -10 -26 -36 -32 -12 18 44 52 32 -10 -56 -80
											-64 -4 74 130 128 48 -86 -215 -263 -168 88 460 854 1153
											1265
	L3	[14]	35	II	4	31	35	1	2.627	7	3 0 -2 -5 -5 0 3 7 8 3 -6 -14 -16 -7 12 40 65 80
	L3	ours	35	II	3	31	34	2	2.60564	7	4 0 -2 -4 -4 0 3 8 8 3 -6 -13 -16 -8 13 40 64 80
	X1	[23], ours 14	I	5	8	13	2	1.6404	10	-4 0 28 0 -113 0 509 840
	G1	[23], ours 15	II	2	15	17	2	2.6338	6	1 2 -1 -7 -7 7 34 56
	Y1 *	[23], ours 29	II	6	23	29	2	2.505	9	-1 -4 0 8 8 -10 -22 0 40 33 -44 -99 0 254 479
	Y1	ours	29	II	7	23	30	2	2.468	10 -2 -8 0 16 15 -20 -44 0 80 65 -88 -195 0 500 944
	Y1	ours	29	II	6	23	29	3	2.472	10 -2 -8 0 16 15 -20 -44 0 80 64 -88 -196 0 501 945
	Y1	[14]	29	II	6	23	29	3	2.5985	10	-2 -8 0 17 16 -21 -46 0 84 68 -92 -205 0 527 994
	Y2	[23]	37	II	9	29	38	3	2.6361	10	-1 0 4 4 -5 -13 0 24 20 -26 -55 0 91 73 -96 -213 0 536 1008
	Y2	ours	37	II	9	29	38	3	2.6259	10	-1 0 4 4

Available as an open-source project at: https://gitlab.com/filteropt/firopt.

https://www.gurobi.com

https://www.ibm.com/analytics/cplex-optimizer