CODEWITHSUNDEEP
pythonnumpyfft
Erica Y
Erica Y

Reputation: 41

FFT in Python: formatting 1-D diffraction Fourier transform

I'm relatively new to Python and the FFT function. I'm trying to use the numpy.fft.fft() function to transform a square pulse (1-D diffraction slit function) to a sinc function (1-D diffraction pattern), and make the output plot identical to the analytical transform of the square pulse, given by the equation:

F(u) = sin(πau)/(πu)

However, I am having trouble formatting the FFT to make it look identical to the analytical function. This is my minimal working code:

a = 40.0  # slit width
N = 5000  # sample points
T = 100.0  # sample spacing
x = np.linspace(-T/2.0,T/2.0,N)

y = np.piecewise(x,[abs(x)>a/2,abs(x)<=a/2],[0,1])  # make 1-D square pulse array

dx = x[1] - x[0]
yf = np.fft.fft(y) * dx  # Fourier transform
xf = np.fft.fftfreq(N,d=dx)
xf = np.fft.fftshift(xf)
yf = np.fft.fftshift(yf)
plt.figure(1)
plt.plot(xf,yf)

f = np.sin(np.pi*a*x)/(np.pi*x)  # analytical transform function
plt.figure(2)
plt.plot(x,f)
plt.show()

I think I need to adjust the frequencies for the FFT, but I'm not sure. Any help is greatly appreciated!

Upvotes: 3

Views: 1752

Answers (1)

Warren Weckesser
Warren Weckesser

Reputation: 114831

You have a discrete signal with finite length. When you use the FFT to compute the Fourier transform of that signal, you are assuming that the signal is periodic. That is, your signal is not a single rectangular pulse; it is a repeating pulse. The analytic result that applies in this case is the periodic sinc function (also known as the aliased sinc function or the Dirichlet function), one form of which is

f(x) = sin(n*x/2) / (n*sin(x/2))

This function is available as scipy.special.diric. Just recently I added an example to the docstring of diric to show exactly what you are asking. Here's the example. (It won't appear in a released version of scipy until version 0.17.)

Standard imports:

>>> import numpy as np
>>> from scipy import special

Create a small rectangular pulse for the demonstration:

>>> m = 8
>>> k = 3
>>> x = np.zeros(m)
>>> x[:k] = 1

Use the FFT to compute the Fourier transform of x, and inspect the magnitudes of the coefficients:

>>> np.abs(np.fft.fft(x))
array([ 3.        ,  2.41421356,  1.        ,  0.41421356,  1.        ,
        0.41421356,  1.        ,  2.41421356])

Now find the same values (up to sign) using diric. We multiply by k to account for the different scaling conventions of numpy.fft.fft and diric:

>>> theta = np.linspace(0, 2*np.pi, m, endpoint=False)
>>> k * special.diric(theta, k)
array([ 3.        ,  2.41421356,  1.        , -0.41421356, -1.        ,
       -0.41421356,  1.        ,  2.41421356])

Upvotes: 2

Related Questions