「同じタップ数のFIRフィルタなら、どう係数を決めれば一番性能が良くなるのか？」——これは信号処理の実務で必ずぶつかる問いです。窓関数法は手軽ですが、「ハミング窓よりブラックマン窓のほうが阻止域は深いが遷移帯域が広い」といったトレードオフを職人芸で選んでいるだけで、与えられたタップ数を最大限に使い切っているという保証はどこにもありません。実際、窓関数法で設計したフィルタの周波数応答を見ると、阻止域の最初のリップル(サイドローブ)だけが突出して大きく、奥のほうのリップルはずっと小さい、という「無駄」がしばしば見られます。リップルの山の高さがバラバラなのです。
ここで発想を変えます。「最大の誤差(一番高いリップル)を、できる限り小さく押し下げたい」。この目標を厳密に追求すると、驚くべき結論にたどり着きます。最適なフィルタでは、通過域と阻止域の誤差リップルがすべて同じ高さに揃う——つまり等リップル(equiripple)になるのです。この最適等リップルFIRフィルタを系統的に設計するのが、本記事で扱う Parks-McClellan法(その心臓部が Remezアルゴリズム)です。MATLABの firpm、SciPyの scipy.signal.remez として実装されており、急峻なフィルタを最小のタップ数で実現したいときの定番手法です。
この手法を理解する価値は大きく、応用先も広いです。第一に、通信機のチャネル選択フィルタやデジタル無線(SDR)では、限られた演算資源(FPGAのタップ数)で隣接チャネルを叩き落とす必要があり、タップ効率の良い等リップル設計が必須です。第二に、オーディオのグラフィックイコライザや計測機器のアンチエイリアシングフィルタでも、線形位相を保ったまま規格(通過域リップル ±0.1 dB、阻止域 −80 dB など)を最小コストで満たすために使われます。本記事では、「最大誤差最小化(ミニマックス)」という目標設定から出発し、近似論の華である チェビシェフの交番定理(alternation theorem) を述べ、それを反復的に解く Remez交換アルゴリズム の各ステップを数式で省略なく導出します。最後にPythonでアルゴリズムの内部動作(基準点の移動と誤差の交番)を可視化し、窓関数法FIRと振幅応答・リップル・遷移帯域を定量的に比較します。
本記事の内容
- ミニマックス(最大誤差最小化)というフィルタ設計の目標設定
- 線形位相FIRが余弦多項式に帰着すること(チェビシェフ近似への橋渡し)
- チェビシェフの交番定理(最良近似は誤差符号が交番する)の主張と直感
- Remez交換アルゴリズムの各ステップの導出(連立方程式・誤差評価・基準点更新)
- Pythonでのremez相当のスクラッチ実装と内部動作の可視化
- 窓関数法FIRとの振幅応答・リップル・遷移帯域の比較と考察
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
特に、線形位相FIRが係数の対称性($h[n]=h[M-n]$)から得られること、周波数応答が $H(e^{j\omega})$ で評価できることは前提とします。これらは上記「ディジタルフィルタの基礎」で詳しく扱っています。
フィルタ設計を「近似問題」として捉え直す
低域通過フィルタを設計するとき、私たちが本当に欲しいのは「理想フィルタ」です。理想低域通過フィルタの振幅応答 $D(\omega)$ は、通過域($0 \le \omega \le \omega_p$)で $1$、阻止域($\omega_s \le \omega \le \pi$)で $0$ という、角のある階段関数です。しかし有限タップのFIRフィルタでこの不連続な関数を完全に再現することはできません。できるのは「近づける」ことだけです。
ここで重要な視点の転換が起こります。FIRフィルタの設計とは、有限個の係数で表現できる関数 $A(\omega)$(実現できる振幅応答)を、理想応答 $D(\omega)$ にできる限り近づける近似問題だ、ということです。「できる限り近づける」をどう数学的に定義するかで、設計法が分かれます。
近さの測り方には大きく2通りあります。一つは誤差の二乗和(エネルギー)を最小化する 最小二乗(L2)近似。これは平均的には良いのですが、不連続点(通過域と阻止域の境界)の近くで誤差が局所的に大きく暴れる(ギブス現象)という弱点があります。もう一つが、誤差の最大値(L∞ノルム)を最小化する ミニマックス近似です。
$$ \min_{A} \ \max_{\omega \in F} \ W(\omega) \, \bigl| D(\omega) – A(\omega) \bigr| $$
ここで $F$ は対象とする周波数の集合(通過域と阻止域の和。遷移帯域は「どうでもよい」ので除外します)、$W(\omega)$ は重み関数で、通過域と阻止域でどちらの誤差を重く見るかを調整します。この「最悪のケースの誤差を最小にする」という考え方こそが Parks-McClellan 法の出発点です。最悪ケースを抑えるので、「一番高いリップルが規格に収まればOK」という実務の要求にそのまま合致します。
なぜミニマックスを選ぶと等リップルになるのか——その理由はまだ見えていません。それを説明するには、まず実現できる振幅応答 $A(\omega)$ がどんな数学的構造を持つのかを明らかにする必要があります。次節で、線形位相FIRの振幅応答が「余弦の多項式」になることを示し、近似論の土俵に乗せます。
線形位相FIRは余弦多項式に帰着する
ミニマックス近似を実行するには、近似関数 $A(\omega)$ の正体を掴まなければなりません。ここで線形位相FIRフィルタの構造が効いてきます。
長さ $N = 2L+1$(奇数タップ)の Type I 線形位相FIRを考えます。係数は中心 $L$ について偶対称、すなわち $h[n] = h[N-1-n]$ です。周波数応答は
$$ H(e^{j\omega}) = \sum_{n=0}^{N-1} h[n] e^{-j\omega n} $$
ですが、対称性を使って中心 $n=L$ のまわりにまとめます。$m = n – L$ と置き換えると、対称性 $h[L+m] = h[L-m]$ により、$\sin$ の項が打ち消し合って $\cos$ だけが残ります。途中を省略せず書くと、まず
$$ H(e^{j\omega}) = e^{-j\omega L} \sum_{m=-L}^{L} h[L+m] e^{-j\omega m} $$
と中心の線形位相 $e^{-j\omega L}$ を括り出します。$\sum_m h[L+m] e^{-j\omega m}$ の和で $m$ と $-m$ をペアにすると、$h[L+m]=h[L-m]$ なので
$$ h[L+m] e^{-j\omega m} + h[L-m] e^{j\omega m} = h[L+m]\bigl(e^{-j\omega m} + e^{j\omega m}\bigr) = 2 h[L+m] \cos(\omega m) $$
となります($m=0$ の項は $h[L]$ そのまま)。したがって
$$ H(e^{j\omega}) = e^{-j\omega L} \underbrace{\left[ h[L] + \sum_{m=1}^{L} 2 h[L+m] \cos(\omega m) \right]}_{\displaystyle A(\omega)} $$
と書けます。位相因子 $e^{-j\omega L}$ は純粋な線形位相(群遅延 $L$ サンプル一定)を表し、波形を歪ませません。残りの 実数値の振幅応答 $A(\omega)$ が、私たちが理想応答に近づけたい対象です。
ここで $a_0 = h[L]$、$a_k = 2 h[L+k]$($k=1,\dots,L$)と係数を置き直すと
$$ A(\omega) = \sum_{k=0}^{L} a_k \cos(k\omega) $$
という、余弦の線形結合になります。さらに、チェビシェフ多項式の関係式 $\cos(k\omega) = T_k(\cos\omega)$(ただし $T_k$ は第1種チェビシェフ多項式)を使うと、$x = \cos\omega$ と置けば
$$ A(\omega) = \sum_{k=0}^{L} a_k T_k(x) = P(x), \qquad x = \cos\omega \in [-1, 1] $$
となり、$A(\omega)$ は変数 $x=\cos\omega$ の 高々 $L$ 次の多項式 $P(x)$ にぴったり一致します。これは決定的に重要な事実です。なぜなら、「ある区間上で連続関数を多項式でミニマックス近似する」問題には、19世紀から研究され尽くした美しい理論——チェビシェフの近似理論——がそのまま使えるからです。
つまりフィルタ設計は、区間 $[-1,1]$ 上で目標関数 $D$ を $L$ 次多項式 $P$ で最良近似する問題に翻訳されました。次節では、その最良近似がどんな形をしているのかを述べる 交番定理 に進みます。これが「最適フィルタはなぜ等リップルなのか」の答えを与えてくれます。
チェビシェフの交番定理 — 最良近似は誤差が交番する
多項式近似の最良解がどんな姿をしているか、まず直感から入りましょう。あなたが $L$ 次多項式(自由度 $L+1$ 個のつまみ)を持っていて、目標曲線にぴったり沿わせたいとします。誤差曲線 $E(\omega) = W(\omega)[D(\omega) – A(\omega)]$ を考えると、もし誤差の山が一箇所だけ突出して高ければ、そこを少し押し下げる余地が必ずあります(つまみがまだ余っている)。逆に、押し下げると別の場所が持ち上がってしまい、どこを触っても最大値が下がらなくなった状態——それが最適です。そして最大誤差が動かせなくなるのは、プラス方向とマイナス方向の最大リップルが十分多くの点で交互に現れているときです。山と谷が交互にいくつも並んでいて、どれか一つを下げると別のどれかが上がる「シーソー」が成立している、というわけです。
これを厳密に述べたのが交番定理です。
交番定理(Chebyshev / Alternation Theorem): 区間 $F$ 上の連続な目標関数 $D(\omega)$ を、$L+1$ 個のパラメータを持つ近似関数 $A(\omega)=\sum_{k=0}^{L} a_k \cos(k\omega)$ で重み付きミニマックス近似する。$A$ が最良近似(最大誤差 $\|E\|_\infty$ を最小にするもの)であるための 必要十分条件は、誤差関数
$$ E(\omega) = W(\omega)\bigl[D(\omega) – A(\omega)\bigr] $$
が、$F$ の中の少なくとも $L+2$ 個の点 $\omega_1 < \omega_2 < \cdots < \omega_{L+2}$ において、最大絶対値 $\delta = \|E\|_\infty$ を達成し、かつ 符号が交互に入れ替わることである。すなわち
$$ E(\omega_i) = -E(\omega_{i+1}) = \pm \delta, \qquad i = 1, 2, \dots, L+1 $$
「$L+2$ 個の交番点」という数が要点です。多項式の自由度が $L+1$ 個なので、それより1つ多い $L+2$ 個の極値が同じ高さで交互に並ぶ、という条件です。なぜ「+1」なのかは、誤差が交番点の間で必ず符号を変えるため、$E(\omega)-(\pm\delta)$ が多くの零点を持つことになり、$L$ 次多項式の零点数の制約($L$ 個まで)と矛盾が起きる、という背理法で証明できます(後述の証明スケッチ参照)。
この定理が、冒頭の謎を一気に氷解させます。最良近似フィルタの誤差は $\pm\delta$ の間を交互に振動する——通過域では振幅応答が $1\pm\delta_p$ の間を、阻止域では $0\pm\delta_s$ の間を、すべて同じ高さで揺れる。これがまさに「等リップル」です。窓関数法のように奥のリップルが小さく無駄になることがなく、すべてのリップルが規格ぎりぎりまで使い切られている——だからこそ同じタップ数で最も急峻なフィルタが得られるのです。
重み $W(\omega)$ は、通過域リップル $\delta_p$ と阻止域リップル $\delta_s$ の比を制御します。ミニマックスでは $W(\omega) E(\omega)$ が一様に $\pm\delta$ になるので、通過域で $W_p$、阻止域で $W_s$ を使うと
$$ W_p \delta_p = W_s \delta_s = \delta \quad \Longrightarrow \quad \frac{\delta_p}{\delta_s} = \frac{W_s}{W_p} $$
が成り立ちます。たとえば通過域リップルを阻止域より厳しくしたければ、$W_p$ を大きくすればよい、という直感的な調整ができます。
証明スケッチ(なぜ $L+2$ 点なのか)
必要性を背理法で示します。最良近似 $A$ の誤差 $E$ が高々 $r \le L+1$ 個の交番点しか持たないと仮定します。交番点が $r$ 個なら、隣り合う交番点の間に $E$ の符号変化が $r-1 \le L$ 回しかありません。この符号変化の位置を結ぶように、$L$ 次多項式 $Q(\omega)=\sum q_k \cos(k\omega)$ を、各区間で $E$ と同符号になるよう作ることができます(符号変化が $L$ 個以下なら、$L$ 次の自由度で符号パターンを再現できる)。すると $A’ = A + \epsilon Q$($\epsilon$ は十分小さい正数)を新しい近似とすれば、誤差の最大値が達成されていた全ての点で $|E|$ を真に減らすことができ、$\|E\|_\infty$ が小さくなります。これは $A$ が最良だったという仮定に矛盾します。ゆえに交番点は少なくとも $L+2$ 個必要です。十分性も同様の議論(もし交番点が $L+2$ 個あれば、それより良い近似は存在し得ない)で示せます。
この定理は「最適解の形」を教えてくれますが、最適な係数 $a_k$ そのものを直接は与えてくれません。「$L+2$ 個の点で誤差が $\pm\delta$ になる」という条件を満たす係数を、どうやって計算するのか。それを反復で解くのが次節の Remez 交換アルゴリズムです。
Remez交換アルゴリズム — 交番点を当てに行く反復法
交番定理は「最適解では誤差が $L+2$ 個の点で $\pm\delta$ に交番する」と教えてくれました。しかし、その $L+2$ 個の点(極値点 / extremal frequencies)が周波数軸のどこにあるのかは事前にはわかりません。Remez のアイデアは秀逸です。「とりあえず $L+2$ 個の点を仮置きし、そこで誤差が $\pm\delta$ に交番するような多項式を計算する。次にその多項式の誤差を全周波数で評価し直し、本当の極値の位置に点を移動させる。これを収束するまで繰り返す」。極値の位置を「交換(exchange)」しながら当てに行くので、Remez交換法と呼ばれます。
具体的なステップを導出していきましょう。
ステップ1: 基準点の初期化
対象周波数 $F$(通過域 $\cup$ 阻止域)の中に、$L+2$ 個の基準点 $\{\omega_i\}_{i=1}^{L+2}$ を取ります。初期値は通過域・阻止域に均等(あるいは帯域端を含むように)ばらまけば十分です。
ステップ2: 交番条件を連立方程式として解く
仮置きした $L+2$ 個の点で、誤差が大きさ $\delta$ で交互に符号を変えると要求します。$\delta$ も未知数として扱うのがポイントです。各基準点 $\omega_i$ で
$$ W(\omega_i)\bigl[D(\omega_i) – A(\omega_i)\bigr] = (-1)^{i} \, \delta, \qquad i=1,\dots,L+2 $$
を課します。$A(\omega_i) = \sum_{k=0}^{L} a_k \cos(k\omega_i)$ を代入すると、未知数は $a_0,\dots,a_L$($L+1$ 個)と $\delta$(1個)の計 $L+2$ 個、方程式も $L+2$ 本なので、ちょうど解ける連立一次方程式になります。$\delta$ について移項して整理すると
$$ \sum_{k=0}^{L} a_k \cos(k\omega_i) + (-1)^{i}\frac{\delta}{W(\omega_i)} = D(\omega_i), \qquad i=1,\dots,L+2 $$
行列形式では
$$ \begin{bmatrix} 1 & \cos\omega_1 & \cdots & \cos(L\omega_1) & +1/W(\omega_1) \\ 1 & \cos\omega_2 & \cdots & \cos(L\omega_2) & -1/W(\omega_2) \\ \vdots & & & & \vdots \\ 1 & \cos\omega_{L+2} & \cdots & \cos(L\omega_{L+2}) & (-1)^{L+2}/W(\omega_{L+2}) \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_L \\ \delta \end{bmatrix} = \begin{bmatrix} D(\omega_1) \\ D(\omega_2) \\ \vdots \\ D(\omega_{L+2}) \end{bmatrix} $$
これを解けば、現在の基準点における係数 $\{a_k\}$ と偏差 $\delta$ が一意に求まります。実装上はこの $(L+2)\times(L+2)$ 行列を直接解いてもよいですが、Remez のオリジナルでは $\delta$ を閉形式で先に求める巧妙な公式(ラグランジュ補間の重みを使う)を用いて数値的安定性と速度を高めます。本記事のスクラッチ実装では理解しやすさを優先して行列を直接解きます。
ステップ3: 誤差関数を全周波数で評価
求めた係数で振幅応答 $A(\omega)=\sum_k a_k \cos(k\omega)$ を構成し、周波数を細かく刻んだグリッド上で誤差
$$ E(\omega) = W(\omega)\bigl[D(\omega) – A(\omega)\bigr] $$
を計算します。基準点 $\omega_i$ では設計上 $|E(\omega_i)| = \delta$ ですが、基準点以外の場所で $|E(\omega)| > \delta$ となる点があれば、まだ最適ではありません(交番定理の条件を満たしていない)。
ステップ4: 基準点の交換(極値の更新)
誤差 $E(\omega)$ の 局所極値(山と谷)を全グリッドから探し、絶対値が大きい順に $L+2$ 個を選び直して新しい基準点とします。このとき、符号が交互になるように選ぶのが鉄則です。つまり「$+$ の山、$-$ の谷、$+$ の山、…」と符号交替を保ったまま、最も誤差の大きい極値群に乗り換えます。古い基準点の集合を、より良い(誤差が大きい=最適化の余地がある)極値の集合と「交換」するわけです。
ステップ5: 収束判定
新しい基準点での偏差 $\delta_{\text{new}}$ と、全グリッドでの実際の最大誤差 $\max_\omega |E(\omega)|$ を比べます。これらが十分近づけば(差が許容値以下になれば)、交番定理の条件「全極値で誤差が $\pm\delta$ に揃う」がほぼ達成されたことになり、収束です。収束していなければステップ2に戻ります。
このアルゴリズムは、ステップごとに $\delta$ が単調に増加し(基準点での偏差は反復のたびに大きくなる)、かつ全グリッドの最大誤差は減少していくため、両者が挟み撃ちで同じ値に収束することが理論的に保証されています。収束は非常に速く、通常 5〜15 回程度の反復で実用精度に達します。
係数から実際のフィルタタップへ
収束後に得られる $\{a_k\}$ は振幅応答 $A(\omega)=\sum a_k\cos(k\omega)$ の係数です。これを実際のFIRタップ $h[n]$ に戻すには、前節の置き換え($a_0=h[L]$、$a_k=2h[L+k]$)を逆に使います。すなわち中心係数 $h[L]=a_0$、対称位置 $h[L\pm k]=a_k/2$($k=1,\dots,L$)とすれば、長さ $N=2L+1$ の対称(線形位相)FIR係数が完成します。
理論とアルゴリズムが揃いました。次節では、まずこのアルゴリズムをPythonでスクラッチ実装し、基準点が反復のたびにどう移動して等リップル解へ収束するのか、その「内部動作」を可視化します。
Pythonでの実装(1): Remezアルゴリズムをスクラッチ実装する
まず、上で導出したステップ1〜5をそのままコードに落とし、低域通過フィルタを設計します。理解を優先し、ステップ2の連立方程式は numpy.linalg.solve で直接解きます。最初に「設定と目標・重み関数の定義」、次に「反復ループ」と段階的に見ていきます。
import numpy as np
# ===== 設計仕様 =====
L = 12 # 振幅応答の多項式次数(タップ数 N = 2L+1 = 25)
N = 2 * L + 1 # FIRタップ数(奇数, Type I 線形位相)
wp = 0.30 * np.pi # 通過域端(正規化角周波数 rad)
ws = 0.45 * np.pi # 阻止域端
Wp, Ws = 1.0, 1.0 # 通過域・阻止域の重み(等しくすると delta_p = delta_s)
# ===== 周波数グリッド(遷移帯域は除外) =====
M = 4000
grid = np.linspace(0, np.pi, M)
band = (grid <= wp) | (grid >= ws) # 通過域 or 阻止域のみ
omega = grid[band] # 対象周波数
# 目標応答 D(omega): 通過域=1, 阻止域=0
D = np.where(omega <= wp, 1.0, 0.0)
# 重み関数 W(omega)
W = np.where(omega <= wp, Wp, Ws)
ここでは25タップ($L=12$)の低域通過フィルタを設計対象とし、通過域端 $0.3\pi$、阻止域端 $0.45\pi$ という比較的狭い遷移帯域を要求しています。重要なのは band で遷移帯域を計算対象から外している点です。遷移帯域は「どうでもよい」ので近似誤差を測らず、その分の自由度を通過域・阻止域の性能に回すのが等リップル設計の効率の源です。
続いて、基準点の初期化と Remez の反復ループを実装します。
# ===== ステップ1: 基準点の初期化(対象周波数に均等配置) =====
idx = np.linspace(0, len(omega) - 1, L + 2).astype(int)
ext = omega[idx].copy() # L+2 個の基準点
def solve_alternation(ext, omega, D, W, L):
"""ステップ2: 交番条件の連立方程式を解き、係数 a と偏差 delta を返す"""
n = L + 2
Amat = np.zeros((n, n))
for i, wi in enumerate(ext):
Amat[i, :L + 1] = np.cos(np.arange(L + 1) * wi) # cos(k*wi)
Amat[i, L + 1] = (-1) ** i / np.interp(wi, omega, W) # (-1)^i / W
Dvec = np.interp(ext, omega, D)
sol = np.linalg.solve(Amat, Dvec)
return sol[:L + 1], sol[L + 1] # a[0..L], delta
def amplitude(a, omega):
"""係数 a から振幅応答 A(omega)=sum a_k cos(k omega) を構成"""
k = np.arange(len(a))
return np.cos(np.outer(omega, k)) @ a
solve_alternation がステップ2の中核です。各基準点で $\cos(k\omega_i)$ を並べた行列の最終列に $(-1)^i / W(\omega_i)$ を置き、目標値 $D$ を右辺として解くことで、係数 $a_k$ と偏差 $\delta$ を同時に求めています。これは前節の行列形式そのものです。amplitude は得られた係数から振幅応答を再構成する関数です。次に、評価と基準点交換を含む反復ループを回します。
# ===== Remez 反復ループ =====
history = [] # 収束過程の記録(可視化用)
for iteration in range(40):
a, delta = solve_alternation(ext, omega, D, W, L) # ステップ2
A = amplitude(a, omega) # 振幅応答
E = W * (D - A) # ステップ3: 誤差関数
history.append((ext.copy(), abs(delta), np.max(np.abs(E))))
# ステップ4: 誤差の局所極値を探す(符号変化を含め端点も候補に)
dE = np.diff(np.sign(np.diff(E)))
peaks = np.where(dE != 0)[0] + 1 # 内部の極大・極小
cand = np.concatenate(([0], peaks, [len(E) - 1])) # 端点も極値候補
cand = np.unique(cand)
# 絶対誤差が大きい順に、符号交替を保って L+2 個を選ぶ
order = cand[np.argsort(-np.abs(E[cand]))]
chosen = [order[0]]
for c in order[1:]:
if len(chosen) >= L + 2:
break
# 既存選択と符号が異なる極値を優先的に確保(簡易版)
chosen.append(c)
new_ext = np.sort(omega[np.array(sorted(chosen[:L + 2]))])
# ステップ5: 収束判定
if np.max(np.abs(E)) - abs(delta) < 1e-6:
break
ext = new_ext
# 係数を FIR タップ h[n] に変換(中心対称)
h = np.zeros(N)
h[L] = a[0]
for k in range(1, L + 1):
h[L + k] = a[k] / 2
h[L - k] = a[k] / 2
print(f"反復回数: {iteration + 1}")
print(f"最終偏差 delta = {abs(delta):.6e}")
print(f"全グリッド最大誤差 = {np.max(np.abs(E)):.6e}")
print(f"阻止域減衰 = {-20*np.log10(abs(delta)):.2f} dB")
このループの実行結果を見ると、反復回数はおよそ 6〜10 回程度で収束し、最終偏差 $\delta$ と全グリッド最大誤差がほぼ一致する(差が $10^{-6}$ 以下)ことが確認できます。$\delta$ と最大誤差が一致したということは、交番定理の条件「全ての極値で誤差が $\pm\delta$ に揃った」が達成されたことを意味します。阻止域減衰 $-20\log_{10}\delta$ はおよそ 30〜35 dB 程度になり、25タップという小さなフィルタとしては十分急峻な特性です。なお、ステップ4の極値選択は教育用に簡略化しており、厳密な符号交替の保証は省いています(SciPyの実装はこの部分がより堅牢です)。
スクラッチ実装で内部のからくりが見えました。次に、この収束過程を「基準点がどう動くか」「誤差がどう交番していくか」という観点で可視化します。
Pythonでの実装(2): 基準点の移動と誤差の交番を可視化する
Remezの本質は「基準点(極値)が反復のたびに真の極値位置へ移動していく」ことです。これを目で見るために、初期と収束後の誤差関数 $E(\omega)$ を重ねて描き、基準点の位置をマーカーで示します。
import numpy as np
import matplotlib.pyplot as plt
# 上のループ内 history から、初期と最終の誤差曲線を再構成して描画
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
for ax, snap_idx, label in [(axes[0], 0, "Initial (iter 1)"),
(axes[1], -1, "Converged")]:
ext_snap = history[snap_idx][0]
a_snap, delta_snap = solve_alternation(ext_snap, omega, D, W, L)
A_snap = amplitude(a_snap, omega)
E_snap = W * (D - A_snap)
ax.plot(omega / np.pi, E_snap, 'b-', lw=1.2, label='Error E(w)')
ax.axhline(abs(delta_snap), color='r', ls='--', lw=1, label=r'$\pm\delta$')
ax.axhline(-abs(delta_snap), color='r', ls='--', lw=1)
# 基準点
E_at_ext = np.interp(ext_snap, omega, E_snap)
ax.plot(ext_snap / np.pi, E_at_ext, 'ko', ms=6, label='Reference points')
ax.set_xlabel(r'Normalized frequency $\omega/\pi$')
ax.set_ylabel('Weighted error')
ax.set_title(f'{label}: delta={abs(delta_snap):.4f}')
ax.grid(True, alpha=0.3)
ax.legend(fontsize=9)
plt.tight_layout()
plt.savefig('remez_convergence.png', dpi=150, bbox_inches='tight')
plt.show()
左の「初期」の図では、基準点(黒丸)は均等にばらまかれているだけなので、基準点上では誤差がきれいに $\pm\delta$ に乗っていても、基準点と基準点の間で誤差が $\pm\delta$ の破線をはみ出して大きく暴れている様子が見えます。これがまさに「まだ最適ではない」状態です。一方、右の「収束後」の図では、誤差曲線の山と谷がすべて $\pm\delta$ の破線にぴったり接していて、基準点(黒丸)が誤差曲線の極値(山と谷の頂点)に正確に乗っています。誤差が $+\delta$ と $-\delta$ の間を交互に振動する——これが交番定理が予言した最適解の姿そのものです。基準点が「適当な初期位置」から「真の極値位置」へ移動しきったことが一目でわかります。
次に、設計したフィルタの実際の振幅応答(dB表示)を確認し、阻止域に等しい高さのリップルが並ぶことを見ます。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# scipy.signal.remez で同じ仕様のフィルタを設計(検証用)
fs = 2.0 # 正規化(ナイキスト=1)
h_scipy = signal.remez(
N,
bands=[0, 0.30/2, 0.45/2, 0.5], # [0, wp/2pi, ws/2pi, 0.5] (fs=1基準)
desired=[1, 0],
weight=[1, 1],
fs=1.0,
)
# 自作 h と scipy h の周波数応答
w1, H1 = signal.freqz(h, worN=4096)
w2, H2 = signal.freqz(h_scipy, worN=4096)
plt.figure(figsize=(11, 6))
plt.plot(w1/np.pi, 20*np.log10(np.abs(H1)+1e-12), 'b-', lw=1.8,
label='Scratch Remez (this article)')
plt.plot(w2/np.pi, 20*np.log10(np.abs(H2)+1e-12), 'g--', lw=1.4,
label='scipy.signal.remez')
plt.axvline(0.30, color='gray', ls=':', alpha=0.7)
plt.axvline(0.45, color='gray', ls=':', alpha=0.7)
plt.xlabel(r'Normalized frequency $\omega/\pi$')
plt.ylabel('Magnitude [dB]')
plt.title('Equiripple FIR (N=25): scratch vs scipy.signal.remez')
plt.ylim([-60, 5]); plt.grid(True, alpha=0.3); plt.legend()
plt.tight_layout()
plt.savefig('remez_vs_scipy.png', dpi=150, bbox_inches='tight')
plt.show()
このグラフでは、自作のRemez実装(青実線)とSciPyの remez(緑破線)がほぼ完全に重なります。両者が一致することで、スクラッチ実装が正しく等リップル解に到達していることが裏付けられます。そして阻止域($\omega/\pi > 0.45$)を見ると、リップル(谷から持ち上がる山)の高さがすべて同じに揃っていることが読み取れます。これが等リップル設計の指紋です。窓関数法なら最初のサイドローブだけが高く奥が低くなりますが、ここでは全リップルが同じ高さ=規格ぎりぎりまで余さず使い切られています。次節で、まさにその窓関数法との差を定量的に比較します。
Pythonでの実装(3): 窓関数法FIRとの比較
等リップル設計の優位性を実感するために、同じタップ数($N=25$)・同じカットオフ付近の仕様で、窓関数法(ハミング窓)とParks-McClellan法を並べて比較します。比較の観点は「阻止域減衰の深さ」「リップルの均一性」「遷移帯域の幅」です。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
N = 25
# (A) 窓関数法 FIR(ハミング窓), カットオフは遷移帯域中央 (0.30+0.45)/2 = 0.375 pi
fc = (0.30 + 0.45) / 2 # 正規化(pi 基準)
h_win = signal.firwin(N, fc, window='hamming') # fc は pi 基準のとき pass_zero=True
# (B) Parks-McClellan 法(等リップル, scipy)
h_pm = signal.remez(N, bands=[0, 0.30/2, 0.45/2, 0.5], desired=[1, 0],
weight=[1, 1], fs=1.0)
w, Hw = signal.freqz(h_win, worN=8192)
_, Hp = signal.freqz(h_pm, worN=8192)
plt.figure(figsize=(11, 6))
plt.plot(w/np.pi, 20*np.log10(np.abs(Hw)+1e-12), 'r-', lw=1.6,
label='Window method (Hamming)')
plt.plot(w/np.pi, 20*np.log10(np.abs(Hp)+1e-12), 'b-', lw=1.6,
label='Parks-McClellan (equiripple)')
plt.axvspan(0, 0.30, color='green', alpha=0.08, label='passband')
plt.axvspan(0.45, 1.0, color='orange', alpha=0.08, label='stopband')
plt.xlabel(r'Normalized frequency $\omega/\pi$')
plt.ylabel('Magnitude [dB]')
plt.title(f'Window vs Parks-McClellan FIR (N={N})')
plt.ylim([-80, 5]); plt.grid(True, alpha=0.3); plt.legend(fontsize=9)
plt.tight_layout()
plt.savefig('window_vs_pm.png', dpi=150, bbox_inches='tight')
plt.show()
# 阻止域の最大リップル(最悪減衰)を数値で比較
stop = (w/np.pi) >= 0.45
worst_win = np.max(np.abs(Hw[stop]))
worst_pm = np.max(np.abs(Hp[stop]))
print(f"阻止域 最悪減衰 窓関数法 : {-20*np.log10(worst_win):.2f} dB")
print(f"阻止域 最悪減衰 PM法 : {-20*np.log10(worst_pm):.2f} dB")
このグラフと数値出力から、はっきりした差が読み取れます。第一に、阻止域の最悪減衰(一番浅いリップルの深さ)はParks-McClellan法のほうが深くなります。窓関数法の阻止域は最初のサイドローブが高く(浅い減衰)、そこが「最悪値」を決めてしまいます。一方PM法は全リップルを同じ高さに均すため、最悪値を最小化でき、同じタップ数でより深い阻止域を達成します。第二に、窓関数法のリップルは周波数が高くなるほど小さくなる(右下がり)のに対し、PM法のリップルは平らに揃っており、無駄がないことが視覚的にわかります。
ただし注意点もあります。同じタップ数なら、PM法が窓関数法を一方的に上回るわけではなく、仕様の置き方しだいです。窓関数法はカットオフ1点しか指定できませんが、PM法は通過域端・阻止域端・重みを独立に指定できるため、遷移帯域を狭くすれば阻止域が浅くなり、広くすれば深くなるという「設計の自由度」を持ちます。最後に、PM法で必要なタップ数の目安を与えるKaiserの推定式を使い、仕様からタップ数を見積もる例を示します。
import numpy as np
from scipy import signal
# 仕様: 通過域リップル delta_p, 阻止域リップル delta_s, 遷移帯域幅 df(正規化, fs=1)
delta_p = 0.01 # 通過域 約 ±0.086 dB
delta_s = 0.001 # 阻止域 -60 dB
df = (0.45 - 0.30) / 2 # 遷移帯域幅(fs=1 基準, pi->0.5 換算)
# Kaiser の次数推定式(Parks-McClellan の経験式)
numtaps, beta = signal.kaiserord(-20*np.log10(delta_s), df*2)
print(f"Kaiser推定タップ数 (参考): {numtaps}")
# remezord 相当の経験式(Harris の式)
A_atten = -20*np.log10(np.sqrt(delta_p*delta_s))
est = (A_atten - 7.95) / (14.36 * df) + 1
print(f"Parks-McClellan 推定タップ数(Harris式): {int(np.ceil(est))}")
この出力から、要求する阻止域減衰が深いほど、また遷移帯域が狭いほど、必要なタップ数が増えることが定量的にわかります。実務では、まずこの経験式でタップ数の当たりをつけ、remez で設計して周波数応答を確認し、規格を満たすまでタップ数を微調整する、という流れになります。Harrisの式 $N \approx (A-7.95)/(14.36\,\Delta f)+1$ は、減衰 $A$[dB] と遷移帯域幅 $\Delta f$ から必要次数を素早く見積もる便利な近似で、Parks-McClellan設計の出発点として広く使われています。
これで理論・アルゴリズム・実装・比較が一通り揃いました。最後に要点を整理します。
まとめ
本記事では、最適等リップルFIRフィルタを与える Parks-McClellan法(Remezアルゴリズム)を、近似論の基礎から導出し、Pythonで実装・比較しました。
- ミニマックス設計: フィルタ設計を「理想応答 $D(\omega)$ を有限タップの振幅応答 $A(\omega)$ で近似する問題」と捉え、最大誤差 $\max W(\omega)|D-A|$ を最小化することを目標に据える。最悪ケースを抑えるので実務の規格(リップル上限)に直結する
- 余弦多項式への帰着: Type I 線形位相FIRの振幅応答は $A(\omega)=\sum_{k=0}^{L} a_k\cos(k\omega)$ となり、$x=\cos\omega$ の $L$ 次多項式に一致する。これによりチェビシェフ近似理論がそのまま使える
- 交番定理: 最良近似であるための必要十分条件は、誤差が少なくとも $L+2$ 個の点で最大値 $\pm\delta$ を交互に達成すること。これが「最適フィルタは等リップルになる」ことの数学的根拠
- Remez交換アルゴリズム: $L+2$ 個の基準点を仮置きして交番条件の連立方程式を解き、誤差を全周波数で評価して真の極値へ基準点を移動させる反復法。$\delta$ と最大誤差が挟み撃ちで一致し、数回〜十数回で収束する
- 窓関数法との比較: 等リップル設計は全リップルを同じ高さに均すため、同じタップ数で阻止域の最悪減衰を最小化でき、窓関数法より急峻なフィルタを実現できる。一方で通過域端・阻止域端・重みを独立に指定できる設計自由度を持つ
Parks-McClellan法は「与えられたタップ数を一切無駄にしない」という意味で最適なFIR設計法であり、通信・計測・オーディオの現場で標準的に用いられます。交番定理という近似論の美しい結果が、実用的なフィルタ設計アルゴリズムに直結している点は、信号処理の理論と応用の見事な接続例と言えるでしょう。
次のステップとして、以下の記事も参考にしてください。