Colouring Noise

Generating coloured noise to simulate physical processes

I’ve been characterising a high performance gyroscope at work recently by measuring the Allan variance of its rate output. From the Allan variance plot, it’s possible to visualise and calculate the various sources of noise that cause error, including quantisation noise, angular random walk, bias instability and rate random walk. Each of these noise sources is caused by a different physical process and can be differentiated by its colour. What is noise colour, you ask?

A dark plot with 5 data series on it for violet, blue, white, pink and brownian noise, each coloured in their respective colours.
The power spectral density of some of the most common colours of noise

The colour of noise refers to the shape of its power spectral density (PSD). Over time, people have assigned colours to particular shapes in the power spectrum.

Ultimately, these colours of noise are often approximations of what we observe over specific frequency bands, and the full PSD of a process can contain multiple colours of noise mixed together. For example, the input referred voltage noise density (which is proportional to the square root of the PSD) of a typical opamp is shown below. This noise could be described as pink noise to around 200–300 Hz, and white noise up to around 50 kHz, but the final measured noise density is more complex.

A typical noise plot from an opamp datasheet showing two curves linearly decreasing in log-log space up to a point where they become horizontal. The two curves represent the voltage noise density for different supply voltages.
The input referred voltage noise density for an ADA4625-1 opamp

Generating Coloured Noise Samples

The process of generating coloured noise is relatively simple:

  1. Start with a representation of white noise signal in the frequency domain
  2. Shape the signal in the frequency domain according to the PSD you’d like to achieve
  3. Take the inverse Fourier transform of the shaped signal

We’ll take a look at how we can easily do this in Python using numpy, scipy and matplotlib.

Base Noise

To begin, we need to generate some white noise in the frequency domain. Ideally this noise will be normalised to have a one-sided PSD of 1.0 at all frequencies which will simplify the shaping process in the next step.

There are couple of ways to do this, with the simplest being to generate the noise in the time domain and take a Fast Fourier Transform (FFT). The PSD of uniformly sampled discrete white noise is equal to $2 \sigma ^ 2 \Delta t$, where $\sigma ^ 2$ is the variance of the noise, and $\Delta t$ is the sampling period. Let’s quickly stub out a class using this knowledge to generate our base noise.

import numpy as np

class NoiseGenerator:
    rng = np.random.default_rng()

    def generate(self, dt, n, colour=None):
        """Generates uniformly sampled noise of a particular colour.

        dt : float
            Sampling period, in seconds.
        n : int
            Number of samples to generate.
        colour : function, optional
            Colouring function that specifies the PSD at a given frequency. If
            not specified, the noise returned will be white Gaussian noise with
            a PSD of 1.0 across the entire frequency range.
            Array of sampled noise, length `n`.

        if colour:
            raise NotImplementedError("We're not ready yet!")
        f, x_f = self._base_noise(dt, n)

        return np.fft.irfft(x_f)

    def _base_noise(self, dt, n):
        """Produces a frequency domain representation of uniformly sampled
        Gaussian white noise.
        # Generate Gaussian white noise in the time domain
        sigma = 1 / np.sqrt(2.0 * dt)
        x_t = rng.normal(0.0, sigma, n)

        # Calculate the frequency bins and associated FFT
        f = np.fft.rfftfreq(n, dt)
        x_f = np.fft.rfft(x)

        return f, x_f

It’s worthwhile checking here to make sure our generated noise has a constant PSD of 1.0. Welch’s method1 is a common technique to estimate the PSD of a time domain signal. Let’s run some code to calculate a PSD estimate and plot it along with the noise signal.

import scipy.signal
import matplotlib.pyplot as plt

def plot_psd(x_t, dt, **kwargs):
    """Plots a signal and its PSD given samples and sampling period."""

    t = np.linspace(0.0, dt * len(x_t), len(x_t), endpoint = False)
    f, s_f = scipy.signal.welch(x_t, fs=1/dt, **kwargs)

    plt.subplot(2, 1, 1)
    plt.plot(t, x_t, linewidth=0.5)
    plt.xlabel("Time (s)")
    plt.ylabel("Value (V)")
    plt.subplot(2, 1, 2)
    plt.loglog(f, s_f, linewidth=0.5)
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("PSD ($V^2/Hz$)")

# Generate white noise and plot the PSD
dt = 1.0e-3
n = 1000000
noise_generator = NoiseGenerator()
x_t = noise_generator.generate(dt, n)
plot_psd(x_t, dt, nperseg=65536)

A figure with two subplots containing a single blue trace each. The top shows a noisy time based signal that is around 1000 seconds long, and has peaks up to a value of 100 V or so. The bottom is a logarithmic plot showing a noisy power spectral density estimate that is relatively close to 1.0 between the frequencies of 0.001 to 500 Hz.
The noise signal and its power spectral density

Judging from the plot above, we do indeed have a white noise signal with a PSD close to 1.0 over most of the frequency band.

Those that looked at the implementation of NoiseGenerator closely would’ve noticed the inefficiency of generating random noise in the time domain, then transforming it to the frequency domain, then transforming it back into the time domain. It’s possible to generate the white noise directly in the frequency domain2. Here’s an alternative implementation of _base_noise() that does exactly that.

    def _base_noise(self, dt, n):
        # Calculate random frequency components
        f = np.fft.rfftfreq(n, dt)
        x_f = self.rng.normal(0,
            0.5, len(f)) + 1j * self.rng.normal(0, 0.5, len(f))
        x_f *= np.sqrt(n / dt)

        # Ensure our 0 Hz and Nyquist components are purely real
        x_f[0] = np.abs(x_f[0])
        if len(f) % 2 == 0:
            x_f[-1] = np.abs(x_f[-1])
        return f, x_f

Shaping the Spectrum

At this point, all that is left is to colour the noise by shaping its spectrum. Updating generate() to handle a colouring function is a suitably generic way to do this. We’ll add some helper static functions to the class for common colours, but of course you’re welcome to do something completely custom. There are really only two things to note here:

  1. The PSD is directly proportional to the square of the FFT, so we need to take the square root of our colouring functions.
  2. A few noise colours have $f$ in the denominator of a fraction, so we need to specifically handle the case that $f = 0$.
    def generate(self, dt, n, colour=None):
        f, x_f = self._base_noise(dt, n)

        if colour:
            x_f *= np.sqrt(colour(f))

        return np.fft.irfft(x_f)
    def pink(scale=1.0):
        """Creates a pink noise colouring function.

        scale : float, optional
            Multiplier to adjust the scale of the pink noise.
            Pink noise colouring function that can be used with `generate()`.
            The function returned will be :math:`y(f) = s / f`, where :math:`s`
            is `scale`. At f = 0, the function will simply return 0.0.

        return lambda f: scale / np.where(f == 0.0, np.inf, f)

I’ve omitted most of the colours in the code snippet above, but the full NoiseGenerator class can be found in this gist. Of particular interest is the piecewise colouring function, which allows creation of some funky noise that can’t be classified by any of the standard colours.

ng = NoiseGenerator()
noise_settings = {
    "Piecewise": {
        "dt": 1e-4,
        "n": 1000000,
        "colour": ng.piecewise_logarithmic(
            [90.0, 100.0, 110.0, 450.0, 500.0, 550.0],
            [0.1, 10.0, 0.01, 0.01, 2.0, 0.001]
    "Brownian": {"dt": 1e-3, "n": 100000, "colour": ng.brownian()},
    "Azure": {"dt": 1e-3, "n": 100000, "colour":}
for name, settings in noise_settings.items():
    x_t = ng.generate(**settings)
    plot_psd(x_t, settings["dt"], nperseg=65536)

A dual plot showing 3 noise signals in the time domain and their power spectral density. A piecewise noise with two distinct spikes in the PSD, brownian noise with its distinct 20 dB / decade falloff and azure noise increasing in PSD from a very low value.
Some examples of coloured noise, including noise with a custom piecewise power spectral density

Questions Raised

Although the coloured noise generation class worked perfectly for my needs, a few questions were raised while I was using it. I’ve left them unanswered for the moment, but if anyone has any answers, please get in touch.

Noise Distribution

The probability density function used in the random number generation above is Gaussian. It’s important to note that the noise distribution is a separate concept entirely to its power spectral density. For our base noise generation in the time domain, we could have used a uniform distribution, or an exponential distribution, or a gamma distribution, and they all would’ve worked.

Interestingly, after taking a Fourier transform, all of these distributions produce something that looks roughly like a 2D Gaussian distribution in the complex plane. Furthermore, after shaping the PSD using a colouring function and taking an inverse Fourier transform, we always seem to end up with noise in the time domain that tends towards a Gaussian distribution.

  1. Is it possible to generate noise with a specific distribution in the frequency domain?
  2. Why does shaping noise in the frequency domain tend to shift the time domain distribution towards a Gaussian distribution? Is it something to do with the central limit theorem?
  3. Is it possible to generate noise with both a specific distribution and a specific colour that isn’t white?

Six separate plots. In the top left, a plot with two traces - one showing gamma distributed noise (which has all samples greater than 0) and one showing pink noise which is distributed evenly each side of the x-axis. In the bottom left, standard looking PSDs of white and pink noise. In the top right are two time domain histograms - one showing the gamma distribution, one showing what looks like a Gaussian distribution. In the bottom right are frequency domain 2D histograms which both look roughly Gaussian.
On the left, time and PSD plots of base white noise with a gamma distribution and pink noise shaped from that base noise. On the right, time and frequency domain histograms of the same. The frequency domain histograms are on the complex plane.

Integration of Noise

I was using the synthetically generated noise to simulate a gyroscope, so I’d integrate the signal to calculate the angular error of the system. In this way, I could test out bias estimation and compensation algorithms without having to acquire long datasets from a real gyroscope. I noticed some strange behaviour in the integrated signal, though. For all coloured noise that had a 0 Hz frequency component of 0.0, the resulting integral would have both a start and end value of 0.0. This of course meant that if I generated 1 million noise samples, my angular error at the millionth sample would always be 0.0.

A plot of six random-walk-like signals that begin with a value of 0.0 and end with a value of 0.0, but are uncorrelated otherwise.
The integral (sum) of six sets of 1 million points of generated pink noise with a PSD of 0.0 at 0 Hz

What is going on here? My current theory is that it is due to the definition of the discrete Fourier transform:

$$X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-\frac {i 2\pi}{N}kn}$$

If the Fourier transform is 0.0 for $k=0$, the above simplifies to:

$$0.0 = \sum_{n=0}^{N-1} x_n$$

Of course, this is exactly the oddity that we see - the sum of our noise samples is 0.0. It does raise the question of whether setting the 0 Hz component in our colour shaping functions to 0.0 is sensible. Perhaps it should be some other value?

Final Words

This small dive into noise generation techniques resulted in a NoiseGenerator class that was useful for gyroscope simulations, as well as other applications. At just 44 lines of functional code, it is testament to the power of the scientific computing libraries in the Python ecosystem. Although generating white noise and filtering it seems to be an accepted technique for coloured noise generation, it did raise some questions on the structure of the generated noise which I hope will be answered in the future.

  1. Welch, P. (1967). The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics, 15(2), pp. 70–73 ↩︎

  2. Timmer, J. and König, M. (1995). On Generating Power Law Noise. Astronomy and Astrophysics, 300, pp. 707–710 ↩︎