{{ message }}

Instantly share code, notes, and snippets.

# endolith/Perfect_FFT.py

Last active Sep 16, 2020
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 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()`.
 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]') plt.tight_layout() plt.show()
 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]') plt.tight_layout() plt.show() ### tysonggraham commented Feb 2, 2016

 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 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.

### joegle commented Sep 30, 2017

 Thanks for sharing this On my version of matplotlib I have to put a final `plt.show()` display the window with the plot.