Differentiated Parabolic Wave

DPW is a simple, efficient and good way to reduce aliasing in a sawtooth wave oscillator. For a proper treatment, please refer to the literature. Otherwise, we'll just indulge ourselves in some sawtooth goodness.

DPW, like many of its brethern, starts with a humble phasor. While the phasor will also produce sawtooth waves, it exhibits some nasty aliasing behavior; the example below is designed to showcase that ugly distortion. Note how some of the harmonics are moving in the "wrong" direction!

Import Gen

; slow modulation to emphasize non-harmonic alias terms
sweep = 4000 + 1000 * Gen:Sin(0.2)

snd = Gen:Phasor(sweep) * 0.1

This is because we're sampling an infinitely sharp sawtooth at the lowly sample rate of your audio hardware. The components at and above the Nyquist frequency, that is, half your sample rate, will be played back at an entirely wrong frequency, a bit like the wheels that seem to spin backwards on a limited frame rate video.

The core of bandlimited synthesis is to remove or cancel out those components before sampling. We could always apply a low-pass filter. However, then we'd need to run the phasor and a filter at a very high internal sample rate.

Turns out we can cheat if we know our source waveform. Integration is a filter, and we can integrate the phasor formula before sampling!

See the paper for the gory details. The integrated shape of a linear segment is a quadratic segment. Indeed, if we just square the phasor, some filtering can be heard.

phasor = Gen:Phasor(440) * 2 - 1
parabolic = phasor * phasor

; fade back and forth between phasor and squared phasor

fade = Gen:Tri(0.2) * 0.5 + 0.5 ; unipolar triangle osc
snd = (phasor + fade * (parabolic - phasor)) * 0.2

So - we filter by 6dB/oct in the continuous domain. This is enough to suppress the aliasing quite a bit. Our waveform became duller though. Turns out we can reverse the filtering in the discrete domain by a simple differentiator.

Differentiator(sig) {
	; just substract a delayed signal from signal!
	sig - z-1(sig)
}

snd = Differentiator(parabolic) * 0.2

Just one problem! The level of the signal varies with frequency, as we readily hear if we plug our original sweep into the DPW:

naive = Gen:Phasor(sweep) * 2 - 1
snd = Differentiator(naive * naive) * 0.1

The reason for that is the fact that the differentiator has a magnitude response that changes with frequency. We want exactly that for the correct harmonic balance, but we must fix the level at the fundamental.

Deriving the differentiator response:

y[t] = x[t] - x[t-1]
Y(Z) = X(Z) - X(Z) Z^-1
H(Z) = Y(Z) / X(Z) = 1 - Z^-1

let Z = exp(i w)
|H(Z)| = sqrt((1 - exp(-i w))(1 - exp(i w)))
|H(Z)| = 2 - 2 cos(w)

So, the response of the differentiator at frequency w (radians) is 2 - 2 cos(w). To fix the level of the oscillator, we just compensate for the response at f0.

Import Math

DPW(freq) {
	w = freq * 2 * Math:Pi / Gen:Rate(nil)
	f0g = 2 - 2 * Math:Cos(w)
	phasor = Gen:Phasor(freq) * 2 - 1

	Differentiator(phasor * phasor) / f0g
}


; steady DPW
snd = DPW(sweep) * 0.1

; demonstrate aliasing by switching back and forth
bad  = Gen:Phasor(sweep) * 2 - 1
good = DPW(sweep)

switch = Gen:Phasor(0.2) > 0.5

snd = Algorithm:Choose(switch good bad) * 0.1

The remaining aliasing we can suppress further with help from oversampling, which synergizes well with the DPW strategy.

; one stage of oversampling for 2x audio rate.
better = Gen:Oversample(#1 { DPW(sweep) } nil)

snd = Algorithm:Choose(switch better good) * 0.1

The nice thing about DPW is that adding phase offsets is easy; we just change the phasor part, the alias supression stays the same. With that in mind, we can build a variable width pulse wave as the difference between two phase-offset sawtooths.

DPW(freq offset) #[Extend] {
	w = freq * 2 * Math:Pi / Gen:Rate(nil)
	f0g = 2 - 2 * Math:Cos(w)

	; add phase offset here and use Fraction to wrap the phasor back to
	; the [0,1] range.
	phasor = Fraction(Gen:Phasor(freq) + offset) * 2 - 1

	Differentiator(phasor * phasor) / f0g
}

width = Gen:Phasor(0.2)

snd = (DPW(55) - DPW(55 width)) * 0.1

Despite its simplicity, DPW is a solid method for virtual analog synthesis. Poly-BLEP would be a worthy contender. Choices, choices!