velipso / sndfilter

Algorithms for sound filters, like reverb, dynamic range compression, lowpass, highpass, notch, etc
BSD Zero Clause License
445 stars 70 forks source link

Low pass filter infinite alpha #9

Closed neevek closed 3 years ago

neevek commented 5 years ago

In the biquad_makeLPF function in reverb.c:

static inline void biquad_makeLPF(sf_rv_biquad_st *biquad, int rate, float freq, float bw){
  freq = clampf(freq, 0, rate / 2);
  float omega = 2.0f * (float)M_PI * freq / (float)rate;
  float cs = cosf(omega);
  float sn = sinf(omega);
  float alpha = sn * sinhf((float)M_LN2 * 0.5f * bw * omega / sn);
  float a0inv = 1.0f / (1.0f + alpha);
  biquad->b0 = a0inv * (1.0f - cs) * 0.5f;
  biquad->b1 = 2.0f * biquad->b0;
  biquad->b2 = biquad->b0;
  biquad->a1 = a0inv * -2.0f * cs;
  biquad->a2 = a0inv * (1.0f - alpha);
  biquad->xn1 = 0;
  biquad->xn2 = 0;
  biquad->yn1 = 0;
  biquad->yn2 = 0;
}

If rate is integral multiple of freq, then omega will be integral multiple of PI, and sn = sinf(M_PI) is 0, sn is used as a divisor for calculating alpha, the result of alpha will be INFINITY, which causes any other values to be NaN when calculated with alpha.

Can I test the value of alpha with isinf(), if it is true, assign a valid value to it? But what is the valid value range of alpha? Or what is the valid value for it if it is INFINITY?

velipso commented 5 years ago

If rate is integral multiple of freq

It isn't for most values of freq... but let's assume it is.

If freq is set to 441Hz, and rate is set to 44100 samples per second, then omega will be:

2 PI 441 / 44100 = 0.06283185307179585

Then sn and cs will be:

Sin(0.06283185307179585) = 0.06279051952931336 Cos(0.06283185307179585) = 0.9980267284282716

That being said, now that I'm looking at the code, if freq <= 0 or >= rate/2, then it looks like your logic is right -- sn will be 0, which gets a divide by zero. If I trace the sources of freq, I can see that a user could pass in a 0 for basslpf in the advanced options, and cause that to happen.

So I think the clamp should probably be fixed.

neevek commented 5 years ago

I just run into the case that freq >= rate/2 with default settings with all presets other than LARGEHALL1 and LARGEHALL2. With the same input, only LARGEHALL1 and LARGEHALL2 work, all other presets will zero out the samples because of the divide by zero issue.

velipso commented 5 years ago

I don't have time to fix at the moment, but a quick fix is to simply change the biquad coefficients to be a passthrough. In order to do that, you set b0 to 1, and the rest (b1, b2, a1, a2) to 0.

neevek commented 5 years ago

I don't have time to fix at the moment, but a quick fix is to simply change the biquad coefficients to be a passthrough. In order to do that, you set b0 to 1, and the rest (b1, b2, a1, a2) to 0.

I set alpha to 1 when sn is 0, which seemed to work as expected.

velipso commented 3 years ago

thanks, sorry it took so long to fix, I don't check here often but I had some free time today