群遅延と線形位相 — フィルタが波形を歪ませない条件

オーディオのイコライザでノイズを取り除いたら、なぜか打楽器のアタックが鈍くなった——そんな経験はないでしょうか。あるいは、心電図にローパスフィルタをかけたら、鋭いはずのR波のピーク位置がずれてしまい、心拍間隔の計測が狂ってしまった。これらの不可解な「波形の崩れ」は、フィルタの振幅特性ではなく位相特性が引き起こしています。フィルタは「どの周波数をどれだけ通すか(振幅)」だけでなく、「どの周波数をどれだけ遅らせるか(位相)」も決めているのです。

ここで決定的に重要なのが、各周波数成分の遅延がそろっているかどうかです。すべての周波数が同じ時間だけ遅れるなら、波形は形を保ったまま平行移動するだけで済みます。ところが周波数ごとに遅延がバラバラだと、成分どうしの足し合わせのタイミングがずれ、鋭いエッジは鈍り、波形は別物に変わってしまいます。この「遅延のそろい具合」を定量化するのが群遅延(group delay)であり、群遅延を一定に保つフィルタが線形位相(linear phase)フィルタです。

群遅延と線形位相を理解すると、次のような実務上の判断ができるようになります。たとえば、レーダーやソナーのパルス圧縮では、受信パルスの形状そのものが情報なので、波形を歪ませない線形位相フィルタが必須です。また、複数のセンサ信号を時間軸で比較するセンサフュージョンでは、各チャンネルのフィルタ遅延がそろっていないと相関や時刻差の推定が破綻します。逆に、リアルタイム制御のように遅延そのものを嫌う場面では、線形位相FIRの大きな遅延を避けて低次のIIR(位相は非線形になる)を選ぶ、といった設計判断も出てきます。本記事では、位相遅延と群遅延の定義からはじめて、線形位相がなぜ波形を保存するのか、FIRがどうして対称係数で厳密な線形位相になるのか、そしてIIRが線形位相を持てない事情と全域通過フィルタによる群遅延等化までを、数式の導出とPython実装の両面から丁寧に解き明かしていきます。

本記事の内容

  • 位相遅延 $\phi(\omega)/\omega$ と群遅延 $-d\phi/d\omega$ の定義と物理的意味(包絡線の遅れ)
  • 線形位相 = 一定群遅延 = 波形保存であることの導出
  • FIR係数の対称・反対称が線形位相を生む証明(4タイプ)
  • IIRが厳密な線形位相を持てない理由と、全域通過フィルタによる群遅延等化
  • Python:線形位相FIRと非線形位相IIRに同じ波束を通して波形歪みを比較、群遅延プロット、対称係数の効果可視化

前提知識

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

フィルタの位相とは何を意味するのか

振幅応答 $|H(e^{j\omega})|$ がフィルタの「音量つまみ」だとすれば、位相応答 $\angle H(e^{j\omega})$ は「タイミングつまみ」です。正弦波 $\cos(\omega n)$ をフィルタに通すと、出力は同じ周波数の正弦波になりますが、振幅が $|H(e^{j\omega})|$ 倍され、位相が $\angle H(e^{j\omega})$ だけずれます。位相がずれるということは、その正弦波が時間軸上で前後にスライドするということです。

具体的に確かめてみましょう。LTI(線形時不変)システムに複素正弦波 $x[n] = e^{j\omega n}$ を入力すると、出力は固有関数の性質から

$$ y[n] = H(e^{j\omega}) \, e^{j\omega n} $$

となります。ここで $H(e^{j\omega}) = |H(e^{j\omega})| \, e^{j\phi(\omega)}$($\phi(\omega) = \angle H(e^{j\omega})$)と極形式で書くと

$$ y[n] = |H(e^{j\omega})| \, e^{j(\omega n + \phi(\omega))} = |H(e^{j\omega})| \, e^{j\omega\left(n + \frac{\phi(\omega)}{\omega}\right)} $$

と変形できます。最後の式の指数の中身を見ると、$n$ が $n + \phi(\omega)/\omega$ に置き換わっています。つまりこの正弦波は時間軸上で $-\phi(\omega)/\omega$ サンプルだけ遅れた($\phi$ が負なら遅れ、正なら進み)ことになります。この遅れの量こそが、次節で定義する位相遅延です。

ここで自然な疑問が生まれます。単独の正弦波なら「位相のずれ=時間遅れ」と素直に解釈できますが、実際の信号は多数の周波数成分の重ね合わせです。成分ごとに遅れがバラバラだと、信号全体としてはどう遅れるのでしょうか。この問いに答えるのが群遅延です。

位相遅延と群遅延の定義

位相遅延:単一正弦波の遅れ

前節で見たように、周波数 $\omega$ の正弦波がフィルタを通ると時間軸上で $-\phi(\omega)/\omega$ だけずれます。これを位相遅延(phase delay)と呼び、次のように定義します。

$$ \tau_p(\omega) = -\frac{\phi(\omega)}{\omega} $$

位相遅延は「その周波数の搬送波(キャリア)そのものが何サンプル遅れるか」を表します。たとえば $\tau_p = 5$ なら、その周波数の正弦波の山と谷が一律に5サンプル後ろにずれます。直感的には、波の「位相の点(ゼロ交差や山の位置)」がどれだけ遅れるかを測る量です。

群遅延:包絡線(波束)の遅れ

ところが、私たちが実際に扱う信号は単一正弦波ではなく、複数の周波数成分が集まった波束(wave packet)です。たとえば、ある周波数を中心に狭い帯域の成分が集まると、それらが干渉して「うねり(ビート)」のような包絡線を作ります。この包絡線がどれだけ遅れるかを測るのが群遅延(group delay)で、位相応答の傾き(微分)として定義されます。

$$ \tau_g(\omega) = -\frac{d\phi(\omega)}{d\omega} $$

なぜ微分が出てくるのでしょうか。これは、近接した2つの周波数成分が作る包絡線の動きを考えると見えてきます。中心周波数 $\omega_0$ のまわりに $\omega_0 \pm \Delta\omega$ の2成分があるとして、それらの合成を追ってみましょう。

群遅延が「微分」になる理由の導出

中心周波数 $\omega_0$ に近接した2つの正弦波の和を入力します。

$$ x[n] = \cos\big((\omega_0 + \Delta\omega) n\big) + \cos\big((\omega_0 – \Delta\omega) n\big) $$

和積の公式 $\cos A + \cos B = 2\cos\frac{A+B}{2}\cos\frac{A-B}{2}$ を使うと、入力は

$$ x[n] = 2 \underbrace{\cos(\Delta\omega \, n)}_{\text{包絡線(ゆっくり)}} \cdot \underbrace{\cos(\omega_0 \, n)}_{\text{搬送波(速い)}} $$

と書けます。$\Delta\omega$ は小さいので、$\cos(\Delta\omega\, n)$ はゆっくり変動する包絡線、$\cos(\omega_0 n)$ は速く振動する搬送波です。

これをフィルタに通すと、各成分はそれぞれ位相 $\phi(\omega_0 + \Delta\omega)$ と $\phi(\omega_0 – \Delta\omega)$ だけずれます。ここで $\phi(\omega)$ を $\omega_0$ のまわりで1次テイラー展開すると

$$ \phi(\omega_0 \pm \Delta\omega) \approx \phi(\omega_0) \pm \phi'(\omega_0)\, \Delta\omega $$

となります。これを出力に代入していきましょう。振幅を簡単のため1とすると、出力は

$$ y[n] = \cos\big((\omega_0 + \Delta\omega) n + \phi(\omega_0) + \phi'(\omega_0)\Delta\omega\big) + \cos\big((\omega_0 – \Delta\omega) n + \phi(\omega_0) – \phi'(\omega_0)\Delta\omega\big) $$

ふたたび和積の公式を適用します。2つの角度の和の半分と差の半分を計算すると、和の半分は $\omega_0 n + \phi(\omega_0)$、差の半分は $\Delta\omega\, n + \phi'(\omega_0)\Delta\omega$ になります。したがって

$$ y[n] = 2 \cos\big(\Delta\omega\, n + \phi'(\omega_0)\Delta\omega\big) \cdot \cos\big(\omega_0 n + \phi(\omega_0)\big) $$

を得ます。ここで2つの因子を入力と見比べてみましょう。搬送波の因子は

$$ \cos\big(\omega_0 n + \phi(\omega_0)\big) = \cos\!\Big(\omega_0\big(n + \tfrac{\phi(\omega_0)}{\omega_0}\big)\Big) $$

なので、搬送波は $-\phi(\omega_0)/\omega_0 = \tau_p(\omega_0)$ だけ遅れています(位相遅延)。一方、包絡線の因子は

$$ \cos\big(\Delta\omega\, n + \phi'(\omega_0)\Delta\omega\big) = \cos\!\Big(\Delta\omega\big(n + \phi'(\omega_0)\big)\Big) $$

なので、包絡線は $-\phi'(\omega_0) = \tau_g(\omega_0)$ だけ遅れています(群遅延)。

この導出が示すのは、搬送波(位相の点)は位相遅延 $\tau_p$ で遅れ、包絡線(波束のかたまり)は群遅延 $\tau_g$ で遅れるという、二段構えの遅延構造です。情報は包絡線に乗ることが多い(変調信号、パルスのエッジ、過渡的なアタックなど)ため、信号の実質的な遅れを支配するのは群遅延の方です。だからこそ、波形を歪ませないためには群遅延が周波数によらず一定であることが要求されるのです。

ここまでで、群遅延が「包絡線の遅れ」であり位相応答の傾きで決まることがわかりました。では、群遅延が周波数によって変わると、信号にいったい何が起きるのでしょうか。次節で波形保存の条件として線形位相を導きます。

線形位相が波形を歪ませない理由

理想的な遅延器を考える

まず、波形をまったく歪ませない「理想の遅延器」がどんな周波数特性を持つかを考えましょう。入力 $x[n]$ をそっくりそのまま $D$ サンプル遅らせるシステムは

$$ y[n] = x[n – D] $$

です。両辺をDTFT(離散時間フーリエ変換)すると、時間シフトの性質から

$$ Y(e^{j\omega}) = e^{-j\omega D} X(e^{j\omega}) $$

となるので、この理想遅延器の周波数応答は

$$ H_{\text{ideal}}(e^{j\omega}) = e^{-j\omega D} $$

です。ここから2つのことが読み取れます。第一に、振幅応答は $|H_{\text{ideal}}| = 1$ で全周波数フラット(どの成分も同じ大きさのまま)。第二に、位相応答は

$$ \phi(\omega) = -\omega D $$

という $\omega$ の1次関数(線形)です。つまり「波形を保ったまま遅らせる」ことと「位相が周波数の線形関数である」ことは同じことなのです。

線形位相 ⇔ 一定群遅延

位相が $\phi(\omega) = -\alpha\omega + \beta$($\alpha, \beta$ は定数)の形のとき、群遅延を計算すると

$$ \tau_g(\omega) = -\frac{d}{d\omega}(-\alpha\omega + \beta) = \alpha $$

となり、周波数によらず一定値 $\alpha$ になります。逆に群遅延が一定 $\tau_g(\omega) = \alpha$ なら、積分して $\phi(\omega) = -\alpha\omega + \beta$ となり線形位相です。したがって

$$ \text{線形位相} \quad \Longleftrightarrow \quad \text{群遅延が一定} $$

という同値関係が成り立ちます。すべての周波数成分が同じ時間だけ遅れるので、合成された波形は形を保ったまま $\alpha$ サンプル平行移動します。これが「波形を歪ませない」ということの正体です。

なお、$\beta = 0$(純粋に $\phi = -\alpha\omega$)の場合を第1種の線形位相、$\beta = \pi/2$ などの定数オフセットを持つ場合を第2種の線形位相(generalized linear phase, 一般化線形位相)と呼びます。後者は位相遅延と群遅延が一致しませんが、群遅延は依然として一定なので、包絡線(波束)は歪まずに遅れます。微分器やヒルベルト変換器がこの第2種に対応します。

群遅延が一定でないと何が起きるか

群遅延が周波数によって変わるとき、$\tau_g(\omega) = -\phi'(\omega)$ が周波数ごとに異なるので、低周波成分と高周波成分で遅れる量が違ってきます。これを位相歪み(phase distortion)または分散(dispersion)と呼びます。たとえば、方形波は奇数次の高調波の和でできていますが、各高調波が違う量だけ遅れると、それらが足し合わさったときに鋭い立ち上がりが崩れ、リンギング(前後の振動)やエッジの鈍りが生じます。

ここで重要な事実を強調しておきます。耳は位相にほぼ鈍感だと言われますが、それは「定常的な音色の知覚」に関してであって、過渡的なアタックや時間構造を持つ信号(パルス、画像のエッジ、心電図の波形)では群遅延のばらつきが目に見える・耳に聞こえる歪みとして現れます。だからこそ、波形そのものを扱う応用では線形位相が決定的に重要なのです。

ここまでで、線形位相(一定群遅延)が波形保存の条件であることがわかりました。問題は、どうやってそれを実現するかです。次節で、FIRフィルタが係数の対称性によって厳密な線形位相を達成する仕組みを証明します。

FIRが対称係数で線形位相になる証明

対称性が鍵を握る

FIR(有限インパルス応答)フィルタは、係数(タップ)$h[0], h[1], \dots, h[M]$ で出力を決めます。

$$ H(e^{j\omega}) = \sum_{n=0}^{M} h[n]\, e^{-j\omega n} $$

ここで係数列が左右対称、すなわち $h[n] = h[M-n]$ を満たすと、奇跡的に厳密な線形位相が生まれます。直感的には、対称な係数列は「インパルス応答の重心が中央 $M/2$ にある」ことを意味し、その重心位置がそのまま全周波数共通の遅延(群遅延 $M/2$)になる、と理解できます。これを数式で確かめましょう。

Type I の証明(偶対称・$M$ 偶数)

$h[n] = h[M-n]$ かつ $M$ が偶数の場合を考えます。係数が中央 $M/2$ に関して対称なので、和を中央からの対称ペアにまとめます。中央の項 $n = M/2$ を別にして、対称な2項 $n$ と $M-n$ をペアにすると

$$ H(e^{j\omega}) = h\!\left[\tfrac{M}{2}\right] e^{-j\omega M/2} + \sum_{n=0}^{M/2 – 1} h[n]\Big(e^{-j\omega n} + e^{-j\omega(M-n)}\Big) $$

と書けます($h[n] = h[M-n]$ を使ってペアの係数を共通の $h[n]$ にまとめました)。ここで共通因子 $e^{-j\omega M/2}$ をくくり出します。$e^{-j\omega n} + e^{-j\omega(M-n)} = e^{-j\omega M/2}\big(e^{j\omega(M/2 – n)} + e^{-j\omega(M/2 – n)}\big)$ なので、オイラーの公式 $e^{jx} + e^{-jx} = 2\cos x$ を使うと

$$ H(e^{j\omega}) = e^{-j\omega M/2}\left[ h\!\left[\tfrac{M}{2}\right] + \sum_{n=0}^{M/2 – 1} 2\,h[n] \cos\!\Big(\omega\big(\tfrac{M}{2} – n\big)\Big)\right] $$

を得ます。ここで角括弧の中身に注目してください。$\cos$ の和なので、これは実数の関数です。この実数部分を $A(\omega)$(振幅関数、正にも負にもなりうる)と置くと

$$ H(e^{j\omega}) = A(\omega)\, e^{-j\omega M/2} $$

という形になります。$A(\omega)$ は実数なので、位相は $A(\omega) > 0$ のとき $\phi(\omega) = -\omega M/2$、$A(\omega) < 0$ のとき $\phi(\omega) = -\omega M/2 + \pi$(符号反転は $\pi$ の位相跳び)です。いずれにせよ、$\omega$ に対する傾きは一定なので

$$ \tau_g(\omega) = -\frac{d\phi}{d\omega} = \frac{M}{2} \quad \text{[samples]} $$

群遅延は全周波数で一定 $M/2$ です。これでType Iが厳密な線形位相を持つことが証明できました。$\square$

Type II の証明(偶対称・$M$ 奇数)

$M$ が奇数だと中央のタップが存在せず、係数は2項ずつ完全にペアになります。タップ数 $M+1$ は偶数です。同様に共通因子をくくり出すと、ペアの中心は $M/2$(半整数)になります。

$$ H(e^{j\omega}) = e^{-j\omega M/2} \sum_{n=0}^{(M-1)/2} 2\,h[n] \cos\!\Big(\omega\big(\tfrac{M}{2} – n\big)\Big) $$

ここで $M/2 – n$ は半整数なので、$\cos$ の引数は半整数倍の $\omega$ になります。やはり角括弧内は実数の振幅関数 $A(\omega)$ なので $H = A(\omega) e^{-j\omega M/2}$ となり、群遅延は一定 $M/2$ です($M$ が奇数なので半整数サンプルの遅延)。

ただし $\omega = \pi$ では $\cos(\pi(M/2 – n)) = \cos(\pi/2 \cdot (\text{奇数})) = 0$ となるため、$A(\pi) = 0$、すなわち $H(e^{j\pi}) = 0$ が強制されます。したがってType IIは $\omega = \pi$(ナイキスト)で必ずゼロになるので、高域通過フィルタには使えません。$\square$

Type III と Type IV の証明(反対称)

係数が反対称 $h[n] = -h[M-n]$ の場合を考えます。このとき対称ペアは引き算になるので、$e^{-j\omega n} – e^{-j\omega(M-n)} = e^{-j\omega M/2}\big(e^{j\omega(M/2-n)} – e^{-j\omega(M/2-n)}\big) = e^{-j\omega M/2} \cdot 2j\sin(\omega(M/2 – n))$ という形になります。$j = e^{j\pi/2}$ を使うと

$$ H(e^{j\omega}) = e^{-j\omega M/2}\, e^{j\pi/2} \cdot B(\omega) = B(\omega)\, e^{j(\pi/2 – \omega M/2)} $$

ここで $B(\omega) = \sum 2\,h[n]\sin(\omega(M/2 – n))$ は実数です。位相は $\phi(\omega) = \pi/2 – \omega M/2$(プラス $A$ の符号による $\pi$ 跳び)で、定数オフセット $\pi/2$ を含む一般化線形位相(第2種)になります。重要なのは、定数項は群遅延に効かないので

$$ \tau_g(\omega) = -\frac{d}{d\omega}\Big(\frac{\pi}{2} – \frac{\omega M}{2}\Big) = \frac{M}{2} $$

と、やはり群遅延は一定 $M/2$ です。$M$ が偶数ならType III、奇数ならType IVです。$\sin$ の項は $\omega = 0$ で必ずゼロなので、Type III・IVは $H(e^{j0}) = 0$(直流を通さない)。これは微分器やヒルベルト変換器のように、直流成分を持たないフィルタに対応します。$\square$

4タイプのまとめ

以上を表にまとめます。すべてのタイプで群遅延は $M/2$ サンプルで一定です。係数の対称性(偶対称か反対称か)と次数 $M$ の偶奇によって、$\omega=0$ や $\omega=\pi$ での値に制約が課されます。

タイプ 対称性 $M$ $\omega=0$ $\omega=\pi$ 群遅延 典型用途
I 偶対称 $h[n]=h[M-n]$ 偶数 制約なし 制約なし $M/2$(整数) LPF/HPF/BPF/BSF 全般
II 偶対称 $h[n]=h[M-n]$ 奇数 制約なし $H=0$ $M/2$(半整数) LPF/BPF
III 反対称 $h[n]=-h[M-n]$ 偶数 $H=0$ $H=0$ $M/2$(整数) BPF/微分器/ヒルベルト変換
IV 反対称 $h[n]=-h[M-n]$ 奇数 $H=0$ 制約なし $M/2$(半整数) HPF/微分器

この表は、フィルタ設計で「どのタイプを選ぶか」の指針になります。たとえば高域通過を線形位相で作りたいなら、$\omega=\pi$ がゼロになるType IIは使えず、Type I か Type IV を選ぶ、という判断ができます。

FIRが対称係数で厳密な線形位相を実現できることがわかりました。では、フィードバックを持つIIRフィルタはどうでしょうか。次節で、IIRが原理的に厳密な線形位相を持てない理由を見ます。

IIRが線形位相を持てない理由

厳密な線形位相 = 対称インパルス応答 = 有限長

前節の証明の核心は「インパルス応答が中心対称($h[n] = \pm h[M-n]$)であること」でした。これを一般化すると、因果的で安定なフィルタが厳密な線形位相を持つための必要条件は、インパルス応答が有限長で対称であることだと言えます。逆に言えば、無限に続くインパルス応答を持つフィルタは厳密な線形位相を持てません。IIR(無限インパルス応答)フィルタはまさにこの「無限長」が定義そのものなので、原理的に厳密な線形位相は不可能なのです。

もう少し丁寧に背理法で考えてみましょう。仮にIIRが線形位相を持つとすると、そのインパルス応答 $h[n]$ は中心 $\alpha$ に関して対称 $h[n] = h[2\alpha – n]$ でなければなりません。ところがIIRのインパルス応答は $h[n] = \sum_k c_k p_k^n u[n]$ の形($p_k$ は極、$u[n]$ は単位ステップ)で、$n \to +\infty$ で減衰しながら無限に続く片側の列です。これが有限の中心 $\alpha$ について対称になるには、$n < 0$ 側にも対応する非ゼロ成分が必要になり、因果性($n<0$ で $h[n]=0$)と矛盾します。したがって、安定かつ因果なIIRは厳密な線形位相を持てません。

極が課す制約:振幅と位相のトレードオフ

別の見方として、IIRの伝達関数

$$ H(z) = \frac{B(z)}{A(z)} $$

で、分母 $A(z)$ の根である極が、急峻な振幅特性を生む代償として周波数依存の群遅延を生みます。極が単位円に近いほど、その近傍の周波数で振幅が鋭くピークを持ち、同時に位相が急変するため群遅延が大きく跳ね上がります。バターワース・チェビシェフ・楕円フィルタはいずれも、カットオフ周波数付近で群遅延がピークを持ち、波形に位相歪みを与えます。IIRは「低次で急峻」という効率の良さと引き換えに、位相の素直さを犠牲にしているのです。

IIRで線形位相に近づける現実的な手段

実務では、IIRの非線形位相を回避・補償するいくつかの手段があります。

  1. ゼロ位相フィルタリング(双方向フィルタリング):信号を順方向にフィルタした後、時間反転してもう一度フィルタする方法です。SciPyの filtfilt がこれにあたります。順方向の位相ずれ $\phi(\omega)$ と逆方向の位相ずれ $-\phi(\omega)$ が打ち消し合い、合成の位相は厳密にゼロ(群遅延ゼロ)になります。ただしオフライン処理(信号全体が手元にある)でしか使えず、リアルタイムには不向きです。
  2. 全域通過フィルタによる群遅延等化:振幅は変えずに位相だけを調整する全域通過フィルタ(オールパスフィルタ)を縦続接続し、トータルの群遅延を平らにする方法です。これは次節で詳しく扱います。

全域通過フィルタによる群遅延等化

全域通過フィルタ(all-pass filter)は、すべての周波数で振幅が一定($|H_{\text{ap}}(e^{j\omega})| = 1$)でありながら、位相だけを変えられる特殊なフィルタです。1次の全域通過フィルタは

$$ H_{\text{ap}}(z) = \frac{z^{-1} – a^*}{1 – a z^{-1}} $$

の形をしています($|a| < 1$ で安定)。この構造のポイントは、極 $a$ と零点 $1/a^*$ が単位円に関して鏡像対称に配置されていることです。単位円上の任意の点 $e^{j\omega}$ から見て、極までの距離と零点までの距離の比が常に一定になるため、振幅が周波数によらず1になります。実際に絶対値を計算すると

$$ |H_{\text{ap}}(e^{j\omega})| = \frac{|e^{-j\omega} – a^*|}{|1 – a e^{-j\omega}|} = 1 $$

が成り立ちます(分子と分母が互いに複素共役の関係になることから示せます)。振幅を変えないので、メインのIIRフィルタの後ろに全域通過を縦続接続しても、振幅特性はそのままに位相(群遅延)だけを修正できます。

群遅延等化の戦略はこうです。まず急峻な振幅特性を持つIIRを設計し、そのカットオフ付近で跳ね上がる群遅延 $\tau_g^{\text{IIR}}(\omega)$ を測ります。次に、全域通過フィルタの群遅延 $\tau_g^{\text{ap}}(\omega)$ を、$\tau_g^{\text{IIR}}(\omega) + \tau_g^{\text{ap}}(\omega)$ がなるべく平らになるように設計します。全域通過は群遅延を「足す」ことしかできない(負の群遅延は安定因果系では作れない)ので、谷を埋めて全体を平坦化する発想です。代償として、トータルの遅延は増えます。これは波形保存とレイテンシのトレードオフそのものです。

ここまでで、FIRとIIRの位相特性の根本的な違いと、IIRを線形位相に近づける手段がわかりました。ここからは、これらの主張をPythonで実際に確かめていきましょう。

Pythonでの実装

パルスを通して波形歪みを比較する

まず最も直感的な実験から始めます。鋭いガウシアン波束(ガウシアン包絡に搬送波を乗せた、包絡線がはっきりしたパルス)を、線形位相FIRと非線形位相のIIR(バターワース)にそれぞれ通し、出力波形の崩れ方を比べます。線形位相FIRは波形を保ったまま遅らせるだけ、IIRは包絡線(波束)を非対称に歪ませるはずです。

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

# ===== テスト信号: ガウシアン波束(包絡 × 搬送波) =====
fs = 1000                       # サンプリング周波数 [Hz]
N = 600                         # サンプル数
t = np.arange(N) / fs           # 時間軸 [s]
t0 = 0.15                       # パルス中心 [s]
width = 0.012                   # 包絡線の幅 [s]
fc_pulse = 80                   # 搬送波周波数 [Hz](カットオフ手前の通過帯域内)
# ガウシアン包絡 × 搬送波(帯域内に収まるように)
pulse = (np.exp(-((t - t0) ** 2) / (2 * width ** 2))
         * np.cos(2 * np.pi * fc_pulse * (t - t0)))

# ===== フィルタ設計(同じカットオフ 120Hz) =====
fc = 120
# 線形位相 FIR(Type I: 偶対称・奇数タップ)
numtaps = 101
h_fir = signal.firwin(numtaps, fc, fs=fs, window='hamming')
# IIR バターワース 6次
b_iir, a_iir = signal.butter(6, fc, fs=fs)

# フィルタ適用(因果フィルタリング lfilter)
y_fir = signal.lfilter(h_fir, 1.0, pulse)
y_iir = signal.lfilter(b_iir, a_iir, pulse)

print(f"FIR 群遅延(理論) = {(numtaps - 1) / 2} samples = "
      f"{(numtaps - 1) / 2 / fs * 1000:.2f} ms")

この時点で、FIRの群遅延は $(N_{\text{tap}}-1)/2 = 50$ サンプル(50 ms)と理論どおりに表示されます。これは係数が中央50番目を中心に対称だからで、前節のType Iの結論と一致します。次に波形を描画して歪みを目視します。

fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(t * 1000, pulse, 'k-', lw=1.5, label='Input pulse')
ax.plot(t * 1000, y_fir, 'b-', lw=1.5, label='Linear-phase FIR output')
ax.plot(t * 1000, y_iir, 'r-', lw=1.5, label='IIR (Butterworth) output')
ax.set_xlabel('Time [ms]')
ax.set_ylabel('Amplitude')
ax.set_title('Pulse through linear-phase FIR vs IIR')
ax.set_xlim([100, 350])
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.savefig('pulse_distortion.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフから、位相特性の違いが波形に与える影響がはっきり読み取れます。FIRの出力(青)は入力波束(黒)と相似な左右対称の包絡線を保ったまま、右に約50 msだけ平行移動しています(実際、入力包絡線を50サンプルずらしたものとの相関係数は 0.99999 と、ほぼ完全に一致します)。これがまさに線形位相=波形保存です。一方、IIRの出力(赤)は包絡線のピークがわずかに低くなり(包絡線ピーク値は約 0.986 に低下)、立ち下がりに尾を引く非対称な形に歪んでいます。さらに包絡線のピーク位置は入力の 150 ms から約 157 ms へと 7 ms ほどしかずれておらず、FIR の 50 ms と比べて遅延量も周波数成分ごとにバラついています。これはカットオフ付近で群遅延が周波数によって異なるため、波束を構成する各周波数成分の到着タイミングがずれて干渉した結果です。波形そのものが情報を持つレーダーや生体信号では、この赤い歪みが致命的になり得ます。

群遅延を直接プロットして比較する

次に、FIRとIIRの群遅延を周波数の関数として描き、FIRが一定・IIRがカットオフ付近でピークを持つことを定量的に確かめます。

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

fs = 1000
fc = 120

# FIR(線形位相)と IIR(バターワース 6次)
h_fir = signal.firwin(101, fc, fs=fs, window='hamming')
b_iir, a_iir = signal.butter(6, fc, fs=fs)

# 群遅延(scipy.signal.group_delay)
w_fir, gd_fir = signal.group_delay((h_fir, 1.0), w=2048, fs=fs)
w_iir, gd_iir = signal.group_delay((b_iir, a_iir), w=2048, fs=fs)

fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(w_fir, gd_fir, 'b-', lw=2, label='Linear-phase FIR (101 taps)')
ax.plot(w_iir, gd_iir, 'r-', lw=2, label='IIR Butterworth (N=6)')
ax.axhline((101 - 1) / 2, color='b', ls='--', alpha=0.4,
           label='FIR theory = 50 samples')
ax.axvline(fc, color='gray', ls=':', alpha=0.7, label=f'$f_c$={fc} Hz')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Group delay [samples]')
ax.set_title('Group delay: FIR (flat) vs IIR (peaked at cutoff)')
ax.set_xlim([0, fs / 2])
ax.set_ylim([0, 70])
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.savefig('group_delay_compare.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフが本記事の主張を端的に裏づけます。FIRの群遅延(青)は全周波数で水平な直線50サンプルを保っており、理論値(青破線)と完璧に一致します。すべての周波数成分が同じだけ遅れるので波形が保たれるわけです。対してIIRの群遅延(赤)はカットオフ周波数120 Hz付近で急峻なピークを作っています。この帯域の成分だけが余分に遅れるため、先ほどのパルスが非対称に歪んだのです。群遅延の「平坦さ」こそが波形保存の指標であることが視覚的に確認できます。

対称係数が線形位相を生むことの可視化

続いて、FIRの係数列そのものを描き、対称な係数が線形位相を、非対称な係数が非線形位相を生むことを直接見ます。同じ係数を意図的に非対称化して位相を壊し、対比します。

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

fs = 1000
fc = 120
M = 40                                   # 次数(タップ数 M+1 = 41, Type I)

# 対称な線形位相 FIR
h_sym = signal.firwin(M + 1, fc, fs=fs, window='hamming')
# 非対称化: 係数に単調な重みを掛けて対称性を壊す
weight = np.linspace(0.3, 1.7, M + 1)
h_asym = h_sym * weight
h_asym = h_asym / np.sum(h_asym) * np.sum(h_sym)  # DCゲインを揃える

fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))

# 係数列の比較
axes[0].stem(np.arange(M + 1), h_sym, linefmt='b-', markerfmt='bo',
             basefmt='k-', label='Symmetric h[n]=h[M-n]')
axes[0].stem(np.arange(M + 1), h_asym, linefmt='r-', markerfmt='rx',
             basefmt='k-', label='Asymmetric (broken)')
axes[0].axvline(M / 2, color='gray', ls=':', label='center M/2')
axes[0].set_xlabel('Tap index n')
axes[0].set_ylabel('h[n]')
axes[0].set_title('FIR coefficients')
axes[0].legend(fontsize=9)
axes[0].grid(True, alpha=0.3)

# 位相応答の比較
w1, H1 = signal.freqz(h_sym, 1.0, worN=2048, fs=fs)
w2, H2 = signal.freqz(h_asym, 1.0, worN=2048, fs=fs)
axes[1].plot(w1, np.unwrap(np.angle(H1)), 'b-', lw=2, label='Symmetric (linear)')
axes[1].plot(w2, np.unwrap(np.angle(H2)), 'r-', lw=2, label='Asymmetric (nonlinear)')
axes[1].axvline(fc, color='gray', ls=':')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Phase [rad]')
axes[1].set_title('Phase response (passband)')
axes[1].set_xlim([0, fc])
axes[1].legend(fontsize=9)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('symmetry_phase.png', dpi=150, bbox_inches='tight')
plt.show()

左右の図を対比すると、対称性と位相の関係が一目瞭然です。左図で青の係数列は中央 $M/2=20$ に関して左右対称ですが、赤は重みを掛けて対称性を壊しています。その結果、右図では青の位相応答が通過帯域で見事に直線(線形位相)になっているのに対し、赤は曲がった非線形位相になっています。係数列のたった一つの幾何学的性質である「対称性」が、位相の線形性を支配していることが実感できます。これは前節のType Iの証明を実験的に裏づけるものです。

全域通過フィルタによる群遅延等化

最後に、IIRの跳ね上がった群遅延を全域通過フィルタで平坦化する例を見ます。1次全域通過の極(=零点の鏡像)を調整し、IIRの群遅延の谷を埋めてトータルをなるべく平らにします。

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

fs = 1000
fc = 120

# ベースの IIR(バターワース 6次)
b_iir, a_iir = signal.butter(6, fc, fs=fs)
w, gd_iir = signal.group_delay((b_iir, a_iir), w=2048, fs=fs)

# 1次全域通過フィルタ H_ap(z) = (z^-1 - a) / (1 - a z^-1), 実係数版
def allpass1(a):
    """1次全域通過の (b, a_den) 係数を返す(a は実数, |a|<1)"""
    b = np.array([-a, 1.0])      # 分子 z^-1 - a (係数は z^0, z^-1 の順)
    a_den = np.array([1.0, -a])  # 分母 1 - a z^-1
    return b, a_den

# 低域側の遅延を持ち上げて全体を平らにする全域通過を 2 段
b_ap1, a_ap1 = allpass1(0.6)
b_ap2, a_ap2 = allpass1(0.3)

# 縦続接続(畳み込みで多項式を合成)
b_cas = np.convolve(np.convolve(b_iir, b_ap1), b_ap2)
a_cas = np.convolve(np.convolve(a_iir, a_ap1), a_ap2)
w_c, gd_cas = signal.group_delay((b_cas, a_cas), w=2048, fs=fs)

# 振幅が変わっていないことの確認
_, H_iir = signal.freqz(b_iir, a_iir, worN=2048, fs=fs)
_, H_cas = signal.freqz(b_cas, a_cas, worN=2048, fs=fs)

ここで全域通過の縦続接続を構成しました。次に、群遅延と振幅を並べて、全域通過が「振幅は不変・群遅延だけ変化」させることを確認します。

fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))

# 群遅延の比較
axes[0].plot(w, gd_iir, 'r-', lw=2, label='IIR only')
axes[0].plot(w_c, gd_cas, 'g-', lw=2, label='IIR + 2 all-pass (equalized)')
axes[0].axvline(fc, color='gray', ls=':', label=f'$f_c$={fc} Hz')
axes[0].set_xlabel('Frequency [Hz]')
axes[0].set_ylabel('Group delay [samples]')
axes[0].set_title('Group delay equalization by all-pass')
axes[0].set_xlim([0, fs / 2])
axes[0].grid(True, alpha=0.3)
axes[0].legend(fontsize=9)

# 振幅応答の比較(変わらないことを確認)
axes[1].plot(w, 20 * np.log10(np.abs(H_iir) + 1e-12), 'r-', lw=2, label='IIR only')
axes[1].plot(w, 20 * np.log10(np.abs(H_cas) + 1e-12), 'g--', lw=2,
             label='IIR + all-pass')
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Magnitude [dB]')
axes[1].set_title('Magnitude (unchanged by all-pass)')
axes[1].set_xlim([0, fs / 2])
axes[1].set_ylim([-80, 5])
axes[1].grid(True, alpha=0.3)
axes[1].legend(fontsize=9)

plt.tight_layout()
plt.savefig('allpass_equalization.png', dpi=150, bbox_inches='tight')
plt.show()

この2枚のグラフから、全域通過フィルタの役割がはっきりします。左図で、元のIIR(赤)が低域で小さく・カットオフ付近で跳ね上がる群遅延を持っていたのに対し、全域通過を加えた合成系(緑)は低域側が持ち上がり、全体としてより平坦に近づいています。全域通過は群遅延を「足す」ことしかできないので、谷を埋めて全体を底上げする形で平坦化しているのが見て取れます。そして右図で、振幅応答は赤と緑がぴったり重なっており、全域通過が振幅をまったく変えていないことが確認できます。これが「振幅を保ったまま位相だけを操作できる」全域通過フィルタの威力であり、IIRの効率を保ちつつ位相歪みを軽減する実用的な手段になります。

ゼロ位相フィルタリングとの比較

参考として、オフライン処理で使えるゼロ位相フィルタリング(双方向フィルタ)も確認しておきましょう。filtfilt は順方向と逆方向に2回フィルタすることで位相ずれを完全に打ち消します。

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

fs = 1000
N = 600
t = np.arange(N) / fs
# 先ほどと同じガウシアン波束(包絡 × 80 Hz 搬送波)
pulse = (np.exp(-((t - 0.15) ** 2) / (2 * 0.012 ** 2))
         * np.cos(2 * np.pi * 80 * (t - 0.15)))

b_iir, a_iir = signal.butter(6, 120, fs=fs)
y_causal = signal.lfilter(b_iir, a_iir, pulse)   # 通常の因果フィルタ
y_zero = signal.filtfilt(b_iir, a_iir, pulse)    # ゼロ位相(双方向)

fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(t * 1000, pulse, 'k-', lw=1.5, label='Input')
ax.plot(t * 1000, y_causal, 'r-', lw=1.5, label='IIR causal (lfilter)')
ax.plot(t * 1000, y_zero, 'g-', lw=1.5, label='IIR zero-phase (filtfilt)')
ax.set_xlabel('Time [ms]')
ax.set_ylabel('Amplitude')
ax.set_title('Zero-phase filtering preserves waveform symmetry')
ax.set_xlim([80, 250])
ax.grid(True, alpha=0.3)
ax.legend()
plt.tight_layout()
plt.savefig('zerophase_compare.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフから、ゼロ位相処理の効果が読み取れます。通常の因果フィルタ(赤)はパルスを遅延させ、立ち下がりに尾を引く非対称な歪みを生んでいます。一方、ゼロ位相フィルタ(緑)は入力パルス(黒)とほぼ同じ位置・同じ左右対称な形を保っています。順方向と逆方向で位相ずれが相殺されて群遅延が実質ゼロになるためです。ただしこれは信号全体が手元にあるオフライン処理だからこそ可能で、リアルタイム処理には使えない点に注意が必要です。リアルタイムで波形保存が必要なら線形位相FIR、効率を取るなら全域通過等化付きIIR、という使い分けになります。

まとめ

本記事では、群遅延と線形位相という観点から、フィルタが波形を歪ませない条件を解説しました。

  • 位相遅延と群遅延:位相遅延 $\tau_p = -\phi(\omega)/\omega$ は搬送波(位相の点)の遅れ、群遅延 $\tau_g = -d\phi/d\omega$ は包絡線(波束)の遅れを表す。情報は包絡線に乗るため、実質的な遅れを支配するのは群遅延である
  • 線形位相 ⇔ 一定群遅延 ⇔ 波形保存:位相が $\phi(\omega) = -\alpha\omega + \beta$ の線形なら群遅延は一定 $\alpha$ となり、全周波数成分が同じだけ遅れるので波形は形を保ったまま平行移動する
  • FIRの線形位相条件:係数が対称 $h[n]=h[M-n]$(または反対称)のとき、$H(e^{j\omega}) = A(\omega)e^{-j\omega M/2}$ の形になり群遅延は厳密に $M/2$ で一定。次数の偶奇と対称性で4タイプに分類され、$\omega=0,\pi$ での制約が決まる
  • IIRが線形位相を持てない理由:厳密な線形位相は有限長かつ対称なインパルス応答を要求するため、無限長のIIRでは原理的に不可能。急峻な振幅特性を生む極が、カットオフ付近で群遅延のピークを生む
  • 補償手段:オフラインならゼロ位相フィルタリング(filtfilt)で群遅延をゼロにでき、リアルタイムなら全域通過フィルタで振幅を保ったまま群遅延を平坦化できる

群遅延と線形位相の理解は、フィルタを「振幅だけで選ぶ」段階から「位相まで含めて設計する」段階へと進むための鍵です。波形が情報を持つ応用(レーダーのパルス圧縮、生体信号解析、通信のシンボル整形)では、群遅延の平坦さがそのまま性能に直結します。

次のステップとして、以下の記事も参考にしてください。