I am trying to use an FFT to numerically obtain the Fourier transform of the function
$ f(t) = \sqrt{\frac{\mu}{\sqrt{\pi}}} \exp\left(-\tfrac{1}{2} (\mu x)^2 \right) $
which is precidely
$F(\omega) = \sqrt{\frac{1}{\mu \sqrt{\pi}}} \exp\left(-\tfrac{1}{2} ( \omega/\mu)^2 \right) $
Note that the input function is real and symmetric and thus I expect the Fourier transform to be real.
To do the Fourier transform I do the following using numpy:
x0 = 20
dx = 0.2
N = 1601
x = np.linspace(-20,20,N)
dx = x[1]-x[0]
mu = 2.0
k0 = 0.0
def norm_gaussian(x,mu=1.0, k0=0.0):
"""Defines a normalized Gaussian"""
return np.exp(-0.5*(x*mu)**2)/np.sqrt(np.sqrt(np.pi)/mu)*np.exp(1j*k0*x)
# Generate the y vector
y = norm_gaussian(x,mu=mu,k0=k0)
# Do the Fourier transform
yt = np.fft.fftshift(np.fft.fft(y,norm="ortho"))
kt = np.fft.fftshift(np.fft.fftfreq(len(x),dx))*2*np.pi
scale = np.sqrt(dx/dk)
ytt = scale*((-1)**np.arange(len(yt)))*yt # Get rid of the (-1)^k
dk = kt[1]-kt[0]
Now here is a plot of the Gaussian in t space Gaussian in t
Here is a plot of what I get from the FFT Gaussian in omega The blue dots are the real part, the green line is the imaginary part, and the orange line is $F(\omega)$.
Note that there is a very significant imaginary part that should not be there. Is there something that I am doing wrong? I would have expected the imaginary part to be zero simply from symmetry considerations.