Hello Thanks to website: http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html I have create a bandstop filter with this parameter: Parameter: ************************************************************************ filtertype = Butterworth passtype = Bandstop ripple = order = 3 samplerate = 500 corner1 = 48 corner2 = 52 ************************************************************************ This is C code: ************************************************************************ #define NZEROS 6 #define NPOLES 6 #define GAIN 1.051555795e+00 static float xv[NZEROS+1], yv[NPOLES+1]; static void filterloop() { for (;;) { xv[0] = xv[1]; xv[1] = xv[2]; xv[2] = xv[3]; xv[3] = xv[4]; xv[4] = xv[5]; xv[5] = xv[6]; xv[6] = next input value / GAIN; yv[0] = yv[1]; yv[1] = yv[2]; yv[2] = yv[3]; yv[3] = yv[4]; yv[4] = yv[5]; yv[5] = yv[6]; yv[6] = (xv[0] + xv[6]) - 4.8556354280 * (xv[1] + xv[5]) + 10.8590651360 * (xv[2] + xv[4]) - 13.9513547570 * xv[3] + ( -0.9043475314 * yv[0]) + ( 4.4647491458 * yv[1]) + (-10.1524477310 * yv[2]) + ( 13.2634555850 * yv[3]) + (-10.4984798610 * yv[4]) + ( 4.7742870210 * yv[5]); next output value = yv[6]; } } *************************************************************************** This is VB.NET code: ************************************************************************ '[x] Vettore_x(1) = Vettore_x(2) Vettore_x(2) = Vettore_x(3) Vettore_x(3) = Vettore_x(4) Vettore_x(4) = Vettore_x(5) Vettore_x(5) = Vettore_x(6) Vettore_x(6) = Vettore_x(7) Vettore_x(7) = Sample_Filtrato / Gain_3_Trappola_50hz '[y] Vettore_y(1) = Vettore_y(2) Vettore_y(2) = Vettore_y(3) Vettore_y(3) = Vettore_y(4) Vettore_y(4) = Vettore_y(5) Vettore_y(5) = Vettore_y(6) Vettore_y(6) = Vettore_y(7) Vettore_y(7) = Vettore_y(7) + (Vettore_x(1) + Vettore_x(7)) Vettore_y(7) = Vettore_y(7) + (-4.855635428 * (Vettore_x(2) + Vettore_x(6))) Vettore_y(7) = Vettore_y(7) + (10.85906513604 * (Vettore_x(3) + Vettore_x(5))) Vettore_y(7) = Vettore_y(7) + (-13.951354757 * Vettore_x(4)) Vettore_y(7) = Vettore_y(7) + (-0.9043475314 * Vettore_y(1)) Vettore_y(7) = Vettore_y(7) + (4.4647491458 * Vettore_y(2)) Vettore_y(7) = Vettore_y(7) + (-10.152447731 * Vettore_y(3)) Vettore_y(7) = Vettore_y(7) + (13.263455585 * Vettore_y(4)) Vettore_y(7) = Vettore_y(7) + (-10.498479861 * Vettore_y(5)) Vettore_y(7) = Vettore_y(7) + (4.774287021 * Vettore_y(6)) Sample_Filtrato = Vettore_y(7) ************************************************************************ It return me overflow on last row --> Sample_Filtrato = Vettore_y(7)[HERE] The range of samples come from 1 to 1024 Thanks
50Hz bandstop filter don't work ....
Started by ●May 21, 2012
Reply by ●May 21, 20122012-05-21
On 5/21/2012 6:57 AM, gpezzella wrote:> Hello > Thanks to website: > http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html > I have create a bandstop filter with this parameter: > > Parameter: > ************************************************************************ > filtertype = Butterworth > passtype = Bandstop > ripple = > order = 3 > samplerate = 500 > corner1 = 48 > corner2 = 52 > ************************************************************************ >The website shows the magnitude and phase response. So it appears to work. What "don't work"? Fred
Reply by ●May 21, 20122012-05-21
On Mon, 21 May 2012 08:57:30 -0500, gpezzella wrote:> Hello > Thanks to website: > http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html I have create > a bandstop filter with this parameter: > > Parameter: > ************************************************************************ > filtertype = Butterworth > passtype = Bandstop > ripple = > order = 3 > samplerate = 500 > corner1 = 48 > corner2 = 52 > ************************************************************************ > > > This is C code: > ************************************************************************ > #define NZEROS 6 > #define NPOLES 6 > #define GAIN 1.051555795e+00 > > static float xv[NZEROS+1], yv[NPOLES+1]; > > static void filterloop() > { for (;;) > { xv[0] = xv[1]; xv[1] = xv[2]; xv[2] = xv[3]; xv[3] = xv[4]; > xv[4] = > xv[5]; xv[5] = xv[6]; > xv[6] = next input value / GAIN; > yv[0] = yv[1]; yv[1] = yv[2]; yv[2] = yv[3]; yv[3] = yv[4]; > yv[4] = > yv[5]; yv[5] = yv[6]; > yv[6] = (xv[0] + xv[6]) - 4.8556354280 * (xv[1] + xv[5]) + > 10.8590651360 * (xv[2] + xv[4]) > - 13.9513547570 * xv[3] > + ( -0.9043475314 * yv[0]) + ( 4.4647491458 * > yv[1]) + (-10.1524477310 * yv[2]) + ( 13.2634555850 > * yv[3]) + (-10.4984798610 * yv[4]) + ( > 4.7742870210 * > yv[5]); > next output value = yv[6]; > } > } >***************************************************************************> > This is VB.NET code: > ************************************************************************ > '[x] > Vettore_x(1) = Vettore_x(2) > Vettore_x(2) = Vettore_x(3) > Vettore_x(3) = Vettore_x(4) > Vettore_x(4) = Vettore_x(5) > Vettore_x(5) = Vettore_x(6) > Vettore_x(6) = Vettore_x(7) > Vettore_x(7) = Sample_Filtrato / Gain_3_Trappola_50hz > > '[y] > Vettore_y(1) = Vettore_y(2) > Vettore_y(2) = Vettore_y(3) > Vettore_y(3) = Vettore_y(4) > Vettore_y(4) = Vettore_y(5) > Vettore_y(5) = Vettore_y(6) > Vettore_y(6) = Vettore_y(7) > > Vettore_y(7) = Vettore_y(7) + (Vettore_x(1) + Vettore_x(7)) > Vettore_y(7) = Vettore_y(7) + (-4.855635428 * (Vettore_x(2) + > Vettore_x(6))) > Vettore_y(7) = Vettore_y(7) + (10.85906513604 * (Vettore_x(3) + > Vettore_x(5))) > Vettore_y(7) = Vettore_y(7) + (-13.951354757 * Vettore_x(4)) > > Vettore_y(7) = Vettore_y(7) + (-0.9043475314 * Vettore_y(1)) > Vettore_y(7) = Vettore_y(7) + (4.4647491458 * Vettore_y(2)) > Vettore_y(7) = Vettore_y(7) + (-10.152447731 * Vettore_y(3)) > Vettore_y(7) = Vettore_y(7) + (13.263455585 * Vettore_y(4)) > Vettore_y(7) = Vettore_y(7) + (-10.498479861 * Vettore_y(5)) > Vettore_y(7) = Vettore_y(7) + (4.774287021 * Vettore_y(6)) > Sample_Filtrato = Vettore_y(7) > ************************************************************************ > > It return me overflow on last row --> Sample_Filtrato = > Vettore_y(7)[HERE] > > The range of samples come from 1 to 1024Without going any deeper into the code, you are implementing a 6-th order filter all in one whack. This is wrong. To implement an IIR filter you want to break it down into the smallest filters that you can, implement each one individually, then cascade them (i.e., the output of filter one feeds filter two, the output of filter two feeds filter three, and the output of filter three is your "final" output). In your case you have three complex pole pairs, so the smallest filter that you can decompose to while keeping all your numbers real is 2nd-order. There's a good chance that if you do this, your overflow problems will go away. Even if they don't, you still need to break the filter up. You're still left with a problem that you have poles at 1/10th the sampling frequency, and you're using single-precision floating point (with an effective precision of 25 bits in the mantissa). This means that you're losing over 7 bits of precision in the filter computation: (500 / 50)^2 = 0.01, 2^(-6) > 0.01 > 2^(-7), and your poles are resonant which exacerbates the problem. So if you're using 16-bit input data then single-precision floating point is barely sufficient, and if you're using more than 16-bit input data you're probably hosed. So: decompose that filter into 2nd-order sections, try it out, and let us know how it went!! -- 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 ●May 21, 20122012-05-21
On Mon, 21 May 2012 08:57:30 -0500, gpezzella wrote:> > It return me overflow on last row --> Sample_Filtrato = > Vettore_y(7)[HERE] > > The range of samples come from 1 to 1024Just out of curiosity -- from what language do the variable names come? Your English is impeccable, but "Sample_Filtrato" clearly isn't English. I'm thinking Italian ("Vettore" somehow says "Italian" to me), but I'm lame enough in the Romance languages that it could just as well be some Spanish language. And -- more pointless curiosity -- what does "Vettore" mean? ("Filtrato" I can guess at!). -- 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 ●May 21, 20122012-05-21
You might want to change your variables to double precision to see if the overflow goes away. If it does, follow the advice given above. If you're trying to remove 50 Hz AC interference, you might want to consider an LMS adaptive filter approach instead of the IIR.
Reply by ●May 21, 20122012-05-21
On Mon, 21 May 2012 08:57:30 -0500, gpezzella wrote:> Hello > Thanks to website: > http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html I have create > a bandstop filter with this parameter: > > Parameter: > ************************************************************************ > filtertype = Butterworth > passtype = Bandstop > ripple = > order = 3 > samplerate = 500 > corner1 = 48 > corner2 = 52 > ************************************************************************ > > > This is C code: > ************************************************************************ > #define NZEROS 6 > #define NPOLES 6 > #define GAIN 1.051555795e+00 > > static float xv[NZEROS+1], yv[NPOLES+1]; > > static void filterloop() > { for (;;) > { xv[0] = xv[1]; xv[1] = xv[2]; xv[2] = xv[3]; xv[3] = xv[4]; > xv[4] = > xv[5]; xv[5] = xv[6]; > xv[6] = next input value / GAIN; > yv[0] = yv[1]; yv[1] = yv[2]; yv[2] = yv[3]; yv[3] = yv[4]; > yv[4] = > yv[5]; yv[5] = yv[6]; > yv[6] = (xv[0] + xv[6]) - 4.8556354280 * (xv[1] + xv[5]) + > 10.8590651360 * (xv[2] + xv[4]) > - 13.9513547570 * xv[3] > + ( -0.9043475314 * yv[0]) + ( 4.4647491458 * > yv[1]) + (-10.1524477310 * yv[2]) + ( 13.2634555850 > * yv[3]) + (-10.4984798610 * yv[4]) + ( > 4.7742870210 * > yv[5]); > next output value = yv[6]; > } > } >***************************************************************************> > This is VB.NET code: > ************************************************************************ > '[x] > Vettore_x(1) = Vettore_x(2) > Vettore_x(2) = Vettore_x(3) > Vettore_x(3) = Vettore_x(4) > Vettore_x(4) = Vettore_x(5) > Vettore_x(5) = Vettore_x(6) > Vettore_x(6) = Vettore_x(7) > Vettore_x(7) = Sample_Filtrato / Gain_3_Trappola_50hz > > '[y] > Vettore_y(1) = Vettore_y(2) > Vettore_y(2) = Vettore_y(3) > Vettore_y(3) = Vettore_y(4) > Vettore_y(4) = Vettore_y(5) > Vettore_y(5) = Vettore_y(6) > Vettore_y(6) = Vettore_y(7) > > Vettore_y(7) = Vettore_y(7) + (Vettore_x(1) + Vettore_x(7)) > Vettore_y(7) = Vettore_y(7) + (-4.855635428 * (Vettore_x(2) + > Vettore_x(6))) > Vettore_y(7) = Vettore_y(7) + (10.85906513604 * (Vettore_x(3) + > Vettore_x(5))) > Vettore_y(7) = Vettore_y(7) + (-13.951354757 * Vettore_x(4)) > > Vettore_y(7) = Vettore_y(7) + (-0.9043475314 * Vettore_y(1)) > Vettore_y(7) = Vettore_y(7) + (4.4647491458 * Vettore_y(2)) > Vettore_y(7) = Vettore_y(7) + (-10.152447731 * Vettore_y(3)) > Vettore_y(7) = Vettore_y(7) + (13.263455585 * Vettore_y(4)) > Vettore_y(7) = Vettore_y(7) + (-10.498479861 * Vettore_y(5)) > Vettore_y(7) = Vettore_y(7) + (4.774287021 * Vettore_y(6)) > Sample_Filtrato = Vettore_y(7) > ************************************************************************ > > It return me overflow on last row --> Sample_Filtrato = > Vettore_y(7)[HERE] > > The range of samples come from 1 to 1024 > > ThanksAha -- you have made a mistake in translating your transfer function into a difference equation. As a consequence you've made a 7-th order unstable system instead of a 6th-order stable one. The really big hint here is that you have seven states where you should only have six -- I suspect that your mistake came in misunderstanding how to handle the leading z^6 term in the denominator polynomial: that coefficient = 1 is implied if you do the difference equation correctly, it DOES NOT have to be explicit. To give an example with a second-order system, if your transfer function is foo * z^2 + bar * z H(z) = ---------------------- z^2 + blah * z + tee then the difference equation is y(t+2) = foo * x(t+2) + bar * x(t+1) - blah * y(t+1) - tee * y(t) Shift this back in time by two, so you're getting today's outputs from today's values: y(t) = foo * x(t) + bar * x(t-1) - blah * y(t-1) - tee * y(t-2) So you have a _2nd_ order transfer function, which gives you a difference equation with _2_ state variables (y(t-1) and y(t-2)). Make sense? Go do some web searching until you find some code for a 2nd-order direct- form filter (either direct form 1 or 2 -- each has its advantages, either will work better than what you're doing). Decompose your transfer function _correctly_, implement the difference equations _correctly_, and I think you'll find joy. -- 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 ●May 21, 20122012-05-21
Tim Wescott wrote:> On Mon, 21 May 2012 08:57:30 -0500, gpezzella wrote: >> >> It return me overflow on last row --> Sample_Filtrato = >> Vettore_y(7)[HERE] >> >> The range of samples come from 1 to 1024 > > Just out of curiosity -- from what language do the variable names come? > Your English is impeccable, but "Sample_Filtrato" clearly isn't English. > I'm thinking Italian ("Vettore" somehow says "Italian" to me), but I'm > lame enough in the Romance languages that it could just as well be some > Spanish language. > > And -- more pointless curiosity -- what does "Vettore" mean? ("Filtrato" > I can guess at!). >I will guess "vector", Victor. Over. -- Les Cargill
Reply by ●May 21, 20122012-05-21
I don't see the transfer function error, but I do see: Vettore_y(7) = Vettore_y(7) + (Vettore_x(1) + Vettore_x(7)) should be: Vettore_y(7) = (Vettore_x(1) + Vettore_x(7)) unless Vettore_y(7) is set to zero by a mechanism I don't understand.
Reply by ●May 22, 20122012-05-22
On May 22, 1:57�am, "gpezzella" <gpezzella@n_o_s_p_a_m.yahoo.com> wrote:> Hello > Thanks to website:http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html > I have create a bandstop filter with this parameter: > > Parameter: > ************************************************************************ > filtertype � � � = � � � Butterworth > passtype � � � � = � � � Bandstop > ripple � = > order � �= � � � 3 > samplerate � � � = � � � 500 > corner1 �= � � � 48 > corner2 �= � � � 52 > ************************************************************************ > > This is C code: > ************************************************************************ > #define NZEROS 6 > #define NPOLES 6 > #define GAIN � 1.051555795e+00 > > static float xv[NZEROS+1], yv[NPOLES+1]; > > static void filterloop() > � { for (;;) > � � � { xv[0] = xv[1]; xv[1] = xv[2]; xv[2] = xv[3]; xv[3] = xv[4]; xv[4] = > xv[5]; xv[5] = xv[6]; > � � � � xv[6] = next input value / GAIN; > � � � � yv[0] = yv[1]; yv[1] = yv[2]; yv[2] = yv[3]; yv[3] = yv[4]; yv[4] = > yv[5]; yv[5] = yv[6]; > � � � � yv[6] = � (xv[0] + xv[6]) - � 4.8556354280 * (xv[1] + xv[5]) + > 10.8590651360 * (xv[2] + xv[4]) > � � � � � � � � � � �- �13.9513547570 * xv[3] > � � � � � � � � � � �+ ( -0.9043475314 * yv[0]) + ( �4.4647491458 * yv[1]) > � � � � � � � � � � �+ (-10.1524477310 * yv[2]) + ( 13.2634555850 * yv[3]) > � � � � � � � � � � �+ (-10.4984798610 * yv[4]) + ( �4.7742870210 * > yv[5]); > � � � � next output value = yv[6]; > � � � } > � } > *************************************************************************** > > This is VB.NET code: > ************************************************************************ > � � � � '[x] > � � � � Vettore_x(1) = Vettore_x(2) > � � � � Vettore_x(2) = Vettore_x(3) > � � � � Vettore_x(3) = Vettore_x(4) > � � � � Vettore_x(4) = Vettore_x(5) > � � � � Vettore_x(5) = Vettore_x(6) > � � � � Vettore_x(6) = Vettore_x(7) > � � � � Vettore_x(7) = Sample_Filtrato / Gain_3_Trappola_50hz > > � � � � '[y] > � � � � Vettore_y(1) = Vettore_y(2) > � � � � Vettore_y(2) = Vettore_y(3) > � � � � Vettore_y(3) = Vettore_y(4) > � � � � Vettore_y(4) = Vettore_y(5) > � � � � Vettore_y(5) = Vettore_y(6) > � � � � Vettore_y(6) = Vettore_y(7) > > � � � � Vettore_y(7) = Vettore_y(7) + (Vettore_x(1) + Vettore_x(7)) > � � � � Vettore_y(7) = Vettore_y(7) + (-4.855635428 * (Vettore_x(2) + > Vettore_x(6))) > � � � � Vettore_y(7) = Vettore_y(7) + (10.85906513604 * (Vettore_x(3) + > Vettore_x(5))) > � � � � Vettore_y(7) = Vettore_y(7) + (-13.951354757 * Vettore_x(4)) > > � � � � Vettore_y(7) = Vettore_y(7) + (-0.9043475314 * Vettore_y(1)) > � � � � Vettore_y(7) = Vettore_y(7) + (4.4647491458 * Vettore_y(2)) > � � � � Vettore_y(7) = Vettore_y(7) + (-10.152447731 * Vettore_y(3)) > � � � � Vettore_y(7) = Vettore_y(7) + (13.263455585 * Vettore_y(4)) > � � � � Vettore_y(7) = Vettore_y(7) + (-10.498479861 * Vettore_y(5)) > � � � � Vettore_y(7) = Vettore_y(7) + (4.774287021 * Vettore_y(6)) > � � � � Sample_Filtrato = Vettore_y(7) > ************************************************************************ > > It return me overflow on last row --> Sample_Filtrato = Vettore_y(7)[HERE] > > The range of samples come from 1 to 1024 > > Thanksxv[0] = xv[1]; xv[1] = xv[2]; xv[2] = xv[3]; xv[3] = xv[4]; xv[4] = xv[5]; xv[5] = xv[6] Should be the other way round ie xv[6]=xv[5] etc
Reply by ●May 24, 20122012-05-24
Hello dear friends, thanks for your help. I have resolved the problem... in a sense that the noise source not was my device but the oscilloscope himself (picoscope toy!!) Now with Professional Agilent MSO-X 2024A no noise is present at output of my high gain preamplifier. Just for curiosity, I have follow one of yours tricks putting two filter in cascade. The filter is a 3° LowPass Filter with Fcut= 200Hz Here the screenshot: https://www.dropbox.com/sh/5r28c1tbh61p0mg/luNus0Po-y P.s. Very good Tim, I live in Naples, Italy. Moreover: Vettore means Array Filtrato means Filtered Vettore_Filtrato = Filtered_Array
