Approximation of Almost Equiripple Low-pass FIR Filters

Report 2 Downloads 39 Views
Approximation of Almost Equiripple Low-pass FIR Filters Miroslav Vlˇcek

Pavel Zahradnik

Czech Technical University in Prague Faculty of Transportation Sciences Department of Applied Mathematics Na Florenci 25, CZ-110 00 Praha 1, Czech Republic Email: [email protected]

Czech Technical University in Prague Faculty of Electrical Engineering Department of Telecommunication Engineering Technick´a 2, CZ-166 27 Praha 6, Czech Republic Email: [email protected]

Abstract—A novel approximation of linear phase almost equiripple low-pass finite impulse response filter is introduced. Its frequency response closely approaches the frequency response of an optimal equiripple low-pass finite impulse response filter. The presented approximation is based on the generating polynomial which it is related to the class of iso-extremal polynomials. The closed form solution for an algebraic evaluation of the impulse response of the filter has been developed on the basis of a generalization of the differential equation suitable for the halfband specifications. No numerical procedures are involved. The practical design procedure based on the developed approximation is presented. An example of the design is included.

I.

I NTRODUCTION

A low-pass (LP) finite impulse response (FIR) filter belongs to the basic filters which are frequently used in digital signal processing. The most efficient LP FIR filters are the equiripple filters. These are usually designed by the McClellanParks [1] numerical optimization, which provides the coefficients of the impulse response of the filter. Despite its practicality it does not solve the problem of approximation theory. The approximation problem of an equiripple LP FIR filter remains still unresolved. This paper is seen as an extension of our previous activity focused on almost equiripple halfband FIR filters [2]. Here we are primarily concerned with the approximation of a linear phase almost equiripple LP FIR filter. The frequency response of such filter closely approaches the frequency response of an optimal (in Chebyshev sense) equiripple LP FIR filter. Here, we present the differential equation for the generating polynomial of the almost equiripple LP FIR filter. Based on the differential equation we have developed an efficient algorithm for algebraic evaluating the impulse response of the filter. II.

BASIC T ERMS

We assume an impulse response h(k) of type I, i.e. with odd length of N = 2n+1 coefficients and with even symmetry a(0) = h(n) , a(k) = 2h(n + k) = 2h(n − k) , k = 1 ...n . (1) The transfer function of the filter is " # n 2n X X −k −n H(z) = h(k)z = z a(0)+ a(k) Tk (w) = z −n Q(w) k=0

k=1

(2)

where Tk (w) is Chebyshev polynomial of first kind. The function Q(w) represents a polynomial in the variable w = 0.5(z + z −1 ) which on the unit circle z = ejωT reduces to the real valued zero phase transfer function Q(w) of the real argument w = cos(ωT ). Consequently we use the standard interval w ∈ [−1, 1] for polynomial approximation.

III.

G ENERATING P OLYNOMIAL

The generating polynomial for the low-pass FIR filter is related to the Zolotarev polynomial [4]. The extremal values of the Zolotarev polynomial Zp,q (w, κ) of degree n = p + q alternate between -1 and +1 (p+1)-times in the interval [wp , 1], and (q + 1)-times in the interval [−1, ws ]. By inspection, the Zolotarev polynomial satisfies the differential equation 2 dZp,q (w, κ) dw  2 2 2 = n (w − wm ) 1 − Zp,q (w, κ) .

(1 − w2 )(w − wp )(w − ws )



(3)

This nonlinear differential equation expresses the fact that ′ the first derivative Zp,q (w, κ) does not vanish at the points w = ±1, wp , ws where Zp,q (w, κ) = ±1 for which the right hand side of eq. (3) vanishes, and that w = wm is a turning point corresponding to the local extrema at which Zp,q (w, κ) 6= ±1. Elliptic modulus κ is the driving parameter of all elliptic functions [4] appearing in the parametric solution of (3). We have derived the linear differential equation of the second order   d2 d (w−wp )(w−ws )(w−wm ) (1 − w2 ) 2 −w Zp,q (w, κ) dw dw   wp + ws 2 + (w−wm )(w− )−(w−wp )(w−ws )(1−w ) 2 d × Zp,q (w, κ) dw +n2 (w − wm )3 (w)Zp,q (w, κ) = 0 . (4) Assuming

the

second

solution

of

(4)

in

form

10

IV.

Z ERO P HASE T RANSFER F UNCTION OF E QUIRIPPLE LP FIR F ILTER

AN

A LMOST

5

The zero phase transfer function Q(w) (Fig. 2) of an almost equiripple LP FIR filter is a normalized primitive function of the generating polynomial Sp,q (w, κ) Z 1 1 Q(w) = −N1 + Sp,q (w, κ)dw = −N1 + S(w) N2 −N1 N2 −N1 (7) where the norming values N1 , N2 are ( S(w = −1) = |P (ωT = π)| for q even N1 = nπ ) for q odd (8) S(w = cos( nnπ )) = P (ωT = +1 n+1

−5

−10

−15

−20

−25

−30 −1

Fig. 1.

−0.8

−0.6

−0.4

−0.2

0

w

0.2

0.4

0.6

0.8

1

Generating polynomial S19,11 (w, 0.55)

√ 1 − w2 Sp,q (w, κ) we obtain the differential equation

1

  d d2 Sp,q (w, κ) (w−wp )(w−ws )(w−wm ) (1−w2 ) 2 −3w dw dw     wp + ws + (w − wm ) w − − (w − wp )(w − ws ) 2 d ×(1 − w2 ) Sp,q (w, κ) dw "   wp +ws + n(n+2)(w−wm)3− wm − (w+wm )(w+2wm ) 2  2  2 !# wp +ws wp +ws + wm wm − +wm wp ws − 2 2 ×Sp,q (w, κ) = 0 .

(5)

√ The function 1 − w2 Sp,q (w, κ) is the iso-extremal solution of the approximation equation (3) with the polynomial Sp,q (w, κ) which forms the generating polynomial (Fig. 1) for the almost equiripple LP FIR filter specifications. With respect to the impulse response of the filter, it is advantageous to express the generating polynomial in the form of an expansion in Chebyshev polynomials of second kind Um (w)

Sp,q (w, κ) =

n X

( S(w = 1) = |P (ωT = 0)| for p even N2 = (9) π π S(w = cos( n + 1 )) = P (ωT = n + 1 ) for p odd and P (ejωT ) is the amplitude frequency response corresponding to the polynomial S(w). The amplitude frequency response |H(ejωT )| corresponding to the zero phase transfer function Q(w) from Fig. 2 is shown in Fig. 3. The algebraic

0.9 0.8 0.7 0.6

Q(w)

S19,11 (w, 0.55)

0

0.5 0.4 0.3 0.2 0.1 0 −1

−0.8

−0.6

−0.4

(6)

m=0

By inserting expansion (6) into equation (5) we can derive the recursive formula for the coefficients a+ (m), see Tab. I. The recursion for the coefficient a+ (m) is closely related to the algorithm for the Zolotarev polynomials [4]. In both algorithms we have used the standard notation [3] for Jacobi’s elliptic functions sn(u0 |κ), cn(u0 |κ), dn(u0 |κ) and cd(u0 |κ), for the complete elliptic integral of the first kind K(κ) of modulus κ and for Jacobi’s Zeta function Z(u0 |κ). The value u0 is a fraction of quarter-period, namely u0 = p K(κ)/n.

0

w

0.2

0.4

0.6

0.8

1

Fig. 2. Zero phase transfer function Q(w) for p = 19, q = 11 and κ = 0.55

algorithm for the impulse response coefficients h(k) is robust and provides easily polynomials with tens of thousands of coefficients. For illustration, the remarkable selectivity of the generating polynomial is shown in Fig. 4. V.

a+ (m)Um (w) .

−0.2

D ESIGN P ROCEDURE

The main objective in filter design is to evaluate the filter degree satisfying filter specification and also to evaluate the impulse response coefficients. Based on the differential equation (5) we have derived a simple procedure for evaluation of the impulse response h(k) of the filter (Tab. I). Its length is 2(p + q) + 3 coefficients. The design procedure consists of few steps: ◦ Specify the minimum attenuation in the stop-band asdB = 20 log(as ) or the maximum attenuation in the pass-band apdB = 20 log(ap ). They are related by ap + as = 1. Further, specify the pass-band frequency ωp T and the stopband frequency ωs T . The maximum width of the transition

TABLE I.

given initialisation

p, q, κ n = p + q,

u0 =

E VALUATION OF THE IMPULSE RESPONSE h(k)

2p + 1 K(κ), 2n + 2

wp = 2 cd2 (u0 |κ) − 1,

ws = 2 cn2 (u0 |κ) − 1

wp + ws sn(u0 |κ)cn(u0 |κ) , wm = ws + 2 Z(u0 |κ) 2 dn(u0 |κ) α(n) = 1, α(n + 1) = α(n + 2) = α(n + 3) = α(n + 4) = α(n + 5) = 0

wq = body (for

m = n + 2 to 3) 8c(1) = n(n + 2) − (m + 3)(m + 5) 4c(2) = 3wm [n(n + 2) − (m + 2)(m + 4)] + (m + 3)(2m + 7)(wm − wq ) 8c(3) = 3[n(n + 2) − (m + 1)(m + 3)] + 12wm [(n + 1)2 wm − (m + 2)2 wq ] −4(m + 2)(m + 3)(wp ws − wm wq ) 2c(4) = 3[(n + 1)2 wm − (m + 1)2 wq ] − (m + 1)2 (wm − wq ) 2 − (m + 1)2 w w ] +2wm [(n + 1)2 wm p s 8c(5) = 3[n(n + 2) − (m − 1)(m + 1)] + 12wm [(n + 1)2 wm − m2 wq ] −4m(m − 1)(wp ws − wm wq ) 4c(6) = 3wm [n(n + 2) − (m − 2)m] + (m − 1)(2m − 3)(wm − wq ) 8c(7) = n(n + 2) − (m − 3)(m − 1) 1 α(m − 3) = c(7)

(end normalisation

loop on m) s(n) =

n X

6 X

(−1)ℓ c(ℓ)α(m + 4 − ℓ)

ℓ=1

α(k)(k + 1)

k=0

(for

k = 0 to n) a+ (k) = (−1)p (n + 1)

(end integration (for

(end impulse response (for

loop on k) k = 0 to n) a+ (k) a(k + 1) = k+1 loop on k) h(n + 1) = −N1 k = 1 to n+1) h(n + 1 ± k) =

(end

α(k) s(n)

a(k) 2(N2 − N1 )

loop on k)

band is ∆ωT = (ωs − ωp )T (Fig. 3). The degree equation relates the degree n of the generating polynomial Sp,q (w, κ), the maximum width of the transition band ∆ωT and the minimum attenuation in the stopband asdB . ◦ Calculate the degree n = p + q of the generating polynomial Sp,q (w, κ) using the approximative degree equation ξ4 (ξ1 n + ξ2 )∆ωT /π + ξ3 + = asdB , (10) (n + ξ5 )ξ6 where the constants ξi are ξ1 = −14.02925485, ξ2 = −32.86119410, ξ3 = −5.80117336, ξ4 = 2.99564719, ξ5 = −21.24188066 and ξ6 = 0.28632078. Equation (10)

approximates the empirical data with excellent accuracy. ◦ Determine the integer values p, q   (ωs + ωp )T p = round n , q =n−p . 2π

(11)

◦ Calculate the elliptic modulus κ κ=

s

1−



1−k 1+k

2

,

(12)

1 0.9

0 −10

ap

0.8

20 log |H(ejωT )| [dB]

−20

|H(ejωT )|

0.7 0.6 0.5 0.4 0.3 0.2

−30 −40 −50 −60 −70 −80

a

s

0.1 0

ω T/π p

0

0.1

0.2

0.3

0.4

0.5

0.6

ωT /π

−90 −100

ω T/π s

0.7

0.8

0.9

1

Fig. 3. Amplitude frequency response |H(ejωT )| for p=19, q=11, κ=0.55

0

0.1

0.2

0.3

0.4

0.5

0.6

ωT /π

0.7

0.8

0.9

1

Fig. 5. Amplitude frequency responses 20 log |H(ejωT )| [dB] for p = 64, q = 160 (11) and κ = 0.467363

1

frequency ωp T = 0.27π, stop-band frequency ωs T = 0.3π and by the minimal attenuation in the stop-band asdB = −100 dB. Based on the presented design procedure, we get the degree of the generating polynomial n = 223.030277 → 224 (10), integer indices p = 64, q = 160 (11) and elliptic modulus κ = 0.467363 (13). The length of the filter is N = 451 coefficients. The actual filter properties are ωp T = 0.271898π, ωs T = 0.301728π and asdB act = −100.36 dB. The amplitude frequency responses 20log|H(ejωT )| [dB] of the designed lowpass filter and of the magnitude complementary high-pass filter are shown in Fig. 5. Their sum is equal 0dB.

0.9 0.8

|H(ejωT )|

0.7 0.6 0.5 0.4 0.3 0.2

VII.

0.1 0 0

0.1

0.2

0.3

0.4

0.5

0.6

ωT /π

0.7

0.8

0.9

1

Fig. 4. Amplitude frequency response |H(ejωT )| for p = 600, q = 200 and κ = 0.25

where k is approximated by the formula    χ2 k= χ1 + n + χ5 p + χ6 wp (13) (p + χ3 )χ4 χ8 1 +χ7 + + χ 10 (p + χ9 ) (n + χ11 p + χ12 )χ13 p+χ14 for wp = cos [(π − ∆ωT )/2] and χ1 = −0.00452871, χ2 = 0.51350112, χ3 = 2.56407699, χ4 = 1.12297611, χ5 = 0.01473844, χ6 = 0.14824220, χ7 = 0.00245539, χ8 = 0.52499043, χ9 = 0.75104615, χ10 = 1.29448910, χ11 = −1.06038228, χ12 = 0.64247743, χ13 = −0.00932499, χ14 = 1.88486768. Equation (13) was obtained by approximating the empirical values. ◦ Evaluate the impulse response h(k) (Tab. I). VI.

E XAMPLE OF

THE

D ESIGN

Let us design a pair of magnitude complementary FIR filters. The low-pass FIR filter is specified by the pass-band

C ONCLUSION

We have introduced an approximation of the almost equiripple low-pass FIR filter based on the generating polynomial which is related to the Zolotarev polynomial. We have presented the differential equation for the generating polynomial. Employing the differential equation we have derived an algebraic procedure for an efficient evaluation of the impulse response of the filter. The practicality of the design procedure has been demonstrated. ACKNOWLEDGMENT The authors wish to thank the Grant Agency of the Czech Republic for supporting this research through project P-102/11/ 1795: Novel selective transforms for non-stationary signal processing. R EFERENCES [1] J.H. McClellan, T. W. Parks and L. R. Rabiner, A Computer Program for Designing Optimum FIR Linear Phase Digital Filters, IEEE Trans. Audio Electroacoust., Vol. AU - 21, Dec. 1973, pp. 506 - 526. [2] P. Zahradnik, M. Vlcek and R. Unbehauen, Almost Equiripple FIR Half-Band Filters, IEEE Transactions on Circuits and Systems - I: Fundamental Theory and Applications, Vol. CAS - 46, No. 6, June 1999, pp. 744 - 748. [3] M. Abramowitz, I. Stegun, Handbook of Mathematical Function, Dover Publication, New York Inc., 1972. [4] M. Vlcek and R. Unbehauen, Zolotarev Polynomials and Optimal FIR Filters, IEEE Transactions on Signal Processing, Vol. 47, No. 3, March 1999, pp. 717-730.