0
$\begingroup$

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:

FT by Wolfram

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:

DFT in Python

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

$\endgroup$
1
$\begingroup$

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.

fft output with proper scaling

$\endgroup$
  • $\begingroup$ 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). $\endgroup$ – S V Apr 29 at 21:22
  • $\begingroup$ @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. $\endgroup$ – LedHead Apr 29 at 22:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.