Designing a FIR lowpass filter in C++ (calculating coefficients?)
-
- KVRer
- Topic Starter
- 5 posts since 8 Jul, 2016
Hey all,
I have a couple use cases in a plug-in I'm writing where having a FIR lowpass filter would either be handy or is necessary.
I've found a couple examples of FIR implementations in C, but they're all relying on coefficients that are pre-calculated. This won't work for me, as I need a filter that can perform on bufferSize # of samples, with a particular cutoffFrequency that's in relation to a sampleRate.
How can I go about calculating FIR coefficients with those three parameters in mind? I've been stumped by this, so your help is much appreciated!
One tentative solution I've come up with is to populate a buffer with a rectangle signal, then to perform an FFT on it to get a Sinc (hopefully straightforward, as I'm working with JUCE which has an FFT class), and then to convolve that Sinc with my signal-to-be-filtered. The thing is, I'm uncertain how to create a rectangle that will communicate to the proper bandwidth for my filter given the above three parameter specifications. I'm also uncertain of if this is the most efficient method of calculating the coefficients or not.
Thanks for any help!
I have a couple use cases in a plug-in I'm writing where having a FIR lowpass filter would either be handy or is necessary.
I've found a couple examples of FIR implementations in C, but they're all relying on coefficients that are pre-calculated. This won't work for me, as I need a filter that can perform on bufferSize # of samples, with a particular cutoffFrequency that's in relation to a sampleRate.
How can I go about calculating FIR coefficients with those three parameters in mind? I've been stumped by this, so your help is much appreciated!
One tentative solution I've come up with is to populate a buffer with a rectangle signal, then to perform an FFT on it to get a Sinc (hopefully straightforward, as I'm working with JUCE which has an FFT class), and then to convolve that Sinc with my signal-to-be-filtered. The thing is, I'm uncertain how to create a rectangle that will communicate to the proper bandwidth for my filter given the above three parameter specifications. I'm also uncertain of if this is the most efficient method of calculating the coefficients or not.
Thanks for any help!
-
- KVRian
- 833 posts since 21 Feb, 2006 from FI
Maybe https://github.com/vinniefalco/DSPFilters can tell you something ?
-
- KVRer
- 15 posts since 14 Jul, 2016
Low-pass FIR filter:
taps[x] = sin(2 * PI * f * x) / (PI * x) * 2 * f;
taps[0] = 2 * f;
where
x = tap index (for example, from -50 to +50)
f = normalized cutoff frequency = cutoff frequency / sample rate frequency
After calculating the coefficients, you will need to apply window function.
taps[x] = sin(2 * PI * f * x) / (PI * x) * 2 * f;
taps[0] = 2 * f;
where
x = tap index (for example, from -50 to +50)
f = normalized cutoff frequency = cutoff frequency / sample rate frequency
After calculating the coefficients, you will need to apply window function.
- KVRAF
- 12555 posts since 7 Dec, 2004
sinc(x) = sin(pi x) / (pi x)
https://en.wikipedia.org/wiki/Sinc_function
Note that you must window the function and so begin the numerous trade-offs made when designing a FIR filter.
Also note that while sinc is the transform of a rectangular function, do you really want a rectangular function? In many cases another "non-ideal" filter may produce preferable results especially where excessive ringing is an issue.
Other foundational functions include the Gaussian function; it is its own transform. (Similar to how sin is its own integral shifted in phase.)
https://en.wikipedia.org/wiki/Gaussian_function
https://en.wikipedia.org/wiki/Gaussian_filter
https://en.wikipedia.org/wiki/Sinc_function
Note that you must window the function and so begin the numerous trade-offs made when designing a FIR filter.
Also note that while sinc is the transform of a rectangular function, do you really want a rectangular function? In many cases another "non-ideal" filter may produce preferable results especially where excessive ringing is an issue.
Other foundational functions include the Gaussian function; it is its own transform. (Similar to how sin is its own integral shifted in phase.)
https://en.wikipedia.org/wiki/Gaussian_function
https://en.wikipedia.org/wiki/Gaussian_filter
Free plug-ins for Windows, MacOS and Linux. Xhip Synthesizer v8.0 and Xhip Effects Bundle v6.7.
The coder's credo: We believe our work is neither clever nor difficult; it is done because we thought it would be easy.
Work less; get more done.
The coder's credo: We believe our work is neither clever nor difficult; it is done because we thought it would be easy.
Work less; get more done.
-
- KVRian
- 1379 posts since 26 Apr, 2004 from UK
I see only two parameters: cut off + steepness of the low pass.
Do you really need a FIR or can you do with an IIR? They are definitely easier to define. Using FFT is one way of defining you filter, but don;t forget that a good filter will require many coefficients.
Do you really need a FIR or can you do with an IIR? They are definitely easier to define. Using FFT is one way of defining you filter, but don;t forget that a good filter will require many coefficients.
-
- KVRian
- 1273 posts since 9 Jan, 2006
Well I think they are considering sample rate a parameter also which could be a valid view point, but it's trivial to arrive at frequency/sample rate and then as you say it's just 2 parameters to design the filter.Miles1981 wrote:I see only two parameters: cut off + steepness of the low pass.
- KVRian
- 519 posts since 12 Apr, 2010 from The Netherlands
-
- KVRian
- 653 posts since 4 Apr, 2010
A couple of tips:
The windowed-sinc method is easy to understand, and to compute on the fly (for instance, if the user can set arbitrary cutoff points, or if you need a cutoff point that's a particular frequency, say 10 kHz, and need to calculate its coefficients based on the current sample rate). For musical purposes, save yourself some head-scratching and go directly to the Kaiser window. It has good properties for audio, including the ease of which you can trade off steepness and stop-band rejection.
If you want a shorter filter for a given performance, look at equiripple designs. Not so easy to design on the fly, but you can find a coefficient calculator on the web. You can trade off, for instance, a 0.1 dB ripple in the passband for a shorter filter. (Windowed sinc filters have a flat passband.)
The windowed-sinc method is easy to understand, and to compute on the fly (for instance, if the user can set arbitrary cutoff points, or if you need a cutoff point that's a particular frequency, say 10 kHz, and need to calculate its coefficients based on the current sample rate). For musical purposes, save yourself some head-scratching and go directly to the Kaiser window. It has good properties for audio, including the ease of which you can trade off steepness and stop-band rejection.
If you want a shorter filter for a given performance, look at equiripple designs. Not so easy to design on the fly, but you can find a coefficient calculator on the web. You can trade off, for instance, a 0.1 dB ripple in the passband for a shorter filter. (Windowed sinc filters have a flat passband.)
My audio DSP blog: earlevel.com
-
- KVRian
- 1273 posts since 9 Jan, 2006
I was wondering how many terms are needed for reasonable accuracy approximating the Bessel function (for the Kaiser window). I also found this paper which gives an alternative method of approximating the function.Tale wrote:Maybe this will help:
http://www.labbookpages.co.uk/audio/firWindowing.html
- KVRian
- 519 posts since 12 Apr, 2010 from The Netherlands
I dunno... Because the factorial gets really large really quick, my JSFX implementation does just 20 rounds, which seems to work pretty well:matt42 wrote:I was wondering how many terms are needed for reasonable accuracy approximating the Bessel function (for the Kaiser window).
Code: Select all
function i0(x)
local(x2, fact, y, i)
(
x2 = x = 0.25*x*x;
fact = y = i = 1;
loop(20,
y += x / (fact*fact);
x *= x2;
fact *= (i += 1);
);
y;
);
Code: Select all
#include <float.h>
double i0(double x)
{
const double x2 = x = 0.25*x*x;
double fact = 1, y = 1;
for (int i = 1;;)
{
const double tmp = x / (fact*fact);
if (tmp < DBL_MIN) break;
y += tmp;
x *= x2;
fact *= (double)++i;
}
return y;
}
Interesting!matt42 wrote:I also found this paper which gives an alternative method of approximating the function.
-
- KVRian
- 1273 posts since 9 Jan, 2006
Thanks for the insight.Tale wrote:I dunno... Because the factorial gets really large really quick, my JSFX implementation does just 20 rounds, which seems to work pretty well
I'll have to compare the two approaches at some point.
- KVRAF
- 12555 posts since 7 Dec, 2004
Before going to such effort you'd better check whether this has any effect on the application you're using it for.
It likely makes a difference with large (>4096) filters but won't make any significant difference above about -160 dB.
range = 1e-10 results in max 27 iterations. The difference between 1e-10 and 1e-20 is 33 (+6) iterations and a tiny fraction of a dB in the result.
I'd rather use Dolph-Chebyshev or sinc windows (produces rectangular lobes) most of the time though both for processing and analysis if I'm going to use a simple window rather than a complex combination.
It likely makes a difference with large (>4096) filters but won't make any significant difference above about -160 dB.
Code: Select all
alpha = 1 to 8
template <int alpha>
kaiser()
{
beta = 1.0 + double(alpha - 1) * 2.5;
}
template <typename T>
static T bessel(T x)
{
const T range = T(1e-10);
double sum = 1.0;
double pw = 1.0;
const double halfx = double(x) / 2.0;
double n = 1.0;
do {
double temp = halfx / n;
pw *= temp * temp;
sum += pw;
n += 1.0;
} while (pw >= range * sum);
return T(sum);
}
template <typename T>
T operator()(T x)
{
double tmp = fabs(x * T(2.0) - T(1.0));
return T((bessel(beta * sqrt(1.0 - tmp*tmp)) / bessel(beta)));
}
I'd rather use Dolph-Chebyshev or sinc windows (produces rectangular lobes) most of the time though both for processing and analysis if I'm going to use a simple window rather than a complex combination.
Free plug-ins for Windows, MacOS and Linux. Xhip Synthesizer v8.0 and Xhip Effects Bundle v6.7.
The coder's credo: We believe our work is neither clever nor difficult; it is done because we thought it would be easy.
Work less; get more done.
The coder's credo: We believe our work is neither clever nor difficult; it is done because we thought it would be easy.
Work less; get more done.
- KVRAF
- 12555 posts since 7 Dec, 2004
Still it would make sense before going to such effort to determine what effect it has on the application.
It is extremely unlikely Kaiser is the best possible window for a particular application, and if so (that's quite a leap!) before looking at the cost of the iterative Bessel function it should be considered as a percentage of the cost of computing the window and what effect adjusting the number of iterations has on the window in the application.
I'm not aware of any reason you wouldn't cache a window and in so doing you might take advantage of the mirror symmetry of the window immediately cutting the cost by 1/2.
I posted the complete code needed to generate the window so people can assess this before looking at such optimizations. I'd say it's quite likely it isn't the best window for an application and if so any other Bessel function implementation would be useless premature optimization. In most cases the most important component of optimization is algorithm selection.
It is extremely unlikely Kaiser is the best possible window for a particular application, and if so (that's quite a leap!) before looking at the cost of the iterative Bessel function it should be considered as a percentage of the cost of computing the window and what effect adjusting the number of iterations has on the window in the application.
I'm not aware of any reason you wouldn't cache a window and in so doing you might take advantage of the mirror symmetry of the window immediately cutting the cost by 1/2.
I posted the complete code needed to generate the window so people can assess this before looking at such optimizations. I'd say it's quite likely it isn't the best window for an application and if so any other Bessel function implementation would be useless premature optimization. In most cases the most important component of optimization is algorithm selection.
Free plug-ins for Windows, MacOS and Linux. Xhip Synthesizer v8.0 and Xhip Effects Bundle v6.7.
The coder's credo: We believe our work is neither clever nor difficult; it is done because we thought it would be easy.
Work less; get more done.
The coder's credo: We believe our work is neither clever nor difficult; it is done because we thought it would be easy.
Work less; get more done.
-
- KVRian
- 1273 posts since 9 Jan, 2006
Hi aciddose,
You must have read my previous comment before I deleted it. I did so, because I felt unsure exactly what you were referring to. Anyway it's clearer now and I see your point.
I haven't had the need to play around with FIR filters, so my knowledge is just what I've read up on. I guess the advantage with Kaiser would be deriving M (filter length) from a set of filter parameter specifications. Where as with Delph-Chebyshev we would need to experiment to find the ideal M. If we need to work on the fly this is a big advantage with Kaiser. If our needs are more static, then we can take the time to test different sizes of M, window types and parameters to get the best performance for a specific application.
This is my hunch anyway. I'm interested in this topic, currently, as I'll soon need to design a high quality resampling algorithm (offline kinda, but still needs to be fast).
You must have read my previous comment before I deleted it. I did so, because I felt unsure exactly what you were referring to. Anyway it's clearer now and I see your point.
I haven't had the need to play around with FIR filters, so my knowledge is just what I've read up on. I guess the advantage with Kaiser would be deriving M (filter length) from a set of filter parameter specifications. Where as with Delph-Chebyshev we would need to experiment to find the ideal M. If we need to work on the fly this is a big advantage with Kaiser. If our needs are more static, then we can take the time to test different sizes of M, window types and parameters to get the best performance for a specific application.
This is my hunch anyway. I'm interested in this topic, currently, as I'll soon need to design a high quality resampling algorithm (offline kinda, but still needs to be fast).