Butterworth

A collection of 2nd-order butterworth filter collections. These will be explained and probably coded up a little bit better later. For now, these are imports from Soundpipe.

The digital Butterworth filter is initially designed in the S-plane like you would an analogue filter, which results in a special transfer function. This transfer function then gets digitized to the Z-plane using a process known as the bilinear transform, or BLT. By the time this filter gets implemented in C code, this transformation has already been applied. What you're left with is a biquad with a bunch of constants from the BLT.

Going through the motions of converting the Butterworth filter from the S-plane to the Z-plane is a little bit beyond the scope of this sndkit. However, Julius Smith has a page that derives a Second-order butterworth lowpass filter.

#ifndef SK_BUTTERWORTH_H
#define SK_BUTLERWORTH_H

#ifndef SKFLT
#define SKFLT float
#endif

typedef struct sk_butterworth sk_butterworth;

void sk_butterworth_init(sk_butterworth *bw, int sr);
void sk_butterworth_freq(sk_butterworth *bw, SKFLT freq);
void sk_butterworth_bandwidth(sk_butterworth *but, SKFLT bw);

SKFLT sk_butlp_tick(sk_butterworth *bw, SKFLT in);
SKFLT sk_buthp_tick(sk_butterworth *bw, SKFLT in);
SKFLT sk_butbp_tick(sk_butterworth *bw, SKFLT in);

#ifdef SK_BUTTERWORTH_PRIV
struct sk_butterworth {
    SKFLT freq, lfreq;
    SKFLT bw, lbw; /* for bandpass filter only */
    SKFLT a[7];
    SKFLT pidsr, tpidsr;
};
#endif

#endif

This code is based off the soundpipe implementation, slightly reworked to make it so that malloc is not used.

Tangled Code

#include <math.h>
#include <stdint.h>
#include <stdlib.h>

#define SK_BUTTERWORTH_PRIV
#include "butterworth.h"

#define ROOT2 1.4142135623730950488

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
<<common>>
<<filters>>

Common

The butterworth filters share some common code, defined below.

The filter function performs the difference equation to compute one sample of audio. in is the input signal, and a refers to the entire filter state: indices 0-4 hold filter coefficients, and 5-6 hold previously computed values. This has been reworked from the cannonical difference equation to be more optimized.

sk_butterworth_init initializes the filter data for all butterworth filters.

sk_butterworth_freq sets the frequency.

static SKFLT filter(SKFLT in, SKFLT *a)
{
    SKFLT t, y;
    /* a5 = t(n - 1); a6 = t(n - 2) */
    t = in - a[3]*a[5] - a[4]*a[6];
    y = t*a[0] + a[1]*a[5] + a[2]*a[6];
    a[6] = a[5];
    a[5] = t;
    return y;
}

void sk_butterworth_init(sk_butterworth *bw, int sr)
{
    int i;
    sk_butterworth_freq(bw, 1000.0);
    sk_butterworth_bandwidth(bw, 1000.0);
    bw->lfreq = -1;
    bw->lbw = -1;
    bw->pidsr = M_PI / (SKFLT)sr;
    for (i = 0; i < 7; i++) bw->a[i] = 0.0;
}

void sk_butterworth_freq(sk_butterworth *bw, SKFLT freq)
{
    bw->freq = freq;
}

void sk_butterworth_bandwidth(sk_butterworth *but, SKFLT bw)
{
    but->bw = bw;
}

Filters

Lowpass

SKFLT sk_butlp_tick(sk_butterworth *bw, SKFLT in)
{
    if (bw->freq != bw->lfreq) {
        SKFLT *a, c;
        a = bw->a;
        bw->lfreq = bw->freq;
        /* derive C constant used in BLT */
        c = 1.0 / tan((SKFLT)(bw->pidsr * bw->lfreq));

        /* perform BLT, store components */
        a[0] = 1.0 / (1.0 + c*ROOT2 + c*c);
        a[1] = 2*a[0];
        a[2] = a[0];
        a[3] = 2.0 * (1.0 - c*c) * a[0];
        a[4] = (1.0 - c*ROOT2 + c*c) * a[0];
    }

    return filter(in, bw->a);
}

Highpass

SKFLT sk_buthp_tick(sk_butterworth *bw, SKFLT in)
{
    if (bw->freq != bw->lfreq) {
        SKFLT *a, c;
        a = bw->a;
        bw->lfreq = bw->freq;
        /* derive C constant used in BLT */
        c = tan((SKFLT)(bw->pidsr * bw->freq));

        /* perform BLT, store components */
        a[0] = 1.0 / (1.0 + c*ROOT2 + c*c);
        a[1] = -2*a[0];
        a[2] = a[0];
        a[3] = 2.0 * (c*c - 1.0) * a[0];
        a[4] = (1.0 - c*ROOT2 + c*c) * a[0];
    }

    return filter(in, bw->a);
}

Bandpass

SKFLT sk_butbp_tick(sk_butterworth *bw, SKFLT in)
{
    if (bw->bw != bw->lbw || bw->freq != bw->lfreq) {
        SKFLT *a, c, d;
        a = bw->a;
        bw->lfreq = bw->freq;
        bw->lbw = bw->bw;

        /* Perform BLT and store components */
        c = 1.0 / tan((SKFLT)(bw->pidsr * bw->bw));
        d = 2.0 * cos((SKFLT)(2.0*bw->pidsr * bw->freq));
        a[0] = 1.0 / (1.0 + c);
        a[1] = 0.0;
        a[2] = -a[0];
        a[3] = - c * d * a[0];
        a[4] = (c - 1.0) * a[0];
    }

    return filter(in, bw->a);
}