0

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?

user1887919
  • 829
  • 2
  • 9
  • 24

1 Answers1

2

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

Severin Pappadeux
  • 18,636
  • 3
  • 38
  • 64
  • Hmmm ok this is v.helpful - but I am still confused: the link effectively states that h(0) = root pi. I think you have explained it slightly re "DFFT in NUmPy is clearly asymmetrical", but it is still not fully clear to me why h(0) does not equal root pi – user1887919 May 15 '19 at 20:23
  • @user1887919 Well, this is convention used in NumPy. There are libraries which have multiple options for scaling forth and back (like likefftw3 in C and mathdotnet in C#), but in NumPy as far as I can see there are only two options - nothing forward and 1/N inverse scaling (norm=None, which you and I used) or 1/sqrt(N) scaling for both forward and inverse transform ("ortho" option). Neither of those options will produce sqrt(pi). You could probably wrap NumPy FFT into your own routines which do proper scaling in and out. – Severin Pappadeux May 15 '19 at 20:56
  • 1
    @user1887919 You could take a look at https://www.mathworks.com/matlabcentral/answers/15770-scaling-the-fft-and-the-ifft, it is Matlab, but discussion about scaling is quite good – Severin Pappadeux May 15 '19 at 21:05