The wavetable oscillator code follows the basic structure of most audio processing blocks, with three types of functions:
- at the start: initialization—create the oscillator and its wavetables
- as needed: control—change the oscillator frequency, and pulse width for PWM
- every sample: processing—get next output sample, update phase
At the start
Initialization—Most of our work for the oscillator is here. Fortunately, this is the part that’s not time-critical. We start by initializing the oscillator’s phase to zero, and creating all of the wavetables we need—starting with the lowest frequency wavetable, followed by its progressively bandwidth-reduced copies for successively higher frequency ranges. In this implementation, each wavetable is accompanied by its table length, and the highest pitch that the table is built to handle without noticeable aliasing.
As needed
Control—You may change the oscillator frequency only for new notes, or it might be as often as every sample if you want modulation with updates at the sample rate (also true of pulse-width if you want to pulse-width modulation). And fortunately this task is the simplest and fastest of all for the oscillator, so it’s no problem supporting it for every sample.
Every sample
Processing—In general you need to update each oscillator’s phase for every sample, though you can choose to suspend updating of any oscillator that’s part of a voice that’s silent. And of course the oscillator won’t do you much good if you don’t get its output for every sample and use it.
Updating the phase is trivial—just add the phase increment, and wrap it back into the range of 0-1 if necessary.
Retrieving the current output is a bit more work. First, we need to determine which table to use, by comparing the current frequency control (phase increment) with the highest-frequency table that can accommodate it. Then we use the current phase value, from the phase accumulator, to determine the wavetable positions and weighting for our linear interpolation to yield the oscillator output.
For instance, if the phase accumulator is 0.51, and the length of the appropriate wavetable is 2048, then the offset into the table is 0.51 x 2048, or 1044.48. So, we use table samples at indices 1044 and 1045, and use the fractional part—0.48—in the interpolation.
Key oscillator components
Variables
The oscillator is implemented in C++. We’ll start with the oscillator’s private member variables, but keep in mind that you won’t access these directly—you’ll interact with the oscillator through the member functions.
First, each oscillator needs a phase accumulator and a phase increment (frequency control):
double phasor; // phase accumulator double phaseInc; // phase increment
And since we’re supporting variable pulse width via the difference between two saws, we need a phase offset control:
double phaseOfs; // phase offset for PWM
And each oscillator needs its list of wavetables and how many there are. We pick a maximum number of wavetables to support—this information doesn’t take up much memory, and it’s more efficient than allocating the list dynamically; we can cover 20 Hz to 20 kHz with ten tables—I picked 32 as the maximum in case you want to try half- or third-octave tables, but you can set numWaveTableSlots to whatever value you want:
// list of wavetables int numWaveTables; waveTable waveTables[numWaveTableSlots];
Here’s the wavetable structure—the length of the wavetable, the top frequency it’s designed for, and the wavetable itself (which is allocated dynamically):
typedef struct {
double topFreq;
int waveTableLen;
float *waveTable;
} waveTable;
Why does each wavetable have its own length variable? Because this oscillator doesn’t require a constant size. For example, for a memory-sensitive environment, you could use longer tables for low frequencies in order to accommodate more harmonics, and shorter ones for higher frequency tables. For constant oversampling, you could cut the table size in half (along with the number of harmonics) for each higher octave.
Functions
Now for the oscillator functions. Besides the constructor—which simply initializes the phase accumulator, phase increment, phase offset, and the empty wavetable list—and the destructor, we have this function to initialize the oscillator:
void addWaveTable(int len, float *waveTableIn, double topFreq);
The addWaveTable function specifies the length of a wavetable to add to the oscillator, its contents, and the top frequency (phase increment) it’s designed for. Call it once for each wavetable you need to add—from lowest frequency tables to highest.
We have these functions for setting frequency (normalized to the sample rate—it’s the phase increment) and the phase offset (for setting pulse width, and perhaps other uses in the future):
void setFrequency(double inc); void setPhaseOffset(double offset);
Here is the function to update the oscillator phase by adding adding the phase increment to the phase accumulator—call this once per sample:
void updatePhase(void);
Finally, the routines that retrieve the current oscillator output. Use getOutput, normally. But for variable pulse width, initialize the oscillator for a sawtooth oscillator, and get the oscillator output via getOutputMinusOffset, where the pulse width is set with setPhaseOffset, with a value between 0.0 and 1.0 corresponding to the pulse width:
float getOutput(void); float getOutputMinusOffset(void);
The code
Here’s the wavetable oscillator source code, in a .zip file—WaveTableOsc.cpp and WaveTableOsc.h, along with a file, main.c, that gives a few oscillator examples and generates .wav file output:
More comments and output examples coming soon.
Update: Be sure to check out the updated WaveTableOsc, in WaveTableOsc optimized.
That’s great thank you! It compiles and works. In the function defineSawtooth() would it work if I read the sawtooth from a file, e.g. if the saw was in a WAV as a single-cycle of say 4096 samples. In that case it would be possible to load any kind of waveform into the wavetable.
Yes, that’s right. I’ll have some different kinds of oscillator examples in an “end notes” article coming soon. One of the examples was made from recording myself saying “ah”, then turning an arbitrary cycle of that into an oscillator, 20 Hz to 20 kHz. You need to start with your single cycle waveform (my “ah”, or your 4096-sample sawtooth), and do an FFT of it. Then you can make the next octave up by setting the upper harmonics to zero and doing an inverse FFT, and repeating the process for higher tables.
Thanks for giving the first confirmation of being able to compile and run the source code!
Nigel
That’s great, looking forward to seeing the end notes for the “ah” example as I need some help understanding how to implement that. Thanks.
OK, I was just planning to post it as an audio example, so I’ll give you some rough notes here. I hacked the code out quickly in a few minutes and didn’t bother to keep it, but the steps are pretty simple: First, you need to grab a single cycle (I “ah’d” a frequency that would be close to an even number of samples, so I could find one cycle quickly in an editor; for instance, play an A two octaves below middle-C on the pain, and “ah” that pitch. That would be 110 Hz, or 400.91 samples—I started at a nice looking zero crossing in Audacity and clipped 401 samples, close enough). I wanted the editor to do some work for me, so I resampled that to 2048 samples to be convenient for the FFT, and saved it as a file of raw float samples. I read it into an array in the program, then used the existing FFT in the code to covert it to the frequency domain. From there, it’s pretty much the same as the oscillator code—I used that to create each octave, each time zeroing out the upper half of the harmonics before calling the FFT to convert to the time-domain waveform.
Hi, it’s been a few months but I finally got around to trying to implement this code in my new synth. Please can you clarify the conversion of the a user sample (like the “ah”). I’m loading the .wav then in defineSawtooth() I’ve replaced it with this where rawwav is the float data from the .wav file. Is this what you mean by zeroing the upper harmonics?
void defineWave(float *d, int len, int numHarmonics, float *ar, float *ai)
{
if (numHarmonics > (len >> 1))
numHarmonics = (len >> 1);
float rawwav[len];
// clear
for (int idx = 0; idx < len; idx++)
{
rawwav[idx] = d[idx + 1];
ai[idx] = 0;
ar[idx] = 0;
}
fft(len, rawwav, ar);
for (int idx = 1, jdx = len – 1; idx <= numHarmonics; idx++, jdx–)
{
ar[jdx] = 0;
}
}
The best way to go about this is to first do an FFT to convert your wave to the frequency domain. The simplest case is if your single-cycle is already a power-of-two in length, so let’s say it’s 1024 samples. Obviously, you can use this as your first subtable. But if the wave covers the full bandwidth, the next subtable, used for the next octave up, will alias. So, for the second table, you’d take the frequency domain version of your wave and zero out the upper half of the harmonics, and do an inverse FFT—now you have something similar to your original waveform, except the upper octave of harmonics have been removed.
You repeat this process, adding subtables to the oscillator, until you’ve added a table with just one harmonic—a sine wave for the top octave. For each successful octave subtable, you zero out half of the remaining harmonics in the frequency domain and do an iFFT to get the next wavetable to add.
You zero a harmonic by setting the corresponding real and imaginary parts to zero. Exactly which elements of the array you zero depends on the FFT length and other details. In my example code, I used a plain-vanilla FFT (I could have used a “real” FFT, for instance—more efficient, but more clutter for the reader), so the upper half of the array mirrors the lower half (but remember that first element is DC and Nyquist—the fundamental is at [1] and mirrored at the top end).
Thanks. I think what I have now is right. There was another problem with the selection of the wavetable from the phaseInc which wasn’t quite helping. It would be great if you have an example to compare against though.
I hacked quick code to do that example, and discarded it immediately after—sorry. But I think it might be more helpful to list the basic steps. I started to write a detailed explanation, but it got so long that it would need to be an article on tis own. I can do that, but I need to finish and post some other articles in the queue first.
In a nutshell, you take your wave and use it for the first table. Then you need to remove the upper half of the harmonics for the next table. The best and easiest way to do that is by converting it to the frequency domain with an FFT. Keep a copy of that, and for each subsequent table, set the upper half of the harmonics to zero and do an inverse FFT to get the wavetable data for the next table. Repeat until you’re down the the last harmonic (the fundamental sine wave).
An example of zeroing the harmonics is that for an FFT of 1024 length, the real and imaginary arrays have the indexed order (where “h” is for harmonic) { DC, h1, h2, …, h511, h512 (Nyquist), -h511, …, -h2, -h1 }. That represents 512 harmonics and the mirrored image above Nyquist—you need to zero the harmonics and its image, for both the real and imaginary arrays. So, to zero out harmonics 257-512, you zero indices [257] through [512] and [1024-257] through [1024-512] (yes, this hits 512 twice, but I’m explaining it conceptually, plus it ends up being the simplest way to do it in a single for-loop). The subsequent table would need 129-256 zeroed, additionally, then 64-128, etc. Keep going till only [1] and [1023] remain un-zeroed—that’s your final table.
I’ve posted a new article with source code to create all the necessary wavetables from a single wavetable here: Replicating wavetables
Do not compile this source code in Microsoft Visual C++ 6.
Too many errors in source code is found.
MVC++ 6 was released in 1998. I tend to use more modern compilers (C++11—2011 usually—but I’m not sure what version I used while developing the oscillator code). The modifications for an older version shouldn’t be difficult, but of course you need to understand the error and the language.
Thank you so much for this!
I made a version compatible with VS 2015, available here:
https://github.com/vberthiaume/waveTableOsc
I essentially replaced some arrays with vectors, and in testSawSweep() from main.cpp, I changed this line:
for (int count = sampleRate * 0.05; count >= 0; –count) {
to this
for (int count = sampleRate * 0.05; count > 0; –count) {
[PS:] oh yeah, and I had to do the same >0 fix in testPWM() and testThreeOsc()
Nice, thanks! (btw, this is nerdy, but the “>= 0” was intentional, figuring there was no need to write a zero at the end, since it would effectively be zero the next sample anyway when it stopped playing 😉 )
Hi Nigel,
thanks for this excellent tutorials. I would love to browse your source code instead of downloading. Github is excellent for that, and I saw you have an account. Would you mind to upload your code to a repository (I don’t mind to do it, but I think is better under you user name).
Thanks!
Yes, I’ll do that—might not happen right away, but “soon”. I was planning another article that will include source code, and already decided it was time to put things into a code repository, for the new article and older stuff.
I am just getting started learning about how to make a synthesizer using the JUCE platform, which have a wavetable tutorial that uses on one table for all frequencies…
createWavetable…
setFrequency…
getNextSample…
Finally coming to my question, what is the advantage of using mutiple wavetables for various frequencies instead of just one?
Probably best to look at my video here first. See how the one-table sawtooth peaks bob up and down as the frequency goes higher? That’s because as you have fewer samples per cycle (to get to higher frequencies), the less chance you have of the corner being where it should be in time. A sine wave, in the JUCE example, is smooth, so it doesn’t have this problem. In general, the higher the harmonics, the shaper the waveform details can be, and the more susceptible to aliasing. Or, to look at in in the frequency domain, a sine wave has only the first harmonic, so you can’t get aliasing as long as it’s played back in the audio range. A sawtooth has many harmonics, so the multiple wavetables are just copies with higher harmonics removed to avoid aliasing at higher fundamental frequencies. That is, at 15 kHz, a sine and a saw sound the same, because the second saw harmonic is at 30 kHz, beyond hearing. But at 20 Hz, they sound very different. We need to preserve harmonics for lower frequencies, but delete unnecessary ones at higher frequencies to avoid aliasing.
I have successfully implemented the code in my synth, however, when I try to read the output, the values in the arrays works only for the lowest notes
How can I scale the oscillator to produce values in the -1…1 range instead of the one used in the wav file?
I don’t understand—in particular, the meaning of “works” here.
I mean to say that when I play a low frequency, i can hear the sounds, however for the highest ones it alias
You’re using the code, unchanged? Are you sure you’re adding the tables with the correct ranges?
It does work, but at this point you might want to use the improved code, WaveTabelOsc optimized and WaveUtils updated.
Just got around to looking at your code, thanks a lot.
I am a bit confused about how to use the “void fft(int N, myFloat *ar, myFloat *ai)” function for my wavetable.
My synth have a “Waveform Creator” that can literally make millions of waveforms, although many will sound very alike. Anyways I had hoped I could then “bandlimit” the table after the before mentioned created it, using your FFT function. My table is a 4096 size array of floats.
Be sure you’ve seen the updated oscillator and WaveUtils. The new WaveUtils has more flexibility and examples. WaveTableOsc optimized is needed because naming conventions changed from the earlier version, and WaveUtils.
Nice article and site
the Wavetable oscillator and waveform examples all work great but one thing is realy bothering me. The “defineSawtooth” function (and also the commented out square and triangle options within it) seems to define the waveform using its frequency content. This content is defined for the “ar” array with its low indices containing what i would expect (harmonics decreasing in magnitude). However the last par of the “ar” array is filled with a mirror image but uses negative values. Why are negative values used here? I dont think I have ever seen a frequency plot that looks this way.
Secondly, once the ar array is filled, the “makeWaveTable” function applies an FFT via
fft(len, ar, ai);
so you are applying the FFT to a frequency spectrum? why is that? I would have expected an IFFT to be used which is very confusing for me. Is this a known method that works due to the “mirror” in the array I mentioned above? Would be great if there is a reference or term for this.
lastly, I assume that “ai” and “ar” refer to IM and RE parts, where the IM parts in the “define” functions are always set to zero.
Thanks for the code and article. Apologies if you have explained this in another post but I couldn’t find anything.
Good questions. First, an amazing and initially counter-intuitive property of the Fourier Transform is that it translates both ways. It’s a little more intuitive if you think about it in terms of my explanation in A gentle introduction to the FFT. In it, sin and cos are used to modulate the time domain samples to yield the frequency domain. But in the other direction, sin and cos modulate the frequency domain “bins”—yielding the time domain signal.
The short story is that the important difference in the two is scaling. There are different scaling conventions (by N, by the square root), and none of it matter if you’d already scaling (such as scaling the waveform to fit a certain amplitude).
As for the other question…recognize that a sawtooth has a progression of harmonic amplitudes, but also phases. In other words, there is no difference in harmonic amplitudes between a ramp (sloping up) and sawtooth (sloping down). And they sound the same, but the shape can be important (certainly for an LFO). I made a sawtooth.
Thanks for the reply. Its appreciated.
Just to be extra clear about what is confusing me, I tried commenting out the following line from the “defineSawtooth” function:
// sawtooth
for (int idx = 1, jdx = len – 1; idx <= numHarmonics; idx++, jdx–) {
myFloat temp = -1.0 / idx;
//ar[idx] = -temp;
ar[jdx] = temp;
}
I then ran the program and fed the signal through a scope. To be honest I cant see or hear any difference. So what is the purpose of using both "idx" and "jdx"?
Incidentally, commenting out the positive "temp" fill also seemed to have no effect on the output.
All the best
First, to be clear, that sort of stuff is just quick & dirty examples of how to use the oscillator and utilities. I could have used a real-only FFT and simplified things, but I didn’t want to add another layer of confusion. idx and jdx count from opposite ends—I could have used a single index variable (idx) and some math to count it down from the top. If you fill just half the array (comment out that line), the result is half amplitude. But when the values are auto-scaled, it doesn’t matter, the final output will be the same.