Silviu-Ioan Filip

A ROBUST AND SCALABLE IMPLEMENTATION OF THE PARKS-MCCLELLAN ALGORITHM FOR DESIGNING FIR FILTERS

With a long history dating back to the beginning of the 1970s, the Parks-McClellan algorithm is probably the most well-known approach for designing finite impulse response filters. Despite being a standard routine in many signal processing packages, it is possible to find practical design specifications where existing codes fail to work. Our goal is twofold. We first examine and present solutions for the practical difficulties related to weighted minimax polynomial approximation problems on multi-interval domains (i.e., the general setting under which the Parks-McClellan algorithm operates). Using these ideas, we then describe a robust implementation of this algorithm. It routinely outperforms existing minimax filter design routines.

Introduction

Digital filtering operations are essential to many engineering applications, which is why, over the last decades, they have enjoyed a sustained interest from researchers in Computer Science and Electrical Engineering. In general, a filtering framework consists of three major steps:

• specify and determine a mathematical representation of the filter (usually in terms of polynomials/rational functions); • quantize the values (i.e., coefficients) found at the previous step, using some imposed numerical representations (be it fixed-point or floating-point formats) [START_REF] Kodek | LLL Algorithm and the Optimal Finite Wordlength FIR Design[END_REF]; • synthesize the obtained filter in hardware/software [START_REF] De Dinechin | Sum-of-product architectures computing just right[END_REF]. The Parks-McClellan exchange algorithm [START_REF] Parks | Chebyshev Approximation for Nonrecursive Digital Filters with Linear Phase[END_REF] solves the first step. It is an iterative procedure converging to the optimal polynomial filter with real valued coefficients satisfying some given constraints.

Current implementations of the Parks-McClellan algorithm suffer from robustness issues. A typical instance where codes such as those found in MATLAB's Signal Processing toolbox, GNURadio or Scipy fail to work, is when computing a minimax filter of degree n = 100 adhering to the specification from Example 2.7. Such a degree is plausible in practice, yet it proves difficult to design. This unsatisfactory situation was the initial motivation of the present study. Together with the three other examples from Section 2.5, it will be used extensively to highlight the various contributions of this article.

The starting point of our work is [START_REF] Pachón | Barycentric-Remez algorithms for best polynomial approximation in the Chebfun system[END_REF], which describes a robust implementation of the Remez exchange algorithm for computing best polynomial approximations on a single compact interval, as part of the Chebfun [START_REF] Driscoll | [END_REF] MATLAB library (the remez command). The main elements that allow Chebfun remez to perform well in practice are: [START_REF] Abu-Al-Saud | Efficient wideband channelizer for software radio systems using modulated PR filterbanks[END_REF] the use of Chebyshev nodes (see Section 5.2.2 for a definition) as starting reference, (2) barycentric Lagrange interpolation [START_REF] Berrut | Barycentric Lagrange Interpolation[END_REF] and (3) a Chebyshev-proxy root-finding method [START_REF] Boyd | Finding the Zeros of a Univariate Equation: Proxy Rootfinders, Chebyshev Interpolation, and the Companion Matrix[END_REF][START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF] based on the recursive subdivision of the approximation domain. The Parks-McClellan algorithm computes minimax approximations in a weighted multi-interval context and can therefore be viewed as an extension of the case addressed by [START_REF] Pachón | Barycentric-Remez algorithms for best polynomial approximation in the Chebfun system[END_REF]. Unfortunately, adapting this implementation to our context is not at all straightforward and raises several significant issues. The present work addresses all of them and, to the best of our knowledge, it describes the first fast, robust and scalable implementation of the Parks-McClellan algorithm. Note also that, despite a central focus on filtering applications (the traditional setting for the Parks-McClellan algorithm), our code can compute more general weighted minimax polynomial approximations over a multi-interval compact subset of R (see Section 7.3).

We begin with an overview of how minimax approximation algorithms are related to digital filter design (Section 2) and practical difficulties in using them (Section 3). The main problem in writing an efficient implementation in our context is finding a suitable initial reference. In particular, there is no simple, closed-form equivalent for Chebyshev nodes in this setting (in the sense that interpolation at those points is close to the optimal approximation). Fast practical convergence of the exchange algorithm is very much dependent on placing the right number of reference points inside each interval of the domain. In Section 4, we describe two complementary heuristic solutions to this initialization problem. Our experiments show that they are quite satisfactory in practice. Section 5 looks at how initialization choices affect execution. We also identify possible sources of numerical instability at runtime, argue how they are related to each other and how to counter them. Central to our study is the Lebesgue constant (see Section 4.1 for a definition) that we use, in an apparently novel way, as a tool for detecting numerical problems. Another important difference with the Chebfun remez code is the fact that we eliminate recursion during root-finding by using problem-specific information about the root positions. Major advantages are a smaller number of computations compared to a recursive search and the possibility of introducing parallelism (see Sections 6 and 7), further speeding up execution on multicore systems.

The resulting code1 allows one to successfully address problematic instances in a scalable manner (see the examples from Sections 7.2 and 7.3). We offer three versions of our implementation. The first one uses IEEE-754 double-precision floating-point arithmetic, while the second one requires long double operations. On x86 machines and compilers, this format corresponds to 80-bit floating-point numbers. The other version uses MPFR [START_REF] Fousse | MPFR: A Multiple-precision Binary Floating-point Library with Correct Rounding[END_REF] multiple precision formats and serves mostly as a verification tool for our double and long double versions. All of the numerical tests in the sequel (and many others) are included with the source code. The double-precision version is sufficient for most practical uses. In our tests, only computations involving Example 2.9 required long double operations to obtain an accurate result.

An overview of minimax FIR filter design

A finite impulse response (FIR) filter with real-valued coefficients {b k } n2 k=-n1 is generally designed using a trigonometric polynomial H(ω) (the frequency response of the filter) of the form:

(1)

H(ω) = n2 k=-n1 b k e -ikω
Many problems in signal processing (for instance in data transmission, audio and image processing) require the b k 's to have real values and be (anti)symmetric around the index k = 0 (so-called linear-phase filters). Depending on the type of symmetry used, there are 4 major types (I to IV) of such filters [START_REF] Oppenheim | Discrete-Time Signal Processing[END_REF]Ch. 5.7.3]. Without any loss of generality, we will focus on type I filters, since type II to IV instances can be designed analogously [START_REF] Antoniou | Digital Signal Processing: Signals, Systems, and Filters[END_REF]]. The solution has a qualitative description: Theorem 2.3 (Alternation theorem). A necessary and sufficient condition for H(ω) to be the unique transfer function of degree at most n that minimizes the weighted approximation error δ Ω = E(ω) ∞,Ω is that E(ω) exhibit at least n + 2 equioscillating extremal frequencies over Ω; i.e., there exist at least n

+ 2 values ω k in Ω such that ω 0 < ω 1 < . . . < ω n+1 and E(ω k ) = -E(ω k+1 ) = λ(-1) k δ Ω , k = 0, . . . , n,
where λ ∈ {±1} is fixed.

Proof. See [ (5) Convergence Parameter Test: Remark 2.4. Algorithm 1 is usually called a multipoint exchange algorithm (or second Remez algorithm). Remarks about its convergence are made in [START_REF] Mcclellan | A personal history of the Parks-McClellan algorithm[END_REF], while in [START_REF] Antoniou | Digital Signal Processing: Signals, Systems, and Filters[END_REF]Ch 15.4] it is mentioned that, in practice, up to 12 iterations (Steps 2 to 4) are usually required for designing two and three-band FIR filters. In particular, [START_REF] Mcclellan | A personal history of the Parks-McClellan algorithm[END_REF] states that convergence is ensured if Ω is a finite set of points.

If max |E(ω k )|-min |E(ω k )| max |E(ω k )| ε t ,
Remark 2.5. Because of the equivalence to polynomial approximation, for much of the paper we will be using the set x = (x k ) 0 k n+1 = (cos(ω k )) 0 k n+1 in place of ω. We denote with X, where X ⊆ [-1, 1], the transformed frequency bands in this case. When talking about x, it will be as x = cos(ω). Also, we consider δ x = δ ω and δ X = δ Ω . [START_REF] Mcclellan | A computer program for designing optimum FIR linear phase digital filters[END_REF], was written in Fortran. It works on a modified version of the filter design problem. At each iteration, the new extremal set ω from Step 5 is chosen from a dense grid G of uniformly spaced points in Ω. The size of G is taken large enough (the initial value [START_REF] Parks | Chebyshev Approximation for Nonrecursive Digital Filters with Linear Phase[END_REF] was 20n), so that the obtained solution will be close to the one over Ω. Computing ω is done by an exhaustive evaluation of E(ω) over G.

Implementations. The first implementation of the Parks-McClellan algorithm

This early Fortran program is a starting point for current codes, like the firpm function from MATLAB's Signal Processing Toolbox and the variations from [START_REF] Ahsan | Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design[END_REF], or similar routines from Scilab, SciPy and GNURadio.

2.4. The role of transition bands. A typical use case for the exchange method is the design of lowpass filters, like the one in Figure 1, with Ω = [0,

ω p ] ∪ [ω s , π], D(ω) = 1, 0 ω ω p , 0, ω s ω π, and
W (ω) = δ2 δ1 , 0 ω ω p , 1, ω s ω π,
where δ 1 and δ 2 denote the maximum errors allowed over [0, ω p ] (the passband ) and [ω s , π] (the stopband ). Between pass and stopbands, we find transition bands. They are crucial here because their absence incurs jump discontinuities in the functions to approximate (at ω = ω p = ω s ), giving rise to the unwanted Gibbs-Wilbraham phenomenon (see [START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF]Ch. 9] for a detailed description).

In practice, the frequency response of a minimax filter inside transition regions is required to be monotonic or at least not present overshoots [START_REF] Rabiner | Some Considerations in the Design of Multiband Finite-Impulse-Response Digital Filters[END_REF]: that is, on a transition band, H(ω) can only take values which are inside the interval determined by its minimum and maximum on the neighboring passband and stopband. Such restrictions are implicitly assured in most practical cases, even though they do not hold in general [START_REF] Rabiner | Some Considerations in the Design of Multiband Finite-Impulse-Response Digital Filters[END_REF].

Figure 2 shows a three band, degree n = 38, type I minimax FIR filter with passband [0, 0.3π], stopband [0.33π, 0.5π], passband [0.6π, π] and corresponding weights 1, 10 and 2. A huge unwanted spike is clearly visible in the transition region between the stopband and the second passband. If we only need to contain it, we can adapt the initial specification using mild constraints on the problematic transition region [START_REF] Ahsan | Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design[END_REF]Sec. 3]. For instance, by adding [0.51π, 0.59π] to Ω and extending D and W to take the values 0.5 and 0.25 on [0.51π, 0.59π], we obtain the second minimax filter from Figure 2, where the transition region response is now between 0 and 1. There is only a minor loss in quality with respect to the initial design (the minimax error δ Ω grows from 0.1172 to 0.1205). If monotonicity is required, extensions of the exchange method that ensure it [START_REF] Shpak | A generalized Remez method for the design of FIR digital filters[END_REF] exist. The cost is, again, a larger δ Ω .

2.5.

Four key examples. Throughout the text, we will be frequently using four filter design problems to highlight our approach. They are taken and/or adapted from the literature and prove to be hard or impossible to solve using current de facto minimax implementations.

We The next problem comes from [START_REF] Zahradnik | Robust analytical design of equiripple comb FIR filters[END_REF] and is related to the design of equiripple comb FIR filters.

Example 2.8. A unit weight type I linear phase filter with passband [0, 0.99π], stopband centered at π and degree n = 520.

Our last example is derived from [START_REF] Abu-Al-Saud | Efficient wideband channelizer for software radio systems using modulated PR filterbanks[END_REF]. There, the exchange algorithm acts as an intermediary step in designing efficient wideband channelizers. The lengths of the filters they use are proportional to the number of desired channels. In order to design a 8192-subchannel system, we can use: Example 2.9. A degree n = 13 • 4096 = 53248 unit weight type II lowpass filter with passband 0, 1 8192 π and stopband 3 8192 π, π .

Practical problems

Over the years, certain difficulties with implementing multi-interval exchange methods have been discussed in the literature. We give short accounts on the most important ones, previous solutions and how we address them in the sequel.

3.1. Robustness issues: initialization is a key step. The need for high-order filters, where n = 100 . . . 500 or even larger, is explored in [START_REF] Ebert | Accelerated design of linear or minimum phase FIR filters with a Chebyshev magnitude response[END_REF][START_REF] Pelloni | On the design of high-order linear phase FIR filters[END_REF]. Applications include the design of transmultiplexers and filter banks for high resolution spectral analysis.

A good choice for the starting reference vector ω = (ω k ) 0 k n+1 is critical for practical convergence. The default method employed in practice [START_REF] Antoniou | Digital Signal Processing: Signals, Systems, and Filters[END_REF]Ch. 15.3.1] is to take the elements of ω uniformly from Ω ⊆ [0, π]. An empirical justification for this comes from a parallel with how one initializes the Remez algorithm on [-1, 1]. There, good choices are usually Chebyshev nodes of the second kind (see Section 4.1), which, according to Remark 2.2, map to a uniform grid on [0, π].

For a large number of high-degree problems this classic choice does not work well. Several alternatives are available. One approach is to pick reference sets according to some empirically computed distributions [START_REF] Adams | On the fast design of high-order FIR digital filters[END_REF][START_REF] Pelloni | On the design of high-order linear phase FIR filters[END_REF][START_REF] Shpak | A generalized Remez method for the design of FIR digital filters[END_REF]. Another tactic [START_REF] Psarakis | A robust initialization scheme for the Remez exchange algorithm[END_REF] is to compute a least-squares optimal filter with the same band and weight constraints as the minimax problem and take the local extrema of its approximation error as the initial reference set ω. Using Carathéodory-Fejer (CF) near-best approximations [START_REF] Pachón | Barycentric-Remez algorithms for best polynomial approximation in the Chebfun system[END_REF]Sec. 3.6] is also an option 2 . Nevertheless, the major shortcoming of these 2 http://www.chebfun.org/examples/approx/FiltersCF.html approaches is that they have been shown to work only in some cases (for example, uniform weighted filters for CF approximation).

Our solution: In Section 4 we introduce two new heuristics for initializing the exchange algorithm. They are quite general and effective in practice, often leading to significant speedups when compared to existing implementations. Moreover, Section 5.2 introduces useful analysis tools for debugging numerical problems, which are, in particular, frequent for high-degree designs.

3.2. Scalability issues: speeding up the extrema search. For most inputs, especially high degree ones, the computational bottleneck of the exchange algorithm is the search for potential extrema. Articles that concentrate on this aspect are [START_REF] Antoniou | Accelerated procedure for the design of equiripple nonrecursive digital filters[END_REF][START_REF] Antoniou | New Improved Method for the Design of Weighted-Chebyshev, Nonrecursive, Digital Filters[END_REF][START_REF] Bonzanigo | Some improvements to the design programs for equiripple FIR filters[END_REF]. The common denominator of all the approaches discussed in these references is that they use information pertaining to the first and/or second derivative of the error function E(ω) for accelerating the extrema search.

Our solution: Section 6 describes a search strategy based on Chebyshev root-finding. Compared to [START_REF] Pachón | Barycentric-Remez algorithms for best polynomial approximation in the Chebfun system[END_REF], our domain subdivision routine is not recursive, making it very easy to introduce parallelism.

Two new heuristics for choosing the initial reference

The first step of the algorithm we stated in Section 2.2 is to choose the initial reference vector ω (or equivalently x).

4.1.

A reminder on Lebesgue constants, Chebyshev nodes and Chebyshev polynomials. Consider a compact set K ⊂ R and the finite-dimensional space of weighted polynomials wP n (R) = span R {w(x), w(x)x, . . . , w(x)x n } of degree at most n, where w ∈ C(K) is a positive weight function. If T = {t 0 , . . . , t n } ⊂ K and {φ 0 , . . . , φ n } is a basis of wP n (R), we introduce the Vandermonde-like matrix

V (t 0 , . . . , t n ) = [v ij ] := [φ j (t i )].
The Lagrange interpolation operator at the nodes T has basis functions

(3) i (x) = det V (t 0 , . . . , t i-1 , x, t i+1 , . . . , t n ) det V (t 0 , . . . , t i-1 , t i , t i+1 , . . . , t n ) , i = 0, . . . , n
and is defined as

L T f (x) = n i=0 f (t i ) i (x).
Its operator norm (also called Lebesgue constant), satisfies [42, Thm. 4.3]

Λ K,T,w = L T = max x∈K n i=0 | i (x)|
and offers a quantitative way to measure the quality of a family of points for performing interpolation over K. Indeed, for every

f ∈ C(K) we have [42, Thm. 3.1] (4) f -L T f K (1 + Λ K,T,w )dist K (f, wP(R)).
For us, Lebesgue constants over the compact X play a major role in the next subsection, while in Section 5, Lebesgue constants over [-1, 1] The n-th Chebyshev polynomial of the first kind is

T n (cos(ω)) = cos(nω), ∀ω ∈ [0, π],
or recursively

T 0 (x) = 1, T 1 (x) = x, T n+2 (x) = 2xT n+1 (x) -T n (x), ∀n ∈ N.
In similar fashion, Chebyshev polynomials of the second kind are given by

U 0 (x) = 1, U 1 (x) = 2x, U n+2 (x) = 2xU n+1 (x) -U n (x), ∀n ∈ N.
Going back to formula (2), we see that the transfer function H(ω) is in fact a linear combination of Chebyshev polynomials, i.e.,

H(ω) = n k=0 h k cos(kω) = n k=0 h k T k (cos(ω)).
The Chebyshev nodes of the first kind are

µ k = cos (n -k + 1 2 )π n + 1 , k = 0, . . . , n,
and correspond to the roots of T n+1 , while the Chebyshev nodes of the second kind are the local extrema of T n over [-1, 1]. They can be given analytically as

ν k = cos (n -k)π n , k = 0, . . . , n.

4.

2. An approximate Fekete points approach. Equation ( 4) suggests that point sets with small Lebesgue constants over X are good choices for starting the exchange algorithm. Fekete points are a natural choice. In the language of Section 4.1, they maximize the Vandermonde determinant in the denominators of the Lagrange basis functions (3), resulting in a value bounded by the dimension of the interpolation space.

As noted in [START_REF] Sommariva | Approximate Fekete points for weighted polynomial interpolation[END_REF][START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF], their exact computation is very challenging in general. Despite this, recent articles [START_REF] Bos | Geometric weakly admissible meshes, discrete least squares approximations and approximate Fekete points[END_REF][START_REF] Bos | On the calculation of approximate Fekete points: the univariate case[END_REF][START_REF] Sommariva | Computing approximate Fekete points by QR factorizations of Vandermonde matrices[END_REF][START_REF] Sommariva | Approximate Fekete points for weighted polynomial interpolation[END_REF]] explore the use of a greedy algorithm for computing approximate Fekete points. It is a QR-based factorization routine which extracts large volume submatrices from rectangular Vandermonde-type matrices defined on an appropriate discretization of X. Numerical examples of weighted polynomial interpolation problems using this approach are discussed in [START_REF] Sommariva | Approximate Fekete points for weighted polynomial interpolation[END_REF]. Here, we only give a very brief statement of this algorithm when applied to our context.

For a degree n minimax approximation in the form of Problem 2.1, if

Y = {y i } ⊂ X, 0 i < N, N > n + 2, is a suitable discretization of X ⊂ [-1, 1], we construct the N × (n + 2) matrix V (y 0 , • • • , y N -1 ) = [v ij ] := [W (arccos(y i )) cos (j • arccos(y i ))], 0 i < N, 0 j n + 1
and aim at extracting an (n + 2) × (n + 2) submatrix of maximum volume. Although this problem is known to be NP-hard [START_REF] Magdon-Ismail | On selecting a maximum volume sub-matrix of a matrix and related problems[END_REF], the following linear algebra-based greedy algorithm gives good results in practice: (2) QR-based Linear System Solver: Using a column pivoting-based QR solver (via an equivalent of LAPACK's DGEQP3 routine [START_REF] Anderson | [END_REF]), find w ∈ R N , a solution to the underdetermined system b = V t w. (3) Subset Selection: Take Y as the set of elements from Y whose corresponding terms inside w are different from zero, that is, if y i ∈ Y and w i = 0, then y i ∈ Y .

Appealing choices for Y are, so-called, (weakly) admissible meshes [START_REF] Calvi | Uniform approximation by discrete least squares polynomials[END_REF]. This is due to the fact that approximate Fekete points chosen from admissible meshes are known to have the same asymptotic behavior as true Fekete points. Using this theory, we can construct very simple weakly admissible meshes for X. For instance, if

X = m k=1 X k , m ∈ N * is a finite disjoint union of closed intervals, then Y = m k=1 Y k ,
where Y k corresponds to the (n + 1)-th order Chebyshev nodes of the second kind mapped to X k , is a suitable discretization of X. Also, it has a relatively small size: O(n) elements if we consider m to be a small constant.

Despite being very robust, the only major downside of this approach is that, because of the O(n 3 ) linear solver, it does not scale particularly well to huge degrees. Hence, in practice, we suggest limiting its use up to degree n = 500 problems and combine it, if needed, with the following complementary approach.

4.3.

A reference scaling idea. The asymptotics of (weighted) minimax polynomial approximation problems is an important subject of Potential Theory [START_REF] Levin | Potential theoretic tools in polynomial and rational approximation[END_REF]. Of particular interest to us is the fact that, asymptotically, the final reference sets of minimax approximations follow a certain distribution over X, called the equilibrium distribution of X (see, for instance, [START_REF] Shen | The Potential Theory of Several Intervals and Its Applications[END_REF]Sec. 5] for the necessary theoretical results). Computing it, in all but the simplest setting (the interval), is an involved procedure [START_REF] Embree | Green's Functions for Multiply Connected Domains via Conformal Mapping[END_REF][START_REF] Shen | The Potential Theory of Several Intervals and Its Applications[END_REF], cumbersome to carry out in practice.

An important element developed in the aforementioned articles, confirmed experimentally in our tests, is that the number of reference values for minimax polynomial approximants inside each interval of X is directly proportional with the degree. As an example, we consider the design specification from Example 2.7 and look at what happens when uniform initialization is used within the exchange algorithm. The a/b entries in Table 1 show the number of reference points inside each band before the first iteration (the a value) and upon convergence (the b value), for three As expected, the b values tend to be proportional with n, while for the second and third band, the corresponding a and b become further apart as we increase the degree. This behavior allowed us to develop a relatively fast and very simple to implement heuristic for choosing an initial reference. The core of the approach is that, if, for a degree n approximation, we first compute the minimax result of degree n/2 , then we have a very good idea on the number of reference values to put inside each band for the degree n problem. The values of the new reference set x are established by taking the n/2 + 2 final references of the smaller result and adding the remaining nn/2 points uniformly between them, thus ensuring a similar distribution of the references inside X for the two degrees. If needed, this strategy can be applied recursively.

Since it is based on asymptotic results, this idea works best for high degree designs, where the AFP-based strategy becomes expensive. 2 shows some results when considering Examples 2.6 to 2.9. For each entry of the Iterations columns, the first value corresponds to the uniform initialization technique, while the second and third ones represent our scaling and AFP approaches. We can observe that, in many cases, there is a considerable reduction in the number of iterations required for convergence. This frequently translates to smaller execution times as well, despite a more involved setup of the starting references. Still, the biggest advantage of these heuristics is that often, they allow us to address, in a simple manner, filter design problems where uniform initialization gives rise to numerical instability issues and the implementation does not converge (NC). We also tried the strategies from [START_REF] Psarakis | A robust initialization scheme for the Remez exchange algorithm[END_REF][START_REF] Shpak | A generalized Remez method for the design of FIR digital filters[END_REF], but encountered numerical problems on the larger designs of Examples 2.8 and 2.9. For [START_REF] Psarakis | A robust initialization scheme for the Remez exchange algorithm[END_REF], the quadratic linear solver used to determine the least-squares optimal filter did not work accurately, while for the heuristic of [START_REF] Shpak | A generalized Remez method for the design of FIR digital filters[END_REF], the starting reference did not contain a good proportion of values inside each band for the algorithm to converge numerically.

In dealing with Example 2.8 we applied reference scaling two times, starting from a degree n = 130 filter and then moving up to n = 260 and finally n = 520. We only needed 4, 3 and 3 iterations for convergence. AFP initialization performed even better, only requiring one iteration. For Example 2.9, we again applied reference scaling two times, for a number of 6, 3 and 3 iterations, whereas AFP proved to be too costly.

Note that, actually, convergence with uniform initialization for Examples 2.7 and 2.8 (the two * -ed cases) occurs in spite of numerical problems: this behavior is fortunate and does not generally happen in such a bad numerical context. See also Remark 5.1 for more insight.

Inside the test files, we have provided code for the design of over 70 different filters that are considered high degree (with n in the order of hundreds and thousands). While the reduction in number of iterations, when computable, was very variable (savings from 0% to over 70% compared to uniform initialization were observable on most of the examples), reference scaling always ensured convergence, whereas for more than 30% of the considered filters, uniform initialization failed to converge. The AFP-based idea also performed very well, with only a small percentage (around 5%) of the tests failing. The examples where this happened had degree n 1000 and were much better suited for the scaling approach.

Computing the current interpolation function

Next, we need to determine the value of the current alternating error δ and derive a formula for computing H(ω) (Step 2). If we apply the Alternation Theorem on the current reference set, our unknowns satisfy the following equations:

E(ω i ) = W (ω i ) [D(ω i ) -H(ω i )] = (-1) i δ x ,
i = 0, . . . , n + 1, where D denotes the ideal filter we are approximating, W is the corresponding weight function for the error and δ x = |δ x |. By taking H(ω) in its cosine expansion from equation (2), we get the following linear system of equations in (h k ) 0 k n and δ x :

          1 cos(ω 0 ) cos(2ω 0 ) • • • cos(nω 0 ) 1 W (ω 0 ) 1 cos(ω 1 ) cos(2ω 1 ) • • • cos(nω 1 ) -1 W (ω 1 ) . . . . . . . . . . . . . . . 1 cos(ω n+1 ) cos(2ω n+1 ) • • • cos(nω n+1 ) (-1) n+1 W (ω n+1 )                  h 0 h 1 . . . h n δ x        =        D(ω 0 ) D(ω 1 ) . . . D(ω n ) D(ω n+1 )       
This system has a unique solution, since the (n + 2) × (n + 2) matrix is non-singular (see [START_REF] Veidinger | On the numerical determination of the best approximation in the Chebyshev sense[END_REF] for a proof idea). According to [START_REF] Antoniou | Digital Signal Processing: Signals, Systems, and Filters[END_REF]Ch. 15.3.3] and [START_REF] Oppenheim | Discrete-Time Signal Processing[END_REF]Ch. 7.7.3], solving it is not encouraged, since it is computationally inefficient and susceptible to numerical ill-conditioning.

In practice, one deduces δ x analytically and uses barycentric Lagrange interpolation to determine H(ω) (i.e. the original choice of [START_REF] Parks | Chebyshev Approximation for Nonrecursive Digital Filters with Linear Phase[END_REF]). The entirety of this section focuses on the numerical stability of this approach and of computing δ x . 5.1. Barycentric Lagrange interpolation. In recent years, barycentric interpolation has become an active research topic in Numerical Analysis. This is mostly due to the nice numerical properties of this interpolation scheme [START_REF] Berrut | Barycentric Lagrange Interpolation[END_REF][START_REF] Floater | Barycentric rational interpolation with no poles and high rates of approximation[END_REF].

The basic setting is the following: if f : [x -, x + ] → R and n ∈ N, let x ∈ R n+2 be a vector of distinct interpolation points x k ∈ [x -, x + ], k = 0, . . . , n + 1, given in increasing order, together with y ∈ R n+2 , where y k = f (x k ), k = 0, . . . , n + 1. We want to find a polynomial p with real coefficients of degree at most n + 1 which interpolates f at x (i.e., p(x k ) = y k , k = 0, . . . , n + 1).

According to [START_REF] Rutishauser | Vorlesungen über Numerische Mathematik[END_REF], the barycentric forms of p are given by ( 5)

p(x) = (x) n+1 k=0 w k x -x k y k ,
known as the first barycentric interpolation formula, and

(x) = n+1 k=0 w k x -x k y k n+1 k=0 w k x -x k , (6) p
the second (or proper ) barycentric formula, where (x) = n+1 i=0 (xx i ) and ( 7)

w k = 1 (x k ) = 1 i =k (x k -x i ) , k = 0, . . . , n + 1.
If the barycentric weights w k have been precomputed, both ( 5) and ( 6) can be evaluated with only O(n) arithmetic operations at an arbitrary point x.

To use barycentric interpolation for the evaluation of H(ω), we first compute δ x using the formula ( 8)

δ x = n+1 k=0 w k D(ω k ) n+1 k=0 (-1) k w k W (ω k ) (see [39, Sec. 3.3] for a proof idea).
We have ( 9)

H(ω) = n+1 k=0 w k x -x k c k n+1 k=0 w k x -x k where x = cos(ω), c k = D(ω k ) -(-1) k δ x W (ω k )
and w ∈ R n+2 is the weight vector whose elements are computed according to [START_REF] Antoniou | Digital Signal Processing: Signals, Systems, and Filters[END_REF]. Such a second barycentric formula is usually preferred for computing the current frequency response. We will elaborate on why this is a good idea in Section 5.2.2.

5.2. Numerical stability issues. The evaluation of ( 8) and ( 9) is split into the computation of the barycentric weights (done once for each new reference vector x) and of the sums in the numerators and the denominators. For a numerically robust evaluation of the w k 's (O(n

= a k (1 + ε k ), k = 0, . . . , n -1 of the input terms, then the relative error E n of computing ŝn = n-1 k=0 âk is upper bounded by max k |ε k | C n . Indeed, E n = |ŝ n -s n | |s n | = n-1 k=0 ε k a k n-1 k=0 a k n-1 k=0 |ε k | |a k | n-1 k=0 a k max k |ε k | n-1 k=0 |a k | n-1 k=0 a k .
Even if the perturbations are as small as possible with the current arithmetic, meaning ε k = O(u), k = 0, . . . , n -1, if there is severe cancellation in computing s n and C n u -1 , then E n will most likely be greater than 1 and give totally inaccurate results. Since p x and q x are obtained using finite precision arithmetic, in certain situations, there is pronounced numerical cancellation when computing p x . This can happen when:

(1) the minimax error at the end of executing the exchange algorithm is very small (in the sense of being relatively close to or smaller than u); (2) the final minimax error is much larger than u, but the starting leveled error δ x is in the order of u or much smaller. In practice, the second bullet point is more likely to happen than the first one, since filters with minimax error δ X < 10 -10 ( u ∼ 10 -16 ) seldom seem to be required.

When using uniform initialization, such issues are not hard to find, as they sometimes occur in the starting iteration(s) of the exchange method. Consider, for instance, the specification from Example 2.8. Although the final alternating error is around 1.6067 • 10 -7 , the first iteration should give δ x 1.521 • 10 -21 u and a condition number for p x of about 6.572 • 10 20 72967 • u -1 . Because of this extreme ill-conditioning, a double-precision computation of δ x gives a wrong value of 5.753 • 10 -19 .

Based on our discussion up to this point, we propose two possible solutions:

(1) use a higher precision arithmetic for doing all computations in the current iteration (barycentric weights, δ x and local extrema of E(ω)); (2) pick a better starting reference (see initialization strategies from Section 4). The first alternative is generally the safest, but costlier choice. From our experience, around 150-200 bits of precision are more than enough for the problems discussed here. Hence, the MPFR-based version of our code uses 165 bits of precision by default. Once enough iterations have passed and the current condition number for p x is sufficiently small (i.e. smaller than say 10 8 ), there is a good chance that using double-precision arithmetic for the remaining iterations will still lead to an acceptable result.

For the second strategy, the basic intuition is that the starting δ x should be much larger and closer to the minimax one. This consequently implies that the condition number for p x is much smaller compared to u -1 and the summations are more stable. To illustrate this, we again consider Example 2.8, but this time, use the reference scaling strategy from Section 4.3. The initial leveled error is δ x 1.456 • 10 -7 u -1 and the condition number for p x is around 6.864 • 10 6 . Double-precision computations are again quite viable. Since it usually results in much smaller execution times, we recommend starting with this second approach.

Numerical problems are much less likely to occur when computing q x . Since we assume the interpolation nodes x are ordered by value, the barycentric weights w i alternate in sign. This means that all the terms of q x have the same sign and hence, its condition number has minimum value one.

5.2.2.

Lebesgue constants as tools for detecting numerical problems. The numerical stability of ( 5) and ( 6) over [-1, 1] was first looked at in a rigorous manner in [START_REF] Rack | The numerical stability of evaluation schemes for polynomials based on the lagrange interpolation form[END_REF] and later refined in [START_REF] Higham | The numerical stability of barycentric Lagrange interpolation[END_REF]. This second reference shows how the modified Lagrange formula ( 5) is backward stable in general, whereas ( 6) is shown to be forward stable for vectors of points x ∈ R n+2 having a small Lebesgue constant 3 . Despite a less favorable numerical behavior of (6) in general, for sets of points having a small Lebesgue constant, like Chebyshev nodes, using the second barycentric formula is preferable [START_REF] Berrut | Barycentric Lagrange Interpolation[END_REF][START_REF] Mascarenhas | The stability of barycentric interpolation at the Chebyshev points of the second kind[END_REF]. More recently, [START_REF] Mascarenhas | On the backward stability of the second barycentric formula for interpolation[END_REF] argues that ( 6) is backward stable if the relevant Lebesgue constant associated with the interpolation vector x is small.

An example of a reference set with small Lebesgue constant is showcased in Figures 3 and4 for the specification given in Example 2.8. Both plots show the relative errors in computed transfer function H(ω) during the first iteration of the exchange algorithm. The 'exact' computations were performed with MPFR using 600 bits of precision (roughly 180 digits). In Figure 3,1],x is quite large and we see that the loss in accuracy is clearly visible around x = -1 (i.e. ω = π). The results are much better in Figure 4 for a smaller Lebesgue constant, reinforcing our earlier statements.

We make the following experimental remark: in general, the Lebesgue constants associated with the references that appear during the execution of the Parks-McClellan algorithm tend to decrease in value, with the final Λ [-1,1],x being not too large.

Typical scenarios of how Λ [-1,1],x changes during the execution of the exchange algorithm are given in Tables 3 to 5 for Examples 2.6 to 2.8. The numerical estimates were computed using the Chebfun routine lebesgue mentioned in [START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF]Ch. 15]. Because of the observed behavior, computations will tend to be more stable during latter iterations. The effect of using the strategies from Section 4 is also clearly visible.

Since a large Lebesgue constant Λ [-1,1],x is symptomatic for a very small δ x , the solutions for countering illconditioning of the previous subsection also apply here.

Remark 5.1. When numerical problems of the types discussed in this section occur (typically because of uniform initialization), our exchange algorithm implementation using double-precision arithmetic can sometimes still find the correct result. What happens in such cases is that, even though our computations of the current interpolant and its approximation error are inaccurate, the next reference we retrieve from them, despite being wrong, can have a small Lebesgue constant. Having such a reference, as we saw, usually improves chances for convergence. It can be viewed as a kind of reinitialization midway through execution. Nevertheless, as our test cases show, convergence in the presence of numerical issues is a completely random behavior. It is not to be expected and, most of all, trusted. One should definitely consider using a higher working precision and/or one of the initialization techniques from Section 4.

= 100 n = 520 Iteration Λ [-1,1],x Iteration Λ [-1,1],x Iteration Λ [-1,1],x 1 4

Extrema search

The final two steps of an exchange algorithm iteration consist of computing the new reference set ω ∈ R n+2 (and equivalently x = (cos(ω k )) 0 k n+1 ). They require finding the local extrema of the current error function E(ω). Chebyshev root-finding is reviewed in the next subsection, while Section 6.2 introduces our non-recursive domain subdivision strategy. For completeness, Sections 6.3 and 6.4 talk about updating the reference set at each iteration and how to retrieve the filter coefficients upon convergence.

6.1. Chebyshev-proxy root-finding. We first consider the slightly different problem of determining the zeros of a function f ∈ C ([x -, x + ]), located inside [x -, x + ]. For simplicity, we again take [x -, x + ] = [-1, 1], and note that by suitable changes of variable, the following results hold for any closed interval. The idea of the Chebyshev-proxy root-finder (CPR) method [START_REF] Boyd | Finding the Zeros of a Univariate Equation: Proxy Rootfinders, Chebyshev Interpolation, and the Companion Matrix[END_REF][START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF] is to replace f (x) by a degree m polynomial proxy p m (x). If the chosen polynomial is an accurate enough approximation of f , then its zeros will very closely match those of f . It interpolates f at the Chebyshev nodes of the second kind ν k = cos((mk)π/m), 0 k m, and is expressed using the basis of Chebyshev polynomials of the first kind. We get The a i coefficients can be evaluated in a numerically stable way with O(m 2 ) operations, by using Clenshaw's recurrence relations [START_REF] Press | Numerical Recipes: The Art of Scientific Computing[END_REF]Ch. 5.4], or faster, in only O(m log m) operations, by means of the Discrete Cosine Transform [START_REF] Strang | The discrete cosine transform[END_REF]. The roots of p m are then computed as the eigenvalues of a generalized companion matrix [START_REF] Boyd | Computing the zeros, maxima and inflection points of Chebyshev, Legendre and Fourier series: solving transcendental equations by spectral interpolation and polynomial rootfinding[END_REF][START_REF] Day | Roots of Polynomials Expressed in Terms of Orthogonal Polynomials[END_REF], also known as a colleague matrix [START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF]Ch. 18]. This approach behaves well in practice [START_REF] Boyd | Numerical experiments on the accuracy of the Chebyshev-Frobenius companion matrix method for finding the zeros of a truncated series of Chebyshev polynomials[END_REF].

p m (x) = m k=0 a k T k (x), x ∈ [-1, 1], a k ∈ R, k = 0, . . . , m, with p m (ν k ) = f (ν k ), k = 0, . . . , m and (10)
a k = 1 m f (ν 0 )T k (ν 0 ) + 2 m m-1 i=1 f (ν i )T k (ν i ) + 1 m f (ν m )T k (ν m ), k = 0, . . . , m.
To obtain good approximations of the local extrema of E(ω), suppose we already have an accurate degree m proxy

E m (ω) = m k=0 a k T k (x), x = cos(ω).
The value of m should be thought of independently from the degree n of the target minimax response H(ω). We look for the roots of e m (ω) = dE m (ω) dx , whose Chebyshev expansion is

e m (ω) = m-1 k=0 b k T k (x)
e m (ω) = m-1 k=0 c k U k (x),
where c k-1 = ka k , k = 1, . . . , m. Although the CPR method is introduced in the context of T m , we found the second formula for e m more natural to use in our setting.

To find the eigenvalues of a m × m colleague matrix, for both types of Chebyshev polynomials, we can use a QR/QZ algorithm (O(m 3 ) operation count). For numerical reasons, we use the colleague matrix form suggested in [START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF]Ex. 18.3].

This cost can be reduced to O(m 2 ) computations by dividing the initial domain [x -, x + ] into several subintervals and taking Chebyshev interpolation polynomials of smaller degree on each of them. The roots can then be determined as the collection of zeros inside all subintervals. More precisely: if N max < m is the maximum degree on which one is willing to use an eigenvalue algorithm, recursively split [x -, x + ] until there are Chebyshev polynomial interpolants of degree at most N max on each subinterval, which approximate f accurately.

There are several general criteria that quantify when we have a good enough approximation of f (x) [START_REF] Boyd | Computing the zeros, maxima and inflection points of Chebyshev, Legendre and Fourier series: solving transcendental equations by spectral interpolation and polynomial rootfinding[END_REF][START_REF] Boyd | Finding the Zeros of a Univariate Equation: Proxy Rootfinders, Chebyshev Interpolation, and the Companion Matrix[END_REF]. The [START_REF] Pachón | Barycentric-Remez algorithms for best polynomial approximation in the Chebfun system[END_REF] Remez algorithm implementation uses such a recursive approach [START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF]Ch. 18] by means of the Chebfun roots command. Note that the intermediary interpolants computed before the final stage subdivisions are not used for the actual root finding.

6.2.

A new subdivision strategy. We use a different, non recursive approach. In general, during each iteration of the Parks-McClellan algorithm, we have estimates about the location and number of extrema of E; if we take a subinterval defined by two consecutive points of the current reference vector x located inside the same band, say [x i , x i+1 ], then we usually expect to have between zero and two potential extrema of E inside it. Figure 5 summarizes this scenario. We then use a degree N max Chebyshev proxy on each subinterval. Values of N max we frequently used in our experiments are 4 and 8. Dynamic range issues [START_REF] Boyd | Finding the Zeros of a Univariate Equation: Proxy Rootfinders, Chebyshev Interpolation, and the Companion Matrix[END_REF]Sec. 7.4] are generally absent because the initial approximations obtained with the approaches from Section 4 are of good quality. Since there is only one subdivision level (no intermediary interpolants), there will be savings in computation time. At each iteration, we use O(n) subintervals of the form [x i , x i+i ]. To cover all of X, we also consider all intervals of the form [a, x i ] ⊆ X and [x j , b] ⊆ X, where a and b are band extremities of X and x i and x j are elements of x which are closest to a and b, respectively. Inside each subinterval, we evaluate the error N max + 1 times in order to construct a suitable proxy for E and its derivative. Because of barycentric interpolation, each evaluation of E will usually require only O(n) operations, while solving each small eigenvalue problem amounts to a constant number of computations. The total operation count will be O(n 2 ).

In addition to requiring the evaluation of E a reasonable number of times (in the range of 4n to 8n), this approach can also be very easily parallelized (the computations on each subinterval are data parallel ). For large degree approximations on multicore systems, the speedup is considerable. We give some examples in Section 7.2. Remark 6.1. This quadratic extrema search scheme gives very reasonable execution times in practice, even for degrees in the order of thousands, as we will see in Section 7.2. Nevertheless, let us mention that potential speedups are still possible; because of apriori knowledge on the exact number of evaluations of E involving formula [START_REF] Bonzanigo | Some improvements to the design programs for equiripple FIR filters[END_REF], it is feasible to use fast multipole methods for evaluating polynomial interpolants [START_REF] Dutt | Fast algorithms for polynomial interpolation, integration, and differentiation[END_REF] at n arbitrary points. If ε is the numerical precision for our computations, such methods need O(n log(1/ε)) operations at each iteration. Parallelization is possible here as well; the O(N max n) polynomial evaluations can be split in O(N max ) independent applications of the fast multipole method on batches of n points. 6.3. The new reference set. Each time we compute a valid eigenvalue y, we add it to the set of potential extrema x only if |E(arccos(y))| δ. We also add to x the extremities of each subinterval if the error they yield is larger or equal in absolute value than the current minimax error δ. This allows us to successfully treat cases where band edges of X are in the final reference set, which is always true for type I stopband and passband filters (see [START_REF] Oppenheim | Discrete-Time Signal Processing[END_REF]Ch. 7.7.1] for a proof).

Assuming the elements of x are ordered, for each subset of consecutive values where the error has the same sign, we only keep one element with largest absolute error. The new reference x is then constructed by taking n + 2 elements of x where the error alternates in sign and has the largest absolute values. A more detailed presentation on how to do this is available in [START_REF] Ahsan | Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design[END_REF]Sec. 4.1]. This strategy is the most robust one we tested. For other discussions, see [START_REF] Antoniou | Digital Signal Processing: Signals, Systems, and Filters[END_REF]Ch. 15.3.4] and [39, Sec. 2.2]. 6.4. Retrieving the filter coefficients. Upon convergence, we compute the coefficients (h k ) 0 k n of the final frequency response H(ω) (eq. ( 2)). As we already saw, we can do this in a numerically stable way with a quadratic number of computations by using formula [START_REF] Bos | Geometric weakly admissible meshes, discrete least squares approximations and approximate Fekete points[END_REF] and Clenshaw's algorithm or in O(n log n) operations with a DCT-based interpolation scheme. Similar to Remark 6.1, the total cost is dictated by how we evaluate the barycentric form of H(ω) at the n-th order Chebyshev nodes: O(n 2 ) operations, since we use formula [START_REF] Bonzanigo | Some improvements to the design programs for equiripple FIR filters[END_REF] at each node. With the evaluation scheme of [START_REF] Dutt | Fast algorithms for polynomial interpolation, integration, and differentiation[END_REF], we would require O(n log(n/ε)) operations (fast polynomial evaluation + DCT-based Chebyshev interpolation).

Figure 6 shows a small sample of the minimax approximation error E, for Example 2.9, after the filter coefficients have been computed using Clenshaw's recurrence relations. Although the degree n = 53248 is quite large, the computations are numerically accurate, with the equioscillations mentioned inside Theorem 2.3 clearly visible.

Implementation details

In this section we give practical information on how someone can use our routines and compare them to equivalent ones available in widely used signal processing packages.

To summarize the previous three sections, each iteration of the exchange algorithm takes O(n 2 ) operations with our code (see https://github.com/sfilip/firpm). Choosing the starting reference takes O(n) operations if uniform initialization is used, but this cost can jump to O(n 2 ) with the approach from Section 4.3 or O(n 3 ) operations for Section 4.2. Since we limit the use of AFP-based initialization to moderate degrees, the cubic cost is not problematic in practice. As discussed in Section 4.4, these more expensive strategies are most of the time worthwhile and sometimes even necessary. Also, the computational bottleneck of our implementation (the extrema search of Section 6.2) is embarrassingly parallel. Exploiting it is very effective in practice (see Table 7).

7.1. User interface. Our code requires a small number of external libraries in order to work. Besides calling MPFR inside the multiple precision version, we also use the Eigen library [27] to perform all the eigenvalue computations when searching for a new reference set. Because it is designed with template metaprogramming techniques, the code calling Eigen requires little to no changes between the different versions. To parallelize the extrema search in Section 6.2, we use OpenMP [START_REF] Dagum | OpenMP: an industry standard API for shared-memory programming[END_REF].

All three versions of our code have an almost identical interface, inspired by the style used for the MATLAB implementation of the Parks-McClellan algorithm. As such, we will focus on describing only the double-precision version.

The function for designing type I and II filters has the following prototype:

PMOutput firpm(std::size_t N, std::vector<double>const& f, std::vector<double>const& a, std::vector<double>const& w, double epsT = 0.01, int Nmax = 4); where • N+1 is the number of coefficients of the designed filter;

• f is the vector of the frequency band edges of Ω ⊆ [0, π], normalized to [0, 1], and given in increasing order;

• a contains the desired amplitudes at each of the points from f; • w is the vector of weights on each band specified by f (its size is half that of f and a);

• epsT convergence parameter threshold from Step 5 in Algorithm 1 from Section 2.2. The default value is taken from [START_REF] Antoniou | Digital Signal Processing: Signals, Systems, and Filters[END_REF]Ch. 15.3]. • Nmax designates the degree of the Chebyshev interpolants used for the extrema search from Section 6.2, with a default value of 4. For instance, designing a degree n = 100, type I filter adhering to the specification given in Example 2.6, results in the following, valid C++11, function call:

PMOutput res = firpm(200, {0, 0.4, 0.5, 1}, {1, 1, 0, 0}, {1, 1});

The returned PMOutput object res has several member variables that can be useful to a filter designer. The vector res.h corresponds to the final coefficients of the filter transfer function from equation (1), while res.x is the final reference vector, belonging to X. The number of iterations required for convergence are stored in the variable res.iter, while the value of the final reference error δ X is denoted with res.delta.

The firpm functions use uniform initialization by default. To use reference scaling or approximate Fekete points, we supply the functions firpmRS and firpmAFP. The firpmRS versions allow the user to choose between uniform and AFP-based initialization at the lowest level.

Timings.

There are multiple implementations of the Parks-McClellan algorithm available, so we compared our approach to those we believe are the most widely used and/or robust in practice. Some results are given in Table 6. Because they are written in different languages, there are bound to be some differences in terms of execution times. To make matters as fair as possible, we disabled parallelization of the extrema search in our code and went for uniform initialization in all the test cases. Default grid size parameters (see the discussion from the beginning of Section 3) were used for the routines in the last four columns and a default value of N max = 4 for our implementation.

The two optimized routines from [START_REF] Ahsan | Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design[END_REF] are at heart efficient rewritings of the original code of [START_REF] Mcclellan | A computer program for designing optimum FIR linear phase digital filters[END_REF], the only difference between the two being how the reference set gets updated at each iteration. Together with our implementation, they were the only ones that were able to converge on 7 out of the 8 test cases. Even so, our code is the most robust one in Table 6. Runtime (real time) comparisons of our IEEE-754 double-precision implementation of the Parks-McClellan algorithm with those available in GNURadio 3.7.8.1 (also written in C++), MATLAB R2014b, the one from SciPy 0.16.1 (python code) and two other MATLAB implementations. The same machine, a quad core 3.6 GHz 64-bit Intel Xeon(R) E5-1620 running Linux 4.1.12 (15.9), was used for all the tests. The average execution times (in seconds) of running each piece of code 50 times, are given. When a routine did not converge to the minimax result, NC was used in place of the execution time. The a/b entries in the last column correspond to the two implementations described in [START_REF] Ahsan | Two Novel Implementations of the Remez Multiple Exchange Algorithm for Optimum FIR Filter Design[END_REF]. terms of execution time. We also tested the corresponding routine from Scilab 5.5.1 (written in Fortran), but did not achieve convergence for any of the test cases. The effects of starting with a better reference (via scaling or approximate Fekete points) and parallelization of the extrema search are showcased in Table 7. As we already emphasized, the greatest gain is for the last two examples, where the filter degrees are larger. In the case of Example 2.9, for which only our routine converged, enabling parallelization on the quad core machine we used for testing, resulted in our code running 3.3 times faster than the purely sequential version. 7.3. Other applications. Despite the focus on filter design problems, the routines we provide are more general. They can be used for arbitrary weighted minimax polynomial approximation problems involving multi-interval domains. Such a scenario arises, for instance, when minimizing the relative error of approximation for a function with zero(s).

Example

As an example, consider the function f (x) = ln(x)e sin(x) . We want to find the polynomial p of degree at most n = 10, which best approximates f over [1 -2 -3 , 1 + 2 -3 ] in terms of relative error (corresponding weight function w(x) = 1/ |f (x)|). Since f has a zero at x = 1, in a double-precision context, we can split the domain into two parts: that is, remove all the values between the largest double-precision value smaller than 1 and the smallest double-precision value greater than 1.

Conclusion

In this article we have presented several ideas that aid in the development of weighted minimax polynomial approximations, with a central focus on optimal linear-phase FIR filters. These contributions amount to the following:

• introduce two new initialization strategies (Sections 4.2 and 4.3); they ensure the convergence of our routine in cases where all other implementations fail to work, while frequently incurring significant speedups as well (consequence of a small number of iterations required for convergence); • present pertinent analysis tools which help diagnose when the exchange algorithm is prone to numerical problems in practice (Section 5.2); • introduce a variation of a well-established root-finding approach [START_REF] Boyd | Finding the Zeros of a Univariate Equation: Proxy Rootfinders, Chebyshev Interpolation, and the Companion Matrix[END_REF][START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF] which allows one to design FIR filters in an efficient way (Section 6.2). Equally noteworthy is the fact that this study has helped us develop an efficient and highly parallel software library. As we saw throughout all of the examples we considered in this text, our routine outperforms other existing codes in terms of scalability and numerical accuracy.

Since at the end of the filter synthesis framework, the final coefficients are quantized, the result of the Parks-McClellan algorithm needs to be adequately modified. We are currently working on efficient ways of solving this problem. We also hope to investigate if similar ideas to the ones we presented here can be applied to other digital signal processing tasks, like those requiring FIR filters with complex coefficients [START_REF] Karam | Complex Chebyshev approximation for FIR filter design[END_REF][START_REF] Karam | Design of optimal digital FIR filters with arbitrary magnitude and phase responses[END_REF] or IIR filters [START_REF] Saramäki | An efficient Remez-type algorithm for the design of optimum IIR filters with arbitrary partially constrained specifications[END_REF][START_REF] Saramäki | Generalizations of classical recursive digital filters and their design with the aid of a Remez-type algorithm[END_REF][START_REF] Selesnick | New exchange rules for IIR filter design[END_REF][START_REF] Selesnick | Magnitude squared design of recursive filters with the Chebyshev norm using a constrained rational Remez algorithm[END_REF].

2. 1 .

 1 Optimal FIR filters and weighted polynomial approximation. Type I filters are characterized by n 1 = n 2 = n in (1), where b -k = b k , k = 1, . . . , n. We then have the equivalent form (2) H(ω) = n k=0 h k cos(ωk), where h 0 = b 0 and h k = 2b k , 1 k n. We address the following: Problem 2.1 (Equiripple (or minimax) FIR filter design). Let Ω be a compact subset of [0, π] and D(ω) an ideal frequency response, continuous on Ω. For a given filter degree n ∈ N, we want to determine H(ω) = n k=0 h k cos(ωk) such that the weighted error function E(ω) = W (ω) (D(ω) -H(ω)) has minimum uniform norm E(ω) ∞,Ω = sup ω∈Ω |E(ω)| , where the weight function W is continuous and strictly positive over Ω. Remark 2.2. By the change of variable x = cos(ω), H(ω) is a polynomial [38, Ch. 7.7] in x of degree at most n. Problem 2.1 is thus equivalent to weighted minimax polynomial approximation over a compact subset of [-1, 1].

( 3 )

 3 Extrema Search: Determine the current error function E(ω) and find its local extrema over Ω, where |E(ω)| δ ω . (4) Reference Set Update: Retrieve a new reference vector ω = (ω k ) 0 k n+1 from the set of potential extremas found at Step 3 by picking n + 2 that include the global extrema of E(ω) and for which the error alternates in sign.

  return the filter characterized by the set of equioscillating frequencies ω . If not, go to Step 2, with ω = ω as the new reference vector.

Figure 1 .

 1 Figure 1. Tolerance scheme and ideal frequency response for a lowpass filter.

Figure 2 .

 2 Figure 2. Design example showcasing a transition band anomaly and how it can be removed.

  with weight function w(x) = 1 are used to justify the numerical stability of barycentric Lagrange interpolation inside the exchange algorithm. Classic examples of excellent interpolation points with almost minimal Lebesgue constants over [-1, 1] are the Chebyshev nodes; they grow as 2 π log n + O(1) when n tends to infinity [58, Thm. 15.2]. To define them, we need Chebyshev polynomials.

Algorithm 2 :

 2 AFP (Approximate Fekete Points) routine Input: appropriate discretized subset Y = {y 0 , . . . , y N -1 } ⊂ X, N > n + 2 Output: Y = y 0 , . . . , y n+1 ⊂ Y s.t. the (n + 2) × (n + 2) Vandermonde-like submatrix constructed from the elements of Y has a large volume (1) Initialization: V = V (y 0 , . . . , y N -1 ), b ∈ R n+2 , b = (1, . . . , 1) t

4. 4 .

 4 Numerical examples. Performance-wise, Table

Figure 3 .

 3 Figure 3. Relative errors (correct significant digits) when computing the starting transfer function for Example 2.8 using both types of barycentric formulas. Uniform initialization (Λ [-1,1],x 4.24973 • 10 14 ) is used.

Figure 4 .

 4 Figure 4. Relative errors (correct significant digits) when computing the starting transfer functions for Example 2.8 using both types of barycentric formulas. Reference scaling (Λ [-1,1],x 1.29773 • 10 5 ) is used.

Figure 5 .

 5 Figure 5. Our interval subdivision strategy for computing the extrema of E(ω).

  , where: b k-1 = 2ka k + b k+1 , k = 1, . . . , m, and b m = b m+1 = 0. Similarly, if we take into account that dT m dx = mU m-1 , for m 1, we have

Figure 6 .

 6 Figure 6. The value of the approximation error E on a small interval [0.247168, 0.24869] of X, for the minimax filter satisfying the specification for Example 2.9.

[ 1 -

 1 2 -3 , 0.99999999999999988897769753748434595763683319091796875] and [1.0000000000000002220446049250313080847263336181640625, 1 + 2 -3 ],

Figure 7 .

 7 Figure 7. Degree n = 10 minimax approximation of f (x) = ln(x)e sin(x) in terms of relative error.

Figure 7

 7 Figure 7 shows the minimax error (f (x)p (x))/ |f (x)| obtained by running our implementation of the exchange algorithm on this new two interval domain, together with the n + 2 = 12 equioscillation points around the extremal values ±6.32824 • 10 -13 .

  Input: the filter degree n, the ideal frequency response D, the weight function W , the frequency set of interest Ω, convergence parameter threshold ε t Output: the (h k ) 0 k n coefficients of the final frequency response H(1) Initialization: Pick a reference vector of frequencies ω = (ω k ) 0 k n+1 from Ω.(2) Finite Set Approximation & Interpolation: Find the current frequency response H(ω) and its associated alternating error δ ω , which correspond to solving Problem 2.1 on the elements of ω.

17, Ch. 3.4]. 2.2. The Parks-McClellan algorithm. Together with a result due to de La Vallée Poussin [17, Ch. 3.4], Theorem 2.3 has been used by Remez [47] to derive iterative procedures which converge to the minimax result [17, Ch. 3.8] [42, Ch. 8-10]. Parks and McClellan [40] were the first to give a successful Remez algorithm implementation for solving Problem 2.1: Algorithm 1: Parks-McClellan (Remez) exchange procedure

  start with two examples used in[START_REF] Psarakis | A robust initialization scheme for the Remez exchange algorithm[END_REF].

	Example 2.6. A unit weight type I lowpass filter with [0, 0.4π] passband and [0.5π, π] stopband. Multiple degrees n
	are considered.
	Example 2.7. A type I, uniformly weighted, bandstop filter with passbands [0, 0.2π], [0.6π, π] and stopband
	[0.3π, 0.5π]. Several degrees n are used.

Table 1 .

 1 Number of reference values in each band for the bandstop specification of Example 2.7 before the first iteration and at the end of execution.

		n = 50 n = 100 n = 200
	[0, 0.2π]	14/13	26/26	51/51
	[0.3π, 0.5π] 14/15	26/31	51/59
	[0.6π, π]	24/24	50/45	100/92

Table 2 .

 2 Iteration count comparison for uniform/reference scaling/ AFP-based initialization

	Example 2.6	Example 2.7	Example 2.8	Example 2.9
	Degree Iterations Degree Iterations Degree Iterations Degree Iterations
	50	11/4/6	50	14/14/4
	80	8/3/4	80	13/3/12	520	12 * /3/1	53248 NC/3/NC
	100	9/8/3	100	23 * /18/16
	different degrees.

  2 ) operations), we use the formulas and ideas of[START_REF] Berrut | Barycentric Lagrange Interpolation[END_REF] Sec. 7] and[START_REF] Pachón | Barycentric-Remez algorithms for best polynomial approximation in the Chebfun system[END_REF] Sec. 3.4]. Notice that inside[START_REF] Bonzanigo | Some improvements to the design programs for equiripple FIR filters[END_REF] we are using n + 2 points to interpolate a polynomial of degree at most n. This is intentional. We could have easily picked only a subset of n + 1 elements of x, but leaving out one element can pose numerical problems, especially if that point is the smallest or largest one from x[START_REF] Webb | Stability of Barycentric Interpolation Formulas for Extrapolation[END_REF].

	only have access to perturbed versions âk	n-1 k=0 |a k | /	n-1 k=0 a k . If we

5.2.1.

The current leveled error. Let p x and q x be the numerator and denominator values of δ x in formula

[START_REF] Berrut | Barycentric Lagrange Interpolation[END_REF]

. Unless otherwise stated, all elementary operations are performed in round-to-nearest double-precision arithmetic with unit roundoff u = 2 -53 .

Central to our analysis is the summation condition number [37, Ch. 6.1]: given n real values a 0 , . . . , a n-1 whose sum s n = n-1 k=0 a k we want to compute, its corresponding condition number is C n =

Table 3 .

 3 Lebesgue constant evolution during the execution of the exchange algorithm (i.e., first, middle and last iteration). The starting reference set is computed using uniform intialization.

	Example 2.6	Example 2.7	Example 2.8
	n = 100	n = 100	n = 520
	Iteration	Λ [-1,1],x	Iteration	Λ [-1,1],x	Iteration	Λ [-1,1],x
	1 5 9	3.72896 • 10 10 8.31678 • 10 6 4.00568 • 10 6	1 12 23	3.00467 • 10 15 7.45988 • 10 7 1.35889 • 10 7	1 6 12	4.24973 • 10 14 6.11242 • 10 5 1.38665 • 10 5

Table 4 .

 4 Lebesgue constant evolution when the reference scaling approach described in Section 4.3 is used.

	Example 2.6	Example 2.7	Example 2.8
	n = 100	n = 100	n = 520
	Iteration	Λ [-1,1],x	Iteration	Λ [-1,1],x	Iteration	Λ [-1,1],x
	1 4 8	5.98098 • 10 6 4.68392 • 10 6 4.00763 • 10 6	1 9 18	4.12072 • 10 8 1.04020 • 10 7 1.35728 • 10 7	1 2 3	1.29773 • 10 5 1.38250 • 10 5 1.38015 • 10 5

Table 5 .

 5 Lebesgue constant evolution when the AFP-based approach described in Section 4.2 is used.

	Example 2.6	Example 2.7	Example 2.8
	n = 100	n

Table 7 .

 7 Timings showing the effect of running our code with different options. As for Table6, the numerical values represent the averages in seconds over 50 executions. For the first seven lines, the double-precision version of our routine was used, while for the last one long double 80-bit operations were carried out. In all cases, our code was compiled using g++5.2.0 with -O3 -DNDEBUG level optimizations.

	(degree)	Uniform	GNURadio MATLAB SciPy	[3]
		(sequential)
	2.6 (n = 50)		0.0034	0.0029	0.0532	0.0711 0.0384/0.0098
	2.6 (n = 80)		0.0045	0.0068	0.1174	0.2587 0.0416/0.0316
	2.6 (n = 100)		0.0073	NC		0.1491	0.3511 0.0451/0.0396
	2.7 (n = 50)		0.0022	0.0023	0.0681	0.0971 0.0395/0.0174
	2.7 (n = 80)		0.0069	NC		0.2002	0.3492 0.0761/0.0293
	2.7 (n = 100)		0.0301	NC		NC	NC	0.2599/0.0521
	2.8 (n = 520)		0.2375	NC		NC	NC	1.1691/2.7011
	2.9 (n = 53248)		NC		NC		NC	NC	NC/NC
	Example (degree)	Uniform	Uniform	Scaling	Scaling	AFP	AFP
	(sequential) (parallel) (sequential) (parallel) (sequential) (parallel)
	2.6 (n = 50)	0.0034	0.0031	0.0021	0.0018	0.0021	0.0011
	2.6 (n = 80)	0.0045	0.0035	0.0029	0.0012	0.0029	0.0015
	2.6 (n = 100)	0.0073	0.0059	0.0099	0.0041	0.0035	0.0023
	2.7 (n = 50)	0.0022	0.0019	0.0024	0.0015	0.0026	0.0014
	2.7 (n = 80)	0.0069	0.0064	0.0058	0.0041	0.0029	0.0028
	2.7 (n = 100)	0.0301	0.0151	0.0123	0.0066	0.0113	0.0071
	2.8 (n = 520)	0.2375	0.1613	0.1632	0.0725	0.0731	0.0693
	2.9 (n = 53248)	NC		NC	537.8	162.6	NC	NC

See https://github.com/sfilip/firpm for the accompanying C++ code

All such constants in this subsection are defined over [-1, 1], x and unit weight function. Thus, we use the simplified notation Λ [-1,1],x when referring to them.