Generating Gaussian Random Image from a Given PSD

# Generate a Gaussian Random 2D Image that has a certain power spectral density

import numpy as np
from numpy.fft import fft2, fftfreq, ifft2
import matplotlib.pylab as plt
from numpy.random import normal

### The PSD that we use to generate the image

def smooth_step(x, w):
    return (np.pi/2+np.arctan(x/w))/np.pi


def ps_model(f, a, f0, alpha, w, e):
    f=np.abs(f)
    a2=a**2
    y=((e + f**2)/(e + f0**2))**(alpha/2)
    s=smooth_step(f-f0, w)
    return a2*y*s+a2*(1-s)

### Image size

nt=32
nf=32

### Frequency value of each pixel

ft=fftfreq(nt)
ff=fftfreq(nf)

### set some parameter of the PSD

fmin_t=ft[1]
a1=3
f01=fmin_t*1
alpha1=-1
fmin_f=ff[1]
a2=1
f02=fmin_f*1
alpha2=-1

### Generate the 2D PSD

ps1=ps_model(ft, a1, f01, alpha1, 1e-9, 1e-9)
ps2=ps_model(ff, a2, f02, alpha2, 1e-9, 1e-9)
#plt.plot(f, ps1)
ps=np.atleast_2d(ps1).T.dot(np.atleast_2d(ps2))
plt.imshow(ps)

<matplotlib.image.AxesImage at 0x7f79621cfcc0>

### Generate the random image

noise_image=ifft2(fft2(normal(size=ps.shape))*ps)

### Check the result, note that the image part is not exactly zero, because the numerical accuracy

plt.figure(figsize=(15,6))
plt.subplot(121)
plt.imshow(noise_image.real, aspect='auto')
plt.colorbar()
plt.subplot(122)
plt.imshow(noise_image.imag, aspect='auto')
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7f7961d44160>