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.
Digital filter woes
Started by ●September 10, 2012
Reply by ●September 10, 20122012-09-10
On Mon, 10 Sep 2012 16:17:52 -0700 "garyr" <garyr@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.
Reply by ●September 10, 20122012-09-10
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 seco=nd> 100 samples and store it as one plot point.=20With only 200 time samples, you may still be within the startup transient o= f the filter. Instead of using this method, substitute z=3Dexp(jwT) into yo= ur transfer function. Then recognize that exp(jwT) =3D 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 filte= r 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.=20Converting from poles/zeroes to difference equations is almost trivially ea= sy. If you're doing that wrong it's probably because you're repeatedly maki= ng 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
Reply by ●September 10, 20122012-09-10
"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
Reply by ●September 11, 20122012-09-11
"Rob Gaddi" <rgaddi@technologyhighland.invalid> wrote in message news:20120910163140.6319a4e0@rg.highlandtechnology.com...> On Mon, 10 Sep 2012 16:17:52 -0700 > "garyr" <garyr@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.
Reply by ●September 11, 20122012-09-11
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 GaryI 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
Reply by ●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.comThe 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.
Reply by ●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 ●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.comThe 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 ●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.> sumx = x - self.x1*self.b1 + self.x2*self.b2 > sumy = self.y1*self.a1 - self.y2*self.a2have 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.006040jare 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.80134685i 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."
