The math of bandlimited oscillators

In the previous post we constructed a bandlimited (non-aliasing) oscillator modeled after the impulse train. In this post we'll look at more aliasing free oscillators from the mathematical point of view.

The previous post looked at the idealized impulse train worked our way down to a useful digital implementation. Let's try it the other way around. We'll start with the discrete impulse train that has one sample with amplitude 1 and (n-1) samples with amplitude 0.
[1, 0, 0, 0, ... , 0], repeated ad infinitum.
The frequency content can be calculated by taking the discrete Fourier transform. It is trivial because each of the partials we're summing over in the transformation starts with 1 + 0i regardless of frequency.
fft([1, 0, ... , 0]) = [1, 1, ... , 1].
So the discrete impulse train contains all the complex frequencies in equal proportion.
( ... + exp(-i x) + exp(0) + exp(i x) + exp(i 2x) + ... )/n
Notice how every complex exponential has a pair except for exp(0). If n is even the Nyquist frequency will also be missing a pair. So writing the series in terms of real functions we have:
( 1 + 2 cos(x) + 2 cos(2x) + ... + 2 cos((n/2-1)x) + cos((n/2)x) )/n, n even
( 1 + 2 cos(x) + 2 cos(2x) + ... + 2 cos(((n-1)/2)x) )/n, n odd.

There is a little trick to get simple expressions for the continuous counterparts of both the single impulse [..., 0, 0, 1, 0, 0, ....] and the impulse train.  The unique non-zero analytic function that has zeros exactly at the sampling points giving [..., 0, 0, 0, ....] is sin(π t). To get the continuous version of the single impulse we just divide out one of the zeros to get
sin(π t)/(π t) =: sinc(t)

as simple and irrigorous as that. For the impulse train we simply divide out every n:th zero with sin(π t/n) to get:
sin(π t)/(n sin(π t/n)) =: windowedsinc(t,n), n odd.
The catch is that this trickery only works when n is odd. But this is the function we'll want to use because the odd case had a more uniform spectrum without any nyquist nuisances. For reference the even case simplifies to:
sin(π t)cot(π t/(2n))/n.
It has a double analytical zero between the impulses. Not something you can easily just deduce by matching zeros with samples. But what happens when you try windowedsinc with an even n?

Oh, so it flips every other impulse around. Well that's cool. We get the difference of two impulse trains 180° apart for free in terms of computational effort. This has the effect of killing every even harmonic leaving us something that sounds like a thin square wave. Actually it is the derivative of the bandlimited square wave, but more on that later.

You can also make bandlimited oscillators by dividing out more zeros per (semi-)period like getting rid of the zeros at indices 0, 1 and 5 mod 10, but they may be more trouble than they're worth. The point of these formulas is to make it easy to calculate a known sum of harmonics without actually doing the summation. Getting the harmonic content of these beasts is somewhat involved.

The computational gain is also marginal because they can always be reduced to a sum of the regular blits.
In the previous post we used the geometric series. Let's look at it in more detail:
1 + z + z² + ... + z^(n-1) = (1-zⁿ)/(1-z).
Let's plug in e^2πit for z and see what happens.

I like the imaginary part better. It's more symmetric. Anyway let's replace z with az, where a is a real number say 0.75.

Hey that imaginary part sort of looks like the saw wave! Well it sounds like it too but much softer. The harmonics fall off more quickly like aⁿ instead of 1/n like in the saw wave. To get something resembling the square wave we'd need to eliminate the even powers of z. However when dealing with complex harmonics there is an easier way. We double the frequency by replacing az with (az)² and then multiply with az to shift the powers to the odd indices.

Again the visual similarity to the square wave is apparent in the imaginary part. The bandlimited gifts of the geometric series do not end here. We can differentiate the series to get:
1 + 2z + 3z² + ... + (n-1)z^(n-2) = (z + (z(n-1)-n)zⁿ)/(z(1-z)²).
A bit messier than before but still much better than actually carrying out the summation on the left hand side. Again plugging in az for z gives a way to control the softness of the sound and for a<1 the spectrum fill have a hump that can be used in combinations with other oscillators to give more color to the sound.

The final trick that can be used to get more sounds out of the geometric series is multiplying it by e^iωt to shift the harmonic frequencies a non-integer amount up or down to get strange aperiodic sounds.
Unfortunately the geometric series cannot be integrated to get cheap formulas for the bandlimited saw and square wave. The integral is:
z + z²/2 + z³/3 + ... + zⁿ/n = -Beta(z; 1+n, 0) - log(z-1),
Where Beta is the incomplete Beta-function. And the right hand side is not defined when z=1. The calculation of Beta is considerably more involved than simple multiplication and division but with the right libraries it might still be computationally cheaper than carrying out the sum naively.