# Help with Butterworth Filter

Started by March 29, 2004
```Hi All,

I've been plowing through pages of theory in dsp books trying to work
out how to create a butterworth filter.  They all seem to use
different symbols and many different forms of the same equation.  Lots
of theory but very little of practical examples.

Basically I have a series of 1500 points taken at even intervals.
Let's call it f(t).  Eg

t    f(t)
1    3412
2    3432
3    3300
4    3433
......... etc

Basically I'm looking to numerically produce g(t) which is f(t) with
all frequencies higher than 0.01 (seconds/data point) filtered out.

Keeping it simple - lets say I want a second order (ie 2 pole)
Butterworth Filter this gives (from a table in a book)

a0= 8.663387E-04
a1= 1.732678E-03 b1= 1.919129E+00
a2= 8.663387E-04 b2= -9.225943E-01

Question is -  what is the formula that I plug these five constants
into to give me g(t), ie the lowpass filter output at a point t???
```
```The constants seem scaled already, should be normalized to 1.
(These are locations of zeros and poles in the complex plane for the
filter).

"Stewart" <windsurfing_stew@yahoo.com.au> wrote in message
> Hi All,
>
> I've been plowing through pages of theory in dsp books trying to work
> out how to create a butterworth filter.  They all seem to use
> different symbols and many different forms of the same equation.  Lots
> of theory but very little of practical examples.
>
> Basically I have a series of 1500 points taken at even intervals.
> Let's call it f(t).  Eg
>
> t    f(t)
> 1    3412
> 2    3432
> 3    3300
> 4    3433
> ......... etc
>
> Basically I'm looking to numerically produce g(t) which is f(t) with
> all frequencies higher than 0.01 (seconds/data point) filtered out.
>
> Keeping it simple - lets say I want a second order (ie 2 pole)
> Butterworth Filter this gives (from a table in a book)
>
> a0= 8.663387E-04
> a1= 1.732678E-03 b1= 1.919129E+00
> a2= 8.663387E-04 b2= -9.225943E-01
>
> Question is -  what is the formula that I plug these five constants
> into to give me g(t), ie the lowpass filter output at a point t???

```
```Stewart wrote:

> Hi All,
>
> I've been plowing through pages of theory in dsp books trying to work
> out how to create a butterworth filter.  They all seem to use
> different symbols and many different forms of the same equation.  Lots
> of theory but very little of practical examples.
>
> Basically I have a series of 1500 points taken at even intervals.
> Let's call it f(t).  Eg
>
> t    f(t)
> 1    3412
> 2    3432
> 3    3300
> 4    3433
> ......... etc
>
> Basically I'm looking to numerically produce g(t) which is f(t) with
> all frequencies higher than 0.01 (seconds/data point) filtered out.
>
> Keeping it simple - lets say I want a second order (ie 2 pole)
> Butterworth Filter this gives (from a table in a book)
>
> a0= 8.663387E-04
> a1= 1.732678E-03 b1= 1.919129E+00
> a2= 8.663387E-04 b2= -9.225943E-01
>
> Question is -  what is the formula that I plug these five constants
> into to give me g(t), ie the lowpass filter output at a point t???

First observation:  that doesn't look right.  The transfer function
should be normalized (see OP) so that the highest-order term in the
denominator is one:

b_2 z^2 + b_1 z
H(z) = -------------------
z^2 + a_1 z + a_0

Then g(t) = a_1 * g(t-1) + a_0 * g(t-2) + b_2 * f(t) + b_1 * f(t-1).
Note that this is recursive, and that g(t) never goes to zero once f(t)
has been wiggled.  That makes it an infinite impulse response filter.

Second observation:  Why in the heck are you messing with IIR filters if
you have a DSP?  Unless you're designing a closed-loop control system
where delay is critical you can get much better attenuation performance
(and better phase matching) by using a FIR (finite impulse response
filter).  And if you _are_ designing a closed-loop control system you
want to avoid the Butterworth and it's ilk like the plague.  Not only
that, but the whole motivation for a DSP having a blazing fast
multiply-accumulate instruction is so you can implement FIR filters
efficiently.

--

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com
```
```Hii Stewart,

look at
http://www.mathworks.com/access/helpdesk/help/toolbox/signal/butter.html
for general information about implementing a butterworth filter.
Then check
http://www.mathworks.com/access/helpdesk/help/toolbox/signal/dfilt.df1.html
and you'll see a graphical representation which explains the use of
the b[] and a[] coefficients.

Be aware: some books name the coefficients and or their signs
differently.
I guess that the source of your coefficients uses a and b the other
way round.

Hope this gives some idea how to go on.

Bernhard

```
```Tim Wescott <tim@wescottnospamdesign.com> wrote in message
news:106i01ldiirmq82@corp.supernews.com...
> Stewart wrote:
>
> > Hi All,
> >
> > I've been plowing through pages of theory in dsp books trying to work
> > out how to create a butterworth filter.  They all seem to use
> > different symbols and many different forms of the same equation.  Lots
> > of theory but very little of practical examples.
> >
> > Basically I have a series of 1500 points taken at even intervals.
> > Let's call it f(t).  Eg
> >
> > t    f(t)
> > 1    3412
> > 2    3432
> > 3    3300
> > 4    3433
> > ......... etc
> >
> > Basically I'm looking to numerically produce g(t) which is f(t) with
> > all frequencies higher than 0.01 (seconds/data point) filtered out.
> >
> > Keeping it simple - lets say I want a second order (ie 2 pole)
> > Butterworth Filter this gives (from a table in a book)
> >
> > a0= 8.663387E-04
> > a1= 1.732678E-03 b1= 1.919129E+00
> > a2= 8.663387E-04 b2= -9.225943E-01
> >
> > Question is -  what is the formula that I plug these five constants
> > into to give me g(t), ie the lowpass filter output at a point t???
>
> First observation:  that doesn't look right.  The transfer function
> should be normalized (see OP) so that the highest-order term in the
> denominator is one:
>
>           b_2 z^2 + b_1 z
> H(z) = -------------------
>          z^2 + a_1 z + a_0
>
> Then g(t) = a_1 * g(t-1) + a_0 * g(t-2) + b_2 * f(t) + b_1 * f(t-1).
> Note that this is recursive, and that g(t) never goes to zero once f(t)
> has been wiggled.  That makes it an infinite impulse response filter.

info.  Note, the author uses beta and alpha instead of a and b.
Here's some C code if that is of use to you:

> Second observation:  Why in the heck are you messing with IIR filters if
> you have a DSP?  Unless you're designing a closed-loop control system
> where delay is critical you can get much better attenuation performance
> (and better phase matching) by using a FIR (finite impulse response
> filter).  And if you _are_ designing a closed-loop control system you
> want to avoid the Butterworth and it's ilk like the plague.  Not only
> that, but the whole motivation for a DSP having a blazing fast
> multiply-accumulate instruction is so you can implement FIR filters
> efficiently.

There are plenty of applications where using IIR filters makes sense on a DSP.
IIR filters can be much more efficient than FIRs, especially when the cut-off
frequency is low compared to the sample rate.  They also have a simplicity and
ease of parameterization that FIRs don't if you need to calculate filter
coefficients "on the fly".

```
```Jon Harris wrote:

> Tim Wescott <tim@wescottnospamdesign.com> wrote in message
> news:106i01ldiirmq82@corp.supernews.com...
>
>>Stewart wrote:
>>
>>
>>>Hi All,
>>>
>>>I've been plowing through pages of theory in dsp books trying to work
>>>out how to create a butterworth filter.  They all seem to use
>>>different symbols and many different forms of the same equation.  Lots
>>>of theory but very little of practical examples.
>>>
>>>Basically I have a series of 1500 points taken at even intervals.
>>>Let's call it f(t).  Eg
>>>
>>>t    f(t)
>>>1    3412
>>>2    3432
>>>3    3300
>>>4    3433
>>>......... etc
>>>
>>>Basically I'm looking to numerically produce g(t) which is f(t) with
>>>all frequencies higher than 0.01 (seconds/data point) filtered out.
>>>
>>>Keeping it simple - lets say I want a second order (ie 2 pole)
>>>Butterworth Filter this gives (from a table in a book)
>>>
>>>a0= 8.663387E-04
>>>a1= 1.732678E-03 b1= 1.919129E+00
>>>a2= 8.663387E-04 b2= -9.225943E-01
>>>
>>>Question is -  what is the formula that I plug these five constants
>>>into to give me g(t), ie the lowpass filter output at a point t???
>>
>>First observation:  that doesn't look right.  The transfer function
>>should be normalized (see OP) so that the highest-order term in the
>>denominator is one:
>>
>>          b_2 z^2 + b_1 z
>>H(z) = -------------------
>>         z^2 + a_1 z + a_0
>>
>>Then g(t) = a_1 * g(t-1) + a_0 * g(t-2) + b_2 * f(t) + b_1 * f(t-1).
>>Note that this is recursive, and that g(t) never goes to zero once f(t)
>>has been wiggled.  That makes it an infinite impulse response filter.
>
>
> info.  Note, the author uses beta and alpha instead of a and b.
> Here's some C code if that is of use to you:
>
>
>>Second observation:  Why in the heck are you messing with IIR filters if
>>you have a DSP?  Unless you're designing a closed-loop control system
>>where delay is critical you can get much better attenuation performance
>>(and better phase matching) by using a FIR (finite impulse response
>>filter).  And if you _are_ designing a closed-loop control system you
>>want to avoid the Butterworth and it's ilk like the plague.  Not only
>>that, but the whole motivation for a DSP having a blazing fast
>>multiply-accumulate instruction is so you can implement FIR filters
>>efficiently.
>
>
> There are plenty of applications where using IIR filters makes sense on a DSP.
> IIR filters can be much more efficient than FIRs, especially when the cut-off
> frequency is low compared to the sample rate.  They also have a simplicity and
> ease of parameterization that FIRs don't if you need to calculate filter
> coefficients "on the fly".
>
>

True, but this guy is obviously a beginner, and his comments make it
sound like he's looking for a brick-wall filter.  Both of these indicate
that his life would be simplified with a FIR.

--

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com
```
```Am 29 Mar 2004 19:30:00 -0800 schrieb Stewart:

> Hi All,
>
> I've been plowing through pages of theory in dsp books trying to work
> out how to create a butterworth filter.  They all seem to use
> different symbols and many different forms of the same equation.  Lots
> of theory but very little of practical examples.
>
> Basically I have a series of 1500 points taken at even intervals.
> Let's call it f(t).  Eg
>
> t    f(t)
> 1    3412
> 2    3432
> 3    3300
> 4    3433
> ......... etc
>
> Basically I'm looking to numerically produce g(t) which is f(t) with
> all frequencies higher than 0.01 (seconds/data point) filtered out.
>
> Keeping it simple - lets say I want a second order (ie 2 pole)
> Butterworth Filter this gives (from a table in a book)
>
> a0= 8.663387E-04
> a1= 1.732678E-03 b1= 1.919129E+00
> a2= 8.663387E-04 b2= -9.225943E-01
>
> Question is -  what is the formula that I plug these five constants
> into to give me g(t), ie the lowpass filter output at a point t???

Hi,
I tried these coefficients in Matlab. They form a perfect Butterworth
lowpass of degree 2 and 3 dB cutoff at exactly 0.01 times the sampling
rate, when used as follows in a transfer function H(z):

a0 + a1*z^(-1) + a2*z^(-2)
H(z) = --------------------------
1 - b1*z^(-1) - b2*z^(-2)

This leads to the recursion formula:

g(t)= b1*g(t-1) + b2*g(t-2)
+ a0*f(t) + a1*f(t-1) + a2*f(t-2)

Try this recursion formual, it will work.

The formula presented by Tim Wescott is different (and it is definitely
wrong).

I agree with Bernhard Holzmeyer: usually, signs of b1 and b2 in H(z) are
positive, and the signs of these 2 coefficients in the recursion formula
are usually negative.

Siegbert
```
```Am 29 Mar 2004 19:30:00 -0800 schrieb Stewart:

> Hi All,
>
> I've been plowing through pages of theory in dsp books trying to work
> out how to create a butterworth filter.  They all seem to use
> different symbols and many different forms of the same equation.  Lots
> of theory but very little of practical examples.
>
> Basically I have a series of 1500 points taken at even intervals.
> Let's call it f(t).  Eg
>
> t    f(t)
> 1    3412
> 2    3432
> 3    3300
> 4    3433
> ......... etc
>
> Basically I'm looking to numerically produce g(t) which is f(t) with
> all frequencies higher than 0.01 (seconds/data point) filtered out.
>
> Keeping it simple - lets say I want a second order (ie 2 pole)
> Butterworth Filter this gives (from a table in a book)
>
> a0= 8.663387E-04
> a1= 1.732678E-03 b1= 1.919129E+00
> a2= 8.663387E-04 b2= -9.225943E-01
>
> Question is -  what is the formula that I plug these five constants
> into to give me g(t), ie the lowpass filter output at a point t???

Hi,
I tried these coefficients in Matlab. They form a perfect Butterworth
lowpass of degree 2 and 3 dB cutoff at exactly 0.01 times the sampling
rate, when used as follows in a transfer function H(z):

a0 + a1*z^(-1) + a2*z^(-2)
H(z) = --------------------------
1 - b1*z^(-1) - b2*z^(-2)

This leads to the recursion formula:

g(t)= b1*g(t-1) + b2*g(t-2)
+ a0*f(t) + a1*f(t-1) + a2*f(t-2)

Try this recursion with your data.

The formula presented by Tim Wescott is different (and it is definitely
wrong).

I agree with Bernhard Holzmayer: usually, signs of b1 and b2 in H(z) are
positive, and the signs of these 2 coefficients in the recursion formula
are usually negative.

Siegbert
```
```Tim Wescott wrote:
> Stewart wrote:
...
> > a0= 8.663387E-04
> > a1= 1.732678E-03 b1= 1.919129E+00
> > a2= 8.663387E-04 b2= -9.225943E-01
> >
> > Question is -  what is the formula that I plug these five constants
> > into to give me g(t), ie the lowpass filter output at a point t???
>
> First observation:  that doesn't look right.  The transfer function
> should be normalized (see OP) so that the highest-order term in the
> denominator is one:
>
>           b_2 z^2 + b_1 z
> H(z) = -------------------
>          z^2 + a_1 z + a_0

Signal processing guys prefer the technically realizable form:

H(z) = (a_0 + a_1 z^{-1} + a_2 z^{-2} )/(1 - (b_1 z^{-1} + b_2
z^{-2}))

and this looks perfectly allright to me. This system has a frequency
response of a bilinear transformed lowpass filter with cutoff
frequency Fc = 0.01 ( = fc / fs). It is realized by this recursion
equation:

g(t) = a_0 f(t) + a_1 f(t-T) + a_2 f(t-2T) + b_1 g(t-T) + b_2 g(t-2T)

where T is the sampling period.

> Second observation:  Why in the heck are you messing with IIR filters if
> you have a DSP?

There are many good reasons for IIR filters. If done right they
perform as expected (btw, the DSP MAC instruction is just as
applicable to IIR as to FIR filters).

Regards,
Andor
```
```Andor wrote:
> Tim Wescott wrote:
>
>>Stewart wrote:
>
> ...
>
>>>a0= 8.663387E-04
>>>a1= 1.732678E-03 b1= 1.919129E+00
>>>a2= 8.663387E-04 b2= -9.225943E-01
>>>
>>>Question is -  what is the formula that I plug these five constants
>>>into to give me g(t), ie the lowpass filter output at a point t???
>>
>>First observation:  that doesn't look right.  The transfer function
>>should be normalized (see OP) so that the highest-order term in the
>>denominator is one:
>>
>>          b_2 z^2 + b_1 z
>>H(z) = -------------------
>>         z^2 + a_1 z + a_0
>
>
> Signal processing guys prefer the technically realizable form:
>
> H(z) = (a_0 + a_1 z^{-1} + a_2 z^{-2} )/(1 - (b_1 z^{-1} + b_2
> z^{-2}))
>
> and this looks perfectly allright to me. This system has a frequency
> response of a bilinear transformed lowpass filter with cutoff
> frequency Fc = 0.01 ( = fc / fs). It is realized by this recursion
> equation:
>
> g(t) = a_0 f(t) + a_1 f(t-T) + a_2 f(t-2T) + b_1 g(t-T) + b_2 g(t-2T)
>
> where T is the sampling period.
>
>
>>Second observation:  Why in the heck are you messing with IIR filters if
>>you have a DSP?
>
>
> There are many good reasons for IIR filters. If done right they
> perform as expected (btw, the DSP MAC instruction is just as
> applicable to IIR as to FIR filters).
>
> Regards,
> Andor

Boy I'm getting a lot of grief for that comment -- which I suppose I
should, since I probably implement more IIR filters than FIR ones
myself.  And yes, the DSP MAC instruction is just as applicable, but if
you're only multiplying a few things together you spend as many clock
ticks in extra setup for the 1-cycle MAC is you save by having it happen
in one cycle.  If _all_ you're doing is one or two IIR filters you may
as well do it on a fast RISC processor as a DSP -- you only see value
from the one-cycle MAC if you can express a bunch of IIR filters as a
matrix multiply or if you're doing FIR filters.

--

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com
```