Specifying mesh spacing for DFT in numpy

Question

I was testing the .fft package of numpy 1.16.1 in Python 3.7.2. In particular I was trying to verify that the transform resembles the analytical one for: $$f(x) = \mathrm{exp}\left[-\left(\frac{x-5}{2}\right)^{2}\right]$$

I get from Wolfram Alpha that $\hat{f} = \mathcal{F}[f]$ looks like this:

Then I tried to replicate this plot with numpy and matplotlib, with the following code:

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(0, 10, 1/1000)
y = np.exp(-((x-5)**2)/4)

y_hat = np.fft.fftshift(np.fft.fft(y))
re_y_hat = np.real(y_hat)
im_y_hat = np.imag(y_hat)

fig, ax = plt.subplots()
ax.plot(x, re_y_hat, "b-", x, im_y_hat, "r-")
plt.show()
plt.close()

But the image I obtain differs greatly from the one Wolfram gives:

In the last image the zero frequecy was shifted to the center by using np.fft.fftshift() so the spike corresponds to frequency zero.

I already figured out that the problem is that nowhere in np.fft.fft() is the $\Delta x$ being specified, so what numpy is interpreting is that I have data that varies very slowly, almost constant$^{1}$, and thus the transform is close to that of a constant function.

I looked at the numpy documentation and other SE posts to see how this can be fixed but found nothing. Does anyone know how to fix this?

---

$^{1}$ We can calculate the average slope of the function numpy sees by $\frac{\mathrm{max}\{f\}-\mathrm{min}\{f\}}{x_{f_{\mathrm{max}}}-x_{f_{\mathrm{min}}}} = \frac{f(5)-f(0)}{n\Delta x} \approx \frac{1}{n\Delta x}$ where $n$ is the number of nodes separating the maximum from the minumum. In this case, since numpy takes $\Delta x = 1$ by default, the slope is about 1/5000=0.0002

Answer 1

I'd love to say that I understand exactly all of the scale factors and shifts that I did below, but I mostly played around with factors until things matched :) I would wait until I've thought it all through, but I wanted to get this solution to you, so forgive me for the incomplete answer.

import numpy as np
import matplotlib.pyplot as plt

fs=1e2
t0=-100
t1=100.0
x = np.arange(t0,t1, 1./fs)
y = np.exp(-((x-5)**2)/4)

y_hat = np.fft.fftshift(np.fft.fft(np.fft.fftshift(y))) / x.shape[0] *(t1-t0)/np.sqrt(2.0*np.pi)
freq=np.arange(-fs/2,fs/2,1.0/(t1-t0))
omega=2*np.pi*freq
y_hat_exact=np.sqrt(2.0)*np.exp(-omega**2-5j*omega)

plt.ion()

ff=plt.figure(1)
ff.clf()
ax=ff.gca()
ax.plot(x,y,'b-')
ax.set_xlabel('x')
ax.set_ylabel('y(x)')

ff=plt.figure(2)
ff.clf()
ax=ff.gca()
ax.plot(omega, y_hat.real, "b.-",label='re fft')
ax.plot(omega, y_hat.imag, "r.-",label='im fft')
ax.plot(omega,y_hat_exact.real,'k--',label='re exact')
ax.plot(omega,y_hat_exact.imag,'k:',label='im exact')
ax.set_xlabel('Radial frequency $\omega$ (rad/s)')
ax.set_ylabel('F[y]')
ax.set_xlim([-1.0,1.0])

ax.legend()

plt.show()

Some comments:

• You can't use x as the independent variable for plotting the FFT, you need frequency (or radial frequency). The code shows how to calculate that.
• You have to divide the result by the number of samples. This is a convention of an FFT.
• Dividing by $\sqrt{2\pi}$ is to match the convention Fourier transform convention used by Wolfram Alpha
• The scaling by the total time t1-t0 was figured out experimentally. If I think about it long enough, I'm sure it has to do with the fact that the analytical Fourier transform is an integral over this time period.
• The input to the FFT needs to start from x=0. If there's non-zero components in negative x, then they need to be wrapped around to the end, hence the fftshift of y
• Your Wolfram Alpha link and corresponding image are for $f(x) = e^{-(x-2)^2/4}$ while your code is for $f(x) = e^{-(x-5)^2/4}$. My example uses the latter.
• Wolfram Alpha uses a $+i\omega t$ convention while this Python fft function uses a $-i\omega t$ convention, which is why the imaginary component in my plot is flipped compared to Wolfram.