#+TITLE: The Oscillator* Introduction.This document will describe an implementation of a classictable-lookup oscillator with linear interpolation. The algorithm for this is an interesting mix of numericalprocessing, with the phasor being implemented in fixedpoint, and everything else being done in floating point.A big advantage to doing it this way is numerical stability:there is no risk of any phase accumulation or truncationlike you'd get with floating-point. Thisis an important characteristic needed for phasor to avoidthings like drift over time. As will be seen later, thefixed-point approach can be a little bit harder tounderstand, especially for those unfamiliar with fixed-pointDSP. The author has put extra effort to explain theseportions as clearly as possible. The implementation also has a famous limitation of onlybeing able to take in tables with sizes that are a power of2.* Generated FilesHeader: #+NAME: osc.h#+BEGIN_SRC c :tangle osc.h#ifndef SK_OSC_H#define SK_OSC_H#ifndef SKFLT#define SKFLT float#endif<<typedefs>>#ifdef SK_OSC_PRIV<<structs>>#endif<<funcdefs>>#endif#+END_SRC C file: #+NAME: osc.c#+BEGIN_SRC c :tangle osc.c#include <stdint.h>#include <stdlib.h>#include <math.h>#define SK_OSC_PRIV#include "osc.h"<<constants>><<funcs>>#+END_SRC* Top-level functionsThe table-lookup oscillator is initialized with=sk_osc_init=. The following arguments must be provided: =osc= is a pre-allocated struct of =sk_osc=. =sr= is the sampling rate. =tab= is a pre-allocated wavetable, an array of =SKFLT=floating-point values. =sz= is the array size of the wavetable =wt=. =iphs= provides the initial phase of the oscillator. It isa value between 0-1. #+NAME: funcdefs#+BEGIN_SRC cvoid sk_osc_init(sk_osc *osc, int sr, SKFLT *wt, int sz, SKFLT iphs);#+END_SRC To compute a sample of audio, use =sk_osc_tick=. #+NAME: funcdefs#+BEGIN_SRC cSKFLT sk_osc_tick(sk_osc *osc);#+END_SRC This oscillator has 2 main parameters: frequency (freq),and amplitude (amp). They can beset with the following functions. #+NAME: funcdefs#+BEGIN_SRC cvoid sk_osc_freq(sk_osc *osc, SKFLT freq);void sk_osc_amp(sk_osc *osc, SKFLT amp);#+END_SRC #+NAME: funcs#+BEGIN_SRC cvoid sk_osc_freq(sk_osc *osc, SKFLT freq){    osc->freq = freq;} void sk_osc_amp(sk_osc *osc, SKFLT amp){    osc->amp = amp;}#+END_SRC* Constants#+NAME: constants#+BEGIN_SRC c#define SK_OSC_MAXLEN 0x1000000L#define SK_OSC_PHASEMASK 0x0FFFFFFL#+END_SRC* Struct and ConstantsThe main struct of this oscillator is called =sk_osc=.#+NAME: typedefs#+BEGIN_SRC ctypedef struct sk_osc sk_osc;#+END_SRC#+NAME: structs#+BEGIN_SRC cstruct sk_osc {    <<sk_osc>>};#+END_SRC #+NAME: sk_osc#+BEGIN_SRC cSKFLT freq, amp;SKFLT *tab;int inc;size_t sz;uint32_t nlb;SKFLT inlb;uint32_t mask;SKFLT maxlens;int32_t lphs;#+END_SRC The oscillator stores it's main parameters =freq= and =amp=as floating point parameters. They are set to be values440 and 0.2 by default. #+NAME: osc_init#+BEGIN_SRC cosc->freq = 440;osc->amp = 0.2;#+END_SRC A reference to the table is stored in the variable =tab=,with its size =sz=. #+NAME: osc_init#+BEGIN_SRC cosc->tab = wt;osc->sz = sz;#+END_SRC A table lookup oscillator indexes through the table usingthe increment rate stored in the integer value =inc=. Thisvalue can be positive or negative. Is is zeroed out atinit-time. #+NAME: osc_init#+BEGIN_SRC cosc->inc = 0;#+END_SRC The variable =lphs= stores the phase position of theprevious sample. The initial phase value =iphs= ismultiplied with the the maximum table value, and thenmasked to keep values in range. #+NAME: osc_init#+BEGIN_SRC cosc->lphs = ((int32_t)(iphs * SK_OSC_MAXLEN)) & SK_OSC_PHASEMASK;#+END_SRC For the fixed point table-lookup, some constants are derivedand stored. Phasor position is stored by splitting the bits of anN-bit integer number into two parts. The upper bits store theinteger portion, while the lower bits store fractionalportion. The maximum number of bits is arbitrary, but theunderlying architecture must be able to accomodate for thewidth. In this implementation, the phasor uses 28 bitsinside of a 32-bit number. This implicitely means the largestvalue can be Split in the phasor position is measured by counting thenumber of lower bits. This value is stored in the variable=nlb=. This value is calculated with the equation @!(fig "osc_log2" "n = log2(M/s)")!@ Where =n= is the number of lower bits, =M= is the maximumwavetable size, and =s= is the size of the wavetable. To calculate =nlb=, and hand-rolled =log2= function iscreated. Smaller values of =s= mean more bits in the fractionalcomponent of the number. #+NAME: osc_init#+BEGIN_SRC c{    uint32_t tmp;    tmp = SK_OSC_MAXLEN / sz;    osc->nlb = 0;    while (tmp >>= 1) osc->nlb++;}<<calculate_mask>><<calculate_inlb>><<calculate_maxlens>>#+END_SRC The =mask= is the lower-bits masking variable. When an ANDoperation is used against this constant, it filters out allthe upper bits, so only the lower bits can pass through.This constant is necessary for being able to extract thelower bits from the fixed-point phase value representation.This sort of value is known in the bit-twiddling world as a=mask=. In binary, all the lower bits up to the number oflower bits are set to be on, with the remaining bits setto be 0. #+NAME: calculate_mask#+BEGIN_SRC cosc->mask = (1<<osc->nlb) - 1;#+END_SRC The inverse of max lower bits value, or $1/2^{nlb}$, isstored as a constant. This cached value is used to replaceand divide operation with a multiply operation, which hastraditionally been a cheaper operation to do on a computer. #+NAME: calculate_inlb#+BEGIN_SRC cosc->inlb = 1.0 / (1<<osc->nlb);#+END_SRC The constant =maxlens= is the maximum table length in unitsof seconds. This is a value used to efficiently convert thefrequency parameter to sample increment value. #+NAME: calculate_maxlens#+BEGIN_SRC cosc->maxlens = 1.0 * SK_OSC_MAXLEN / sr;#+END_SRC* InitializationIn addition to setting variables, the init function willalso set the starting phase of oscillator.#+NAME: funcs#+BEGIN_SRC cvoid sk_osc_init(sk_osc *osc, int sr, SKFLT *wt, int sz, SKFLT iphs){    <<osc_init>>}#+END_SRC* ComputationThe meat of the algorithm is here. Here outlines the=tick= function, where a single sample of an oscillatoris computed. #+NAME: funcs#+BEGIN_SRC cSKFLT sk_osc_tick(sk_osc *osc){    SKFLT out;    SKFLT fract;    SKFLT x1, x2;    int32_t phs;    int pos;     out = 0;    <<update_increment_amount>>    <<lookup_values>>    <<obtain_fractional_component>>    <<interpolate_values>>    <<update_the_state>>    return out;}#+END_SRC First, the increment amount =inc= is updated.The increment amount tells how much further to move the readpointer in the table. This increment amount is based on thecurrent oscillator frequency =freq= and the variable=maxlens=. =lrintf= is used to round to the nearest integer. This is kind of a baffling operation. How could multiplyingthe frequency by some arbitrary duration yield a phasorincrement amount? And where is the sampling rate in allthis? The thing that throws everything off is this fixed pointbusiness. Things make a lot more sense if you wanted todo this the floating-point way. Recall that a phasor is a repeating line that ramps from0 to 1 over a given period of time. If we call this periodof time =t=, the increment value =I= needs to work so that=t * sr * I = 1=. In other words, it's the slope of a linediscretised. Linear slope is found using good-ol rise over run, change invalue over change in time. Inverting the frequency =1/F=gives it's period length in seconds. Multiplying by thesampling rate =sr= converts that value to samples. Thisgives us =sr/f=, our change in time. Because of thenormalized range, change in value is just 1. Putting it alltogether This gives us a slope of =1/(sr/f)=, or =f/sr=.That is the normalized increment value. If we wanted this to work with our fixed point phasor,we'd need to scale it by the maximum length of the phasor.That looks like =(f/sr)*maxlen)=, or =(f * maxlen)/sr)=.The frequency =f= can be pulled out of the numerator toget =f * (maxlen/sr)=, which can be reduced to the operationsimilarly seen below of =f * maxlens=. Before frequency units were measured in Hertz, they weremeasured in cycles-per-second. If you ever took a highschoolchemistry or physics class, you may recall that units can"cancel out" one another like a fraction. Whencycles-per-second (cycles / seconds) gets multiplied bya value in seconds, the seconds cancel. What you are leftwith is a unit called *cycles*. #+NAME: update_increment_amount#+BEGIN_SRC cosc->inc = floor(lrintf(osc->freq * osc->maxlens));#+END_SRC Look up values A =n= and B =n + 1= samples from wavetable.Perform table lookup. Both the current position, and it'sneighor are needed.This position is found by looking at the upper bits ofthe current phase. #+NAME: lookup_values#+BEGIN_SRC cphs = osc->lphs;pos = phs >> osc->nlb;x1 = osc->tab[pos];x2 = osc->tab[(pos + 1) % osc->sz];#+END_SRC Now, it's time to interpolate between points A and B. Thisoscillator uses linear interpolation, which can be thoughtof as a crossfade between two values. The equation forlinear interpolation is commonly shown as: @!(fig "osc_lerp" "y = (1 - \\alpha) x_1 + \\alpha x_2")!@ Where $x_1$ and $x_2$ are two values, and $\alpha$ is afractional value between 0 and 1. The $\alpha$ valuesdetermines the distribution balance ofthe two values. When $\alpha = 0$, it is entirely $x_1$, andwhen $\alpha = 1$, the value is $x_2$. The fractional (alpha) value is obtained by taking the lowerbits portion of the current fixed-point phase position, andnormalizing it to be a floating-point value between 0 and 1. #+NAME: obtain_fractional_component#+BEGIN_SRC cfract = (phs & osc->mask) * osc->inlb;#+END_SRC There are now all the parts to do the interpolation. Itturns out the equation above can be simplified further toshave off a multiply operation. @!# $$# $$(fig "osc_lerp_simp"``\eqalign{y &= (1 - \alpha)x_1 + \alpha x_2 \cr  &= x_1 - \alpha x_1 + \alpha x_2 \cr  &= x_1 + \alpha x_2 - \alpha x_1 \cr  &= x_1 + (x_2 - x_1)\alpha}``)!@ Which then gets translated to the following C code below.In this step, the output is also being scaled by theamplitude =amp=. #+NAME: interpolate_values#+BEGIN_SRC cout = (x1 + (x2 - x1) * fract) * osc->amp;#+END_SRC And now the sample has been computed! To wrap up, theinternal phase amount =lphs= is updated and masked toprevent overflow. #+NAME: update_the_state#+BEGIN_SRC cphs += osc->inc;phs &= SK_OSC_PHASEMASK;osc->lphs = phs;#+END_SRC