Prefer SVF over DF1/DF2 biquads

Thanks for the pointer. I will check this out, because it's very much the kind of core low-level thing I'm trying to learn about.

And I did do DF2 initially because I was lured in by the "canonical" part and the smaller state. The numerical instability I learnt the hard way, although it was probably triggered by me doing other stupid things. But I'll check out those SVF filters as an alternative to DF1.

And no, I'm not interested in mapping the raw DF1 coefficients to SVF ones - I'd much rather do the coefficients based on the filter specs, since that's how I generate all the DF1 ones anyway. So I'll take a look at the papers.

From a very quick look, it looked like your default coefficients were based on tan(𝜔) in the core papers, but there is one pdf about using just sin() instead. It looks like the preferred form to me (because I'm lazy and I already have 'sin()' in my microcontroller library). That paper seems to also say that the sin version is more accurate with low precision (which is what I have), so it seems a win-win situation, but maybe I misunderstand from my quick scan.

Thanks for the pointer. I will check this out, because it's very much the kind of core low-level thing I'm trying to learn about.

And I did do DF2 initially because I was lured in by the "canonical" part and the smaller state. The numerical instability I learnt the hard way, although it was probably triggered by me doing other stupid things. But I'll check out those SVF filters as an alternative to DF1.

And no, I'm not interested in mapping the raw DF1 coefficients to SVF ones - I'd much rather do the coefficients based on the filter specs, since that's how I generate all the DF1 ones anyway. So I'll take a look at the papers.

From a very quick look, it looked like your default coefficients were based on tan(𝜔) in the core papers, but there is one pdf about using just sin() instead. It looks like the preferred form to me (because I'm lazy and I already have 'sin()' in my microcontroller library). That paper seems to also say that the sin version is more accurate with low precision (which is what I have), so it seems a win-win situation, but maybe I misunderstand from my quick scan.

Thanks for the pointer. I will check this out, because it's very much the kind of core low-level thing I'm trying to learn about.

And I did do DF2 initially because I was lured in by the "canonical" part and the smaller state. The numerical instability I learnt the hard way, although it was probably triggered by me doing other stupid things. But I'll check out those SVF filters as an alternative to DF1.

And no, I'm not interested in mapping the raw DF1 coefficients to SVF ones - I'd much rather do the coefficients based on the filter specs, since that's how I generate all the DF1 ones anyway. So I'll take a look at the papers.

From a very quick look, it looked like your default coefficients were based on tan(𝜔) in the core papers, but there is one pdf about using just sin() instead. It looks like the preferred form to me (because I'm lazy and I already have 'sin()' in my microcontroller library). That paper seems to also say that the sin version is more accurate with low precision (which is what I have), so it seems a win-win situation, but maybe I misunderstand from my quick scan.

Thanks for the pointer. I will check this out, because it's very much the kind of core low-level thing I'm trying to learn about.

And I did do DF2 initially because I was lured in by the "canonical" part and the smaller state. The numerical instability I learnt the hard way, although it was probably triggered by me doing other stupid things. But I'll check out those SVF filters as an alternative to DF1.

And no, I'm not interested in mapping the raw DF1 coefficients to SVF ones - I'd much rather do the coefficients based on the filter specs, since that's how I generate all the DF1 ones anyway. So I'll take a look at the papers.

From a very quick look, it looked like your default coefficients were based on tan(𝜔) in the core papers, but there is one pdf about using just sin() instead. It looks like the preferred form to me (because I'm lazy and I already have 'sin()' in my microcontroller library). That paper seems to also say that the sin version is more accurate with low precision (which is what I have), so it seems a win-win situation, but maybe I misunderstand from my quick scan.

Thanks for the pointer. I will check this out, because it's very much the kind of core low-level thing I'm trying to learn about.

And I did do DF2 initially because I was lured in by the "canonical" part and the smaller state. The numerical instability I learnt the hard way, although it was probably triggered by me doing other stupid things. But I'll check out those SVF filters as an alternative to DF1.

And no, I'm not interested in mapping the raw DF1 coefficients to SVF ones - I'd much rather do the coefficients based on the filter specs, since that's how I generate all the DF1 ones anyway. So I'll take a look at the papers.

From a very quick look, it looked like your default coefficients were based on tan(𝜔) in the core papers, but there is one pdf about using just sin() instead. It looks like the preferred form to me (because I'm lazy and I already have 'sin()' in my microcontroller library). That paper seems to also say that the sin version is more accurate with low precision (which is what I have), so it seems a win-win situation, but maybe I misunderstand from my quick scan.

Can't say I understand this well, but here's code from my discussion with Gemini about filters and it uses a sine coefficient.

#include <stdio.h>
#include <math.h>

typedef struct {
    float low, band; // The "State" variables
    float f, q;      // Coefficients
} SVF;

// Setup: freq is 0.0 to 1.0 (normalized to Nyquist), res is 0.0 to 1.0
void svf_init(SVF *s, float freq, float res) {
    s->low = 0;
    s->band = 0;
    // Simple tuning coefficients
    s->f = 2.0f * sinf(M_PI * freq); 
    s->q = 1.0f - res;
}

void svf_process(SVF *s, float input, float *lp, float *hp, float *bp) {
    // The "Magic" of the SVF: Three lines of code for three filters
    float high = input - s->low - s->q * s->band;
    s->band = s->f * high + s->band;
    s->low = s->f * s->band + s->low;

    *lp = s->low;
    *hp = high;
    *bp = s->band;
}

int main() {
    SVF filter;
    svf_init(&filter, 0.1f, 0.5f); // 10% of Nyquist, medium resonance

    float input = 1.0f; // A single impulse
    float lp, hp, bp;

    printf("Step | Lowpass | Bandpass | Highpass\n");
    printf("-----|---------|----------|----------\n");

    for (int i = 0; i < 10; i++) {
        svf_process(&filter, input, &lp, &hp, &bp);
        input = 0; // Turn off input after first sample (Impulse Response)
        printf(" %d   |  %.3f  |  %.3f   |  %.3f\n", i, lp, bp, hp);
    }

    return 0;
}

https://g.co/gemini/share/222287d499f5

It's great to see your interested in audio circuits and DSP :) It's a fascinating field, and one I'm lucky enough to be able to make my full time job. It would be my pleasure to help out with any questions you have on various low level DSP type topics. I've got a lot of experience in this area, and can hopefully save you a lot of time. Back to the SVF vs DF1...

In advance, sorry if I made some typos in the following, please consider this a summary and check my technical paper pdfs for the full details because that is where I've double checked the DSP in the Mathematica notebooks by running it and doing various plots.

Although the sin only coefficient version is a bit more accurate for the coefficient part, but the difference is pretty minimal, and even for "pro audio" the tan version is just fine (this is what I got Ableton to use for Live and their other products). The sin only version is ok for basic responses, but as soon as you throw in shelving shapes it gets a little uglier. For easier coefficient calculation, and the same state increment implementation, I recommend using the trapezoidal integrated SVF as outlined in this paper:

https://cytomic.com/files/dsp/SvfLinearTrapAllOutputs.pdf

It does use tan, but your ARM chip for this project looks ok for the some limited 32-bit float divisions, which you would have to do for coefficient calculation of the DF1 anyway. There are two divs all up for a "set" of the new cutoff frequency, one for the tan approximation, and one for the normalisation of the coefficients. Please feel free to use one of the following rational approximations for the tan I designed for use with this sort of SVF coefficient computation:

// tanh approximations for x in range +-1.54

// relative error of around +-1e-3, or +-2 cents
tan(x) ~= (1.00111 x - 0.075021 x^3)/(1 - 0.405097 x^2)

// relative error of around +-1e-7, or +-2e-4 cents
tan(x) ~= (x - 0.112218 x^3 + 0.00113748 x^5)/(1 - 0.445551 x^2 + 0.0163193 x^4)

The basic psuedo code for for the trapezoidal integrated SVF from the paper is as follows:

init:
  w = pi*cutoff/samplerate
  g0 = tan(w)
  k0 = damp = 1/Q = 2 - 2*resonance, for resonance in the range [0, 1)
  A = 10^(gaindb/40)

  switch response:
    case high:
      g = g0
      k = k0
      m0 = 1
      m1 = 0
      m2 = 0
    case band:
      g = g0
      k = k0
      m0 = 0
      m1 = 1
      m2 = 0
    case low:
      g = g0
      k = k0
      m0 = 0
      m1 = 0
      m2 = 1
    case notch:
      g = g0
      k = k0
      m0 = 1
      m1 = 0
      m2 = 1
    case peak:
      g = g0
      k = k0
      m0 = 1
      m1 = 0
      m2 = -1
    case all:
      g = g0
      k = k0
      m0 = 1
      m1 = -k0
      m2 = 1
    case bell:
      g = g0
      k = k0/A
      m0 = 1
      m1 = k0*A
      m2 = 1
    case low shelf:
      g = g0/sqrt(A)
      k = k0
      m0 = 1
      m1 = k0*A
      m2 = A*A
    case high shelf:
      g = g0*sqrt(A)
      k = k0
      m0 = A*A
      m1 = k0*A
      m2 = 1

  gk = g + k
  gt0 = 1/(1 + g*gk)
  gk0 = gk*gt0

  // linear interpolate these between blocks
  gt1 = g*gt0
  gk1 = g*gk0
  gt2 = g*gt1
  m0x = m0
  m1x = m1 - k*m0
  m2x = m2 - m0

clear:
  ic1eq = 0
  ic2eq = 0

tick:
  t0 = input - ic2eq
  t1 = gt1*t0 - gk1*ic1eq
  t2 = gt2*t0 + gt1*ic1eq
  v1 = t1 + ic1eq
  v2 = t2 + ic2eq
  ic1eq = ic1eq + 2*t1
  ic2eq = ic2eq + 2*t2
  output = m0x*input + m1x*v1 + m2x*v2

This form is great since you can do linear interpolation on the 6 computed coefficients {gt1, gk1, gt2, m0x, m1x, m2x}, and get a smooth filter in between these endpoints. This way your coefficient calcuation only has to be done once per block of say 32 samples, then you += on these terms each sample so after the block is done they have reached their respective endpoints.

I don't recommend using the Chamberlain SVF (quasi forward Euler SVF) as posted by @freealise as a replacement for a DF1 since it has very limited frequency range before being numerically unstable.

It's great to see you're interested in audio circuits and DSP :) It's a fascinating field, and one I'm lucky enough to be able to make my full time job. It would be my pleasure to help out with any questions you have on various low level DSP type topics. I've got a lot of experience in this area, and can hopefully save you a lot of time. Back to the SVF vs DF1...

In advance, sorry if I made some typos in the following, please consider this a summary and check my technical paper pdfs for the full details because that is where I've double checked the DSP in the Mathematica notebooks by running it and doing various plots.

Although the sin only coefficient version is a bit more accurate for the coefficient part, but the difference is pretty minimal, and even for "pro audio" the tan version is just fine (this is what I got Ableton to use for Live and their other products filter/eq, and I got Valve to use for their Source 2 audio engine filters/eq). The sin only version is ok for basic responses, but as soon as you throw in shelving shapes it gets a little uglier. For easier coefficient calculation, and the same state increment implementation, I recommend using the trapezoidal integrated SVF as outlined in this paper:

https://cytomic.com/files/dsp/SvfLinearTrapAllOutputs.pdf

It does use tan, but your ARM chip for this project looks ok for the some limited 32-bit float divisions, which you would have to do for coefficient calculation of the DF1 anyway. There are two divs all up for a "set" of the new cutoff frequency, one for the tan approximation, and one for the normalisation of the coefficients. Please feel free to use one of the following rational approximations for the tan I designed for use with this sort of SVF coefficient computation:

// tanh approximations for x in range +-1.54

// relative error of around +-1e-3, or +-2 cents
tan(x) ~= (1.00111 x - 0.075021 x^3)/(1 - 0.405097 x^2)

// relative error of around +-1e-7, or +-2e-4 cents
tan(x) ~= (x - 0.112218 x^3 + 0.00113748 x^5)/(1 - 0.445551 x^2 + 0.0163193 x^4)

The basic psuedo code for for the trapezoidal integrated SVF from the paper is as follows:

init:
  w = pi*cutoff/samplerate
  g0 = tan(w)
  k0 = damp = 1/Q = 2 - 2*resonance, for resonance in the range [0, 1)
  A = 10^(gaindb/40)

  switch response:
    case high:
      g = g0
      k = k0
      m0 = 1
      m1 = 0
      m2 = 0
    case band:
      g = g0
      k = k0
      m0 = 0
      m1 = 1
      m2 = 0
    case low:
      g = g0
      k = k0
      m0 = 0
      m1 = 0
      m2 = 1
    case notch:
      g = g0
      k = k0
      m0 = 1
      m1 = 0
      m2 = 1
    case peak:
      g = g0
      k = k0
      m0 = 1
      m1 = 0
      m2 = -1
    case all:
      g = g0
      k = k0
      m0 = 1
      m1 = -k0
      m2 = 1
    case bell:
      g = g0
      k = k0/A
      m0 = 1
      m1 = k0*A
      m2 = 1
    case low shelf:
      g = g0/sqrt(A)
      k = k0
      m0 = 1
      m1 = k0*A
      m2 = A*A
    case high shelf:
      g = g0*sqrt(A)
      k = k0
      m0 = A*A
      m1 = k0*A
      m2 = 1

  gk = g + k
  gt0 = 1/(1 + g*gk)
  gk0 = gk*gt0

  // linear interpolate these between blocks
  gt1 = g*gt0
  gk1 = g*gk0
  gt2 = g*gt1
  m0x = m0
  m1x = m1 - k*m0
  m2x = m2 - m0

clear:
  ic1eq = 0
  ic2eq = 0

tick:
  t0 = input - ic2eq
  t1 = gt1*t0 - gk1*ic1eq
  t2 = gt2*t0 + gt1*ic1eq
  v1 = t1 + ic1eq
  v2 = t2 + ic2eq
  ic1eq = ic1eq + 2*t1
  ic2eq = ic2eq + 2*t2
  output = m0x*input + m1x*v1 + m2x*v2

This form is great since you can do linear interpolation on the 6 computed coefficients {gt1, gk1, gt2, m0x, m1x, m2x}, and get a smooth filter in between these endpoints. This way your coefficient calcuation only has to be done once per block of say 32 samples, then you += on these terms each sample so after the block is done they have reached their respective endpoints.

I don't recommend using the Chamberlin SVF (quasi forward Euler SVF) as posted by @freealise as a replacement for a DF1 since it has very limited frequency range before being numerically unstable.

It's great to see you're interested in audio circuits and DSP :) It's a fascinating field, and one I'm lucky enough to be able to make my full time job. It would be my pleasure to help out with any questions you have on various low level DSP type topics. I've got a lot of experience in this area, and can hopefully save you a lot of time. Back to the SVF vs DF1...

In advance, sorry if I made some typos in the following, please consider this a summary and check my technical paper pdfs for the full details because that is where I've double checked the DSP in the Mathematica notebooks by running it and doing various plots.

Although the sin only coefficient version is a bit more accurate for the coefficient part, but the difference is pretty minimal, and even for "pro audio" the tan version is just fine (this is what I got Ableton to use for Live and their other products filter/eq, and I got Valve to use for their Source 2 audio engine filters/eq). The sin only version is ok for basic responses, but as soon as you throw in shelving shapes it gets a little uglier. For easier coefficient calculation, and the same state increment implementation, I recommend using the trapezoidal integrated SVF as outlined in this paper:

https://cytomic.com/files/dsp/SvfLinearTrapAllOutputs.pdf

It does use tan, but your ARM chip for this project looks ok for the some limited 32-bit float divisions, which you would have to do for coefficient calculation of the DF1 anyway. There are two divs all up for a "set" of the new cutoff frequency, one for the tan approximation, and one for the normalisation of the coefficients. Please feel free to use one of the following rational approximations for the tan I designed for use with this sort of SVF coefficient computation:

// tanh approximations for x in range +-1.54

// relative error of around +-1e-3, or +-2 cents
tan(x) ~= (1.00111 x - 0.075021 x^3)/(1 - 0.405097 x^2)

// relative error of around +-1e-7, or +-2e-4 cents
tan(x) ~= (x - 0.112218 x^3 + 0.00113748 x^5)/(1 - 0.445551 x^2 + 0.0163193 x^4)

The basic psuedo code for for the trapezoidal integrated SVF from the paper is as follows:

init:
  w = pi*cutoff/samplerate
  g0 = tan(w)
  k0 = damp = 1/Q = 2 - 2*resonance, for resonance in the range [0, 1)
  A = 10^(gaindb/40)

  switch response:
    case high:
      g = g0
      k = k0
      m0 = 1
      m1 = 0
      m2 = 0
    case band:
      g = g0
      k = k0
      m0 = 0
      m1 = 1
      m2 = 0
    case low:
      g = g0
      k = k0
      m0 = 0
      m1 = 0
      m2 = 1
    case notch:
      g = g0
      k = k0
      m0 = 1
      m1 = 0
      m2 = 1
    case peak:
      g = g0
      k = k0
      m0 = 1
      m1 = 0
      m2 = -1
    case all:
      g = g0
      k = k0
      m0 = 1
      m1 = -k0
      m2 = 1
    case bell:
      g = g0
      k = k0/A
      m0 = 1
      m1 = k0*A
      m2 = 1
    case low shelf:
      g = g0/sqrt(A)
      k = k0
      m0 = 1
      m1 = k0*A
      m2 = A*A
    case high shelf:
      g = g0*sqrt(A)
      k = k0
      m0 = A*A
      m1 = k0*A
      m2 = 1

  gk = g + k
  gt0 = 1/(1 + g*gk)
  gk0 = gk*gt0

  // linear interpolate these between blocks
  gt1 = g*gt0
  gk1 = g*gk0
  gt2 = g*gt1
  m0x = m0
  m1x = m1 - k*m0
  m2x = m2 - m0

clear:
  ic1eq = 0
  ic2eq = 0

tick:
  t0 = input - ic2eq
  t1 = gt1*t0 - gk1*ic1eq
  t2 = gt2*t0 + gt1*ic1eq
  v1 = t1 + ic1eq
  v2 = t2 + ic2eq
  ic1eq = ic1eq + 2*t1
  ic2eq = ic2eq + 2*t2
  output = m0x*input + m1x*v1 + m2x*v2

This form is great since you can do linear interpolation on the 6 computed coefficients {gt1, gk1, gt2, m0x, m1x, m2x}, and get a smooth filter in between these endpoints. This way your coefficient calcuation only has to be done once per block of say 32 samples, then you += on these terms each sample so after the block is done they have reached their respective endpoints.

I don't recommend using the Chamberlin SVF (quasi forward Euler SVF) as posted by @freealise as a replacement for a DF1 since it has very limited frequency range before being numerically unstable.

Wonderful manner of discussion, truly in the spirit of open source. Cheers everyone!

Wonderful manner of discussion, truly in the spirit of open source. Cheers everyone!

Fascinating :)

Fascinating :)

Fascinating :)

You might also want to try implementing the biquads in Transposed Direct Form II. This form is more numerically robust than normal direct forms and is as easy and as efficient to implement. Though, an SVF-based model still trumps in terms of numerical stability and will yield better improvements.

You might also want to try implementing the biquads in Transposed Direct Form II. This form is more numerically robust than normal direct forms and is as easy and as efficient to implement. Though, an SVF-based model still trumps in terms of numerical stability and will yield better improvements.

Isn't the ladder structure the usual go to for numerical accuracy in LF? I've not looked in quantifiying the amount of computation to recompute parameters on the fly I must admit.

Question to @andy-cytomic, in your SVF-vs-DF1.pdf, "lad" means ladder structure?

Question to @andy-cytomic, in your SVF-vs-DF1.pdf, "lad" means ladder structure?

yes, lad is a ladder form df1, and it does look pretty good in terms of noise performance in the LTI case, but for musical applications where the cutoff / resonance change the SVF is the best structure to use, especially for all fp32 computation.