DSP Audio Programming Series: Part 2

It’s been quite some time since part 1 of this series (all the way back in June of last year) where we looked at and implemented  a basic delay effect. In this part, I will be introducing filters, another staple of the audio effect and processing toolkit. Audio filtering is a very large and complex topic, spanning many different types and designs, so we’ll be looking at a fairly basic resonating low-pass and high-pass filter in this part to keep things simple and manageable.

Like part 1, the code will be using Portaudio for cross-platform audio input/output. The code will also be made available on Github for both Mac and Windows to do with as you please. We will be sticking to a command-line based application as well, which means it is less practical to adjust parameters in real-time unfortunately. This wasn’t a big problem when looking at a simple delay effect, but having real-time control over the parameters of a filter is far more convenient. Fortunately it’s fairly simple to set up a basic GUI application that connects controls to the parameters of the filter (the README in the Github projects contains some tips of setting up a GUI in both OS X and Windows).

Ok, let’s get to it!

Digital filters can be categorized into two basic types: finite impulse response (FIR), and infinite impulse response (IIR). The former works by mixing the samples of the original signal with delayed copies of the input, which is essentially a feed-forward network, whereas the latter combines the original signal with delayed copies of both the input and the output, making it a feedback network. IIR filters have the advantage of requiring fewer delay elements than FIR filters to achieve comparable results, and are usually less computationally expensive, but they introduce a non-linear phase shift in the signal in contrast with FIR filters that can be designed with linear phase.

Here we will be looking at IIR filters because they are more commonly used for real-time audio processing, though FIR filters have many important applications in digital audio as well (such as anti-aliasing filters and noise reduction).

The most basic implementation of an IIR filter is the biquad direct form I, with the equation:

y(n) = a₀x(n) + a₁x(n-1) + a₂x(n-2) – b₁y(n-1) – b₂y(n-2)

where x is the input and y is the output. Thus, x(n) is the current input sample, x(n-1) is the 1-sample input delay, and x(n-2) is the 2-sample input delay. It follows then that y(n-1) is the previous output and y(n-2) is the 2-sample output delay. The a and b terms of the equation are the filter coefficients, which determine the response and characteristics of the filter (which includes its type, such as low-pass, high-pass, etc.).

This topology is illustrated below.

Direct form I filter topology.

The corresponding code in CSFilter that calculates the outgoing sample according to the equation above is straightforward:

inline float processSample (const float inSample, const int channel)
{
    double outSample;
    int idx0 = mNumChannels * channel;
    int idx1 = idx0 + 1;
 
    outSample = (mCoeffs.a0 * inSample) + (mCoeffs.a1 * mZx[idx0]) + (mCoeffs.a2 * mZx[idx1]) - (mCoeffs.b1 * mZy[idx0]) - (mCoeffs.b2 * mZy[idx1]);
 
    mZx[idx1] = mZx[idx0];
    mZx[idx0] = inSample;
    mZy[idx1] = mZy[idx0];
    mZy[idx0] = outSample;
 
    return (float)(outSample * mGain);
}

The last block of code in the method above updates the filter state delay variables before returning the filtered sample. This is similar to dealing with buffers in the delay effect from part 1, except here we have two delay buffers (one for input and one for output) of only two samples each.

Now that we’ve covered the filter equation, let’s look at calculating the filter coefficients.

w = 2πfc/fs // angular frequency
d = 1 / Q // temporary variable
beta = ( (1 - (d/2)sin(w)) / (1 + (d/2)sin(w)) ) / 2
gamma = (0.5 + beta)cos(w)

Low-pass coefficients:

a0 = (0.5 + beta - gamma) / 2
a1 = 0.5 + beta - gamma
a2 = a0
b1 = -2 * gamma
b2 = 2 * beta

High-pass coefficients:

a0 = (0.5 + beta + gamma) / 2
a1 = -(0.5 + beta + gamma)
a2 = a0
b1 = -2 * gamma
b2 = 2 * beta

(Source: “Designing Audio Effect Plug-Ins in C++”, by Will Pirkle)

Where above, fc = cutoff frequency, fs = sampling rate, and Q = quality factor. With Q set to a value of 1 / sqrt(2), there will be no resonant peaking in the filter, while greater values of Q introduce a resonant peak at the cutoff frequency.

It is very helpful to plot the frequency response of a filter in order to visualize how it will affect the audio signal, so below we see plots of both the low-pass and high-pass filters at a cutoff frequency of 2000Hz and varying Q (the frequency axis is logarithmic).

LPF with cutoff frequency at 2000Hz and Q = 0.707.

LPF with cutoff frequency at 2000Hz and Q = 2.3.

HPF with cutoff frequency at 2000Hz and Q = 0.707.

HPF with cutoff frequency at 2000Hz and Q = 2.3.

The cutoff frequency for low-pass and high-pass filters is defined as the frequency at which attenuation is -3dB, and we can see with a magnified view of the plots that this is the case with this particular filter.

Attenuation is -3dB at the 2000Hz cutoff frequency.

Attenuation is -3dB at the 2000Hz cutoff frequency.

We can now proceed to the implementation by filling in the callback function we supply to the Portaudio engine.  This is where processing happens; Portaudio provides us with a buffer of audio data (input) and a place to store our output that is sent to the audio hardware.  In order to apply the filter effect to the incoming audio, we pass in the CSFilter instance to the userData parameter, and this gives us access to the class methods that implement the filter.

int audioCallback (const void* input, void* output, 
        unsigned long samples, const PaStreamCallbackTimeInfo* timeInfo, 
        PaStreamCallbackFlags statusFlags, void* userData)
{
    const float *in = (const float*)input;
    float *out = (float*)output;
    CSFilter *filter = (CSFilter*)userData;
    int channels = filter->getNumChannels();
 
    for (int i = 0; i < samples; ++i) {
        for (int chan = 0; chan < channels; ++chan) {
            *out++ = filter->processSample(*in++, chan);
        }
    }
 
    return paContinue;
}

The processSample method as we saw above, executes the filter’s difference equation, so the only main thing left is for the filter class to calculate the coefficients.

void CSFilter::calculateCoeffs ()
{
    double theta = 2. * M_PI * mNormalizedFreq; // normalized frequency has been precalculated as fc/fs
    double d = 0.5 * (1. / mQuality) * sin(theta);
    double beta = 0.5 * ( (1. - d) / (1. + d) );
    double gamma = (0.5 + beta) * cos(theta);
 
    if (mFilterType == FILTER_TYPE_LOW_PASS) {
        mCoeffs.a0 = 0.5 * (0.5 + beta - gamma);
        mCoeffs.a1 = 0.5 + beta - gamma;
    }
    else {
        mCoeffs.a0 = 0.5 * (0.5 + beta + gamma);
        mCoeffs.a1 = -(0.5 + beta + gamma);
    }
 
    mCoeffs.a2 = mCoeffs.a0;
    mCoeffs.b1 = -2. * gamma;
    mCoeffs.b2 = 2. * beta;
}

To finish off, let’s briefly examine how the filter equation and coefficients work to produce the signal the it does. As we know from the Fourier theorem, any complex signal is made up of an infinite (at least in the continuous, analog world) number of frequency components, which can be decomposed into its individual frequencies; i.e. sine waves. So let’s see what happens when we send a single sine wave through the filter with its samples passing through the direct form I network illustrated above; first at 200Hz, then at 2000Hz. The filter has the same parameters as in the above plots — 2000Hz cutoff, and a Q of 0.707.

Top: original 200Hz sine wave. Bottom: filtered 200 Hz sine wave.

We can see from comparing the two that aside from a modest phase shift, the sine wave remained unaffected by the filter. When we pass a 2000Hz sine wave through, however, it’s a different story.

Top: original 2000Hz sine wave. Bottom: filtered 2000Hz sine wave.

The phase shift is obviously more pronounced and the amplitude of the wave has been reduced from 1.0 to 0.7, which is equal to -3dB, exactly what we would expect at the 2000Hz cutoff frequency. Extrapolating from this, it’s clear that higher frequency sine waves will be attenuated more and more, resulting in the low-pass filter effect. The high-pass filter works the same way of course, just reversed.

There’s really no limit to the use of filters in audio processing, from the practical to the artful. I hope this has been an interesting introduction to digital filters. Here are the links to the Xcode and Visual Studio projects on Github if you’re interested in exploring more on what we’ve built here.

CSFilter for Mac on Github: csfilter-mac
CSFilter for Windows on Github: csfilter-win