Fast Fourier Transform in Python?

Question

I have this dataset with oscillations:

y = array([ 9.88706879e-05, -1.80853647e-05,  2.42572582e-05, 
1.12215205e-04,
        1.32105126e-04,  1.13424614e-05, -1.58262175e-04, -2.62013276e-04,
       -2.58070932e-04, -1.53975865e-04, -8.19357356e-05, -1.55734157e-04,
       -2.90791620e-04, -3.70294471e-04, -3.46855608e-04, -2.23495910e-04,
       -1.35441615e-04, -2.11411786e-04, -4.21891416e-04, -6.77753516e-04,
       -8.09657243e-04, -6.97948704e-04, -5.01935670e-04, -4.20075723e-04,
       -5.28464040e-04, -8.14942203e-04, -1.03669983e-03, -9.76604755e-04,
       -7.50889655e-04, -5.34882634e-04, -4.06928662e-04, -3.96093220e-04,
       -4.31306957e-04, -4.25399844e-04, -3.26933980e-04, -1.32440493e-04,
        5.40550849e-06, -4.87299567e-05, -2.04672372e-04, -3.15870097e-04])

x = array([-25, -20, -15, -10,  -5,   0,   5,  10,  15,  20,  25,  30, 35,
        40,  45,  50,  55,  60,  65,  70,  75,  80,  85,  90,  95, 100,
       105, 110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165,
       170])

How can I measure the Fast Fourier Transform in Python?

That is how my data look:

and this is the output of the FFT:

N = np.shape(x)
T = (200/N)*10e-12 #time step in sec
xf = fftfreq(N,1/T)[:N//2]
yf = fft(y)
plt.figure()
plt.plot(xf, 2.0/N*np.abs(yf[0:N//2]))
plt.xlabel('freq.(Hz)')
plt.ylabel('f(t)')

I am wondering if the script is correct.

Answer 1

There are two practical mistakes in your original code:

• 10e-12 is not 10^-12; it's 10^-11. Use 1e-12 for 10^-12

• N = np.shape(x) results in N being a tuple, which shouldn't work in the other parts of the code. So the given code fails with a TypeError.

But that won't fix your problem. I'm not too familiar with FFTs, but simplifying your code a bit, combined with some intuition (which might be incorrect), I can get the following with seemingly correct results:

N = len(x)
dt = 5e-12
xf = fftfreq(N, dt)[:N//2]
yf = fft(y)
plt.figure()
plt.plot(xf, 2.0/N * np.abs(yf)[0:N//2])
plt.xlabel('freq.(Hz)')
plt.ylabel('f(t)')

(The "magic" 5e-12 can be obtained in numerous ways. One way is dt = (x[1] - x[0]) * 1e-12, assuming equal spacing for x (which it should have anyway) and converting to seconds.)

That yields the image below. There is a small peak around 3e10 Hz, which equals a period of about 1 / 3e10 = 3.3e-11 seconds or 33 picoseconds. Which seems about right, looking at the x-y data graph in your question.