Skip to content
Create a gist now

Instantly share code, notes, and snippets.

Perfect FFT

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, 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 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.
  • If the result of an IFFT has some complex residue, use real() to get rid of it, not abs().
from __future__ import division
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]')
from __future__ import division
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]')
@maxwillzq

The function stem does not work? It is from which module? Thanks

@endolith
Owner
endolith commented Mar 2, 2012

Sorry:

from matplotlib.pyplot import stem

or just

from pylab import *

which is what I do

@maxwillzq

Thanks a lot, endolith.

@tysonggraham

Is there an equation to get the frequency in hz from the sample output? I thought is was (fs * k/N) but in implementing that I don't get 8 with the numbers here. So that would be (64*56)/64 and that just gives 56...
So to recap:
N = sample width
fs = sample frequency
k = index that we get our spike on.
I am definitely confused here so sorry if I am unclear. Let me know if you need more info concerning the question and thank you very much for your work.

@endolith
Owner
endolith commented May 20, 2016 edited

sorry @tysonggraham this doesn't notify me when people leave comments. scipy has functions fftfreq and rfftfreq to give you the frequencies in Hz. You can look at their code to see how.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.