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???
Help with Butterworth Filter
Started by ●March 29, 2004
Reply by ●March 30, 20042004-03-30
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 news:3e8a4d01.0403291930.2b02e550@posting.google.com...> 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???
Reply by ●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
Reply by ●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 ●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 ●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 ●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 ●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 ●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_0Signal 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 ●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, > AndorBoy 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
