Calculating Frequency Response of BiQuad filters

Development
1

So, on a new module I’m working on I want to display its frequency response. It is seven filters that can be in any combination of parallel and serial. That part of the display I have working.

The issue is the underlying frequency response of the the Biquad filters I am using. I’m using the standard BiQuad filter code that rack offers, and I found this code for calculating the frequency response:

float w = 2.0*M_PI*frequency;  
float numerator = a0*a0 + a1*a1 + a2*a2 + 2.0*(a0*a1 + a1*a2)*cos(w) + 2.0*a0*a2*cos(2.0*w);
float denominator = 1.0 + b1*b1 + b2*b2 + 2.0*(b1 + b1*b2)*cos(w) + 2.0*b2*cos(2.0*w);
float magnitude = sqrt(numerator / denominator);

The issue is that the above code works great except for calculating the response for frequencies below about 100hz - and the Q makes a big difference. Above 100hz or so, the graph looks great and just like what you would expect, but it gets real noisy below.

Any ideas what I’m doing wrong here? The filters themselves sound great…just the display

1 Like
2

so turns out issue was one of precision. I read the BiQuads are very sensitive to the values of the coefficients - I switched from floats to doubles and the problem went away. Of course now I can’t use simd until it supports double_4s or whatever their called :confused:

3

The formula for the amplitude of a given frequency is |H(z)| with H defined at https://en.wikipedia.org/wiki/Digital_biquad_filter and z = e^{i 2\pi f / f_n}. If you’re getting low numerical stability due to the numerator and denomenator being close to zero, maybe you can approach it like a limit problem and apply L’Hôpital’s rule? Idk. Almost always there’s a way to rewrite a formula as simple as that to be numerically stable.