Perspective on dither

Recently, I’ve had lengthy discussions on the topic of dither with a couple of different people—of opposite views. One believes that everything should be dithered, including truncations to 24-bit. The other feels that dither is a waste of time even for 16-bit, and is needed only for shorter word lengths, rarely used today.

This got me thinking about how to give perspective on the truncation distortion and dither levels we’re talking about for 16-bit and 24-bit files. The usual demonstration would be to find or create a recording that is not unrealistic, yet produces noticeable truncation distortion at the bit levels we’re interested in—which is mainly near the floor of 16-bit and 24-bit samples.

But it occurred to me that people don’t have a good idea of how loud the distortion levels are by themselves, and that should be the place to start. I thought of a way to generate signals of those levels with no distortion of their own.

Perfect audio

Are you ready to listen to perfect digital audio files? By perfect I mean that they are exactly what they are designed to be. If I’d tried for low-level sine waves, for instance, they would be crude, distorted approximations. Instead, I created square waves at precise bit levels to avoid quantization errors, synchronized to the sample rate to avoid aliasing. The sample values are exact—they are not values that are “close” to expected values—and the harmonic content is representative of truncation effects at those levels.

The signal begins with the LSB set (we’ll call this “+1”), then the next sample is negated (-1). This repeats (+1, -1, +1, -1, …) for a short duration, then the period increases by one sample (+1, +1, -1, -1, +1, +1, -1, -1, …) for a short duration, and the process continues, increasing the period by one sample each time. This is a classic “divide by n” oscillator. You’ll notice that the pitch resolution is very poor at the beginning, as the second tone is an octave down from the first (which is at Nyquist and will be swallowed by your reconstruction filter), but gets better and better as pitch drops and each additional sample is a small percentage of the period. The waveform amplitude is not the smallest possible—toggling between +1 and 0 would do that—but is representative of the amplitude of the smallest truncation effects for bipolar signals.

So, the sound starts dropping at discrete frequencies, and moves increasingly towards a smooth frequency sweep.

Listening tests

Listen to these files on a quality monitoring system. First, the generated test signal toggling at the level of the fifth bit, moderately loud, so that you can become familiar with what you’ll be listening for at the lower levels:

Sweep at the 5th bit level (-24.1 dB)

All files are 24-bit. The difference is that the +1 level for the 24-bit version is the 24th bit, for the 16-bit version it’s the 16th bit, and for the 5-bit version it’s the 5th bit (making it 2048 times, or 66 dB, louder than the 16-bit version).

Now for 16-bit. Start by setting your monitoring level to as loud as you would normally play music, by running a song through it. Crank it up, but don’t hurt your ears. Then play this file by itself:

Sweep at the 16th bit level (-90.3 dB)

On a quality system at a relatively high monitoring level, you’ll have no problem hearing this sweep. After listening a few times, you might want to play a song again to remind yourself of the relative level of the test signal. Think about how difficult it would be to hear the test signal during the chunking guitars of a rock tune, but how it might be heard in the fading of quiet piano notes in classical music. (Yes, truncation noise can resulting in more of a tearing sound that the tone of the test sweep, but the relative levels are still representative of the relationship.)

Now, 24-bit. Warning: Do NOT try to compensate by cranking your volume level. You are NOT going to hear this signal anyway, and you risk making a mistake and assaulting your ears horrifically, or blowing speakers. No matter how great you think your 24-bit converters are, it and all of your other gear generate higher levels of noise than the 24th-bit level.

Sweep at the 24th bit level (-138.5 dB)

If you want to try other levels, use a sample editor (such as the free Audacity, or your favorite DAW) to boost the level 2x, or 6.02 dB, for every additional bit that you want to move up.

Think about the levels, and at what bit levels you’d care to conceal truncation distortion with dither, and I’ll comment further in a future article.

Posted in Dither | 4 Comments

Replicating wavetables

My WaveTableOsc object uses a series of wavetables that are copies of single-cycle waveforms pre-rendered at a progressively smaller bandwidth. By “replicating wavetables” I mean taking a full-bandwidth single cycle waveform, and using it to make a complete set of progressively bandwidth-reduced wavetables. In this article, I’ll demonstrate principles and source code that will let you take any source waveform suitable for the lowest octave, and build the rest of the tables automatically.

The source waveform might be any time-domain single-cycle wave—something you drew, created with splines, recorded with a microphone, calculated—or it can be a frequency-domain description—harmonics and phases. Fundamentally, we’ll use a frequency-domain description—the array pair we use in an FFT—that holds amplitude and phase information for all possible harmonics. If we start with a time-domain wave instead, we simple convert it to frequency domain as the first step. The source code builds the bandwidth-reduced copies, converts them to the time domain, and adds them to the oscillator until it contains all necessary tables to cover the audio range without noticeable aliasing.

We’ll build a new table for each octave higher, by removing the upper half of the previous table’s harmonics. If it’s not clear why, consider that all harmonics get doubled in frequency when shifting pitch. A 100 Hz sawtooth wave has harmonics spaced 100 Hz apart; shifting it up an octave to 200 Hz shifts the spacing to 200 Hz apart, so half of the harmonics that were below Nyquist move above.

The code

In principle, the code is simple:

repeat while the wavetable has harmonic content:
	add wavetable to oscillator
	remove upper half of the harmonics

This is one of those places where I’d like to be more sophisticated, and use a faster “real” FFT and avoid the silliness of filling in a mirror-image array, but it’s more code and clutter—I’m sticking with the generic FFT from the WaveTableOsc source code. I’ve also reused the makeWaveTable function, so the new code here is minimal. The end product is WaveUtils.h and WaveUtils.cpp, plain C utility functions. I’d probably wrap the functions into a manager class that dealt with WaveTableOsc and wavetable classes, if I were focusing on a specific implementation and not an instructive article.

The two reused functions are used internally, and we need just one new function that we call directly—fillTables:

void fillTables(WaveTableOsc *osc, double *freqWaveRe, double *freqWaveIm, int numSamples) {
    int idx;
    
    // zero DC offset and Nyquist
    freqWaveRe[0] = freqWaveIm[0] = 0.0;
    freqWaveRe[numSamples >> 1] = freqWaveIm[numSamples >> 1] = 0.0;
    
    // determine maxHarmonic, the highest non-zero harmonic in the wave
    int maxHarmonic = numSamples >> 1;
    const double minVal = 0.000001; // -120 dB
    while ((fabs(freqWaveRe[maxHarmonic]) + fabs(freqWaveIm[maxHarmonic]) < minVal)
        && maxHarmonic) --maxHarmonic;

    // calculate topFreq for the initial wavetable
    // maximum non-aliasing playback rate is 1 / (2 * maxHarmonic), but we allow
    // aliasing up to the point where the aliased harmonic would meet the next
    // octave table, which is an additional 1/3
    double topFreq = 2.0 / 3.0 / maxHarmonic;
    
    // for subsquent tables, double topFreq and remove upper half of harmonics
    double *ar = new double [numSamples];
    double *ai = new double [numSamples];
    double scale = 0.0;
    while (maxHarmonic) {
        // fill the table in with the needed harmonics
        for (idx = 0; idx < numSamples; idx++)
            ar[idx] = ai[idx] = 0.0;
        for (idx = 1; idx <= maxHarmonic; idx++) {
            ar[idx] = freqWaveRe[idx];
            ai[idx] = freqWaveIm[idx];
            ar[numSamples - idx] = freqWaveRe[numSamples - idx];
            ai[numSamples - idx] = freqWaveIm[numSamples - idx];
        }
        
        // make the wavetable
        scale = makeWaveTable(osc, numSamples, ar, ai, scale, topFreq);

        // prepare for next table
        topFreq *= 2;
        maxHarmonic >>= 1;
    }
}

Defining an oscillator in the frequency domain

Here’s an example that creates and returns a sawtooth oscillator by specifying the frequency spectrum (note that we must include the mirror of the spectrum, since we’re using a complex FFT):

WaveTableOsc *sawOsc(void) {
    int tableLen = 2048;    // to give full bandwidth from 20 Hz
    int idx;
    double *freqWaveRe = new double [tableLen];
    double *freqWaveIm = new double [tableLen];
    
    // make a sawtooth
    for (idx = 0; idx < tableLen; idx++) {
        freqWaveIm[idx] = 0.0;
    }
    freqWaveRe[0] = freqWaveRe[tableLen >> 1] = 0.0;
    for (idx = 1; idx < (tableLen >> 1); idx++) {
        freqWaveRe[idx] = 1.0 / idx;                    // sawtooth spectrum
        freqWaveRe[tableLen - idx] = -freqWaveRe[idx];  // mirror
    }
    
    // build a wavetable oscillator
    WaveTableOsc *osc = new WaveTableOsc();
    fillTables(osc, freqWaveRe, freqWaveIm, tableLen);

    return osc;
}

Defining an oscillator in the time domain

I you have a cycle of a waveform in the time domain, you can create and oscillator this way—just do an FFT to convert to the frequency domain, and pass the result to fillTable to complete the oscillator:

WaveTableOsc *waveOsc(double *waveSamples, int tableLen) {
    int idx;
    double *freqWaveRe = new double [tableLen];
    double *freqWaveIm = new double [tableLen];
    
    // convert to frequency domain
    for (idx = 0; idx < tableLen; idx++) {
        freqWaveIm[idx] = waveSamples[idx];
        freqWaveRe[idx] = 0.0;
    }
    fft(tableLen, freqWaveRe, freqWaveIm);

    // build a wavetable oscillator
    WaveTableOsc *osc = new WaveTableOsc();
    fillTables(osc, freqWaveRe, freqWaveIm, tableLen);
    
    return osc;
}

Here’s a zip file containing the source code, including the two examples:

WaveUtils source code

Caveat

As with our wavetable oscillator, this source code requires wavetable arrays—whether time domain or frequency domain—of a length that is a power of 2. As examined in the WaveTableOsc articles, 2048 is an excellent choice. The source code does no checking to ensure powers of 2—it’s up to you.

Posted in Oscillators, Source Code, Wavetable Oscillators | 43 Comments

About source code examples

My goal is to teach audio DSP principles in a way that is more intuitive than most available material. And part of that goal is to help you to think about the your goals and how to solve them by showing you my thought process.

So at one extreme, source code shouldn’t be part of what I’m doing here—ideas over implementations. However, examples are a powerful learning tool, and I want to show that the principles introduced here are practical. But the other extreme is to develop a full open-source library here, and that’s something different entirely. Foremost, it conflicts with what I’m trying to do here—even if I had time to do both. My goal with source code is to give you minimal, but commercial-grade examples in source code. Something lean enough for you to dissect and see the principles applied.

For example, the WaveTableOscillator object is a high-quality DSP component—good stuff. But it’s not exactly how I would write it for myself for a major project. In addition to the WaveTableOscillator class, which doesn’t create the wavetables themselves, I might write a wavetable class that does, and a manager class that would use the two classes, keeping track of oscillators and sharing common wavetables between them to save space and time. And I might use C++ templates. And asserts for trapping bad parameters at development time. But in the context of this blog, these things tend to obscure the functionality that I’m trying to focus on.

I’m getting ready to lay more source code on you in coming articles. I think that working C++ code, that you can experiment with, is better than using pseudo code, and I’ll continue in that direction for now.

Posted in Source Code | 2 Comments

Ventilator adapter in a mint tin

I don’t eat mints often, but when I do, they are sugar free…

A Leslie speaker simulation

Players of Hammond organs and their clones know how important the Leslie speaker is to the classic sound. Modern clones have electronic simulations built-in, with varying degrees of success. The Neo Instruments Ventilator is a particularly convincing DSP-based Leslie simulation in a foot pedal. I was never satisfied with the Leslie simulation in my Korg CX-3 (the modern DSP variety—they also made an electronic-component-based model of the same name circa 1980). Bypassing the internal Leslie simulation in favor of the Ventilator greatly improved my enjoyment of the instrument. (I have a Leslie 122, but it’s in a state of disrepair—besides, an electronic simulation makes it much quicker to record a sudden inspiration.)

But something’s missing

But the Ventilator has a serious omission: no MIDI control. A great advantage of Hammond clones over the real thing is that you can capture your performance via MIDI, then audition alternate drawbar settings and try different Leslie speeds later, in playback, without the need to repeat the performance. Without MIDI control, the Ventilator has no way to capture and play back the simulation of switching Leslie speeds.

To fix this shortcoming, I bought a MIDI-controlled relay (made by MIDI Solutions), as the Ventilator has an external footswitch input. But when I went to integrate it, I found what I deem a design failure in the Ventilator.

The (optional) external controller for the Ventilator has two footswitches, one to control fast or slow modulation, and the other for “brake” (stop). It connects with the Ventilator via a tip-ring-sleeve plug. I don’t use the brake function—I only want to switch fast/slow. but the way the Ventilator is designed, a single relay can’t do the job. To switch to fast speed, one pin goes high and the other low, and to go slow, the two pins switch to the opposites states. This would require two single relays configured to switch in opposite directions.

These MIDI relays aren’t cheap, so the thought of buying another was not a desirable choice. And it was irksome that Neo Instruments could have easily designed the Ventilator to allow a single relay to do the job with no change of other functionality. This would also allow a user to locate the Ventilator on top of the organ, for access to the knobs during performance, and use a simple and cheap single footswitch to switch speed instead of their relatively expensive external dual switch.

The truth

Here’s how they could have have done it (it’s a matter of software, so they could still do it): The tip and ring are each “pulled up” to +5 volts, individually, and are sensed by the processor so that each yield a logical “1” when disconnected, and a logical “0” when connected to ground (the sleeve). This give us two bits, or four possible states. Three states are used—brake, fast, and slow. Here’s the truth table (“0” indicates that the pin is connected to the sleeve—common—and “1” indicates no connection—open):

Tip Ring State
0 0 (undefined/unused)
0 1 fast
1 0 slow
1 1 brake

Note that there is one left over state—an undefined state (not possible to attain with their dual footswitch). However, if the Ventilator had been designed to use this state as a redundant speed selector, the device would work with either a single switch (SPST—single-pole single-throw), or their dual footswitch with no change.

The fix

OK, how do I fix this slight, and use my MIDI relay for speed switching? One choice would be to modify the MIDI relay. I popped it open, and found that the output jack is already a TRS connector, so I had all the pins I needed. I could try to find a small dual relay (SPDT—single-pole, double-throw), but even if I could, it’s possible there would be switching problems if I didn’t have enough current available (the MIDI relay box is powered from the MIDI current loop). While this would let the MIDI relay switch speeds, it wouldn’t let me use a simple single footswitch as an alternative, and I’d like that option. Alternatively, I could build a mod into the Ventilator—this is probably the easiest and most flexible choice, since I’d have access to the Ventilator’s power supply. After some thought, I decided that I’d prefer building an external adapter instead, leaving the Ventilator unaltered.

The biggest issue of an external adapter is that I didn’t want to use a battery. Since the Ventilator pulls the tip and ring pins up to +5 volts through resistors, internally, I needed a solution that tapped that current for power. Fortunately, the solution requires only a single transistor and resistor. When the input switch is closed, the NPN transistor switch opens, and when the switch opens, the NPN turns on and closes the connection.

Great, a small, simple circuit. But not small enough to fit in the shell of a connector, so I needed to put it in a box…

Building it

About the mints…Eclipse mints are sugar-free and come in these cool little metal boxes with a hinged, snap-shut lid. I’ve kept them instead of tossing them in the recycle bin—it seemed they could be handy for something. The only use I’d found so far is as a temporary holder of screws when I’m fixing gear, but I took a close look at one for this project. Yes, the dimensions are ideal for two 1/4″ phone jacks, with enough room left over for my circuit board and wiring inside.

Here’ s a circuit diagram—click to see a larger view, or go here to run a circuit simulation:

Adapter schematic

One transistor, one resistor, and a few jumpers on a prototyping board is all it takes. Trim the circuit board—easily under a half square inch. Drill a couple of 3/8″ holes for the jacks, and feed them through the open end with needle nose pliers. The metal’s thin, but once you tighten the nuts everything firms up again. Cost: $1.59 for the prototyping board, 69 cents per jack. I had the transistor (you can substitute other common small-signal NPN transistors) and resistor on hand, but you can find them for about $1.50 for packs of four or five.

Here’s a view before securing the board to the inside of the box with adhesive-backed foam tape:

Adapter, open view

The input connected to the relay or footswitch with a mono cable, and the output to the Ventilator with a stereo cable. I installed the output jack next to the twin mints of the box graphics as a reminder.

My CX-3 accepts a simple momentary footswitch to toggle Leslie simulation speeds, in turn outputting a MIDI control change that can toggle the MIDI relay. Since the MIDI relay box has a MIDI output that passes the MIDI notes and control changes, I can record notes and speed switching for playback later.

Posted in Uncategorized | Leave a comment

A one-pole filter

Here’s a very simple workhorse of DSP applications—the one-pole filter. By comparison, biquads implement two zeros and two poles. You can see that our one-pole simply discards the zeros (the feed-forward delay paths) and the second pole (feedback path):

We keep the input gain control, a0, in order to compensate for the gain of the feedback path due to b1, which sets the corner frequency of the filter. Specifically, b1 = -e-2πFc and a0 = 1 – |b1|. (In the implementation, we’ll roll the minus sign in the summation into b1 and make it positive, for convenience.)

Filter types

A pole is a single-frequency point of pushing gain up, and a zero is a single-frequency point of pulling the gain down; with a single pole, you are not going to get complex response curves such as bandpass, peak, and shelving filters that you can get with the two poles and zeros of a biquad. That leave us with lowpass and highpass, which have a single point of change.

However, a one-pole makes a poor highpass filter for cases in which we might be likely to use it—in particular, a DC blocker. That’s because it makes a highpass by pushing response up at high frequencies—we really need a zero to pull the response down at very low frequencies. So, we’ll only implement coefficient calculation for lowpass response here. However, we can still make a highpass filter suitable for a DC blocker by subtracting the output of a one-pole lowpass filter, set to a low frequency, from the direct signal.

Typical uses

What are typical uses of this one-pole lowpass filter? Sometimes, we do need the gentle slope of a first order frequency roll-off. And often, we use this simple filter on parameters instead of one the audio directly—some gentle smoothing of a calculated signal in an automatic gain control, for instance.

The 6 dB per octave slope of the one-pole filter—a halving of amplitude for every doubling of frequency—is gentle and natural. It’s extremely cheap, computationally, so it’s the perfect choice when you need just a bit of smoothing. And it doesn’t “ring” (overshoot), so it’s an excellent choice for filtering knob changes. Run each of your user interface sliders and knobs through a one-pole lowpass set to a low (sub-audio) frequency, and your glitchy, zippered control changes turn into smooth parameter changes.

Note that if you feed the one-pole lowpass an impulse, it will yield a perfect exponential decay—the same decay that we use on typical ADSR envelope generators on synthesizers. To look at it another way, the feedback path is an iterative solution to calculating an exponential curve.

The source code

The one-pole filter is so simple that we’ll make it inline, entirely in the header file:

//
//  OnePole.h
//

#ifndef OnePole_h
#define OnePole_h

#include <math.h>

class OnePole {
public:
    OnePole() {a0 = 1.0; b1 = 0.0; z1 = 0.0;};
    OnePole(double Fc) {z1 = 0.0; setFc(Fc);};
    ~OnePole();
    void setFc(double Fc);
    float process(float in);
    
protected:    
    double a0, b1, z1;
};

inline void OnePole::setFc(double Fc) {
    b1 = exp(-2.0 * M_PI * Fc);
    a0 = 1.0 - b1;
}

inline float OnePole::process(float in) {
    return z1 = in * a0 + z1 * b1;
}

#endif

Examples

Here’s a DC blocker:

OnePole *dcBlockerLp = new OnePole(10.0 / sampleRate);
...
// for each sample:
sample -= dcBlockerLp->process(sample);

Need highpass?

It’s trivial to add one-pole highpass response, if you’d like—just negate the lowpass calculation; this inverts the spectrum, so you’ll need to flip the frequency range by subtracting Fc from 0.5:

b1 = -exp(-2.0 * M_PI * (0.5 - Fc));
a0 = 1.0 + b1;

Pole position

You can experiment with how the location of the pole affects frequency response by using the pole-zero placement app. Click the Pair checkbox near the Poles sliders to disable it. Move both Poles sliders to center to move the poles to the origin. (The zeros are already at the origin by default when you open the page. At the origin, poles and zeros have no effect.) Now, move just one pole slider; you’ll see highpass response when the pole is to the left of the origin, and lowpass to the right.

Posted in DC Blocker, Digital Audio, IIR Filters, Source Code | 19 Comments

A note about de-normalization

One ugly thing that we need to be aware of, especially for filters, is the de-normalization of numbers. Basically, computer processors try to keep floating point numbers normalized—they try to keep them in the binary form 1.xxxxxx (where each x is 0 or 1) times 2 raised to a positive or negative power. When numbers become extremely small—smaller than can be expressed by the largest negative exponent—the processor tries to maintain precision by allowing the number to become de-normalized—to slip to 0.xxxxxx times 2 raised to the lowest exponent, in order to allow smaller but less precise numbers. The catch is that most processors become much slower when performing computations with de-normalized numbers (“denormals”).

Denormals happen when you multiply a small number by a tiny number. Typically, denormals are created in filters that ring out with no new input. Almost any new input, including background noise when recording, will keep denormals from appearing. However, if you process some sound, then follow it with samples of zero, a filter or reverb can decay to denormals; the denormals are then passed to the next processing block and cause it to slow as well.

Some modern processors have specific instructions to treat denormals as zero, solving our problem—if your code enforces that setting. You can repair a denormal by checking flags in a floating point number and substituting zero. But that’s a lot of work. A simple way to avoid denormals is to add a tiny but normalized number to calculations that are at risk—1E-18, for instance is -200 dB down from unity, and won’t affect audio paths, but when added to denormals will result in 1E-18; this won’t affect audible samples at all, due to the way floating point works.

The Pentium 4 in particular becomes devastatingly slow—so slow that it can’t keep up with real time audio processing. For audio purposes, denormals might as well be zero. They are so many dB down from unity that their only effect is the potential of bogging down our processing.

While denormal protection could be built into the filter code, it’s more economical to use denormal protection on a more global level, instead of in each module. As an example of one possible solution, you could add a tiny constant (1E-18) to every incoming sample at the beginning of your processing chain, switching the sign of the constant (to -1E-18) after every buffer you process to avoid the possibility of a highpass filter in the chain removing the constant offset.

The best solution depends on your exact circumstances, so you should make a point of understanding the threat. For instance, in the biquad code, there is a potential to generate denormals in the calculations of z1 and z2; if the filter input goes to zero, the output alone recirculates through z1 and z2, potentially getting closer and closer to zero until the values become de-normalized and continue to recirculate in slower de-normalized calculations, as well as passing to the next processing block and slowing its calculations.

Lest you think that you need only be concerned about denormals when fine-tuning your code, take a look at the performance hit for various processors (special thanks to RCL for permission to reproduce his chart):

Yes, that’s more than 100 times slower execution of test code on a Pentium 4 when the calculations become denormalized. I hit this one myself, years ago: After working on DSP chips, primarily, which are unaffected by this issue, I was developing a commercial product for Mac and PC that ran on the host CPU. I did the initial development on a Power PC processor, then gave it a run on my PC—a P4. Audio processing came to a dead stop. Fortunately, I was aware of the issue with denormals, and went straight to a solution that fixed it.

Please read RCL’s excellent article on denormals, “Denormal floats across architectures“, in which he goes into more detail, including comparisons of various CPU architectures.

Posted in Digital Audio | 1 Comment

Biquad C++ source code

I don’t want to get into the business of teaching people how to code—there are a huge number of free resources available on the internet to that do that. But I’ll give a small taste for those trying to get started.

First, I want to be clear that biquads, especially the direct forms, get a lot of coverage simply because they are a straight-forward and basic building block—not because they are always the best choice in filters. The direct forms are terrible choices if you want a sweepable synthesizer filter, for instance. But they are often an adequate and computationally efficient choice for many types of audio and control signal filtering. And they are simple and easy to understand. And since beginners need a place to start, they serve as a good example for coding a DSP example.

C++ is a popular language for DSP coding, and coding in general. It’s object-oriented, and that makes for a good model in developing our audio processing building blocks. Here’s how we approach a biquad filter object.

Filter object requirements

Our filter object should carry it’s own settings and states. That way, we can have multiple instances of the filter, each remaining independent, with its own filter type, frequency, and internal state. We define the state and settings variables (.h) file, giving them a name and reserving space for them.

For this filter, we’ll have a filter type (lowpass, highpass, bandpass, notch, peak, lowshelf, and highshelf), frequency, Q, and peak gain (used only for the peak and shelf types). I’ll use normalized frequency, where 1.0 is the sample rate—divide the frequency you want to set, in Hertz or cycles per second, by your sample rate to get normalized frequency (2438.0/44100 for 2438 Hz at a sample rate of 44100 Hz). We’ll set peak gain in dB, where negative numbers are for cut, and positive numbers for boost.

Biquads need to convert the frequency, Q, and peak gain parameters to the a0, a1, a2, b1, and b2 factors. Since that takes a fair amount of computation, it makes sense to calculate the factors only when the parameters change, which is almost always much lower than the sample rate at which the filter processes audio. So, we’ll want to reserve space to hold those factors. And because we want to be able to set the parameters individually, and because they rely on each other, we’ll keep copies of the current parameter settings as well.

For this example, we’ll use the transposed direct form II architecture (which has good floating point characteristics, and requires only two unit delays), so we’ll need to reserve two memory locations for the delay elements.

To specify the filter type, we assign a number for each type, using enumerated types so that we can attach a more meaningful name to the type. Here’s a look at that definition:

enum {
    bq_type_lowpass = 0,
    bq_type_highpass,
...
    bq_type_highshelf
};

And here are the variables we’re allocating for each filter instance:

    int type;
    double a0, a1, a2, b1, b2;
    double Fc, Q, peakGain;
    double z1, z2;

Filter functions

Next, we need to think about the software functions. We’ll need functions to set the type, frequency, Q, and peak gain individually. For convenience, we’ll also include a function that accepts all of the parameters at once, particularly handy for setting fixed filters. Also, we always need a constructor, called to create a filter instance, and a destructor, called (sometimes implicitly) to remove a filter instance and free its memory. For convenience we’ll include a bare-bones constructor, and a second one that includes setting the filter—especially handy for filters that need to be set just once.

And we’ll need to function that processes and input sample and returned the filtered result—this is called for each sample, and we want it to be as efficient as possible. If our processing routine is small, we’d like to make it an inline function, which saves some overhead of the function call by letting the compiler place the code contained in the function inline. To do that, we define the function as inline, and place the code in the header file. Since all all of these functions need to be accessible by other code, they need to be defined as public; since all of the variables are internal to the filter, we define them as protected—external code has no access. (We make the code a little more tidy by splitting off the computation biquad factors into a separate function; since it’s only accessed internally, we make it protected.)

We’ll use the double type (64-bit floating point) for our computations, and float (32-bit floating point) as the input and output sample size. (We could define 64-bit samples, or use templates to accommodate different sample types, but that’s overdoing it for this example.) Here’s our completed header file (with a common technique to avoid nested includes from being a problem for the compiler):

//
//  Biquad.h
//
//  Created by Nigel Redmon on 11/24/12
//  EarLevel Engineering: earlevel.com
//  Copyright 2012 Nigel Redmon
//
//  For a complete explanation of the Biquad code:
//  http://www.earlevel.com/main/2012/11/26/biquad-c-source-code/
//
//  License:
//
//  This source code is provided as is, without warranty.
//  You may copy and distribute verbatim copies of this document.
//  You may modify and use this source code to create binary code
//  for your own purposes, free or commercial.
//

#ifndef Biquad_h
#define Biquad_h

enum {
    bq_type_lowpass = 0,
    bq_type_highpass,
    bq_type_bandpass,
    bq_type_notch,
    bq_type_peak,
    bq_type_lowshelf,
    bq_type_highshelf
};

class Biquad {
public:
    Biquad();
    Biquad(int type, double Fc, double Q, double peakGainDB);
    ~Biquad();
    void setType(int type);
    void setQ(double Q);
    void setFc(double Fc);
    void setPeakGain(double peakGainDB);
    void setBiquad(int type, double Fc, double Q, double peakGain);
    float process(float in);
    
protected:
    void calcBiquad(void);

    int type;
    double a0, a1, a2, b1, b2;
    double Fc, Q, peakGain;
    double z1, z2;
};

inline float Biquad::process(float in) {
    double out = in * a0 + z1;
    z1 = in * a1 + z2 - b1 * out;
    z2 = in * a2 - b2 * out;
    return out;
}

#endif // Biquad_h

Now, for the implementation of the rest of our functions in the .cpp file:

//
//  Biquad.cpp
//
//  Created by Nigel Redmon on 11/24/12
//  EarLevel Engineering: earlevel.com
//  Copyright 2012 Nigel Redmon
//
//  For a complete explanation of the Biquad code:
//  http://www.earlevel.com/main/2012/11/26/biquad-c-source-code/
//
//  License:
//
//  This source code is provided as is, without warranty.
//  You may copy and distribute verbatim copies of this document.
//  You may modify and use this source code to create binary code
//  for your own purposes, free or commercial.
//

#include <math.h>
#include "Biquad.h"

Biquad::Biquad() {
    type = bq_type_lowpass;
    a0 = 1.0;
    a1 = a2 = b1 = b2 = 0.0;
    Fc = 0.50;
    Q = 0.707;
    peakGain = 0.0;
    z1 = z2 = 0.0;
}

Biquad::Biquad(int type, double Fc, double Q, double peakGainDB) {
    setBiquad(type, Fc, Q, peakGainDB);
    z1 = z2 = 0.0;
}

Biquad::~Biquad() {
}

void Biquad::setType(int type) {
    this->type = type;
    calcBiquad();
}

void Biquad::setQ(double Q) {
    this->Q = Q;
    calcBiquad();
}

void Biquad::setFc(double Fc) {
    this->Fc = Fc;
    calcBiquad();
}

void Biquad::setPeakGain(double peakGainDB) {
    this->peakGain = peakGainDB;
    calcBiquad();
}
    
void Biquad::setBiquad(int type, double Fc, double Q, double peakGainDB) {
    this->type = type;
    this->Q = Q;
    this->Fc = Fc;
    setPeakGain(peakGainDB);
}

void Biquad::calcBiquad(void) {
    double norm;
    double V = pow(10, fabs(peakGain) / 20.0);
    double K = tan(M_PI * Fc);
    switch (this->type) {
        case bq_type_lowpass:
            norm = 1 / (1 + K / Q + K * K);
            a0 = K * K * norm;
            a1 = 2 * a0;
            a2 = a0;
            b1 = 2 * (K * K - 1) * norm;
            b2 = (1 - K / Q + K * K) * norm;
            break;
            
        case bq_type_highpass:
            norm = 1 / (1 + K / Q + K * K);
            a0 = 1 * norm;
            a1 = -2 * a0;
            a2 = a0;
            b1 = 2 * (K * K - 1) * norm;
            b2 = (1 - K / Q + K * K) * norm;
            break;
            
        case bq_type_bandpass:
            norm = 1 / (1 + K / Q + K * K);
            a0 = K / Q * norm;
            a1 = 0;
            a2 = -a0;
            b1 = 2 * (K * K - 1) * norm;
            b2 = (1 - K / Q + K * K) * norm;
            break;
            
        case bq_type_notch:
            norm = 1 / (1 + K / Q + K * K);
            a0 = (1 + K * K) * norm;
            a1 = 2 * (K * K - 1) * norm;
            a2 = a0;
            b1 = a1;
            b2 = (1 - K / Q + K * K) * norm;
            break;
            
        case bq_type_peak:
            if (peakGain >= 0) {    // boost
                norm = 1 / (1 + 1/Q * K + K * K);
                a0 = (1 + V/Q * K + K * K) * norm;
                a1 = 2 * (K * K - 1) * norm;
                a2 = (1 - V/Q * K + K * K) * norm;
                b1 = a1;
                b2 = (1 - 1/Q * K + K * K) * norm;
            }
            else {    // cut
                norm = 1 / (1 + V/Q * K + K * K);
                a0 = (1 + 1/Q * K + K * K) * norm;
                a1 = 2 * (K * K - 1) * norm;
                a2 = (1 - 1/Q * K + K * K) * norm;
                b1 = a1;
                b2 = (1 - V/Q * K + K * K) * norm;
            }
            break;
        case bq_type_lowshelf:
            if (peakGain >= 0) {    // boost
                norm = 1 / (1 + sqrt(2) * K + K * K);
                a0 = (1 + sqrt(2*V) * K + V * K * K) * norm;
                a1 = 2 * (V * K * K - 1) * norm;
                a2 = (1 - sqrt(2*V) * K + V * K * K) * norm;
                b1 = 2 * (K * K - 1) * norm;
                b2 = (1 - sqrt(2) * K + K * K) * norm;
            }
            else {    // cut
                norm = 1 / (1 + sqrt(2*V) * K + V * K * K);
                a0 = (1 + sqrt(2) * K + K * K) * norm;
                a1 = 2 * (K * K - 1) * norm;
                a2 = (1 - sqrt(2) * K + K * K) * norm;
                b1 = 2 * (V * K * K - 1) * norm;
                b2 = (1 - sqrt(2*V) * K + V * K * K) * norm;
            }
            break;
        case bq_type_highshelf:
            if (peakGain >= 0) {    // boost
                norm = 1 / (1 + sqrt(2) * K + K * K);
                a0 = (V + sqrt(2*V) * K + K * K) * norm;
                a1 = 2 * (K * K - V) * norm;
                a2 = (V - sqrt(2*V) * K + K * K) * norm;
                b1 = 2 * (K * K - 1) * norm;
                b2 = (1 - sqrt(2) * K + K * K) * norm;
            }
            else {    // cut
                norm = 1 / (V + sqrt(2*V) * K + K * K);
                a0 = (1 + sqrt(2) * K + K * K) * norm;
                a1 = 2 * (K * K - 1) * norm;
                a2 = (1 - sqrt(2) * K + K * K) * norm;
                b1 = 2 * (K * K - V) * norm;
                b2 = (V - sqrt(2*V) * K + K * K) * norm;
            }
            break;
    }
    
    return;
}

Download the Biquad source code.

Usage examples

Job done. Now how to we use it? When we need a biquad, we just instantiate the Biquad object, request its initial setting, then call the process function for every sample we want to filter.

Biquad *lpFilter = new Biquad();	// create a Biquad, lpFilter;

lpFilter->setBiquad(bq_type_lowpass, Fc / sampleRate, 0.707, 0);

// filter a buffer of input samples, in-place
for (int idx = 0; idx < bufSize; idx++) {
    in[idx] = lpFilter->process(in[idx]);
}

Here’s an example of cascading two peak filters to yield a 3 dB dip at 200 Hz and a 6 dB bump at 6 kHz:

Biquad *filter1 = new Biquad(bq_type_peak, 200.0 / sampleRate, 1.0, -3.0);
Biquad *filter2 = new Biquad(bq_type_peak, 6000.0 / sampleRate, 5.0, 6.0);

// filter a buffer of input samples, in-place
for (int idx = 0; idx < bufSize; idx++) {
    in[idx] = filter2->process(filter1->process(in[idx]));
}
Posted in Biquads, Digital Audio, Filters, IIR Filters, Source Code | 52 Comments

A wavetable oscillator—end notes

I wrote the wavetable oscillator code from scratch as I wrote this series of articles, using just the thought process I gave. Hence, no references to other articles on wavetable oscillators. The concept of a phase accumulator is fairly obvious; one of the first books I recall reading about it was Musical Applications of Microprocessors, by Hal Chamberlin, about 30 years ago (I have an autographed first edition), though the book didn’t explore the concept of multiple tables per oscillator to control aliasing.

These end notes are an opportunity to touch on various aspects of the oscillator.

Oscillator examples

First, a few examples. Here’s a pulse-width modulation example, using a sine-wave instance of an oscillator to modulate the pulse width of another oscillator:

PWM 110 Hz mod 0.3 Hz sine, 4s

Three detuned and summed sawtooth oscillators:

Saw times 3 detuned 55 Hz, 4s

But we aren’t limited to computing our wavetables. Here, I recorded myself singing “aahh…” (originally, 110 Hz, the second A below middle C). I cut one cycle of that in an audio editor and stretched (resampled) it to 2048 samples, and did an FFT of that. From there it was easy to create higher octaves by eliminating harmonics in the upper octave and doing an inverse FFT to make the next higher wavetable—the same process as we did with the computed waves. Here it is, swept 20 Hz to 20 kHz (but keep in mind that the typical useful range is in the lower octaves):

Ahh sweep 20-20k, 20s

Sawtooth revisited

Because I wanted to show the harmonic series of sine waves building a sawtooth in the earlier articles, I didn’t mention that sawtooth waves traditionally ramp up (remember the charging capacitor in the analog oscillator?). Creating this wave from sines make it ramp down—often called an inverted sawtooth. It really doesn’t make a difference in the audio range; in the sub-audio range for a control signal, you might prefer one or the other depending on the effect you’re after. But either way, it’s trivial to build the up-ramping sawtooth, and I did that in the final version just for fun. I left out the 20-40 Hz octave in the tutorial to make a point, but I’m including it here—a non-inverted sawtooth sweep, and you’ll notice that the lowest octave is brighter with a true 20 Hz table:

Saw sweep 2048 lin 20-20k (base=20Hz), 20s

More on low-end brightness

Though we don’t often listen to the audio output of an oscillator in the sub-audio frequency range, the harmonic content is apparent for the brighter waveforms such as sawtooth. Here’s what it sounds like when we sweep from audio range down to 1 Hz with our wavetables built for the audio range:

Saw sub-audio test (20-20k 2048), 20s

Note how the harmonic content drops as the ticks slow.

Because the lower the frequency, the less that aliasing is a factor (since the sample resolution for a single cycle increases the lower we go), we could switch to a more straight-forward sawtooth for sub-audio. Or, because our oscillator allows a variable table size for each subtable, we could simply place a long wavetable to handle all sub-audio. Here’s the same oscillator, but with a 32768-sample sawtooth table to handle everything below 20 Hz:

Saw sub-audio test (20-20k 2048 + 32768 ramp), 20s

Optimizations

My goal with the accompanying code was to make it simple and understandable. It’s pretty quick, though—the full audio sweep takes about 0.9 seconds to generate 1000 seconds of audio on my computer. And there are certainly ways to make it faster, though in general a little bit faster makes it a lot less readable, and often requires some loss of generality. For instance, we could make oscillator table selection faster by mandating octave spacing, then using a fast-log2 floating-point trick to to get to the proper table in a single step—but that would put a limitation on the oscillator, the relevant code would become unreadable to many people, and in practice the speed gain is small.

Some people might prefer a fixed-point phase accumulator. The math is fast and precise, and you can get a free wrap-around when you increment. But the double-precision floating point phase accumulator of this oscillator is extremely precise, is plenty fast, and the code is very easy to follow.

Interpolation

If you search the web for wavetable oscillators, you’ll find terrific educational papers and forum threads discussing various types of interpolation and their relative performance. It may leave you wondering why you should trust the lowly and relatively pathetic linear interpolation presented here. Good point.

It’s about table size. If your table is big enough, you don’t need any interpolation at all—truncation is fine. But linear interpolation is only a bit more costly, and gives a boost that’s probably worth the effort. Higher forms of interpolation require more sample points in the calculation (for instance, needing two points on each side of the point to be calculated, instead of one point on each side for linear interpolation). They are great for squeezing better quality out of a small table, but you need to ask yourself why you are using the small tables. If you’re implementing this on a modern computer, what is the point of increasing your calculations so that you can use a 512-sample table instead of a 2048-sample table—on a systems that has gigabytes of free memory?

Constant versus variable table size

If we use a constant subtable size, as each subtable moves up to cover the next octave, with half of the harmonic content, they become another factor of 2 x oversampled. The top octave will always be a sine wave, and at the 2048 samples we’re using for our subtables, it’s extremely oversampled—way more than we need or will be able to discern.

And while I said that we want to be at least 2 x oversampled at our lowest subtable, for linear interpolation to have good fidelity, at the same time it’s unlikely that we’ll hear it in the lowest octave. The reason is that low frequency wavetables for a sawtooth will be very crowded in the high harmonics—half of them are in the top octave of our hearing, spaced closely. So, we could go to 1024-sample tables for 40 Hz and up, or 2048 for 20 Hz and up (as in the “Saw sweep 2048 lin 20-20k” example above), and you probably won’t hear the difference. All of the higher tables will still be increasingly oversampled. Here’s the 1024 equivalent of the 2048 example in Part 3 of the wavetable oscillator series:

Saw sweep 1024 trunc 20-20k (base=20Hz), 20s

In a nutshell, we might suspect that it will be easier to hear the benefits of oversampling in the upper octaves, but at the same time the minimum table size to hold their required harmonic content is smaller.

We could just as easily go to a constant oversampling ratio by dropping the subtable length by a factor of two as we go to the next octave’s subtable. It’s easy enough to make all of these things variable in your code, so that you can play with tradeoffs.

Lesser interpolation experiments

Linear interpolation is referred to as “first order” interpolation. Truncation is “zero order” interpolation, implying none. So let’s hear how awful the oscillator would sound with no interpolation:

Saw sweep 2048 trunc 20-20k, 20s

Um, that’s not as bad as expected—or should we not have expected it to be bad? Remember, we’re using fixed table sizes so far, which makes the higher, more sensitive tables progressively more oversampled. And oversampling helps improve any interpolation, even none.

Let’s explore variable table sizes, to keep a constant oversampling ratio per octave. We’ll start with the worst-case of 1 x oversampling, and with truncation:

Saw sweep variable 1x trunc 20-20k, 20s

As expected, the aliasing is bad, and gets worse in the higher octaves. Here it is again, with linear interpolation:

Saw sweep variable 1x lin 20-20k, 20s

Even though 1 x oversampling is still inadequate, the improvement with linear interpolation is obvious. Let’s try again, but with 4 x oversampling; first with no interpolation (and let’s focus on the higher octaves):

Saw sweep trunc 1280-20k, constant 4x oversample, 20s

Now with linear interpolation:

Saw sweep lin 1280-20k, constant 4x oversample, 20s

You can still hear a little aliasing at the top, but 4 x is adequate for much of the range. We might consider 8 x, but a problem with this approach is that we’re wasting a lot of table space for the lower octaves, and not achieving the objective of saving memory by using variable tables. Instead of keeping a constant oversampling ratio, we can go with a variable table size, but impose a minimum length, so that just the top octaves are progressively more oversampled, where we need it most.

Here are constant 1 x oversampled tables, but with a minimum of 64 samples so that the top octaves are oversampled progressively higher; first truncated, then with linear interpolation:

saw sweep 2048 trunc 20-20k (base=20, variable limit 64), 20s

saw sweep 2048 lin 20-20k (base=20, variable limit 64), 20s

A couple of things we’ve learned from these exercises: The upper octaves need progressively more oversampling to sound as good as the lower octaves—at least for a waveform like the sawtooth. And linear interpolation is a noticeable improvement over truncation.

But I hope you’ve also realized that truncation is a viable alternative, as long as the table size is big enough. And the “big enough” qualification is not difficult if you’re using fixed table sizes, which yield progressively higher oversampling in higher octaves.

While variable table size is helpful in memory-resticted environments, it’s not very important in most modern computing environments. I included the ability to handle variable table sizes within an oscillator for educational experimentation, mainly.

Note: Some people are thinking…if truncation might be acceptable, why not use rounding instead? The answer is…there is no difference in noise/distortion between truncation and rounding—use whichever is most convenient (which may depend on the computer language you use or rounding mode of the processor). If you doubt that truncation and rounding are equivalent for our use, consider that rounding is equivalent to adding 0.5 and truncating; this means that rounding in our lookup table is equivalent to truncation except for a phase shift of half a sample—something we can’t hear (because there’s no reference to make it matter, but if it really bothers you, you could always create the wavetables pre-shifted by half a sample!).

Setting oscillator frequency

We set the pitch by setting the increment that gets added to our phasor for every new sample lookup. An increment of 0.5 steps half-way through the table. For a single-cycle table, that would correspond to a frequency of half the sample rate—22.05 kHz, half of our 44.1 kHz sample rate. Another way to say that is that the increment for a frequency of 22.05 kHz is 22050 / 44100, or 0.5.

So, the increment for any frequency is that frequency divided by the sample rate. For 440 Hz, it’s 440 / 44100, or 0.00997732426303855.

On the subject of frequency, one simple and useful addition to the oscillator might be the ability to adjust the phase, for quadrature and phase modulation effects.

Hearing

OK, I hear better than I thought. 15k loud and clear, 16k is plenty clear too, but I could tell it wasn’t quite as loud. 17k is getting pretty tough, needing plenty of volume. But in this case I was listening to naked sine waves at pretty decent volume—I stand by my reasoning that it would be awfully tough to hear those aliased harmonics—even if we were playing ridiculously high notes with no synth filter, at high volume, and soloed. But, if you disagree, just go with closer wavetable spacing!

Posted in Digital Audio, Oscillators, Wavetable Oscillators | 7 Comments

A wavetable oscillator—the code

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:

WaveTableOsc source code

More comments and output examples coming soon.

Posted in Digital Audio, Oscillators, Source Code, Wavetable Oscillators | 15 Comments

A wavetable oscillator—Part 3

In Part 2, we looked at how to use subtables designed to cover one octave each, removing upper harmonics for each higher octave to avoid aliasing. Here’s a view of the nine wavetables for our sawtooth oscillator, with the starting frequency and number of harmonics for each octave:

The table size for each subtable is 2048 samples of single precision floating point. We’re using a double precision floating point phase accumulator. Here’s the result, with an exponential sweep from 20 Hz to 20 kHz:

Sawtooth sweep 20Hz-20kHz, 20s

That sounds…good! Some observations:

  • The wave starts out a little short of harmonics—we knew it would, because we built our lowest table for 40 Hz and up.
  • You probably can’t hear the very end of the sweep, because the frequency is too high.
  • If you play it back on a good audio system and listen carefully, you can hear at least two “ticks” in the audio, nearing the end of the sweep, about 2 seconds apart.

The ticks are due to the abrupt change in energy as we switch to the next table with fewer harmonics. It’s much easier to notice it on in the high range, since the relative power lost is larger (stronger harmonics are dropped). This might be especially noticeable with a strong vibrato that spans tables on a very high note. You can see the transitions in our sweep:

It would be easy to fix—with a higher sample rate, or by crossfading the switch. Or, you can ask yourself how often you’ve ever played sustained notes that high…with strong vibrato…soloed in a mix…with the synth filter missing. The fact is that all synthesizer oscillators make tradeoffs such as this. Either way, we’ll carry on and keep this oscillator simple and easy to understand.

Other waveforms

So far, we’re building the waves computationally. We could sample a single cycle of anything, resample it into a wavetable length of our choice, then use an FFT to get the spectrum, and modify the spectrum as needed for our additional octave tables. This is where the wavetable oscillator has a huge advantage over analog oscillators—any combination of harmonics is possible. For now, we’ll just deal with completing the classic synthesizer waveforms. For most analog synthesizers, that means a sawtooth and a variable-width pulse (rectangle) wave.

The Minimoog had just three choices of pulse width, but most synthesizers have voltage control of pulse width for pulse width modulation (PWM) effects. At first it might seem that we’d need to generate many variations of pulses, but fortunately we can create the same thing by subtracting one sawtooth from another at the same frequency but different phase. Then we get our pulse width modulation by controlling the phase offset. Here’s a square wave (50% pulse) generated using this method:

Square wave sweep 20Hz-20kHz, 20s

The triangle wave was also found in the Minimoog, and is routine for modular analog synths, but was sacrificed in most classic non-modular analog synths because it’s not used as often as sawtooth and pulse. Compared with a sawtooth, a triangle has only odd harmonics, and the harmonics fall off as the inverse square of the harmonics number (and they alternate in sign as well), giving it a hollow sound like the square wave (also odd harmonics), but much more mellow sounding due to the fast rolloff of harmonic energy.

Triangle wave sweep 20Hz-20kHz, 20s

We’ll look more closely at how the oscillator is implemented, as well as give additional observations about wavetable oscillators, in upcoming articles.

Posted in Digital Audio, Oscillators, Wavetable Oscillators | 6 Comments