IIR filter reference designs

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."