Introduction to Continuous Time Filter Design

MASSACHUSETTS INSTITUTE OF TECHNOLOGY DEPARTMENT OF MECHANICAL ENGINEERING 2.161 Signal Processing – Continuous and Discrete Introduction to Continuous Time Filter Design 1 1 Introduction The design techniques described here are based on the creation of a prototype low-pass filter, with subsequent transformation to other filter forms (high-pass, band-pass, band-stop) if necessary. 2 Low-Pass Filter Design The prototype low-pass filter is based upon the magnitude-squared of the frequency response function |H(jΩ)|2 , or the frequency response power function. The phase response of the filter is not considered. We begin by defining tolerance regions on the power frequency response plot, as shown in Fig. 1. | H ( jW ) |2 1 1 2 1 + e 1 1 + l 0 2 0 p a s s b a n d W c tr a n s itio n b a n d W r s to p b a n d W (ra d /s e c ) Figure 1: Tolerance regions in the frequency response plot. The filter specifications are that 1 ≥ |H(jΩ)|2 > 1 1 1 + 2 D. Rowell, October 2, 2008 1 for |Ω| < Ωc (1) and |H(jΩ)|2 < 1 1 + λ2 for |Ω| > Ωr , (2) where Ωc is the cut-off frequency, Ωr is the rejection frequency, and  and λ are design parameters that select the filter attenuation at the two critical frequencies. For example, if  = 1, at Ωc the power response response |H(jΩc )|2 = 0.5, the -3 dB response frequency. In general we expect the response function to be monotonically decreasing in the transition band. The filter functions examined in this document will be of the form |H(jΩ)|2 = 1 . 1 + f 2 (Ω) (3) where f (Ω) → 0 as Ω → 0, and f (Ω) → ∞ as Ω → ∞ to generate a low-pass filter action. 2.1 The Butterworth Filter The Butterworth filter is defined by the power gain |H(jΩ)|2 = 1 1 + 2 (Ω/Ωc )2N (4) where N is a positive integer defining the filter order. Note that λ does not appear in this formulation, but clearly N and λ are interrelated, since at Ω = Ωr 1 1 ≥ 2 1+λ 1 + 2 (Ωr /Ωc )2N (5) which may be solved to show log(λ/) (6) log(Ωr /Ωc ) The power response function of Butterworth filters for N = 1 . . . 5 is shown in Fig. 2. Butterworth filters are also known as “maximally flat” filters because the response has the maximum number of vanishing derivatives at Ω = 0 and Ω = ∞ for filters of the form of Eq. 3. N≥ 2.1.1 The Transfer Function of the Butterworth Filter The poles of the power gain transfer function may be found from the characteristic equation � 1+ 2 s jΩc �2N =0 (7) which yields 2N roots (poles) that lie on a circle: pn = Ωc −1/N ejπ(2n+N −1)/2N n = 1 . . . 2N (8) with radius r = Ωc −1/N , and angular separation of π/N rad. Notice that if N is odd a pair of poles will lie on the real axis at s = ±Ωc −1/N , while if N is even the roots will form complex conjugate pairs. 2 |H(jω)| 2 1 Power response 0.9 0.8 0.7 0.6 0.5 0.4 0.3 N increasing 0.2 N=1 0.1 N=5 0 0 0.5 1 1.5 2 2.5 3 Normalized frequency 3.5 4 4.5 5 ω/ω c Figure 2: Power response of Butterworth filters for N = 1 . . . 5 ( = 1). Figure 3a shows the six poles of |H(s)|2 for a third-order (N = 3) Butterworth filter. For a causal system we must have |H(jΩ)|2 = H(jΩ)H(jΩ) = H(s)H(−s)|s=jΩ which allows us to take take the N poles of |H(jΩ)|2 in the left half-plane as the poles of the filter H(s), that is the poles specified by (n = 1 . . . N ) in Eq. (8) above, pn = Ωc −1/N ejπ(2n+N −1)/2N n = 1...N as is shown in Fig. 3b. If the filter is to have unity gain at low frequencies lim |H(jΩ)| = 1 Ω→0 we require the complete Butterworth transfer function to be (−p1 )(−p2 ) . . . (−pN ) (s − p1 )(s − p2 ) . . . (s − pN ) (−1)N p1 p2 . . . pN = , (s − p1 )(s − p2 ) . . . (s − pN ) H(s) = where only the N stable left-half-plane poles are included. 3 (9) r = W n = 2 c e -1 /N n = 1 jW x s - p la n e x n = 6 2 p /6 x n = 3 x x r = W x n = 2 s n = 5 c e -1 /N n = 1 jW s - p la n e x p /3 x s n = 3 x n = 4 (a ) (b ) Figure 3: The poles of (a) |H(s)|2 , and (b) H(s) for a third-order (N = 3) Butterworth filter. jW r = W c s - p la n e x p /4 s x Figure 4: Poles of a second-order Butterworth filter (N = 2,  = 1). 4 Example Consider a second-order (N = 2,  = 1) Butterworth filter with cut-off frequency Ωc . Equation (9) generates a pole pair as shown shown in the pole-zero of Fig. 4. The transfer function is that of a second-order system with damping ratio ζ = 0.707 and undamped natural frequency Ωn = Ωc . Show that the power frequency response satisfies the Butterworth specification of Eq. (4). The transfer function for the second-order filter is H(s) = s2 + √ Ω2c 2Ωc s + Ω2c (10) and the frequency response is H(jΩ) = H(s)|s=jΩ = |H(jΩ)|2 = H(jΩ)H(jΩ) = Ω2c √ (Ω2c − Ω2 ) + j 2Ωc Ω (11) Ω4c 1 = 4 4 (Ωc + Ω ) 1 + (Ω/Ωc )4 (12) which is of the Butterworth form of Eq. (4). Example Design a Butterworth low-pass filter to meet the power gain specifications shown in Fig. 5. Comparing Figs. 1 and 5 1 = 0.9 −→  = 0.3333 1 + 2 1 = 0.05 −→ λ = 4.358 1 + λ2 Then N≥ log(λ/) = 3.70 log(Ωr /Ωc ) we therefore select N=4. The 4 poles (p1 . . . p4 ) lie on a circle of radius r = Ωc −1/N = 13.16 and are given by |pn | = 13.16  pn = π(2n + 3)/8 for n = 1 . . . 4, giving a pair of complex conjugate pole pairs p1,4 = −5.04 ± j12.16 p2,3 = −12.16 ± j5.04 5 | H ( jW ) |2 0 .9 1 0 .0 5 0 0 1 0 p a s s b a n d tr a n s itio n b a n d 2 0 W s to p b a n d Figure 5: Filter specification for Butterworth design. Bode Diagram 50 Magnitude (dB) 0 −50 −100 −150 0 Phase (deg) −90 −180 −270 −360 −1 10 0 1 10 10 2 10 3 10 Frequency (rad/sec) Figure 6: Bode plots for fourth-order Butterworth design example. 6 The transfer function, normalized to unity gain, is H(s) = (s2 29993 + 10.07s + 173.2)(s2 + 24.32s + 173.2) and the filter Bode plots are shown in Fig. 6. 2.2 Chebyshev Filters The order of a filter required to met a low-pass specification may often be reduced by relaxing the requirement of a monotonically decreasing power gain with frequency, and allowing “ripple” to occur in either the pass-band or the stop-band. The Chebyshev filters allow these conditions: 1 Type 1 |H(jΩ)|2 = (13) 2 1 + � TN2 (Ω/Ωc ) 1 Type 2 |H(jΩ)|2 = (14) 2 2 1 + � (TN (Ωr /Ωc )/TN2 (Ωr /Ω)) Where TN (x) is the Chebyshev polynomial of degree N . Note the similarity of the form of the Type 1 power gain (Eq. (13)) to that of the Butterworth filter, where the function TN (Ω/Ωc ) has replaced (Ω/Ωc )N . The Type 1 filter produces an all-pole design with slightly different pole placement from the Butterworth filters, allowing resonant peaks in the pass­ band to introduce ripple, while the Type 2 filter introduces a set of zeros on the imaginary axis above Ωr , causing a ripple in the stop-band. The Chebyshev polynomials are defined recursively as follows T0 (x) T1 (x) T2 (x) T3 (x) .. . TN (x) = = = = 1 x 2x2 − 1 4x3 − 3x = 2xTN −1 (x) − TN −2 (x), N >1 (15) with alternate definitions TN (x) = cos(N cos−1 (x)) = cosh(N cosh−1 (x)) The Chebyshev polynomials have the min-max property: Of all polynomials of degree N with leading coefficient equal to one, the polynomial TN (x)/2N −1 has the smallest magnitude in the interval |x| ≤ 1. This “minimum maximum” amplitude is 21−N . 7 (16) (17) In low-pass filters given by Eqs. (13) and (14), this property translates to the following characteristics: Filter Butterworth Chebyshev Type 1 Chebyshev Type 2 2.2.1 Pass-Band Characteristic Maximally flat Ripple between 1 and 1/(1 + �2 ) Maximally flat Stop-Band Characteristic Maximally flat Maximally flat Ripple between 1 and 1/(1 + λ2 ) The Chebyshev Type 1 Filter With the power response from Eq. (13) |H(jΩ)|2 = 1 1+ �2 TN2 (Ω/Ωc ) and the filter specification from Fig. 1, the required filter order may be found as follows. At the edge of the stop-band Ω = Ωr |H(jΩr |2 = 1 1+ �2 TN2 (Ωr /Ωc ) so that ≤ 1 1 + λ2 � � λ ≤ �TN (Ωr /Ωc ) = � cosh N cosh−1 (Ωr /Ωc ) and solving for N N≥ cosh−1 (λ/�) cosh−1 (Ωr /Ωc ) (18) The characteristic equation of the power transfer function is � 1+� 2 TN2 s jΩc � � = 0 or TN s jΩc � =± j � Now TN (x) = cos(N cos−1 (x)), so that � cos N cos � −1 s jΩc �� =± j � (19) If we write cos−1 (s/jΩc ) = γ + jα, then s = Ωc (j cos (γ + jα)) = Ωc (sinh α sin γ + j cosh α cos γ) (20) which defines an ellipse of width 2Ωc sinh(α) and height 2Ωc cosh(α) in the s-plane. The poles will lie on this ellipse. Substituting into Eq. (16) � TN s jΩc � = cos (N (γ + jα)) = cos N γ cosh N α − j sin N γ sinh N α, 8 the characteristic equation becomes j cos N γ cosh N α − j sin N γ sinh N α = ± . � (21) � 0 for real x we require Equating the real and imaginary parts in Eq. (21), (1) since cosh x = cos N γ = 0, or (2n − 1)π n = 1, . . . , 2N (22) γn = 2N and, (2) since at these values of γ, sin N γ = ±1 we have α=± 1 1 sinh−1 N � (23) As in the Butterworth design procedure, we select the left half-plane poles as the poles of the filter frequency response. The design procedure is: 1. Determine the filter order cosh−1 (λ/�) N≥ cosh−1 (Ωr /Ωc ) 2. Determine α α=± 1 1 sinh−1 N � 3. Determine γn , n = 1 . . . N γn = (2n − 1)π 2N n = 1, . . . , N 4. Determine the N left half-plane poles pn = Ωc (sinh α sin γn + j cosh α cos γn ) n = 1, . . . , N 5. Form the transfer function (a) If N is odd H(s) = (b) If N is even H(s) = −p1 p2 . . . pN (s − p1 )(s − p2 ) . . . (s − pN ) 1 p1 p2 . . . pN 2 1 + � (s − p1 )(s − p2 ) . . . (s − pN ) The difference in the gain constants in the two cases arises because of the ripple in the pass-band. When N is odd, the response |H(j0)|2 = 1, whereas if N is even the value of |H(j0)|2 = 1/(1 + �2 ). 9 Example Repeat the previous Butterworth design example using a Chebyshev Type 1 design. From the previous example we have Ωc = 10 rad/s., Ωr = 20 rad/s., � = 0.3333, λ = 4.358. The required order is N≥ cosh−1 (λ/�) cosh−1 13.07 = = 2.47 cosh−1 (Ωr /Ωc ) cosh−1 2 Therefore take N = 3. Determine α: � � 1 1 1 α= sinh−1 = sinh−1 (3) = 0.6061 N � 3 and sinh α = 0.6438, and cosh α = 1.189. Also, γn = (2n − 1)π/6 for n = 1 . . . 6 as follows: n: γn : sin γn : cos γn : 1 π/6 √1/2 3/2 2 π/2 1 0 3 5π/6 1/2 √ − 3/2 4 5 7π/6 3π/2 −1/2 -1 √ − 3/2 0 6 11π/6 −1/2 √ 3/2 Then the poles are pn = Ωc (sinh α sin γn + j cosh α cos γn ) √ � � 1 3 p1 = 10 0.6438 × + j1.189 × = 3.219 + j10.30 2 2 p2 = 10 (0.6438 × 1 + j1.189 × 0) = 6.438 √ � � 1 3 p3 = 10 0.6438 × − j1.189 × = 3.219 − j10.30 2 2 √ � � 1 3 p4 = 10 −0.6438 × − j1.189 × = −3.219 − j10.30 2 2 p5 = 10 (−0.6438 × 0 − j1.189 × 0) = −6.438 √ � � 1 3 p6 = 10 −0.6438 × + j1.189 × = −3.219 + j10.30 2 2 and the gain adjusted transfer function of the resulting Type 1 filter is H(s) = (s2 750 + 6.438s + 116.5)(s + 6.438) The pole-zero plot for the Chebyshev Type 1 filter is shown in Fig. 7. 10 jW X X -6 .4 3 8 j1 0 .3 0 s - p la n e 0 -3 .2 1 9 s - j1 0 .3 0 X Figure 7: Pole-zero plot for third-order Chebyshev Type 1 design example. 2.2.2 The Chebyshev Type 2 Filter The Chebyshev Type 2 filter has a monotonically decreasing magnitude function in the pass­ band, but introduces equi-amplitude ripple in the stop-band by the inclusion of system zeros on the imaginary axis. The Type 2 filter is defined by the power gain function: |H(jΩ)|2 = 1 1+ T 2 (Ω /Ω ) �2 TN2 (Ωrr /Ω)c N (24) If we make the substitutions ν= Ωr Ωc Ω and �ˆ = 1 �TN (Ωr /Ωc ) Eq. 24 may be written in terms of the modified frequency ν |H(jν)|2 = �ˆ2 TN2 (ν/Ωc ) 1 + �ˆ2 TN2 (ν/Ωc ) (25) which has a denominator similar to the Type 1 filter, but has a numerator that contains a Chebyshev polynomial, and is of order 2N . We can use a method similar to that used in the Type 1 filter design to find the poles as follows: 1. First define a complex variable, say τ = µ + jν (analogous to the Laplace variable s = σ + jΩ used in the type 1 design) and write the power transfer function: |H(τ )|2 = �ˆ2 TN2 (τ /jΩc ) 1 + �ˆ2 TN2 (τ /jΩc ) 11 The poles are found using the method developed for the Type 1 filter, the zeros are found as the roots of the polynomial TN (τ /jΩc ) on the imaginary axis τ = jν. From the definition TN (x) = cos (N cos−1 (x)) it is easy to see that the roots of the Chebyshev polynomial occur at � (n − 1/2)π x = cos N � n = 1...N and from Eq. (25) the system zeros will be at � (n − 1/2)π τn = jΩc cos N � n = 1 . . . N. 2. The poles and zeros are mapped back to the s-plane using s = Ωr Ωc /τ and the N left half-plane poles are selected as the poles of the filter. 3. The transfer function is formed and the system gain is adjusted to unity at Ω = 0. Example Repeat the previous Chebyshev Type 1 design example using a Chebyshev Type 2 filter. From the previous example we have Ωc = 10 rad/s., Ωr = 20 rad/s., � = 1/3, λ = 4.358. The procedure to find the required order is the same as before, and we conclude that N = 3. Next, define Ωr Ωc 200 = Ω Ω 1 3 �ˆ = = = 0.1154 �TN (Ωr /Ωc ) T3 (2) ν = Determine α: � � α= 1 1 1 sinh−1 = sinh−1 (8.666) = 0.9520 N �ˆ 3 and sinh α = 1.1024, and cosh α = 1.4884. The values of γn = (2n − 1)π/6 for n = 1 . . . 6 are the same as the design for the Type 1 filter, so that the poles of |H(τ )|2 are pn = Ωc (sinh α sin γn + j cosh α cos γn ) √ � � 1 3 τ1 = 10 1.1024 × + j1.4884 × = 5.512 + j12.890 2 2 τ2 = 10 (1.1024 × 1 + j1.4884 × 0) = 11.024 12 √ � 1 3 10 1.1024 × − j1.488 × = 5.512 − j12.890 2 2 √ � � 1 3 10 −1.1024 × − j1.4884 × = −5.512 − j12.890 2 2 � � 1 10 −1.1024 × − j1.488 × 0 = −11.024 2 √ � � 1 3 10 −1.1024 × + j1.4884 × = −5.512 + j12.890 2 2 � τ3 = τ4 = τ5 = τ6 = The three left half-plane poles (τ4 , τ5 , τ6 ) are mapped back to the s-plane using s = Ωr Ωc /τ giving three filter poles p1 , p2 = −5.609 ± j13.117 p3 = −18.14 The system zeros are the roots of T3 (ν/jΩc ) = 4(ν/jΩc )3 − 3(ν/jΩc ) = 0 from the definition of TN (x), giving ν1 = 0 and ν2 , ν3 = ±j8.666. Mapping these back to the s-plane gives two finite zeros z1 , z2 = ±j23.07, z3 = ∞ (the zero at ∞ does not affect the system response) and the unity gain transfer function is −p1 p2 p3 (s − z1 )(s − z2 ) z1 z2 (s − p1 )(s − p2 )(s − p3 ) 6.9365(s2 + 532.2) = (s + 18.14)(s2 + 11.22s + 203.5) H(s) = The pole-zero plot for this filter is shown in Fig. 8. Note that the poles again lie on ellipse, and the presence of the zeros in the stop-band. 2.3 Comparison of Filter Responses Bode plot responses for the three example filters are shown in Fig. 9. While all filters meet the design specification, it can be seen that the Butterworth and the Chebyshev Type 1 filters are all-pole designs and have an asymptotic high-frequency magnitude slope of −20N dB/decade, in this case -80 dB/decade for the Butterworth design and -60 dB/decade for the Chebyshev Type 1 design. The Type 2 Chebyshev design has two finite zeros, with the result that its asymptotic high frequency response has a slope of only -20 dB/decade. Note also the singularity in the phase response of the Type 2 Chebyshev filter, caused by the two purely imaginary zeros. The pass-band and stop-band power responses are shown in Fig. 10. Notice that the design method developed here guarantees that the response will meet the specification at the cut-off frequency (in this case |H(jΩ)|2 = 0.9 at Ωc = 10. Other design methods (such as used by Matlab) may not use this criterion. 13 jW s - p la n e j2 0 o j2 3 .0 7 X s to p -b a n d j1 3 .1 2 j1 0 p a s s -b a n d X -1 8 .1 4 0 -5 .6 1 - j1 3 .1 2 X Magnitude (dB) 50 0 s o - j2 3 .0 7 Chebyshev Type 2 (N=3) Figure 8: Pole-zero plot for third-order Chebyshev Type 2 design example. −50 −100 Chebyshev Type 1 (N=3) Butterworth (N=4) Phase (deg) −150 360 180 Chebyshev Type 2 (N=3) 0 −180 −360 0 10 Chebyshev Type 1 (N=3) Butterworth (N=4) 1 2 10 10 Angular frequency (rad/sec) Figure 9: Comparison of Bode plots for the three design examples. 14 3 10 1.02 Butterworth (N=4) Power response 1 Chebyshev Type 2 (N=3) 0.98 0.96 0.94 0.92 Chebyshev Type 1 (N=3) 0.9 0.88 0 1 2 3 4 5 6 7 8 9 10 Angular frequency (rad/sec) 0.035 Power response 0.03 0.025 0.02 0.015 Butterworth (N=4) Chebyshev Type 1 (N=3) Chebyshev Type 2 (N=3) 0.01 0.005 0 20 30 40 50 60 70 80 90 100 Angular frequency (rad/sec) Figure 10: Comparison of the pass-band (top) and stop-band (bottom) power responses for the three design examples. 2.4 Other Filter Designs While the Butterworth and Chebyshev filters are perhaps the most common designs, there are other filter groups that have additional advantages. 2.4.1 Bessel Filters The Bessel (or sometimes known as the Thomson) filter has a phase response that most closely approximates a pure time delay, which is often an advantage if phase distortion is to be avoided. In addition the Bessel filters have a very small overshoot in their step response, which make them attractive for the transmission of pulse-like waveforms. The improvement in the phase response comes at the expense of the of the gain response. Although the gain response is monotonic, Bessel filters are not maximally flat, and do not have as steep a transition-band as other filter types. Bessel filters are all-pole filters. 2.4.2 Legendre Filters The filters, also known as optimally monotonic, combine properties of the Butterworth and Chebyshev designs. The Legendre filter has no ripples in its magnitude response, but it is not as flat as the Butterworth design. It has better transition-band cut-off rate than the Butterworth filter, and no ripple in either the passband or the stop-band. 2.4.3 The Elliptic Filter Also known as the Chebyshev-Cauer filter, allows a very sharp cut-off in the transition-band by allowing ripples in both the pass-band and in the stop-band. 15 2.5 MATLAB Low-Pass Filter Design Functions MATLAB filter design functions are contained in the Signal Processing Toolbox. See MAT­ LAB’s help for complete descriptions. • Continuous low-pass filter order calculation. The following return the required order, and equivalent -3dB cut-off frequency (Wn ) or stop-band edge (Ws ) to meet the given specifications: [N, [N, [N, [N, Wn] Wn] Ws] Wn] = = = = buttord(Wc, Wr, Rp, Rs, ’s’) - Butterworth filter. cheb1ord(Wc, Wr, Rp, Rs, ’s’) - Chebyshev Type I filter. Cheb2ord(Wc, Wr, Rp, Rs, ’s’) - Chebyshev Type II filter. ellipord(Wc, Wr, Rp, Rs, ’s’) - Elliptic filter. • Continuous low-pass filter prototypes. Return the zeros, poles, and gain for an nth order filter with a -3dB cut-off(� = 1), and cut-off frequency of Ωc = 1: [Z,P,K] [Z,P,K] [Z,P,K] [Z,P,K] = = = = buttap(N) - Butterworth filter prototype. cheb1ap(N) - Chebyshev Type I filter prototype. cheb2ap(N) - Chebyshev Type II filter prototype. ellipap(N) - Elliptic filter prototype. • Continuous filter design. Note that these functions will design both continuous and discrete filters, and that continuous filters require the last argument to be ’s’, for example [z,p,k] = butter(3,2*pi*150,’s’) will design as 3rd-order Butterworth filter with a -3 dB cut-off frequency of 150 Hz and return the zeros, poles and gain. A filter “object” may then be created using filt = zpk(z,p,k) Similarly, [num,den] = butter(3,2*pi*150,’s’) will design the same filter and return the transfer function numerator and denominator polynomial coefficients. A filter “object” may then be created using filt = tf(num, den) butter(N,Wn,’s’) - Butterworth filter design. cheby1(N,Rp,Wn,’s’) - Chebyshev Type I filter design. cheby2(N,Rs,Ws,’s’) - Chebyshev Type II filter design. ellip(N,Rp,Rs,Wn,’s’) - Elliptic filter design. besself(N,Wo) - Bessel analog filter design (continuous). Example Use Matlab to design a Chebyshev Type 2 low-pass filter with Ωc = 50 Hz, Ωs = 60 Hz, with maximum attenuation in the pass-band of 1 dB, and minimum attenuation in the stop-band of 40 dB. 16 The following MATLAB commands design and plot the required filter. (Note that cheb2ord() returns ws. the stop-band edge frequency, which is then passed to cheby2()): [n,ws] = cheb2ord(2*pi*50, 2*pi*60, 1, 40, ’s’); [num,den] = cheby2(n,40,ws,’s’); cheby2_lpf = tf(num,den); f=[0:1:100]; w=2*pi*f; [mag,phase]=bode(cheby2_lpf); plot(f, squeeze(mag)) generates the following plot Chebyshev Type 2 Filter 1 0.9 Frequency response magnitude 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 3 0 20 40 60 Frequency (Hz) 80 100 Transformation to Other Filter Classes Filter specification tolerance bands for high-pass, band-pass and band-stop filters are shown in Fig. 11. The most common procedure for the design of these filters is to design a prototype low-pass filter using the methods described above, and then to transform the low-pass filter to the desired form by a substitution in the transfer function, that is we substitute a function g(s) for s in the low-pass transfer function H(s), so that the new transfer function is H � (s) = H(g(s)). The effect is to modify the filter poles and zeros to produce the desired frequency response characteristic. The critical frequencies used in the design are as shown in Fig. 11. For band-pass and band-stop filters it is convenient to define a center frequency Ωo as the geometric mean of the pass-pand edges, and a band width ΔΩ: � Ωcu Ωcl Ωo = ΔΩ = Ωcu − Ωcl . 17 | H ( jW ) |2 1 R c R s 0 W c W r 0 (ra d /s e c ) W ( a ) H ig h - p a s s | H ( jW ) |2 1 R c R s 0 W 0 rl W W c l (b ) B a n d -p a s s W c u W ru (ra d /s e c ) | H ( jW ) |2 1 R R 0 c s 0 W c l W rl (c ) B a n d -s to p W ru W c u W (ra d /s e c ) Figure 11: Specification tolerance bands for (a) high-pass, (b) band-pass, and (c) band-stop filters. 18 The transformation formulas for a low-pass filter with cut-off frequency Ωc are given in Table 1. The band-pass and band-stop transformations both double the order of the filter, since Ω c1 s Ωc2 Ω2 Low-pass (Ωc ) → High-pass (Ωc ) g(s) = c s s2 + Ω2o Low-pass (Ωc = ΔΩ ) → Band-pass (Ωcl , Ωcu ) g(s) = s sΩ2c Low-pass (Ωc = ΔΩ ) → Band-stop (Ωcl , Ωcu ) g(s) = 2 s + Ω2o Low-pass (Ωc1 ) → Low-pass (Ωc2 ) g(s) = Table 1: Transformation to other filter forms from a prototype low-pass filter with cut-off frequency Ωc . s2 is involved it the transformation. In the above table Ω1 and Ω2 are the edges of the pass/stop-band, and the low-pass filter is designed to have a cut-off frequency Ωc = Ω2 − Ω1 . The above transformations will create an ideal gain characteristic from an ideal lowpass filter. For practical filters, however, the “skirts” of the pass-bands will be a warped representation of the low-pass prototype filter. This does not usually cause problems. Example Show the effect of the low-pass to high-pass conversion by examining the poles and zeros of the transformed first-order filter H(s) = Ωc . s + Ωc The transformation Ω2c /s for s in H(s) gives H � (s) = s . s + Ωc which generates a zero at s = 0 (creating the high-pass action) and leaves the pole unchanged. It is easy to show that the low-pass to high-pass transformation on an nth order all-pole filter will create n zeros at the origin. 3.1 MATLAB Filter Transformation Functions MATLAB’s Signal Processing Toolbox contains filter transformation functions lp2lp(), lp2hp(), lp2bp() and lp2bs() that transform a prototype low-pass filter with a cut-off frequency of 1 rad/s to another low-pass, a high-pass, a band-pass, or a band-stop filter respectively: 19 [NUMT,DENT] [NUMT,DENT] [NUMT,DENT] [NUMT,DENT] = = = = lp2lp(NUM,DEN,Wc) - Low-pass to low-pass. lp2hp(NUM,DEN,Wc) - Low-pass to high-pass. lp2bp(NUM,DEN,Wo,Bw) - Low-pass to band-pass. lp2bs(NUM,DEN,Wo,Bw) - Low-pass to band-stop. The requirement that the prototype low-pass filter have a cut-off frequency of Ωc = 1 rad/s changes the required substitutions in the transfer functions from Table 1 to those given in Table 2. Let the prototype low-pass filter have frequency response H(jΩ) with cut-off Low-pass (Ωc = 1) → Low-pass (Ωc ) Low-pass (Ωc = 1) → High-pass (Ωc ) Low-pass (Ωc = 1) → Band-pass (Ωcl , Ωcu ) Low-pass (Ωc = 1) → Band-stop (Ωcl , Ωcu ) s Ωc Ωc g(s) = s s2 + Ω2o g(s) = ΔΩ s ΔΩ s g(s) = 2 s + Ω2o g(s) = Table 2: Transformation to other filter forms from a prototype low-pass filter with cut-off frequency Ωc = 1 as required by MATLAB functions. frequency Ωc = 1. Then Table 2 allows the stop-band limit Ωr to be computed from Ωr = g(s)|s=jΩr as shown in Table 3. Because band-pass and band-stop filters have two stop-band edges Ωrl and Ωru , we should choose the edge that imposes the the more stringent filter specification, that is the edge that defines the lowest value of Ωr . Ωc Low-pass 1 High-pass 1 Band-pass 1 Band-stop 1 Ωr Ωr Ωc Ωc Ωr � � |Ω2o − Ω2rl | |Ω2o − Ω2ru | min , Δ ΔΩ .Ωru � .Ω Ω rl � ΔΩ .Ωru ΔΩ .Ωrl min , 2 2 2 |Ωo − Ωrl | |Ωo − Ω2ru | Table 3: Mapping of critical frequencies as defined in Fig. 11 to a prototype low-pass filter H(jΩ) with Ωc = 1 for transformation using MATLAB functions lp2lp(), lp2hp(),lp2bp() and lp2bs() Design Procedure: 20 1. Determine the filter specifications and choose band edge frequencies and attenuation values as in Fig. 11. 2. Use Table 3 to define Ωr for the prototype filter. 3. Design the prototype filter using Ωc = 1, Ωr , Rc , and Rs . 4. Transform the prototype filter to the desired form. Example A band-stop filter is to be used to minimize 60 Hz ac interference superimposed on experimental data. A suitable set of specifications was selected as Ωo Ωcu Ωru Ωr1 Rc Rs = = = = = = 2π60 rad/s 2π70 rad/s 2π63 rad/s 2π57 rad/s 1 dB 40 dB (60 (70 (63 (57 Hz) Hz) Hz) Hz) The following MATLAB script was used to design the band-stop filter: % Design a band-stop filter % First define some ctitical frequencies % Passband edges wo=2*pi*60; wcu=2*pi*70; wcl = wo^2/wcu; BW=(wcu-wcl); % Stop band edges wsu=2*pi*63; wsl=2*pi*57; % Rc = 1; Rs = 40; % Determine the stop-band edge in the lp prototype % W1 = BW*wsu/(wo^2-wsu^2); W2 = BW*wsl/(wo^2-wsl^2); Wr = min(abs(W1),abs(W2)); % % design the prototype low-pass filter % [N,Wn] = buttord(1,Wr, Rc, Rs, ’s’); [num,den] = butter(N, Wn, ’s’); % Convert to a band-stop filter 21 [num_stop,den_stop] = lp2bs(num,den,wo,BW); filt = tf(num_stop,den_stop); % % Plot the frequency response magnitude % f=[30:1:90]; [mag,phase]=bode(filt,2*pi*f); plot(f,(squeeze(mag))) grid; xlabel(’Frequency (Hz)’); ylabel(’Response Magnitude’) which produces the following plot 1 0.9 Response Magnitude 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 30 40 50 60 Frequency (Hz) 22 70 80 90 A Sample Low-Pass Filter Design Routines The following MATLAB “.m” routines are included to demonstrate the Butterworth and Chebyshev pole placement methods described in this handout. A.1 Butterworth Low-Pass Filter % *** 2.161 Signal Processing - Continuous and Discrete *** % % lpbutter - Butterworth continuous lowpass filter design based on % derivation in the 2.161 class notes. % % [z,p,k,n] = lpbutter(wc, Rc, ws, Rs) or % [z,p,k,n] = lpbutter(wc, Rc, ws, Rs,’zpk’) returns the zeros, poles, and gain % and order for a Butterworth continuous lowpass filter. % [b,a,n] = lpbutter(wc, Rc, ws, Rs,’tf’) returns the transfer function % numerator and denominator coefficients, and order for a % Butterworth continuous lowpass filter. % Arguments: % wc - passband cut-off frequency (rad/s) % Rc - attenuation at passband cut-off (dB) (Rc > 0) % ws - stopband frequency (rad/s) % Rs - attenuation at stopband (dB) (Rs > 0). % % Author: D. Rowell % Revision: 2.0 9-23-2007 % %--------------------------------------------------------------­ function [varargout] = lpbutter(wc,Rc,ws,Rs,varargin) % % Compute the required order n % epsilon = sqrt(10^(Rc/10)-1); lambda = sqrt(10^(Rs/10)-1); n = ceil(log(lambda/epsilon)/log(ws/wc)); % % There are no zeros in the lpf - create an empty zeros vector: % z = []; % % Create the complex conjugate poles on a circle % p = wc*epsilon^(-1/n)*exp(i*(pi*(1:2:n-1)/(2*n) + pi/2)); p = [p; conj(p)]; p = p(:); % If n is odd, add a real pole if rem(n,2)==1 % n is odd 23 p = [p; -wc*epsilon^(-1/n)]; end % % For unity low freq gain, the k is the product of the negative of the poles. % k = real(prod(-p)); % % Return the appropriate system representation % if (nargin == 5 ) if (varargin{1}(1:2) == ’tf’) sys = zpk(z,p,k); [b,a] = tfdata(sys); varargout(1) = b; varargout(2) = a; varargout(3) = {n}; elseif (varargin{1}(1:3) == ’zpk’) varargout(1)={z}; varargout(2)={p}; varargout(3)={k}; varargout(4)={n}; end else varargout(1)={z}; varargout(2)={p}; varargout(3)={k}; varargout(4)={n}; end 24 A.2 Chebyshev Type 1 Low-pass Filter % *** 2.161 Signal Processing - Continuous and Discrete *** % % lpcheby1 - Chebyshev Type 1 continuous lowpass filter design based on % derivation in the 2.161 class notes. % % [z,p,k,n] = lpcheby1(wc, Rc, ws, Rs) or % [z,p,k,n] = lpcheby1(wc, Rc, ws, Rs,’zpk’) returns the zeros, poles, and gain % and order for a Chebyshev Type 1 continuous lowpass filter. % [b,a,n] = lpcheby1(wc, Rc, ws, Rs,’tf’) returns the transfer function % numerator and denominator coefficients, and order for a % Chebyshev Type 1 continuous lowpass filter. % Arguments: % wc - passband cut-off frequency (rad/s) % Rc - attenuation at passband cut-off (dB) (Rc > 0) % ws - stopband frequency (rad/s) % Rs - attenuation at stopband (dB) (Rs > 0). % % Author: D. Rowell % Revision: 2.0 9-23-2007 % %--------------------------------------------------------------­ function [varargout] = lpcheby1(wc,Rc, ws, Rs,varargin) % % Compute the required order n % epsilon = sqrt(10^(Rc/10)-1); lambda = sqrt(10^(Rs/10)-1); n = ceil(acosh(lambda/epsilon)/acosh(ws/wc)); % % There are no zeros in the lpf - create an empty zeros vector: % z = []; % % Determine the poles - choosing first only complex conjugates in the % left-half plane % alpha = asinh(1/epsilon)/n; p = []; for j=1:floor(n/2) gamma = (2*j - 1)*pi/(2*n); newp = -wc*(sinh(alpha)*sin(gamma)+ i*cosh(alpha)*cos(gamma)); p = [p; newp; conj(newp)]; end % If n is odd - add a real pole if rem(n,2)==1 % n is odd gamma = (2*(floor(n/2)+1) - 1)*pi/(2*n); 25 p = [p; -wc*(sinh(alpha)*sin(gamma))]; end % Compute the gain k = real(prod(-p)); % % If n is even adjust the gain to unity at zero frequency % if ~rem(n,2) % n is even so patch k k = k/sqrt((1 + epsilon^2)); end % % Return the appropriate system representation % if (nargin == 5 ) if (varargin{1}(1:2) == ’tf’) sys = zpk(z,p,k); [b,a] = tfdata(sys); varargout(1) = b; varargout(2) = a; varargout(3) = {n}; elseif (varargin{1}(1:3) == ’zpk’) varargout(1)={z}; varargout(2)={p}; varargout(3)={k}; varargout(4)={n}; end else varargout(1)={z}; varargout(2)={p}; varargout(3)={k}; varargout(4)={n}; end 26 A.3 Chebyshev Type 2 Low-pass Filter % *** 2.161 Signal Processing - Continuous and Discrete *** % % lpcheby2 - Chebyshev Type 2 continuous lowpass filter design based on % derivation in the 2.161 class notes. % % [z,p,k,n] = lpcheby2(wc, Rc, ws, Rs) or % [z,p,k,n] = lpcheby2(wc, Rc, ws, Rs,’zpk’) returns the zeros, poles, and gain % and order for a Chebyshev Type 2 continuous lowpass filter. % [b,a,n] = lpcheby2(wc, Rc, ws, Rs,’tf’) returns the transfer function % numerator and denominator coefficients, and order for a % Chebyshev Type 2 continuous lowpass filter. % Arguments: % wc - passband cut-off frequency (rad/s) % Rc - attenuation at passband cut-off (dB) (Rc > 0) % ws - stopband frequency (rad/s) % Rs - attenuation at stopband (dB) (Rs > 0). % % Author: D. Rowell % Revision: 2.0 9-23-2007 % %--------------------------------------------------------------­ function [varargout] = lpcheby2(wc,Rc, ws, Rs,varargin) % % Compute the required order n % epsilon = sqrt(10^(Rc/10)-1); lambda = sqrt(10^(Rs/10)-1); n = ceil(acosh(lambda/epsilon)/acosh(ws/wc)); wprod = wc*ws; epsilonhat= 1/(epsilon*cos(n*acos(ws/wc))); % % Create the finite (imaginary) zeros % Handle even and odd orders separately to eliminate the zero at infinity % for n odd if (rem(n,2)) m = n - 1; z = cos([1:2:n-2 n+2:2:2*n-1]*pi/(2*n))’; else m = n; z = cos((1:2:2*n-1)*pi/(2*n))’; end z = (z - flipud(z))./2; % Transform back to the s-plane z = i*ws./z; % Organize zeros in complex pairs: j = [1:m/2; m:-1:m/2+1]; 27 z = z(j(:)); % % Determine the poles - choosing first only complex conjugates in the % left-half plane % alpha = asinh(1/epsilonhat)/n; p = []; for j=1:floor(n/2) gamma = (2*j - 1)*pi/(2*n); newp = -wc*(sinh(alpha)*sin(gamma)+ i*cosh(alpha)*cos(gamma)); p = [p; newp; conj(newp)]; end % If n is odd - add a real pole if rem(n,2)==1 % n is odd gamma = (2*(floor(n/2)+1) - 1)*pi/(2*n); p = [p; -wc*(sinh(alpha)*sin(gamma))]; end p=wprod./p; % Compute the gain k = real(prod(-p)/prod(-z)); % % Return the appropriate system representation % if (nargin == 5 ) if (varargin{1}(1:2) == ’tf’) sys = zpk(z,p,k); [b,a] = tfdata(sys); varargout(1) = b; varargout(2) = a; varargout(3) = {n}; elseif (varargin{1}(1:3) == ’zpk’) varargout(1)={z}; varargout(2)={p}; varargout(3)={k}; varargout(4)={n}; end else varargout(1)={z}; varargout(2)={p}; varargout(3)={k}; varargout(4)={n}; end 28 A.4 Demonstration of Filter Design Routines % Test program for low-pass filter design. % Designs lp filters using the example in the class handout % D.Rowell 9/20/2007 % % Define the filter specifications: % cut-off frequency: 10 Hz, max pass-band attenuation: 1 dB wc = 2*pi*10; passbandatten = 1; % stop-band: 20Hz: min stop-band attenuation: 30 dB ws = 2*pi*20; stopbandatten = 30; % Design the filters [z,p,k,nbutt] = lpbutter(wc,passbandatten,ws,stopbandatten); buttfilt=zpk(z,p,k); [z,p,k,ncheby1] = lpcheby1(wc,passbandatten,ws,stopbandatten); cheby1filt=zpk(z,p,k); [z,p,k,ncheby2] = lpcheby2(wc,passbandatten,ws,stopbandatten); cheby2filt=zpk(z,p,k); % plot the three power responses on a linear scale f=[0:.5:30]; w=2*pi*f; [buttmag, buttphase] = bode(buttfilt,w); [cheby1mag, cheby1phase] = bode(cheby1filt,w); [cheby2mag, cheby2phase] = bode(cheby2filt,w); plot(f,squeeze(buttmag).^2, f,squeeze(cheby1mag).^2, f,squeeze(cheby2mag).^2); title(’2.161 Demo Filter Design Software’); xlabel(’Frequency (Hz)’); ylabel(’Power Response’); grid; 2.161 Demo Filter Design Software 1 Chebyshev Type1, N=4 Power Response 0.8 0.6 Chebyshev Type2, N=4 0.4 Butterworth, N=6 0.2 0 0 5 10 15 Frequency (Hz) 29 20 25 30