Perfect FFT

readme.md

This is just a note to self, for remembering the little details about NumPy's FFT implementation.

• To get the FFT bins to line up perfectly with a single frequency bin, without any "skirts" or spectral leakage, you need to make a perfect cycle, where the next sample after this chunk lines up with the first. (In other words, the first and last samples should not be the same.)
• To get a sinusoid of amplitude 1 (= 0 dBFS = -3 dBov) to produce 2 complex exponentials of amplitude 0.5, you need to divide the fft() results by the number of samples.
• The fft() output is from 0 Hz to Nyquist frequency to sampling rate. To plot the spectrum from negative Nyquist frequency to positive Nyquist frequency, with 0 in the center, use fftshift() on both the freqs and ampl variables. You can also just use fftfreq() to generate a horizontal axis for plotting, but the plot will have extraneous lines on it.
  • fftshift() shifts the DC component from the bottom to the center of the spectrum.
  • ifftshift() shifts the DC component from the center to the bottom of the spectrum.
  • (They do the same thing for even-length, but not for odd-length.)
• If the result of an IFFT has some complex residue, use real() to get rid of it, not abs().

Perfect_FFT.py

from numpy import linspace, cos, pi, absolute
from numpy.fft import fft, fftfreq, fftshift
import matplotlib.pyplot as plt

# Sampling rate
fs = 64  # Hz

# Time is from 0 to 1 seconds, but leave off the endpoint, so
# that 1.0 seconds is the first sample of the *next* chunk
length = 1  # second
N = fs * length
t = linspace(0, length, num=N, endpoint=False)

# Generate a sinusoid at frequency f
f = 8  # Hz
a = cos(2 * pi * f * t)

# Plot signal, showing how endpoints wrap from one chunk to the next
plt.subplot(2, 1, 1)
plt.plot(t, a, '.-')
plt.plot(1, 1, 'r.')  # first sample of next chunk
plt.margins(0.1, 0.1)
plt.xlabel('Time [s]')

# Use FFT to get the amplitude of the spectrum
ampl = 1/N * absolute(fft(a))

# FFT frequency bins
freqs = fftfreq(N, 1/fs)

# Plot shifted data on a shifted axis
plt.subplot(2, 1, 2)
plt.stem(fftshift(freqs), fftshift(ampl))
plt.margins(0.1, 0.1)
plt.xlabel('Frequency [Hz]')
plt.tight_layout()
plt.show()

Perfect_RFFT.py

from numpy import linspace, cos, pi, absolute
from numpy.fft import rfft, rfftfreq
import matplotlib.pyplot as plt

# Sampling rate
fs = 64  # Hz

# Time is from 0 to 1 seconds, but leave off the endpoint, so
# that 1.0 seconds is the first sample of the *next* chunk
length = 1  # second
N = fs * length
t = linspace(0, length, num=N, endpoint=False)

# Generate a sinusoid at frequency f
f = 8  # Hz
a = cos(2 * pi * f * t)

# Plot signal, showing how endpoints wrap from one chunk to the next
plt.subplot(2, 1, 1)
plt.plot(t, a, '.-')
plt.plot(1, 1, 'r.')  # first sample of next chunk
plt.margins(0.1, 0.1)
plt.xlabel('Time [s]')

# Use RFFT to get the amplitude of the one-sided spectrum
ampl = 1/N * absolute(rfft(a))

# RFFT frequency bins
freqs = rfftfreq(N, 1/fs)

# Plot spectrum
plt.subplot(2, 1, 2)
plt.stem(freqs, ampl)
plt.margins(0.1, 0.1)
plt.xlabel('Frequency [Hz]')
plt.tight_layout()
plt.show()

signal and spectrum.png

@joegle Thanks, I added it.