calculating the coefficients of a butterworth iir filter

Hallo everyone!

I wanted to ask your advice on the following problem:

i have to generate a butterworth iir bandpass filter with cutoff
frequencies of 0.008 and 0.05 Hertz. This will be used in filtering a time
signal sampled at 1 Hz (so the sampling frequency is high enough to work
with the time signal) 

My initial problem is to calculate the coefficients. Ive thought of a few
ways to calculate them which could be (if im not wrong):

my input parameters are : Filter Order and Gain. The programming language
used is c++ with it++ as a signal processing library.

   1. generate a frequency vector (size : Filter Order + 1) running from 0
to  1.0 and a corresponding Magnitude vector and calculate the coefficients
using the Yule Walker equations. this however is a bit inaccurate because i
have to modify the magnitude vector using trial and error.

   2. taking the coefficients from a table, which is easy and would work,
but others have calculated these, so I would like to learn it.

  3. Billinear transform. But this appears somewhat complicated for a
filter of higher order and looking at my extremely narrow passband, i
wonder if this is the right way. :-/

ive looked at the sourcecode for some java applets on the net that
calculate filter coefficients, but found it quite hard to follow. 

so my question is: what would be an effective method for calculating the
butterworth filter coefficients especially in c++?

many thanks for your replies in advance!

dilpreet

dilpreet06 skrev:
> Hallo everyone!
>
> I wanted to ask your advice on the following problem:
>
> i have to generate a butterworth iir bandpass filter with cutoff
> frequencies of 0.008 and 0.05 Hertz. This will be used in filtering a time
> signal sampled at 1 Hz (so the sampling frequency is high enough to work
> with the time signal)
>
> My initial problem is to calculate the coefficients. Ive thought of a few
> ways to calculate them which could be (if im not wrong):
>
> my input parameters are : Filter Order and Gain.

The most effective input comprise sampling frequency, stop band(s)
and pass-band(s) cut-off frequencies, stop band attenuation and
passband ripple.

The gain is merely a sacling factor and the filter order is computed
along the way based on the above spec.

> The programming language
> used is c++ with it++ as a signal processing library.
>
>    1. generate a frequency vector (size : Filter Order + 1) running from 0
> to  1.0 and a corresponding Magnitude vector and calculate the coefficients
> using the Yule Walker equations. this however is a bit inaccurate because i
> have to modify the magnitude vector using trial and error.

It is possible to design filters that way, but you are asking
specifically for Butterworth IIRs. There are far more efficient
ways to design those.

>    2. taking the coefficients from a table, which is easy and would work,
> but others have calculated these, so I would like to learn it.

Use this approach if you need immediate results. If you have
more time on your hands, find a decent text on filter design that
contains the recipe for how to do this. One such text is suggested
below.

>   3. Billinear transform. But this appears somewhat complicated for a
> filter of higher order and looking at my extremely narrow passband, i
> wonder if this is the right way. :-/

This is the way to go, if you wnat to learn.

The key is to see that higher-order IIRs are implemented as a
sequence of 2-pole primitives, the "biquads." You only have to
work out the bilinear transform for orders 1 and 2. When you
design your 73-order IIR filter you represent it as 36 biquads
and one 1st order section, and apply the 1st and 2nd order
BLTs to each of those.

> ive looked at the sourcecode for some java applets on the net that
> calculate filter coefficients, but found it quite hard to follow.

Don't bother. You have to dig deep into the textbook literature to
find these recipes. No one have written textbooks on these
matters for more than 30 years. The one recent book I know of,
that contains the recipes, is

http://www.amazon.com/gp/product/0071454241/sr=8-1/qid=1154767557/ref=sr_1_1/103-3640117-4140620?ie=UTF8

It contains recipes for Butterworths, Chebychevs and elliptic
filters, among others. It also contain frequency tranforms and
shows how to apply pre-warping to fine-tune the designs.

This is a new version of some text that first was written in the
60s. It is good on filter design recipes, but some details in
the notation are at odds with what has emerged as the standard
in DSP these days. Usually, the coefficients in the numerator
polynomial in transfer functions H(w) are designated b_i,
and the coefficients in the denominator polynomial are
denoted as a_i. In this book, they are the oposite.

Not a big deal, but maybe worth keeping a note of if you
read other DSP texts as well, and want to compare
approaches.

> so my question is: what would be an effective method for calculating the
> butterworth filter coefficients especially in c++?

It depends on what you mean by "effective". Designing the
filter is not very computationally expensive; they are done
with in a matter of milliseconds. The main problem in this
application is to keep track of all the details while preserving
the sanity of whoever maintains the code.

Starting with the spec as indicated above (normalized corner
frequencies and attenuations/ripple), one has to do a lot
of tasks that are trivial per se, but that need to be fitted into
a rather complex picture:

- Transform the spec to an LP prototype. This step depends on
  if one wants an LP, HP BP or BS filter.
- Pre-warp the prototype parameters from discrete-time to
  continuous time.
- Compute the order of the filter. This step depends on type
  of filter (Butterworth, Chebyshev, Elliptic).
- Compute the coefficients on the biquads.This step depends
  on type of filter (Butterworth, Chebyshev, Elliptic).
- Apply the BLT to the biquads
- Transform the dicrete-time prototype filter to the desired
  discrete-time filter. This step depends on if one wants an
  LP, HP BP or BS filter.

In some of these stages one might want to use polynomial
root solvers to separate 4-pole sections into two biquads.

There is lots of book-keeping to be done here. Implement
this wrong and you immediately have a maintenance night-mare
on your hands. I used this very task as an oportunity to program
"real" OO code with C++. I managed to keep track of all these
details without using a single "case switch" statement. Everything
is taken care of only using classes and argument type matching
on method calls.

It can be done, but the starting point needs to be that you
understand the nuts'n bolts of designing filters.
If you want to learn how do design filters, design both
a Butterworth and a Chebyshev filter according to the
spec above (and a LP spec and a HP spec), so you get a
feel for what is similar and what is different. Don't get mixed
up with program organization until you have done that; the
computational routines are almost trivial.

When you get a sense for the big picture, start playing
around with ways of merging all this into a C++/Java
setting. At this stage you ought to have all the computational
details available, so you can focus your attention to
program organization.

Rune

Rune, thanks for the depthful and carefully explained thread! 

You've more or less confirmed what i suspected earlier. I shall use the
tables for now, but try out the biquad and billinear transform in the
future when i have more time. 

Ive seen some of your (and other peoples') replies to other threads here.
Thanks for sharing your wealth of info with others. Knowledge is power.

Peace.

dilpreet06 wrote:
...
> Knowledge is power.

Others here would counter: "Imagination is more important than
knowledge".

:-)

dilpreet06 wrote:

> Rune, thanks for the depthful and carefully explained thread!

Y're welcome.

> You've more or less confirmed what i suspected earlier. I shall use the
> tables for now, but try out the biquad and billinear transform in the
> future when i have more time.

That's the trick, isn't it: Use the right approach at the right time.
Get results when needed, expand your knowledge base when you
have the oportunity.

> Ive seen some of your (and other peoples') replies to other threads here.
> Thanks for sharing your wealth of info with others.

I've learned more DSP here on comp.dsp in the last few years, than
during a decade in universities and R&D organizations. My experience
is that one haven't learned the material well until one is able to
explain it to others in their own terms.

> Knowledge is power.

Opinions may be divided on the matter...

Rune