Checkerboard pattern after FFT

Question

I'm trying to create a Gaussian random field, by creating a grid in Fourier space and then inverse Fourier transorming it to get the random field. For this, the inverse Fourier transformed image needs to be real valued. I seem to be getting residuals in the imaginary part of the grid of the order 10^-18 - -22, so I expected this to be numerical errors in the FFT. The real part of the image displays a weird checkerboard pattern on pixelscale though, where the pixels jump from positive to negative. To see if the FFT functions correctly I tried transforming a Gaussian, which should give back another Gaussian and again the checkerboard pattern is present in the image. When taking the absolute value of the image, it looks fine, but I also need it to allow for negative values for my Gaussian random field.

For the Fourier transformation of the Gaussian I use the following code:

#! /usr/bin/env python

import numpy as n
import math as m
import pyfits


def fourierplane(a):
  deltakx = 2*a.kxmax/a.dimkx #stepsize in k_x
  deltaky = 2*a.kymax/a.dimky #stepsize in k_y

  plane = n.zeros([a.dimkx,a.dimky]) #empty matrix to be filled in for the Fourier grid

  for y in range(n.shape(plane)[0]):
    for x in range(n.shape(plane)[1]):
      #Defining coordinates centred at x = N/2, y = N/2
      i1 = x - a.dimkx/2 
      j1 = y - a.dimky/2

      #creating values to fill in in the grid:    
      kx = deltakx*i1  #determining value of k_x at gridpoint
      ky = deltaky*j1  #determining value of k_y at gridpoint
      k = m.sqrt(kx**2 + ky**2) #magnitude of k-vector


      plane[y][x] = m.e**(-(k**2)/(2*a.sigma_k**2)) #gaussian
  return plane

def substruct():

  class fougrid:
    pass

  grid = fougrid()

  grid.kxmax = 2.00 #maximum value k_x
  grid.kymax = 2.00 #maximum value k_y

  grid.sigma_k = (1./20.)*grid.kxmax #width of gaussian

  grid.dimkx = 1024
  grid.dimky= 1024

  fplane = fourierplane(grid) #creating the Fourier grid

  implane = n.fft.ifftshift(n.fft.ifft2(fplane)) #inverse Fourier transformation of the grid to get final image

  ##################################################################
  #seperating real and imaginary part of the Fourier transformed grid
  ##################################################################

  realimplane = implane.real
  imagimplane = implane.imag

  #taking the absolute value:
  absimplane = n.zeros(n.shape(implane))
  for a in range(n.shape(implane)[0]):
    for b in range(n.shape(implane)[1]):
      absimplane[a][b] = m.sqrt(implane[a][b].real**2 + implane[a][b].imag**2)

  #saving images to files:
  pyfits.writeto('randomfield.fits',realimplane) #real part of the image grid
  pyfits.writeto('fplane.fits',fplane) #grid in fourier space
  pyfits.writeto('imranfield.fits',imagimplane) #imaginary part of the image grid
  pyfits.writeto('absranfield.fits',absimplane) #real part of the image grid

substruct() #running the script

Does anyone have any idea how this pattern is created and how to solve this problem?

Answer 1

Whenever you see unexpected alternating signs in one DFT domain, it could mean the data in the other DFT domain was rotated halfway through the array (similar to an fftshift). If you have a symmetric "hump" of real values in one domain, then centering that hump on array element 0 (instead of array element n/2) will be the arrangement that most likely won't produce alternating signs in the transform domain.