# Specifying mesh spacing for DFT in numpy

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

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 *(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()


• 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. • It does clarify several issues, but the code you gave does not work if I change for example the interval. Note that you are taking the transform from -100 to 100, whereas I was trying to do it from 0 to 10 (the function is centered at 5 and at 0,10 it takes only 0.193% of its value, and inside the interval 99.959% of the total area, so the approximation within this interval should be pretty good).
– user30058
Apr 29, 2019 at 21:22
• @SalvadorVillarreal If you start at x=0, you need to not do the fftshift on y. How good the approximation is obviously depends on what you're doing with it. The approximation I show here is accurate to several digits. Also, note that the larger the initial range t1-t0, the finer the resolution in frequency even if there's a lot of zeros or near zeros. Apr 29, 2019 at 22:20