Forums

IIR filter reference designs

Started by Rune Allnor June 9, 2005

robert bristow-johnson wrote:
> in article 1118655971.796709.5760@g47g2000cwa.googlegroups.com, Rune Allnor > at allnor@tele.ntnu.no wrote on 06/13/2005 05:46: > > > No, I don't think so. The sentence "it's left with a previous emplyer" > > is a turn-off; I don't want to get you into trouble. > > i understand. it's really just textbook sorta stuff (BLT, nothing more),
Somebody once told me to 'never "never mind"', and that goes double when attorneys, patents and commercial property is involved, no matter how trivial the stuff might turn out to be.
> but if you're having trouble with your textbook about it, then maybe it > isn't.
Well, the textbooks I have, are Oppenheim & Schafer's 1975 and 1999 editions, Proakis & Manolakis 1996 edition, and some old stuff by Jackson. All of these are DSP texts that mention IIR filter design by means of the BLT merely in passing. Some (O&S, P&M) work through a Butterworth LP design example, but that's all. P&M mentions some formulae for the Cheb 1 case that I can't get to work. Van Valkenburg's 1982 book is a text particularly aimed at analog filter design, which means he takes stuff for granted that I, who was trained in college with DSP as the main material and analog electronics as the esoteric stuff, may have missed. Or at least forgotten. Whether that is a weakness with the text or with my preparation, is probably a matter of definition.
> i have this 30 year old book, "Active Network Design" by Lindquist that has > closed form s-plane formulae. nobody can get after me for that. do you want > those or does your Van Valkenberg book have those?
If I could get those for the Cheb 1 and 2 filters, I might perhaps be able to make sense of van Valkenburg... To be fair, van Valkenburg's Cheb 1 material was brilliant. My filters worked 30 minutes after I had read his chapter on Cheb 1 filters. Now I have trouble with the Cheb 2 filters.
> one thing to do, to compare both types of Tchebychev to each other and to > Butterworth, is to fix the -3 dB corner frequency to the same place. there > is a simple expression that relates the -3 dB frequency to the passband edge > (Tcheb type 1) or stopband edge (type 2).
Yep, that much I've got from van Valkenburg. I've even been able to reproduce some pole and zero locations van Valkenburg list for the Cheb 2 filter. The hard part is to understand where to put the numbers in the formulae: I have a vague recollection from Control Systems 101 (which was the last time I used S plane formulae) that causal, stable continuous-time systems have an s transform of the type b_0 + b_1s + ... + b_(n-1)s^(n-1) H(s) = ------------------------------------------ a_0 + a_1s + ... + a_(n-1)s^(n-1) + s^n i.e. that the denominator must be a polynomial in s that is of strictly higher degree than the numerator. For the Cheb 2 systems, I find that the zeros come in pairs, i.e. that the 2nd order "filter primitives" become b_0 + b_2 s^2 H'(s) = ----------------------- a_0 + a_1 s + a_2 s^2 which contradicts my very vague recollection from an intro class of some 12-15 years ago. Oh well, I suspect I have revealed enough of my ignorance in such matters for one post... ;)
> by fudging those frequencies a > little, you can make sure that you compare apples to apples.
That's always a good starting point. Rune
in article 1118686950.000754.48610@g49g2000cwa.googlegroups.com, Rune Allnor
at allnor@tele.ntnu.no wrote on 06/13/2005 14:22:

> I have a vague recollection from Control Systems 101 > (which was the last time I used S plane formulae) that causal, stable > continuous-time systems have an s transform of the type > > b_0 + b_1s + ... + b_(n-1)s^(n-1) > H(s) = ------------------------------------------ > a_0 + a_1s + ... + a_(n-1)s^(n-1) + s^n > > i.e. that the denominator must be a polynomial in s that is of strictly > higher degree than the numerator.
no, they can be equal order b_0 + b_1 s + ... + b_(n-1) s^(n-1) + b_n s^n H(s) = ------------------------------------------------ a_0 + a_1 s + ... + a_(n-1) s^(n-1) + s^n which is equal to c_0 + c_1 s + ... + c_(n-1) s^(n-1) H(s) = b_n + ------------------------------------------------ a_0 + a_1 s + ... + a_(n-1) s^(n-1) + s^n where c_k = b_n*a_k - b_k
> For the Cheb 2 systems, I find that > the zeros come in pairs, i.e. that the 2nd order "filter primitives" > become > > b_0 + b_2 s^2 > H'(s) = ----------------------- > a_0 + a_1 s + a_2 s^2 > > which contradicts my very vague recollection from an intro class of > some 12-15 years ago.
it will for even order, anyway. for odd order, the number of poles exceeds the number of zeros by 1.
> Oh well, I suspect I have revealed enough of my > ignorance in such matters for one post... ;)
don't sweat it. i found that lot'sa current DSP people know little to nothing about classical s-plane filter design. i don't know diddley about half or more of the stuff posted here. here's what i see in the textbooks about Tchebyshev Type I and II filters: Tchebyshev polynomials: { (-1)^n * cosh(n*arccosh(-x)) x <= -1 { T_n(x) = { cos(n*arccos(x)) -1 <= x <= 1 { { cosh(n*arccosh(x)) 1 <= x even though it doesn't look like an n^th order polynomial, it is. it's obvious that T_0(x) = 1 T_1(x) = x and you can work out that T_(n+1)(x) = 2*x*T_n(x) - T_(n-1)(x) even order Tchebyshevs have even symmetry, odd order has odd symmetry. note that whatever order, |T_n(x)| <= 1 for |x| <= 1 and that, given higher orders, n, that T_n(x) takes off like a raped ape after |x| exceeds 1 (no other such constrained n^th order polynomial will exceed T_n(x) for x >= 1). that differential of magnitude (with constraint of magnitude for |x| <= 1) is what is used to make a "filter" to effectively discriminate between what is kept and what gets reduced. Type I LPF: |H(jw)|^2 = 1/(1 + (e*T_n(w/w_e))^2) which is a degenerate case of |H(s)|^2 = 1/(1 + (e*T_n(s/j*w_e))^2) it can be easily shown that in the passband (|w| < w_e) that the magnitude ripples as much as dB_ripple = 10*log10(1 + e^2) and that the -3 dB frequency is w0 = w_e * cosh( 1/n * arccosh(1/e) ) where e <= 1 or dB_ripple <= 3.01 dB . it's also easy to see there are no analog zeros to H(s) (there's just a 1 in the numerator of |H(s)|^2). actually there are zeros, but they're all at infinity (this distinction is salient when mapping them to z = -1 using the BLT). what isn't easy to show (but it's possible) is that the poles (whatever makes the denominator of |H(s_p)|^2 go to zero) are n even, 0 <= k <= n/2 - 1 s_p/w_e = -cos(pi*(k+0.5)/n) * sinh(1/n*arcsinh(1/e)) +/- j*sin(pi*(k+0.5)/n) * cosh(1/n*arcsinh(1/e)) n odd, 0 <= k <= (n-1)/2 s_p/w_e = -cos(pi*k/n) * sinh(1/n*arcsinh(1/e)) +/- j*sin(pi*k/n) * cosh(1/n*arcsinh(1/e)) this is all textbook. so solve for e, substitute w0 in for w_e, do the bilinear transforms on your poles and zeros and you should have something closed form. for the Type II (Inverse Tchebyshev) LPF: |H(jw)|^2 = 1/(1 + 1/(e*T_n(w_e/w))^2) = (e*T_n(w_e/w))^2 / (1 + (e*T_n(w_e/w))^2) which is a degenerate case of |H(s)|^2 = = (e*T_n(j*w_e/s))^2 / (1 + (e*T_n(j*w_e/s))^2) this time the stopband edge is w_e. it can be easily shown that in the stopband (|w| > w_e) that the maximum gain is dB_stopband = 10*log10(1/(1 + 1/e^2)) = 10*log10(e^2/(1 + e^2)) and the -3 dB frequency is w0 = w_e / cosh( 1/n * arccosh(1/e) ) where e <= 1 or dB_stopband <= -3.01 dB . the zeros are whatever makes the numerator of |H(s)|^2 go to zero and the poles are whatever makes the denominator of |H(s)|^2 go to zero. zeros (even or odd n): s_z/w_e = +/- j/cos(pi*(k+0.5)/n) ( when k = (n-1)/2, then there is a zero at infinity and the gain of the LPF will die out to 0 at very high frequencies. for even n, the gain is sqrt(e^2/(1 + e^2)). ) poles: for n even, 0 <= k <= n/2 - 1 s_p/w_e = [ -cos(pi*(k+0.5)/n) * cosh(1/n*arcsinh(1/e)) +/- j*sin(pi*(k+0.5)/n) * sinh(1/n*arcsinh(1/e)) ]^(-1) for n odd, 0 <= k <= (n-1)/2 s_p/w_e = [ -cos(pi*k/n) * cosh(1/n*arcsinh(1/e)) +/- j*sin(pi*k/n) * sinh(1/n*arcsinh(1/e)) ]^(-1) with the same song and dance as before, you get your closed form solution in the z-domain. the cool thing is, if you fix w0 to be constant and let your "e" go to zero, they both go to the n^th order Butterworth. i might have mistranscribed this (i tried to simplify some expressions and symbols). if anyone sees mistakes, please correct them here. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."

robert bristow-johnson wrote:
> in article 1118686950.000754.48610@g49g2000cwa.googlegroups.com, Rune Allnor > at allnor@tele.ntnu.no wrote on 06/13/2005 14:22: > > > I have a vague recollection from Control Systems 101 > > (which was the last time I used S plane formulae) that causal, stable > > continuous-time systems have an s transform of the type > > > > b_0 + b_1s + ... + b_(n-1)s^(n-1) > > H(s) = ------------------------------------------ > > a_0 + a_1s + ... + a_(n-1)s^(n-1) + s^n > > > > i.e. that the denominator must be a polynomial in s that is of strictly > > higher degree than the numerator. > > no, they can be equal order > > b_0 + b_1 s + ... + b_(n-1) s^(n-1) + b_n s^n > H(s) = ------------------------------------------------ > a_0 + a_1 s + ... + a_(n-1) s^(n-1) + s^n > > which is equal to > > c_0 + c_1 s + ... + c_(n-1) s^(n-1) > H(s) = b_n + ------------------------------------------------ > a_0 + a_1 s + ... + a_(n-1) s^(n-1) + s^n > > > where c_k = b_n*a_k - b_k > > > For the Cheb 2 systems, I find that > > the zeros come in pairs, i.e. that the 2nd order "filter primitives" > > become > > > > b_0 + b_2 s^2 > > H'(s) = ----------------------- > > a_0 + a_1 s + a_2 s^2 > > > > which contradicts my very vague recollection from an intro class of > > some 12-15 years ago. > > it will for even order, anyway. for odd order, the number of poles exceeds > the number of zeros by 1. > > > Oh well, I suspect I have revealed enough of my > > ignorance in such matters for one post... ;) > > don't sweat it. i found that lot'sa current DSP people know little to > nothing about classical s-plane filter design. i don't know diddley about > half or more of the stuff posted here. > > here's what i see in the textbooks about Tchebyshev Type I and II filters:
[... lots of good stuff snipped ...] Hi RBJ, and thanks for the elaborate post. Your post gave basically the same message I got from van Valkenburg's book, so I got a bit more confidence that I actually had the big picture; that there were no hidden insights I had missed. So I took a closer look at my code. It turned out that I had misinterpreted the term "zero at infinity" and introduced a b_1*s term in the first-order section of the filter. Once that was corrected (and when I trusted that the pair of zeros actually give a b_2*s^2 term in the 2nd order sections), everything came together. Thanks for your help. I'd probably spent days or weeks sorting this out without your assistance. Rune
in article 1118921914.707555.55380@o13g2000cwo.googlegroups.com, Rune Allnor
at allnor@tele.ntnu.no wrote on 06/16/2005 07:38:

> robert bristow-johnson wrote: > in article 1118686950.000754.48610@g49g2000cwa.googlegroups.com, Rune Allnor > at allnor@tele.ntnu.no wrote on 06/13/2005 14:22: > > > I have a vague recollection from Control Systems 101 > > (which was the last time I used S plane formulae) that causal, stable > > continuous-time systems have an s transform of the type > > > > b_0 + b_1s + ... + b_(n-1)s^(n-1) > > H(s) = ------------------------------------------ > > a_0 + a_1s + ... + a_(n-1)s^(n-1) + s^n > > > > i.e. that the denominator must be a polynomial in s that is of strictly > > higher degree than the numerator. > > no, they can be equal order > > b_0 + b_1 s + ... + b_(n-1) s^(n-1) + b_n s^n > H(s) = ------------------------------------------------ > a_0 + a_1 s + ... + a_(n-1) s^(n-1) + s^n > > which is equal to > > c_0 + c_1 s + ... + c_(n-1) s^(n-1) > H(s) = b_n + ------------------------------------------------ > a_0 + a_1 s + ... + a_(n-1) s^(n-1) + s^n > > > where c_k = b_n*a_k - b_k
screwup. i think it's c_k = b_k - b_n*a_k
> > For the Cheb 2 systems, I find that > > the zeros come in pairs, i.e. that the 2nd order "filter primitives" > > become > > > > b_0 + b_2 s^2 > > H'(s) = ----------------------- > > a_0 + a_1 s + a_2 s^2 > > > > which contradicts my very vague recollection from an intro class of > > some 12-15 years ago. > > it will for even order, anyway. for odd order, the number of poles exceeds > the number of zeros by 1.
...
> [... lots of good stuff snipped ...] > > Hi RBJ, and thanks for the elaborate post. > > Your post gave basically the same message I got from van Valkenburg's > book, so I got a bit more confidence that I actually had the big > picture; that there were no hidden insights I had missed. > > So I took a closer look at my code. It turned out that I had > misinterpreted the term "zero at infinity" and introduced a b_1*s term > in the first-order section of the filter.
yeah, instead of thinking of it as (s-q1)(s-q2)...(s-qn) H(s) = G ----------------------- (s-p1)(s-p2)...(s-pn) think of it as (1-s/q1)(1-s/q2)...(1-s/qn) H(s) = H(0) ----------------------------- (1-s/p1)(1-s/p2)...(1-s/pn) then the concept of "zero at infinity" or "pole at infinity" has more tractable meaning. either way, under the BLT, they get mapped to z = -1.
> Once that was corrected (and > when I trusted that the pair of zeros actually give a b_2*s^2 term in > the 2nd order sections), everything came together.
the zeros of the Tchebyshev Type II (or "Inverse Tchebyshev") low-pass filter are on the jw axis. (except for that one zero at infinity for an odd-order Type II LPF.)
> Thanks for your help. I'd probably spent days or weeks sorting this > out without your assistance.
now, just plug and chug. the BLT maps poles and zeros just like it maps any other s to the z-plane. for each conjugate pair of poles/zeros, you get coefs for a biquad. lemme know (email if you want) if there are any holes left. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."