Short One Pole Low-pass

I often use this filter for parameter smoothing.

You can get an intuitive feeling for the filter coefficient c by observing that when it is near 1.0, the out slowly moves to the new value (low-pass). That is very little of the new input is passed, and mostly the previous state is preserved. If c is near 0.0, the out follows the new value closely (pass-thru)

out = new * (1 - c) + out * c;

c can be calculated for a given frequency by:

c = exp(-2 * PI * frequency/samplingRate);

You can obtain a high-pass by subtracting the low-pass output from the input.

hpout = new - out;

Coefficient calculation and sound processing

Most filters have one section of code for setting the frequency, resonance and gain (the coefficients calculation) and another section for filtering the audio (using those coefficients and saving the internal state in one or more delay lines). It’s typical for the coefficient calculation to be more CPU intensive than the filtering of sound. This is because coefficient calculation involves trigonometric functions.

Low and high pass filters derived from an all-pass filter

All filters have both a defined amplitude and phase response over frequency. An allpass filter is a filter that has a flat amplitude response. The particular allpass filter we will be using is one which causes the phase to shift by π or 180 degrees at the nyquist frequency, and has 0 degrees phase shift at 0 Hz. It uses a single coefficient c, which will vary the transition from 0 degrees to 180 degrees, either causing it to happen more quickly or more slowly. If it happens more quickly, lower frequencies will have a phase shift of near 180.

out = c * (in - out) + in1;
in1 = in;

Note that the previous value of in is copied to in1 to be used in the next sample calculation, and out is used in the next sample calculation as well. These two variables need to be preserved as state variables.

The coefficient can be calculated for a given frequency as follows:

omega = PI * frequency/samplerate;
tf = tan(omega);
c = (tf -1.0)/(tf + 1.0);

Now as the phase is shifted by 0 degrees at 0 Hz and 180 degrees at the nyquist frequency (and lower frequencies, dependent on c), adding the in signal to the allpass out will cause high frequencies to cancel for a 0 amplitude, and low frequencies to reinforce for an amplitude of 2. With one additional line of code, we have a simple low-pass filter.

lpout = (in + out) * 0.5;

And the high-pass is just as simple.

hpout = (in - out) * 0.5;

Notch and band-pass filters derived from a second order all-pass filter

By putting two all-pass filters in series, we can double the phase shift to 2π or 360 degrees. With this second order all-pass filter we can create a notch or bandpass filter by adding the input to the output of the all-pass. This is because the 180 degree phase shift will now appear between 0 Hz and Nyquist. Bandwidth of this filter is controlled by how quickly the phase goes from 0 to 360 phase shift.

The all-pass is calculated as follows:

out = c * (out2 - in) + d * (1.0 - c) * (in1 - out1) + in2;
out2 = out1;
out1 = out;
in2 = in1; 
in1 = in;

It shouldn’t be any surprise that we now have four state variables (out1, out2, in1, in2) – twice as many as in the first order filter. We now have two coefficients c and d. These are calculated so that we have independent control of frequency and bandwidth.

tf = tan(PI * bandwidth/sampleRate);
c = (tf -1.0)/(tf + 1.0);  
d = -cos(2 * PI * frequency/sampleRate);

And like the first-order low-pass and high-pass filters, we can create notch and band-pass filters by adding and subtracting the output of the all-pass to the input.

notch = (in + out) * 0.5;
bandpass = (in + out) * 0.5;

Externals for Pure Data

Just to show this code working, here are two filters programmed as Pure Data externals (fof~ and sof~). Each can be switched between addition and subtraction for low-pass and high-pass or notch and band-pass.

Source for the fof~ external – download

The sample processing (fof_perform) function from the fof~ external:

static t_int *fof_perform(t_int *w)
{
    t_float tf;
    t_fof *x = (t_fof *)(w[1]);
    t_float *in = (t_float *)(w[2]);
    t_float *out = (t_float *)(w[3]);
    int n = (int)(w[4]);

    // calculate coefficient only if frequency changes
    if(x->frequency != x->frequencyold)
    {
       if(x->frequency > (x->samplerate * 0.48f))
            x->frequencyold = x->frequency = (x->samplerate * 0.48f);
        else if(x->frequency < 0.0001f)
            x->frequencyold = x->frequency = 0.0001f;
        else
            x->frequencyold = x->frequency;
        tf = tanf(x->PI * x->frequency/x->samplerate);
        x->c = (tf -1.0f)/(tf + 1.0f);
    }

    while (n--)
    {
        float input = *(in++);
        x->out = x->c * (input - x->out) + x->in1;
        x->in1 = input;

        switch(x->type)
        {
            case 0:
                *out++ = (input + x->out) * 0.5f; // lp
                break;
            case 1:
                *out++ = (input - x->out) * 0.5f; // hp
                break;
            case 2:
                *out++ = x->out; // ap
                break;
            case 3:
            default:
                *out++ = input; // bypass
                break;
        }
    }
    return (w+5);
}

Source for the sof~ external – download

The sample processing (sof_perform) and bandwidth setting (sof_bandwidth) function from the sof~ external:

void sof_bandwidth(t_sof *x, t_floatarg g)
{
    t_float tf;
    x->bandwidth = g;
    tf = tanf(x->PI * x->bandwidth/x->samplerate);
    x->c = (tf - 1.0f)/(tf + 1.0f);
}

 

static t_int *sof_perform(t_int *w)
{
    t_float tf;
    t_sof *x = (t_sof *)(w[1]);
    t_float *in = (t_float *)(w[2]);
    t_float *out = (t_float *)(w[3]);
    int n = (int)(w[4]);

    // calculate coefficient only if frequency changes
    if(x->frequency != x->frequencyold)
    {
       if(x->frequency > (x->samplerate * 0.48f))
            x->frequencyold = x->frequency = (x->samplerate * 0.48f);
        else if(x->frequency < 0.0001f)
            x->frequencyold = x->frequency = 0.0001f;
        else
            x->frequencyold = x->frequency;
        x->d = -1.0f * cosf(2.0f * x->PI * x->frequency/x->samplerate);    }
    while (n--)
    {
       float input = *(in++);
        x->out = x->c * (x->out2 - input)
            + x->d * (1.0f - x->c) * (x->in1 - x->out1)
            + x->in2;
        x->out2 = x->out1;
        x->out1 = x->out;
        x->in2 = x->in1;
        x->in1 = input;
        switch(x->type)
        {
            case 0:
                *out++ = (input + x->out) * 0.5f;
                break;
            case 1:
                *out++ = (input - x->out) * 0.5f;
                break;
            case 2:
                *out++ = x->out;
                break;
            case 3:
            default:
                *out++ = input;
                break;
        }
    }
    return (w+5);
}