Fast Synthesis of Basic Oscillators
I’ve been working on synthesizers for the GBA and rp2040 the past couple years, and I’ve made some fast oscillators.
So first of all, this all applies to doing synthesis with integers. I think for oscillators, even if you do have an FPU, ints are much easier to deal with anyway because of the wrapping.
We’ll be doing math entirely with 32-bit values in this post, for two reasons. One is, it’s the native register size on the GBA, and what I’ve been doing with for the past couple years. The other: it’s a very nice sweet spot for pitch resolution. You can represent all pitches within the audible range quite accurately when dealing with a 32-bit phase. As you reduce the bit-depth of your phase, you start losing pitch accuracy in the bass. You can work with that, and many demosceners do, but it becomes something you have to actively think about when composing your music. With 32-bit, its a non-issue. Anything beyond 32-bit doesn’t really provide much benefit.
Most of the waveforms I’m showing here will have aliasing problems. The sharp transitions create audible artifacts that don’t sound super great. One way to fix this is with PolyBLEP.
I’ve never implemented that, and anyway the point of this code is to be fast at all cost. This is for super constrained environments, where you have to accept the aliasing is there and move on. For alias-free sounds on a budget, you can use phase modulation between sine waves and you’ll get super clean sounds for very cheap, though not as cheap as these waves. I’ll do another post about that. If I ever do a fast PolyBLEP I’ll post about that too, but it’s not on my to-do right now.
Saw
All oscillators start life as a saw wave, and a saw wave has two components: the current phase, and a value added to the phase each audio sample.
phase: u32 = 0
step: u32 = (frequency_in_hz * 0x100000000u64 / sample_rate) as u32
outbuf: [i32; num_samples]
Each audio sample, we add the step to our phase, and the output of our oscillator is the phase. This gives us a ramp-up saw.
let phase: u32 = 0;
let step: u32 = (frequency_in_hz * 0x100000000u64 / sample_rate) as u32;
let mut outbuf: [i32; num_samples];
for i in 0..num_samples {
phase = phase.wrapping_add(step);
outbuf[i] = phase as i32;
}
// Arm assembly
// loop:
// add phase, step
// str phase, [outbuf], 4
// subs count, 1
// bne loop
1 / /
/| /|
/ | / |
/ | / |
0 / | / |
| / | /
| / | /
|/ |/
-1 / /
Each time phase overflows, we’ve completed one oscillation.
If you’re familiar with the concept of the nyquist
frequency, it’s quite obvious
here: when our phase step is 0x80000000, our wave oscillates once every
two-samples, at exactly 1/2 our sample rate. If our phase step is higher than
0x80000000, we start overflowing more than once every two samples, but our
output actually starts slowly rolling backwards. This also flips our saw wave
to a ramp-down saw instead of a ramp-up. Taken as a signed integer, phases >=
0x80000000 are negative. Naturally, as our phase step grows beyond
0x80000000 we start oscillating slower and slower.
Wave from memory
If you shift down your saw, you can use it to index a wave stored in memory, and do all the wave table synthesis your heart desires. You’ll get the best performance with 8-bit waves, because it takes an extra instruction to mask off the low address bits when accessing 16 or 32-bit waves. So I usually use 8-bit waves, as I don’t find waveforms benefit from the dynamic range of 16-bits the way samples do.
This example uses 1024-sample waves, adjust the bit shifts for other sizes. Make sure to store your wave somewhere with fast access times.
let phase: u32 = 0;
let step: u32 = (frequency_in_hz * 0x100000000u64 / sample_rate) as u32;
let mut outbuf: [i32; num_samples];
let wavedata: [i8; 1024] = /* ... */
for i in 0..num_samples {
phase = phase.wrapping_add(step);
let out = wavedata[(phase >> 22) as usize];
// Scale up the output if you need to, or
// leave it 8-bit to save some cycles cycles
// going into a fixed-point volume control :)
outbuf[i] = out;
}
// Arm assembly
// loop:
// add phase, step
// lsr out, phase, 22
// ldrsb out, [wavedata, out]
// str out, [outbuf], 4
// subs count, 1
// bne loop
Square
There’s two approaches you can take depending on personal taste, and what’s fast for your hardware.
Mathematic calculation
Calculating mathematically is best if you don’t have branchless conditionals. This is true of thumbv1, like you find on common Cortex-M0+ processors. It’s fine everywhere else too, so reach for this first.
let phase: u32 = 0;
let step: u32 = (frequency_in_hz * 0x100000000u64 / sample_rate) as u32;
let mut outbuf: [i32; num_samples];
for i in 0..num_samples {
phase = phase.wrapping_add(step);
// out = 0x00000000 when phase < 0x80000000
// out = 0xFFFFFFFF when phase >= 0x80000000
let mut out = (phase as i32) >> 31;
// out = INT_MIN when phase < 0x80000000
// out = INT_MAX when phase >= 0x80000000
out = out + (-0x80000000);
outbuf[i] = out;
}
// arm assembly
// loop:
// add phase, step
// asr out, phase, 31
// sub out, 0x80000000
// str out, [outbuf], 4
// subs count, 1
// bne loop
// thumbv1 assembly (cannot subtract large immediates like arm)
// movs dc_offset, 0x80
// lsls dc_offset, 24
// loop:
// adds phase, step
// asrs out, phase, 31
// subs out, dc_offset
// str out, [outbuf], 4
// subs count, 1
// bne loop
1 _________ _________
| | | |
| | | |
| | | |
0 | | | |
| | | |
| | | |
| | | |
-1 ________| |________| |
For thumbv1, you can sacrifice a bit of precision and a cycle for each instruction to get rid of the dc offset register:
// loop:
// adds phase, step
// asrs out, phase, 31
// subs out, 0x80
// lsls out, 24
// str out, [outbuf], 4
// subs count, 1
// bne loop
Conditional calculation
This second algorithm is only fast if you have branchless-conditionals. arm32
can make any instruction conditional, and thumbv2 has the IT instruction to
form a small branchless conditional block, so this algorithm is good on both of
those, and it’s what I use for the GBA. The main advantage over the mathematic
calculation is it’s perfectly centered around 0, and it’s easy to control the
amplitude of the square coming out of it by adjusting the constants. It may or
may not fit better into the signal chain around it depending on what you’ve got
going on.
let phase: u32 = 0;
let step: u32 = (frequency_in_hz * 0x100000000u64 / sample_rate) as u32;
let mut outbuf: [i32; num_samples];
for i in 0..num_samples {
phase = phase.wrapping_add(step);
// We use 0x7F000000 instead of 0x7FFFFFFF,
// because it takes fewer cyles to load a
// constant with <= 8 sigfigs on arm/thumb
if ((phase & 0x80000000) == 0) {
outbuf[i] = -0x7F000000;
else {
outbuf[i] = 0x7F000000;
}
outbuf[i] = phase as i32;
}
// Arm assembly. thumbv2 would be similar, but use `IT`
// loop:
// adds phase, step
// movpl out, -0x7F000000
// movmi out, 0x7F000000
// str out, [outbuf], 4
// subs count, 1
// bne loop
Triangle
To explain the theory here I’m going to lead with an implementation using conditionals, but then I’ll give you a conditional-free variant that works the same way but is a little less clear as to what’s going on.
Uni-polar triangle wave with Conditionals
A triangle wave ramps up and then back down again. This is rather annoying to express in floating point, but here in integer land it’s quite easy:
let phase: u32 = 0;
let step: u32 = (frequency_in_hz * 0x100000000u64 / sample_rate) as u32;
let mut outbuf: [i32; num_samples];
for i in 0..num_samples {
phase = phase.wrapping_add(step);
if (phase & 0x80000000 == 0) {
outbuf[i] = phase as i32;
} else {
outbuf[i] = (!phase) as i32;
}
}
// Arm assembly
// loop:
// adds phase, step
// movpl out, phase /* count-up in +phase */
// mvnmi out, phase /* count-down in -phase */
// str out, [outbuf], 4
// subs count, 1
// bne loop
1 - -
/ \ / \
/ \ / \
/ \ / \
0 - - -
-1
Note that ! is the one’s complement operator, meaning all the bits in phase
get flipped.
This ramps up normally from 0 to 0x7FFFFFFF. Once we hit 0x80000000, we
start inverting the number. !0x80000000 = 0x7FFFFFFF, !0x80000001 =
0x7FFFFFFE, and so on, counting all the back down to !0xFFFFFFFF = 0.
Note that this means our triangle has very slightly flat tips, since 0 appears
twice, and so does 0x7FFFFFFF. When dealing with 32-bit phase, you will never
actually hear this artifact, but it might become more apparent with smaller
phases.
Taken as a signed integer, !x = (-x) - 1. You might wonder, why don’t we use
the actual - operator to take the twos complement, or use the .abs() absolute value function? It seems like
it would get rid of our flat tips! Well, it almost does, except there is a big
problem: INT_MIN.abs() == INT_MIN. The twos-complement of 0x80000000 is 0x80000000. When we output a signed audio sample, our wave would have a very nasty wrap-around negative spike at the tip, like this:
1 | |
/|\ /|
/ | \ / |
/ | \ / |
0 / | \/ |
| |
| |
| |
| |
-1 | |
Bi-polar triangle wave with Conditionals
Outputting a uni-polar wave can be fine sometimes, but usually you want a wave that oscillates between negative and positive. Our solution to that is crude. We just subtract the DC-offset from our triangle, and there you have it.
let phase: u32 = 0;
let step: u32 = (frequency_in_hz * 0x100000000u64 / sample_rate) as u32;
let mut outbuf: [i32; num_samples];
for i in 0..num_samples {
phase = phase.wrapping_add(step);
if (phase & 0x80000000 == 0) {
outbuf[i] = (phase as i32) - 0x40000000;
} else {
outbuf[i] = ((!phase) as i32) - 0x40000000;
}
}
// Arm assembly
// loop:
// adds phase, step
// movpl out, phase
// mvnmi out, phase
// sub out, 0x40000000
// str out, [outbuf], 4
// subs count, 1
// bne loop
1
- -
/ \ / \
0 / \ / \
/ \ / \
- - -
-1
This produces a triangle wave which does not span the full i32 range, instead spanning from -0.5 to +0.5. That’s never an issue in practice, since you’re always going into some sort of volume control after generating the oscillator, and you can account for the loss in volume somewhere in that signal chain.
Mathematic calculation of the triangle
The conditional implementation is easy to understand, but you can do the exact same thing mathematically, in fewer instructions.
let phase: u32 = 0;
let step: u32 = (frequency_in_hz * 0x100000000u64 / sample_rate) as u32;
let mut outbuf: [i32; num_samples];
for i in 0..num_samples {
phase = phase.wrapping_add(step);
// Will be all-zeroes if phase is positive or all-ones if phase is negative
let mut out = (phase as i32) >> 31;
// When phase is positive, out = phase
// When phase is negative, out = !phase
out = out ^ (phase as i32);
// Remove DC-offset
outbuf[i] = out - 0x40000000;
}
// arm assembly
// loop:
// add phase, step
// eor out, phase, phase, asr 31
// sub out, 0x40000000
// str out, [outbuf], 4
// subs count, 1
// bne loop
Parabola approximation of a sine wave
I’m stretching the limits of “basic” here, but it’s an approximation of a basic wave so we’ll do it for fun.
If you have a fast 32-bit multiply instruction, you can get an approximation of a sine wave by taking a parabola for the first half of the phase, and then the negation of that parabola for the second half of the phase.
This approximation is nowhere close to accurate enough for physics or graphics, but for audio it will get you a sound that is closer to a pure tone than a triangle wave, without being as pure as a proper sine.
The usecase for this is very niche. On the GBA, I don’t do this, because it’s slower than just accessing a sine lookup table in ROM. On the rp2040 I don’t do this, because while ROM access is slow, there’s plenty of SRAM to keep a sine table in. On a high end system, there’s so much cache that you can pre-load a sine table into cache if you’re chasing perf instead of doing this.
Best I can tell the only real reason to use this is if you’re doing some size coding and don’t have space for a lookup table, or you want the unique sound of it.
let phase: u32 = 0;
let step: u32 = (frequency_in_hz * 0x100000000u64 / sample_rate) as u32;
let mut outbuf: [i32; num_samples];
for i in 0..num_samples {
phase = phase.wrapping_add(step);
// Saw wave at 2x frequency
let mut out: i32 = (phase << 1) as i32;
// Saw -> Unipolar triangle (see Triangle section)
// Range 0 - 0x7FFFFFFF
let mut out = out ^ (out >> 31);
// We are going to calculate X - X^2 to get a
// parabolic arc from 0 to half-phase. This is
// can also be thought of as
// lerp(from: X, to: 0, alpha: X)
let x = out;
// Scale to 0 - 0x7FFF
let mut out = phase >> 16;
// Output arc from 0 - 0x4000FFFE
out = out * out;
out = x - out
// Negate the second arc in the second half of
// the phase
out = out ^ ((phase as i32) >> 31);
// Final output: -0x4000FFFF - 0x4000FFFE
outbuf[i] = phase as i32;
}
// Arm assembly
// loop:
// add phase, step
// lsl out, phase, 1
// eor out, out, out, asr 31
// mov x, out
// lsr out, 16
// mov tmp, out /* arm limitation: cant square
// mul out, tmp with single reg */
// sub out, x, out
// eor out, out, phase, asr 31
// str out, [outbuf], 4
// subs count, 1
// bne loop
It’s hard to do this one justice with ascii-art, so here’s a proper render:
