双一次変換（バイリニア変換）— アナログフィルタをデジタル化する設計法

「アナログフィルタの設計技術をそのままデジタルで使えたら、どれだけ楽だろう」—— デジタルフィルタ設計を学ぶとき、多くの人がそう考えます。実は、その夢を叶える変換が存在します。双一次変換（バイリニア変換）です。

Butterworth フィルタや Chebyshev フィルタは、アナログ回路の世界で100年以上かけて磨き上げられた設計法です。これらのアナログプロトタイプ $H(s)$ をそのまま流用し、1つの代数的な置き換えでデジタルフィルタ $H(z)$ を得る——それが双一次変換の核心です。この変換は安定性を厳密に保存し、折り返し歪み（エイリアシング）を一切生じません。

双一次変換が活躍する応用場面は非常に広い範囲にわたります。音声処理 では、オーディオイコライザや音声強調フィルタの設計に使われ、アナログ回路のロールオフ特性をデジタル DSP 上で忠実に再現します。生体信号計測 では、心電図（ECG）の基線変動除去や脳波（EEG）の特定周波数帯抽出に使われ、医療機器の FIR では表現しにくい急峻な IIR 特性が威力を発揮します。scipy の butter・cheby1・ellip 関数の内部でも、まさにこの双一次変換が動いています。

本記事の内容

• 双一次変換の変換式 $s = \frac{2}{T}\cdot\frac{z-1}{z+1}$ の由来——台形積分近似から自然に導出する
• $s$ 平面から $z$ 平面への写像と安定性保存の仕組み
• 周波数ワーピング $\Omega = \frac{2}{T}\tan\!\left(\frac{\omega_d}{2}\right)$ の導出と意味
• プリワーピング——カットオフ周波数を事前に補正する技術
• RC ローパスフィルタを例にした手計算と scipy による確認
• インパルス不変法との比較

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

双一次変換とは何か

「アナログフィルタの伝達関数 $H(s)$ に、$s = f(z)$ という式を代入して $H(z)$ を得る」——それが双一次変換の全体像です。ただし、どんな代入式でもよいわけではありません。良い変換は次の3つの条件を満たす必要があります。

1. 安定性保存: アナログフィルタが安定（左半面に極）なら、デジタルフィルタも安定（単位円内に極）であること
2. エイリアシング不発生: デジタル変換の折り返し歪みがないこと
3. 周波数特性の対応: アナログの通過域と遮断域の特性がデジタル側でも保たれること

双一次変換はこの3条件をすべて満たす唯一の有理関数変換です。その式は

$$ \boxed{s = \frac{2}{T}\cdot\frac{z – 1}{z + 1}} $$

と表されます。ここで $T$ はサンプリング周期です。この式を逆に解いて $z$ を $s$ の関数として表すと

$$ z = \frac{1 + \frac{T}{2}s}{1 – \frac{T}{2}s} $$

となります。「$s = 0$ のとき $z = 1$」「$s \to +\infty$ のとき $z \to -1$」——この対応関係が周波数ワーピングの原因でもあり、安定性保存の源でもあります。

なぜこの式が上記3条件を満たすのかは、導出と写像の解析を通じて明らかになります。まずはこの式がどこから来るのかを見ていきましょう。

変換式の導出——台形積分近似から

積分の離散化という視点

双一次変換の式 $s = \frac{2}{T}\cdot\frac{z-1}{z+1}$ は、天下り的に与えられたように見えますが、実は台形積分近似という自然なアイデアから導出できます。

連続時間システムにおいて、積分器の伝達関数は $H(s) = \frac{1}{s}$ です。「アナログの積分器をデジタルに変換する」という問いを解くと、双一次変換の式が自然に現れます。

連続時間の積分

$$ y(t) = \int_{0}^{t} x(\tau)\, d\tau $$

を離散化するとき、最も簡単な方法は「前の時刻での長方形近似（後退差分）」です。しかし精度が低いため、代わりに台形近似を使います。

上の図を見ると、後退差分は時刻 $nT$ での値だけで長方形を作るのに対し、台形近似は両端の平均を使ってより正確な面積を計算しています。具体的には、時刻 $t = nT$ での積分値の近似式は

$$ y[n] = y[n-1] + \frac{T}{2}(x[n] + x[n-1]) $$

となります。右辺の第2項が「現時刻と1サンプル前の値の平均 × 時間幅 $T$」という台形の面積です。

この差分方程式に $z$ 変換を適用します。$z$ 変換の遅延則 $y[n-1] \leftrightarrow z^{-1}Y(z)$ を使うと、

$$ Y(z) = z^{-1}Y(z) + \frac{T}{2}(1 + z^{-1})X(z) $$

左辺に $Y(z)$ をまとめると、

$$ Y(z)(1 – z^{-1}) = \frac{T}{2}(1 + z^{-1})X(z) $$

デジタル積分器の伝達関数は

$$ H(z) = \frac{Y(z)}{X(z)} = \frac{T}{2}\cdot\frac{1 + z^{-1}}{1 – z^{-1}} = \frac{T}{2}\cdot\frac{z + 1}{z – 1} $$

となります（分子分母に $z$ を掛けて整理しました）。

変換式の導出

アナログ積分器の伝達関数は $H(s) = \frac{1}{s}$ でした。「アナログ積分器がデジタル積分器に対応する」という要請から、次の等式を立てます。

$$ \frac{1}{s} \longleftrightarrow \frac{T}{2}\cdot\frac{z + 1}{z – 1} $$

この対応関係から、$\frac{1}{s}$ を $\frac{T}{2}\cdot\frac{z+1}{z-1}$ に置き換えれば、任意のアナログフィルタをデジタル化できます。これは逆に言うと、$s$ を

$$ \boxed{s = \frac{2}{T}\cdot\frac{z – 1}{z + 1}} $$

に置き換えることと等価です。これが双一次変換の変換式の正体です。台形積分という「精度の高い離散化の選択」から自然に導出されました。

ここで重要なのは、台形近似の代数的な構造です。後退差分なら $s \approx \frac{1-z^{-1}}{T}$ となり（インパルス不変法の基礎）、前進差分なら $s \approx \frac{z-1}{T}$ となります。双一次変換はこれら2つの「平均」とも解釈でき、安定性に関して後者が不安定を引き起こすのに対し、台形近似は安定性を保存するという優れた性質を持っています。

s 平面から z 平面への写像

安定性がなぜ保存されるのか

双一次変換の最大の強みは「安定性の保存」です。アナログフィルタが安定（全ての極が $s$ 平面の左半面にある）なら、変換後のデジタルフィルタも必ず安定（全ての極が $z$ 平面の単位円内にある）——これを数式で確認しましょう。

$s = \sigma + j\Omega$（$\sigma$ は実部、$\Omega$ は虚数角周波数）に対し、$z$ の絶対値を計算します。

変換式 $z = \frac{1 + \frac{T}{2}s}{1 – \frac{T}{2}s}$ に $s = \sigma + j\Omega$ を代入すると、

$$ z = \frac{1 + \frac{T}{2}(\sigma + j\Omega)}{1 – \frac{T}{2}(\sigma + j\Omega)} $$

分子と分母の絶対値を計算します。分子 $|1 + \frac{T}{2}(\sigma + j\Omega)| = \sqrt{\left(1 + \frac{T\sigma}{2}\right)^2 + \left(\frac{T\Omega}{2}\right)^2}$、分母 $|1 – \frac{T}{2}(\sigma + j\Omega)| = \sqrt{\left(1 – \frac{T\sigma}{2}\right)^2 + \left(\frac{T\Omega}{2}\right)^2}$ となるので、

$$ |z|^2 = \frac{\left(1 + \frac{T\sigma}{2}\right)^2 + \left(\frac{T\Omega}{2}\right)^2}{\left(1 – \frac{T\sigma}{2}\right)^2 + \left(\frac{T\Omega}{2}\right)^2} $$

分子から分母を引くと、

$$ |z|^2 – 1 = \frac{\left(1 + \frac{T\sigma}{2}\right)^2 – \left(1 – \frac{T\sigma}{2}\right)^2}{\left(1 – \frac{T\sigma}{2}\right)^2 + \left(\frac{T\Omega}{2}\right)^2} $$

分子を展開すると $\left(1 + \frac{T\sigma}{2}\right)^2 – \left(1 – \frac{T\sigma}{2}\right)^2 = 2 \cdot 2 \cdot \frac{T\sigma}{2} = 2T\sigma$ となるので、

$$ |z|^2 – 1 = \frac{2T\sigma}{\left(1 – \frac{T\sigma}{2}\right)^2 + \left(\frac{T\Omega}{2}\right)^2} $$

分母は常に正なので、$|z|^2 – 1$ の符号は $\sigma$ の符号に一致します。つまり、

• $\sigma < 0$（左半面）$\Rightarrow$ $|z| < 1$（単位円内）: 安定
• $\sigma = 0$（虚軸上）$\Rightarrow$ $|z| = 1$（単位円上）: 境界
• $\sigma > 0$（右半面）$\Rightarrow$ $|z| > 1$（単位円外）: 不安定

図に示すように、$s$ 平面の虚軸（$\sigma = 0$）上の点は $z$ 平面の単位円上に写像され、$s$ 平面の左半面全体が $z$ 平面の単位円内に写像されます。アナログの安定条件（左半面）とデジタルの安定条件（単位円内）がこの変換によって厳密に対応しているのです。

これはインパルス不変法と大きく異なる点です。インパルス不変法では、折り返し誤差によって高域に歪みが生じます。一方、双一次変換では $s$ 平面のすべての周波数（$\Omega: -\infty \to +\infty$）が $z$ 平面の単位円（$\omega_d: -\pi \to +\pi$）に一対一に対応します。この全単射性が、エイリアシングをゼロにする理由です。

この写像の代償として、周波数の対応関係が非線形（ワーピング）になります。次のセクションでその詳細を見ていきましょう。

周波数ワーピングの導出

デジタル周波数とアナログ周波数の非線形な対応

$s$ 平面の虚軸 $s = j\Omega$ を双一次変換で $z$ 平面に写像すると、$z = e^{j\omega_d}$ という単位円上の点に対応します（$\omega_d$ はデジタル角周波数）。この対応関係を数式で求めましょう。

まず $z = e^{j\omega_d}$ を変換式 $s = \frac{2}{T}\cdot\frac{z-1}{z+1}$ に代入します。

$$ s = \frac{2}{T}\cdot\frac{e^{j\omega_d} – 1}{e^{j\omega_d} + 1} $$

指数を半角に分離するため、分子分母に $e^{-j\omega_d/2}$ を掛けます。

分子: $e^{j\omega_d/2} – e^{-j\omega_d/2} = 2j\sin(\omega_d/2)$

分母: $e^{j\omega_d/2} + e^{-j\omega_d/2} = 2\cos(\omega_d/2)$

これを代入すると、

$$ s = \frac{2}{T}\cdot\frac{2j\sin(\omega_d/2)}{2\cos(\omega_d/2)} = j\cdot\frac{2}{T}\tan\!\left(\frac{\omega_d}{2}\right) $$

$s = j\Omega$ と比較することで、アナログ角周波数 $\Omega$ とデジタル角周波数 $\omega_d$ の対応関係が得られます。

$$ \boxed{\Omega = \frac{2}{T}\tan\!\left(\frac{\omega_d}{2}\right)} $$

これが周波数ワーピングの式です。逆を求めると

$$ \omega_d = 2\arctan\!\left(\frac{\Omega T}{2}\right) $$

となります。

ワーピングの意味と影響

上の図が周波数ワーピング曲線です。グレーの破線が「理想的な線形対応」、青の実線が実際のワーピング曲線を示しています。

2つの重要な特徴が読み取れます。まず、低域では線形に近い——$\omega_d$ が小さい範囲では $\tan(\omega_d/2) \approx \omega_d/2$ なので $\Omega \approx \omega_d/T$ となり、ほぼ比例します。次に、高域（ナイキスト周波数 $\omega_d \to \pi$）に近づくにつれて急激に圧縮される——$\tan(\pi/2) \to \infty$ なので、アナログの無限大の周波数がデジタルのナイキスト周波数 $\pi$ にまとめて押し込まれます。

この非線形な対応が意味することは、「アナログプロトタイプの周波数スケールとデジタルフィルタの周波数スケールは非線形に対応する」ということです。低域フィルタの通過域・遮断域の境界（カットオフ周波数）の位置が、単純にスケールを変えても一致しません。ここで「プリワーピング」が必要になります。

プリワーピング——設計周波数の事前補正

なぜプリワーピングが必要か

双一次変換を使って「カットオフ 1 kHz のデジタルローパスフィルタを設計したい」とします。直感的には「アナログで 1 kHz のフィルタを設計し、双一次変換でデジタル化する」と考えるでしょう。

しかしここに落とし穴があります。ワーピングが非線形なので、アナログでカットオフ $\Omega_c$ に設計しても、双一次変換後のデジタルフィルタのカットオフは $\omega_d = 2\arctan(\Omega_c T/2)$ となり、目標の値からずれてしまいます。

図から明確に読み取れるように、プリワーピングなし（赤線）では、アナログで設計した周波数が双一次変換後にデジタル域でずれてしまいます。一方、プリワーピングあり（緑線）では、設計後のデジタルフィルタが正確に目標のカットオフ周波数に対応しています。

プリワーピングの計算

解決策は逆方向から考えることです。

目標: デジタルフィルタのカットオフが $\omega_c$（rad/sample）になってほしい。

手順: ワーピング式を使って、このデジタル周波数に対応するアナログ周波数を先に計算し、そのアナログ周波数でプロトタイプを設計する。

$$ \Omega_c = \frac{2}{T}\tan\!\left(\frac{\omega_c}{2}\right) $$

この $\Omega_c$ を使ってアナログプロトタイプを設計し、双一次変換すると、変換後のデジタルフィルタのカットオフがちょうど $\omega_c$ になります。

具体的な数値例を考えましょう。サンプリング周波数 $f_s = 8000$ Hz、目標カットオフ $f_c = 1000$ Hz の場合、

デジタル角周波数: $\omega_c = 2\pi \cdot \frac{f_c}{f_s} = 2\pi \cdot \frac{1000}{8000} = \frac{\pi}{4}$ rad/sample

サンプリング周期: $T = 1/f_s = 1.25 \times 10^{-4}$ s

プリワーピング後のアナログ角周波数:

$$ \Omega_c = \frac{2}{T}\tan\!\left(\frac{\pi/4}{2}\right) = \frac{2}{T}\tan\!\left(\frac{\pi}{8}\right) \approx \frac{2}{T} \times 0.4142 = 6628\, \text{rad/s} $$

これに対し、素直にスケール変換した場合は $\Omega_{\text{naive}} = \omega_c / T = \pi/4 \times 8000 = 6283$ rad/s となり、345 rad/s（実周波数で 55 Hz 相当）のずれが生じます。フィルタの次数が高くなるほど、このわずかなずれが特性に響くため、プリワーピングは必須の補正です。

プリワーピングの効果は低域では小さいですが、カットオフ周波数がナイキスト周波数 $f_s/2$ に近づくほど補正量が大きくなります。

RC ローパスフィルタを例にした設計手順

アナログ伝達関数の確認

最もシンプルなアナログフィルタとして、RC ローパスフィルタ（1次バターワース）を使って双一次変換の全手順を追いましょう。

RC ローパスフィルタのアナログ伝達関数は

$$ H(s) = \frac{\Omega_c}{s + \Omega_c} $$

です。ここで $\Omega_c$ はカットオフ角周波数です。$\omega = \Omega_c$ のとき $|H(j\Omega_c)| = 1/\sqrt{2} \approx 0.707$（$-3$ dB）になることは、分子に $j\Omega_c$ を代入して確認できます。

双一次変換の手順（手計算）

仕様: カットオフ $f_c = 1000$ Hz、サンプリング周波数 $f_s = 8000$ Hz。

ステップ 1: デジタル角周波数を求める

$$ \omega_c = 2\pi \cdot \frac{f_c}{f_s} = 2\pi \cdot \frac{1000}{8000} = \frac{\pi}{4} \approx 0.7854 \text{ rad/sample} $$

ステップ 2: プリワーピングでアナログ設計周波数を求める（$T = 1/8000$ s）

$$ \Omega_c = \frac{2}{T}\tan\!\left(\frac{\omega_c}{2}\right) = \frac{2}{1/8000}\tan\!\left(\frac{\pi}{8}\right) = 16000 \times 0.41421 \approx 6628 \text{ rad/s} $$

ステップ 3: アナログプロトタイプを設計する

$$ H(s) = \frac{\Omega_c}{s + \Omega_c} = \frac{6628}{s + 6628} $$

ステップ 4: 双一次変換を適用する。$s = \frac{2}{T}\cdot\frac{z-1}{z+1}$ を代入すると

$$ H(z) = \frac{6628}{\dfrac{2}{T}\cdot\dfrac{z-1}{z+1} + 6628} $$

分母を通分するため、分子分母に $\frac{z+1}{2/T}$ を掛けます。

$$ H(z) = \frac{6628 \cdot \frac{T}{2}(z+1)}{(z-1) + 6628 \cdot \frac{T}{2}(z+1)} $$

$T = 1/8000$ s、$\frac{T}{2} = 1/16000$ s を代入し、$6628 \cdot \frac{T}{2} = 6628/16000 \approx 0.4143$ と定義します（$k = \Omega_c \cdot T/2$）。

$$ H(z) = \frac{k(z+1)}{(z-1) + k(z+1)} = \frac{k(z+1)}{(1+k)z + (k-1)} $$

分子分母を $(1+k)$ で割り、$z^{-1}$ 形式に変換します。

$$ H(z) = \frac{\dfrac{k}{1+k}(1 + z^{-1})}{1 + \dfrac{k-1}{1+k}z^{-1}} $$

$k = 0.4143$ のとき $b_0 = b_1 = \frac{k}{1+k} = \frac{0.4143}{1.4143} \approx 0.2929$、$a_1 = \frac{k-1}{1+k} = \frac{-0.5857}{1.4143} \approx -0.4142$ となります。

これが双一次変換で得られたデジタルフィルタの係数です。差分方程式で書くと

$$ y[n] = 0.2929\,x[n] + 0.2929\,x[n-1] + 0.4142\,y[n-1] $$

となります。

scipy による確認

手計算の係数が正しいか、Python で検証します。

import numpy as np
from scipy import signal

# 設計仕様
fs = 8000    # サンプリング周波数 [Hz]
fc = 1000    # カットオフ周波数 [Hz]
T = 1.0 / fs

# プリワーピング（scipy が自動でやってくれるが、原理を確認）
omega_d = 2 * np.pi * fc / fs           # デジタル角周波数
Omega_c = (2 / T) * np.tan(omega_d / 2) # プリワーピング後のアナログ周波数
print(f"プリワーピング後の Ω_c = {Omega_c:.2f} rad/s")
print(f"素直なスケール変換: {omega_d / T:.2f} rad/s")

# アナログプロトタイプ係数 (scipy)
b_a = np.array([Omega_c])       # 分子: Ω_c
a_a = np.array([1.0, Omega_c])  # 分母: s + Ω_c

# 双一次変換
b_d, a_d = signal.bilinear(b_a, a_a, fs=fs)
print(f"\n双一次変換後の係数:")
print(f"  b = {b_d}")
print(f"  a = {a_d}")

このコードを実行すると次のような出力が得られます。

プリワーピング後の Ω_c = 6627.42 rad/s
素直なスケール変換: 6283.19 rad/s

双一次変換後の係数:
  b = [0.29289322 0.29289322]
  a = [ 1.        -0.41421356]

手計算で求めた $b_0 = b_1 \approx 0.2929$、$a_1 \approx -0.4142$ とぴったり一致しています。scipy の bilinear 関数が内部でプリワーピングを含む双一次変換を正確に実行していることが確認できます。

手計算の値と scipy の値が数値的に一致したことで、変換手順が正しいことが裏付けられました。次は、この係数で得られるデジタルフィルタの周波数応答を可視化してみましょう。

振幅応答と極零配置の確認

アナログとデジタルの振幅応答比較

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

for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

fs = 8000
fc = 1000
T = 1.0 / fs
omega_d = 2 * np.pi * fc / fs
Omega_c = (2 / T) * np.tan(omega_d / 2)

# アナログフィルタの周波数応答
f = np.logspace(1, np.log10(fs / 2), 1000)
omega = 2 * np.pi * f
H_analog = Omega_c / (1j * omega + Omega_c)
amp_analog = 20 * np.log10(np.abs(H_analog))

# デジタルフィルタの周波数応答
b, a = signal.bilinear([Omega_c], [1, Omega_c], fs=fs)
f_dig, H_dig = signal.freqz(b, a, worN=1000, fs=fs)
amp_dig = 20 * np.log10(np.abs(H_dig) + 1e-12)

fig, ax = plt.subplots(figsize=(9, 5))
ax.semilogx(f, amp_analog, lw=2.5, color="#1d4ed8", label="アナログ $H(s)$")
ax.semilogx(f_dig, amp_dig, "--", lw=2.5, color="#dc2626", label="デジタル $H(z)$（双一次変換）")
ax.axvline(fc, color="#374151", lw=1.5, ls=":")
ax.axhline(-3, color="#d97706", lw=1.5, ls=":")
ax.plot(fc, -3, "s", color="#d97706", markersize=10, label=f"カットオフ {fc} Hz")

ax.set_xlabel("周波数 [Hz]", fontsize=11)
ax.set_ylabel("振幅特性 [dB]", fontsize=11)
ax.set_title("RC ローパスフィルタ: アナログ vs デジタル振幅応答", fontsize=12, fontweight="bold")
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, which="both")
ax.set_xlim(10, fs / 2)
ax.set_ylim(-60, 5)
plt.tight_layout()
plt.show()

# カットオフ確認
idx = np.argmin(np.abs(f_dig - fc))
print(f"デジタルフィルタのカットオフ付近の振幅: {amp_dig[idx]:.3f} dB（目標 -3 dB）")

デジタルフィルタのカットオフ付近の振幅: -3.011 dB（目標 -3 dB）

グラフから2つのことが読み取れます。まず、1000 Hz のカットオフ点で両フィルタがほぼ完全に一致しており、プリワーピングの効果が確認できます（差は 0.011 dB にすぎません）。次に、高域（数 kHz 以上）でデジタルフィルタの遮断がアナログより若干急峻になっていますが、これはワーピングによる高域圧縮の効果です。主要な通過域では両者の応答はほぼ同一であることが分かります。

極零写像の確認

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

for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

fs = 8000
fc = 1000
T = 1.0 / fs
omega_d = 2 * np.pi * fc / fs
Omega_c = (2 / T) * np.tan(omega_d / 2)

# デジタルフィルタ係数
b, a = signal.bilinear([Omega_c], [1, Omega_c], fs=fs)

# 極と零点
zeros, poles, gain = signal.tf2zpk(b, a)
print(f"零点: {zeros}")
print(f"極:   {poles}")
print(f"ゲイン: {gain:.4f}")

零点: [-1.+0.j]
極:   [0.41421356+0.j]
ゲイン: 0.29289322

極零マップから、双一次変換の特徴的な構造が読み取れます。アナログフィルタの極 $s = -\Omega_c \approx -6628$ rad/s（左半面）は、デジタルフィルタの極 $z \approx 0.4142$（単位円内の正実軸上）に写像されています。そして、アナログにはなかった零点 $z = -1$（ナイキスト周波数） が双一次変換によって自動的に生成されます。

この $z = -1$ の零点は双一次変換の普遍的な性質です。変換式 $H(z) = \frac{k(z+1)}{(1+k)z + (k-1)}$ の分子に $(z+1)$ があることから、ナイキスト周波数（$z = e^{j\pi} = -1$）での利得は常に 0（$-\infty$ dB）になります。ローパスフィルタにとってこれは都合が良い特性ですが、ハイパスやバンドパスを設計する際は注意が必要です。

ワーピングによる高域圧縮

双一次変換の周波数ワーピングは、高域において「無限大の周波数を有限の $\pi$ に押し込む」という効果を持ちます。これはどういう意味でしょうか。

アナログフィルタの遮断域（高域）がデジタルフィルタではナイキスト周波数 $\omega_d = \pi$ 付近に圧縮されます。その結果、アナログよりもデジタルの方が高域の遮断が急峻に見えます。

左のグラフは低域〜カットオフ付近の応答で、プリワーピングによってカットオフ周波数が正確に一致しています。右のグラフは高域を拡大したものです。青の実線（双一次変換）がグレーの破線（ワーピングなしの理想的な線形対応）と比べて、ナイキスト周波数付近でより急激に落ちていることが分かります。

この高域圧縮効果はデジタルフィルタの設計においては概ね有益です。遮断域での減衰がより急峻になるため、同じ次数でより鋭いロールオフが得られます。一方で、高域の位相応答はアナログプロトタイプとは異なる挙動を示すため、位相特性が重要な用途（音声コーデックの位相一致など）では注意が必要です。

インパルス不変法との比較

双一次変換の理解を深めるため、もう一つの代表的な変換手法であるインパルス不変法と比較します。

インパルス不変法はアナログフィルタのインパルス応答 $h(t)$ を $T$ ごとにサンプリングして $h[n] = h(nT)$ とすることでデジタルフィルタを作る方法です。時間領域でのサンプリングなので、周波数域ではエイリアシングが生じます。

たとえば 1次ローパス $H(s) = \Omega_c / (s + \Omega_c)$ のインパルス応答は $h(t) = \Omega_c e^{-\Omega_c t} u(t)$ です。これをサンプリングすると

$$ h[n] = \Omega_c e^{-\Omega_c nT} u[n] $$

$z$ 変換すると

$$ H_{\text{II}}(z) = \frac{\Omega_c}{1 – e^{-\Omega_c T}z^{-1}} $$

となります。この式にはインパルス不変法の特徴が表れています。第一に、零点がありません（$z = -1$ の零点が存在しない）。第二に、高い周波数（$\omega_d$ が大きい領域）でアナログとデジタルの特性がずれます。このずれがエイリアシングの現れです。

実際に数値で比較すると（fc = 1 kHz、fs = 8 kHz）：

import numpy as np
from scipy import signal

fs = 8000
fc = 1000
T = 1.0 / fs
omega_d = 2 * np.pi * fc / fs
Omega_c = (2 / T) * np.tan(omega_d / 2)

# 双一次変換（プリワーピングあり）
b_blt, a_blt = signal.bilinear([Omega_c], [1, Omega_c], fs=fs)

# インパルス不変法: h[n] = h(nT) = Ω_c * exp(-Ω_c * n * T) をサンプリング
# DC ゲインが 1 になるよう正規化: b = (1 - exp(-Ω_c * T))
Omega_c_ii = 2 * np.pi * fc        # プリワーピングなし
alpha = np.exp(-Omega_c_ii * T)    # 極: z = e^{-Ω_c T}
b_ii = np.array([1 - alpha])       # DC 正規化
a_ii = np.array([1.0, -alpha])

# 周波数応答を計算
f_blt, H_blt = signal.freqz(b_blt, a_blt, worN=2000, fs=fs)
f_ii, H_ii = signal.freqz(b_ii, a_ii, worN=2000, fs=fs)

idx_blt = np.argmin(np.abs(f_blt - fc))
idx_ii = np.argmin(np.abs(f_ii - fc))

amp_blt = 20 * np.log10(np.abs(H_blt[idx_blt]))
amp_ii = 20 * np.log10(np.abs(H_ii[idx_ii]))
amp_ii_nyq = 20 * np.log10(np.abs(H_ii[-1]))

print(f"双一次変換: カットオフ付近の振幅 = {amp_blt:.3f} dB")
print(f"インパルス不変法: カットオフ付近の振幅 = {amp_ii:.3f} dB")
print(f"インパルス不変法: ナイキスト {fs//2} Hz での振幅 = {amp_ii_nyq:.3f} dB")

双一次変換: カットオフ付近の振幅 = -3.010 dB
インパルス不変法: カットオフ付近の振幅 = -2.793 dB
インパルス不変法: ナイキスト 4000 Hz での振幅 = -8.550 dB

結果から2点が読み取れます。カットオフ付近（1 kHz）ではインパルス不変法も約 -2.8 dB と比較的近い値ですが、これはエイリアシングによる折り返しが通過域の利得を持ち上げているためで、正しく -3 dB を実現しているわけではありません。より顕著な違いはナイキスト周波数（4 kHz）付近で、双一次変換は -70 dB 近い急峻な遮断を示すのに対し、インパルス不変法ではエイリアシングにより高域が -8.6 dB にとどまります（十分に遮断できていない）。高域をしっかり抑圧したい用途では、双一次変換が明確に優位です。カットオフ周波数の正確な制御と高域遮断の確実性の両面から、双一次変換が実用の IIR フィルタ設計で広く選ばれる理由はここにあります。

完全な設計フロー

ここまでの内容を整理して、双一次変換による IIR フィルタ設計の完全なフローを示します。

設計フローは6ステップで構成されます。このフローにより、アナログ設計のノウハウがデジタル実装に直接活かせます。

scipy での実装（Butterworth ローパスの例）

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

for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False

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

# scipy の butter は内部でプリワーピング込みの双一次変換を実行
b, a = signal.butter(order, fc, btype="low", fs=fs)

# 周波数応答
f, H = signal.freqz(b, a, worN=2000, fs=fs)
amp = 20 * np.log10(np.abs(H) + 1e-12)
phase = np.angle(H) * 180 / np.pi

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 7), sharex=True)

ax1.semilogx(f, amp, color="#1d4ed8", lw=2.5)
ax1.axvline(fc, color="#dc2626", lw=1.5, ls=":")
ax1.axhline(-3, color="#d97706", lw=1.5, ls=":")
ax1.set_ylabel("振幅 [dB]", fontsize=11)
ax1.set_title(f"4次 Butterworth LPF（双一次変換、カットオフ {fc} Hz）", fontsize=12, fontweight="bold")
ax1.set_ylim(-80, 5)
ax1.grid(True, alpha=0.3, which="both")
ax1.legend([f"$H(z)$（{order}次）", f"カットオフ {fc} Hz", "-3 dB"], fontsize=10)

ax2.semilogx(f, phase, color="#16a34a", lw=2.5)
ax2.axvline(fc, color="#dc2626", lw=1.5, ls=":")
ax2.set_xlabel("周波数 [Hz]", fontsize=11)
ax2.set_ylabel("位相 [deg]", fontsize=11)
ax2.set_ylim(-200, 10)
ax2.grid(True, alpha=0.3, which="both")

plt.tight_layout()
plt.show()

# カットオフ確認
idx = np.argmin(np.abs(f - fc))
print(f"カットオフ {fc} Hz での振幅: {amp[idx]:.3f} dB")
print(f"3 kHz での減衰: {amp[np.argmin(np.abs(f - 3000))]:.1f} dB")
print(f"係数 b: {b}")
print(f"係数 a: {a}")

カットオフ 1500 Hz での振幅: -3.010 dB
3 kHz での減衰: -44.6 dB
係数 b: [0.03797304 0.15189217 0.22783825 0.15189217 0.03797304]
係数 a: [ 1.         -0.97836878  0.79008574 -0.24188218  0.03773388]

4次 Butterworth フィルタでカットオフ 1500 Hz が正確に -3 dB に実現されています。3 kHz（カットオフの 2 倍）では -44.6 dB の減衰が得られており、4次フィルタとして十分に急峻な特性が確認できます（ワーピングによる高域圧縮で、線形スケールでの2倍周波数よりも急峻になっています）。b 係数が対称構造（$b_0 = b_4$、$b_1 = b_3$）を持つのは 4次バターワースの特徴で、この係数をそのまま差分方程式に代入して実装できます。

設計上の注意点

次数の選択

双一次変換は n 次アナログフィルタを n 次デジタルフィルタに変換します。次数が上がると急峻なロールオフが得られますが、係数の量子化誤差の影響も大きくなります（特に固定小数点実装）。

実用的な目安として、4〜8次が多くの用途でのバランスが良い選択です。高精度な設計では、2次セクション（バイカッド）の直列接続として実装することで、数値安定性が大幅に向上します。

通過域リプルと遷移帯域幅の調整

フィルタの種類によって特性が大きく変わります。

• Butterworth: 通過域が最大平坦（リプルなし）。遷移が緩やか
• Chebyshev I型: 通過域に等リプル。同次数でより急峻
• Chebyshev II型: 遮断域に等リプル。通過域は平坦
• 楕円（Cauer）: 通過域・遮断域両方に等リプル。最も急峻

双一次変換はどのアナログプロトタイプにも同様に適用できます。scipy の signal.cheby1・signal.ellip も内部で双一次変換を使用しています。

ハイパス・バンドパスへの拡張

ローパスプロトタイプから他のフィルタ形状（ハイパス・バンドパス・バンドストップ）への変換も周波数変換と双一次変換を組み合わせることで実現できます。

たとえばハイパスフィルタは $s \to \Omega_c / s$ という置き換えをアナログ段階で行ってから双一次変換します。または、scipy の btype="high" 引数で直接指定することもできます。

from scipy import signal

fs = 8000
fc = 1000
b_hp, a_hp = signal.butter(4, fc, btype="high", fs=fs)
b_bp, a_bp = signal.butter(4, [800, 2000], btype="band", fs=fs)

バンドパス（band）では 2倍の次数（4次指定で8次のデジタルフィルタ）になる点に注意が必要です。

まとめ

本記事では、双一次変換（バイリニア変換）の原理から実践的な設計手順まで解説しました。

• 変換の由来: $s = \frac{2}{T}\cdot\frac{z-1}{z+1}$ は台形積分近似から自然に導出される
• 安定性保存: $s$ 平面の左半面が $z$ 平面の単位円内に厳密に写像される（$|z|^2 – 1$ の符号は $\sigma$ の符号に等しい）
• エイリアシングなし: $s$ 平面の全域が単位円上に一対一対応するため、折り返し誤差は生じない
• 周波数ワーピング: $\Omega = \frac{2}{T}\tan(\omega_d/2)$ という非線形な対応が生じる
• プリワーピング: カットオフ周波数を $\Omega_c = \frac{2}{T}\tan(\omega_c/2)$ で事前補正することで、目標のデジタルカットオフを正確に実現できる
• RC 例での確認: $b_0 = b_1 = 0.2929$、$a_1 = -0.4142$ という係数が scipy と完全一致

双一次変換を習得すると、Butterworth・Chebyshev・楕円フィルタといったアナログの古典設計資産をデジタル DSP で完全に活用できるようになります。次のステップとして、IIR フィルタの安定実装（バイカッド構造）や、より高次な設計仕様（等リプル帯域、群遅延補正）への応用に進んでみてください。