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


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

Jerry Avins <jya@ieee.org> writes:

> Error. Typesetter's or van Valkenburg's we'll never know, so we can't blame the
> proofreader.

"Mac" van Valkenburg has been dead for a few years now, so we may 
never know.  I will check my copy at home for any penciled in marks 
from his class.

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),
but if you're having trouble with your textbook about it, then maybe it
isn't.

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?

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).  by fudging those frequencies a
little, you can make sure that you compare apples to apples.

-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."


robert bristow-johnson wrote:
> in article 1118556798.292050.210830@g47g2000cwa.googlegroups.com, Rune
> Allnor at allnor@tele.ntnu.no wrote on 06/12/2005 02:13:
>
> > I just wanted the coefficients that pop out from a given spec and
> > algorithm,
>
> sounds like closed form.
>
> > to check that the code I implement myself actually works.
> > I believe my code for the Cheb 1 case works (I haven't tested it
> > with pre-warping and frequency transforms) and I thought I would have
> > a try at the Cheb 2 case today.
>
> i think i have closed form solutions to the s-plane poles as zeros of the
> Tchebyshev (type 1) and Inverse Chebyshev (type 2) low-pass filters.  i used
> to have formulae for the complete digital design (z plane poles/zeros) but
> it's left at a previous employer.
>
> i might be able to rederive it.  do you want me to post it?

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.

Rune

in article 1118556798.292050.210830@g47g2000cwa.googlegroups.com, Rune
Allnor at allnor@tele.ntnu.no wrote on 06/12/2005 02:13:

> I just wanted the coefficients that pop out from a given spec and
> algorithm,

sounds like closed form.

> to check that the code I implement myself actually works.
> I believe my code for the Cheb 1 case works (I haven't tested it
> with pre-warping and frequency transforms) and I thought I would have
> a try at the Cheb 2 case today.

i think i have closed form solutions to the s-plane poles as zeros of the
Tchebyshev (type 1) and Inverse Chebyshev (type 2) low-pass filters.  i used
to have formulae for the complete digital design (z plane poles/zeros) but
it's left at a previous employer.

i might be able to rederive it.  do you want me to post it?

-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

robert bristow-johnson wrote:
> Rune,
>
> i haven't figgered out the gist of the thread.  do you want closed form
> expressions for the poles/zeros or coefficients for Tchebyshev type 1 or 2
> filters?  (mapping to z domian with BLT is not a big deal, unless there is
> prewarping.)

I just wanted the coefficients that pop out from a given spec and
algorithm, to check that the code I implement myself actually works.
I believe my code for the Cheb 1 case works (I haven't tested it
with pre-warping and frequency transforms) and I thought I would have
a try at the Cheb 2 case today.

Naively, I thought it would be a piece of cake to implement a recipe
for Cheb 1 and Cheb 2 filters as found in e.g. Proakis & Manolakis,
but I was wrong.

While working on the Butterworth case, I found two typos in the
expression for the lowpass -> bandstop transform. First, the sign
of the substitute for z^{-1} is wrong, it should be the opposite of
the sign in the lowpass -> bandpass case. Second, there is an
internal variable on the form

    q = 2*alpha /(K-1)

(sorry, I don't have my copy of P&M here, so the details in the
formula above may be wrong, and I can't tell you the table number
and page where the typos are) that has the wrong sign. The fact
that these typos exist in the 3rd edition indicate to me that no
one have used the P&M book as basis for an implementation of these
transforms. And so there may be further undetected typos in the
examples. My confidence in P&M have faded a bit over the last few
days...

With trusted reference data available I can track down my own typos,
misunderstandings and blunders. The recipes and tables I found in
van Valkenburgs book were invaluable. I had working code half an
hour after I sat down with that book. The algorithm was explained
very well in the text, succintly distilled in a design recipe,
and the tables of pole locations in prototype filters were very
useful during debugging. For comparision, I spent a whole day on
the P&M Cheb 1 recipe and got only garbage. 

Rune