Digital Signal Processing (DSP)

Audio programming

Audio programming can be very fun, but it also come with a particular set of challenges that are not often found in other programming fields. Audio processing is typically a real-time problem (unless you are performing offline rendering), in particular we can classify it as a firm real-time system, meaning that if we fail to meet the processing “deadline”, we won’t normally cause a full system failure, but audible glitches can be very very noticeable by humans and can be annoying or outright insufferable.

In real-time audio programming, we will usually have a callback function that will be periodically be called by the audio thread when new samples are needed. This function will be tasked to fill up N samples in a buffer (or buffers in stereo). For example, if the audio sampling rate is 44.1KHz and the buffer size is of 256 samples, the callback can be called each (256 / 44100 * 1000) = 5.80ms. Furthermore, the longer the audio buffers are, the longer the latency/responsiveness of our audio. Anything longer than 5-15ms will be quite noticeable when trying to play in real time, but it’s not so problematic otherwise (for example, a 100ms audio latency in a DAW may be annoying, but it’s still usable, but good luck trying to play an instrument with that kind of delay).

Some important things to consider when working on the audio callback function are:

Math

Complex numbers

A lot of signal processing math makes use of complex numbers in the rectangular form (a + j * b) or polar form (A = sqrt(a^2 + b^2); phi = atan2(b/a)). The later is commonly used as the vector representation of the polar form makes more intuitive sense when multiplying or adding vectors.

The Euler’s formula gives us a mapping between the rectangular and polar forms:

exp(j * phi) = cos(phi) + j * sin(phi)
cos(phi) = (exp(j * phi) + exp(-j * phi)) / 2
sin(phi) = (exp(j * phi) - exp(-j * phi)) / (2 * j)

Sinusoidal waves

The most fundamental of audio waves is the sine wave. In a discrete context, a real sinewave takes the form of:

y[n] = A * cos(w * n * Ts + phi) =
     = A * cos(2 * pi * f * n * Ts + phi) =
     = A * cos(2 * pi * n * f / fs + phi)

where:  A := amplitude
        w := angular frequency (radians/second)
        f = w / (2 * pi) := frequency in Hertz (cycles/second)
        phi := initial phase in radians
        n := time index
        Ts := sampling period (= 1 / fs)
        fs := sampling frequency

Here is how to generate a real sine wave in Python:

A = .9
f = 440
fs = 44100
phi = np.pi / 2
n = np.arange(-0.002, 0.002, 1.0 / fs)
y = A * np.cos(2 * np.pi * f * n + phi)

In audio DSP we often work with complex sinewaves of the form:

s[n] = exp(j * 2 * pi * k * n / N) = cos(2 * pi * k * n / N) + j * sin(2 * pi * k * n / N)
where: N := Number of samples
       k := Number of wave cycles
N = 500
k = 3
n = np.arange(-N/2, N/2)
s = np.exp(2j * np.pi * k * n / N)

Convolution

A very important operation in DSP is the convolution, which in a discrate domain can be represented as:

c[n] = convolve(a, b, n)

def convolve(a, b, n):
    c = 0
    for m in range(0, N - 1):
        c += a[m] * b[n - m]
    return c

Discrete Fourier Transform (DFT)

The DFT transforms a discrete time domain signal into its spectral representation.

X[k] = sum_n(x[n] * exp(-j2 * pi * k * n / N)), k = 0,..., N-1, n = 0, ... N-1

The counterpart inverse DFT (IDFT) takes a spectral signal and recovers the corresponding time domain data:

x[n] = 1/N * sum_k(X[k] * exp(j * 2 * pi * k * n / N)), k = 0,..., N-1, n = 0, ... N-1

In Python, these can be implemented as follows:

def dft(x):
    N = len(x)
    X = np.zeros(N, dtype = np.complex_)
    for k in range(0, N):
        n = np.arange(N)
        s = np.exp(-2j * np.pi * k * n / N)
        X[k] = np.sum(x * s)
    return X

def idft(X):
    N = len(X)
    x = np.zeros(N, dtype = np.complex_)
    for n in range(0, N):
        k = np.arange(N)
        s = np.exp(2j * np.pi * k * n / N)
        x[n] = np.sum(X * s) / N
    return x

Each bin n in the DFT output corresponds to bin = n * fs / N. Where N is the number of samples used to calculate the DFT. Note that higher FFT windows result in narrower spectral peaks for detected harmonics and thus higher spectral resolution.

While the spectrum is a complex signal, we normally visualize it using its magnitude and phase graphs. The phase graph can look quite messy, but using phase unwrapping (adding 2 * pi at discontinuity points) we have a more visually clear representation of the phase.

If we use zero-padding (adding zeroes at the end of a signal) on a time domain signal before DFT, the computed spectrum will be smoother. It is akin to an interpolation in the spectral domain.

Properties

Fast Fourier Transform (FFT)

The DFT can be quite demanding to calculate, but a more efficient algorithm called the FFT is used in practice. The FFT takes advantage of the symetry properties to dramatically reduce the computation time, going from an O(N^2) (DFT) to O(N log(N)) (FFT).

To use the FFT, our input signal must be a power of 2. To calculate the FFT of an arbitrary length signal, we use a combination of zero-padding and zero-phase windowing. The way it works, is that we split the original signal in the middle, the rightmost part will be located at the beginning of the fft buffer, followed by the zero padding and lastly by the leftmost part of the original signal.

Resources

Short-time Fourier Transform (STFT)

If we work with real time audio, we don’t always have the entirety of the audio signal when we start processing, but rather, we operate on chunks/frames of audio. Furthermore we can define the STFT as a way of working with these chunks of audio. In short, it involves the multiplication of the chunk with a window function and then performing a normal fourier transform of the signal for that chunk:

X[k] = sum(w[n] * x[n+iH] * exp(-j * 2 * pi * k * n / N))
    where:
        n := time index (-N/2 <= n <= N/2 - 1)
        w := analysis window
        i := frame number
        H := hop-size

The effect of windowing on the spectra of a sinusoid is the same as shifting the spectra of the window to both the negative and positive frequencies of the sinusoid.

There are different types of window, with different spectral properties. Usually we evaluate windows by measuring the width of the main lobe and the amplitude of the side lobes.

Some of the most common windows used in audio processing are:

In addition to the window type, we also need to consider the effect of window size, odd vs even windows, FFT size and hop size. Smaller window sizes have poor spectral resolution compared with larger ones. Odd size windows have their phases centered around zero, which is preferable. FFT size can be chosen independently from the window size (this involves zero-padding the signal). Large hopping size can produce modulation artifacts of the spectra, which can be audible after the IDFT (for example if M = 201 and H = 100, this isn’t audible when M = 201 and H = 50).

For analysis purposes we should be aware of the time-frequency tradeoff: While larger window sizes give us better spectral resolution, we lose time domain information, and viceversa.

The inverse STFT for a single frame is the same as the regular IDFT, though when multiple frames are involved, they should be shifted and added accordingly. For a full sound it bouls down to: y[n] = x[n] * sum(w[n - i * H]). Depending on the window type and if the amplitude modulation is to be avoided, different overlap sizes can be tried. For example, for a Blackman window, a 4 * N overlap avoids any distortion. A similar result is obtained with a 2 * N overlap when using a non-symmetric Hann window.

For more information about FFT windows, check out this article

The sinusoidal model

For a given sound, we can detect the main harmonics by using a sinusoidal model, where the harmonic peaks produced by the STFT can be detected on the spectral domain by looking at the peaks. In essence we model a sound as a sum of sinusoids:

y[n] = sum(As[n] * cos(2 * pi * fs[n] * n))

To obtain enough frequency resolution for this purpose we need to consider the effect of the main lobe width of the window function.

If B = Bw * fs / M and delta = abs(f[k+1] - f[k]):

M >= Bw * fs / delta

where: B := Bandwidth of the main-lobe of a window
       Bw := Number of samples of the main-lobe
       fs := Sampling frequency (Hz)
       M  := Window size
       f[k] and f[k+1] := Frequency of sinusoids in Hz.
       delta := Frequencies we want to resolve as distinct

In many cases we want to detect the harmonics of a signal, so we consider delta = f0 (the fundamental frequency / initial harmonic of a sound) and thus M >= Bw * fs / f0. A good window size is twice the minimum number, in other words: M = 2 * Bw * fs / f0

The peaks in a spectrum can be detected as a local maxima by looking at the derivative sign or for each point checking if the neighbours at either side are smaller. This method is subject to noise, luckily an easier method to obtain better resolution in a signal spectrum is to use zero-padding, for example by making the FFT size N much larger than the window size M. Other smoothing methods could also be used depending on your needs. This method may still not give us enough resolution unless the N increse is very large however, which would be computationally very expensive. An alternative method is to try to fit a parabola x[n] = a * (n - p) ^ 2 + b where p is the center of the parabola, a is the concavity measure and b is the offset. We can combine the zero-padding with the parabolic interpolation to obtain a more precise peak center.

So far the method described would work for a single frame, but in order to use a real signal, we need to make sure the peaks are stable across multiple frames. We can use a spectrogram for this purpose and define stability in terms of frequency and amplitude as well as the phase derivative in time/frequency.

The sinusoidal model can be used for additive synthesis in addition to be an analysis tool. Adding the values of multiple sinusoids can be quite an expensive operation when working on the time domain, but this can be done much faster by working with the spectral model and then using the IDFT to obtain the final sound. Using one of the windows we previously discuss, we can generate arbitrary sinusoids. Even better, we only need to focus on the main lobe of the window for this purpose if we use a window with a high enough SNR, like the Blackman-Harris.

The harmonic model

Following from the sinusoidal model, where a sound is represented by a sum of sinusoids. The harmonic model restricts which harmonics/sinusoids can be part of the model by making them a multiple of the fundamental frequency:

y[n] = sum(A[n] * cos(2 * pi * r * f0[n] * n))

where: A[n] := instantaneous amplitude
       f0[n] := fundamental frequency (Hz)
       r := harmonic number

Y[k] = sum(A * W[k - r * f0'])

where: A := instantaneous amplitude
       W := spectrum of the analysis window
       f0' := normalized fundamental frequency
       r := harmonic number

The main challenge of this model is the identification of the fundamental frequence (f0), specially when considering not only monophonic signals, but also polyphonic or enharmonic ones. Polyphonic detection is only really feasible in the spectral domain, but we can perform f0 detection for monophonic signals on the time domain.

One of the tools we have at our disposal is the autocorrelation function in the time domain:

r[l] = sum(x[n] * x[n + l]), where l := lag = 0, 1, ..., N - 1

We can calculate the autocorrelation for different lag values to study the peak distribution, in some signals is easier than others.

Another useful function for this purpose is the YIN algorithm, often use in speech monophonic signals. It tries to find the minimum value of the square difference of samples:

d[l] = sum((x[n] - x[n + l])^2), where l := lag = 0, 1, ..., N - 1

In the frequency domain we can detect harmonic peaks as previously described, and then heuristically select the common divisor(s) of the harmonic series that best explain the spectral peaks. The two-way mismatch (TWM) algorithm by Maher and Beuchamp is one method of performing pattern matching between the predicted and measured peaks. This algorithm will work with polyphonic signals as well as monophonic, but clearly the former are more challenging to deal with.

Similarly Salamon and Gómez (2012) propose and algorithm to find the prominent pitch in polyphonics signals.

These are quite involved algorithms that go a bit too far from me to summarize right now, but the original research articles should serve as a good baseline for their understanding. Frequency detection is still an ongoing research field, particularly with complex signals.

One useful thing to note is that if we know the fundamental frequency, we can calculate the maximum number of harmonics that can appear in a signal with:

n_harm = fs / (2 * f0)

When talking about harmonics and the harmonic series, also known as the overtone series it’s worth noting that when working with the western chromatic scale of 12 semitones per scale, the difference between one pitch and the next is of:

pitch[p] = pitch[p - 1] * 2 ^ (1 / 12)

Between two subsequent pitches we divide the frequency range in 100 “cents”. The difference in “cents” between two frequencies f0 and f1 (where f0 < f1) is:

n_cents = 1200 * log2(f1 / f0)

Since the pitch distribution is logarithmic, the multiple harmonics of a fundamental frequency for a given pitch have an musical intervalic interpretation. The fundamental frequency is considered the root of the interval, the second harmonic is an octave, the 3rd is a 5th. These intervals are not perfect when compared to equal temperament tuning. In this table you can find the initial 16 harmonics and the intervalic representation for each of them.

Harmonic Interval from root Detune (Cents) Cents from root
1 I 0 ct
2 VIII 1200 ct
3 V +2 ct 1902 ct
4 VIII 2400 ct
5 III -14 ct 2786 ct
6 V +2 ct 3102 ct
7 VIIb -31 ct 3369 ct
8 VIII 3600 ct
9 II +4 ct 3804 ct
10 III -14 ct 3986 ct
11 Vb -49 ct 4151 ct
12 V +2 ct 4302 ct
13 VIb +41 ct 4441 ct
14 VIIb -31 ct 4569 ct
15 VII -12 ct 4688 ct
16 VIII 4800 ct

The stocastic model

Sinusoidal and harmonic models tell a story about the periodic sounds of a signal. We can analyze a signal using stochastic models with tools like autocorrelation, power spectral density, mean, variance and probability distributions. Complex sounds, the attack of instruments and noise can be studied with stochastic models, and we can use these in addition to the harmonic or sinusoidal models. Think of sounds like the ocean waves breaking on the shore, the wind hitting a microphone, etc.

There are many stochastic models available, one of which consist onthe convolution of white noise with the filtered approximation of our signal:

y[n] = sum(u[n] * h[n-k])
    where: u[n] := white noise
           h[n] := impulse response of filter approximating input signal x[n]

The phase of white noise is essentially composed of random numbers within the phase range, for this reason we can see the model in the spectral domain as the magnitude of our filtered input sound and random numbers as the phase.

One approach to obtain the filtered magnitude spectra is to use the LPC approximation, where we minimize the error function such that:

h[n] = sum(a[k] * x[n - k])
error = infsum((x[n] - sum(a[k] * x[n - k]) ^ 2)) = infsum((x[n] - h[n]) ^ 2)
    where: a[k] := filter coefficients

This method is very useful to obtain the formant (main resonances) of a vocal sound. Read more about formant synthesis here.

A simpler approach is the use of an envelope approximation by using a low pass filter and zero padding the dft before the idft transform.

a[n] = IDFT(LP(DFT(x[n])))
b[n] = IDFT(DFT(ZP(a[n])))
    where: LP := low pass filter
           ZP := zero-padding

We can also use residuals for stochastic models. The residual refers to the the stochastic model obtained from the substraction of the modeled sinusoids (from the sinusoidal or harmonic model) to the original signal:

y[n] = sum(A[n] * cos(2 * pi * f[n] * n)) + xr[n] = ys[n] + xr[n]
    where: ys[n] := sinusoidal model composed of R sinusoids
           xr[n] := x[n] - ys[n]
           x[n] := our original signal

Y[k] = sum(A * W[k - f]) + Xr[k] = Ys[k] + Xr[k]

Residuals are difficult to model, however, since they are formed of large amounts of data without many knobs to adjust, which is why other stochastic models may be more useful for our purposes if that knob tweaking is the goal. What we can do, however, is to take the extracted residuals, and then use it for modeling the stochastic components as previously described.

Transformations

We can apply transformations to our signal in several ways, but it is very convenient to do these on the frequency domain or to modify one of our models. For example, filtering via convolution is an operation that can be expensive to compute in the time domain, but it’s just a multiplication on the frequency domain. Better yet, if our magnitude spectra is in logarithmic scale (dB), instead of multiplications we use additions (log(a * b) = log(a) + log(b)), which are generally much cheaper to compute. In the STFT context, a filter doesn’t change with time and it will be applied to each frame as a multiplication in the spectral magnitude and addition of the phases, though we generally only care about the magnitude of a filter, and if the filter is zero-phase windowed it will have no effect on the end result:

Y[l][k] = |X[l][k]| * |H[k]| * exp(j * phase_x[l][k]  + j * phase_h)

Another operation we may want to perform on the spectral domain is “morphing”. It’s similar to filtering, except the filter also changes in time. For example if we filtered our original x[n] signal with a z[n] signal that is transformed to the spectral domain and whose magnitude is smoothed out. For example we can make an orchestral sample x[n] sound as if being played by a vocal sound z[n], think a Vocoder effect.

If we have a sinusoidal model, we only want to modify the amplitude values of the sinusoids and later regenerate the phase values based on the modified model, since the phase is very sensitive to modifications. With this, we can resythesize samples, setting onsets in different locations. Likewise, we can modulate the frequencies to change in a variety of ways. The similar concept can be applied to harmonic models and compound harmonic/stochastic models. The possibilities are endless with regards to transformation and morphing operations. Some examples include frequency transposition (fh[l] = s * f[l] * f[l]), frequency shifting (fh[l] = s * f[l] + f[l]), frequency stretching (fh[l] = fh[l] / h * h^(s * f[l])).

Audio features

Sound analysis has a large number of tools in addition to the models previously described. Analyzing an audio file can result in a variety of metrics/features that can be used for studying sound characteristics or as a jumping point for more creative transformations.

Feature extraction can be done in the time and spectral domain, and we can also differentiate algorithms that work on a single frame or multiple frames. Some examples:

Tools

Sampling

Sampling is easy enough when the sampling rate of the sample matches the playing audio rate (e.g. 48KHz samples on a 48KHz audio stream), wrapping around as needed in case of looping samples:

typedef struct Sample {
    float *data;
    int len;
    int pos;
} Sample;

int
audio_callback(float *output, int len) {
    for (int x = 0; x < len; x++) {
        if (sample->data == NULL) {
            continue;
        }

        // Looping.
        if (sample->pos >= sample->len) {
            sample->pos = 0;
        }

        output[x] = sample->data[sample->pos++];
    }
}

When “resampling”, for increasing/decreasing the pitch or changing the sample duration, we will run into a number of issues. Resampling means we will be playing the sample at a different speed that the one it was originally recorded at. We generally want a fine control over this speed, meaning we want a fractional number for our positional increment (floating or fixed point numbers). For now I’ll stick to floating point numbers, but fixed precision integers have a number of advantages we will go over at a later time.

We already know that an increment of 1 will play the sound at the original rate, but what should the sample increment be to play the note at a desired pitch? Well, a rate of 2 corresponds to the original pitch shifted one octave up, and conversely a rate of 0.5 will be one octave down instead.

For all the semitones in between, they are separated from the previous semitone by 2 ^ (1 / 12), so:

...
B  = C  / pow(2, (1/12))
C# = C  * pow(2, (1/12))
D  = C# * pow(2, (1/12))
...

Furthermore, we can pre-calculate the frequencies and fractional increment amounts for each note, for example, using a Python script to cover C0-C9:

C4 = 261.625580
interval = 2 ** (1 / 12)
for i, octave in enumerate([1./16., 1./8., 1./4., 1./2., 1., 2., 4., 8., 16., 32.]):
    base = C4 * octave
    freqs = [base]
    increments = [octave]
    for j in range(1, 12):
        freqs += [freqs[j - 1] * interval]
        increments += [increments[j - 1] * interval]

    print(list(zip(freqs, increments)))

If we run this script and look at the C4-B4 octave we get the following note/increment pairs:

Note Frequency Increment
C4 261.62558 1.0
C#4 277.1826465503454 1.0594630943592953
D4 293.6647844169278 1.122462048309373
D#4 311.12700120271364 1.1892071150027212
E4 329.6275754329552 1.2599210498948734
F4 349.2282510543508 1.3348398541700346
F#4 369.9944434997273 1.4142135623730954
G4 391.99545800596655 1.498307076876682
G#4 415.3047209137905 1.5874010519682
A4 440.00002472134804 1.6817928305074297
A#4 466.16378770944584 1.7817974362806792
B4 493.8833290048991 1.8877486253633877

We can store this table as an array and use a lookup to get the increment for any given note. Our sample playing code can now be something like this, making also sure that when looping, we keep the remainder of the fractional part to avoid glitches:

typedef struct Sample {
    float *data;
    float len;
    float pos;
} Sample;

int
audio_callback(float *output, int len) {
    for (int x = 0; x < len; x++) {
        if (sample->data == NULL) {
            continue;
        }

        // Looping.
        while (sample->pos >= sample->len) {
            sample->pos -= sample->len;
        }

        int pos = sample->pos;
        output[x] = sample->data[pos];
        sample->pos += sample->inc;
    }
}

Note that this doesn’t quite do the trick just yet. If you try using this function, you will notice some high pitch harmonics creeping in, even if you use a pure sine sample. One of the things causing this effect is the fact that we are using a sample and hold method for selecting samples. By truncating the fractional part of the position, we create a “staircase” effect, adding a bunch of harmonics to our sample. Instead, you probably want to use some form of sample interpolation, the simplest of which is linear interpolation. This method is already a huge improvement over no interpolation, and can be quite fast to calculate:

float
sample_lerp(float *data, size_t len, float pos) {
    float x = pos;
    int x0 = (int)x;
    int x1 = (x0 + 1);
    float y0 = data[x0];
    float y1 = data[x1 % len];
    x = x - x0;
    float y = y0 + x * (y1 - y0);
    return y;
}

Another popular interpolation method is the cubic interpolation, more specifically the Hermite interpolation. This method requires 4 points, where the point we want to interpolate is located between the two middle ones. I was scratching my head trying to understand the derivation for this polynomial, I’m not gonna go into a lot of detail here, but essentially, on a 3rd order polynomial you set some constrains so that your derivative behave in a smoother manner. More details can be found in this StackExchange answer and The Continuity of Splines video by Freya Holmer. The important thing about this method is that it presents a smoother path across samples, and with a smoother wave, you have less added artificial harmonics.

float
sample_hermite(float *data, size_t len, float pos) {
    float x = pos;
    int x0 = x - 1;
    int x1 = x;
    int x2 = x + 1;
    int x3 = x + 2;
    float y0 = data[x0 % len];
    float y1 = data[x1];
    float y2 = data[x2 % len];
    float y3 = data[x3 % len];
    x = x - x1;
    float c0 = y1;
    float c1 = 0.5f * (y2 - y0);
    float c2 = y0 - 2.5f * y1 + 2.0f * y2 - 0.5f * y3;
    float c3 = 1.5f * (y1 - y2) + 0.5f * (y3 - y0);
    return ((c3 * x + c2) * x + c1) * x + c0;
}

In practice, linear interpolation is good enough for many purposes, but it’s good to have a better method available in case we have the processing power to spare.

Resources

Filtering

EMA

Formula:

EMA LP: y[n] = alpha * x[n] + (1 - alpha) * y[n - 1]
        where:     x := input
                   y := output
               alpha := filter amount [0.0 .. 1.0]

EMA HP: y[n] = 1/2 * (2 - beta) * (x[n] - x[n - 1]) + (1 - beta) * y[n - 1]
        where:     x := input
                   y := output
                beta := filter amount [0.0 .. 1.0]

Polyphase filtering

A common task in DSP is to perform resampling via upsampling and downsampling of a signal. Without getting too much into the math (if you are reading this, you probably already know what the Nyquist limit is right?), we need to apply a low pass filter in addition to the resampling to avoid aliasing.

For upsampling, we can zero stuff a signal with zeroes. For example, for an upsamping factor L = 2, we put one zero in between each pair of samples, for L = 3 we put 3 and so on. After this zero stuffing, we apply a low pass filter to interpolate between the original samples.

         ┌─────┐    ┌─────┐
x[n] ───>│  L  │───>│ LPF │───> y[n], fs = fs * L
         └─────┘    └─────┘

In the case of downsampling, we do it the other way around, with a low pass filter first to avoid frequencies that would otherwise alias, and then we apply decimation, meaning we keep only each sample multiple of M. So for M = 2, we would remove every other sample.

         ┌─────┐    ┌─────┐         
x[n] ───>│ LPF │───>│  M  │───> y[n], fs = fs / M
         └─────┘    └─────┘

You can already see a problem if you pay attention to this, we end up having to filter a bunch of samples that we end up throwing away. What a waste!

Here is where Polyphase filters come into play. The math it’s a bit too involved to cover here, but I’ll link some resources below with more in depth explanations. I’m going to focus on the practical part here. For a filter H[n] we can switch our signal path so that the downsampler becomes:

                          ┌──────┐    
              ┌──────────>│  H1  │───┐
              │           └──────┘   │ 
              │                      │
              │           ┌──────┐   │ 
              │ ┌────────>│  H2  │───┤
           ┌─────┐        └──────┘   │
 x[n]  ─── │ MUX │                   ├── + ──> y[n]
           └─────┘        ┌──────┐   │
              │ └────────>│  H3  │───┤
              │           └──────┘   │
              │              ...     │
              │           ┌──────┐   │
              └──────────>│  Hn  │───┘
                          └──────┘

Where MUX is a sample multiplexer that will push sample n = 0 to H1, sample n = 1 to H2, etc. And each of the subfilters can be easily obtained by separating the coefficients in a similar way, each N intervals. This sounds a bit abstract so let’s look at an example.

Say we have an FIR filter of 4 coefficients in the form:

b = c0 + c1*z-1 + c2*z-2 + c3*z-3

If we want to do decimation by a factor of 2, the subfilters H0 and H1 are split as such:

H0 = c0 + c2*z-1
H1 = c1 + c3*z-1

So far so good right? For me, I’ve found at this point it’s a bit difficult to understand how the signal flow actually work, should the multiplexing work in sync with the inner filters or not? How does the summation work? In summary, we process all the samples, but each of the input samples get sent into a different filter, so x[0] goes to H0, x[1] to H1, x[2] to H0, etc. The trick is that the output is half the rate of the input, so for y[0] we need to send the sample x[0] to the first filter and gather the output of all of the filters for adding them together. Still confused? So was I! So I decided to experiment a bit in python. Here is the code for an M=2 decimator using a stateful FIRFilter class that keep track of the internal delayed samples.

class FIRFilter:
    def __init__(self, b):
        self.b = b
        self.xz = np.zeros(len(b))
        self.y = 0.0

    def filter_sample(self, xn):
        self.xz[0] = xn
        acc = 0.0
        for n in range(len(self.b)):
            acc += self.xz[n] * self.b[n]

        # Shift delay line.
        self.xz = np.roll(self.xz, 1)
        self.y = acc
        return acc

    def filter(self, x):
        y = np.zeros(len(x))
        for i, x in enumerate(x):
            y[i] = self.filter_sample(x)
        return y

def filter_polyphase(b, x):
    N = len(x)
    y = np.zeros(N//M)
    h0 = FIRFilter(b[0::M])
    h1 = FIRFilter(b[1::M])
    for i, n in enumerate(range(0, N, M)):
        acc = h0.filter_sample(x[n]) + h1.y
        h1.filter_sample(x[n + 1])
        y[i] = acc
    return y

We can exted this to an M bank polyphase decimator easily:

def filter_polyphase(b, x, M):
    N = len(x)
    y = np.zeros(N//M)
    bank = []
    for i in range(M):
        bank.append(FIRFilter(b[i::M]))
    for i, n in enumerate(range(0, N, M)):
        acc = 0.0
        for j in range(0, M):
            acc += bank[j].filter_sample(x[n + M - 1 - j])
        y[i] = acc
    return y

Let’s prove that this works! Assuming we are already oversampled by a factor of 2, we would like to downsample to 48kHz. We can generate a simple filter:

width = 1000
fs = 48000 * 2
fc = 24000 - width / 2
numtaps = 256
b = signal.remez(numtaps, [0, fc, fc + width, fs/2], [1, 0], fs=fs)

We will test a simple implementation of the combination of filtering and decimation with an impulse response.

def decimate(x, M):
    y = np.zeros(len(x)//M + 1)
    for i, n in enumerate(range(0, len(x), M)):
        y[i] = x[n]
    return y

x = np.zeros(1024)
x[0] = 1.0
fir = FIRFilter(b)
y0 = fir.filter(x)
y1 = decimate(y0, 2)

Giving us the following time and frequency response:

Now let’s see how it compares to our polyphase filter:

y2 = filter_polyphase(b, x, 2)

We can see that the output of the polyphase filter fully overlaps with the filtered and decimated signal.

Amazing! And now each sample just have to do half of the filtering operations. The magic of this technique is that if we were to downsample by a factor of 4, for the same filter coefficients each sample only has to do 1/4 of of the original multiplies.

Resources

Antialiasing

So we know that we may want to use a polyphase filter for an antialiasing filter with downsampling, and this may work well for M = 2. In the context of audio applications, we oversample and downsample our signal to avoid aliasing when harmonics are generated, for example by some saturation or waveshaping function. Once the aliasing harmonics are baked into a signal, they are impossible to remove, however if we operate at a higher sample rate, the harmonics will take longer and have more time to dissipate before bouncing back from the Nyquist limit. The more we oversample, the less reflections we will have, though our processing becomes much more expensive as we would have to process many more samples.

Let’s show an example of what this looks like by creating different chirp signals (sine sweeps) and distorting them with a tanh waveshaper. Note that we apply the distortion to the “oversampled” signal:

gain = 1.0
x0 = np.tanh(signal.chirp(np.arange(48000) / 48000, f0=100, f1=22000, t1=1) * gain)
x2 = np.tanh(signal.chirp(np.arange(48000 * 2) / 48000 / 2, f0=100, f1=22000 * 2, t1=2) * gain)
x4 = np.tanh(signal.chirp(np.arange(48000 * 4) / 48000 / 4, f0=100, f1=22000 * 4, t1=4) * gain)
x8 = np.tanh(signal.chirp(np.arange(48000 * 8) / 48000 / 8, f0=100, f1=22000 * 8, t1=8) * gain)
x16 = np.tanh(signal.chirp(np.arange(48000 * 16) / 48000 / 16, f0=100, f1=22000 * 16, t1=16) * gain)
x32 = np.tanh(signal.chirp(np.arange(48000 * 32) / 48000 / 32, f0=100, f1=22000 * 32, t1=32) * gain)
x64 = np.tanh(signal.chirp(np.arange(48000 * 64) / 48000 / 64, f0=100, f1=22000 * 64, t1=64) * gain)
x128 = np.tanh(signal.chirp(np.arange(48000 * 128) / 48000 / 128, f0=100, f1=22000 * 128, t1=128) * gain)

Before we start, why don’t we take a look at an optimized FIR filter for this purpose? We would like a pretty steep filter with a flat passband and a low stopband attenuation, say of -100dB. A good candidate for this design is to use a blackman-harris window FIR filter. We can check a few other windows like the blackman or hamming for comparison.

As you can see, for the same number of taps, the attenuation is much better on the stopband on the blackmanharris. For this window we can study the frequency response for different number of taps.

Until n = 16 we get sufficient attenuation, even if it’s progressively less steep. If we were to downsample in multiple steps, we could use a 128 tap filter for the last pass, a 32 for the second to last pass and 16 for everything before that. Here is how you can get the coefficients for those filters for a frequency cut slightly below fs/2.

b2 = signal.firwin(128, 0.455, window='blackmanharris')
b4 = signal.firwin(32, 0.455, window='blackmanharris')
b8 = signal.firwin(16, 0.455, window='blackmanharris')

However, with our polyphase filters, we can actually just downsample directly with a factor M. Bear in mind though, that to maintain the same filter performance at lower fractions of the frequency range, we need to progressively increase the number of taps. The good news is that the performance of the filter will be the same no matter what, since the work is divided, though the increse in number of samples also will affect the amount of work needed.

ntap2 = 128
ntap4 = ntap2 * 2
ntap8 = ntap4 * 2
ntap16 = ntap8 * 2
ntap32 = ntap16 * 2
ntap64 = ntap32 * 2
ntap128 = ntap64 * 2
b2 = signal.firwin(ntap2, 0.90/2, window=window)
b4 = signal.firwin(ntap4, 0.90/4, window=window)
b8 = signal.firwin(ntap8, 0.90/8, window=window)
b16 = signal.firwin(ntap16, 0.90/16, window=window)
b32 = signal.firwin(ntap32, 0.90/32, window=window)
b64 = signal.firwin(ntap64, 0.90/64, window=window)
b128 = signal.firwin(ntap128, 0.90/128, window=window)

Time of the truth, we will up to x128 oversampling and see the effect in our spectrograms.

y2 = filter_polyphase(b2, x2, 2)
y4 = filter_polyphase(b4, x4, 4)
y8 = filter_polyphase(b8, x8, 8)
y16 = filter_polyphase(b16, x16, 16)
y32 = filter_polyphase(b32, x32, 32)
y64 = filter_polyphase(b64, x64, 64)
y128 = filter_polyphase(b128, x128, 128)

Starting at x8 oversampling, our aliasing is mostly gone, hurray!

But wait, what happens if you increase the gain in the signal, hence adding more distortion? Let’s start by bumping it up to 2.0

Uh oh, we see now some aliasing starting to creep in at x8. May still be fine for most applications, but here is what happens if we increase the gain to 10.0

A lot more audible harmonics there, but we may just increase the oversampling factor if we have the cpu to spare. What if we go all the way to a gain of 100 times?

Wowza, that’s a lot of aliasing, even at x64 (and arguably x128). Of course that’s quite a lot of extra gain and distortion being applied, but it shows that oversampling is not always the solution for all our aliasing woes. Things like hard clipping or square waves will inevitably add a lot of aliasing to our signal if we are not careful, though there are other techniques that can be used to produce alias free oscillators.

Comb filters

Synthesis

Karplus–Strong

Wavetable

Effects

Reverb

Other Resources

Audio APIs