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)}{nDelta x} approx frac{1}{nDelta 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

LedHead · Answer

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{2pi}$ 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 $+iomega t$ convention while this Python fft function uses a $-iomega t$ convention, which is why the imaginary component in my plot is flipped compared to Wolfram.

Specifying mesh spacing for DFT in numpy

One Answer

Add your own answers!

Ask a Question