scipy.signal.freqz: Plot a Filter’s Frequency Response

scipy.signal.freqz computes the frequency response of a digital filter from its coefficients. It returns two arrays:

from scipy import signal

b, a = signal.butter(4, 0.2)
w, h = signal.freqz(b, a)

# w -> 512 frequencies, 0 to pi radians per sample
# h -> the COMPLEX response at each one

h is complex, which is the detail that trips most people up. You take abs(h) for magnitude and np.angle(h) for phase.

Every number and plot below came from running the code on Python 3.12.5, SciPy 1.18.1, NumPy 2.5.3, Matplotlib 3.11.2.

Plotting a frequency response with freqz

The standard plot is magnitude in decibels against frequency, and it’s three lines of work:

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

# a 4th-order Butterworth low-pass, cutoff at 0.2 of the Nyquist frequency
b, a = signal.butter(4, 0.2)

w, h = signal.freqz(b, a)          # w: frequencies, h: complex response

plt.figure(figsize=(8, 4.2))
plt.plot(w, 20 * np.log10(abs(h)), color="#0b6bcb", lw=2)
plt.axhline(-3, color="#c0392b", ls="--", lw=1, label="-3 dB")
plt.title("Frequency response of a 4th-order Butterworth filter")
plt.xlabel("Frequency (radians / sample)")
plt.ylabel("Magnitude (dB)")
plt.ylim(-80, 5)
plt.grid(alpha=.3); plt.legend(); plt.tight_layout()
plt.show()

20 * np.log10(abs(h)) converts magnitude to decibels. Without that conversion the interesting part of the curve is squashed against zero.

The -3 dB line matters because that is the conventional definition of a filter’s cutoff, not the point where the response reaches zero.

What freqz returns: the w and h arrays

Print them before you plot anything. It saves a lot of guessing:

import numpy as np
from scipy import signal

b, a = signal.butter(4, 0.2)
w, h = signal.freqz(b, a)

print("w shape:", w.shape, "| h shape:", h.shape)
print("h dtype:", h.dtype)               # complex, not float
print()
print("w starts at:", w[0])
print("w ends at  :", round(float(w[-1]), 6))
print("pi is      :", round(float(np.pi), 6), "<- the grid stops just short of it")
print()
print("magnitude at DC   :", round(float(abs(h[0])), 4))
print("phase at DC (rad) :", round(float(np.angle(h[0])), 4))

Output:

w shape: (512,) | h shape: (512,)
h dtype: complex128

w starts at: 0.0
w ends at  : 3.135457
pi is      : 3.141593 <- the grid stops just short of it

magnitude at DC   : 1.0
phase at DC (rad) : 0.0

Returned	Contains	Use
w	Frequencies, radians per sample by default	The x axis
h	Complex response, one value per frequency	abs(h) and np.angle(h)

Notice where the grid ends: 3.135457, not 3.141593. The frequency points are evenly spaced but the final endpoint is excluded, so you never quite reach π.

That rarely matters for a plot. It does matter if you’re indexing h[-1] expecting the response exactly at Nyquist.

The scipy freqz fs parameter: frequencies in hertz

Radians per sample is awkward to read. Pass fs and you’ll get the frequencies back in hertz instead:

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

b, a = signal.butter(4, 0.2)

w_rad, h1 = signal.freqz(b, a)                 # radians per sample, 0 to pi
w_hz, h2 = signal.freqz(b, a, fs=1000)         # hertz, 0 to fs/2

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
ax1.plot(w_rad, 20 * np.log10(abs(h1)), color="#0b6bcb", lw=2)
ax1.set_title("Default: radians / sample")
ax1.set_xlabel("Frequency (rad/sample)")

ax2.plot(w_hz, 20 * np.log10(abs(h2)), color="#0a7d32", lw=2)
ax2.axvline(100, color="#c0392b", ls="--", lw=1, label="100 Hz cutoff")
ax2.set_title("fs=1000: hertz")
ax2.set_xlabel("Frequency (Hz)")
ax2.legend()

for ax in (ax1, ax2):
    ax.set_ylabel("Magnitude (dB)"); ax.set_ylim(-80, 5); ax.grid(alpha=.3)
plt.tight_layout(); plt.show()

print("without fs, w ends at:", round(float(w_rad[-1]), 4), "rad/sample")
print("with fs=1000, w ends at:", round(float(w_hz[-1]), 4), "Hz")

Output:

without fs, w ends at: 3.1355 rad/sample
with fs=1000, w ends at: 499.0234 Hz

With fs=1000, the axis runs from 0 to just under 500 Hz. That upper limit is the Nyquist frequency, always half the sampling rate.

Design the filter with the same fs and you can specify the cutoff in hertz too, which is far easier to reason about than a fraction of Nyquist.

The scipy freqz worN parameter: an integer or an array

worN does two different jobs depending on what you hand it. An integer sets how many points to compute:

import numpy as np
from scipy import signal

b, a = signal.butter(4, 0.2)

# worN as an integer: how many points on the grid
for n in (8, 512, 2048):
    w, h = signal.freqz(b, a, worN=n)
    print(f"worN={n:<5} -> {len(w)} points")

print()
# worN as an ARRAY: evaluate at exactly the frequencies you name
wanted = [0, 50, 100, 200, 400]
w, h = signal.freqz(b, a, worN=wanted, fs=1000)

print(f"{'Hz':>6} {'|h|':>9} {'dB':>9}")
for f, mag in zip(w, abs(h)):
    print(f"{f:>6.0f} {mag:>9.4f} {20 * np.log10(mag):>9.2f}")
print()
print("100 Hz gives |h| = 0.7071, which is exactly -3 dB: the cutoff.")

Output:

worN=8     -> 8 points
worN=512   -> 512 points
worN=2048  -> 2048 points

    Hz       |h|        dB
     0    1.0000     -0.00
    50    0.9984     -0.01
   100    0.7071     -3.01
   200    0.0400    -27.97
   400    0.0001    -78.12

100 Hz gives |h| = 0.7071, which is exactly -3 dB: the cutoff.

• worN=2048 — a finer grid, useful for a sharp filter whose notch falls between the default points.
• worN=[0, 50, 100] — evaluate at those frequencies and nothing else.
• worN=8 — too coarse to plot, but handy when you just want a few numbers.

The array form is the one worth remembering. It answers questions like “how much does this filter attenuate 200 Hz” in one line, with no searching through a 512-point array.

The table above confirms the design: at 100 Hz the magnitude is 0.7071, which is -3.01 dB. That is the cutoff landing exactly where butter was told to put it.

Plotting the freqz phase response

The same h carries the phase. Take the angle, then unwrap it:

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

b, a = signal.butter(4, 0.2)
w, h = signal.freqz(b, a, fs=1000)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6), sharex=True)

ax1.plot(w, 20 * np.log10(abs(h)), color="#0b6bcb", lw=2)
ax1.set_ylabel("Magnitude (dB)")
ax1.set_title("Magnitude and phase response")
ax1.set_ylim(-80, 5); ax1.grid(alpha=.3)

phase = np.unwrap(np.angle(h))          # unwrap removes the jumps at +/- pi
ax2.plot(w, np.degrees(phase), color="#8e44ad", lw=2)
ax2.set_ylabel("Phase (degrees)")
ax2.set_xlabel("Frequency (Hz)")
ax2.grid(alpha=.3)

plt.tight_layout(); plt.show()

raw = np.angle(h)
print("raw phase jumps between", round(float(raw.min()), 2), "and", round(float(raw.max()), 2), "radians")
print("unwrapped runs smoothly from", round(float(phase[0]), 2), "to", round(float(phase[-1]), 2))

Output:

raw phase jumps between -3.13 and 3.12 radians
unwrapped runs smoothly from 0.0 to -6.28

np.angle(h) returns values confined to -π..π, so the real phase curve arrives chopped into pieces with vertical jumps.

np.unwrap adds the missing multiples of 2π back and gives you a continuous line. Always unwrap before plotting phase.

Phase tells you the delay the filter introduces at each frequency. If that delay matters for your application, look at Butterworth filter design or at filtfilt, which cancels phase distortion by filtering twice.

Why does scipy freqz fail on second-order sections?

Because it isn’t the right function for them. If you designed the filter with output="sos", you need sosfreqz:

import numpy as np
from scipy import signal

sos = signal.butter(4, 0.2, output="sos")
print("sos shape:", sos.shape, "(sections x 6 coefficients)")

try:
    w, h = signal.freqz(sos)
except ValueError as err:
    print("freqz(sos) raised ValueError:", err)

# sosfreqz is the function for second-order sections
w, h = signal.sosfreqz(sos, fs=1000)
print()
print("sosfreqz worked:", w.shape, h.shape)

b, a = signal.butter(4, 0.2)
_, h_ba = signal.freqz(b, a, fs=1000)
print("largest difference vs freqz(b, a):", float(np.max(np.abs(abs(h) - abs(h_ba)))))

Output:

sos shape: (2, 6) (sections x 6 coefficients)
freqz(sos) raised ValueError: Array shapes are incompatible for broadcasting.

sosfreqz worked: (512,) (512,)
largest difference vs freqz(b, a): 4.884981308350689e-15

The error is Array shapes are incompatible for broadcasting, which is not an obvious clue. It happens because freqz reads the first argument as numerator coefficients and an SOS array is two-dimensional.

sosfreqz takes the same keyword arguments, including fs and worN. The output matches freqz(b, a) to about 5e-15, so the two are interchangeable in accuracy.

Prefer SOS for filters above roughly order 8. Transfer-function coefficients become numerically unstable at high orders, and SOS avoids that.

Comparing filter orders with freqz

Because freqz just returns arrays, comparing designs is an ordinary loop:

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

plt.figure(figsize=(8, 4.5))
for order, colour in ((2, "#0b6bcb"), (4, "#0a7d32"), (8, "#c0392b")):
    b, a = signal.butter(order, 100, fs=1000)
    w, h = signal.freqz(b, a, fs=1000)
    plt.plot(w, 20 * np.log10(abs(h)), color=colour, lw=2, label=f"order {order}")

plt.axvline(100, color="#555", ls=":", lw=1)
plt.axhline