"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)
Digital filter woes
Started by ●September 10, 2012
Reply by ●September 12, 20122012-09-12
Reply by ●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 ●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.80134685Your 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 ●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.comIf you will post the s-plane poles I'll compute the z-plane coefficients using the bilinear transform.
Reply by ●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 ●September 13, 20122012-09-13
"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.comThe 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
