C ++ Filter Design Butterworth Nth order

Question

I am looking for a function that calculates Butterworth Nth filter design coefficients, such as the Matlab function:

[bl,al]=butter(but_order,Ws); 

and

 [bh,ah]=butter(but_order,2*bandwidth(1)/fs,'high'); 

I found many examples of calculating 2nd order, but not Nth order (for example, I work with order 18 ...). - Unfortunately, I do not know about DSP.

Do you know any library or a way to easily implement this method? When I only know the order, the cutoff frequency and sample rate. I just need to get the vectors B (numerator) and A (denominator).

There is also a requirement that this method work on different platforms - Windows, Linux, ...

Thanks in advance.

Answer 1

It can be easily found (on Debian or Ubuntu):

 $ aptitude search ~dbutterworth | grep lib 

Which gives you the answer right away:

 p librtfilter-dev - realtime digital filtering library (dev) p librtfilter1 - realtime digital filtering library p librtfilter1-dbg - realtime digital filtering library (dbg) 

So you need a library called rtfilter . Description:

rtfilter is a library that provides a set of routines that implement a real-time digital filter for multi-channel signals (i.e. filtering multiple signals with the same filter parameters). It implements FIR, IIR and downsampler filters for float and dual data type (for both real and complex-valued signals). Additional functions are also provided for the development of several conventional filters: Butterworth , Chebyshev, window sinc, analytical filter ...

This library is cross-platform, that is, it works under Linux, MacOS and Windows. From the official website :

rtfilter should compile and run on any POSIX platform (GNU / Linux, MacOSX, BSD ...) and Windows.

You can install it like this:

 $ sudo aptitude install librtfilter-dev librtfilter1 

After installing the -dev package -dev you can even find an example (using the Butterworth filter) in /usr/share/doc/librtfilter1/examples/butterworth.c . This example (along with the corresponding Makefile ) can also be found here .

In particular, you are interested in the rtf_create_butterworth() function. You can access the documentation for this function with the command:

 $ man rtf_create_butterworth 

or you can read it here .

You can specify any filter order passing it as num_pole param for the rtf_create_butterworth() function (as far as I remember the number of poles, this is the same as the order of the filters).

UPDATE

This library does not provide an external API for calculating coefficients. It provides only the actual filtering capabilities, so you can use rtf_filter() to get samples (data) after filtering.

But you can find code for calculating coefficients in library sources. See the compute_cheby_iir () function . This function is static , so it can only be used inside the library itself. But you can copy this function code as is into your project and use it. In addition, you should not confuse the name of this function: this is the same algorithm for calculating the Butterworth filter and Chebyshev filter coefficients.

Let's say you prepared the parameters for the rtf_create_butterworth() function:

 const double cutoff = 8.0; /* cutoff frequency, in Hz */ const double fs = 512.0; /* sampling rate, in Hz */ unsigned int nchann = 1; /* channels number */ int proctype = RTF_FLOAT; /* samples have float type */ double fc = cutoff / fs; /* normalized cut-off frequency, Hz */ unsigned int num_pole = 2; /* filter order */ int highpass = 0; /* lowpass filter */ 

Now you want to calculate the numerator and denominator for your filter. I wrote a wrapper for you:

 struct coeff { double *num; double *den; }; /* TODO: Insert compute_cheby_iir() function here, from library: * https://github.com/nbourdau/rtfilter/blob/master/src/common-filters.c#L250 */ /* Calculate coefficients for Butterworth filter. * coeff: contains calculated coefficients * Returns 0 on success or negative value on failure. */ static int calc_coeff(unsigned int nchann, int proctype, double fc, unsigned int num_pole, int highpass, struct coeff *coeff) { double *num = NULL, *den = NULL; double ripple = 0.0; int res = 0; if (num_pole % 2 != 0) return -1; num = calloc(num_pole+1, sizeof(*num)); if (num == NULL) return -2; den = calloc(num_pole+1, sizeof(*den)); if (den == NULL) { res = -3; goto err1; } /* Prepare the z-transform of the filter */ if (!compute_cheby_iir(num, den, num_pole, highpass, ripple, fc)) { res = -4; goto err2; } coeff->num = num; coeff->den = den; return 0; err2: free(den); err1: free(num); return res; } 

You can use this shell as follows:

 int main(void) { struct coeff coeff; int res; int i; /* Calculate coefficients */ res = calc_coeff(nchann, proctype, fc, num_pole, highpass, &coeff); if (res != 0) { fprintf(stderr, "Error: unable to calculate coefficients: %d\n", res); return EXIT_FAILURE; } /* TODO: Work with calculated coefficients here (coeff.num, coeff.den) */ for (i = 0; i < num_pole+1; ++i) printf("num[%d] = %f\n", i, coeff.num[i]); for (i = 0; i < num_pole+1; ++i) printf("den[%d] = %f\n", i, coeff.den[i]); /* Don't forget to free memory allocated in calc_coeff() */ free(coeff.num); free(coeff.den); return EXIT_SUCCESS; } 

If you are interested in the mathematical background for calculating these coefficients, see the DSP Handbook, chapter 33 .