0

So I wrote a short Python program to estimate the accuracy of the Python's FFT method.

import numpy as np
import matplotlib.pyplot as plt

#Aufgabe 1
x0=0
a=2.5
k0=3
X=np.linspace(-4,4,100)
timestep=0.1
k=np.fft.fftfreq(X.size,d=timestep)
psi_analytical=[(2/(np.pi*a**2))**(1/4)*np.exp(-((i-x0)**2)/a**2)*np.exp(1j*k0*(i-x0)) for i in X]
psi_tilde_numerical=np.fft.fft(psi_analytical)
psi_tilde_analytical=[(2/(np.pi*a**2))**(1/4)*(a/2)*np.exp(-(a*(i-k0))**2/4)*np.exp(-1j*i*x0) for i in k]
psi_numerical=np.fft.ifft(psi_tilde_analytical)


#plt.plot(k,np.abs(psi_tilde_numerical),label='numerical psi tilde')
#plt.plot(k,np.abs(psi_tilde_analytical),'--',color='tab:orange', label='analytical psi tilde')

plt.plot(X,np.abs(psi_analytical),label='analytical psi, real')
plt.plot(X,np.abs(psi_numerical),'--',color='tab:orange',label='numerical psi, real')
plt.legend()
plt.show()

The analytical function is as follows:

Analytical functions

To my surprise, the numerical and analytical functions are totally different. However, I'm not sure why this is the case.

The normalisation constant N is (2/(np.pi*a**2))**(1/4)

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
Jung
  • 145
  • 6
  • It looks like the analytical functions are continuous while FFT is discrete in nature. There is also the issue of only looking at a window of the infinite function. You are effectively multiplying the function with a box which results in convolution in the frequency domain. – mandulaj Jan 06 '21 at 20:38
  • If you plot your sampled input function and FFT result, and add them to the question, it'll be easier for us to point at possible issues. These answers might help you find your problem: https://stackoverflow.com/a/49331862/7328782 , https://stackoverflow.com/a/49142862/7328782 . Especially the first one, since I see you have a time axis with 0 in the middle, but haven't used `ifftshift` in your code. – Cris Luengo Jan 06 '21 at 20:47

1 Answers1

0

Doing a bit more research, I think I might have an answer for you.

Discrete Fourier Transform (DFT), which is computed efficiently using the Fast Fourier Transform algorithm (FFT), operates on discrete time domain signals. The Fourier Transform (FT) operates on function in continuous time domain.

DFT will approximate the FT under certain condition. One of those conditions is that the signal has to be band limited. This means that the FT for the function has to be zero for all frequencies above a certain frequency threshold α and the DFT has to have a sample rate that is at least 2*α This goes back to Nyquist-Shannon sampling theorem.

In your case, you are trying to approximate a Gaussian function exp(-x²) which is not band limited. This is because as you can see from your formulas, FT of a Gaussian is also a Gaussian. This means that it has negligible but non-zero components for frequencies all the way to infinity. As a result, you won't be able to approximate the FT using a DFT since you would need to have infinite sampling rate.

In conclusion, its important to realize that the DFT and FT are vastly different transforms and thus can not just be compared.

mandulaj
  • 733
  • 3
  • 10