"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
Reply by Tim Wescott●September 13, 20122012-09-13
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
Reply by garyr●September 13, 20122012-09-13
>
> 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.
Reply by Tim Wescott●September 13, 20122012-09-13
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
>>>>> news:cb3b8189-0a4a-4234-aec9-0859fc4b6019@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
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
Reply by robert bristow-johnson●September 12, 20122012-09-12
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."
Reply by garyr●September 12, 20122012-09-12
"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)
Reply by robert bristow-johnson●September 12, 20122012-09-12
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.
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?
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."
Reply by garyr●September 11, 20122012-09-11
"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
>>>> news:cb3b8189-0a4a-4234-aec9-0859fc4b6019@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
Reply by Tim Wescott●September 11, 20122012-09-11
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
>>> news:cb3b8189-0a4a-4234-aec9-0859fc4b6019@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
Reply by garyr●September 11, 20122012-09-11
"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
>> news:cb3b8189-0a4a-4234-aec9-0859fc4b6019@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.