C implementation of http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

biquad.h

#ifndef BIQUAD_H
#define BIQUAD_H

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

#ifndef M_PI
#   define M_PI 3.14159265358979323846
#endif

typedef struct Biquad_Params Biquad_Params;
typedef struct Biquad Biquad;

typedef enum {
    BIQUAD_LPF,
    BIQUAD_HPF,
    BIQUAD_BPF_CSG,
    BIQUAD_BPF_CZPG,
    BIQUAD_NOTCH,
    BIQUAD_APF,
    BIQUAD_PEAKING_EQ,
    BIQUAD_LOW_SHELF,
    BIQUAD_HIGH_SHELF
} Biquad_Type;

struct Biquad_Params {
    Biquad_Type type;

    double b0;
    double b1;
    double b2;
    double a1;
    double a2;
};

struct Biquad {
    Biquad_Params* params;

    float x1;
    float x2;
    float y1;
    float y2;
};

static void biquad_setup (Biquad_Params* self, Biquad_Type type, float gain, float freq, float Q, float rate) {
    self->type = type;

    double A;
    double w0;
    double cos_w0;
    double sin_w0;
    double alpha;

    double b0;
    double b1;
    double b2;
    double a0;
    double a1;
    double a2;

    double shelf;

    if (type < BIQUAD_PEAKING_EQ) {
        A = sqrt (pow (10, gain / 20.0));
    } else {
        A =       pow (10, gain / 40.0);
    }
    A = gain;

    w0 = 2 * M_PI * (freq / rate);
    cos_w0 = cos (w0);
    sin_w0 = sin (w0);

    switch (self->type) {
    case BIQUAD_LPF:
    case BIQUAD_HPF:
    case BIQUAD_APF:
        alpha = sin_w0 / (2 * Q);
        break;

    case BIQUAD_BPF_CSG:
    case BIQUAD_BPF_CZPG:
    case BIQUAD_NOTCH:
    case BIQUAD_PEAKING_EQ:
        alpha = sin_w0 * sinh ((log (2.0) / 2.0) * Q * (w0 / sin_w0));
        break;

    case BIQUAD_LOW_SHELF:
    case BIQUAD_HIGH_SHELF:
        Q = Q > 1.0 ? 1.0 : Q;
        alpha = (sin_w0 / 2.0) * sqrt (((A + (1 / A)) * ((1 / Q) - 1)) + 2);
        break;
    }

    shelf = 2 * sqrt (A) * alpha;

    switch (self->type) {
    case BIQUAD_LPF:
        b0 =  (1 - cos_w0) / 2.0;
        b1 =   1 - cos_w0;
        b2 =  (1 - cos_w0) / 2.0;
        a0 =   1 + alpha;
        a1 =  -2 * cos_w0;
        a2 =   1 - alpha;
        break;

    case BIQUAD_HPF:
        b0 =  (1 + cos_w0) / 2.0;
        b1 = -(1 + cos_w0);
        b2 =  (1 + cos_w0) / 2.0;
        a0 =   1 + alpha;
        a1 =  -2 * cos_w0;
        a2 =   1 - alpha;
        break;

    case BIQUAD_BPF_CSG:
        b0 =   Q * alpha;
        b1 =   0.0;
        b2 =  -Q * alpha;
        a0 =   1 + alpha;
        a1 =  -2 * cos_w0;
        a2 =   1 - alpha;
        break;

    case BIQUAD_BPF_CZPG:
        b0 =   alpha;
        b1 =   0;
        b2 =  -alpha;
        a0 =   1 + alpha;
        a1 =  -2 * cos_w0;
        a2 =   1 - alpha;
        break;

    case BIQUAD_NOTCH:
        b0 =   1;
        b1 =  -2 * cos_w0;
        b2 =   1;
        a0 =   1 + alpha;
        a1 =  -2 * cos_w0;
        a2 =   1 - alpha;
        break;

    case BIQUAD_APF:
        b0 =   1 - alpha;
        b1 =  -2 * cos_w0;
        b2 =   1 + alpha;
        a0 =   1 + alpha;
        a1 =  -2 * cos_w0;
        a2 =   1 - alpha;
        break;

    case BIQUAD_PEAKING_EQ:
        b0 =   1 + (alpha * A);
        b1 =  -2 * cos_w0;
        b2 =   1 - (alpha * A);
        a0 =   1 + (alpha / A);
        a1 =  -2 * cos_w0;
        a2 =   1 - (alpha / A);
        break;

    case BIQUAD_LOW_SHELF:
        b0 =      A * ((A + 1) - ((A - 1) * cos_w0) + shelf);
        b1 =  2 * A * ((A - 1) - ((A + 1) * cos_w0));
        b2 =      A * ((A + 1) - ((A - 1) * cos_w0) - shelf);
        a0 =           (A + 1) + ((A - 1) * cos_w0) + shelf;
        a1 =     -2 * ((A - 1) + ((A + 1) * cos_w0));
        a2 =           (A + 1) + ((A - 1) * cos_w0) - shelf;
        break;

    case BIQUAD_HIGH_SHELF:
        b0 =      A * ((A + 1) + ((A - 1) * cos_w0) + shelf);
        b1 = -2 * A * ((A - 1) + ((A + 1) * cos_w0));
        b2 =      A * ((A + 1) + ((A - 1) * cos_w0) - shelf);
        a0 =           (A + 1) - ((A - 1) * cos_w0) + shelf;
        a1 =      2 * ((A - 1) - ((A + 1) * cos_w0));
        a2 =           (A + 1) - ((A - 1) * cos_w0) - shelf;
        break;
    }

    self->b0 = b0 / a0;
    self->b1 = b1 / a0;
    self->b2 = b2 / a0;
    self->a1 = a1 / a0;
    self->a2 = a2 / a0;
}

static void biquad_init (Biquad* self) {
    self->x1 = 0;
    self->x2 = 0;
    self->y1 = 0;
    self->y2 = 0;
}

static void biquad_run (Biquad* self, uint32_t nframes, const float* input, float* output) {
    Biquad_Params* p = self->params;

    for (uint32_t i = 0; i < nframes; ++i) {
        float out = (p->b0 * input[i])
        + (p->b1 * self->x1)
        + (p->b2 * self->x2)
        - (p->a1 * self->y1)
        - (p->a2 * self->y2);

        self->x2 = self->x1;
        self->x1 = input[i];

        self->y2 = self->y1;
        self->y1 = out;

        output[i] = out;
    }
}

#endif

This is beautiful code to reference. Thank you!