# Digital filter woes

Started by September 10, 2012
```"robert bristow-johnson" <rbj@audioimagination.com> wrote in message
news:k2ov9g\$l50\$1@dont-email.me...
>
> 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?

What possible difference could that make? sumx & sumy are just temporary
variables.

>
> 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                  rbj@audioimagination.com
>
> "Imagination is more important than knowledge."
>
>
I described the filter specs in my first post but had the type as Besses
which was wrong. It is Butterworth.
The specs are: pass, stop & sampling frequency are 4.8 KHz, 8.5 KHz and 30
KHz; maximum passband ripple 0.5 dB, minimum stopband attenuation 30 dB.
With those parameters the Scipy function iirdesign returns 7 Z-plane poles
and zeros, six of which form three sets of conjugate pole/zero pairs like
the above.

In the above:
a1 =  (0.742966,  0.499348j) + (0.742966,  -0.499348j)
a2 =  (0.742966,  0.499348j) * (0.742966,  -0.499348j)

```
```On 9/12/12 10:34 AM, garyr wrote:
> "robert bristow-johnson"<rbj@audioimagination.com>  wrote in message
..
>>>           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?
>
> What possible difference could that make? sumx&  sumy are just temporary
> variables.
>

i misread it.  i was thinking "sumx += ..."  in the C language.

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

now this question needs answering.  where did you get these poles and
zeros from?  did some design program give you these?  i happen to know
that Butterworth is minimum phase, and it appears that the zeros lie
outside the unit circle.  i wonder if they were meant to be a -1 + 0j
(that would be Butterworth LPF mapped to z-plane with the bilinear
transform).

>>
>>> 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."
>>

this sign issue is important.  what is "self.a1"?

> I described the filter specs in my first post but had the type as Bessel
> which was wrong. It is Butterworth.
> The specs are: pass, stop&  sampling frequency are 4.8 KHz, 8.5 KHz and 30
> KHz; maximum passband ripple 0.5 dB, minimum stopband attenuation 30 dB.
> With those parameters the Scipy function iirdesign returns 7 Z-plane poles
> and zeros, six of which form three sets of conjugate pole/zero pairs like
> the above.
>
> In the above:
> a1 =  (0.742966,  0.499348j) + (0.742966,  -0.499348j)
> a2 =  (0.742966,  0.499348j) * (0.742966,  -0.499348j)
>

how does scipy map the Butterworth s-plane poles (there are no zeros, or
the zeros live out at infinity) to z-plane.  if it's bilinear transform
(which, i think, is the most logical because the frequency warping can
be easily accounted for), then i can tell there is something wrong with
the zeros.

--

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."

```
```On Tue, 11 Sep 2012 19:36:44 -0700, garyr wrote:

> "Tim Wescott" <tim@seemywebsite.com> wrote in message
> news:EoidndGxxu0YE9LNnZ2dnUVZ_vWdnZ2d@web-ster.com...
>> On Tue, 11 Sep 2012 07:50:48 -0700, garyr wrote:
>>
>>> "Tim Wescott" <tim@seemywebsite.please> wrote in message
>>> news:mIWdnZU-B7K7RdPNnZ2dnUVZ_sSdnZ2d@web-ster.com...
>>>> On Mon, 10 Sep 2012 19:38:18 -0700, garyr wrote:
>>>>
>>>>> "Greg Berchin" <gjberchin@charter.net> wrote in message
>>>>> 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
>>>>> 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.
>>>
>>
>> --
>> 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

Your given pole values don't look right for a filter with the given
passband -- are you sure that's not where your factor of two problem is?

A really hack job of this in Scilab* gives me pole values of

0.3824408 + 0.6852941i
0.3341172 + 0.3814747i
0.3337718 + 0.1706111i
0.3365241 + 4.488D-17i
0.3337718 - 0.1706111i
0.3341172 - 0.3814747i
0.3824408 - 0.6852941i

I wouldn't expect them to be exactly right -- but the order of magnitude
should be so, and doing a Bode plot of H = 1 / prod(z - pols) shows the
right frequency range.

* I just found the poles of an analog Butterworth filter, and found their
equivalents in the z plane with exp(pole / 30000).

--
Tim Wescott
Control system and signal processing consulting
www.wescottdesign.com
```
```>
> Your given pole values don't look right for a filter with the given
> passband -- are you sure that's not where your factor of two problem is?
>
> A really hack job of this in Scilab* gives me pole values of
>
>    0.3824408 + 0.6852941i
>    0.3341172 + 0.3814747i
>    0.3337718 + 0.1706111i
>    0.3365241 + 4.488D-17i
>    0.3337718 - 0.1706111i
>    0.3341172 - 0.3814747i
>    0.3824408 - 0.6852941i
>
> I wouldn't expect them to be exactly right -- but the order of magnitude
> should be so, and doing a Bode plot of H = 1 / prod(z - pols) shows the
> right frequency range.
>
> * I just found the poles of an analog Butterworth filter, and found their
> equivalents in the z plane with exp(pole / 30000).
>
> --
> Tim Wescott
> Control system and signal processing consulting
> www.wescottdesign.com

If you will post the s-plane poles I'll compute the z-plane coefficients
using the bilinear transform.

```
```On Thu, 13 Sep 2012 09:14:05 -0700, garyr wrote:

>> Your given pole values don't look right for a filter with the given
>> passband -- are you sure that's not where your factor of two problem
>> is?
>>
>> A really hack job of this in Scilab* gives me pole values of
>>
>>    0.3824408 + 0.6852941i
>>    0.3341172 + 0.3814747i
>>    0.3337718 + 0.1706111i
>>    0.3365241 + 4.488D-17i
>>    0.3337718 - 0.1706111i
>>    0.3341172 - 0.3814747i
>>    0.3824408 - 0.6852941i
>>
>> I wouldn't expect them to be exactly right -- but the order of
>> magnitude should be so, and doing a Bode plot of H = 1 / prod(z - pols)
>> shows the right frequency range.
>>
>> * I just found the poles of an analog Butterworth filter, and found
>> their equivalents in the z plane with exp(pole / 30000).
>>
>> --
>> Tim Wescott
>> Control system and signal processing consulting www.wescottdesign.com
>
> If you will post the s-plane poles I'll compute the z-plane coefficients
> using the bilinear transform.

Be sure to pre-warp.  These aren't exactly the ones I gave before -- it's
not paying work, so I didn't take notes.  In Scilab I used

[pols, gain] = zpbutt(7, 2 * %pi * 0.6 * 8500)

to get these:

- 7130.5153 + 31240.829i
- 19979.26 + 25053.2i
- 28870.867 + 13903.477i
- 32044.245 + 3.924D-12i
- 28870.867 - 13903.477i
- 19979.26 - 25053.2i
- 7130.5153 - 31240.829i

I did double-check the response -- it's pretty good.

Here's the z-plane poles I get by direct frequency translation this time,
using

exp(pols / 30000) .'

to get:

0.3982046 + 0.6805068i
0.3447918 + 0.3808943i
0.3416957 + 0.1707631i
0.3436466 + 4.495D-17i
0.3416957 - 0.1707631i
0.3447918 - 0.3808943i
0.3982046 - 0.6805068i

--
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" <tim@seemywebsite.com> wrote in message
news:Cr2dnaF9j9Jsk8_NnZ2dnUVZ_gOdnZ2d@web-ster.com...
> On Thu, 13 Sep 2012 09:14:05 -0700, garyr wrote:
>
>
>>> Your given pole values don't look right for a filter with the given
>>> passband -- are you sure that's not where your factor of two problem
>>> is?
>>>
>>> A really hack job of this in Scilab* gives me pole values of
>>>
>>>    0.3824408 + 0.6852941i
>>>    0.3341172 + 0.3814747i
>>>    0.3337718 + 0.1706111i
>>>    0.3365241 + 4.488D-17i
>>>    0.3337718 - 0.1706111i
>>>    0.3341172 - 0.3814747i
>>>    0.3824408 - 0.6852941i
>>>
>>> I wouldn't expect them to be exactly right -- but the order of
>>> magnitude should be so, and doing a Bode plot of H = 1 / prod(z - pols)
>>> shows the right frequency range.
>>>
>>> * I just found the poles of an analog Butterworth filter, and found
>>> their equivalents in the z plane with exp(pole / 30000).
>>>
>>> --
>>> Tim Wescott
>>> Control system and signal processing consulting www.wescottdesign.com
>>
>> If you will post the s-plane poles I'll compute the z-plane coefficients
>> using the bilinear transform.
>
> Be sure to pre-warp.  These aren't exactly the ones I gave before -- it's
> not paying work, so I didn't take notes.  In Scilab I used
>
> [pols, gain] = zpbutt(7, 2 * %pi * 0.6 * 8500)
>
> to get these:
>
>  - 7130.5153 + 31240.829i
>  - 19979.26 + 25053.2i
>  - 28870.867 + 13903.477i
>  - 32044.245 + 3.924D-12i
>  - 28870.867 - 13903.477i
>  - 19979.26 - 25053.2i
>  - 7130.5153 - 31240.829i
>
> I did double-check the response -- it's pretty good.
>
> Here's the z-plane poles I get by direct frequency translation this time,
> using
>
> exp(pols / 30000) .'
>
> to get:
>
>    0.3982046 + 0.6805068i
>    0.3447918 + 0.3808943i
>    0.3416957 + 0.1707631i
>    0.3436466 + 4.495D-17i
>    0.3416957 - 0.1707631i
>    0.3447918 - 0.3808943i
>    0.3982046 - 0.6805068i
>
> --
> 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 pre-warp didn't cause much change: 2*30*tan(0.5*4.8/30) = 4.81
Applying the bilinear transform to your s-plane values I get:
-7130.5153   31240.8290j :      0.469342,    0.683794
-19979.2600   25053.2000j :      0.366321,    0.427995
-28870.8670   13903.4770j :      0.318015,    0.206198
-32044.2450       0.0000j :      0.303721,    0.000000
-28870.8670  -13903.4770j :      0.318015,   -0.206198
-19979.2600  -25053.2000j :      0.366321,   -0.427995
-7130.5153  -31240.8290j :      0.469342,   -0.683794
Plus 7 zeros at z = -1.
I computed the pole positions as (A + s)/(A - s) with A = (2/ts) = 2*30e3.
The plot had the -3 dB point at 4.7 KHz and -30 at 6.9 KHz
The difference equation coefficients are:
a1              a2              b1                 b2
0     0.938684     0.687857   -2.000000     1.000000
1     0.732642     0.317371   -2.000000     1.000000
2     0.636030     0.143651   -2.000000     1.000000
And for the single pole:   b1=  0.303721   b2= -1.000000
The digital filter response was: -3 dB at 4.6 KHz and -30 at 6.8 KHz. An
excellent filter for my application. I still don't know exactly what the
problem is but at least I do know somewhat more now than I did.

Thanks very much for your interest and help with my problem.
Gary Richardson

```