Designing a FIR lowpass filter in C++ (calculating coefficients?)

[TristanaWild]

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!

[juha_p]

Maybe https://github.com/vinniefalco/DSPFilters can tell you something ?

[Chang-soo Kim]

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.

[aciddose]

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

[Miles1981]

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.

[matt42]

I see only two parameters: cut off + steepness of the low pass.

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.

[Tale]

Maybe this will help:
http://www.labbookpages.co.uk/audio/firWindowing.html

[earlevel]

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.)

[matt42]

Maybe this will help:
http://www.labbookpages.co.uk/audio/firWindowing.html

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]

I was wondering how many terms are needed for reasonable accuracy approximating the Bessel function (for the Kaiser window).

I dunno... Because the factorial gets really large really quick, my JSFX implementation does just 20 rounds, which seems to work pretty well:

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;
);

And here is a more precise C version, which continues until it runs out of precision:

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;
}

I also found this paper which gives an alternative method of approximating the function.

Interesting!

[matt42]

I dunno... Because the factorial gets really large really quick, my JSFX implementation does just 20 rounds, which seems to work pretty well

Thanks for the insight.

I'll have to compare the two approaches at some point.

[aciddose]

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.

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)));
}

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.

[matt42]

... deleted post. Never mind.

[aciddose]

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.

[matt42]

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).