Pythonでディジタルフィルタの周波数特性を求める(scipy.signal.freqz)

Python でディジタルフィルタの周波数特性（振幅特性・位相特性）を求めるには scipy.signal.freqz 関数を使います。主な使用例を説明したいと思います。

パラメータと返り値

scipy.signal.freqz の主なパラメータと返り値は以下です。

主な使用例

周波数特性を求める

バターワースローパスフィルタ（サンプリング周波数 48000Hz、カットオフ周波数 1000Hz、4 次）を設計し、freqz で振幅特性と位相特性を求めるソースコードは以下です。

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

fs = 48000       # サンプリング周波数 [Hz]
fc = 1000        # カットオフ周波数 [Hz]
order = 4        # フィルタ次数

b, a = sg.butter(order, fc, btype="low", fs=fs)
w, h = sg.freqz(b, a, worN=2048, fs=fs)

mag_db = 20 * np.log10(np.abs(h))    # 振幅特性[dB]
phase = np.unwrap(np.angle(h))       # 位相特性[rad]（アンラップ済み）

# 振幅特性と位相特性をプロット
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 9))

ax1.plot(w, mag_db)
ax1.axhline(-3, color="gray", linestyle="--")
ax1.axvline(fc, color="gray", linestyle="--")
ax1.set_xlim(0, fs / 2)
ax1.set_ylim(-80, 5)
ax1.set_xlabel("Frequency [Hz]")
ax1.set_ylabel("Magnitude [dB]")
ax1.grid(True)

ax2.plot(w, phase, color="tab:orange")
ax2.set_xlim(0, fs / 2)
ax2.set_xlabel("Frequency [Hz]")
ax2.set_ylabel("Phase [rad]")
ax2.grid(True)

fig.tight_layout()
fig.savefig("freqz_basic.png")

freqz に b と a を与えたとき、以下のディジタルフィルタの周波数特性 \( H(e^{j\omega}) \) を求めます。

\[ H(e^{j\omega}) = \dfrac{\sum_{k=0}^{M} b[k]\, e^{-j\omega k}}{\sum_{k=0}^{N} a[k]\, e^{-j\omega k}} \]

振幅特性は複素数 h の絶対値を dB 変換したもの、位相特性は h の偏角として求められます。プロットした結果は以下です。

グラフを見ると、カットオフ周波数 1000Hz のところでちょうど -3dB になっており、設計どおりの特性が得られていることが確認できます。

fsを指定してHz単位で結果を得る

fs 引数を省略すると w は rad/sample 単位（0〜π）で返ってきますが、fs にサンプリング周波数を渡すと w が直接 Hz 単位（0〜fs/2）で返ってきます。

b, a = sg.butter(4, 1000, btype="low", fs=48000)

w1, h1 = sg.freqz(b, a, worN=5)            # fsを指定しない場合
w2, h2 = sg.freqz(b, a, worN=5, fs=48000)  # fsを指定した場合

print(w1)
# [0.         0.62831853 1.25663706 1.88495559 2.51327412]
print(w2)
# [    0.  4800.  9600. 14400. 19200.]

fs を指定しない場合は rad/sample、指定した場合は Hz で得られていることが分かります。fs を指定しておくと正規化の計算をせずにそのままグラフの横軸として使えるので便利です。

worNに周波数を直接指定する

worN には評価点数だけでなく、調べたい周波数の配列を直接渡すこともできます。

freqs = np.array([500, 1000, 2000, 5000])    # 特性を知りたい周波数[Hz]
w, h = sg.freqz(b, a, worN=freqs, fs=48000)

mag_db = 20 * np.log10(np.abs(h))
print(mag_db)
# [ -0.02  -3.01 -24.25 -57.14]

500Hz ではほぼ減衰なし、カットオフの 1000Hz で -3dB、2000Hz・5000Hz と離れるほど大きく減衰していることが確認できます。全帯域のグラフを描かなくても、特定の周波数の特性だけをピンポイントで知りたいときに便利です。

フィルタ次数を変えて比較する

同じカットオフ周波数のまま次数を 2 次・4 次・8 次と変えて振幅特性を比較すると、次数によるロールオフ（減衰の急峻さ）の違いを確認できます。

for order in [2, 4, 8]:
    b, a = sg.butter(order, fc, btype="low", fs=fs)
    w, h = sg.freqz(b, a, worN=2048, fs=fs)
    mag_db = 20 * np.log10(np.abs(h))

おまけ

wholeパラメータ

whole=True を指定すると、0〜π ではなく0〜2π の全域で周波数応答を評価します。通常の実数係数フィルタでは 0〜π 側と π〜2π 側が対称になるため使う場面は少ないですが、複素数係数のフィルタを扱う場合などに使用します。

その他のパラメータ

● include_nyquist
whole=False かつ worN が整数のとき、True にするとナイキスト周波数（π または fs/2）を評価点に含めます。

● plot
結果をプロットする関数を渡せるパラメータですが、非推奨（deprecated）扱いになっているため、matplotlib で別途プロットするのがおすすめです。

おわりに

本記事では、ディジタルフィルタの周波数特性を求める関数 scipy.signal.freqz について紹介しました。フィルタ係数から特性を直接可視化できるので、設計したフィルタが意図どおりの特性になっているかを確かめる際にとても便利です。

■参考文献
[1] The SciPy community. “scipy.signal.freqz — SciPy Reference Guide”. (参照 2026-08-25). https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.freqz.html