Efficient algorithms for the design of finite impulse response - PowerPoint PPT Presentation

Efficient algorithms for the design of finite impulse response digital filters Silviu Filip under the supervision of N. Brisebarre and G. Hanrot (AriC, LIP, ENS Lyon) Journées Nationales de Calcul Formel (JNCF) CIRM, Luminy, November 3-7, 2014 1 / 19

Digital Signal Processing became increasingly relevant over the past 4 decades: ANALOG → DIGITAL 2 / 19

Digital Signal Processing became increasingly relevant over the past 4 decades: ANALOG → DIGITAL think of: data communications (ex: Internet, HD TV and digital radio) audio and video systems (ex: CD, DVD, BD players) many more 2 / 19

Digital Signal Processing became increasingly relevant over the past 4 decades: ANALOG → DIGITAL think of: data communications (ex: Internet, HD TV and digital radio) audio and video systems (ex: CD, DVD, BD players) many more What are the ’engines’ powering all these? 2 / 19

Digital filters y [ n ] = x [ n ] ⋆ h [ n ] → we get two categories of filters finite impulse response ( FIR ) filters infinite impulse response ( IIR ) filters 3 / 19

Digital filters y [ n ] = x [ n ] ⋆ h [ n ] → we get two categories of filters finite impulse response ( FIR ) filters infinite impulse response ( IIR ) filters → natural to work in the frequency domain 3 / 19

Digital filters Y ( ω ) = X ( ω ) H ( ω ) , ω ∈ [0 , π ] → we get two categories of filters finite impulse response ( FIR ) filters infinite impulse response ( IIR ) filters → natural to work in the frequency domain H is the transfer function of the filter 3 / 19

Digital filters Y ( ω ) = X ( ω ) H ( ω ) , ω ∈ [0 , π ] → we get two categories of filters finite impulse response ( FIR ) filters H is a polynomial infinite impulse response ( IIR ) filters H is a rational fraction → natural to work in the frequency domain H is the transfer function of the filter 3 / 19

The filtering toolchain Steps: 1. derive a concrete mathematical representation of the filter → use theory of minimax approximation 4 / 19

The filtering toolchain Steps: 1. derive a concrete mathematical representation of the filter → use theory of minimax approximation 2. quantization of the filter coefficients using fixed-point or floating-point formats → use tools from algorithmic number theory (euclidean lattices) 4 / 19

The filtering toolchain Steps: 1. derive a concrete mathematical representation of the filter → use theory of minimax approximation 2. quantization of the filter coefficients using fixed-point or floating-point formats → use tools from algorithmic number theory (euclidean lattices) 3. hardware synthesis of the filter 4 / 19

The filtering toolchain Steps: 1. derive a concrete mathematical representation of the filter → use theory of minimax approximation 2. quantization of the filter coefficients using fixed-point or floating-point formats → use tools from algorithmic number theory (euclidean lattices) 3. hardware synthesis of the filter Today’s focus: first step for FIR filters 4 / 19

Finite Impulse Response (FIR) filters large class of filters, with a lot of desirable properties Usual representation: H ( ω ) = � L k =0 a k cos( ωk ) 5 / 19

Finite Impulse Response (FIR) filters large class of filters, with a lot of desirable properties Usual representation: H ( ω ) = � L k =0 a k cos( ωk ) = � L k =0 a k T k (cos( ω )) → if x = cos( ω ) , view H in the basis of Chebyshev polynomials 5 / 19

Finite Impulse Response (FIR) filters large class of filters, with a lot of desirable properties Usual representation: H ( ω ) = � L k =0 a k cos( ωk ) = � L k =0 a k T k (cos( ω )) → if x = cos( ω ) , view H in the basis of Chebyshev polynomials Specification: 5 / 19

Finite Impulse Response (FIR) filters large class of filters, with a lot of desirable properties Usual representation: H ( ω ) = � L k =0 a k cos( ωk ) = � L k =0 a k T k (cos( ω )) → if x = cos( ω ) , view H in the basis of Chebyshev polynomials Specification: 8 � H ( ω ) = a k cos( ωk ) k =0 5 / 19

Optimal FIR design with real coefficients The problem: Given a closed real set F , find an approximation H ( ω ) = � L k =0 a k cos( ωk ) of degree L for a continuous function D ( ω ) , ω ∈ F such that δ = � E ( ω ) � ∞ ,F = max ω ∈ F | H ( ω ) − D ( ω ) | is minimal . 6 / 19

Optimal FIR design with real coefficients The solution : characterized by the Alternation Theorem Theorem The unique solution H ( ω ) = � L k =0 a k cos( ωk ) has an error function E ( ω ) , for which there exist L + 2 values ω 0 < ω 1 < · · · < ω L +1 , belonging to F , such that E ( ω i ) = − E ( ω i +1 ) = ± δ, for i = 0 , . . . , L. 7 / 19

Optimal FIR design with real coefficients The solution : characterized by the Alternation Theorem Theorem The unique solution H ( ω ) = � L k =0 a k cos( ωk ) has an error function E ( ω ) , for which there exist L + 2 values ω 0 < ω 1 < · · · < ω L +1 , belonging to F , such that E ( ω i ) = − E ( ω i +1 ) = ± δ, for i = 0 , . . . , L. → well studied in Digital Signal Processing literature 1972: Parks and McClellan → based on a powerful iterative approach from Approximation Theory: 1932: Remez 7 / 19

The Parks-McClellan design method: Example 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0 8 / 19

The Parks-McClellan design method: Example 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0 8 / 19

The Parks-McClellan design method: Example 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0 8 / 19

The Parks-McClellan design method: Example 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0 8 / 19

The Parks-McClellan design method: Example 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0 8 / 19

The Parks-McClellan design method: Example 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0 8 / 19

The Parks-McClellan design method: Example 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0 8 / 19

The Parks-McClellan design method: Example 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 0 8 / 19

The Parks-McClellan design method: Steps 9 / 19

Step 1: Choosing the L + 2 initial references Traditional approach: take the L + 2 references uniformly from F → can lead to convergence problems 10 / 19

Step 1: Choosing the L + 2 initial references Traditional approach: take the L + 2 references uniformly from F → can lead to convergence problems → want to start from better approximations Existing approaches: not general enough and/or costly to execute 10 / 19

Step 1: Choosing the L + 2 initial references Traditional approach: take the L + 2 references uniformly from F → can lead to convergence problems → want to start from better approximations Existing approaches: not general enough and/or costly to execute Our approach : extrema position extrapolation from smaller filters → although empirical, it is rather robust in practice 10 / 19

Step 2: Computing the current error function E ( ω ) and δ Amounts to solving a linear system in a 0 , . . . , a L and δ .       1 cos( ω 0 ) · · · cos( Lω 0 ) 1 D ( ω 0 ) a 0 . . . . . . . . . . . .       . . . . . .  =             1 cos( ω L ) · · · cos( Lω L ) ( − 1) L D ( ω L ) a L      ( − 1) L +1 1 cos( ω L +1 ) · · · cos( Lω L +1 ) D ( ω L +1 ) δ → solving system directly: can be numerically unstable 11 / 19

Step 2: Computing the current error function E ( ω ) and δ Amounts to solving a linear system in a 0 , . . . , a L and δ .       1 cos( ω 0 ) · · · cos( Lω 0 ) 1 D ( ω 0 ) a 0 . . . . . . . . . . . .       . . . . . .  =             1 cos( ω L ) · · · cos( Lω L ) ( − 1) L D ( ω L ) a L      ( − 1) L +1 1 cos( ω L +1 ) · · · cos( Lω L +1 ) D ( ω L +1 ) δ → solving system directly: can be numerically unstable → use barycentric form of Lagrange interpolation [Berrut&Trefethen2004] 11 / 19

Barycentric Lagrange interpolation Problem: p polynomial with deg p � L interpolates f at points x j , i.e., p ( x j ) = f j , j = 0 , . . . , L 12 / 19

Barycentric Lagrange interpolation Problem: p polynomial with deg p � L interpolates f at points x j , i.e., p ( x j ) = f j , j = 0 , . . . , L → the barycentric form of p is: L w j � f j x − x j j =0 p ( x ) = , L w j � x − x j j =0 1 where w j = k � = j ( x j − x k ) . � Cost: O ( L 2 ) for computing all w j , O ( L ) for evaluating p ( x ) . 12 / 19

Barycentric Lagrange interpolation → numerically stable if the family of interpolation nodes used has a small Lebesgue constant [Mascarenhas&Camargo2014] The Lebesgue constant: specific to each grid of points; quantifies the convergence/divergence properties of polynomial interpolants using those nodes 13 / 19