C++ Butterworth Nth order filter design

Question

I am looking for a function which calculate a Butterworth Nth filter design coefficients like a 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 haven't any knowledge about DSP.

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

There is also a requirement that the method works under different platforms - Windows, Linux, ...

Answer 1

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

$ aptitude search ~dbutterworth | grep lib

Which gives you answer immediately:

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 library called rtfilter. Description:

rtfilter is a library that provides a set of routines implementing realtime digital filter for multichannel signals (i.e. filtering multiple signals with the same filter parameters). It implements FIR, IIR filters and downsampler for float and double data type (both for real and complex valued signal). Additional functions are also provided to design few usual filters: Butterworth, Chebyshev, windowed sinc, analytical filter...

This library is cross-platform, i.e. works under Linux, MacOS and Windows. From official site:

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

You can install it like this:

$ sudo aptitude install librtfilter-dev librtfilter1

After -dev package is installed, you can even find an example (with Butterworth filter usage) at /usr/share/doc/librtfilter1/examples/butterworth.c. This example (along with corresponding Makefile) also can be found here.

Particularly you are interested in rtf_create_butterworth() function. You can access documentation for this function via command:

$ man rtf_create_butterworth

or you can read it here.

You can specify any filter order passing it as num_pole param to rtf_create_butterworth() function (as far as I remember the number of poles it's the same thing as filter order).

UPDATE

This library doesn't provide external API for coefficients calculation. It only provides actual filtering capabilities, so you can use rtf_filter() to obtain samples (data) after filtering.

But, you can find the code for coefficients calculation in library sources. See compute_cheby_iir() function. This function is static, so it's only can be used inside the library itself. But, you can copy this function code as is to your project and use it. Also, don't let the name of this function confuse you: it is the same algorithm for Butterworth filter and Chebyshev filter coefficients calculation.

Let's say, you have prepared parameters for 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 numerator and denominator for your filter. I have written the 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 wrapper like this:

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 math background for those coefficients calculation, look at DSP Guide, chapter 33.