## filter design - code problem

If you require help or assistance with anything then please post here

Moderators: electrogear, exonerate

### filter design - code problem

Im trying to implement the legendre (optimum) filter described here.

the coefficient formulae come from a z-transform of the transfer function and seem to be correct (well at least theyre returning realistic numbers )

the problem seems to be coming from the actual filter itself - it jumps directly to infinity when the feedback is introduced

Code: Select all
streamin in;
streamout out;
streamin k;
float k2, k3;

float z0,z1,z2,z3;
float a0,a1,a2,a3;
float b0,b1,b2,b3;

float c0P5=0.57735;float c1P31=1.3107;
float c1P35=1.35897;float c1P7=1.73205;

//precalc
k2=k*k;
k3=k2*k;

//calculate coefficients
a0=c0P5*k3;
a1=c1P7*k3;
a2=c1P7*k3;
a3=c0P5*k3;

b0=c0P5*k3 + c1P35*k2 + c1P31*k + 1;
b1=c1P7*k3 + c1P35*k2 - c1P31*k - 3;
b2=c1P7*k3 - c1P35*k2 - c1P31*k + 3;
b3=c0P5*k3 + c1P35*k2 + c1P31*k - 1;

//filter
z0 = in + b1*z1 - b2*z2 - b3*z3;
out= a0*z0 + a1*z1 + a2*z2 + a3*z3;

z3=z2;
z2=z1;
z1=z0;
firmament
essemer

Posts: 16
Joined: Mon May 23, 2011 9:07 am

### Re: filter design - code problem

invert the feedback?
Last edited by MegaHurtz on Fri Aug 26, 2011 1:50 pm, edited 1 time in total.
Visit my website at: http://www.schlukhash.nl

MegaHurtz
smaniac

Posts: 1515
Joined: Mon Aug 11, 2008 5:29 pm
Location: Eindhoven/Netherlands

### Re: filter design - code problem

MegaHurtz wrote:invert the feedback?

more specific ? i did try some switching of signs but didnt get anywhere

im starting to think it may be the coefficients, i tried it with arbitrary numbers (all 1s, all 0.5s etc) and the filter itself seems to be fine . .
firmament
essemer

Posts: 16
Joined: Mon May 23, 2011 9:07 am

### Re: filter design - code problem

Problem is I can recognise filters by sound,create and calculate them digitally and analog.. but i never related to filters in butterworth polynominals. To me an "optimal filter" is a LR4.
Visit my website at: http://www.schlukhash.nl

MegaHurtz
smaniac

Posts: 1515
Joined: Mon Aug 11, 2008 5:29 pm
Location: Eindhoven/Netherlands

### Re: filter design - code problem

firmament wrote:
Code: Select all
//filter
z0 = in + b1*z1 - b2*z2 - b3*z3;
out= a0*z0 + a1*z1 + a2*z2 + a3*z3;

I don't have time to look at the math right now, but I can tell you just from a quick glance that your transfer function is still in the frequency domain. It needs to be transformed into the time domain to work correctly.

-cyto

ps. if you leave your "c" values rounded to such high numbers (x.xxxx), I can almost guarantee that if you do get this thing working, it will be very unstable. Use all the bits you can!

cyto
essemilian

Posts: 317
Joined: Sun Nov 28, 2010 4:36 am
Location: CIN | OH | USA

### Re: filter design - code problem

have a look at the two pole example here:
https://ccrma.stanford.edu/~jos/fp/Direct_Form_II.html
if it is not compensated somewhere else in the coef calculations i guess the "in + b1*" should read "in - b1*"
i know some times people switch their a and b's in biquads. i have not yet read the paper you refered to so i cant tell atm in your case.

this is how i would write a DF2 filter:
Code: Select all
wn = in - wn1*a1 - wn2*a2 - wn3*a3;
out = wn*b0 + wn1*b1 + wn2*b2 +  wn3*b3;
wn3 = wn2;
wn2 = wn1;
wn1 = wn;
Any sufficiently advanced technology is indistinguishable from magic.
Arthur C. Clarke, "Profiles of The Future", 1961 (Clarke's third law)

http://www.audioteknikk.net

tor
essemilian

Posts: 462
Joined: Wed Apr 14, 2010 8:52 pm

### Re: filter design - code problem

Hello, I have scilab's s-z transform calculator, which transforms the H(s) to H(z)

8.414e-16 + 2.524e-15*z^-1 + 2.524e-15*z^-2 + 8.414e-16*z^-3
----------------------------------------------------------------------------------------
1 - 2.9999703*z^-1 + 2.9999406*z^-2 - 0.9999703*z^-3

where T (sample time) = 1/44100.

so,

Code: Select all
streamin in;
streamout out;

float x0, x1, x2, x3;
float y0, y1, y2, y3;

x0 = in;
y0 = 2.9999703*y1 - 2.9999406*y2 + 0.9999703*y3 + 8.414e-16*x0 + 2.524e-15*x1 + 2.524e-15*x2 + 8.414e-16*x3;
y3 = y2;
y2 = y1;
y1 = y0;
x3 = x2;
x2 = x1;
x1 = x0;
out = y0;

I think that the problem is maybe that SM's code component does not provide the type of DOUBLE float...
kaj
essemer

Posts: 40
Joined: Sun Jan 02, 2011 9:02 am

### Re: filter design - code problem

kaj wrote:I think that the problem is maybe that SM's code component does not provide the type of DOUBLE float...
Did you say you were going to vote on this thread?

viewtopic.php?f=5&t=3815
Last edited by infuzion on Sat Aug 27, 2011 11:17 pm, edited 1 time in total.
Need help? First search the forum & WiKi, then post in the help forum with a clear topic, request, & OSM. Then please WiKi the correct solution. If you want my personal assistance, I charge by the hour or for an exchange of services.
infuzion
smstar

Posts: 6163
Joined: Wed May 04, 2005 8:02 pm
Location: Earth, USA, CO, Denver

### Re: filter design - code problem

thanks for the help guys!

MegaHurtz wrote:Problem is I can recognise filters by sound,create and calculate them digitally and analog.. but i never related to filters in butterworth polynominals. To me an "optimal filter" is a LR4.

this uses legendre polynomials rather than butterworth ones. theres already code around for l-r filters.

cyto wrote:I don't have time to look at the math right now, but I can tell you just from a quick glance that your transfer function is still in the frequency domain. It needs to be transformed into the time domain to work correctly.

the z values are delays, its configured like this

except with another delay in. pretty sure ive translated it to time domain - mathematically i took the H(s) transfer function and use thed bilinear transform to get discrete time H(z). i read this could be used directly for filter coefficients? or do i need a laplace transform for impulse response or something

cyto wrote:ps. if you leave your "c" values rounded to such high numbers (x.xxxx), I can almost guarantee that if you do get this thing working, it will be very unstable. Use all the bits you can!

thanks haha ive fixed this. ive also corrected one of the signs in my calculations for the b coefficients (interestingly they show some sort of symmetry) and the filter corner is now at the frequency i set. the response also looks about right - however its upside down / rotated somehow . . any ideas?
firmament
essemer

Posts: 16
Joined: Mon May 23, 2011 9:07 am

### Re: filter design - code problem

tor wrote:have a look at the two pole example here:
https://ccrma.stanford.edu/~jos/fp/Direct_Form_II.html
if it is not compensated somewhere else in the coef calculations i guess the "in + b1*" should read "in - b1*"
i know some times people switch their a and b's in biquads. i have not yet read the paper you refered to so i cant tell atm in your case.

this is how i would write a DF2 filter:
Code: Select all
wn = in - wn1*a1 - wn2*a2 - wn3*a3;
out = wn*b0 + wn1*b1 + wn2*b2 +  wn3*b3;
wn3 = wn2;
wn2 = wn1;
wn1 = wn;

thanks, ive changed the "in + a1" to a - and my code matches up with yours. the link is useful as well, this stuff is making more and more sense!
kaj wrote:Hello, I have scilab's s-z transform calculator, which transforms the H(s) to H(z)

8.414e-16 + 2.524e-15*z^-1 + 2.524e-15*z^-2 + 8.414e-16*z^-3
----------------------------------------------------------------------------------------
1 - 2.9999703*z^-1 + 2.9999406*z^-2 - 0.9999703*z^-3

where T (sample time) = 1/44100.

no idea whether this is right or not but id say its fixed cutoff and samplerate so i cant do much with it . . what cutoff frequency did you give it? i could check whether my coeffs are giving the same results . .
kaj wrote:so,

Code: Select all
streamin in;
streamout out;

float x0, x1, x2, x3;
float y0, y1, y2, y3;

x0 = in;
y0 = 2.9999703*y1 - 2.9999406*y2 + 0.9999703*y3 + 8.414e-16*x0 + 2.524e-15*x1 + 2.524e-15*x2 + 8.414e-16*x3;
y3 = y2;
y2 = y1;
y1 = y0;
x3 = x2;
x2 = x1;
x1 = x0;
out = y0;

I think that the problem is maybe that SM's code component does not provide the type of DOUBLE float...

yeah this isnt working, sm cant seem to handle the 8.414e-16 numbers.

wait isnt there a double primitive in the toolbox? does that get full 64bit precision?
firmament
essemer

Posts: 16
Joined: Mon May 23, 2011 9:07 am

### Re: filter design - code problem

to firmament

Let me know the transfer function H(s) that you used.
kaj
essemer

Posts: 40
Joined: Sun Jan 02, 2011 9:02 am

### Re: filter design - code problem

1/C*F*PI2

So.. its getting integrated at the plus signs?
Visit my website at: http://www.schlukhash.nl

MegaHurtz
smaniac

Posts: 1515
Joined: Mon Aug 11, 2008 5:29 pm
Location: Eindhoven/Netherlands

### Re: filter design - code problem

firmament wrote:wait isnt there a double primitive in the toolbox? does that get full 64bit precision?

you can do plus minus and multiplication in 64 bit but you can only feed them with 32 bit values. i designed a biquad that takes advantage of this. you find it in the RBJ ported thread. though it is not truly 64bit in the whole filter due to the delay part of the scematic in 32 bit, the multiplications of the feedbacked signals results in more presicion in the bass and less rounding noise in general. problem is when you want to any thing but plus, minus & multiplication you have to use the stream2double and double2stream prims. i do belive this is not so cpu friendly and it prevent you for having a end to end 64bit design
Any sufficiently advanced technology is indistinguishable from magic.
Arthur C. Clarke, "Profiles of The Future", 1961 (Clarke's third law)

http://www.audioteknikk.net

tor
essemilian

Posts: 462
Joined: Wed Apr 14, 2010 8:52 pm

### Re: filter design - code problem

i made a quick 3-pole version for you to test. I have not tested it myself. But im pretty sure it should work right. Now you should also do the precalculations and coeffisient calculations in double to get even bettet results. At the cost of some more cpu of cource.

If you do coef maths in double you must delete the stream2double prims for the coefs in my scematic. but i guess you would understand that

3poleDF2double.osm
Any sufficiently advanced technology is indistinguishable from magic.
Arthur C. Clarke, "Profiles of The Future", 1961 (Clarke's third law)

http://www.audioteknikk.net

tor
essemilian

Posts: 462
Joined: Wed Apr 14, 2010 8:52 pm

### Re: filter design - code problem

kaj wrote:Let me know the transfer function H(s) that you used.

tfn.png (14.52 KiB) Viewed 1080 times

see what the calculator comes up with
MegaHurtz wrote:1/C*F*PI2

So.. its getting integrated at the plus signs?

sorry not really sure what you mean here ? i have a similar part in the schematic [ tan(pi*f/sr) ] to calculate the k coefficient in the code, its to do with making the frequency setting samplerate independent
tor wrote:you can do plus minus and multiplication in 64 bit but you can only feed them with 32 bit values. i designed a biquad that takes advantage of this. you find it in the RBJ ported thread. though it is not truly 64bit in the whole filter due to the delay part of the scematic in 32 bit, the multiplications of the feedbacked signals results in more presicion in the bass and less rounding noise in general. problem is when you want to any thing but plus, minus & multiplication you have to use the stream2double and double2stream prims. i do belive this is not so cpu friendly and it prevent you for having a end to end 64bit design

the 64bit stuff seems . . useless? if you can only feed it 32 bit values then round them back off to 32 bit at the end . .
firmament
essemer

Posts: 16
Joined: Mon May 23, 2011 9:07 am

Next