Reply by axlq July 16, 20042004-07-16
In article <3e8a4d01.0403291930.2b02e550@posting.google.com>,
Stewart <windsurfing_stew@yahoo.com.au> wrote:
>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???
I found the following article easy to understand, useful for writing code to test, and as a bonus, included a description of BOTH a 2-pole Butterworth and a 2-pole Critically-Damped filter, both adjusted to have a 3dB cutoff anywhere you specify, accounting for the sampling rate. It contains the formulas for the coefficients, not just tables of constants. http://people.umass.edu/exsci735/Robertson&Dowling.pdf Beware of combinations of cutoff and sampling rate that result in the adjusted angular cutoff frequency to be Wc >= tan(pi/2) which blows up (this is true for both kinds of filters described). It starts ringing at Wc=tan(pi/4), and I've found it best to keep Wc<tan(pi/8). -A
Reply by John Smith March 30, 20042004-03-30
"Butterworth" is used a lot in radio because of very low loss and nice phase
change characteristics, (real HW filter) however Stewart's problem, is
already digitized, so brick wall etc, something else easy to pick up should
work, just cutting out the highs.

"Tim Wescott" <tim@wescottnospamdesign.com> wrote in message
news:106k00ui3k5hj4a@corp.supernews.com...
> 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
Reply by Tim Wescott March 30, 20042004-03-30
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
Reply by Andor March 30, 20042004-03-30
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
Reply by Siegbert Steinlechner March 30, 20042004-03-30
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
Reply by Siegbert Steinlechner March 30, 20042004-03-30
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
Reply by Tim Wescott March 30, 20042004-03-30
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. > > > http://ccrma-www.stanford.edu/~jos/filters/BiQuad_Section.html has some more > 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: > http://www.dspguru.com/sw/lib/biquad.c > > >>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
Reply by Jon Harris March 30, 20042004-03-30
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.
http://ccrma-www.stanford.edu/~jos/filters/BiQuad_Section.html has some more 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: http://www.dspguru.com/sw/lib/biquad.c
> 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".
Reply by Bernhard Holzmayer March 30, 20042004-03-30
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





Reply by Tim Wescott March 30, 20042004-03-30
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