[sdiy] BLIT/BLEP virtual analogue synthesis

[Olivier Gillet]

Another idea from two observations: 1/ you're probably computing the
samples at a faster rate than the control signal modulating
phase_increment. 2/ What really matters is the scaled reciprocal of
phase_increment. Scaled as in 65536 / phase_increment

So here's a way of doing it: whenever you do something that ends up
modify phase_increment, set a "dirty" bit to true.

Now in your blep code:

if (/*wrapping*/  ... ) {
  phase = ...
  if (dirty) {
    phase_increment_reciprocal = ... // Smart code to compute 65536 /
phase_increment  efficiently
    dirty = 0;
  }
  fractional_wrap_point = phase * phase_increment_reciprocal >> 8;  //
Assuming you want the fractional value in the 0-255 scale.

  // Add a blep parametrized by fractional_wrap_point to the playlist
}

The reciprocal computation can be done with a Newton method, or using
the following shift + lookup table idea:

# 256 bytes LUT with reciprocal of scaled large values
reciprocal_lut = [0] * 128
for i in xrange(128):
  reciprocal_lut[i] = int(65536 / ((i + 128.5)))

# 512 bytes LUT with reciprocal of small values
small_values_lut = [0] * 256
for i in xrange(256):
  small_values_lut[i] = 65536 / i if i > 0 else 1

# The pseudo-code for the reciprocal itself
def FastScaledReciprocal(x):
  if x < 256:
    return small_values_lut[x]
  else:
    r = 1
    p = x
    while p < 32768:
      r <<= 1
      p <<= 1
    return r * reciprocal_lut[(p >> 8) - 128] >> 8

I gave it a test run and the max error is 1, average is 0.0066, and
there's at most 7 iterations in the loop. It probably beats a
division.


Now that I've written it, I recognized the pattern: this is the same
flavour of code as the shift+antilog trick used to convert a
fractional MIDI pitch into a phase increment. For instance:
http://github.com/pichenettes/shruthi-1/blob/master/hardware/shruthi/synthesis_engine.cc#L893

So you already have this kind of lookup table and loop somewhere, and
you use it to compute the phase increment. Thus, in your MIDI >
increment loops, do the exact dual thing, in the opposite direction,
and you'll get the reciprocal for free. With the reciprocal, the
fractional wrapped_phase / phase_increment ratio is only a multiply
and shift away.

Olivier