Tutorials

FFT Signal Analysis with SciPy

Every periodic signal — a heartbeat, a sound, a vibration — can be described as a sum of sine waves at different frequencies. The Fast Fourier Transform (FFT) computes this decomposition efficiently, converting a time-domain signal into its frequency-domain representation. You'd use it to identify the dominant frequencies in a sensor reading, filter out noise at specific frequencies, detect anomalies in vibration data, or analyze audio. For smoothing and peak detection in the time domain, see [signal filtering and smoothing](/tutorials/signal-filtering-and-smoothing-with-scipy) and [finding peaks with SciPy](/tutorials/find-peaks-with-scipy).

### Computing the Frequency Spectrum

`scipy.fft.fft` converts a signal array into complex amplitudes at each frequency. `fftfreq` gives the corresponding frequency values in Hz. The magnitude `|FFT|` tells you how much of each frequency is present in the signal.

from scipy.fft import fft, fftfreq
import numpy as np
import matplotlib.pyplot as plt

fs = 500               # sampling rate (samples per second)
t = np.linspace(0, 1, fs, endpoint=False)
signal = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t)

freqs = fftfreq(fs, d=1/fs)
spectrum = np.abs(fft(signal))

pos = freqs > 0  # keep only positive frequencies
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot(t[:80], signal[:80], color="steelblue", linewidth=1.5)
axes[0].set_title("Signal (first 0.16 s)")
axes[0].set_xlabel("Time (s)")
axes[0].set_ylabel("Amplitude")
axes[1].plot(freqs[pos], spectrum[pos], color="crimson", linewidth=1.5)
axes[1].set_title("Frequency Spectrum")
axes[1].set_xlabel("Frequency (Hz)")
axes[1].set_ylabel("|FFT|")
plt.tight_layout()
plt.show()
- `fftfreq(fs, d=1/fs)` produces frequency bins from 0 to fs/2 (positive) and −fs/2 to 0 (negative mirror); `d=1/fs` sets the sample spacing in seconds.
- The two peaks in the spectrum appear at exactly 50 Hz and 120 Hz — the frequencies we built into the signal.
- The height of each peak is proportional to the amplitude of that frequency component: the 50 Hz component (amplitude 1) is twice as tall as the 120 Hz component (amplitude 0.5).

### Finding Dominant Frequencies

Once you have the spectrum, you can automatically identify the top frequencies by sorting the magnitudes. Normalizing by `N/2` converts the raw FFT magnitude into the actual signal amplitude.

from scipy.fft import fft, fftfreq
import numpy as np

fs = 500
t = np.linspace(0, 1, fs, endpoint=False)
signal = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t)

freqs = fftfreq(fs, d=1/fs)
amplitudes = np.abs(fft(signal)) / (fs / 2)  # normalize to physical amplitude

pos = freqs > 0
top_idx = np.argsort(amplitudes[pos])[-3:][::-1]
print("Top 3 frequency components:")
for i in top_idx:
    print(f"  {freqs[pos][i]:.1f} Hz  amplitude ≈ {amplitudes[pos][i]:.3f}")
Top 3 frequency components:
  50.0 Hz  amplitude ≈ 1.000
  120.0 Hz  amplitude ≈ 0.500
  119.0 Hz  amplitude ≈ 0.000
- Dividing by `fs/2` normalizes the FFT output so that a pure sine wave of amplitude A shows up as A in the spectrum.
- `np.argsort(...)[-3:][::-1]` finds the indices of the three largest values in descending order — a compact way to rank components.
- In a real signal with noise, many small peaks appear; typically you'd threshold by amplitude to focus on the significant ones.

### FFT-Based Filtering

You can filter in the frequency domain by zeroing out unwanted frequencies in the spectrum and then converting back with the inverse FFT (`ifft`). This is equivalent to a perfect brick-wall filter: all frequencies outside the passband are removed exactly.

from scipy.fft import fft, ifft, fftfreq
import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)
fs = 500
t = np.linspace(0, 1, fs, endpoint=False)
signal = np.sin(2 * np.pi * 5 * t) + 0.6 * rng.standard_normal(fs)

spectrum = fft(signal)
freqs = fftfreq(fs, d=1/fs)

# Remove all frequencies above 15 Hz
filtered_spectrum = spectrum.copy()
filtered_spectrum[np.abs(freqs) > 15] = 0
filtered = np.real(ifft(filtered_spectrum))

plt.figure(figsize=(9, 3.5))
plt.plot(t, signal, color="steelblue", alpha=0.5, linewidth=1, label="Noisy signal (5 Hz)")
plt.plot(t, filtered, color="crimson", linewidth=2, label="FFT low-pass (cutoff = 15 Hz)")
plt.title("FFT-Based Low-Pass Filtering")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.legend()
plt.tight_layout()
plt.show()
- Setting `filtered_spectrum[np.abs(freqs) > 15] = 0` zeros every frequency component above 15 Hz, in both the positive and negative halves of the spectrum.
- `np.real(ifft(...))` converts back to the time domain; the small imaginary part from floating-point rounding is discarded.
- Brick-wall FFT filtering can introduce ringing (Gibbs phenomenon) at sharp transitions. For smooth filters see [signal filtering and smoothing](/tutorials/signal-filtering-and-smoothing-with-scipy).

### Spectrogram: Frequency Content Over Time

A spectrogram shows how the frequency content of a signal changes over time by computing a short-time FFT in a sliding window. `scipy.signal.spectrogram` handles this directly and returns the result ready to plot.

from scipy.signal import spectrogram
import numpy as np
import matplotlib.pyplot as plt

fs = 1000
t = np.linspace(0, 4, 4 * fs, endpoint=False)
# Frequency sweeps from 10 Hz to 200 Hz (chirp signal)
freq_inst = 10 + 47.5 * t
signal = np.sin(2 * np.pi * np.cumsum(freq_inst) / fs)

f, t_spec, Sxx = spectrogram(signal, fs=fs, nperseg=256)

plt.figure(figsize=(9, 4))
plt.pcolormesh(t_spec, f[:60], 10 * np.log10(Sxx[:60] + 1e-12),
               shading="gouraud", cmap="viridis")
plt.colorbar(label="Power (dB)")
plt.title("Spectrogram of a Chirp Signal (10–200 Hz)")
plt.xlabel("Time (s)")
plt.ylabel("Frequency (Hz)")
plt.tight_layout()
plt.show()
- `np.cumsum(freq_inst) / fs` integrates the instantaneous frequency to get the phase, producing a chirp that sweeps upward in frequency over 4 seconds.
- `spectrogram` divides the signal into overlapping windows of `nperseg=256` samples and FFTs each one; `f`, `t_spec`, and `Sxx` are the frequency bins, time bins, and power values.
- `10 * np.log10(Sxx + 1e-12)` converts power to decibels; the `+ 1e-12` prevents `log10(0)`. Slicing `[:60]` limits the y-axis to the lowest 60 frequency bins.

The FFT is one of the most widely used algorithms in science and engineering. Once you have the spectrum, you can [find peaks](/tutorials/find-peaks-with-scipy) in it to identify the dominant frequency components automatically.