Instantly share code, notes, and snippets.

# endolith/Perfect_FFT.py

Last active June 16, 2023 08:25
Show Gist options
• Save endolith/236567 to your computer and use it in GitHub Desktop.
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()`.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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()
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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.