A quick and accurate method to generate pink noise that I would like to share. The code is available on GitHub here, if you are curious how a few binary operations yield pink noise - or if you just want to hear it - read on.
Pink noise, or 1⁄f noise, has a frequency rolloff of 3 dB per octave, which for audio purposes is often the preferred choice because it is perceived as wide band noise. A reasonable way for generating pink noise digitally is to filter white noise, but unfortunately a 3 dB per octave filter design is not straightforward, so it comes with a compromise between performance and accuracy.
Another common approach found here is to sum a bunch of white noise sources that run on the nominal sample rate divided by powers of two and use zero order hold interpolation for upsampling. This produces lots of aliasing, but as the goal is noise, who cares? But still the spectrum of this emits some unpleasant bumps and valleys.
If linear interpolation, also known as first order hold, is used for upsampling, the spectrum comes much closer to the desired 3 dB rolloff. So this new method adds up twelve independent white noise sources at octave spaced sample rates upsampled by linear interpolation, added to white noise passed through a 12-tap finite impulse response (FIR) filter to fine tune the frequency response, especially close to fs/2. The illustration below shows the twelve upsampled noise sources, the FIR filtered white noise in bright purple, an below the sum of all in a nice shade of pink:
This might look too complicated to be effective, but with white noise coming in the right form, things get surprisingly easy.
Consider white noise consisting of uncorrelated 16-bit numbers, if quantised down to 8 bit, is it still white? Unsurprisingly the answer is yes, the numbers won't get correlated by quantisation. This also holds for quantisation down to a single bit. So once all assumptions that white noise needs to have uniform, gaussian, or whatever distribution are flushed, all that is needed is a source of uncorrelated bits: Welcome linear feedback shift register! A proper LFSR generates a stream of uncorrrelated bits, spectrally flat like white noise if you interpret those as +1 and -1. If the feedback term is chosen properly, the result is a maximum length sequence (MLS) of pseudorandom bits that repeats after 2n-1 cycles. With a sample rate of 44.1 kHz, a 32-bit LFSR runs about a day before repeating, which is good enough for me.
FIR thing first
It is easy to see that a 12-tap FIR filter has 4096 possible outputs if the input signal is just a single bit wide. This can be precomputed and put into a lookup table, but to avoid cache misses, two tables of 64 entries each might be better. The LFSR feedback term is chosen so that twelve successive bits can be used as an 12-bit index. Split into two 6-bit indices, the whole FIR calculation boils down to two table lookups. The purpose of this filter is to minimise the error introduced by the upsampled noise signals. Finding the coefficients requires a moderate amount of math:
The zero order hold sequence g of a sequence s upsampled by m to the nominal rate, noted s↑ here, is
||∑||s↑(n) − s↑(n − m)|
The first order hold / linear interpolation is
|y(n) = m−1||∑||∑||s↑(n) − s↑(n − m)|
With the transfer function
|Y(z) =||m−1||(1 − z−m)2|
|(1 − z−1)2|
For adding up and interpolating the independent bit-noise sequences sk↑ with mk= 2k, the complete signal is
|y(n) =||∑||2−k||∑||∑||sk↑(n) − sk↑(n − 2k)|
As each sk↑ is upsampled by inserting 2k−1 zeros, the gain of sk↑ is 2-k and the squared magnitude for the sum of all is
||H(z)|2 = (1 − z−1)−4||∑||2−3k (1 − z−2k)4|
This squared magnitude, or power as you might call it, is used for calculating the FIR coefficients, which should have the squared magnitude response of
||F(z)|2 ≈ |H(z)|2 - a |z|−1|
There is surely some clever way to calculate the gain value a, iterative guessing worked well for me.
Generating and summing up all noise sources can be made quite effective, for this we start with the zero order hold versions gk(n). We don't really perform the time expansion with zero stuffing and integration to get there, so we rewrite the summation as
|y(n) =||∑||2−k||∑||gk(n) − gk(n − 2k)|
As all gk are single bit signals and any 2-k can be represented by a fixed point number with just a single bit set, both terms 2−k gk(n) and 2−k gk(n − 2k) can be combined into one single word respectively, which simplifies integration a lot. For simplicity we call these two words inc and dec, where each bit represents one of the random bit sequences. Bit bK-k of inc needs to be updated every 2k times, and with a clever scheme like this, only one bit position has to be altered per step. Before updating, the old bit value of inc has to be moved to the same position in dec for the 2k long delay.
So how does it compare? Below I plotted the magnitude error of some pink noise generators I found. The horizontal position of graphs is arbitrarily chosen just to make it readable. The error of the pink noise generator found on simplynoise.com, the self-proclaimed "best free color noise generator on the Internet", is unfortunately too large for this graph. All other algorithms can be found here. The grid is 0.1 dB for magnitude or 0.2 db for squared magnitude.
Below some more information about the specific generators. Error is magnitude ripple or ± deviation from ideal spectral power. Sampling rate is assumed to be 44100 Hz.
|rbj||20 Hz - 21000 Hz||0.61 dB||6||6|
|pke||10 Hz - 22050 Hz||0.91 dB||7||7|
|pk3||10 Hz - 22050 Hz||0.06 dB||14||13||excellent|
|voss||5 Hz - 22050 Hz||1.83 dB||~2||>=11|
|simplynoise.com||15 Hz - 18000 Hz||~2.4 dB||?||?||measured|
|new shade||9 Hz - 22050 Hz||0.04 dB||0||13|
- More octaves can be added at almost no performance penalty, perhaps this can be useful for other applications than audio.
- This is all for discrete time, I wonder if for continuous time, an infinite sum of similar shaped noise sources would yield perfect pink noise.
- Similar schemes for higher order interpolation exist, but those won't bring down the error.
- The Web Audio is wonderful, but here it generates a constant CPU load even if the Play button is not active in the above demo. On a battery-powered device, this might drain the battery. If you know a way using the ScriptProcessor without this issue, please let me know.