Fast Digital Sine Oscillator

I originally intended this blog to be about math and digital signal processing as well and not just about C++. Since I recently answered a programming forum question regarding efficient generation of sine waves and took the time to code one in — surprize — C++, I thought to myself, why not sharing the code along with explanations about how it works.

If you’re interested in software synthesizers in general or generating digital representations of DTMF tones and don’t want to waste many CPU cycles when generating sine waves, keep on reading.

Generating sampled sine waves without invoking the math library’s sin function for every sample is actually fairly simple. I will show two approaches of doing that with differing pros and cons and finally a version that combines the advantages of both.

Approach1: Imagine the 2D plane, a unit circle, and a point on that circle that rotates around the origin (see below). If you take the y-coordinate of that point, send this to the soundcard as one sample, rotate the point by an angle of a and repeat these steps you’ll hear a sine wave with a frequency of a*fs/(2*pi) Hz where fs is the sampling rate in Hertz. We only need to compute the rotation matrix in advance depending on the desired frequency. After that we need 4 multiplications and two additions per sample. Not so bad. But it’s also not really that fast.

Since we are only interested in one component of this point coordinate but we’re computing the other one as well, one might ask “Can we improve this and make it faster?”. And yes, there’s a way to do it.

Approach2: Instead of rotating the point by computing the product of a rotation matrix and a vector we can write the next point as a linear combination of the previous two points, see the picture below.

In this example, p2 can be computed using p2 = q + (q-p0) = 2 q – p0 where the point q is simply a scaled version of p1, namely q = cos(a) p1. If we put this together, we’ll arrive at the following linear recurrence equation

The nice thing about it is that to compute the y-coordinate of the new point we don’t need any previous x-coordinates. In addition, we saved a multiplication per component. The factor 2 cos(a) can be computed in advance, so, we’re able to compute additional sine wave samples with only one multiplication and one addition per sample. Actually, this variant is an oscillating IIR filter with a pair of complex conjugated poles on the unit circle. Usually, this is something you would avoid because such filters are unstable or “self-oscillating” to use another term. But in this case we exploit this oscillation property.

Are we done yet? Almost. What if we want to reset the phase of the oscillator or change the frequency without resetting the phase? How is that supposed to work? Admittedly, this is not that easy if we only remember the last two samples and the weighting factor 2 cos(a). Also, we have to keep in mind that calculations with floating point numbers involves rounding errors. It is quite possible that the amplitude of the sine wave will increase or decrease due to rounding errors. So, once in a while we should normalize it. Normalization is also not that trivial if you rely on the recurrence equation alone.

I think the most elegant solution is a combination of both approaches. We can use the current point including x-coordinate and the rotation matrix as the oscillator’s state for easier phase, frequency and amplitude manipulations. But when the user is interested in generating many samples at a time, we unroll the loop to produce a block of samples using the linear recurrence equation. Also, instead of vectors and rotation matrices we can simply use complex numbers. Finally, here’s the C++0x source code of my sine generator. To keep it short but still flexible, I made use of C++0x’s new lambda feature:

#ifndef SINE_GENERATOR_HPP_INCLUDED
#define SINE_GENERATOR_HPP_INCLUDED

#include <cmath>
#include <complex>

class sine_generator
{
    std::complex<double> phase;
    std::complex<double> rotation;

    template<typename Sink>
    void internal_generate(int count, Sink sink);

public:
    sine_generator() : phase(1),rotation(1) {}

    void set_frequency(double rad)
    { rotation = std::polar(double(1),rad); }

    void reset_phase(double rad)
    { phase = std::polar(double(1),rad); }

    template<typename OutputIter>
    void generate(int count, OutputIter begin)
    {
        internal_generate(count,[&](double smp){
            *begin = smp; ++begin;
        });
    }

    template<typename ForwardIter, typename BinOp>
    void generate(int count, ForwardIter begin, BinOp binop)
    {
        internal_generate(count,[&](double smp){
           *begin = binop(*begin,smp); ++begin;
        });
    }
};

template<typename Sink>
void sine_generator::internal_generate(int count, Sink sink)
{
    std::complex<double> rot_forx = pow(rotation,12);
    phase /= abs(phase); // prevent drifting off due to rounding errors
    // use linear recurrence for generating the samples
    double const fy = 2 * real(rotation);
    double const fx = 2 * real(rot_forx);
    double y0 = imag(phase);
    double x0 = real(phase);
    double y1 = imag(phase*conj(rotation));
    double x1 = real(phase*conj(rot_forx));
    double tt; // temporary variable
    while (count >= 12) {
        // 4*3 = 12 samples in a row ...
        for (int k=0; k<4; ++k) {
            // output       | compute next sample
            sink(y0); tt = fy*y0-y1;
            sink(tt); y1 = fy*tt-y0;
            sink(y1); y0 = fy*y1-tt;
        }
        // update x as well ...
        tt = fx*x0-x1; x1 = x0; x0 = tt;
        count -= 12;
    }
    phase = std::complex<double>(x0,y0);
    while (count-- > 0) {
        sink(imag(phase));
        phase *= rotation;
    }
}

#endif // SINE_GENERATOR_HPP_INCLUDED 

Here’s an example program that uses the sine generator:

#include <iostream>
#include <vector>
#include <functional>

#include "sine.hpp"

using std::vector;
using std::cout;
using std::plus;

int main()
{
    const int buffsize = 160;
    vector<double> buffer(buffsize);

    sine_generator sg;

    sg.set_frequency(440 * 2*3.141592653589793 / 8000);
    sg.reset_phase(0);
    // overwrite buffer with sine wave ...
    sg.generate(buffsize,buffer.begin());

    sg.set_frequency(880 * 2*3.141592653589793 / 8000);
    sg.reset_phase(0);
    // mix another sine into the buffer ...
    sg.generate(buffsize,buffer.begin(),plus<double>());

    for (int i=0; i<buffsize; ++i) {
        cout << buffer.at(i) << '\n';
    }
}

– P