Fast Phasor ala R. Hoelderich

Written by Nick Donaldson, 2013

The following is a description of the "fast phasor" bit-twiddling algorithm (supposedly) described by R. Hoelderich in the 1995 ICMC Banff proceedings. I have as of yet been unable to find the original paper from these proceedings or any other documentation on the algorithm, but this is an explanation of how it works.

Purpose

In digital audio, there are a number of techniques that make use of a periodic "phasor", which is essentially a sampled version of a perfectly angular rising sawtooth (slopes linearly upward to a maximum value, then drops back to zero, and repeats indefinitely).

(image of phasor signal might go here)

This type of phasor is commonly used in audio DSP for indexing a lookup table for oscillators, as an aliasing sawtooth, and for various other applications.

For a general-purpose phasor, the range of the phasor is typically [0-1) (fractional values only). The phasor can be linearly scaled as necessary for a particular application. For table lookup oscillators, where it is useful to have both an integral value to index an array as well as a fractional value for interpolation, the phasor represents a sample index ranging from zero to the length of the table. Lookup tables are conventionally created with a power-of-two + 1 number of samples, where the extra sample makes linear interpolation more efficient when the base sample index is at the end of the table.

A functional yet somewhat inefficient approach that is commonly used in audio DSP applications is to simply increment a floating-point phase value on every subsequent sample, check if the phase has exceeded the maximum, and reset it back to zero if that is the case:

float phase;
while (samplesRemaining--){
  phase += increment;
  if (phase >= max){
    phase -= max;
  }
  /* Do something with phase value */
}

The problem here is that the branch instructions for the "if" statement in a tight loop that repeats 44.1k+ times per second are very inefficient, as the branch will only evaluate "true" a relatively small percentage of the time. The wrapping can be made faster for lookup tables if the int floor of the phasor is taken first, and its modulo is evaluated against the length of the table, but this still requires a somewhat costly cast to an integer, as explained below.

For use as a table lookup oscillator, the inefficiency is compounded by the fact that interpolating lookup tables need to derive the integer portion of the phase for indexing into an array. A typical approach is to floor the phase:

float phase;
while (samplesRemaining--){

  phase += increment;
  if (phase >= max){
    phase -= max;
  }
  
  int idx = (int)phase;
  float frac = phase - (float)idx;

  /* Do something with phase/idx/frac value */
}

The implicit floor by casting to int is often faster than calling floor(), but this is still somewhat inefficient.

The "fast phasor" algorithm described below is a way to use the bits of a double-precision floating point number to quickly and efficiently perform "wrapping" on a unit-range phasor and also derive both the integral and fractional parts of the output of a phasor with a power-of-two range.

The Algorithm

(TODO)

References

Mention on linux-audio-dev mailing list archive
Implementation in pd source code