In [1]:
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
from scipy.special import spence
from IPython.core.display import SVG, Image, display


# Practical Considerations for Antiderivative Anti-Aliasing¶

In audio signal processing, it is often useful to use nonlinear functions for things like saturation, wavefolding, and more. While nonlinear functions can create interesting sounds, they often create difficulty in the world of signal processing, particularly by invalidating many of the handy mathematical theorems that hold up in the world of linear signal processing. One particularly nasty issue that arises when using nonlinear functions on digital audio is aliasing.

In this article, we'll give a brief discussion of what aliasing is, and how it is typically dealt with. Then, we'll introduce a relatively new anti-aliasing technique called Antiderivative Anti-Aliasing (ADAA), and discuss some practical considerations for implementing various nonlinear systems with ADAA.

## Aliasing¶

First, what is aliasing? This question brings us to the heart of digital signal processing. A digital system with sample rate $f_s$ can only faithfully reproduce signals up to the frequency $f_s/2$, often referred to as the Nyquist frequency, named after electrical engineer Harry Nyquist of Bell Labs. Any signal above the Nyquist frequency will be reflected back over it. For example, if digital system with sample rate $f_s = 48$ kHz attempts to reproduce a signal at 50 kHz, the signal will be reflected back to $48 - (50 - 48) = 46$ kHz. This reflected signal is known as an "aliasing artifact" and is typically considered undesirable, particularly since the artifact is not harmonically related to the original signal. For a more in-depth explanation of aliasing, see here.

Typically, digital systems use filters to supress any signal above the Nyquist frequency, so that aliasing artifacts don't occur. However, when using a nonlinear function, the digital system can create signal above the Nyquist frequency, thus creating aliasing artifacts. As an example let's look at a specific type of distortion known as "hard-clipping" distortion. A hard-clipper is a waveshaping distortion with a waveshaping curve that looks as follows:

In [2]:
adsp.plot_static_curve(adsp.hard_clipper, gain=5)
plt.grid()
plt.title('Hard Clipper Response')

Out[2]:
Text(0.5,1,&#39;Hard Clipper Response&#39;)

If I now create a 1.2 kHz sine wave, and process it through a hard-clipping distortion with no aliasing, I would see the following frequency response:

In [3]:
def plot_fft(x, fs, sm=1.0/24.0):
fft = 20 * np.log10(np.abs(np.fft.rfft(x) + 1.0e-9))
freqs = np.fft.rfftfreq(len(x), 1.0 / fs)
return freqs, fft

def process_nonlin(fc, FS, nonlin, gain=10):
N = 200000
sin = np.sin(2 * np.pi * fc / FS * np.arange(N))
y = nonlin(gain * sin)
freqs, fft = plot_fft(y, FS)
return freqs, fft

FC = 1244.5
FS = 1920000

freqs_analog, fft_analog = process_nonlin(FC, FS, adsp.hard_clipper)
peak_idxs = signal.find_peaks(fft_analog, 65)[0]
plt.plot(freqs_analog, fft_analog)
plt.scatter(freqs_analog[peak_idxs], fft_analog[peak_idxs], c='r', marker='x')
plt.xlim(0, 20000)
plt.ylim(10)
plt.title('Hard Clipping Distortion w/ No Aliasing')
plt.ylabel('Magnitude [dB]')
plt.xlabel('Frequency [Hz]')
plt.grid()


The spike farthest to the left represents the original sine wave, while all the other spikes represent the higher order harmonics generated by the hard-clipper. Note that all of the harmonics are evenly spaced.

Now, if we were to sample the sine wave at $f_s = 48$ kHz before processing it with the hard clipper, any generated harmonics above $f_s / 2 = 24$ kHz would be reflected back to produce aliasing artifacts. Let's see what that looks like:

In [4]:
FC = 1244.5
FS = 48000

freqs_alias, fft_alias = process_nonlin(FC, FS, adsp.hard_clipper)
plt.plot(freqs_alias, fft_alias)
plt.scatter(freqs_analog[peak_idxs], fft_analog[peak_idxs], c='r', marker='x')
plt.xlim(0, 20000)
plt.ylim(10)
plt.title('Hard Clipping Distortion w/ Aliasing')
plt.ylabel('Magnitude [dB]')
plt.xlabel('Frequency [Hz]')
plt.grid()


In this case, we see that all the blue spikes not marked with the red x's are aliasing artifacts. This level of aliasing will almost certainly be audible for people listening to signals that have been processed with this effect, and probably won't sound very pleasant.

## Anti-Aliasing with Oversampling¶

Now that we have some understanding of what aliasing is, and how nonlinear signal processing can create aliasing, let's look at the most common method for supressing aliasing artifacts: oversampling.

The idea behind oversampling is very simple: if we run our nonlinear process at a higher sample rate, then any harmonics produced above the Nyquist frequency can be filtered out before the audio is downsampled back to its original sample rate. Let's see how well 4x oversampling works for out hard-clipping distortion:

In [5]:
freqs, fft = process_nonlin(FC, FS*4, adsp.hard_clipper)
plt.plot(freqs, fft)
plt.scatter(freqs_analog[peak_idxs], fft_analog[peak_idxs], c='r', marker='x')
plt.xlim(0, 20000)
plt.ylim(10)
plt.title('Hard Clipping Distortion w/ 4x Oversampling')
plt.ylabel('Magnitude [dB]')
plt.xlabel('Frequency [Hz]')
plt.grid()


The result is pretty good! While aliasing artifacts are still visible they are at a much lower magnitude than the original signal. More oversampling can be used to supress aliasing artifacts even further.

However, using oversampling has some drawbacks, notably with performance. Specifically, the processing time for an effect will be multiplied by the oversampling factor. In other words, if a process is done with 4x oversampling, the oversampled process will be $1/4$th as efficient as the original process. With that in mind, let us look for other anti-aliasing methods, that hopefully introduce less performance overhead.

## Antiderivative Antialiasing¶

Antiderivative anti-aliasing (abbreviated as ADAA), was first introduced in a 2016 DAFx paper by Julian Parker, Vadim Zavalishin, and Efflam Le Bivic from Native Instruments. The technical background for ADAA is then further developed in an IEEE paper by Stefan Bilbao, FabiÃ¡n Esqueda, Parker, and Vesa VÃ¤limÃ¤ki. I won't go to much into the mathematical details of ADAA, but the basic idea is that instead of applying a nonlinear function to a signal, we instead apply the anti-derivative of that function, and then use discrete-time differentiation resulting in a signal with supressed aliasing artifacts.

ADAA can be implemented as follows: Say we have a nonlinear function, $y[n] = f(x[n])$, with an antiderivative $F_1(x) = \int_0^x f(t) dt$. A first-order ADAA version of the function, can be written as follows: $$y[n] = \frac{F_1(x[n]) - F_1(x[n-1])}{x[n] - x[n-1]}$$

Unfortunately, when $|x[n] - x[n-1]|$ is very small, this equation becomes ill-conditioned, similar to dividing by zero. To remedy this issue, we define a tolerance, below which $|x[n] - x[n-1]|$ is treated as if it were zero, and then use the following equation to express first-order ADAA: $$y[n] = \begin{cases} \frac{F_1(x[n]) - F_1(x[n-1])}{x[n] - x[n-1]}, & |x[n] - x[n-1]| > TOL \\ f\left( \frac{x[n] + x[n-1]}{2} \right), & \text{else} \end{cases}$$

For the hard-clipper, we can write the original nonlinear function as follows: $$f(x) = \begin{cases} x, & -1 \leq x \leq 1 \\ \text{sgn}(x), & \text{else} \end{cases}$$

And the first antiderivative as: $$F_1(x) = \begin{cases} \frac{1}{2}x^2, & -1 \leq x \leq 1 \\ x \; \text{sgn}(x) - \frac{1}{2}, & \text{else} \end{cases}$$

Let's see how the aliasing artifacts look using first-order ADAA:

In [6]:
class ADAA_1:
def __init__(self, f, F1, TOL=1.0e-5):
self.TOL = TOL
self.f = f
self.F1 = F1

def process(self, x):
y = np.copy(x)
x1 = 0.0
for n, _ in enumerate(x):
if np.abs(x[n] - x1) < self.TOL: # fallback
y[n] = self.f((x[n] + x1) / 2)
else:
y[n] = (self.F1(x[n]) - self.F1(x1)) / (x[n] - x1)
x1 = x[n]
return y

def signum(x):
return int(0 < x) - int(x < 0)

def hardClip(x):
return x if np.abs(x) < 1 else signum(x)

return x * x / 2.0 if np.abs(x) < 1 else x * signum(x) - 0.5

freqs, fft = process_nonlin(FC, FS, hardClip_ADAA.process)
plt.plot(freqs_alias, fft_alias, '--', c='orange', label='No ADAA')
plt.legend()
plt.scatter(freqs_analog[peak_idxs], fft_analog[peak_idxs], c='r', marker='x')
plt.xlim(0, 20000)
plt.ylim(10)
plt.title('Hard Clipping Distortion w/ 1st-order ADAA')
plt.ylabel('Magnitude [dB]')
plt.xlabel('Frequency [Hz]')
plt.grid()


Clearly, the aliasing artifacts are significantly more supressed than with no anti-aliasing, however they are still of significant magnitude. One option would be to use first-order ADAA in tandem with a modest amount of oversampling, maybe 2x:

In [7]:
freqs_2x, fft_2x = process_nonlin(FC, FS*2, hardClip_ADAA.process)
plt.legend()
plt.scatter(freqs_analog[peak_idxs], fft_analog[peak_idxs], c='r', marker='x')
plt.xlim(0, 20000)
plt.ylim(10)
plt.title('Hard Clipping Distortion w/ 1st-order ADAA')
plt.ylabel('Magnitude [dB]')
plt.xlabel('Frequency [Hz]')
plt.grid()


A second option is to use higher-order antiderivatives. Let's say that our nonlinear function $y[n] = f(x[n])$ has a second antiderivative, which we'll call $F_2(x)$. Then second-order ADAA can be written as: $$y[n] = \frac{2}{x[n] - x[n-2]} \left(\frac{F_2(x[n]) - F_2(x[n-1])}{x[n] - x[n-1]} - \frac{F_2(x[n-1]) - F_2(x[n-2])}{x[n-1] - x[n-2]} \right)$$

Since this equation has more divisions where a divide-by-zero error is possible, several "fallback" computations are necessary. I won't go through them in detail here, but full derivations can be found in Stefan Bilbao's IEEE paper (linked above), and implementations can be seen in the attached code (see below). Finally, we need the second antiderivative of our hard clipping function: $$F_2(x) = \begin{cases} \frac{1}{6}x^3, & -1 \leq x \leq 1 \\ \left(\frac{1}{2}x^2 + \frac{1}{6} \right) * \text{sgn}(x) - \frac{x}{2}, & \text{else} \end{cases}$$

Now we can examine the response of second-order ADAA:

In [8]:
class ADAA_2:
def __init__(self, f, F1, F2, TOL=1.0e-5):
self.TOL = TOL
self.f = f
self.F1 = F1
self.F2 = F2

def process(self, x):
y = np.copy(x)

def calcD(x0, x1):
if np.abs(x0 - x1) < self.TOL:
return self.F1((x0 + x1) / 2.0)
return (self.F2(x0) - self.F2(x1)) / (x0 - x1)

def fallback(x0, x2):
x_bar = (x0 + x2) / 2.0
delta = x_bar - x0

if delta < self.TOL:
return self.f((x_bar + x0) / 2.0)
return (2.0 / delta) * (self.F1(x_bar) + (self.F2(x0) - self.F2(x_bar)) / delta)

x1 = 0.0
x2 = 0.0
for n, _ in enumerate(x):
if np.abs(x[n] - x1) < self.TOL: # fallback
y[n] = fallback(x[n], x2)
else:
y[n] = (2.0 / (x[n] - x2)) * (calcD(x[n], x1) - calcD(x1, x2))
x2 = x1
x1 = x[n]
return y

return x * x * x / 6.0 if np.abs(x) < 1 else ((x * x / 2.0) + (1.0 / 6.0)) * signum(x) - (x/2)