JSFX Butterworth HPF noises problem

[timboid]

I made this small JSFX mostly to entertain myself. It works reasonably well, although with my limited knowledge of DSP-math I couldn't figure out how to calculate the recurrence coefficients, so instead, I got them calculated by the awesome http://www-users.cs.york.ac.uk/~fish...-bin/mkfscript and then run polynomial approximation (and in fact, I borrowed the code structure too...). Amazingly it works, the precision of math in Jesusonic is outstanding!
But I've run into problems with plugin generating LOUD noises when I change filter frequency. Is there a way to alleviate that?

Here's the code, but BE WARNED IF YOU TRY IT OUT!!!!! IT CAN POSSIBLY BLOW UP YOUR SPEAKERS!!!!

desc: 3rd order HPF

/////WARNING!!!! LOUD NOISES WHEN CHANGING FREQUENCY!!!

slider1:20<20,400,20>HPF (Hz)

@init
  cDcAdd = 10^-30;
  cDenorm = 10^-30;

@slider
 hpfreq = slider1*96000/srate;
 polyA1 = -1.30899690866119/10^4;
 polyA2 = 8.56732268461399/10^9;
 polyA3 = -3.96942694457829/(10^13);
 polyA4 = 1.46879322765979/(10^17);
 yc0 = 1 + polyA1*hpfreq + polyA2*(hpfreq^2) + polyA3*(hpfreq^3) + polyA4*(hpfreq^4);
 polyB1 = 2.61799388802009/10^4;
 polyB2 = -8.56738348538136/10^9;
 polyB3 = 9.35691781968719/10^14;
 polyB4 = 2.77054841032207/10^18;
 yc1 =-3 + polyB1*hpfreq + polyB2*(hpfreq^2) + polyB3*(hpfreq^3) + polyB4*(hpfreq^4);
 polyC1 = -1.30899694128233/10^4;
 polyC2 = 7.29823379873708/10^15;
 polyC3 = 2.33329646402345/10^14;
 yc2 = 3 + polyC1*hpfreq + polyC2*(hpfreq^2) + polyC3*(hpfreq^3);
 polyD1 = 6.54498827934359/10^5;
 polyD2 = 2.14147812558925/10^9;
 polyD3 = 5.97209136623587/10^14;
 GHPF = 1 + polyD1*hpfreq + polyD2*(hpfreq^2) + polyD3*(hpfreq^3);
 invGHPF = 1/GHPF;

@sample
  xv00 = xv10;
  xv10 = xv20;
  xv20 = xv30;
  xv30 = spl0 * invGHPF;
  yv00 = yv10;
  yv10 = yv20;
  yv20 = yv30;
  yv30 = (xv30 - xv00) + 3*(xv10 - xv20) + yc0*yv00 + yc1*yv10 + yc2*yv20;
  spl0 = yv30;
 
  xv01 = xv11;
  xv11 = xv21;
  xv21 = xv31;
  xv31 = spl1 * invGHPF;
  yv01 = yv11;
  yv11 = yv21;
  yv21 = yv31;
  yv31 = (xv31 - xv01) + 3*(xv11 - xv21) + yc0*yv01 + yc1*yv11 + yc2*yv21;
  spl1 = yv31;

[SaulT]

Very cool!

Some filter structures are very sensitive to rapid changes in coefficients (I'm thinking of an all-pass or two, but whatever). Basically either choose an algorithm that doesn't have this issue (e.g. the tpt zdf) or interpolate your coefficients, that is, change them slowly from one value to the other.

I went ahead and implemented the interpolation using single pole filters. This way, the processing variables (e.g. yc0, yc1) are set to the initial defaults, then when an updated value is calculated they slide to the new value. For me it looked like the fluctuations need at least a 5 ms time constant to be manageable, you can pick a value and hard code it, of course.

Can you speak a little more specifically as to how you used the polynomial approximation to get these coefficients?

EDIT: I think I see. Using a polynomial approximation to ballpark the coefficient values that would go with certain frequencies. Very nifty!

desc: timboid - 3rd order HPF

slider1:200<20,400,20>HPF (Hz)
slider2:20<1,100,1>glide (ms)

@init
  cDcAdd = 10^-30;
  cDenorm = 10^-30;

// since these don't change they can go in @init

 polyA1 = -1.30899690866119/10^4;
 polyA2 = 8.56732268461399/10^9;
 polyA3 = -3.96942694457829/(10^13);
 polyA4 = 1.46879322765979/(10^17);
 polyB1 = 2.61799388802009/10^4;
 polyB2 = -8.56738348538136/10^9;
 polyB3 = 9.35691781968719/10^14;
 polyB4 = 2.77054841032207/10^18;
 polyC1 = -1.30899694128233/10^4;
 polyC2 = 7.29823379873708/10^15;
 polyC3 = 2.33329646402345/10^14;
 polyD1 = 6.54498827934359/10^5;
 polyD2 = 2.14147812558925/10^9;
 polyD3 = 5.97209136623587/10^14;

function singlepole(in,target,coeff) ( in*coeff + target*(1-coeff); );

function calculate_var() (
 nhpfreq = slider1*96000/srate;
 nyc0 = 1 + polyA1*nhpfreq + polyA2*(nhpfreq^2) + polyA3*(nhpfreq^3) + polyA4*(nhpfreq^4);
 nyc1 =-3 + polyB1*nhpfreq + polyB2*(nhpfreq^2) + polyB3*(nhpfreq^3) + polyB4*(nhpfreq^4);
 nyc2 = 3 + polyC1*nhpfreq + polyC2*(nhpfreq^2) + polyC3*(nhpfreq^3);
 nGHPF = 1 + polyD1*nhpfreq + polyD2*(nhpfreq^2) + polyD3*(nhpfreq^3);
 ninvGHPF = 1/nGHPF;
  );

function update_var() (
 yc0 = singlepole(yc0,nyc0,glide_coeff);
 yc1 = singlepole(yc1,nyc1,glide_coeff);
 yc2 = singlepole(yc2,nyc2,glide_coeff);
 GHPF = singlepole(GHPF,nGHPF,glide_coeff);
 invGHPF = singlepole(invGHPF,ninvGHPF,glide_coeff);
 );

// starting with 0's in coefficients creates big spike when plugin is
// first started, so let's put some starting values in the processing vars

 calculate_var();
 hpfreq = nhpfreq;
 yc0 = nyc0;
 yc1 = nyc1;
 yc2 = nyc2;
 GHPF = nGHPF;
 invGHPF = ninvGHPF;

@slider
 glide_coeff = exp(-1/(0.001 * slider2 * srate));
 calculate_var();

@sample
  update_var();

  xv00 = xv10;
  xv10 = xv20;
  xv20 = xv30;
  xv30 = spl0 * invGHPF;
  yv00 = yv10;
  yv10 = yv20;
  yv20 = yv30;
  yv30 = (xv30 - xv00) + 3*(xv10 - xv20) + yc0*yv00 + yc1*yv10 + yc2*yv20;
  spl0 = yv30;
 
  xv01 = xv11;
  xv11 = xv21;
  xv21 = xv31;
  xv31 = spl1 * invGHPF;
  yv01 = yv11;
  yv11 = yv21;
  yv21 = yv31;
  yv31 = (xv31 - xv01) + 3*(xv11 - xv21) + yc0*yv01 + yc1*yv11 + yc2*yv21;
  spl1 = yv31;

  // this keeps your output from blowing up (speaker protection)

  spl0 = min(max(-1,spl0),1);
  spl1 = min(max(-1,spl1),1);

[timboid]

Whoa, thanks man! Works great!

Originally Posted by SaulT

Can you speak a little more specifically as to how you used the polynomial approximation to get these coefficients?

EDIT: I think I see. Using a polynomial approximation to ballpark the coefficient values that would go with certain frequencies. Very nifty!

Yep, you are correct, I generated a number of coefficients via webscript and got them approximate with a polynomial. I thought of having a LUT, but that would be valid only for a particular samplerate. I could do this for most used ones, but maybe someone would wish to do something crazy, so having a polynomial felt better. I think I may have to add a couple of points below 20Hz, to get better behavior outside of the range I initially planned.
I will post the update here if it does something useful.


EDIT: OK, here it is, I approximated a wider range of coefficients, so more samplerates are covered.

desc: timboid - 3rd order HPF

slider1:200<20,400,20>HPF (Hz)
slider2:20<1,100,1>glide (ms)

@init
  cDcAdd = 10^-30;
  cDenorm = 10^-30;

// since these don't change they can go in @init

 polyA1 = -1.30899694187114/10^4;
 polyA2 = 8.56736585100884/10^9;
 polyA3 = -3.97180419540762/(10^13);
 polyA4 = 1.52635392228233/(10^17);
 polyA5 = -5.04700087875157/10^22;
 polyB1 = 2.6179938821757/10^4;
 polyB2 = -8.5673694152914/10^9;
 polyB3 = 9.3473089577013/10^14;
 polyB4 = 3.02591121470417/10^18;
 polyB5 = -2.32200010482097/10^22;
 polyC1 = -1.30899688523548/10^4;
 polyC2 = -4.81938237476238/10^14;
 polyC3 = 2.35163494605937/10^14;
 polyC4 = -1.97596635235226/10^19;
 polyD1 = 6.54498438435183/10^5;
 polyD2 = 2.14186778172279/10^9;
 polyD3 = 5.83223951666622/10^14;
 polyD4 = 1.64188114174327/10^18;


function singlepole(in,target,coeff) ( in*coeff + target*(1-coeff); );

function calculate_var() (
 nhpfreq = slider1*96000/srate;
 nyc0 = 1 + polyA1*nhpfreq + polyA2*(nhpfreq^2) + polyA3*(nhpfreq^3) + polyA4*(nhpfreq^4) + polyA5*(nhpfreq^5);
 nyc1 =-3 + polyB1*nhpfreq + polyB2*(nhpfreq^2) + polyB3*(nhpfreq^3) + polyB4*(nhpfreq^4) + polyB5*(nhpfreq^5);
 nyc2 = 3 + polyC1*nhpfreq + polyC2*(nhpfreq^2) + polyC3*(nhpfreq^3) + polyC4*(nhpfreq^4);
 nGHPF = 1 + polyD1*nhpfreq + polyD2*(nhpfreq^2) + polyD3*(nhpfreq^3) + polyD4*(nhpfreq^4);
 ninvGHPF = 1/nGHPF;
  );

function update_var() (
 yc0 = singlepole(yc0,nyc0,glide_coeff);
 yc1 = singlepole(yc1,nyc1,glide_coeff);
 yc2 = singlepole(yc2,nyc2,glide_coeff);
 GHPF = singlepole(GHPF,nGHPF,glide_coeff);
 invGHPF = singlepole(invGHPF,ninvGHPF,glide_coeff);
 );

// starting with 0's in coefficients creates big spike when plugin is
// first started, so let's put some starting values in the processing vars

 calculate_var();
 hpfreq = nhpfreq;
 yc0 = nyc0;
 yc1 = nyc1;
 yc2 = nyc2;
 GHPF = nGHPF;
 invGHPF = ninvGHPF;

@slider
 glide_coeff = exp(-1/(0.001 * slider2 * srate));
 calculate_var();

@sample
  update_var();

  xv00 = xv10;
  xv10 = xv20;
  xv20 = xv30;
  xv30 = abs(spl0) < cDenorm ? 0 : spl0;
  xv30 *= invGHPF;
  yv00 = yv10;
  yv10 = yv20;
  yv20 = yv30;
  yv30 = (xv30 - xv00) + 3*(xv10 - xv20) + yc0*yv00 + yc1*yv10 + yc2*yv20;
  spl0 = yv30;
 
  xv01 = xv11;
  xv11 = xv21;
  xv21 = xv31;
  xv31 = abs(spl1) < cDenorm ? 0 : spl1;
  xv31 *= invGHPF;
  yv01 = yv11;
  yv11 = yv21;
  yv21 = yv31;
  yv31 = (xv31 - xv01) + 3*(xv11 - xv21) + yc0*yv01 + yc1*yv11 + yc2*yv21;
  spl1 = yv31;

  spl0 += cDcAdd;
  spl1 += cDcAdd;

  spl0 = min(max(-1,spl0),1);
  spl1 = min(max(-1,spl1),1);