I am trying to understand Discrete Fourier Transform (DFT), after only having experience with the continuous transformation. The natural idea was to try to understand DFT on the simplest function, being an exponential, which should transform into an exponential too. To check this in practice, I've drafted up a Python script to perform DFT on a 6D Gaussian.
Surprisingly, the result seems to have both real an imaginary parts. Based on arguments from the continuous Fourier transform, I would only expect an imaginary part if the function is not even - which the Gaussian should be. The grid is centered symmetrically around zero, so there is no issue with numerical asymmetry either.
Should imaginary numbers appear in the DFT of a six-D Gaussian, and, if yes, what is the significance of this observation?
I believe my Python script is correct, but just in case, here is the implementation:
import numpy as np
from scipy.fft import fftn
def sixdgauss(coords):
return np.exp(-np.pi * (coords[0] ** 2 + coords[1] ** 2 + coords[2] ** 2 +
coords[3] ** 2 + coords[4] ** 2 + coords[5] ** 2))
nx, ny, nz = 8, 8, 8
xx, yy, zz = np.linspace(-2, 2, nx), np.linspace(-2, 2, ny), np.linspace(-2, 2, nz)
nxp, nyp, nzp = 8, 8, 8
xxp, yyp, zzp = np.linspace(-2, 2, nxp), np.linspace(-2, 2, nyp), np.linspace(-2, 2, nzp)
XX, YY, ZZ, XXP, YYP, ZZP = np.meshgrid(xx, yy, zz, xxp, yyp, zzp)
gauss = sixdgauss((XX, YY, ZZ, XXP, YYP, ZZP))
print(fftn(gauss))
Edit:
As it had been pointed out to me in the comments below, the issue is already present for the 1D case; running xx = np.linspace(-4,4,16); fft( np.exp(-np.pi * xx**2) ) results in imaginary result too.
You should use simpler examples when you’re looking into such things. The same thing already happens in one dimension.
The transform is real if the input is symmetric about the origin. The origin is at $0$, not in the middle of the interval. The transform doesn’t know that you used the coordinate value $0$ to generate the function value in the middle of the interval.