How To Get Frequency Axis From An Fft Function?

So, I am probably missing something obvious, but I have searched through lots of tutorials and documentation and can't seem to find a straight answer. How do you find the frequency

Solution 1:

Here is an example.

First create a sine wave with sampling interval pre-determined. we will combine two sine waves with frequencies 20 and 40. Remember high frequencies might be aliased if the time interval is large.

#Import the necessary packages
from scipy import fftpack
import matplotlib.pyplot as plt
import numpy as np

# sampling freq in herts 20Hz, and40Hz
freq_sampling1 =10
freq_sampling2 =20
amplitude1 =2 # amplitude offirst sine wave
amplitude2 =4 # amplitude ofsecond sine wave
time= np.linspace(0, 6, 500, endpoint=True) # timerangewith total samples of500from0to6withtimeintervalequals6/500
y = amplitude1*np.sin(2*np.pi*freq_sampling1*time) + amplitude2*np.sin(2*np.pi*freq_sampling2*time)

plt.figure(figsize=(10, 4))
plt.plot(time,y, 'k', lw=0.8)
plt.xlim(0,6)
plt.show()

Notice in the figure that two sine waves are superimposed. One with freq. 10 and amplitude 2 and the other with freq. 20 and amplitude 4.

# apply fft function
yf = fftpack.fft(y, time.size)

amp = np.abs(yf) # get amplitude spectrum 
freq = np.linspace(0.0, 1.0/(2.0*(6/500)), time.size//2) # get freq axis# plot the amp spectrum

plt.figure(figsize=(10,6))
plt.plot(freq, (2/amp.size)*amp[0:amp.size//2])
plt.show()

Notice in the amplitude spectrum the two frequencies are recovered while amplitude is zero at other frequencies. the Amplitude values are also 2 and 4 respectively.

you can use instead fftpack.fftfreq to obtain frequency axis as suggested by tom10 Therefore, the code changes to

yf = fftpack.fft(y, time.size)
amp = np.abs(yf) # get amplitude spectrum 
freq = fftpack.fftfreq(time.size, 6/500)
plt.figure(figsize=(10,6))
plt.plot(freq[0:freq.size//2], (2/amp.size)*amp[0:amp.size//2])
plt.show()

We are only plotting the positive part of the amplitude spectrum [0:amp.size//2]

Solution 2:

Once you feed your window of samples into the FFT call it will return an array of imaginary points ... the freqency separation between each element of returned array is determined by

freq_resolution =  sampling_freq / number_of_samples

the 0th element is your DC offset which will be zero if your input curve is balanced straddling the zero crossing point ... so in your case

freq_resolution = 1000 / 5378

In general, for efficiency, you will want to feed an even power of 2 number of samples into your FFT call, important if you are say sliding your window of samples forward in time and repeatedly calling FFT on each window

To calculate the magnitude of a frequency in a given freq_bin (an element of the returned imaginary array)

X = A + jB

    A on real axis
    B on imag axis

for above formula its

mag = 2.0 * math.Sqrt(A*A+B*B) / number_of_samples

phase  == arctan( B / A )

you iterate across each element up to the Nyquist limit which is why you double above magnitude

So yes its a linear increment with same frequency spacing between each freq_bin