Using Numpy's FFT in Python
There are numerous ways to call FFT libraries both in Numpy, Scipy or standalone packages such as PyFFTW. In this post, we will be using Numpy's FFT implementation.
Definition and normalization
As usual, we want to make sure we get the definition right, as the normalization coefficients or the sign of the exponent can be arbitrary.
As per the official documentation, the DFT is defined as
while the inverse DFT is defined as
We can check that this matches the definition we gave in the previous section.
fftshift and fftfreq
In many applications of the FFT, we are not really interested in the raw indices or , but rather, on the true frequencies that these indices represent.
In order to obtain those true frequencies, we must shift the zero-frequency component to the center of the spectrum and reorder the coefficients, as we discussed on the previous section along with the almost identical topic of aliasing and negative frequencies, which are a common source of confusion when implementing algorithms that rely on the FFT.
Let’s see how the fftshift
function reorders a vector of size N
by considering two simple examples, one for an even value of N
, and an other for an odd value of N
:
import numpy as np
from numpy.fft import fft, fftshift, fftfreq
fftshift(np.arange(0,10))
array([5, 6, 7, 8, 9, 0, 1, 2, 3, 4])
fftshift(np.arange(0,11))
array([ 6, 7, 8, 9, 10, 0, 1, 2, 3, 4, 5])
Another useful function, fftfreq
gives the correct frequencies (that is, below Nyquist) that the FFT algorithm is producing. By default, this function will return the frequencies divided by N
, so we multiply by N
in order to get the shifted indices we desire (see the docs on fftfreq).
N = 10
fftfreq(N)*N
array([ 0., 1., 2., 3., 4., -5., -4., -3., -2., -1.])
Now lets combine these two functions to compare the original with the shifted indices and coefficients.
import numpy as np
from numpy.fft import fft, fftshift, fftfreq
from matplotlib import pyplot as plt
N = 21
xj = np.arange(0,N)*2*np.pi/N
f = 2*np.exp(17*1j*xj) + 3*np.exp(6*1j*xj) + np.random.rand(N)
original_k = np.arange(0,N)
shifted_k = fftshift(fftfreq(N)*N)
original_fft = fft(f)
shifted_fft = fftshift(fft(f))
fig, axs = plt.subplots(2,1)
axs[0].plot(original_k, np.abs(original_fft))
axs[0].plot([0,6,17], np.abs(original_fft[[0,6,17]]),"o")
axs[0].set_title("Original FFT coefficients")
axs[0].set_xticks(original_k)
axs[1].plot(shifted_k, np.abs(shifted_fft))
axs[1].plot([-4,0,6], np.abs(shifted_fft[[6,10,16]]),"o")
axs[1].set_title("Shifted FFT coefficients")
axs[1].set_xticks(shifted_k)
plt.tight_layout()
plt.show()
Sampling Rate and Frequency Spectrum Example
Using the functions fft
, fftshift
and fftfreq
, let’s now create an example using an arbitrary time interval and sampling rate.
We will use a sampling rate of 44100 Hz, and measure a simple sinusoidal signal for a total of 0.1 seconds.
import numpy as np
from numpy.fft import fft, fftshift, fftfreq
from matplotlib import pyplot as plt
t0 = 0 # Start time
fs = 44100 # Sampling rate (Hz)
tmax = 0.1 # End time
t = np.arange(t0, tmax, 1/fs)
signal = np.sin(2 *np.pi * 60 * t)
F = fftshift(fft(signal))
freqs = fftshift(fftfreq(len(t), 1/fs))
fig, axs = plt.subplots(2,1)
axs[0].plot(t, signal)
axs[0].set_title("Signal")
axs[1].plot(freqs, np.abs(F))
axs[1].set_title("Spectrum")
axs[1].set_xlim((-100, 100))
axs[1].set_xticks(np.arange(-100,100,20))
plt.tight_layout()
plt.show()