Discrete filters by sin substitution

[enslaver101]

AUTO-ADMIN: Non-MP3, WAV, OGG, SoundCloud, YouTube, Vimeo, Twitter and Facebook links in this post have been protected automatically. Once the member reaches 5 posts the links will function as normal.

Here is a method to create a digital filter from an analog prototype, by replacing the s term with j*sin(x*pi/2)/(w*pi/2) and solving 2 quadratic equations to get the coefficients a1 and a2, then for the lowpass we put (1 + a1 + a2) in the b0 coefficient. The technique is similar to the bilinear transform, except that with BLT the coefficients are pretty simple, whereas here you need to calculate a nested square root and the coefficients are prone to numerical precision loss.
Here's some code and a Desmos graph. The blue line represents our desired analog filter, the other one is the approximation.

https://www.desmos.com/calculator/xhqqysafgx (https://www.desmos.com/calculator/xhqqysafgx)

Code: Select all (#)

Coefficient calculations:

float w = 2*pi * freq / samplerate;

float ww = w*w;
float asd = (w - 2)*(w + 2);
float wq2 = 2*w/q;
float alpha = sqrt(asd*asd + wq2*wq2);
	
float a = 0.5 * (ww - alpha);
float b = 0.5 * (ww + alpha);
float c = 0.5 * (b - sqrt((b - 2)*(b + 2)));

float a1 = a*c;
float a2 = c*c;

float b0 = 1 + a1 + a2;
float b1 = 0, b2 = 0;
if (false) //bandpass
{
	b0 = b0/w;
	b1 = -b0;
	b2 = 0;
}
else if (false) //hp
{
	b0 = b0/(w*w);
	b1 = -2*b0;
	b2 = b0;
}

Somewhere in the processing loop:

float x = in[chan][i];
float tmp = x*b0 + x1*b1 + x2*b2 - a1*y1 - a2*y2;
x2 = x1; x1 = x; y2 = y1; y1 = tmp;
out[chan][i] = tmp;

A big problem is that the lower the filter frequency and the higher the Q are, the worse it performs. Using doubles helps (although if you go even lower it'll probably still blow up), but the BLT, for example does not even need double precision. If anyone has an idea how to get these coefficients in a better way, it would be appreciated.

[mystran]

A big problem is that the lower the filter frequency and the higher the Q are, the worse it performs.

This applies to (somewhat) to direct form filters in general, no matter how you compute the coefficients... because for a biquad P(z)/(1+a1*z^-1+a2*z^-2) the a1 coefficient is -2*cos(w) [edit: well, technically -2*cos(w)/a0, with a0=1+sin(w)/(2*Q), but that doesn't help] where w is the angular pole frequency.. and cos is very close to 1 for small values of w, which means precision goes bye-bye...

...but in your formula with small w, at least (w-2)*(w+2) loses a bunch of precision, because 2 is large in comparison and if we expand this to w*w-4 that's even worse, because w*w is smaller and 4 is bigger. Then we square that value, we lose even more of what's left of w.

That seems like the biggest loss of precision. Essentially you're computing sqrt(16 + wq2*wq2) and now 2*w/q is also tiny, so wq2*wq2 is tinier... so we have sqrt(16 + some rounding errors).

I'm not convinced that I understand what you're actually computing, so I'm not going to try and suggest how to compute it more accurately... but like it's probably the "alpha" term where things go bad.

[enslaver101]

AUTO-ADMIN: Non-MP3, WAV, OGG, SoundCloud, YouTube, Vimeo, Twitter and Facebook links in this post have been protected automatically. Once the member reaches 5 posts the links will function as normal.

This version will not blow up at low frequencies:

Code: Select all (#)

float w = (pi * freq / samplerate);

float ww = w*w;
float asd = (w - 1)*(w + 1);
float wq = w/q,  wq2 = wq*wq;
float alpha = sqrt(asd*asd + wq2);

float a = ww - alpha;
//float b = asd + alpha;
float b = -wq2/(a - 1);
float c = 1 + b - sqrt((b)*(b + 2));

a1 = 2*a*c;
a2 = c*c;

('b' is calculated using alternative quadratic formula c/(-b - sqrt(b^2 - 4ac))
Now it gets more inaccurate towards the high end (unnoticeable - I was printfing the numbers at this time, so I could tell). If set to 22579200hz, it blows up (this boundary lowers as Q is increased). Also there's an additional division, which is undesirable.

it's probably the "alpha" term where things go bad.

It definitely is, at least part of it. But pretty much any sort of factoring I applied made things worse, which was disappointing. Ideally I'd like not to have sqrt at all in there.

Turns out it's possible to use other formulas, with more frequency warping - prewarp w so that it cancels the square root. For example, this formula results in a similar filter (basically it's the same filter, but shifted; the shift depends on Q so it makes things much less predictable):

Code: Select all (#)

float w = (pi * freq / samplerate);

float ww = w*w;
float c = 1 - w/q;
a1 = 2*c*(ww*(1 - 1/(2*q*q)) - c)/(c + ww);
a2 = c*c;

The additional warping may be compensated in some way; I haven't experimented with it so far.