はじめに
**ボード線図(Bode plot)**は、線形時不変システムの周波数応答 \(H(j\omega)\) を 振幅[dB] と 位相[deg] に分け、横軸を 対数スケール(log scale)の周波数 で表した2段グラフです。Hendrik Bodeが1930年代に制御工学のために体系化したこの可視化手法は、現在ではフィルタ設計、伝達関数解析、サーボ系の安定性評価など、信号処理と制御の現場で標準ツールになっています。
特にフィルタ設計では、 バターワースフィルタ の「最大平坦特性」、 チェビシェフフィルタ の「等リップル」、 ハイパスフィルタ や ローパスフィルタ の「ロールオフ率」など、すべての特性がボード線図上で視覚化されます。本記事では、ボード線図の数学的基礎・読み方・Python実装を整理し、各種フィルタの周波数応答を横断的に理解できるハブとして使えるようにまとめます。
数学的基礎
伝達関数と周波数応答
連続時間 LTI システムの伝達関数を \(H(s)\) とします。\(s = j\omega\) を代入して得られる 周波数応答 は、複素数値関数として振幅と位相を持ちます。
\[H(j\omega) = |H(j\omega)| \, e^{j\angle H(j\omega)} \tag{1}\]ボード線図はこの \(H(j\omega)\) を 2 つのプロットに分解します。
- 振幅プロット: \(20\log_{10}|H(j\omega)|\) を縦軸に、\(\log_{10}\omega\) を横軸にとる
- 位相プロット: \(\angle H(j\omega)\) (度)を縦軸に、\(\log_{10}\omega\) を横軸にとる
デシベル定義
振幅をデシベル(dB)で表す理由は、乗算が加算に変わることにあります。\(H_1\) と \(H_2\) を直列接続したシステムの振幅は
\[20\log_{10}|H_1 H_2| = 20\log_{10}|H_1| + 20\log_{10}|H_2| \tag{2}\]となり、ボード線図上では 個別のプロットの和 として描けます。これにより、複数次の極・零点を持つフィルタも 1 次系の重ね合わせで近似的に読み解けるのです。
位相遅延と群遅延
位相プロットから直接読める量として、位相遅延 と 群遅延 が重要です。
\[\tau_p(\omega) = -\frac{\angle H(j\omega)}{\omega}, \quad \tau_g(\omega) = -\frac{d\angle H(j\omega)}{d\omega} \tag{3}\]群遅延 \(\tau_g\) は位相プロットの 傾き に直結し、波形歪みの大きさを表します。FIR フィルタの線形位相は群遅延が定数になることを意味し、ボード線図上では位相が周波数の一次関数で減少します。
1 次・2 次系のボード線図
1 次ローパス
\[H(s) = \frac{1}{1 + s/\omega_c} \tag{4}\]の振幅は
\[|H(j\omega)| = \frac{1}{\sqrt{1 + (\omega/\omega_c)^2}} \tag{5}\]であり、漸近線で読むと
- \(\omega \ll \omega_c\) : \(|H|_{dB} \approx 0\) dB(通過域)
- \(\omega = \omega_c\) : \(|H|_{dB} = -3\) dB(カットオフ点)
- \(\omega \gg \omega_c\) : \(|H|_{dB} \approx -20\log_{10}(\omega/\omega_c)\) dB(-20 dB/dec のロールオフ)
位相は \(\omega_c\) で \(-45^\circ\) 、高域漸近で \(-90^\circ\) になります。1 次ハイパスはこの 鏡像(低域で \(-20\) dB/dec の立ち上がり、\(+90^\circ \to 0^\circ\) への位相回転)として読めます。
2 次系の共振
2 次系
\[H(s) = \frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2} \tag{6}\]の振幅プロットは、減衰係数 \(\zeta < 1/\sqrt{2}\) で 共振ピーク を持ちます。ピーク値は
\[M_p = \frac{1}{2\zeta\sqrt{1-\zeta^2}} \tag{7}\]であり、\(\zeta\) が小さいほど鋭いピークが立ちます。チェビシェフフィルタの通過域リップルや楕円フィルタの両域リップルは、この 2 次系の共振が複数組み合わさった結果と見なせます。高域漸近のロールオフは -40 dB/dec で、1 次系の倍の急峻さです。
典型的フィルタとボード線図の対応
| フィルタ | 通過域の振幅特性 | ロールオフ | 位相特性 |
|---|---|---|---|
| バターワース \(N\) 次 | 最大平坦 | \(-20N\) dB/dec | 比較的滑らか |
| チェビシェフ I 型 | 等リップル | \(-20N\) dB/dec | 通過域端で急峻に変化 |
| 楕円(Cauer) | 両域リップル | 最急峻 | 最も急峻 |
| ハイパス | \(f > f_c\) で平坦 | 低域 \(-20N\) /dec | \(+90N^\circ \to 0^\circ\) |
| バンドパス | 帯域内で平坦 | 両側 \(-20N\) /dec | 帯域中心で 0° |
| ノッチ | \(f_0\) で深いディップ | 急峻 | \(f_0\) で \(\pm 90^\circ\) 回転 |
このようにボード線図は 「フィルタの種類を一目で見分ける視覚的辞書」 として機能します。次数 \(N\) はロールオフ率から、減衰特性は通過域・阻止域のリップル形状から、位相歪みは位相プロットの曲がりから読み取れます。
Python による作成
scipy.signal.bode による標準実装
scipy.signal.bode は連続時間伝達関数からボード線図のデータを直接生成できます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# --- 1次ローパス: H(s) = 1 / (1 + s/wc) ---
wc = 2 * np.pi * 100 # カットオフ角周波数 [rad/s] (100 Hz)
num = [1.0]
den = [1.0 / wc, 1.0]
system = signal.TransferFunction(num, den)
# bode は 振幅[dB] と 位相[deg] を返す
w, mag, phase = signal.bode(system, w=np.logspace(0, 4, 1000))
f = w / (2 * np.pi) # Hz に変換
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7), sharex=True)
ax1.semilogx(f, mag, linewidth=2)
ax1.axhline(-3, color='gray', linestyle=':', label='-3 dB')
ax1.axvline(wc / (2 * np.pi), color='r', linestyle='--', label=f'$f_c$ = {wc/(2*np.pi):.0f} Hz')
ax1.set_ylabel('Magnitude [dB]')
ax1.set_title('Bode Plot - 1st-order Lowpass')
ax1.grid(True, which='both', alpha=0.3)
ax1.legend()
ax2.semilogx(f, phase, linewidth=2, color='orange')
ax2.axvline(wc / (2 * np.pi), color='r', linestyle='--')
ax2.axhline(-45, color='gray', linestyle=':', label='-45°')
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('Phase [degrees]')
ax2.grid(True, which='both', alpha=0.3)
ax2.legend()
plt.tight_layout()
plt.show()
semilogx で横軸を対数化することがポイントです。リニア軸では低域の構造が潰れてしまい、ボード線図本来の 「数十年(decade)にわたる広帯域挙動を一度に俯瞰する」 という長所が失われます。
デジタルフィルタ(IIR)のボード線図
バターワース
や
チェビシェフ
のような IIR デジタルフィルタには scipy.signal.freqz を使います。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
fs = 1000.0 # サンプリング周波数 [Hz]
fc = 100.0 # カットオフ周波数 [Hz]
orders = [2, 4, 8]
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
for N in orders:
sos = signal.butter(N, fc, btype='low', fs=fs, output='sos')
w, h = signal.sosfreqz(sos, worN=4096, fs=fs)
# 振幅プロット(対数周波数 + dB)
ax1.semilogx(w, 20 * np.log10(np.abs(h) + 1e-12), label=f'N={N}', linewidth=2)
# 位相プロット(unwrap で 360° 跳びを除去)
ax2.semilogx(w, np.degrees(np.unwrap(np.angle(h))), label=f'N={N}', linewidth=2)
ax1.axhline(-3, color='gray', linestyle=':', label='-3 dB')
ax1.axvline(fc, color='r', linestyle='--', alpha=0.5)
ax1.set_ylabel('Magnitude [dB]')
ax1.set_title('Bode Plot - Butterworth Lowpass (digital)')
ax1.set_ylim(-100, 5)
ax1.grid(True, which='both', alpha=0.3)
ax1.legend()
ax2.axvline(fc, color='r', linestyle='--', alpha=0.5)
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('Phase [degrees]')
ax2.set_xlim(1, fs / 2)
ax2.grid(True, which='both', alpha=0.3)
ax2.legend()
plt.tight_layout()
plt.show()
このプロットからは、次数 \(N\) を 1 段上げるごとに高域のロールオフが 20 dB/dec ずつ急峻になる こと、位相回転が \(-90^\circ\) ずつ深まることが確認できます。これが「次数から特性を予測する」感覚を養う最良のトレーニングです。
手動計算による複数フィルタの比較
bode / freqz を使わず手動で計算すると、各フィルタの違いを並列に可視化できます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
fs = 1000.0
fc = 100.0
N = 4
f = np.logspace(0, np.log10(fs / 2), 2000)
w = 2 * np.pi * f
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=True)
filters = {
'Butterworth': signal.butter(N, fc, btype='low', fs=fs, output='ba'),
'Chebyshev I': signal.cheby1(N, 1.0, fc, btype='low', fs=fs, output='ba'),
'Chebyshev II': signal.cheby2(N, 40.0, fc, btype='low', fs=fs, output='ba'),
'Elliptic': signal.ellip(N, 1.0, 40.0, fc, btype='low', fs=fs, output='ba'),
}
for name, (b, a) in filters.items():
_, h = signal.freqz(b, a, worN=f, fs=fs)
ax1.semilogx(f, 20 * np.log10(np.abs(h) + 1e-12), label=name, linewidth=2)
ax2.semilogx(f, np.degrees(np.unwrap(np.angle(h))), label=name, linewidth=2)
ax1.axvline(fc, color='r', linestyle='--', alpha=0.5, label=f'$f_c$ = {fc} Hz')
ax1.axhline(-3, color='gray', linestyle=':', alpha=0.5)
ax1.set_ylabel('Magnitude [dB]')
ax1.set_title(f'Bode Plot Comparison (N={N})')
ax1.set_ylim(-80, 5)
ax1.grid(True, which='both', alpha=0.3)
ax1.legend()
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('Phase [degrees]')
ax2.grid(True, which='both', alpha=0.3)
ax2.legend()
plt.tight_layout()
plt.savefig("bode_plot_comparison.png", dpi=150)
plt.show()
![バターワース・チェビシェフI型・チェビシェフII型・楕円フィルタ(いずれも4次、カットオフ100Hz)のボード線図比較。上段が振幅[dB]、下段が位相[deg]](/posts/20260521_bode_plot/bode_plot_comparison.png)
同じ次数・同じカットオフで並べると、バターワースの平坦な通過域、チェビシェフ I 型の通過域リップル、チェビシェフ II 型の阻止域リップル、楕円フィルタの両域リップルが ボード線図上で直接比較 できます。
設計指針
カットオフから次数を決める
通過域 \(f_p\) の最大減衰 \(A_p\) [dB] と阻止域 \(f_s\) の最小減衰 \(A_s\) [dB] が与えられたとき、バターワース次数は
\[N \geq \frac{\log_{10}\!\left(\dfrac{10^{A_s/10} - 1}{10^{A_p/10} - 1}\right)}{2\log_{10}(f_s / f_p)} \tag{8}\]で求まります。これは ボード線図上の「2 点」から傾きを逆算して必要なロールオフ率を決める 操作です。scipy.signal.buttord がこの計算を自動化します。
数値例: \(f_p = 80\)
Hz で \(A_p = 1\)
dB以内、\(f_s = 150\)
Hz で \(A_s = 40\)
dB以上の減衰を要求する場合、式(8)を手計算すると \(N \geq 8.40\)
となり、切り上げて \(N=9\)
が必要です。scipy.signal.buttord(80, 150, 1.0, 40.0, analog=True) を実行すると同じく \(N=9\)
(\(\omega_n=86.24\)
) を返し、式(8)の手計算と一致することを確認できます。一方、analog=False(デジタルフィルタ、\(f_s=1000\)
Hz でバイリニア変換を使用)で同じ条件を計算すると \(N=8\)
に減ります。これはバイリニア変換の周波数ワーピングにより高周波側の遷移帯域が実質的に広がるためで、アナログ設計の式(8)をそのままデジタルフィルタの次数見積もりに使うと過大に見積もる場合がある点に注意が必要です。
位相回りに注意
ロールオフ率を急峻にするほど位相回転が深くなり、群遅延の周波数依存性が増します。
- 音声処理・医療信号: 位相歪みを避けたい →
FIR で線形位相を確保
、または
filtfiltでゼロ位相化 - リアルタイム制御: 群遅延を最小化したい → ベッセルフィルタや低次バターワース
- 通過域の平坦性が最優先: バターワース
- 遷移帯域の急峻さが最優先: チェビシェフ・楕円
ボード線図の位相プロットは、この 「振幅性能と位相性能のトレードオフ」 を一望で示してくれる重要なツールです。
適応フィルタとの違い
ボード線図は 時不変な伝達関数を持つ固定フィルタ にのみ意味を持ちます。 LMS/RLS などの適応フィルタ は係数が時間で変化するため、ある瞬間のスナップショットとしてのみボード線図を描けます。逆に、ボード線図で見た「目標の周波数応答」を適応フィルタが学習で近づける、という関係になります。
FFT / PSD との関係
ボード線図はフィルタ自体の特性を表しますが、実信号 がフィルタを通った結果を見るには FFT や 窓関数 + PSD と組み合わせます。具体的には、入力スペクトル \(X(f)\) と出力スペクトル \(Y(f) = H(f) X(f)\) から \(|H(f)| = |Y(f)/X(f)|\) を推定すれば、実測ボード線図 が得られます。これがシステム同定の基本原理です。
極・零点のロールオフ:漸近線の厳密な導出
これまでは「\(-20\) dB/dec になる」「位相は \(-45°\) を通る」という結果だけを提示してきました。ここでは、なぜそうなるのかを式変形だけから導出し、複数の極・零点が本当に単純な足し算で重ね合わされるのかを数値で検証します。
1次極が -20 dB/dec を生む理由
式(5)の振幅 \(|H(j\omega)| = 1/\sqrt{1+(\omega/\omega_c)^2}\) をデシベルに直すと、
\[ |H(j\omega)|_{dB} = -20\log_{10}\sqrt{1+(\omega/\omega_c)^2} = -10\log_{10}\!\left[1+\left(\frac{\omega}{\omega_c}\right)^{2}\right] \]となります。高域 \(\omega \gg \omega_c\) では括弧内の \(1\) が無視できるため、
\[ |H(j\omega)|_{dB} \ \xrightarrow{\ \omega \gg \omega_c\ }\ -20\log_{10}\!\left(\frac{\omega}{\omega_c}\right) \tag{9} \]という \(\log_{10}\omega\) の一次関数に漸近します。\(\omega\) を10倍(1 decade)したときの変化量を計算すると、
\[ -20\log_{10}\!\left(\frac{10\omega}{\omega_c}\right) - \left[-20\log_{10}\!\left(\frac{\omega}{\omega_c}\right)\right] = -20\log_{10}10 = -20\ \text{dB} \tag{10} \]となり、周波数が1桁上がるごとに振幅が厳密に20 dB下がることが式変形だけから導けます。これは近似ではなく、\(\omega \gg \omega_c\) の極限で厳密に成り立つ漸近線の傾きです。
位相は \(\angle H(j\omega) = -\arctan(\omega/\omega_c)\) であり、\(\omega=\omega_c\) でちょうど \(-45°\) (\(\arctan 1 = 45°\) )、\(\omega\to\infty\) で \(-90°\) に漸近します。教科書でよく使われる「\(0.1\omega_c\) で \(0°\) 、\(10\omega_c\) で \(-90°\) 、その間を直線で結ぶ」という折れ線近似が実際の \(\arctan\) 曲線とどの程度ずれるかを数値的に確認します。
import numpy as np
from scipy import signal
wc = 10.0
sys_pole = signal.TransferFunction([1.0], [1.0 / wc, 1.0])
for wt in [0.1 * wc, wc, 10 * wc]:
_, mag, phase = signal.bode(sys_pole, w=np.array([wt]))
print(f"w={wt:6.1f}: mag={mag[0]:7.3f} dB, phase={phase[0]:7.2f} deg")
実行結果:
w= 1.0: mag= -0.043 dB, phase= -5.71 deg
w= 10.0: mag= -3.010 dB, phase= -45.00 deg
w= 100.0: mag=-20.043 dB, phase= -84.29 deg
折れ線近似は \(0.1\omega_c\) で \(0°\) 、\(10\omega_c\) で \(-90°\) と予測しますが、実際にはそれぞれ \(-5.71°\) 、\(-84.29°\) であり、decade の端でちょうど \(5.71°\) (\(=\arctan(0.1)\times180/\pi\) )の誤差が生じます。\(\omega_c\) ちょうどでは近似・厳密解ともに \(-45°\) で完全に一致します。
1次零点は鏡像
零点 \(H(s) = 1+s/\omega_c\) は極と分子・分母の関係が逆転するため、上記の議論の符号がそのまま反転します。
\[ |H(j\omega)|_{dB} \ \xrightarrow{\ \omega \gg \omega_c\ }\ +20\log_{10}\!\left(\frac{\omega}{\omega_c}\right) \tag{11} \]位相は \(+\arctan(\omega/\omega_c)\) で \(0° \to +90°\) に単調増加します。極が「振幅は下向きに20 dB/dec、位相は \(0°\) から \(-90°\) 」であるのに対し、零点は「振幅は上向きに20 dB/dec、位相は \(0°\) から \(+90°\) 」と、振幅・位相の両方が符号反転した鏡像関係にあります。
複数の極・零点は対数軸上でそのまま足し算できる
式(2)のデシベル加法性は、極・零点の個数によらず成り立ちます。カスケード接続された系全体の振幅は、各要素の \(dB\) 値の単純な和になるはずです。これを実際に数値で検証します。
import numpy as np
from scipy import signal
wc1, wc2, wcz = 5.0, 50.0, 500.0
sys1 = signal.TransferFunction([1.0], [1.0 / wc1, 1.0]) # 極 (wc1 = 5 rad/s)
sys2 = signal.TransferFunction([1.0], [1.0 / wc2, 1.0]) # 極 (wc2 = 50 rad/s)
sysz = signal.TransferFunction([1.0 / wcz, 1.0], [1.0]) # 零点 (wcz = 500 rad/s)
num_casc = np.polymul([1.0 / wcz, 1.0], [1.0])
den_casc = np.polymul([1.0 / wc1, 1.0], [1.0 / wc2, 1.0])
sys_casc = signal.TransferFunction(num_casc, den_casc) # 2極+1零点のカスケード全体
w = np.array([1000.0])
_, mag1, _ = signal.bode(sys1, w=w)
_, mag2, _ = signal.bode(sys2, w=w)
_, magz, _ = signal.bode(sysz, w=w)
_, mag_casc, _ = signal.bode(sys_casc, w=w)
print(f"pole1={mag1[0]:.4f} dB, pole2={mag2[0]:.4f} dB, zero={magz[0]:.4f} dB")
print(f"sum of individual dB = {mag1[0] + mag2[0] + magz[0]:.4f} dB")
print(f"cascade actual dB = {mag_casc[0]:.4f} dB")
実行結果(\(\omega=1000\) rad/s):
pole1=-46.0207 dB, pole2=-26.0314 dB, zero=6.9897 dB
sum of individual dB = -65.0625 dB
cascade actual dB = -65.0625 dB
個別要素の \(dB\) 値を単純に足した \(-65.0625\) dBは、実際にカスケード接続した伝達関数を直接計算した値と小数点以下4桁まで完全に一致します。2極1零点(正味の相対次数1)の高域ロールオフが1次極1本分(\(-20\) dB/dec)に収束することも、この足し算から直接読み取れます。
発展トピック:非最小位相系とゲイン・位相余裕
ここまでは「与えられた伝達関数のボード線図をどう読むか」が主題でした。ここからは、ボード線図を制御系の安定性解析に使う際に見落としやすい3つの落とし穴を、実際に scipy.signal で数値検証しながら見ていきます。
非最小位相系(RHPゼロ):振幅は同じでも位相はまったく違う
零点が左半平面(LHP, \(\mathrm{Re}(s)<0\) )にある系を最小位相系、右半平面(RHP, \(\mathrm{Re}(s)>0\) )にある系を非最小位相系と呼びます。同じ極配置に対して零点をLHPからRHPへ移しても、振幅 \(|H(j\omega)|\) はまったく変化しません。零点 \(s=-a\) (LHP)なら \(|-a+j\omega|=\sqrt{a^2+\omega^2}\) 、零点 \(s=+a\) (RHP)なら \(|a-j\omega|=\sqrt{a^2+\omega^2}\) と、絶対値がそもそも一致するためです。
\[ H_{\mathrm{mp}}(s) = \frac{s+1}{(s+2)(s+5)}, \qquad H_{\mathrm{nmp}}(s) = \frac{-s+1}{(s+2)(s+5)} \tag{12} \](\(H_{\mathrm{mp}}\) は零点 \(s=-1\) の最小位相系、\(H_{\mathrm{nmp}}\) は零点 \(s=+1\) の非最小位相系。極配置は共通。)
import numpy as np
from scipy import signal
num_mp = [1, 1] # 零点 s=-1 (LHP)
num_nmp = [-1, 1] # 零点 s=+1 (RHP)
den = [1, 7, 10] # (s+2)(s+5)
sys_mp = signal.TransferFunction(num_mp, den)
sys_nmp = signal.TransferFunction(num_nmp, den)
w = np.logspace(-2, 2, 2000)
w, mag_mp, phase_mp = signal.bode(sys_mp, w=w)
_, mag_nmp, phase_nmp = signal.bode(sys_nmp, w=w)
print("max|mag_mp - mag_nmp| =", np.max(np.abs(mag_mp - mag_nmp)))
print(f"phase_mp(w=100) = {phase_mp[-1]:.2f} deg")
print(f"phase_nmp(w=100) = {phase_nmp[-1]:.2f} deg")
実行結果:
max|mag_mp - mag_nmp| = 7.105427357601002e-15
phase_mp(w=100) = -86.56 deg
phase_nmp(w=100) = -265.42 deg
振幅の差は浮動小数点誤差レベル(\(7\times10^{-15}\) )で、事実上完全に一致します。一方で位相は \(\omega=100\) rad/s において、最小位相系が \(-86.56°\) (相対次数1なので\(-90°\) に漸近)であるのに対し、非最小位相系は \(-265.42°\) (\(-270°\) に漸近)と、\(180°\) 近くもの差が生じます。RHPゼロは高域で余分な \(-180°\) の位相遅れ(LHPゼロと比べて2倍の位相回転)を追加するためです。

この違いが実務で重要な理由: フィードバック制御において非最小位相ゼロは、制御帯域を零点周波数より大幅に低く抑えないと不安定化しやすいという制約を課します。振幅のボード線図だけを見ていると最小位相系と区別がつかないため、位相プロットを必ず確認することが非最小位相系を見逃さないための鉄則です。倒立振子やロケットのピッチ制御(操作直後に一瞬逆方向へ動く現象)が典型例として知られています。
ゲイン余裕と位相余裕:似て非なる2つの安定性指標
開ループ伝達関数 \(L(j\omega)\) に対して、安定余裕は2つの独立した量として定義されます。
\[ \mathrm{PM} = 180° + \angle L(j\omega_{gc}), \quad \text{where}\ |L(j\omega_{gc})| = 0\ \mathrm{dB} \] \[ \mathrm{GM} = -20\log_{10}|L(j\omega_{pc})|, \quad \text{where}\ \angle L(j\omega_{pc}) = -180° \]ここで \(\omega_{gc}\) はゲイン交差周波数(振幅が0dBを切る点)、\(\omega_{pc}\) は位相交差周波数(位相が\(-180°\) を切る点)です。両者は定義上異なる周波数で評価される独立な指標であり、片方が大きくても他方が小さいことは十分ありえます。
3次系の開ループ伝達関数 \(L(s)=K/[s(s+1)(s+5)]\) (積分器1つ + 実極2つ)で \(K=8\) とした場合を数値的に確認します。
import numpy as np
from scipy import signal
def find_margins(num, den, w):
sys = signal.TransferFunction(num, den)
w, mag, phase = signal.bode(sys, w=w)
gc_idx = np.where(np.diff(np.sign(mag)))[0][0]
w0, w1 = w[gc_idx], w[gc_idx + 1]
m0, m1 = mag[gc_idx], mag[gc_idx + 1]
wgc = w0 + (0 - m0) * (w1 - w0) / (m1 - m0)
p0, p1 = phase[gc_idx], phase[gc_idx + 1]
pm = 180 + (p0 + (wgc - w0) * (p1 - p0) / (w1 - w0))
pc_idx = np.where(np.diff(np.sign(phase + 180)))[0][0]
w0, w1 = w[pc_idx], w[pc_idx + 1]
p0, p1 = phase[pc_idx], phase[pc_idx + 1]
wpc = w0 + (-180 - p0) * (w1 - w0) / (p1 - p0)
m0, m1 = mag[pc_idx], mag[pc_idx + 1]
gm = -(m0 + (wpc - w0) * (m1 - m0) / (w1 - w0))
return wgc, pm, wpc, gm
K = 8.0
den = np.polymul([1, 0], np.polymul([1, 1], [1, 5]))
w = np.logspace(-2, 2, 400000)
wgc, pm, wpc, gm = find_margins([K], den, w)
print(f"wgc={wgc:.4f} rad/s, PM={pm:.2f} deg")
print(f"wpc={wpc:.4f} rad/s, GM={gm:.2f} dB")
実行結果:
wgc=1.0689 rad/s, PM=31.02 deg
wpc=2.2361 rad/s, GM=11.48 dB
実務では「\(\mathrm{GM}>6\) dB かつ \(\mathrm{PM}>45°\) 」がよく使われる目安ですが、この系はGM=11.48 dBで目安を十分満たす一方、PM=31.02°で目安(45°)を下回り、ステップ応答のオーバーシュートや振動が目立つ水準にあります。ゲイン余裕だけを見て「6 dB以上あるから大丈夫」と判断すると、位相余裕が示す過渡応答の質の悪化を見逃すことになります。
![L(s)=8/[s(s+1)(s+5)] のボード線図。ゲイン余裕11.48dBは目安を満たすが、位相余裕31.02°は目安(45°)を下回る](/posts/20260521_bode_plot/gain_phase_margin_bode.png)
実務ではGM・PMの両方を必ず併記し、どちらか一方だけで安定余裕を判断しないことが重要です。特に共振ピークを持つ系や高次系では、ゲイン交差点が複数存在する「条件付き安定」系も現れ、単純なGM/PMの定義自体が破綻する場合もあります(発展的な扱いは ナイキスト線図・根軌跡 を参照)。 PID制御 のゲインチューニングでも、このGM/PMの併記は欠かせません。
共振ピークと減衰係数 \(\zeta\) :\(M_p\) 公式の数値検証
式(7)の共振ピーク公式 \(M_p = 1/(2\zeta\sqrt{1-\zeta^2})\) (\(\zeta<1/\sqrt2\approx0.7071\) で有効)を、複数の \(\zeta\) について数値的に検証します。
import numpy as np
from scipy import signal
wn = 1.0
zetas = [0.05, 0.1, 0.2, 0.3, 0.5, 0.6, 1 / np.sqrt(2), 0.8, 1.0]
for zeta in zetas:
num = [wn ** 2]
den = [1, 2 * zeta * wn, wn ** 2]
sys = signal.TransferFunction(num, den)
w = np.logspace(-2, 1, 200000)
w, mag, phase = signal.bode(sys, w=w)
idx = np.argmax(mag)
if zeta < 1 / np.sqrt(2):
Mp_db_theory = 20 * np.log10(1 / (2 * zeta * np.sqrt(1 - zeta ** 2)))
else:
Mp_db_theory = 0.0 # zeta >= 1/sqrt(2) では共振ピークが消失
print(f"zeta={zeta:.4f}: Mp_theory={Mp_db_theory:8.4f} dB, Mp_numeric={mag[idx]:8.4f} dB")
実行結果:
zeta=0.0500: Mp_theory= 20.0109 dB, Mp_numeric= 20.0109 dB
zeta=0.1000: Mp_theory= 14.0230 dB, Mp_numeric= 14.0230 dB
zeta=0.2000: Mp_theory= 8.1361 dB, Mp_numeric= 8.1361 dB
zeta=0.3000: Mp_theory= 4.8466 dB, Mp_numeric= 4.8466 dB
zeta=0.5000: Mp_theory= 1.2494 dB, Mp_numeric= 1.2494 dB
zeta=0.6000: Mp_theory= 0.3546 dB, Mp_numeric= 0.3546 dB
zeta=0.7071: Mp_theory= 0.0000 dB, Mp_numeric= -0.0000 dB
zeta=0.8000: Mp_theory= 0.0000 dB, Mp_numeric= -0.0002 dB
zeta=1.0000: Mp_theory= 0.0000 dB, Mp_numeric= -0.0009 dB
\(\zeta<1/\sqrt2\) の範囲では理論値と数値計算値が小数点以下4桁まで完全に一致します。\(\zeta=1/\sqrt2\approx0.7071\) を境に共振ピークは消失し(\(M_p=1=0\) dB)、それ以上の \(\zeta\) では振幅は共振を持たず単調減少します(\(\zeta=0.8, 1.0\) の「数値値」がわずかに負なのは、ピークが存在しないため探索窓内の最大値が実質0dB地点に一致することによる数値誤差で、理論と矛盾しません)。\(\zeta\) が小さいほどピークは鋭く高くなり、\(\zeta=0.05\) では実に \(20\) dB(電圧比で10倍)ものピークが立つことが確認できます。

チェビシェフフィルタの通過域リップルや楕円フィルタの両域リップルは、この共振ピークが複数組み合わさった結果と見なせます。減衰係数 \(\zeta\) が小さい2次系のカスケードほど、フィルタの阻止域除去性能と引き換えに通過域近くでピーキング(リップル)が増大するのは、この \(M_p\) 公式が本質的な原因です。
最近の研究動向
ボード線図に基づく周波数領域解析は古典的な手法ですが、データ駆動制御やモデルフリー制御が主流になりつつある現在でも活発に研究が続いています。
- Isoshima, Tanemura, & Chida (2023) は、モデルを用いずに入出力データだけからゲイン余裕・位相余裕の下界を推定する手法を提案しました。古典的なボード線図解析がモデルベースであるのに対し、実測周波数応答データのみから安定余裕の保証を与える点が特徴です。
- Ren, Quan, Xu, Wang, & Cai (2024) は、多入力多出力(MIMO)系に対する安定余裕をデータ駆動で定義・推定する枠組みを提案しました。SISO系のGM/PMがそのままでは使えないMIMO系に対して、周波数領域データから実験的に安定余裕を求める方法を体系化しています。
- Broens, Butler, & Tóth (2024) は、周波数応答関数(FRF)の推定値に基づいて線形パラメータ変動(LPV)構造化制御器を自動チューニングする手法を提案し、精密位置決め系での高精度モーション制御に適用しています。ボード線図の自動読み取りとループ整形の自動化という方向性を示す研究です。
いずれも「ボード線図・周波数応答という古典的な表現形式」と「データ駆動・自動化という現代的な要求」を橋渡しする研究であり、本記事で導出した振幅・位相の基礎理論が、最新の制御研究でも変わらず土台になっていることがわかります。
まとめ
- ボード線図は周波数応答 \(H(j\omega)\) を 振幅[dB] + 位相[deg] の 2 段、対数周波数軸 で描く可視化
- デシベルにより直列接続が加算に、対数軸により広帯域挙動が等間隔で読める
- 1 次系のロールオフは \(\pm 20\) dB/dec、2 次系は \(\pm 40\) dB/dec、\(N\) 次は \(\pm 20N\) dB/dec で、これは \(\log_{10}\omega\) の一次関数への漸近と数値検証(式変形一致・カスケード加法性が小数点以下4桁一致)から厳密に導かれる
- バターワース・チェビシェフ・楕円・ハイパス・バンドパス・ノッチなど すべてのフィルタの個性がボード線図に現れる
- 非最小位相系(RHPゼロ)は最小位相系と振幅が完全に一致するが、位相は高域で \(180°\) 近くも余分に遅れる(数値検証: \(\omega=100\) で\(-86.56°\) vs \(-265.42°\) )
- ゲイン余裕と位相余裕は定義上異なる周波数で評価される独立な指標で、GM=11.48 dBと目安を満たしてもPM=31.02°と目安を下回る系が実在する(\(L(s)=8/[s(s+1)(s+5)]\) )
- 2次系の共振ピーク \(M_p=1/(2\zeta\sqrt{1-\zeta^2})\) (\(\zeta<1/\sqrt2\) )は理論値と数値計算値が小数点以下4桁一致し、\(\zeta\) が小さいほど鋭いピークになる
- Python では
scipy.signal.bode(連続時間)とscipy.signal.freqz(離散時間)が標準
各種フィルタの設計理論については以下の関連記事を参照してください。
関連記事
- ナイキスト線図・根軌跡・安定余裕の読み方とPython実装 - ボード線図で読み取る位相余裕・ゲイン余裕を、ナイキスト線図と根軌跡から閉ループ安定性として体系化した発展編です。制御系設計における本記事の直接の続きになります。
- PID制御の基礎理論と各要素の役割 - 本記事のゲイン余裕・位相余裕の議論は、PID制御のゲインチューニングにおける安定限界の評価に直接応用できます。
- PID制御のPython実装:シミュレーションとチューニング - PIDゲインを上げたときの過渡応答の悪化を、本記事のGM/PMの数値例と対応づけて理解できます。
- H∞制御とロバスト安定性の理論とPython実装 - ボード線図・GM/PMという古典的な周波数領域指標を、不確かさを明示的に扱うロバスト制御の枠組みへ一般化した発展編です。
- バターワースフィルタの設計とPython実装 - 最大平坦特性のロールオフをボード線図上で確認できる代表的 IIR フィルタです。
- チェビシェフフィルタの設計原理とPython実装 - 通過域リップルがボード線図でどう現れるかを学べます。
- ハイパスフィルタの設計とPython実装 - 低域ロールオフと \(+90^\circ\) への位相回転をボード線図で読み解けます。
- ローパスフィルタの設計と比較 - ボード線図比較の最も基本的な題材です。
- バンドパスフィルタの設計とPython実装 - 中心周波数まわりの帯域特性と位相の対称性が観察できます。
- ノッチフィルタの設計とPython実装 - \(f_0\) における深いディップと急峻な位相回転をボード線図で確認できます。
- FIRフィルタとIIRフィルタの比較 - 線形位相 FIR の位相プロットが直線になることをボード線図で理解できます。
- 適応フィルタ(LMS/RLS)の理論とPython実装 - 時不変のボード線図と対比される、時変フィルタの設計手法を解説しています。
- 高速フーリエ変換(FFT)の仕組みとPython実装 - ボード線図の実測に用いる FFT の理論を解説しています。
- 窓関数と PSD 推定の実践 - 実測ボード線図におけるリーケージ抑制と PSD 推定を解説しています。
- ベッセルフィルタの理論とPython実装 - 群遅延が最大限平坦なベッセルは、位相プロットがほぼ直線として可視化されます。バターワース・チェビシェフとの位相比較に最適です。
- 時間周波数解析の選び方ハブ - ボード線図で評価する時不変フィルタの周波数応答を出発点に、時変・非定常信号を扱う FFT・STFT・Wavelet・Hilbert変換の使い分けを整理した姉妹ハブです。
参考文献
- Bode, H. W. (1945). Network Analysis and Feedback Amplifier Design. Van Nostrand.
- Ogata, K. (2010). Modern Control Engineering (5th ed.). Prentice Hall.
- Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
- Isoshima, K., Tanemura, M., & Chida, Y. (2023). “Data-driven estimation of the lower bounds of gain and phase margins”. Automatica, 153, 111008.
- Ren, J., Quan, Q., Xu, B., Wang, S., & Cai, K.-Y. (2024). “Data-driven stability margin for linear multivariable systems”. International Journal of Robust and Nonlinear Control, 34(13), 8844-8862.
- Broens, Y., Butler, H., & Tóth, R. (2024). “Frequency Domain Auto-tuning of Structured LPV Controllers for High-Precision Motion Control”. arXiv:2403.05878.
- scipy.signal.bode — SciPy documentation
- scipy.signal.freqz — SciPy documentation
関連ツール
- DevToolBox - 開発者向け無料ツール集 - JSON整形、正規表現テスターなど85種類以上の開発者向けツール
- CalcBox - 暮らしの計算ツール - 統計計算、周波数変換など61種類以上の計算ツール