filter design - code problem

[firmament]

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;

[MegaHurtz]

invert the feedback?

[firmament]

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 . .

[MegaHurtz]

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.

[cyto]

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!

[tor]

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;

[kaj]

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...

[infuzion]

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

[firmament]

thanks for the help guys!

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.

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

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]

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!

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 . .

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?

[kaj]

to firmament

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

[MegaHurtz]

1/C*F*PI2

So.. its getting integrated at the plus signs?

[tor]

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

[tor]

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

(780 Bytes) Downloaded 88 times

[firmament]

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

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

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

I'll try and find the ported rbj thread as the example you gave me isnt loading.
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 . .