I'm trying to design an LP IIR digital filter but not having much luck.
The bandwidth of my filter is too low. The filter specifications are: pass,
stop and sampling frequencies 4.8, 8.5 and 30 KHz respectively; maximum
pass band ripple 0.5 dB and minimum stop band attenuation 30 dB;
filter type Bessel. I obtained the poles and zeros with the Python Scipy
function iirdesign. There are seven poles and zeros, six of which form three
sets of conjugate pairs.

I  have written Python code to simulate the filter using floating-point
arithmetic. I implemented the filter as a cascade of three Direct Form I
second-order sections and one first-order section. I compute the frequency
response in the following manner: for each frequency I apply 100 samples of
a sinusoidal signal to the input and discard the outputs. I apply another
100 samples, save the output values and then compute the rms of this second
100 samples and store it as one plot point. This process indicates that the
bandwidth of my filter is about one-half of what it should be. I get similar
results if I apply 200 samples in each case. I also get similar results when
I plot the Fourier transform of the impulse response.

When I compute the filter response by evaluating the ratio of the product of
the zeros and the product of the poles the results are very close to the
specifications given above which indicates the pole/zero values obtained
from iirdesign() are correct, as expected. I believe I have computed the
coefficients of the difference equations correctly. I must have something
wrong but I can't find it. If I increase the values of start and stop
frequency I pass to iirdesign() I can generate coefficients that will cause
my filter to have an acceptable bandwidth but I didn't think that would be
necessary. Am I wrong in thinking that?

My apologies for the length of this post. I've tried to provide enough
detail to allow someone to make a reasonable guess as to where I might have
screwed up. Thanks in advance for your replies.

On Mon, 10 Sep 2012 16:17:52 -0700
"garyr" <g...@fidalgo.net> wrote:

> I'm trying to design an LP IIR digital filter but not having much luck.
> The bandwidth of my filter is too low. The filter specifications are: pass,
> stop and sampling frequencies 4.8, 8.5 and 30 KHz respectively; maximum
> pass band ripple 0.5 dB and minimum stop band attenuation 30 dB;
> filter type Bessel. I obtained the poles and zeros with the Python Scipy
> function iirdesign. There are seven poles and zeros, six of which form three
> sets of conjugate pairs.
> 

This is without actually simulating your design, but I think your
filter can't exist.  Bessel filters by definition have a very soft knee
in the frequency response.  It winds up working a bit like a
roundabout as you add poles; you always enter the same radius knee from
the same direction in the same place, but leave it later and later at a
steeper and steeper angle.

If you're really asking for a Bessel filter that's 30dB down under an
octave away from being less than 0.5dB down, just off the top of my head
that's impossible by about an order of magnitude.

For the kind of sharp knee you're talking about, you want a
Chebyschev.  If the linear phase of the Bessel matters, then you want
an FIR with a whole mess of taps.  If you specifically need the lack of
overshoot from the Bessel, it probably can't be achieved with a filter
with a knee that sharp.

-- 
Rob Gaddi, Highland Technology -- www.highlandtechnology.com
Email address domain is currently out of order.  See above to fix.

On Monday, September 10, 2012 6:21:24 PM UTC-5, garyr wrote:

> I compute the frequency
> response in the following manner: for each frequency I apply 100 samples of
> a sinusoidal signal to the input and discard the outputs. I apply another
> 100 samples, save the output values and then compute the rms of this second
> 100 samples and store it as one plot point. 

With only 200 time samples, you may still be within the startup transient of the filter. Instead
of using this method, substitute z=exp(jwT) into your transfer function. Then recognize that
exp(jwT) = cos(wT) + jsin(wT), compute the complex response at your frequencies of interest and
convert to magnitude/phase. That will give you the steady-state response of the filter instead
of something that is potentially corrupted by filter transients.

> When I compute the filter response by evaluating the ratio of the product of
> the zeros and the product of the poles the results are very close to the
> specifications given above which indicates the pole/zero values obtained
> from iirdesign() are correct, as expected. I believe I have computed the
> coefficients of the difference equations correctly. I must have something
> wrong but I can't find it. 

Converting from poles/zeroes to difference equations is almost trivially easy. If you're doing
that wrong it's probably because you're repeatedly making a simple mistake and not noticing it.
Take a break, find someone else to look over your work, and you'll probably find it.

Greg

"Greg Berchin" <g...@charter.net> wrote in message
news:c...@googlegroups.com...
On Monday, September 10, 2012 6:21:24 PM UTC-5, garyr wrote:

> I compute the frequency
> response in the following manner: for each frequency I apply 100 samples
> of
> a sinusoidal signal to the input and discard the outputs. I apply another
> 100 samples, save the output values and then compute the rms of this
> second
> 100 samples and store it as one plot point.

With only 200 time samples, you may still be within the startup transient of
the filter. Instead of using this method, substitute z=exp(jwT) into your
transfer function. Then recognize that exp(jwT) = cos(wT) + jsin(wT),
compute the complex response at your frequencies of interest and convert to
magnitude/phase. That will give you the steady-state response of the filter
instead of something that is potentially corrupted by filter transients.

> When I compute the filter response by evaluating the ratio of the product
> of
> the zeros and the product of the poles the results are very close to the
> specifications given above which indicates the pole/zero values obtained
> from iirdesign() are correct, as expected. I believe I have computed the
> coefficients of the difference equations correctly. I must have something
> wrong but I can't find it.

Converting from poles/zeroes to difference equations is almost trivially
easy. If you're doing that wrong it's probably because you're repeatedly
making a simple mistake and not noticing it. Take a break, find someone else
to look over your work, and you'll probably find it.

Greg
=============================================
The impulse response of the filter is quite small after 50 samples. Also,
the plot didn't change much in going from 100 to 200 samples.

I did substitute z = exp(jwT) in the ratio of the product of the poles &
zeros obtained from the iirdesign() function. The resulting plot looks
exactly as it should. The -3 dB and the -30 dB values are very close to the
specified values.

I've been over the conversion from poles/zero to the difference equations
many times; I agree is pretty straightforward.

As I mentioned, I can generate coefficients that result in my digital
version meeting my requirements so I have pushed on and generated a Python
version employing integer arithmetic; it behaves in the same manner. I've
also completed a Verilog version. I guess I'll find out when I get the
hardware working. This is a one-of project so if it doesn't work I've had
some fun fooling around with Gnuplot and MyHDL.

Thanks for your suggestions
Gary

"Rob Gaddi" <r...@technologyhighland.invalid> wrote in message 
news:2...@rg.highlandtechnology.com...
> On Mon, 10 Sep 2012 16:17:52 -0700
> "garyr" <g...@fidalgo.net> wrote:
>
>> I'm trying to design an LP IIR digital filter but not having much luck.
>> The bandwidth of my filter is too low. The filter specifications are: 
>> pass,
>> stop and sampling frequencies 4.8, 8.5 and 30 KHz respectively; maximum
>> pass band ripple 0.5 dB and minimum stop band attenuation 30 dB;
>> filter type Bessel. I obtained the poles and zeros with the Python Scipy
>> function iirdesign. There are seven poles and zeros, six of which form 
>> three
>> sets of conjugate pairs.
>>
>
> This is without actually simulating your design, but I think your
> filter can't exist.  Bessel filters by definition have a very soft knee
> in the frequency response.  It winds up working a bit like a
> roundabout as you add poles; you always enter the same radius knee from
> the same direction in the same place, but leave it later and later at a
> steeper and steeper angle.
>
> If you're really asking for a Bessel filter that's 30dB down under an
> octave away from being less than 0.5dB down, just off the top of my head
> that's impossible by about an order of magnitude.
>
> For the kind of sharp knee you're talking about, you want a
> Chebyschev.  If the linear phase of the Bessel matters, then you want
> an FIR with a whole mess of taps.  If you specifically need the lack of
> overshoot from the Bessel, it probably can't be achieved with a filter
> with a knee that sharp.
>
> -- 
> Rob Gaddi, Highland Technology -- www.highlandtechnology.com
> Email address domain is currently out of order.  See above to fix.

I added the filter type to my email at the last moment and got it wrong. I'm 
actually trying to design a Butterworth filter.

Thanks for your reply.

On Mon, 10 Sep 2012 19:38:18 -0700, garyr wrote:

> "Greg Berchin" <g...@charter.net> wrote in message
> news:c...@googlegroups.com...
> On Monday, September 10, 2012 6:21:24 PM UTC-5, garyr wrote:
> 
>> I compute the frequency response in the following manner: for each
>> frequency I apply 100 samples of a sinusoidal signal to the input and
>> discard the outputs. I apply another 100 samples, save the output
>> values and then compute the rms of this second 100 samples and store it
>> as one plot point.
> 
> With only 200 time samples, you may still be within the startup
> transient of the filter. Instead of using this method, substitute
> z=exp(jwT) into your transfer function. Then recognize that exp(jwT) > cos(wT) +
jsin(wT), compute the complex response at your frequencies of
> interest and convert to magnitude/phase. That will give you the
> steady-state response of the filter instead of something that is
> potentially corrupted by filter transients.
> 
>> When I compute the filter response by evaluating the ratio of the
>> product of the zeros and the product of the poles the results are very
>> close to the specifications given above which indicates the pole/zero
>> values obtained from iirdesign() are correct, as expected. I believe I
>> have computed the coefficients of the difference equations correctly. I
>> must have something wrong but I can't find it.
> 
> Converting from poles/zeroes to difference equations is almost trivially
> easy. If you're doing that wrong it's probably because you're repeatedly
> making a simple mistake and not noticing it. Take a break, find someone
> else to look over your work, and you'll probably find it.
> 
> Greg =============================================> 
> The impulse response of the filter is quite small after 50 samples.
> Also, the plot didn't change much in going from 100 to 200 samples.
> 
> I did substitute z = exp(jwT) in the ratio of the product of the poles &
> zeros obtained from the iirdesign() function. The resulting plot looks
> exactly as it should. The -3 dB and the -30 dB values are very close to
> the specified values.
> 
> I've been over the conversion from poles/zero to the difference
> equations many times; I agree is pretty straightforward.
> 
> As I mentioned, I can generate coefficients that result in my digital
> version meeting my requirements so I have pushed on and generated a
> Python version employing integer arithmetic; it behaves in the same
> manner. I've also completed a Verilog version. I guess I'll find out
> when I get the hardware working. This is a one-of project so if it
> doesn't work I've had some fun fooling around with Gnuplot and MyHDL.
> 
> Thanks for your suggestions Gary

I think you're not giving it enough samples for settling.  With a 
sampling/passband ratio of 4000:1, you should be letting it settle for 
about 4000 samples before you believe what's coming out.  One or two 
hundred just isn't going to be representative.

When you calculate your response with the z = exp(jwT) method, are you 
using the transfer functions from which you derived your filters, or are 
you using the actual filter difference equations?  If your error is in 
the difference equations or the software, and you're looking at the 
source transfer function, you're not going to catch it.

You may also have a problem with coefficient rounding, if you're doing 
this all in single-precision arithmetic.  It may be a good thing to try 
it with double-precision to see if that improves things.

-- 
Tim Wescott
Control system and signal processing consulting
www.wescottdesign.com

"Tim Wescott" <t...@seemywebsite.please> wrote in message 
news:m...@web-ster.com...
> On Mon, 10 Sep 2012 19:38:18 -0700, garyr wrote:
>
>> "Greg Berchin" <g...@charter.net> wrote in message
>> news:c...@googlegroups.com...
>> On Monday, September 10, 2012 6:21:24 PM UTC-5, garyr wrote:
>>
>>> I compute the frequency response in the following manner: for each
>>> frequency I apply 100 samples of a sinusoidal signal to the input and
>>> discard the outputs. I apply another 100 samples, save the output
>>> values and then compute the rms of this second 100 samples and store it
>>> as one plot point.
>>
>> With only 200 time samples, you may still be within the startup
>> transient of the filter. Instead of using this method, substitute
>> z=exp(jwT) into your transfer function. Then recognize that exp(jwT) >> cos(wT) +
jsin(wT), compute the complex response at your frequencies of
>> interest and convert to magnitude/phase. That will give you the
>> steady-state response of the filter instead of something that is
>> potentially corrupted by filter transients.
>>
>>> When I compute the filter response by evaluating the ratio of the
>>> product of the zeros and the product of the poles the results are very
>>> close to the specifications given above which indicates the pole/zero
>>> values obtained from iirdesign() are correct, as expected. I believe I
>>> have computed the coefficients of the difference equations correctly. I
>>> must have something wrong but I can't find it.
>>
>> Converting from poles/zeroes to difference equations is almost trivially
>> easy. If you're doing that wrong it's probably because you're repeatedly
>> making a simple mistake and not noticing it. Take a break, find someone
>> else to look over your work, and you'll probably find it.
>>
>> Greg =============================================>>
>> The impulse response of the filter is quite small after 50 samples.
>> Also, the plot didn't change much in going from 100 to 200 samples.
>>
>> I did substitute z = exp(jwT) in the ratio of the product of the poles &
>> zeros obtained from the iirdesign() function. The resulting plot looks
>> exactly as it should. The -3 dB and the -30 dB values are very close to
>> the specified values.
>>
>> I've been over the conversion from poles/zero to the difference
>> equations many times; I agree is pretty straightforward.
>>
>> As I mentioned, I can generate coefficients that result in my digital
>> version meeting my requirements so I have pushed on and generated a
>> Python version employing integer arithmetic; it behaves in the same
>> manner. I've also completed a Verilog version. I guess I'll find out
>> when I get the hardware working. This is a one-of project so if it
>> doesn't work I've had some fun fooling around with Gnuplot and MyHDL.
>>
>> Thanks for your suggestions Gary
>
> I think you're not giving it enough samples for settling.  With a
> sampling/passband ratio of 4000:1, you should be letting it settle for
> about 4000 samples before you believe what's coming out.  One or two
> hundred just isn't going to be representative.
>
> When you calculate your response with the z = exp(jwT) method, are you
> using the transfer functions from which you derived your filters, or are
> you using the actual filter difference equations?  If your error is in
> the difference equations or the software, and you're looking at the
> source transfer function, you're not going to catch it.
>
> You may also have a problem with coefficient rounding, if you're doing
> this all in single-precision arithmetic.  It may be a good thing to try
> it with double-precision to see if that improves things.
>
> -- 
> Tim Wescott
> Control system and signal processing consulting
> www.wescottdesign.com

The sampling/passband ratio is 30/4.8 or 6.25.

Using the z = exp(jwT) method I calculated the filter response using the 
pole/zero information returned from the Scipy function iirdesign. This 
information was also used to calculate the difference equation coefficients.

There is no double-precision arithmetic in Python. Basic floating point 
arithmetic is done with 64-bit values, I believe. There is a Decimal class 
that allows arbitrary precision but I did not find that necessary. In my 
fixed-point version of the filter I used 16-bit coefficients and 12-bit 
data. The resulting frequency response plot was very similar to that 
obtained with the floating point version.

Thanks for your reply. 

On Tue, 11 Sep 2012 07:50:48 -0700, garyr wrote:

> "Tim Wescott" <t...@seemywebsite.please> wrote in message
> news:m...@web-ster.com...
>> On Mon, 10 Sep 2012 19:38:18 -0700, garyr wrote:
>>
>>> "Greg Berchin" <g...@charter.net> wrote in message
>>> news:c...@googlegroups.com... On
>>> Monday, September 10, 2012 6:21:24 PM UTC-5, garyr wrote:
>>>
>>>> I compute the frequency response in the following manner: for each
>>>> frequency I apply 100 samples of a sinusoidal signal to the input and
>>>> discard the outputs. I apply another 100 samples, save the output
>>>> values and then compute the rms of this second 100 samples and store
>>>> it as one plot point.
>>>
>>> With only 200 time samples, you may still be within the startup
>>> transient of the filter. Instead of using this method, substitute
>>> z=exp(jwT) into your transfer function. Then recognize that exp(jwT) >>>
cos(wT) + jsin(wT), compute the complex response at your frequencies
>>> of interest and convert to magnitude/phase. That will give you the
>>> steady-state response of the filter instead of something that is
>>> potentially corrupted by filter transients.
>>>
>>>> When I compute the filter response by evaluating the ratio of the
>>>> product of the zeros and the product of the poles the results are
>>>> very close to the specifications given above which indicates the
>>>> pole/zero values obtained from iirdesign() are correct, as expected.
>>>> I believe I have computed the coefficients of the difference
>>>> equations correctly. I must have something wrong but I can't find it.
>>>
>>> Converting from poles/zeroes to difference equations is almost
>>> trivially easy. If you're doing that wrong it's probably because
>>> you're repeatedly making a simple mistake and not noticing it. Take a
>>> break, find someone else to look over your work, and you'll probably
>>> find it.
>>>
>>> Greg =============================================>>>
>>> The impulse response of the filter is quite small after 50 samples.
>>> Also, the plot didn't change much in going from 100 to 200 samples.
>>>
>>> I did substitute z = exp(jwT) in the ratio of the product of the poles
>>> & zeros obtained from the iirdesign() function. The resulting plot
>>> looks exactly as it should. The -3 dB and the -30 dB values are very
>>> close to the specified values.
>>>
>>> I've been over the conversion from poles/zero to the difference
>>> equations many times; I agree is pretty straightforward.
>>>
>>> As I mentioned, I can generate coefficients that result in my digital
>>> version meeting my requirements so I have pushed on and generated a
>>> Python version employing integer arithmetic; it behaves in the same
>>> manner. I've also completed a Verilog version. I guess I'll find out
>>> when I get the hardware working. This is a one-of project so if it
>>> doesn't work I've had some fun fooling around with Gnuplot and MyHDL.
>>>
>>> Thanks for your suggestions Gary
>>
>> I think you're not giving it enough samples for settling.  With a
>> sampling/passband ratio of 4000:1, you should be letting it settle for
>> about 4000 samples before you believe what's coming out.  One or two
>> hundred just isn't going to be representative.
>>
>> When you calculate your response with the z = exp(jwT) method, are you
>> using the transfer functions from which you derived your filters, or
>> are you using the actual filter difference equations?  If your error is
>> in the difference equations or the software, and you're looking at the
>> source transfer function, you're not going to catch it.
>>
>> You may also have a problem with coefficient rounding, if you're doing
>> this all in single-precision arithmetic.  It may be a good thing to try
>> it with double-precision to see if that improves things.
>>
>> --
>> Tim Wescott
>> Control system and signal processing consulting www.wescottdesign.com
> 
> The sampling/passband ratio is 30/4.8 or 6.25.

Ah.  You meant 4.8 _kilo_Hz, not 4.8Hz.  That does make a difference.

> Using the z = exp(jwT) method I calculated the filter response using the
> pole/zero information returned from the Scipy function iirdesign. This
> information was also used to calculate the difference equation
> coefficients.

Did you go back and calculate the transfer function _by hand_ from the 
difference equation?  If you haven't, you haven't checked that step.

http://www.wescottdesign.com/articles/zTransform/z-transforms.html

Your two most likely points of error are going from the pole/zero values 
to the difference equation, and going from the difference equation to 
code -- and these are the steps you haven't checked in any way but 
running samples through.  

Normally I discourage posting code ('cause I don't like plowing through 
it), but if you post a set of pole-zero values for one stage of the 
filter, and your difference equation, and the code (which will hopefully 
be three or four lines, and certainly should be less than ten), then 
someone will take a look.  Maybe even me.

> There is no double-precision arithmetic in Python. Basic floating point
> arithmetic is done with 64-bit values, I believe.

IEEE double-precision floating point _is_ 64 bits.

> There is a Decimal
> class that allows arbitrary precision but I did not find that necessary.
> In my fixed-point version of the filter I used 16-bit coefficients and
> 12-bit data. The resulting frequency response plot was very similar to
> that obtained with the floating point version.
> 
> Thanks for your reply.

-- 
My liberal friends think I'm a conservative kook.
My conservative friends think I'm a liberal kook.
Why am I not happy that they have found common ground?

Tim Wescott, Communications, Control, Circuits & Software
http://www.wescottdesign.com

"Tim Wescott" <t...@seemywebsite.com> wrote in message
news:E...@web-ster.com...
> On Tue, 11 Sep 2012 07:50:48 -0700, garyr wrote:
>
>> "Tim Wescott" <t...@seemywebsite.please> wrote in message
>> news:m...@web-ster.com...
>>> On Mon, 10 Sep 2012 19:38:18 -0700, garyr wrote:
>>>
>>>> "Greg Berchin" <g...@charter.net> wrote in message
>>>> news:c...@googlegroups.com... On
>>>> Monday, September 10, 2012 6:21:24 PM UTC-5, garyr wrote:
>>>>
>>>>> I compute the frequency response in the following manner: for each
>>>>> frequency I apply 100 samples of a sinusoidal signal to the input and
>>>>> discard the outputs. I apply another 100 samples, save the output
>>>>> values and then compute the rms of this second 100 samples and store
>>>>> it as one plot point.
>>>>
>>>> With only 200 time samples, you may still be within the startup
>>>> transient of the filter. Instead of using this method, substitute
>>>> z=exp(jwT) into your transfer function. Then recognize that exp(jwT)
>>>> cos(wT) + jsin(wT), compute the complex response at your frequencies
>>>> of interest and convert to magnitude/phase. That will give you the
>>>> steady-state response of the filter instead of something that is
>>>> potentially corrupted by filter transients.
>>>>
>>>>> When I compute the filter response by evaluating the ratio of the
>>>>> product of the zeros and the product of the poles the results are
>>>>> very close to the specifications given above which indicates the
>>>>> pole/zero values obtained from iirdesign() are correct, as expected.
>>>>> I believe I have computed the coefficients of the difference
>>>>> equations correctly. I must have something wrong but I can't find it.
>>>>
>>>> Converting from poles/zeroes to difference equations is almost
>>>> trivially easy. If you're doing that wrong it's probably because
>>>> you're repeatedly making a simple mistake and not noticing it. Take a
>>>> break, find someone else to look over your work, and you'll probably
>>>> find it.
>>>>
>>>> Greg =============================================>>>>
>>>> The impulse response of the filter is quite small after 50 samples.
>>>> Also, the plot didn't change much in going from 100 to 200 samples.
>>>>
>>>> I did substitute z = exp(jwT) in the ratio of the product of the poles
>>>> & zeros obtained from the iirdesign() function. The resulting plot
>>>> looks exactly as it should. The -3 dB and the -30 dB values are very
>>>> close to the specified values.
>>>>
>>>> I've been over the conversion from poles/zero to the difference
>>>> equations many times; I agree is pretty straightforward.
>>>>
>>>> As I mentioned, I can generate coefficients that result in my digital
>>>> version meeting my requirements so I have pushed on and generated a
>>>> Python version employing integer arithmetic; it behaves in the same
>>>> manner. I've also completed a Verilog version. I guess I'll find out
>>>> when I get the hardware working. This is a one-of project so if it
>>>> doesn't work I've had some fun fooling around with Gnuplot and MyHDL.
>>>>
>>>> Thanks for your suggestions Gary
>>>
>>> I think you're not giving it enough samples for settling.  With a
>>> sampling/passband ratio of 4000:1, you should be letting it settle for
>>> about 4000 samples before you believe what's coming out.  One or two
>>> hundred just isn't going to be representative.
>>>
>>> When you calculate your response with the z = exp(jwT) method, are you
>>> using the transfer functions from which you derived your filters, or
>>> are you using the actual filter difference equations?  If your error is
>>> in the difference equations or the software, and you're looking at the
>>> source transfer function, you're not going to catch it.
>>>
>>> You may also have a problem with coefficient rounding, if you're doing
>>> this all in single-precision arithmetic.  It may be a good thing to try
>>> it with double-precision to see if that improves things.
>>>
>>> --
>>> Tim Wescott
>>> Control system and signal processing consulting www.wescottdesign.com
>>
>> The sampling/passband ratio is 30/4.8 or 6.25.
>
> Ah.  You meant 4.8 _kilo_Hz, not 4.8Hz.  That does make a difference.
>
>> Using the z = exp(jwT) method I calculated the filter response using the
>> pole/zero information returned from the Scipy function iirdesign. This
>> information was also used to calculate the difference equation
>> coefficients.
>
> Did you go back and calculate the transfer function _by hand_ from the
> difference equation?  If you haven't, you haven't checked that step.
>
> http://www.wescottdesign.com/articles/zTransform/z-transforms.html
>
> Your two most likely points of error are going from the pole/zero values
> to the difference equation, and going from the difference equation to
> code -- and these are the steps you haven't checked in any way but
> running samples through.
>
> Normally I discourage posting code ('cause I don't like plowing through
> it), but if you post a set of pole-zero values for one stage of the
> filter, and your difference equation, and the code (which will hopefully
> be three or four lines, and certainly should be less than ten), then
> someone will take a look.  Maybe even me.
>
>> There is no double-precision arithmetic in Python. Basic floating point
>> arithmetic is done with 64-bit values, I believe.
>
> IEEE double-precision floating point _is_ 64 bits.
>
>> There is a Decimal
>> class that allows arbitrary precision but I did not find that necessary.
>> In my fixed-point version of the filter I used 16-bit coefficients and
>> 12-bit data. The resulting frequency response plot was very similar to
>> that obtained with the floating point version.
>>
>> Thanks for your reply.
>
> -- 
> My liberal friends think I'm a conservative kook.
> My conservative friends think I'm a liberal kook.
> Why am I not happy that they have found common ground?
>
> Tim Wescott, Communications, Control, Circuits & Software
> http://www.wescottdesign.com

The following is the code for a second-order section. x is the current
input, x1 & x2 are the two previous values. y is the output, y1 & y2 are the
two previous output values. g2 is the reciprocal of the filter gain.

   def next(self, x):
       # y(n) = b0*x(n) - b1*x(n-1) + b2*x(n-2) + a1*y(n-1) - a2*y(n-2)
        sumx = x - self.x1*self.b1 + self.x2*self.b2
        sumy = self.y1*self.a1 - self.y2*self.a2
        y = (sumx*self.g2) + sumy
        self.x2 = self.x1
        self.x1 = x
        self.y2 = self.y1
        self.y1 = y
        return y

Given the conjugate pair of zeros (z - z0), (z - z0*) and poles (z - p0),
(z - p0*) the coefficients are given by: b0 = 1.0; b1 = z0 + z0*; b2 = the
product of z0 and z0*. Similarly, a1 = p0 + p0*, a2 = the product of p0 and
p0*. The gain coefficient is g2 = (1.0 - a1 + a2)/(1.0 - b1 + b2)

Here is a set of typical pole/zero values for a second-order section:

 p:  0.742966    0.499348j           z:  -1.013548    0.006040j
 p:  0.742966   -0.499348j           z:  -1.013548   -0.006040j

The corresponding coefficients are:
b1: -2.02709587   b2: 1.02731590
a1:  1.48593139   a2:  0.80134685

i was about to suggest that gary post his butterworth filter specs and 
the resulting coefficients from cooking those specs.



On 9/11/12 10:36 PM, garyr wrote:
...
>
> The following is the code for a second-order section. x is the current
> input, x1&  x2 are the two previous values. y is the output, y1&  y2 are the
> two previous output values. g2 is the reciprocal of the filter gain.
>
>     def next(self, x):
>         # y(n) = b0*x(n) - b1*x(n-1) + b2*x(n-2) + a1*y(n-1) - a2*y(n-2)

why alternate the signs with b1 and with a1?  the bees should have 
positive sign and the ayes should have negative.  i s'pose it's just 
another convention, but i don't see what motivates it.

>          sumx = x - self.x1*self.b1 + self.x2*self.b2
>          sumy = self.y1*self.a1 - self.y2*self.a2

have you guaranteed that sumx and sumy are initialized to zero?

and it appears that

    self.g2 = b0
    self.b1 = b1/b0     as you have defined b1.
    self.b2 = b2/b0


>          y = (sumx*self.g2) + sumy
>          self.x2 = self.x1
>          self.x1 = x
>          self.y2 = self.y1
>          self.y1 = y
>          return y
>
> Given the conjugate pair of zeros (z - z0), (z - z0*) and poles (z - p0),
> (z - p0*) the coefficients are given by: b0 = 1.0; b1 = z0 + z0*; b2 = the
> product of z0 and z0*. Similarly, a1 = p0 + p0*, a2 = the product of p0 and
> p0*. The gain coefficient is g2 = (1.0 - a1 + a2)/(1.0 - b1 + b2)
>
> Here is a set of typical pole/zero values for a second-order section:
>
>   p:  0.742966    0.499348j           z:  -1.013548    0.006040j
>   p:  0.742966   -0.499348j           z:  -1.013548   -0.006040j

are these z-plane poles and zeros or are they s-plane poles and zeros?

so it's two identical 2nd-order sections?


> The corresponding coefficients are:
> b1: -2.02709587   b2: 1.02731590
> a1:  1.48593139   a2:  0.80134685

i believe that a2 = |p|^2 .

i don't think that 0.80134685 = 0.742966^2 + 0.499348^2

or maybe it does.

now does self.a1 = +1.48593139  or  -1.48593139 ?

"inquiring minds want to know."

-- 

r b-j                  r...@audioimagination.com

"Imagination is more important than knowledge."