Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

LowPassFilter coefficients incorrectly calculated #144

Closed
mzero opened this issue Jan 11, 2022 · 4 comments
Closed

LowPassFilter coefficients incorrectly calculated #144

mzero opened this issue Jan 11, 2022 · 4 comments

Comments

@mzero
Copy link

mzero commented Jan 11, 2022

The coefficients used in filter are incorrectly calculated.

The fb coefficient is calculated, effectively as:

fb = q + q * (1 - cutoff)

However, the proper equation, which is the comments at line 34, is:

fb = q + q / (1 - f)

The referenced algorithm at musicdsp.org correctly shows the division. The referenced fixed point code at dave's blog (pawfal.org) made the error of replacing the divide with a multiply. The error is noted in the comments to his post, and is now fixed the code in the post itself.

The effect of the error is that you can't get consistent resonance results with varying frequency.


The f coefficient is simply set to the input value cutoff, where cutoff maps the 8 bit fixed point value 0 to 255 to frequencies of 0 to ½ sample rate.

However, it turns out that f in the DSP algorithm is computed as follows:

f = 2 * sin(pi * cutoff / samplerate)

The effect of setting f to the input parameter directly, is that it maps 0 to 255 to frequencies of 0 to 1/6 sample rate.

The approximation of sin(x) = x for small x is good enough for this purpose, but the range of the filter is much less than expected. Rather than being able to set the cutoff of this filter between 0 and 8kHz, you can only set it between 0 and 2,7kHz.

@mzero
Copy link
Author

mzero commented Jan 11, 2022

The fix for q is straightforward, requiring a fix-point divide to be defined and used (and be sure to use f, not cutoff in the calculation).

The fix for f is depends on what the API should achieve:

  1. Leave it as it is: The filter has a cutoff range of 0Hz to AUDIO_RATE / 2π - or about 2.6kHz for the 16384kHz sample rate, 5.2kHz with 32768kHz sample rate.

  2. Push the available range by scaling the input by 1.5.

f = 1.5 * cutoff

This still makes use of the sin x = x approximation, though pushes into less accurate territory. This raises the range to 3.9kHz for the lower sampling rate, 7.8kHz for the higher.

  1. The range is still rather low for a low pass filter - there is no way to set the center frequency high enough that it effectively goes away. You can solve this by increasing the range again:
f = 2.0 * cutoff

Now the input parameter, cutoff, in the range 0 ~ 1 does drive the filter from 0Hz to Nyquist. However, it is neither linear nor exponential in frequency, and more than half its range is in in the high end of the spectrum

For all of these, adjustments need to be made for the fixed point arithmetic, since values can be more than 1.

@tomcombriat
Copy link
Collaborator

Hi,

Thanks for your input, and sorry for the late reply!

If I understand well, there are basically two problems with the current implementation:
(1) fb = q + q * (1 - f) should be fb = q + q / (1 - f)
(2) f = cutoff should be (ideally) f=2 * sin(pi * cutoff / samplerate)

I think the LowPassFilter is not the only one with this problem, StateVariable is using the same equations.

Weirdly, I think (1) is an harder fix than (2): division a notoriously costly on small MCU so just replacing the multiplication by a division won't do the trick on the performance side I think.
Actually, (2) might be solved by calculating the good value (with the sin), maybe using one of the wavetable already embedded in Mozzi. This might increase memory consumption a bit…

I will start digging into it and will keep you updated!
Best,

@tomcombriat
Copy link
Collaborator

Okay, I have an idea, will try to implement and test it soon.

Basically, the idea would be to use the limited development (if someone can enlighten me with the correct English translation, I'll be grateful…) of 1/(1-f):

Close to f=0 we have: 1/(1-f) \approx 1+f (+f^2). Hence, fb=q+q/(1-f) \approx q + q (1+f (+ f^2)). At the first order (without the square), this will not be very correct for large f, but will be closer to the ideal value already, without added computational cost.

@sensorium
Copy link
Owner

On behalf of Mozzi users, thanks to @mzero and @tomcombriat for your insights and work on this.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants