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 \(f\), we'll use the equation \(\sin(2 \pi f t)\), where \(t\) is a discrete notion of the current time (which we can calculate using the sample rate \(f_s\) 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 can be seen in analog telephones, where each button generates a unique tone and the algorithm determines which button was pressed. Internally, Goertzel uses a second-order IIR filter that resonates at the target frequency. While its purpose is detection, we can repurpose this resonant structure for generation: if we initialize the filter with the right starting values and feed it no input, it will continue to oscillate at the target frequency, producing a sine wave.
The key insight is a recurrence relation that computes \(\sin(a + nb)\) using only the previous two values in the series. With some algebra and trigonometry, we can derive this equation and use it to build a fixed-frequency sine wave generator.
So, here's the algorithm:
$$ \sin(a + nb) = x \sin(a + n'' b) + y \sin(a + n' b) $$
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 \(sin(a + b) = sin(a) cos(b) + cos(a)sin(b)\), we can rearrange the equation like so:
$$ \begin{aligned} \sin&(a + nb) \\ &= x \sin(a + n'' b) + y \sin(a + n' b) \\ &= x \sin(a + nb - 2b) + y \sin(a + nb - b) \\ &= \sin(a + nb) (x \cos 2b + y \cos b) - \cos(a + nb) (x \sin 2b + y \sin b) \\ \end{aligned} $$
Therefore:
$$ \sin(a + nb)(x \cos 2b + y \cos b - 1) = \cos(a + nb)(x \sin 2b + y \sin b) $$
In order for the equation to hold for all n, we need the following to be true.
$$ \begin{aligned} x \cos 2b + y \cos b &= 1 \\ x \sin 2b + y \sin b &= 0 \end{aligned} $$
Solving for x and y yields \(x = -1\) and \(y = 2 \cos b\). When we apply that to the original equation, we get:
$$ \sin(a + nb) = 2 \cos b \sin(a + n' b) - \sin(a + n'' b) $$
Now, assuming that the initial phase offset for our sine wave generator is 0 (\(a\)):
$$ \sin(b n) = 2 \cos b \sin(b n') - \sin(b n'') $$
Hooray! Given that our sine wave generator has a fixed frequency, \(f\) and therefore \(w\) will not change for each calculation of the next sample. This allows us to cache \(w\) 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 \(S\):
$$ \begin{aligned} S(0) &= 0 \\ S(1) &= \sin(\omega h) \\ S(n) &= c \times S(n - 1) - S(n - 2) \\ \text{where} \\ c &= 2 \cos(\omega h) \\ \omega &= 2 \pi f \\ h &= 1 / f_s \end{aligned} $$
We'll initialize our delay lines with \(S(0)\) and \(S(1)\) along with \(c\). 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!
§ A note on numerical stability
Recursive oscillators like this one can accumulate floating-point error over time, causing the amplitude to slowly drift. However, with f64 precision at typical audio sample rates (44.1kHz or 48kHz), this drift is negligible for any practical use case. You would need to run the oscillator continuously for days before the error becomes measurable, and longer still before it becomes audible.
For f32, or for applications that truly run indefinitely, you could add periodic normalization. For typical audio synthesis applications, f64 without correction is perfectly fine.
§ References
- Eric Cheng and Paul Hudak. Audio Processing and Sound Synthesis in Haskell. 2009.