Robust tools for weighted Chebyshev approximation and applications to digital filter design

136 Pages·2017·3.46 MB·English

Robust tools for weighted Chebyshev approximation and
applications to digital filter design
Silviu-Ioan Filip
To cite this version:
Silviu-IoanFilip. RobusttoolsforweightedChebyshevapproximationandapplicationstodigitalfilter
design. Signal and Image Processing. Université de Lyon, 2016. English. ￿NNT: 2016LYSEN063￿.
￿tel-01447081￿
HAL Id: tel-01447081
https://theses.hal.science/tel-01447081
Submitted on 26 Jan 2017
HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est
archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
Num´ero National de Th`ese : 2016LYSEN063
` ´
THESE de DOCTORAT DE L’UNIVERSITE DE LYON
op´er´ee par
´
l’Ecole Normale Sup´erieure de Lyon
E´cole Doctorale No 512
´
Ecole Doctorale en Informatique et Math´ematiques de Lyon
Sp´ecialit´e de doctorat :
Informatique
Soutenue publiquement le 07/12/2016, par :
Silviu-Ioan FILIP
Robust tools for weighted Chebyshev approximation and
applications to digital filter design
Outils robustes pour l’approximation de Chebyshev pond´er´ee et applications `a la synth`ese de filtres num´eriques
Devant le jury compos´e de :
Beckermann, Bernhard Professeur Universit´e de Lille Rapporteur
Directeur de Universit´e de Rennes
Sentieys, Olivier Rapporteur
recherche ENSSAT Lannion
Belfiore, Jean-Claude Professeur T´el´ecom ParisTech Examinateur
Directrice de
Leblond, Juliette Centre de recherche INRIA Examinatrice
recherche
Trefethen, Lloyd Nicholas Professeur University of Oxford Examinateur
Brisebarre, Nicolas Charg´e de recherche ENS de Lyon Co-directeur
Directeur
Hanrot, Guillaume Professeur ENS de Lyon
de th`ese
ii
Contents
1 Introduction 1
1.1 Minimax approximation - outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 Technical prerequisites 7
2.1 Chebyshev technology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1.1 Basic facts about Chebyshev polynomials . . . . . . . . . . . . . . . . . . . . . . . . 7
2.1.2 Chebyshev series and Chebyshev interpolants . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Digital filter design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.1 Frequency domain representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.2 Discrete-time systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.3 Minimax FIR filter design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3 Computer arithmetic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3.1 Floating-point arithmetic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.3.2 Fixed-point arithmetic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3 Efficient tools for weighted minimax approximation 21
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.2 Properties of best approximations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2.1 Unique best approximations: the Haar condition . . . . . . . . . . . . . . . . . . . . 23
3.2.2 When uniqueness is no longer guaranteed . . . . . . . . . . . . . . . . . . . . . . . . 25
3.3 Minimax exchange algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.3.1 The second Remez algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.3.2 The first Remez algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.4 Our setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.4.1 The role of multi-interval domains . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.4.2 Four key examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.5 Practical problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.5.1 Robustness issues: initialization is a key step . . . . . . . . . . . . . . . . . . . . . . 33
3.5.2 Scalability issues: speeding up the extrema search . . . . . . . . . . . . . . . . . . . 35
3.6 Initialization: two new heuristic approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.6.1 Lagrange interpolation problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.6.2 Polynomial meshes and an AFP-based approach . . . . . . . . . . . . . . . . . . . . 39
3.6.3 A reference scaling idea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.6.4 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.7 Computing the current approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
iii
iv
3.7.1 Barycentric interpolation formulas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.7.2 Numerical stability issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.8 Extrema search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.8.1 Chebyshev-proxy root-finding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.8.2 A new subdivision strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.8.3 Updating the current reference set . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.8.4 Retrieving the minimax coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.9 Implementation details . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.9.1 User interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.9.2 Timings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.10 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4 Machine-efficient approximations and filter design 59
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.1.1 State of the art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.2 The finite wordlength design problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.2.1 Exact approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.2.2 A heuristic approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.3 Euclidean lattices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.3.1 Some basic results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.3.2 The shortest vector problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.3.3 Lattice basis reduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.3.4 The closest vector problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.4 FIR filter design using the FPminimax algorithm . . . . . . . . . . . . . . . . . . . . . . . . 74
4.4.1 Suitable interpolation points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.4.2 An approximate solution for the resulting CVP instance . . . . . . . . . . . . . . . . 81
4.4.3 A theoretical estimate of the quality of our solution . . . . . . . . . . . . . . . . . . 81
4.5 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.6 Conclusion and future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5 A complete design chain for FIR filter synthesis on FPGAs 89
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.1.1 Why FPGAs? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.1.2 Digital filters: from specification to implementation . . . . . . . . . . . . . . . . . . 90
5.1.3 Implementation parameter space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
5.1.4 Contributions and outline of the present work . . . . . . . . . . . . . . . . . . . . . . 92
5.2 Background and state of the art. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
5.3 Integrated hardware FIR filter design. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
5.3.1 Computing the optimal filter out of the I/O format specification . . . . . . . . . . . 93
5.4 Examples and results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
5.5 Conclusion and future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6 Rational minimax approximation problems 99
6.1 The setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
6.2 The rational exchange algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
6.2.1 Initialization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
6.2.2 Determining the current approximant . . . . . . . . . . . . . . . . . . . . . . . . . . 103
6.2.3 Updating the reference set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
6.3 More examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
6.4 Conclusion & perspectives for future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
List of Figures
2.1 T to T on [ 1,1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
0 3
−
3.1 Error curve for the minimax approximation of x2 over [0,2] using the basis (x,ex). . . . . . 29
3.2 Tolerance scheme and ideal frequency response for a lowpass filter. . . . . . . . . . . . . . . 30
3.3 Design example showcasing a transition band anomaly and how it can be removed. . . . . . 31
3.4 Degree n=10 minimax approximation of f(x)=ln(x)esin(x) in terms of relative error. . . . 32
3.5 Relativeerrors(correctsignificantdigits)whencomputingthestartingpforExample3.20with
both types of barycentric formulas. Uniform initialization (ΛRn[x],[−1,1](x)(cid:39)4.24973·1014) is
used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.6 Relative errors (correct significant digits) when computing the starting p for Example 3.20
using both types of barycentric formulas. Reference scaling (ΛRn[x],[−1,1](x)(cid:39)1.29773·105)
is used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.7 Our interval subdivision strategy for computing the extrema of e(x). . . . . . . . . . . . . . 51
3.8 The value of the approximation error e(cid:63) on a small interval [0.247168,0.24869] of X, for the
minimax result satisfying the specification for Example 3.21. . . . . . . . . . . . . . . . . . 54
4.1 The Z2 lattice with unit basis vectors e1 = [1 0]T and e2 = [0 1]T and its associated
fundamental parallelepiped. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.2 The Z2 lattice, now with basis vectors b1 = [1 1]T and b2 = [1 2]T and its fundamental
parallelepiped. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.3 A lattice for which the (cid:96) and (cid:96) SVP solutions (b and c respectively) are different. . . . . 68
2 ∞
4.4 Successive minima of a two-dimensional lattice (in the (cid:96) norm). . . . . . . . . . . . . . . . 68
2
4.5 The relative lengths of successive Gram-Schmidt vectors before (left column) and after (right
column) LLL-reduction, for Example 4.2. The top row uses the final reference set of p(cid:63), the
middle row is for the zeros of e(cid:63) and endpoints of X, while the last row corresponds to an
n+2-point AFP set. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.6 TherelativelengthsofsuccessiveGram-Schmidtvectorsforthedesignofadegreen=62type
I filter adhering to specification B from Table 4.2, with 22-bit (21 bits for the fractional part)
filter taps. The top row corresponds to the situation before any basis reduction algorithm
was applied, whereas the middle row shows the effect of LLL reduction. BKZ (with block
size β =8) and HKZ reductions give the same result (bottom line). An AFP grid containing
n+2 elements was used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
v
vi LIST OF FIGURES
4.7 TherelativelengthsofsuccessiveGram-Schmidtvectorsforthedesignofadegreen=62type
I filter adhering to specification B from Table 4.2, with 22-bit (21 bits for the fractional part)
filter taps. The top row corresponds to the situation before any basis reduction algorithm
was applied, whereas the bottom row shows the effect of LLL reduction. BKZ (with block
size β =8) and HKZ reductions give very similar results, which is why they are not shown.
A grid of n+2 equispaced points inside the convex hull of X was used. . . . . . . . . . . . 80
4.8 Stopband attenuation for the minimax design of a given lowpass specification, alongside the
16-bit fixed-point designs of the error shaped approach and our LLL-based algorithm. . . . 85
5.1 Simplified depiction of a generic FPGA layout. . . . . . . . . . . . . . . . . . . . . . . . . . 90
5.2 An ideal FIR filter architecture. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
5.3 An actual FIR filter architecture. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
5.4 MATLAB design flow with all the parameters. The dashed ones can be computed by the
tool or input by the user. The external minimizecoeffwl function can help compute the
coefficient formats, while the quantized coefficients themselves can be computed using one of
minimizecoeffwl, constraincoeffwl or maximizestopband. . . . . . . . . . . . . . . . . 91
5.5 Proposed design flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
6.1 Minimax approximation errors for Example 6.3 when a = 0.8 (top left), a = 0.9 (top
− −
right), a= 0.99 (bottom left) and a= 1 (bottom right). In all four cases the abscissa is
− −
linearly scaled to [ 1,1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
−
6.2 Approximation errors after the first iteration of the exchange algorithm has been executed
(left) and at the end of execution (right) for the degree (8,8) best approximation of f(x)=
x 0.4 sin(x 0.2),x [ 1,1]. The extrema of the degree (8,8) CF approximation to f
| − | | − | ∈ −
were used as starting reference vector. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
6.3 The error function for the best minimax approximation of f(x)= 1 over [ 1,1].106
R4,4 log(x+2.005) −
6.4 The error function for the best minimax approximation of f(x)=√x over [0,1]. . . 107
40,40
R
List of Tables
2.1 Parameters of the binary floating-point formats from the 754-2008 standard [107]. . . . . . 17
3.1 Link between the two representations used interchangeably in the rest of the text. . . . . . 30
3.2 Default weighting values for linear-phase FIR filters (x=cos(ω)). . . . . . . . . . . . . . . . 42
3.3 NumberofreferencevaluesineachbandforthebandstopspecificationofExample3.19before
the first iteration and at the end of execution. . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.4 Iteration count comparison for uniform/reference scaling/ AFP-based initialization . . . . . 43
3.5 Statistics for the exchange algorithm execution on Example 3.18. . . . . . . . . . . . . . . . 43
3.6 Statistics for the exchange algorithm execution on Example 3.19. . . . . . . . . . . . . . . . 43
3.7 Lebesgueconstantevolutionduringtheexecutionoftheexchangealgorithm(i.e.,first,middle
and last iteration). The starting reference set is computed using uniform intialization. . . . 48
3.8 Lebesgue constant evolution when the reference scaling approach described in Section 3.6 is
used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.9 Lebesgue constant evolution when the AFP-based approach described in Section 3.6 is used. 48
3.10 Relative error of h with respect to δ for Example 3.18, n=100. . . . . . . . . . . . . . . 54
k X
3.11 Relative error of h with respect to δ for Example 3.19, n=100. . . . . . . . . . . . . . . 54
k X
3.12 Runtime (real time) comparisons of our IEEE-754 double-precision implementation of the
second Remez 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 [6]. . . . . . . . . . . . . . 56
3.13 Timings (real time) showing the effect of running our code with different options. As for
Table 3.12, 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. . . . . . . . . . . . . . . . . . . . . . . . 57
4.1 Scaled basis functions for finite word length FIR design. . . . . . . . . . . . . . . . . . . . . 62
4.2 Type I filter specifications considered in [124]. . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.3 Lattice basis reduction algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
vii
viii LIST OF TABLES
4.4 Approximation error comparison for the three different discretization choices. LLL basis
reduction is used. We have taken t =2,t =n for n<40 and t =1,t = n/2 for n(cid:62)40
1 2 1 2
(cid:98) (cid:99)
to limit the vicinity search runtime. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.5 Approximation error and execution times for LLL, BKZ and HKZ reduction.. . . . . . . . 83
4.6 Comparison with optimal results and the heuristic approach from [124]. . . . . . . . . . . . 84
4.7 Adaptive weighting improvements on the examples in Table 4.6. . . . . . . . . . . . . . . . 85
4.8 High degree type I FIR specifications. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.9 High degree quantization results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.1 Virtex6 synthesis results and accuracy measurements for generated FIR filters. . . . . . . . 96
6.1 Rational minimax approximation errors for ex on ( ,0] of type (n,n),n=1,...,60. . . . 108
−∞
Notation
N The set of natural numbers 0,1,... .
{ }
N(cid:63) The set of natural numbers, without 0: 1,2,... .
{ }
Z The set of integers ..., 2, 1,0,1,2,... .
{ − − }
R The set of real numbers.
C The set of complex numbers.
R[x] The set of univariate polynomials with real coefficients.
Rn[x] The set of univariate polynomials with real coefficients of degree at most n.
Rn[x1,...,xd] The set of d-variate polynomials with real coefficients of degree at most n.
(X) The set of all continuous functions over the set X.
C
sgn(x) The sign function: returns 1 if x>0, -1 if x<0 and is undefined if x=0.
a Small bold letters are used to represent column vectors.
A Capital bold letters are used to represent matrices.
At The transpose of the matrix A.
A(cid:63) The conjugate transpose of the matrix A (i.e., At).
diag(a) The diagonal matrix whose diagonal is formed by taking the elements of a.
ker(A) The kernel of the matrix A, that is, the set of vectors x for which Ax=0.
([a,b],w(x),dx) The space of square-integrable functions on [a,b] (i.e., if f ([a,b],w(x),dx),
2 2
L ∈L
then b f(x)2w(x)dx< ).
a | | ∞
(cid:82)
ix