Adaptive Filter | ADPCM | ADSP | ADSP-2181 | Aliasing | AMR | Anti-Aliasing | ARMA | Autocorrelation | AutoCovariance | Beamforming | Bessel | Blackfin | Butterworth | C6713 | CCS | Chebyshev | CIC Filter | Circular Convolution | Code Composer Studio | Comb Filter | Compression | Convolution | Cross Correlation | DCT | Decimation | Deconvolution | Demodulation | DM642 | DSP Boards | DSP/BIOS | DTMF | Echo Cancellation | Equalization | Equalizer | ETSI | EZLITE (Ez-kit Lite) | FFT | FFTW | FIR Filter | Fixed Point | FSK | G.711 | G.723 | G.729 | Gaussian Noise | Goertzel | GPIO | Hilbert Transform | IFFT | IIR Filter | Interpolation | Invariance | JTAG | Kalman | Laplace Transform | Levinson | LPC | McBSP | MIPS | Modulation | MPEG | Multirate | Notch Filter | Nyquist | OFDM | Oversampling | Pink Noise | Pitch | PLL | Polyphase | QAM | QDMA | Quantization | Quantizer | Radar | Random Noise | Reed Solomon | Remez | Resampling | RTDX | Sampling | Sharc | TI C6711 | Undersampling | Viterbi | Wavelets | White Noise | Wiener Filter | Windowing | XDS510PP | Z Transform

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

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

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" <w...@yahoo.com.au> wrote in message
news:3...@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???

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

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

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

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

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Tim Wescott <t...@wescottnospamdesign.com> wrote in message
news:1...@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.

www.stanford.edu/~jos/filters/BiQuad_Section.html" target=_blank rel="nofollow">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".

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

Jon Harris wrote:

> Tim Wescott <t...@wescottnospamdesign.com> wrote in message
> news:1...@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.
> 
> 
> www.stanford.edu/~jos/filters/BiQuad_Section.html" target=_blank rel="nofollow">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

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

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

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

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

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

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

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

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

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.