音楽プレイヤーのイコライザ、無線受信機のチャネル選択、心電図からの基線変動の除去 — これらはすべて「特定の周波数だけを通し、それ以外を遮断する」というフィルタの仕事です。ディジタル信号処理では、こうしたフィルタを差分方程式の係数 $\{b_m\}, \{a_k\}$ として実装します。では、その係数はどうやって決めるのでしょうか。「カットオフ 1 kHz の急峻なローパスが欲しい」という仕様から、いきなり係数が出てくるわけではありません。
実は、ディジタルIIRフィルタの設計には100年近い歴史を持つ強力な道筋があります。まずアナログフィルタの世界で、理論的に最適な「プロトタイプ（原型）」を設計し、それを双一次変換という数学的な写像でディジタルの世界へ移すのです。アナログ側で Butterworth・Chebyshev・楕円という3つの古典的プロトタイプが完成されており、それぞれ「最大平坦」「等リプルで急峻」「最も急峻」という明確な個性を持っています。この道筋を理解すると、SciPy の butter や ellip がブラックボックスではなくなり、自分の仕様に応じて方式と次数を根拠を持って選べるようになります。
この設計法は、レーダー受信機のアンチエイリアシングフィルタ、無線通信のチャネルフィルタ、生体信号計測のノイズ除去など、信号処理が登場するあらゆる現場で使われています。本記事では、3つのプロトタイプの振幅2乗特性 $|H(j\Omega)|^2$ を出発点に、極の配置まで数式で導出し、双一次変換と周波数プリワーピングを丁寧に導いたうえで、SciPy と自前実装で BPF/HPF を設計し、同一仕様で3方式を比較します。
本記事の内容
- アナログプロトタイプという設計思想（なぜアナログを経由するのか）
- Butterworth の振幅2乗特性と単位円上の極配置の導出
- Chebyshev の等リプルと Chebyshev 多項式 $T_n$ による特性の導出
- 楕円フィルタの等リプル特性と楕円関数の役割
- 双一次変換 $s = \frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}$ と周波数プリワーピングの導出
- ローパスプロトタイプから LPF/HPF/BPF への周波数変換
- SciPy と自前双一次変換による設計、3方式の比較プロット
設計は「仕様決定 → 正規化ローパスプロトタイプ → 周波数変換 → 双一次変換」という4段階のパイプラインで進みます。各ステップが独立して体系化されているため、アナログ理論の豊富な遺産を丸ごとディジタル設計に活かせます。プリワーピングはステップ③と④の境界で必要になる補正で、カットオフ周波数が変換後もぴったり一致するよう事前に$\Omega_c$を調整します。
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
特に、伝達関数 $H(z)$・極と零点・単位円上での周波数応答という概念は本記事の前提になります。アナログ側では複素変数 $s = \sigma + j\Omega$（ラプラス変数）を扱いますが、考え方は z 領域とよく似ています。
アナログプロトタイプとは — なぜわざわざアナログを経由するのか
ディジタルフィルタを設計したいのに、なぜいったんアナログフィルタを設計するのでしょうか。直感的なイメージとしては「すでに完成された名作レシピを翻訳して使う」ことに似ています。アナログフィルタの理論は1930年代から発展し、Butterworth・Chebyshev・楕円（Cauer）という、ある意味で「数学的に最適」なフィルタ特性が解析的に求められています。これらは連続周波数 $\Omega$（rad/s）の関数として、極と零点の位置が閉じた式で与えられます。ディジタル設計でこれを一から作るのは大変ですが、完成品を借りてきて写像するだけなら容易です。
もう一つの理由は、仕様の表現が自然であることです。フィルタの仕様は通常、「通過帯域のリプルは何 dB 以内」「阻止帯域では何 dB 以上減衰」「遷移帯域の幅」といった振幅特性で与えられます。アナログプロトタイプはこうした仕様を満たす最小次数を解析的に決められるため、設計が体系化されています。
ここでまず「プロトタイプ」という言葉を定義しておきます。正規化ローパスプロトタイプとは、カットオフ角周波数を $\Omega_c = 1$ rad/s に正規化したローパスフィルタのことです。実際に欲しいのが 1 kHz のハイパスであっても、まず $\Omega_c = 1$ のローパスを設計し、後述する周波数変換で目的のフィルタに変形します。設計の核心部分を一度だけ作り、あとは変換で使い回す、という発想です。
設計全体の流れは次のとおりです。
- 仕様（通過帯域端・阻止帯域端・リプル・減衰量）を決める
- 必要次数 $N$ を求め、正規化ローパスプロトタイプ $H_p(s)$ を得る（極・零点の式）
- 周波数変換でアナログの目的フィルタ $H_a(s)$（LPF/HPF/BPF/BSF）にする
- 双一次変換 $s \to z$ でディジタル $H(z)$ にする
- プリワーピングでカットオフ周波数のずれを補正する
まずはステップ2の主役である3つのプロトタイプを、振幅2乗特性から順に導いていきましょう。最初は最もシンプルな Butterworth です。
Butterworth プロトタイプの導出
振幅2乗特性と「最大平坦」の意味
この図から、Butterworthの振幅2乗特性の本質がはっきりわかります。次数$N$が大きいほど、通過帯域（$\Omega<1$）ではゲインが1.0に貼り付いたまま平坦で、$\Omega=1$のカットオフを超えると急激に落ちる「崖」の形が鮮明になります。$N=1$では緩やかな肩しかありませんが、$N=8$では理想矩形に近づき、「最大平坦（maximally flat）」という名称が視覚的に納得できます。
理想のローパスフィルタは、通過帯域でゲイン1、阻止帯域でゲイン0という「崖」のような特性を持ちますが、有限次数の有理関数ではこれを実現できません。そこで Butterworth は「通過帯域をできるだけ平らに保つ」という方針を取りました。具体的には、ゲインを $\Omega=0$ のまわりでテイラー展開したとき、できるだけ多くの微分係数がゼロになるように設計します。これが最大平坦（maximally flat）特性です。
Butterworth フィルタの振幅2乗特性は次式で定義されます。
$$ |H(j\Omega)|^2 = \frac{1}{1 + \left(\dfrac{\Omega}{\Omega_c}\right)^{2N}} $$
ここで $N$ は次数、$\Omega_c$ はカットオフ角周波数です。正規化プロトタイプでは $\Omega_c = 1$ として
$$ |H(j\Omega)|^2 = \frac{1}{1 + \Omega^{2N}} $$
とします。この式の振る舞いを見ると、$\Omega = 0$ で $|H|^2 = 1$、$\Omega = 1$ で $|H|^2 = 1/2$（すなわち $-3$ dB）、$\Omega \to \infty$ で $|H|^2 \to 0$ となり、ローパスらしい単調減少特性になっています。$\Omega^{2N}$ という高次の項のおかげで、$\Omega < 1$ ではゲインがほぼ1に貼りつき、$\Omega > 1$ で急速に落ちます。次数 $N$ を上げるほど崖が急になります。
極の導出 — 解析接続と $s$ 平面
振幅特性だけでは伝達関数 $H(s)$ は決まりません。$H(s)$ を求めるために、$|H(j\Omega)|^2$ を複素平面へ拡張します。実係数フィルタでは $|H(j\Omega)|^2 = H(j\Omega)H^*(j\Omega) = H(j\Omega)H(-j\Omega)$ が成り立ち、$s = j\Omega$ すなわち $\Omega = s/j = -js$ と置き換えると $\Omega^2 = -s^2$ なので、
$$ H(s)H(-s) = \frac{1}{1 + (-s^2)^N} = \frac{1}{1 + (-1)^N s^{2N}} $$
と書けます。この式の極、つまり分母が0になる $s$ を求めます。分母を0とおくと
$$ 1 + (-1)^N s^{2N} = 0 \quad\Longrightarrow\quad s^{2N} = -(-1)^{-N} = (-1)^{N+1} $$
となります。ここで $(-1)^{N+1} = e^{j\pi(N+1)}$ と指数表示し、さらに $2\pi$ の整数倍 $k$ を加えても同じ値であることを使うと、
$$ s^{2N} = e^{j\pi(N+1)}\, e^{j 2\pi k} = e^{j\pi(N + 1 + 2k)}, \quad k = 0, 1, \dots, 2N-1 $$
両辺の $2N$ 乗根を取ると、極は
$$ s_k = \exp\!\left(j\,\frac{\pi(N + 1 + 2k)}{2N}\right), \quad k = 0, 1, \dots, 2N-1 $$
と求まります。指数の中身を見ると、すべての極は絶対値が1、つまり $s$ 平面の単位円上に等間隔で並びます。これが Butterworth の美しい性質です。
安定な極の選択
$2N$ 個の極のうち、安定なフィルタ $H(s)$ には左半平面（$\mathrm{Re}(s) < 0$）にある $N$ 個だけを割り当てます。右半平面の極は $H(-s)$ の側に属します。連続系では左半平面の極が安定（インパルス応答が減衰）の条件だからです。安定極の偏角を整理すると
$$ \theta_k = \frac{\pi}{2} + \frac{(2k+1)\pi}{2N}, \quad k = 0, 1, \dots, N-1 $$
の形で、虚軸を中心に左半平面に対称配置されます。たとえば $N=2$ なら極は $s = e^{j 3\pi/4}, e^{j 5\pi/4}$、すなわち $-\frac{1}{\sqrt2}\pm j\frac{1}{\sqrt2}$ です。プロトタイプの伝達関数は零点を持たず（全極型）、
$$ H_p(s) = \frac{1}{\prod_{k=0}^{N-1}(s – s_k)} $$
となります。分母を展開すると $N$ 次の Butterworth 多項式が得られ、係数表としても知られています。
Butterworth は通過帯域が最も平坦で位相特性も比較的素直ですが、遷移が緩やかなので同じ急峻さを得るには次数が大きくなります。「通過帯域のなめらかさを少し犠牲にして、遷移をもっと急にできないか」と考えると、次の Chebyshev にたどり着きます。
Chebyshev 多項式 $T_N(x)$ の振る舞いを見ると、設計思想がそのまま図に現れています。左図の通過帯域（$|x|\leq 1$）では、$T_N$ が$-1$から$1$の間を$N+1$回振動し、等リプルを実現します。右図の阻止帯域（$x>1$）では $\cosh$ で急増し、次数が上がるほど急峻に増大します。これが「等リプルを許す代わりに遷移が急峻になる」Chebyshev の核心です。
Chebyshev プロトタイプの導出
等リプルという発想と Chebyshev 多項式
Butterworth は通過帯域を平らにすることに全力を注いだ結果、遷移が緩慢になりました。発想を変えて、「通過帯域に一定の小さな波（リプル）を許す代わりに、その分のエネルギーを遷移の急峻さに回す」のが Chebyshev フィルタです。許したリプルを通過帯域全体に均等にばらまく（等リプルにする）ことで、与えられた次数とリプル量に対して最も急峻な遷移が得られます。これはミニマックス的な最適性として知られています。
この等リプルを実現する道具が Chebyshev 多項式 $T_N(x)$ です。定義は
$$ T_N(x) = \begin{cases} \cos(N \arccos x), & |x| \le 1 \\ \cosh(N \,\mathrm{arccosh}\, x), & |x| > 1 \end{cases} $$
です。$|x|\le 1$ の領域で $T_N(x)$ は $\cos$ なので $[-1, 1]$ の間を $N+1$ 回振動し、$|x|>1$ では $\cosh$ なので急激に増大します。この「内側で振動し外側で爆発する」性質が、通過帯域の等リプルと阻止帯域の急峻な減衰の両方を生みます。漸化式 $T_{N+1}(x) = 2x\,T_N(x) – T_{N-1}(x)$、$T_0=1$、$T_1=x$ で計算できます。
振幅2乗特性
Chebyshev I 型（通過帯域等リプル型）の振幅2乗特性は
$$ |H(j\Omega)|^2 = \frac{1}{1 + \varepsilon^2\, T_N^2\!\left(\dfrac{\Omega}{\Omega_c}\right)} $$
で定義されます。ここで $\varepsilon$ はリプルの大きさを決めるパラメータです。正規化プロトタイプ（$\Omega_c=1$）では引数が $\Omega$ になります。通過帯域 $|\Omega|\le 1$ では $T_N(\Omega)\in[-1,1]$ なので $T_N^2\in[0,1]$ をとり、ゲインは
$$ \frac{1}{1+\varepsilon^2} \le |H|^2 \le 1 $$
の間で $N$ 回振動します。これが等リプルです。リプルの深さ（dB）は
$$ R_p = 10\log_{10}(1+\varepsilon^2) \quad\Longrightarrow\quad \varepsilon = \sqrt{10^{R_p/10} – 1} $$
で $\varepsilon$ と結びつきます。たとえば $R_p = 1$ dB なら $\varepsilon \approx 0.5088$ です。一方 $\Omega > 1$ では $T_N$ が $\cosh$ で急増するので、ゲインは Butterworth より速く落ちます。
極の導出 — 楕円上の配置
極を求めるため、Butterworth と同様に $s = j\Omega$ で解析接続し、分母 $1 + \varepsilon^2 T_N^2(\Omega) = 0$ を解きます。$\Omega = s/j$ を入れ、$T_N(x) = \cos(N\arccos x)$ を用います。$w = \arccos(s/j)$ と置くと条件は
$$ \cos(N w) = \pm \frac{j}{\varepsilon} $$
となります。$w = u + jv$ と実部虚部に分け、$\cos(N w) = \cos(Nu)\cosh(Nv) – j\sin(Nu)\sinh(Nv)$ を代入すると、右辺が純虚数であることから実部は0、虚部は $\mp 1/\varepsilon$ でなければなりません。実部が0となる条件
$$ \cos(Nu)\cosh(Nv) = 0 \quad\Longrightarrow\quad \cos(Nu) = 0 \quad\Longrightarrow\quad u_k = \frac{(2k+1)\pi}{2N} $$
が得られます（$\cosh \ge 1 > 0$ なので $\cos(Nu)=0$ が必要）。虚部の条件 $\sin(Nu_k)\sinh(Nv) = 1/\varepsilon$ から、$\sin(Nu_k)=\pm1$ なので
$$ \sinh(Nv) = \frac{1}{\varepsilon} \quad\Longrightarrow\quad v = \frac{1}{N}\,\mathrm{arcsinh}\frac{1}{\varepsilon} \equiv v_0 $$
と一意に定まります。最後に $s = j\cos w = j\cos(u + jv)$ を展開すると、安定（左半平面）な極は
$$ s_k = -\sinh(v_0)\sin\!\left(\frac{(2k+1)\pi}{2N}\right) + j\,\cosh(v_0)\cos\!\left(\frac{(2k+1)\pi}{2N}\right), \quad k=0,\dots,N-1 $$
と求まります。この式をよく見ると、実部は $\sinh v_0$ で、虚部は $\cosh v_0$ でスケールされており、極は $s$ 平面上で横軸半径 $\sinh v_0$・縦軸半径 $\cosh v_0$ の楕円の上に並びます。Butterworth の円が楕円に押しつぶされた形です。$\varepsilon$ を小さくすると（リプルを浅くすると）$v_0$ が大きくなって楕円が円に近づき、Butterworth に漸近します。
Chebyshev は通過帯域にリプルを許す代わりに遷移を急峻にしました。Chebyshev II 型（逆 Chebyshev）は逆に通過帯域を平坦に保ち阻止帯域を等リプルにします。では「通過帯域も阻止帯域も両方リプルを許せば、もっと急峻にできるのでは」という究極の発想が、次の楕円フィルタです。
楕円フィルタプロトタイプの導出
両側等リプルと最小次数
楕円フィルタ（Cauer フィルタ）は、通過帯域も阻止帯域も両方とも等リプルにしたフィルタです。許容できるリプルと減衰の予算を、通過帯域と阻止帯域の両方に均等にばらまくことで、所定の仕様を満たす最小次数を達成します。同じ遷移帯域幅・同じリプル・同じ減衰の仕様なら、楕円フィルタが最も低い次数で実現できます。代償は位相特性が大きく歪むこと（群遅延のばらつきが大きい）です。
振幅2乗特性は、Chebyshev の $T_N$ をChebyshev 有理関数 $R_N$ に置き換えた形
$$ |H(j\Omega)|^2 = \frac{1}{1 + \varepsilon^2\, R_N^2\!\left(\Omega,\,\xi\right)} $$
で与えられます。ここで $R_N$ は楕円有理関数（Chebyshev rational function）で、選択係数 $\xi$ をパラメータに持ちます。$R_N$ は通過帯域 $|\Omega|\le 1$ で $[-1,1]$ を等リプルに振動し、阻止帯域 $|\Omega|\ge \Omega_s$ では $[1/L, \infty)$ の領域で振動するように作られた有理関数です。$T_N$ が多項式（零点だけ）だったのに対し、$R_N$ は分子と分母の両方に零点を持つ点が決定的に異なります。
楕円関数による零点・極の配置
$R_N$ の零点 $z_i$ と極 $p_i$ は、ヤコビの楕円関数 $\mathrm{sn}, \mathrm{cn}, \mathrm{dn}$ を使って与えられます。直感的には、Chebyshev のときに使った $\cos$（円関数、一重周期）を、二重周期を持つ楕円関数に置き換えることで、通過帯域と阻止帯域の二つの領域で同時に等リプルを実現しているのです。具体的には、楕円有理関数の零点は
$$ \Omega_{z,i} = \mathrm{cd}\!\left(\frac{(2i-1)}{N}K(k),\; k\right), \quad i=1,\dots, \lfloor N/2\rfloor $$
の形で与えられ（$\mathrm{cd} = \mathrm{cn}/\mathrm{dn}$、$K(k)$ は第一種完全楕円積分、$k$ は選択係数 $\xi$ と関連するモジュラス）、阻止帯域の零点は $\Omega_s/\Omega_{z,i}$ に対応して配置されます。これにより、阻止帯域の特定周波数で振幅が完全にゼロになる伝達零点（虚軸上の零点）が生じます。
楕円フィルタの伝達関数 $H_p(s)$ は、左半平面の極（楕円関数で与えられる）と虚軸上の零点を持ち、
$$ H_p(s) = H_0 \frac{\prod_i (s^2 + \Omega_{0,i}^2)}{\prod_k (s – p_k)} $$
の形になります。分子の $s^2 + \Omega_{0,i}^2$ が阻止帯域の伝達零点で、これが「阻止帯域に深いノッチを刻む」働きをします。導出の全体は楕円積分の逆関数や次数方程式（degree equation）$\frac{N K'(k)}{K(k)} = \frac{K'(k_1)}{K(k_1)}$ を要し相当に込み入るため、本記事では特性と配置の意味を押さえることに留め、数値計算は SciPy の ellipap に委ねます。重要なのは「両側等リプル＝最小次数」という設計上の価値と、「阻止帯域に虚軸上の零点（ノッチ）を持つ」という構造的特徴です。
ここまでで、3つのアナログプロトタイプ $H_p(s)$ が（全極型あるいは零点付きで）得られました。しかしこれらはまだ連続時間のフィルタです。これをディジタルの $H(z)$ にするための橋渡しが、いよいよ双一次変換です。
3方式の $s$ 平面での極零配置を見ると、理論が視覚的に確認できます。Butterworth の極は完全な単位円上に等間隔で並び（零点なし）、Chebyshev I の極は虚軸方向に伸びた楕円上に配置されています。楕円フィルタだけが虚軸上に零点（図中の○印）を持ち、これが阻止帯域に深いノッチを刻む伝達零点です。すべての極が左半平面（Re(s) < 0）に収まっており、安定なプロトタイプであることも確認できます。
同一次数$N=5$の3プロトタイプを正規化周波数で比較すると、遷移の急峻さの違いが一目瞭然です。楕円フィルタが最も急峻に落ち、Chebyshev がそれに続き、Butterworth が最もなだらかです。ただし楕円は通過帯域にリプルが見え、Butterworth は完全に平坦です。「急峻さと通過帯域の平坦さはトレードオフ」という設計原則を、このグラフが定量的に裏づけています。
双一次変換の導出
なぜ単純な置き換えではダメなのか
アナログ伝達関数 $H_a(s)$ をディジタル $H(z)$ にするには、$s$ と $z$ を結ぶ写像が必要です。理論的には $z = e^{sT}$（$T$ はサンプリング周期）が厳密な関係ですが、これは超越関数で有理関数になりません。そこで近似が要ります。素朴に微分を差分で置き換える（$s \approx (1-z^{-1})/T$、後退差分）と、$s$ 平面の左半平面が $z$ 平面の単位円内にきれいに収まらず、安定性が保証されません。双一次変換は、左半平面を単位円内に過不足なく写すように設計された写像で、安定なアナログフィルタが必ず安定なディジタルフィルタになることを保証します。
台形則からの導出
双一次変換は、積分を台形則で離散化することから導けます。1次の系 $\dot{y}(t) = x(t)$（積分器、$Y(s) = X(s)/s$）を考え、$t = nT$ と $t=(n-1)T$ の間で積分を台形近似します。
$$ y(nT) = y((n-1)T) + \int_{(n-1)T}^{nT} x(\tau)\,d\tau \approx y((n-1)T) + \frac{T}{2}\bigl(x(nT) + x((n-1)T)\bigr) $$
両辺を z 変換します。$y[n] – y[n-1] = \frac{T}{2}(x[n] + x[n-1])$ なので、
$$ Y(z)(1 – z^{-1}) = \frac{T}{2}X(z)(1 + z^{-1}) $$
積分器の伝達関数 $Y/X$ は $z$ 領域で
$$ \frac{Y(z)}{X(z)} = \frac{T}{2}\,\frac{1 + z^{-1}}{1 – z^{-1}} $$
となります。一方アナログの積分器は $Y/X = 1/s$ ですから、両者を等しいとみなして $1/s$ をこの式に対応させると、$s$ は
$$ \boxed{\; s = \frac{2}{T}\,\frac{1 – z^{-1}}{1 + z^{-1}} \;} $$
と置き換えればよいことがわかります。これが双一次変換です。$H(z) = H_a(s)\big|_{s = \frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}}$ として代入すれば、有理関数のまま $z$ の有理関数が得られます。
安定性の保証 — 左半平面が単位円内へ
この図は、双一次変換が安定性を保存する仕組みを直感的に示しています。左図（$s$ 平面）の青色の左半平面にある極は、右図（$z$ 平面）の単位円の内側に過不足なく収まります。単位円の内側が安定領域ですから、安定なアナログフィルタから設計した場合、必ず安定なディジタルフィルタが得られます。右半平面への飛び出しが起こらないことも確認できます。
この写像が安定性を保つことを確認します。$z = e^{j\omega}$（単位円上）を双一次変換に入れると、
$$ s = \frac{2}{T}\,\frac{1 – e^{-j\omega}}{1 + e^{-j\omega}} $$
分子分母に $e^{j\omega/2}$ を掛けてオイラーの公式で整理すると、
$$ s = \frac{2}{T}\,\frac{e^{j\omega/2} – e^{-j\omega/2}}{e^{j\omega/2} + e^{-j\omega/2}} = \frac{2}{T}\,\frac{2j\sin(\omega/2)}{2\cos(\omega/2)} = j\,\frac{2}{T}\tan\frac{\omega}{2} $$
となり、$s$ は純虚数になります。つまり単位円（$|z|=1$）はちょうど虚軸（$s=j\Omega$）に写るのです。さらに $z$ 平面の単位円内（$|z|<1$）は $s$ 平面の左半平面に、単位円外は右半平面に写ることも示せます。したがってアナログの安定極（左半平面）は必ずディジタルの安定極（単位円内）に写り、安定性が保たれます。
周波数プリワーピング
上の式 $\Omega = \frac{2}{T}\tan\frac{\omega}{2}$ は、アナログ角周波数 $\Omega$ とディジタル角周波数 $\omega$（$= 2\pi f/f_s$）の対応関係そのものです。ここに重大な注意点があります。この対応は線形ではなく $\tan$ で歪んでいるのです。低い周波数では $\tan(\omega/2)\approx \omega/2$ なのでほぼ線形ですが、$\omega$ がナイキスト（$\pi$）に近づくと $\Omega$ が急激に伸び、無限大の $\Omega$ が有限の $\omega=\pi$ に押し込まれます。これが周波数ワーピング（warping）です。
何も対策しないと、設計したカットオフ周波数が双一次変換後にずれてしまいます。そこで、アナログプロトタイプを設計する段階で、あらかじめ目的のディジタルカットオフ $\omega_c$ に対応するアナログ周波数
$$ \Omega_c = \frac{2}{T}\tan\frac{\omega_c}{2} $$
を使ってプロトタイプを設計します。これをプリワーピング（prewarping）と呼びます。あらかじめ $\tan$ の歪みを織り込んでおくことで、双一次変換後にカットオフがぴったり $\omega_c$ に来ます。具体的には、設計仕様の各エッジ周波数を上式で $\Omega$ に変換してからプロトタイプ設計・周波数変換を行い、最後に双一次変換します。
左図のワーピング曲線を見ると、低周波（50 Hz 以下）では線形（破線）とほぼ重なりますが、ナイキスト（500 Hz）に近づくにつれてアナログ周波数が急激に伸びているのがわかります。右図はプリワーピングの必要性を示し、補正なしでは 100 Hz に設定したはずのカットオフがデジタル変換後にずれてしまいます。プリワーピング済みの $\Omega_c$ を使うと、変換後にカットオフがぴったり 100 Hz に来ることが確認できます。
双一次変換で $s\to z$ の橋ができました。あとは、正規化ローパスプロトタイプ $H_p(s)$ を目的の LPF/HPF/BPF に変える「アナログ周波数変換」を組み合わせれば設計が完結します。
ローパスから所望特性への周波数変換
正規化ローパスプロトタイプ（カットオフ $\Omega=1$）の変数 $s$ を別の式で置き換えることで、目的のフィルタに変形できます。代表的な変換は次のとおりです。$\Omega_c$ は目的のカットオフ（プリワーピング済み）です。
- ローパス→ローパス（カットオフ $\Omega_c$ へ）: $s \to \dfrac{s}{\Omega_c}$
- ローパス→ハイパス（カットオフ $\Omega_c$ へ）: $s \to \dfrac{\Omega_c}{s}$
- ローパス→バンドパス（中心 $\Omega_0=\sqrt{\Omega_1\Omega_2}$、帯域幅 $B=\Omega_2-\Omega_1$）: $s \to \dfrac{s^2 + \Omega_0^2}{B\,s}$
- ローパス→バンドストップ: $s \to \dfrac{B\,s}{s^2 + \Omega_0^2}$
ハイパス変換 $s\to \Omega_c/s$ の意味を確認します。プロトタイプで $\Omega=1$ がカットオフでしたが、$s\to\Omega_c/s$ では $\Omega$ が大きいほど変換後の引数が小さくなり、低周波と高周波の役割が入れ替わります。結果として通過帯域が高域側に移り、ハイパスになります。バンドパス変換は $s$ の次数を2倍にするため、$N$ 次プロトタイプから $2N$ 次のバンドパスが生まれます。中心周波数 $\Omega_0$ で $s^2+\Omega_0^2=0$ となり、プロトタイプの $s=0$（最大ゲイン点）に対応するため、中心で通過します。
これで、(1) 仕様をプリワーピングして $\Omega$ に直す、(2) プロトタイプ $H_p(s)$ を得る、(3) 周波数変換で $H_a(s)$ にする、(4) 双一次変換で $H(z)$ にする、という設計パイプラインが完成しました。次はこれを Python で実装し、3方式を比較します。
Python での実装
プロトタイプの極配置を確認する
まず Butterworth と Chebyshev I 型の極が、それぞれ単位円・楕円の上に並ぶことを自前計算と SciPy で確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# ===== 自前計算: Butterworth と Chebyshev I の極（正規化プロトタイプ） =====
N = 6 # 次数
# --- Butterworth: 左半平面の極（単位円上） ---
k = np.arange(N)
theta = np.pi/2 + (2*k + 1)*np.pi/(2*N)
poles_butter = np.exp(1j*theta) # 単位円上
# --- Chebyshev I: 楕円上の極 ---
Rp = 1.0 # 通過帯域リプル [dB]
eps = np.sqrt(10**(Rp/10) - 1.0)
v0 = np.arcsinh(1.0/eps)/N
ang = (2*k + 1)*np.pi/(2*N)
poles_cheby = -np.sinh(v0)*np.sin(ang) + 1j*np.cosh(v0)*np.cos(ang)
# SciPy のプロトタイプ取得関数と照合
zb, pb, kb = signal.buttap(N)
zc, pc, kc = signal.cheb1ap(N, Rp)
print("Butterworth |pole|（自前）:", np.round(np.abs(poles_butter), 4))
print("Butterworth |pole|（SciPy）:", np.round(np.sort(np.abs(pb)), 4))
print("Chebyshev 実部範囲（自前）:", np.round(poles_cheby.real.min(), 4),
"〜", np.round(poles_cheby.real.max(), 4))
print("Chebyshev |pole|（SciPy）:", np.round(np.sort(np.abs(pc)), 4))
出力を見ると、Butterworth の極は自前計算・SciPy ともに絶対値がすべて 1.0 で、確かに単位円上に乗っています。Chebyshev の極は絶対値が 1 より小さくばらつき、楕円上に押しつぶされていることが数値で確認できます。自前の式と SciPy の buttap/cheb1ap が一致することから、本文の導出が正しいこともわかります。
極配置を可視化する
数値だけでなく、$s$ 平面での極配置を図にして「円か楕円か」を目で確かめます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
N = 6
Rp, Rs = 1.0, 40.0
# 各プロトタイプの極・零点
zb, pb, kb = signal.buttap(N)
zc, pc, kc = signal.cheb1ap(N, Rp)
ze, pe, ke = signal.ellipap(N, Rp, Rs)
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
theta = np.linspace(0, 2*np.pi, 300)
for ax, (name, z, p) in zip(
axes,
[('Butterworth', zb, pb), ('Chebyshev I', zc, pc), ('Elliptic', ze, pe)]):
ax.plot(np.cos(theta), np.sin(theta), 'k--', lw=0.6) # 単位円（参照）
ax.scatter(p.real, p.imag, marker='x', s=80, color='red', label='Poles')
if len(z) > 0:
ax.scatter(z.real, z.imag, marker='o', s=80,
facecolors='none', edgecolors='blue', label='Zeros')
ax.axvline(0, color='gray', lw=0.5)
ax.set_title(f'{name} prototype (N={N})')
ax.set_xlabel('Re(s)'); ax.set_ylabel('Im(s)')
ax.set_aspect('equal'); ax.grid(alpha=0.3); ax.legend()
plt.tight_layout()
plt.savefig('prototype_poles.png', dpi=150, bbox_inches='tight')
plt.show()
この図から3方式の構造的な違いが一目でわかります。Butterworth の極は破線の単位円にぴったり乗り、零点はありません（全極型）。Chebyshev I の極は虚軸方向に伸びた楕円上に並び、やはり零点を持ちません。楕円フィルタだけは虚軸上に零点を持っており、これが阻止帯域のノッチ（深い減衰点）に対応します。すべての極が左半平面（Re(s) < 0）にあり、安定なアナログプロトタイプであることも確認できます。
自前の双一次変換でディジタル化する
SciPy の bilinear に頼らず、プリワーピングを含む双一次変換を自前で実装し、ローパスをディジタル化します。アナログプロトタイプを lp2lp でカットオフ変換し、bilinear 相当の写像を多項式の係数レベルで行います。
import numpy as np
from scipy import signal
def my_bilinear(b_s, a_s, fs):
"""アナログ伝達関数 (b_s, a_s) を双一次変換でディジタル化する自前実装。
s = 2*fs*(1 - z^-1)/(1 + z^-1) を代入し、z^-1 の多項式に展開する。"""
K = 2.0 * fs
# s -> K*(z-1)/(z+1) と見て、(z+1)^M, (z-1)^k の積で展開する
M = max(len(b_s), len(a_s)) - 1
# 係数を次数 M にそろえる（高次=左）
b_s = np.r_[np.zeros(M + 1 - len(b_s)), b_s]
a_s = np.r_[np.zeros(M + 1 - len(a_s)), a_s]
def s_power_to_z(coeffs):
# coeffs: s^M, s^(M-1), ..., s^0 の係数
out = np.zeros(M + 1)
for i, c in enumerate(coeffs):
p = M - i # s の指数
# s^p = (K)^p * (z-1)^p * (z+1)^(M-p) を z^-? 多項式へ
term = (K**p) * np.polymul(np.poly1d([1, -1])**p if p > 0 else np.poly1d([1]),
np.poly1d([1, 1])**(M - p) if (M - p) > 0 else np.poly1d([1]))
out += c * np.r_[np.zeros(M + 1 - len(term.coeffs)), term.coeffs]
return out
b_z = s_power_to_z(b_s)
a_z = s_power_to_z(a_s)
# a_z[0] で正規化
b_z = b_z / a_z[0]
a_z = a_z / a_z[0]
return b_z, a_z
fs = 1000.0
fc = 100.0
# プリワーピング: アナログ角周波数 Omega_c
wc = 2*fs*np.tan(np.pi*fc/fs) # = 2*fs*tan(omega_c/2), omega_c = 2*pi*fc/fs
# 4次 Butterworth プロトタイプ -> カットオフ wc のアナログLPF
z, p, k = signal.buttap(4)
b_s, a_s = signal.zpk2tf(z, p, k)
b_s, a_s = signal.lp2lp(b_s, a_s, wo=wc)
# 自前双一次 と SciPy を比較
b_my, a_my = my_bilinear(b_s, a_s, fs)
b_sp, a_sp = signal.bilinear(b_s, a_s, fs)
print("自前 b:", np.round(b_my, 6))
print("SciPy b:", np.round(b_sp, 6))
print("最大係数差:", np.max(np.abs(np.r_[b_my - b_sp, a_my - a_sp])))
出力では自前実装と SciPy の bilinear の係数がほぼ完全に一致し（差は浮動小数点誤差レベル）、本文で導いた $s = 2f_s(1-z^{-1})/(1+z^{-1})$ の写像が正しく実装できていることが確認できます。プリワーピング $\Omega_c = 2f_s\tan(\pi f_c/f_s)$ を入れた点が重要で、これにより双一次変換後のカットオフがちょうど 100 Hz に来ます。
3方式の振幅特性を同一仕様で比較する
同じ仕様（通過帯域端・阻止帯域端・リプル・減衰）で3方式を設計し、必要次数の違いを見ます。SciPy の buttord/cheb1ord/ellipord で最小次数を自動算出させます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
fs = 1000.0
fp = 100.0 # 通過帯域端 [Hz]
fst = 150.0 # 阻止帯域端 [Hz]
Rp = 1.0 # 通過帯域リプル [dB]
Rs = 50.0 # 阻止帯域減衰 [dB]
wp, ws = fp/(fs/2), fst/(fs/2) # 正規化（ナイキスト=1）
# 各方式の最小次数を算出
Nb, wnb = signal.buttord(wp, ws, Rp, Rs)
Nc, wnc = signal.cheb1ord(wp, ws, Rp, Rs)
Ne, wne = signal.ellipord(wp, ws, Rp, Rs)
print(f"必要次数 Butterworth={Nb}, Chebyshev I={Nc}, Elliptic={Ne}")
# 設計（双一次変換は SciPy 内部で実行される）
bb, ab = signal.butter(Nb, wnb)
bc, ac = signal.cheby1(Nc, Rp, wnc)
be, ae = signal.ellip(Ne, Rp, Rs, wne)
plt.figure(figsize=(11, 6))
for (b, a), name, N in [((bb, ab), 'Butterworth', Nb),
((bc, ac), 'Chebyshev I', Nc),
((be, ae), 'Elliptic', Ne)]:
w, H = signal.freqz(b, a, worN=4096, fs=fs)
plt.plot(w, 20*np.log10(np.abs(H) + 1e-12), lw=2, label=f'{name} (N={N})')
plt.axvline(fp, color='gray', ls='--', alpha=0.6)
plt.axvline(fst, color='gray', ls=':', alpha=0.6)
plt.axhline(-Rs, color='red', ls=':', alpha=0.4)
plt.xlim(0, fs/2); plt.ylim(-90, 5)
plt.xlabel('Frequency [Hz]'); plt.ylabel('Magnitude [dB]')
plt.title('Same spec, 3 prototypes: magnitude response')
plt.legend(); plt.grid(alpha=0.3)
plt.savefig('three_proto_mag.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから設計上のトレードオフが鮮明に読み取れます。同じ仕様（通過 100 Hz・阻止 150 Hz・リプル 1 dB・減衰 50 dB）を満たすのに、Butterworth は最も高い次数（$N=15$）を要し、Chebyshev I はそれより少なく（$N=7$）、楕円フィルタが最小次数（$N=5$）で済みます。楕円フィルタの阻止帯域には深いノッチ（伝達零点）が並び、$-50$ dB ラインを最小限の次数で割り込んでいます。Butterworth は通過帯域が完全に平坦、Chebyshev は通過帯域にリプルが見え、楕円は通過・阻止の両方にリプルが見える点も、本文の振幅2乗特性の式どおりです。
群遅延と遷移帯域のトレードオフを定量化する
急峻さの代償として現れる群遅延（位相歪み）を比較します。急峻なフィルタほど通過帯域端で群遅延が暴れる傾向を確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
fs = 1000.0
fp, fst, Rp, Rs = 100.0, 150.0, 1.0, 50.0
wp, ws = fp/(fs/2), fst/(fs/2)
designs = {}
Nb, wnb = signal.buttord(wp, ws, Rp, Rs); designs['Butterworth'] = signal.butter(Nb, wnb)
Nc, wnc = signal.cheb1ord(wp, ws, Rp, Rs); designs['Chebyshev I'] = signal.cheby1(Nc, Rp, wnc)
Ne, wne = signal.ellipord(wp, ws, Rp, Rs); designs['Elliptic'] = signal.ellip(Ne, Rp, Rs, wne)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
for name, (b, a) in designs.items():
# 群遅延
w, gd = signal.group_delay((b, a), w=4096, fs=fs)
ax1.plot(w, gd, lw=2, label=name)
# 位相
w2, H = signal.freqz(b, a, worN=4096, fs=fs)
ax2.plot(w2, np.unwrap(np.angle(H))*180/np.pi, lw=2, label=name)
for ax in (ax1, ax2):
ax.axvline(fp, color='gray', ls='--', alpha=0.6)
ax.set_xlim(0, fs/2); ax.set_xlabel('Frequency [Hz]')
ax.grid(alpha=0.3); ax.legend()
ax1.set_ylabel('Group delay [samples]'); ax1.set_title('Group delay')
ax2.set_ylabel('Phase [deg]'); ax2.set_title('Phase response')
plt.tight_layout()
plt.savefig('three_proto_delay.png', dpi=150, bbox_inches='tight')
plt.show()
この図から、急峻さと位相のトレードオフが定量的に見えます。Butterworth の群遅延は通過帯域内で最もなだらかで、位相も比較的素直です。楕円フィルタは通過帯域端（100 Hz 付近）で群遅延が鋭いピークを示し、波形を最も歪ませます。Chebyshev はその中間です。つまり「次数を下げ遷移を急峻にした分のツケが位相歪みとして現れる」という関係が読み取れ、波形保存が重要な用途（オーディオ・パルス伝送）では Butterworth、選択度最優先の用途（無線チャネル分離）では楕円、という選択指針が裏づけられます。
HPF と BPF を周波数変換で設計する
最後に、同じプロトタイプから周波数変換でハイパスとバンドパスを設計し、ローパス以外への展開を確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
fs = 1000.0
# HPF: カットオフ 150 Hz、楕円 4次
b_hp, a_hp = signal.ellip(4, 1, 50, 150, btype='highpass', fs=fs)
# BPF: 通過 100-200 Hz、Chebyshev I 4次（=8次相当の帯域通過）
b_bp, a_bp = signal.cheby1(4, 1, [100, 200], btype='bandpass', fs=fs)
plt.figure(figsize=(11, 6))
for (b, a), name in [((b_hp, a_hp), 'Elliptic HPF (fc=150Hz)'),
((b_bp, a_bp), 'Chebyshev BPF (100-200Hz)')]:
w, H = signal.freqz(b, a, worN=4096, fs=fs)
plt.plot(w, 20*np.log10(np.abs(H) + 1e-12), lw=2, label=name)
plt.xlim(0, fs/2); plt.ylim(-90, 5)
plt.xlabel('Frequency [Hz]'); plt.ylabel('Magnitude [dB]')
plt.title('HPF / BPF via frequency transform + bilinear')
plt.legend(); plt.grid(alpha=0.3)
plt.savefig('hpf_bpf.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフから、ローパスプロトタイプが周波数変換と双一次変換を経て、確かにハイパス・バンドパス・バンドストップとして機能していることが確認できます。楕円 HPF は 150 Hz より下を急峻に減衰させ、阻止帯域（低域側）にノッチを持ちます。Chebyshev BPF は 100〜200 Hz だけを通し、通過帯域内に等リプルが見えます。バンドストップは BPF の補として 80〜120 Hz を深く阻止しています。バンドパスは内部で次数が2倍になるため、低域側・高域側の両方の遷移帯域が形成されている点も、$s\to(s^2+\Omega_0^2)/(Bs)$ という変換が次数を倍にする性質と整合します。
まとめ
本記事では、アナログプロトタイプから双一次変換を経てディジタルIIRを設計する道筋を解説しました。
- アナログプロトタイプ経由: 正規化ローパス $\Omega_c=1$ のプロトタイプを設計し、周波数変換と双一次変換でディジタル目的フィルタに写すのが古典的IIR設計の王道
- Butterworth: 振幅2乗 $|H|^2 = 1/(1+\Omega^{2N})$。極は $s$ 平面の単位円上に等間隔配置。最大平坦で位相が素直、ただし遷移が緩やか
- Chebyshev I: $|H|^2 = 1/(1+\varepsilon^2 T_N^2(\Omega))$。通過帯域等リプル。極は楕円上に配置され、Butterworth より低次で急峻
- 楕円: 楕円有理関数 $R_N$ による両側等リプル。阻止帯域に虚軸上の零点（ノッチ）を持ち、同一仕様で最小次数を達成。代償は位相歪み
- 双一次変換: $s = \frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}$。台形則から導かれ、左半平面を単位円内に写すので安定性を保存。$\Omega = \frac{2}{T}\tan\frac{\omega}{2}$ の歪みをプリワーピングで補正
- トレードオフ: 急峻さ（楕円）↔ 位相の素直さ・低群遅延（Butterworth）。用途で選ぶ
3方式の比較から、波形保存重視なら Butterworth、選択度最優先なら楕円、という設計指針が定量的に裏づけられました。次のステップとして、以下の記事も参考にしてください。