Fast Sine

Implementing a fast sine wave generator in Rust using Goertzel's algorithm

Thursday, May 14, 2020

In audio synthesis, the sine wave is a fundamental element in many algorithms, sounds, and effects. A sine wave can be used as an LFO (low-frequency oscillator) to control various elements of a sound like its amplitude or frequency. Sine waves can also be used for additive synthesis, such as generating band-limited versions of saw, triangle, or square waves. Many audio processing applications are also real-time applications, so the computation needs to happen quickly and consistently.

While doing research into audio processing and functional reactive programming, I stumbled upon the work of Paul Hudak. His work on the The Haskell School of Music and Euterpea are inspiring. In his paper Audio Processing and Sound Synthesis in Haskell (2009) with Eric Cheng, a section is devoted to the optimizing of signal functions, sine waves in particular.

While the sample code in the paper is in Haskell, we're going to be working with Rust as it's what I use for audio processing.

§ Initial implementation

A straight-forward method for generating a sine wave is to use the standard sin function. In order to generate a sine wave of the frequency $𝑓$ , we'll use the equation $s i n (2 𝜋 𝑓 𝑡)$ , where $𝑡$ is a discrete notion of the current time (which we can calculate using the sample rate $𝑓 𝑠$ and a count of the nth sample).

struct Sine {
    n: usize,
    frequency: f64,
    fs: f64,
}

impl Sine {
    fn step(&mut self) -> f64 {
        let &mut Sine {
            ref mut n,
            frequency,
            fs,
        } = self;

        let t = *n as f64 / fs;
        let value = (2.0 * PI * frequency * t).sin();
        *n += 1;
        value
    }
}

§ The Goertzel algorithm & fixed rate oscillators

The Goertzel algorithm is used to detect a single frequency/tone in a waveform. A common application of the algorithm can be seen in analog telephones where each button generates a unique tone and the Goertzel algorithm is used to detect from the tone which button was pressed. The Goertzel algorithm computes a value of $s i n (𝑎 + 𝑛 𝑏)$ using the previous two values in the series. With some algebra and trigonometry, we can find an equation which holds for all n and can be used to develop a fixed-frequency sine wave generator that creates the next value in the series from the two previous.

So, here's the algorithm:

$s i n (𝑎 + 𝑛 𝑏) = 𝑥 s i n (𝑎 + 𝑛 ″ 𝑏) + 𝑦 s i n (𝑎 + 𝑛' 𝑏)$

In order to use it for our sine wave generator, we need to solve for an x and y that hold for all n. Using the identity $𝑠 𝑖 𝑛 (𝑎 + 𝑏) = 𝑠 𝑖 𝑛 (𝑎) 𝑐 𝑜 𝑠 (𝑏) + 𝑐 𝑜 𝑠 (𝑎) 𝑠 𝑖 𝑛 (𝑏)$ , we can rearrange the equation like so:

$s i n (𝑎 + 𝑛 𝑏) = 𝑥 s i n (𝑎 + 𝑛 ″ 𝑏) + 𝑦 s i n (𝑎 + 𝑛' 𝑏) = 𝑥 s i n (𝑎 + 𝑛 𝑏 - 2 𝑏) + 𝑦 s i n (𝑎 + 𝑛 𝑏 - 𝑏) = s i n (𝑎 + 𝑛 𝑏) (𝑥 c o s 2 𝑏 + 𝑦 c o s 𝑏) - c o s (𝑎 + 𝑛 𝑏) (𝑥 s i n 2 𝑏 + 𝑦 s i n 𝑏)$

Therefore:

$s i n (𝑎 + 𝑛 𝑏) (𝑥 c o s 2 𝑏 + 𝑦 c o s 𝑏 - 1) = c o s (𝑎 + 𝑛 𝑏) (𝑥 s i n 2 𝑏 + 𝑦 s i n 𝑏)$

In order for the equation to hold for all n, we need the following to be true.

$𝑥 c o s 2 𝑏 + 𝑦 c o s 𝑏 = 1 𝑥 s i n 2 𝑏 + 𝑦 s i n 𝑏 = 0$

Solving for x and y yields $𝑥 = - 1$ and $𝑦 = 2 c o s 𝑏$ . When we apply that to the original equation, we get:

$s i n (𝑎 + 𝑛 𝑏) = 2 c o s 𝑏 s i n (𝑎 + 𝑛' 𝑏) - s i n (𝑎 + 𝑛 ″ 𝑏)$

Now, assuming that the initial phase offset for our sine wave generator is 0 ( $𝑎$ ):

$s i n (𝑏 𝑛) = 2 c o s 𝑏 s i n (𝑏 𝑛') - s i n (𝑏 𝑛 ″)$

Hooray! Given that our sine wave generator has a fixed frequency, $𝑓$ and therefore $𝑤$ will not change for each calculation of the next sample. This allows us to cache $𝑤$ along with the previous two samples and therefore compute the next sample of the sine wave with only one multiplication and one subtraction!

§ Fast implementation

Now, let's turn this into code. We can break the above equation down a bit more to simplify the creation of a sine wave generator $𝑆$ :

$𝑆 (0) = 0 𝑆 (1) = s i n (𝜔 ℎ) 𝑆 (𝑛) = 𝑐 \times 𝑆 (𝑛 - 1) - 𝑆 (𝑛 - 2) w h e r e 𝑐 = 2 c o s (𝜔 ℎ) 𝜔 = 2 𝜋 𝑓 ℎ = 1 / 𝑓 𝑠$

We'll initialize our delay lines with $𝑆 (0)$ and $𝑆 (1)$ along with $𝑐$ . Then, for each step of the wave, we'll calculate the next sample s, update the delay lines, and return s.

pub struct Sine {
    c: f64,
    z1: f64,
    z2: f64,
}

impl Sine {
    pub fn with_frequency(frequency: f64, fs: f64) -> Self {
        let (s, c) = (TAU * frequency / fs).sin_cos();
        Self {
            c: 2. * c,
            z1: 0.,
            z2: s,
        }
    }

    pub fn step(&mut self) -> f64 {
        let Self { c, z1, z2 } = *self;

        self.z2 = c.mul_add(z2, -z1);
        self.z1 = z2;

        z1
    }
}

And there we have it!

§ References

Eric Cheng and Paul Hudak. Audio Processing and Sound Synthesis in Haskell. 2009.