How to generate 2D colored noise

Question

I'm unsuccessfully trying to generate a 2D array of pink noise. Can anyone explain me how to do this? Otherwise, is there a python library that allows to generate 2D (or higher dimensionality) colored noise (1/f noise)?

Answer 1

The answers given both here and in a related question (How to generate a pink noise image?) give you most of the picture - but I thought I'd lay out a full step-by-step to show how everything comes together.

Especially as it seems people run into issues understanding the "frequency matrix".

To convert white noise to pink noise

1. Generate some white noise: whitenoise = np.random.uniform(0, 1, (256, 256))
2. Fourier transform and shift low-frequencies to centre (this makes generating the "frequency matrix" more intuitive): ft_arr = np.fft.fftshift(np.fft.fft2(whitenoise))
3. Generate the "frequency matrix". In the 1-D case this would be an array of the actual frequencies which correspond to the amplitudes given by the transform. In the 2-D case, this is the distance from the center of our fftshifted fourier space, as the further out to the edges we go, the greater the frequency captured at that point. This can be generated via np.meshgrid and np.hypot :
   
   _x, _y = np.mgrid[0:ft_arr.shape[0], 0:ft_arr.shape[1]]
f = np.hypot(_x - ft_arr.shape[0] / 2, _y - ft_arr.shape[1] / 2)
   
   We then divide the fourier space by the frequencies.
   pink_ft_arr = ft_arr / f
4. Remove singularities that show up due to having a zero-frequency component. (Probably there's a smarter way to do this, but this gave me the expected output regardless):
   
   pink_ft_arr = np.nan_to_num(pink_ft_arr, nan=0, posinf=0, neginf=0)
   
5. Convert the fourier space pink noise back into image space:
   pinknoise = np.fft.ifft2(np.fft.ifftshift(pink_ft_arr)).real

You can then plot this using matplotlibs plt.imshow, or scale it into some sensible range and write it out as an image.

White noise, from np.random.uniform
Pink noise generated by the above.

Answer 2

1. generate white noise in 2D, e.g. using np.random.randn

2. calculate FFT_2D of it (numpy.fft.fft2)

3. multiply result (2D spectrum) by 1/f**2 matrix, calculated in such a way: 1/f_along_x * 1/f_along_y. That's how I interpret this definition of pink noise in N-dimensions, but I am not sure I am correct! EDIT: Cris Luengo version (see the comment below) - 1/sqrt(f_x**2+f_y**2) - seems better for me, but its for you to decide which definition of 2d pink noise to use.

4. Use IFFT_2D to get pink noise image (numpy.fft.ifft2)

Answer 3

Here's some code derived from the answer by ckyleda which plots the landscape in 3D using matplotlib. I use normal white noise, rather than uniform white noise, and I use the modern NumPy random interface.

In my tests, plots using the frequencies from np.hypot are way too spiky. Using the squares of those frequencies, as suggested in Power-law spectra of pink noise from Wikipedia appears to give much nicer results.

""" Pink noise surface in matplotlib
    Written by PM 2Ring 2023.10.20
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

def pink_noise(size, rand, roughness):
    white = rand(size=(size, size))
    white_ft = np.fft.fftshift(np.fft.fft2(white))

    hsize = size / 2
    y, x = np.ogrid[-hsize:hsize, -hsize:hsize]
    freq = (x**2 + y**2) ** (1 / roughness)
    pink_ft = np.divide(white_ft, freq,
      out=np.zeros_like(white_ft), where=freq!=0)

    pink = np.fft.ifft2(np.fft.ifftshift(pink_ft)).real
    lo, hi = np.min(pink), np.max(pink)
    return (pink - lo) / (hi - lo)

def test(size=128, seed=None, roughness=1, scale=1/4):
    rng = np.random.default_rng(seed)
    Z = pink_noise(size, rng.normal, roughness)
    Y, X = np.ogrid[0:size, 0:size]

    #plt.style.use('dark_background')
    fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, figsize=(8, 8), constrained_layout=True)
    ax.plot_surface(X, Y, Z, rstride=1, cstride=1, antialiased=False, cmap=cm.terrain)
    ax.set_box_aspect((1, 1, scale))
    plt.axis('off')
    plt.show()

test(size=128, seed=137, roughness=1, scale=1/4)

Output

Here's the result for the same size and seed, but with roughness=2, so it uses the square roots of the frequencies used in the previous plot.

---

Here's a version with a few more parameters to twiddle. It runs on the SageMath server, but the code is mostly plain Python. It just uses Sage to get the parameters and to do the interactive 3D rendering (which is actually done in the browser by three.js).

Here are the camera controls

//    Orbit - left mouse / touch: one-finger move
//    Zoom - middle mouse, or mousewheel / touch: two-finger spread or squish
//    Pan - right mouse, or left mouse + ctrl/meta/shiftKey, or arrow keys / touch: two-finger move

There's an info menu in the bottom right corner, which can be used to save screenshots, or to save the HTML so you can view the landscape without Sage.

To generate cloudscapes, select the 'Blues_r' colormap, and set the scale to something small, eg 1/20 (make the scale negative if you want to view the clouds from below). Increase the gamma parameter to 2 or 3 for bluer skies.