Fast Fourier Transform (FFT) in Python with NumPy

A fast Fourier transform in Python takes one NumPy call: np.fft.rfft(signal) turns a signal into its frequency components, and np.fft.rfftfreq() tells you which frequency each value belongs to:

import numpy as np
spectrum = np.fft.rfft(signal)
freqs = np.fft.rfftfreq(len(signal), d=1 / sample_rate)

All the numbers and the plot below come from Python 3.12.5, NumPy 2.5.3 and SciPy 1.18.1, run on a synthetic signal so you can see exactly what goes in and what comes out.

We’ll find the hidden frequencies in a noisy signal, plot the spectrum, filter out noise, spot a weekly cycle in daily website traffic, and see why it’s called fast.

Finding frequencies with numpy.fft

The signal you’ll analyze mixes a 50 Hz wave, a weaker 120 Hz wave and a lot of random noise. Looking at it, you’d never guess what’s inside; the FFT pulls both frequencies straight out:

import numpy as np

fs = 1000                                  # samples per second
t = np.arange(0, 1, 1 / fs)                # one second of data
rng = np.random.default_rng(0)
signal = 1.0 * np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t) + rng.normal(0, 0.8, t.size)

spectrum = np.fft.rfft(signal)                    # FFT for real-valued input
freqs = np.fft.rfftfreq(signal.size, d=1 / fs)    # the frequency of each bin, in Hz
amplitude = 2 * np.abs(spectrum) / signal.size    # scale back to the signal's units

top = np.argsort(amplitude)[-2:][::-1]
for i in top:
    print(f"{freqs[i]:>6.1f} Hz  amplitude {amplitude[i]:.2f}")
print("bins:", spectrum.size, "| resolution:", freqs[1] - freqs[0], "Hz")

Output:

  50.0 Hz  amplitude 1.01
 120.0 Hz  amplitude 0.52
bins: 501 | resolution: 1.0 Hz

You’ll want rfft() for real-valued data like audio or sensor readings; it returns only the non-negative frequencies. Multiplying by 2 and dividing by the length scales the result back to the original amplitudes.

Your frequency resolution is the sample rate divided by the number of samples. One second of data gives 1 Hz bins; record longer to separate frequencies that are closer together.

Plotting the FFT spectrum

A spectrum plot lets you see the same thing at a glance. The time-domain signal looks like noise, but the frequency domain has two clear spikes:

import numpy as np

fs = 1000                                  # samples per second
t = np.arange(0, 1, 1 / fs)                # one second of data
rng = np.random.default_rng(0)
signal = 1.0 * np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t) + rng.normal(0, 0.8, t.size)

import matplotlib.pyplot as plt

amplitude = 2 * np.abs(np.fft.rfft(signal)) / signal.size
freqs = np.fft.rfftfreq(signal.size, d=1 / fs)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 5.5))
ax1.plot(t[:200], signal[:200], lw=0.8)
ax1.set_title("Noisy signal (first 0.2 s)")
ax1.set_xlabel("Time (s)")
ax2.plot(freqs, amplitude, lw=0.8)
ax2.set_title("FFT amplitude spectrum")
ax2.set_xlabel("Frequency (Hz)")
ax2.set_xlim(0, 250)
plt.tight_layout()
plt.show()

How do you filter noise with an FFT?

Zero out the weak frequency bins, then run the inverse transform, irfft(), to get back to a cleaned-up signal:

import numpy as np

fs = 1000                                  # samples per second
t = np.arange(0, 1, 1 / fs)                # one second of data
rng = np.random.default_rng(0)
signal = 1.0 * np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t) + rng.normal(0, 0.8, t.size)

clean = 1.0 * np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t)

spectrum = np.fft.rfft(signal)
freqs = np.fft.rfftfreq(signal.size, d=1 / fs)
amplitude = 2 * np.abs(spectrum) / signal.size

spectrum[amplitude < 0.25] = 0              # keep only the strong frequencies
filtered = np.fft.irfft(spectrum, n=signal.size)

rms = lambda x: np.sqrt(np.mean(x ** 2))
print(f"error before filtering: {rms(signal - clean):.3f}")
print(f"error after filtering:  {rms(filtered - clean):.3f}")
print("frequencies kept:", freqs[spectrum != 0])

Output:

error before filtering: 0.782
error after filtering:  0.032
frequencies kept: [ 50. 120.]

I use a hard cutoff like this for demos. For real audio or sensor filtering, SciPy’s filter design tools give smoother results without ringing artifacts.

Using FFT to find a weekly cycle in business data

FFT isn’t just for audio, and you can use it on business data too. A year of daily website visits has a weekday/weekend pattern buried under growth and noise. Remove the trend first, then look for the strongest cycle:

import numpy as np

rng = np.random.default_rng(3)
days = np.arange(364)                                   # 52 weeks of daily visits
weekly = 1200 * np.cos(2 * np.pi * days / 7)             # weekday/weekend swing
visits = 8000 + 4 * days + weekly + rng.normal(0, 500, days.size)

detrended = visits - np.polyval(np.polyfit(days, visits, 1), days)
amp = np.abs(np.fft.rfft(detrended))
freqs = np.fft.rfftfreq(days.size, d=1)                  # cycles per day
peak = np.argmax(amp[1:]) + 1                            # skip the zero-frequency bin
print(f"strongest cycle: every {1 / freqs[peak]:.1f} days")

Output:

strongest cycle: every 7.0 days

The trend line comes from np.polyfit. Removing it keeps the slow growth from swamping the lowest frequencies. For other ways to explore patterns in data, see principal component analysis in Python.

Why is it called the fast Fourier transform?

The direct formula compares every sample with every frequency, which is O(n²). The FFT reorganizes that work into O(n log n). On just 2,048 samples you’ll see a dramatic difference:

import timeit
import numpy as np
from scipy import fft as sfft

x = np.random.default_rng(0).normal(size=2048)

def dft(x):                                    # the direct O(n^2) formula
    n = np.arange(x.size)
    return np.exp(-2j * np.pi * n[:, None] * n / x.size) @ x

print("same result:", np.allclose(dft(x), np.fft.fft(x)), np.allclose(np.fft.fft(x), sfft.fft(x)))
t_dft = timeit.timeit(lambda: dft(x), number=3) / 3
t_np = timeit.timeit(lambda: np.fft.fft(x), number=200) / 200
t_sp = timeit.timeit(lambda: sfft.fft(x), number=200) / 200
print(f"direct DFT {t_dft * 1000:.1f} ms | numpy.fft {t_np * 1000:.3f} ms | scipy.fft {t_sp * 1000:.3f} ms")

Output:

same result: True True
direct DFT 168.7 ms | numpy.fft 0.022 ms | scipy.fft 0.012 ms

NumPy and SciPy give identical results. scipy.fft adds extras such as multithreading with workers= and more transform types, so I use it for heavier work.

numpy.fft functions you’ll use most

Function	What it does
np.fft.rfft(x)	FFT of real input (non-negative frequencies only)
np.fft.rfftfreq(n, d)	Frequency for each rfft bin, given sample spacing d
np.fft.irfft(X, n)	Inverse of rfft, back to the time domain
np.fft.fft(x)	Full complex FFT, positive and negative frequencies
np.fft.fftshift(X)	Move zero frequency to the center for plotting

Keep exploring NumPy and SciPy:

• SciPy interpolate in Python
• Monte Carlo simulation in Python
• Inverse of a matrix in Python
• Standard deviation in Python

Frequently asked questions