Search code examples
pythonfftcontinuous-fourier

Fourier transform of a Gaussian function in Python

I want to calculate the Fourier transform of some Gaussian function. Consider the simple Gaussian g(t) = e^{-t^2}. The Fourier transform of g(t) has a simple analytical expression , such that the 0th frequency is simply root pi.

If I try to do the same thing in Python:

N = 1000
t = np.linspace(-1,1,N)
g = np.exp(-t**2)

h = np.fft.fft(g) #This is the Fourier transform of expression g

Simple enough. Now as per the docs h[0] should contain the zero frequency term, which we know from the analytical expression is root pi. But instead it gives 746.444?!

Why the discrepancy between the analytical solution and the computational one?

Solution

  • Not sure why do you think you should get analytical expression. DFFT in NUmPy is clearly asymmetrical, and if you look at formula for Ak here, you could clearly see that for A0 you ought to get sum of input. Also, having gaussian from [-sigma...sigma] interval is not right.

    Here is modified example

    import numpy as np
import matplotlib.pyplot as plt

N = 4001
t = np.linspace(-4.0, 4.0, N)
print((t[0], t[2000], t[4000]))
g = np.exp(-t*t)
print(np.sum(g)) # sum of input

h = np.fft.fft(g, norm=None)
print(h[0]) # should be the same as sum of input

    and it prints

    (-4.0, 0.0, 4.0)
886.2269119018041
(886.226911901804+0j)

    You could do inverse transform and plot it

    q = np.fft.ifft(h, norm=None)

plt.plot(t, g, label = "Gauss")
plt.show()
plt.plot(t, np.abs(q), label = "dFFT Gauss")
plt.show()
f = np.fft.fftfreq(N)
plt.plot(f, np.angle(h), f, np.abs(h))
plt.show()

    and get

    enter image description here