crossover_filter

何かあれば GitHub のリポジトリに issue を作るか ryukau@gmail.com までお気軽にどうぞ。

インデックスに戻る

Update: 2025-01-07

Table of Contents

周波数帯域の分割に用いるクロスオーバーフィルタを実装します。

かなりややこしいので概要を設けました。ページ上部の Table of Contents も多少は役に立つかと思います。

以下の 3 つのフィルタを実装します。

素朴な FIR フィルタ
IIR の Linkwitz-Riley フィルタ
FIR の Linkwitz-Riley フィルタ

IIR の Linkwitz-Riley フィルタは以下の属性の組み合わせによって実装の詳細が 6 つに分かれます。

フィルタが偶数次あるいは奇数次。
2 バンド、 3 バンド以上かつ汎用的、 3 バンド以上かつ効率重視。

偶数次、奇数次によってフィルタ係数の計算が変わります。偶数次は 2 次セクションだけで構築できますが、奇数次では 1 次のフィルタが現れます。

2 バンドでは補正オールパスフィルタが不要ですが、 3 バンド以上では補正オールパスフィルタが必要となってきます。また補正オールパスフィルタを計算するタイミングで汎用性と計算効率のトレードオフがあります。

FIR の Linkwitz-Riley フィルタの節では小さい部品から順にクロスオーバーフィルタを組み立てています。つまり:

1 次の複素数フィルタの FIR 近似から 2 次の複素共役フィルタの FIR 近似が組み立てられる。
2 次の複素共役フィルタの FIR 近似から Butterworth フィルタの FIR 近似が組み立てられる。
Butterworth フィルタの FIR 近似から線形位相 FIR Linkwitz-Riley フィルタが組み立てられる。

という流れになっています。

フィルタ係数が左右対称な線形位相の FIR フィルタはクロスオーバーフィルタとして使えます。

素朴な線形位相 FIR フィルタの利点は多様な設計方法が確立していることと実装が簡単なことです。欠点はカットオフ周波数の変更時の計算量が多く、レイテンシが生じることです。カットオフ周波数を変更するとフィルタ係数をすべて書き換える必要があります。また線形位相フィルタは出力の遅れを補正するためにレイテンシが生じます。

窓関数法 (windowed sinc) で使うローパスの線形位相 FIR フィルタ係数は以下の式で計算できます。

$\begin{matrix} (1) & fir_lp (n) = \frac{\sin (2 π n f_{0} / f_{s})}{π n}, n \in [- ceil (N_{FIR} / 2), floor (N_{FIR} / 2)] . \end{matrix}$

$f_{s}$ : サンプリング周波数 Hz 。
$f_{0}$ : カットオフ周波数 Hz 。
$N_{FIR}$ : フィルタ係数の数。タップ数。

さらに窓関数をかけることで周波数特性を変更できます。以下はパラメータ $α$ で手軽に特性を変更できる Kaiser 窓の近似計算式です。

$\begin{matrix} (2) & w_{n} = \frac{\cosh [α \sqrt{1 - {(\frac{n - K}{N_{FIR}})}^{2}}]}{\cosh (α)} . \end{matrix}$

式 $2$ は Oboe の Resampler で参照されている dsp.stackexchange.com の回答に掲載されていました。 Windows の VC++ か、 Linux の GCC 環境なら、式 $2$ の近似を使わなくても C++ 17 で追加された <cmath> の cyl_bessel_k で直接 Kaiser 窓を実装することができます。しかし、この記事を書いた時点では Mac の Xcode (Apple Clang) では数学特殊関数が使えなかったので式 $2$ の近似を掲載しています。

FIR の畳み込みの計算については「レイテンシのない畳み込み」に詳細を掲載しています。ここでは省略します。

以下は C++ 20 による実装へのリンクです。

C++ による素朴な FIR の実装を読む (github.com)

フィルタ係数が短いときは上のリンク先のように std::inner_product や単純な for によって簡単に実装できます。フィルタ係数が長いときは「レイテンシのない畳み込み」で紹介しているような FFT を用いた畳み込みで効率よく計算しなければリアルタイムの締め切りに間に合いません。

Linkwitz-Riley フィルタは Butterworth を 2 つ直列につないだ IIR フィルタです。以下は Linkwitz-Riley フィルタのブロック線図です。

任意の次数の Butterworth フィルタについて、ローパスを LP 、ハイパスを HP とすると Linkwitz-Riley フィルタによるクロスオーバーは (LP -> LP) + (HP -> HP) = AP と表せます。ここで -> は直列接続、 AP はオールパスです。つまり、分割した (LP -> LP) と (HP -> HP) を足し合わせても振幅特性が平坦なまま、ということです。

Linkwitz-Riley フィルタの次数は内部で使われる Butterworth フィルタの 2 倍となります。以降では Linkwitz-Riley フィルタの次数を $M$ 、内部で使われる Butterworth フィルタの次数を $N = M / 2$ とします。 2 つの Butterworth を直列につなぐため $M$ は常に偶数です。

3 バンド以上の Linkwitz-Riley フィルタの実装は $M = 4 n$ のときと $M = 4 n + 2$ で場合分けされます。これは位相の補正に使うオールパスフィルタの実装が奇数次と偶数次で変わるためです。

分割した信号を投げっぱなしにせず、併合時に再び処理を行えるときは補正オールパスフィルタの数を減らして効率よく計算できます。

$N$ 個の 2 次セクション (second order sections, sos) を直列につなぐことで 2 つ直列につないだ Butterworth フィルタを実装できます。すべてのセクションのカットオフ周波数はクロスオーバー周波数をそのまま使います。 $Q$ 値は earlevel.com の “Cascading filters” という記事で紹介されていた以下の式で計算します。

$\begin{matrix} (3) & Q_{k} = \frac{1}{2 \sin ((k + 0.5) \frac{π}{N})} = \frac{1}{2 \cos (\frac{k π}{2 N})}, k \in [0, N) . \end{matrix}$

2 次セクションは biquad フィルタのローパスとハイパスが使えます。「Biquad フィルタの比較」で実装を紹介しています。

以下は C++ による実装へのリンクです。テンプレートパラメータの filterType を lowpass あるいは highpass とすると 2 つ直列につないだ Butterworth フィルタとなります。

C++ による 2 つ直列につないだ Butterworth フィルタの実装を読む (github.com)

$M = 4 n$ と $M = 4 n + 2$ の場合で実装が変わります。どちらの場合でも、すべてのカットオフ周波数にはクロスオーバー周波数を使います。

2 次セクションは biquad オールパスを使います。 TPT (topology-preserving transform) という離散化に基づいた biquad オールパスの実装を以下のリンク先に掲載しています。 Audio EQ Cookbook のレシピでも動きます。

C++ による TPT biquad フィルタの実装を読む (github.com)

2 次セクションの数は $N_{AP} = N / 2$ 個です。 $Q$ 値は 2 つ直列につないだ Butterworth フィルタと同じ式が使えますが、 $N$ が $N_{AP}$ に置き換わります。以下は $Q$ 値を計算する式の再掲です。

$\begin{matrix} (4) & Q_{k} = \frac{1}{2 \sin ((k + 0.5) \frac{π}{N_{AP}})} = \frac{1}{2 \cos (\frac{k π}{2 N_{AP}})}, k \in [0, N_{AP}) . \end{matrix}$

以下は C++ による実装へのリンクです。テンプレートパラメータの filterType を allpass とすると $M = 4 n$ のときの補正オールパスフィルタとなります。

C++ による $M = 4 n$ のときの補正オールパスフィルタの実装を読む (github.com)

1 つの 1 次オールパスと $floor (N / 2)$ 個の 2 次セクションの組み合わせとなります。 $M = 2$ のとき 2 次セクションは 0 個です。

1 次オールパスは以下の伝達関数で表されます。

$\begin{matrix} (5) & H_{AP 1} (z) = \frac{a + z^{- 1}}{1 + z^{- 1}} . \end{matrix}$

$a$ はカットオフ周波数 $f_{c}$ とサンプリング周波数 $f_{s}$ から以下の式で計算できます。厳密にはオールパスの場合はカットオフ周波数とは言わず、ブレーク周波数 (break frequency) と言うそうです。ブレーク周波数は位相シフトの量が最大の半分となる周波数のことです。

$\begin{matrix} (6) & a = \frac{t - 1}{t + 1}, t = \tan (π f_{c} / f_{s}) . \end{matrix}$

差分方程式にします。

y0 = a * (x0 - y1) + x1;

x0: 現在の入力値。
x1: 1 サンプル前の入力値。
y0: 現在の出力値。
y1: 1 サンプル前の出力値。

これで計算できる形になりました。

2 次セクションの $Q$ 値は以下の正規化された Butterworth 多項式 (normalzied Butterworth polynomial) で計算します。

$\begin{matrix} (7) & Q = - 2 \cos (\frac{2 k + N - 1}{2 N} π), k \in [1, N_{AP}] . \end{matrix}$

$k$ の範囲が 1 から始まり、終端は含めていることに注意してください。

以下は Python 3 による実装へのリンクです。

Python 3 による $M = 4 n + 2$ のときの補正オールパスフィルタの実装を読む (github.com)

3 バンド以上の分割と併合について 2 つの方法を紹介します。原則としては汎用的な実装を行い、クロスオーバーがボトルネックとなっていれば効率のいい実装を使うことが考えられます。

2 バンドでは補正オールパスフィルタが不要なので以下の処理は実装しなくていいです。

汎用的な帯域分割は補正オールパスの総数が分割するバンド数に応じた三角数 (triangular number) となりますが、単純な加算で併合を行える特長があります。つまりバンド数が多いと分割の計算が重たくなりますが、分割後の使い勝手はこちらのほうがいいです。

Block diagram of 4 band crossover using Linkwitz-Riley filter (LR). Input first goes to LR, then splits to high and low. High goes to Butterworth allpass, and low goes to another LR. Repeat this until desired number of bands are obtained.

処理を書き下します。

入力を Linkwitz-Riley に通過させて high と low に分割。
ここまでで得られた high をそれぞれ新たなオールパスに通過させる。
Low を Linkwitz-Riley に通過させて high と low に分割。
所望のバンド数が得られるまで 1 に戻って分割を繰り返す。
所望のバンド数が得られたら、ここまでで得られたすべての high について以下の係数を乗算する。
- Linkwitz-Riley の次数が $4 m + 0$ なら $+ 1$ を乗算。 (何もしない)
- Linkwitz-Riley の次数が $4 m + 2$ なら $- 1$ を乗算。

以下は Python 3 による実装へのリンクです。クロスオーバー周波数のリスト cutoffsHz の長さがバンド数となります。

Python 3 による汎用的な Linkwitz-Riley クロスオーバーの分割と併合の実装を読む (github.com)

分割後の信号を再取得できるときは以下のブロック線図のようにオールパスの数を (バンド数) - 1 に減らすことができます。例えば 1 つの音のプラグインの内部で分割と併合の処理を両方行うときに使えます。この方法は出力を投げっぱなしにできないため汎用性は落ちます。

Block diagram of efficient 4 band crossover.

処理を書き下します。まずは分割します。

Low に入力を代入する。
Low を Linkwitz-Riley に通過させて high と low に分割。
High を出力のキューに加える。
現在の出力のキューの長さに low の 1 を加えた数が所望のバンド数でなければ手順 1 に戻る。

キューに格納された各帯域に任意の処理を施したあと、以下の併合の処理を行います。

変数 y を用意して 0 に初期化。
キューの先頭を取り出し、新たなオールパスに通過させ、出力を y に加算。
キューの長さが 1 で無ければ手順 1 に戻る。
y に以下の係数を乗算。
- Linkwitz-Riley の次数が $4 m + 0$ なら $+ 1$ を乗算。 (何もしない)
- Linkwitz-Riley の次数が $4 m + 2$ なら $- 1$ を乗算。
キューの最後の要素を取り出し、 y に加算して出力とする。

Martin Vicanek さんによる “A New Reverse IIR Filtering Algorithm” で Linkwitz-Riley フィルタを線形位相の FIR フィルタとして効率よく近似する方法が紹介されています。

Vicanek の論文では 2 次の複素共役フィルタの FIR を効率よく近似して、時間反転フィルタに変換する方法が紹介されています。ここで時間反転 (time reversal) とはインパルス応答の前後を逆向きにすることです。元の前向きフィルタと時間反転フィルタを直列につなぎ合わせることで線形位相フィルタを構築することができます。

Linkwitz-Riley フィルタは、次数の同じ 2 つの Butterworth フィルタを直列につないだフィルタです。つまり 2 つの Butterworth フィルタの片方を前向きフィルタ、もう片方を逆向きフィルタとすることで線形位相 FIR にできます。また Butterworth フィルタは 2 次の複素共役フィルタの組み合わせとして実装できます。

まずは以降で部品として使う 1 次の複素数フィルタの FIR 近似を実装します。以下は “A New Reverse IIR Filtering Algorithm” の式 1 として記載されている、 1 次の複素数フィルタの伝達関数です。 $c$ は複素数です。

$\begin{aligned} (8) & H_{c} (z) = \frac{1}{1 - c z^{- 1}} & = 1 + c z^{- 1} + c^{2} z^{- 2} + c^{3} z^{- 3} \dots \\ (9) & = (1 + c z^{- 1}) (1 + c^{2} z^{- 2}) (1 + c^{4} z^{- 4}) (1 + c^{8} z^{- 8}) \dots \end{aligned}$

何をしているのかというと、 $\frac{1}{1 - c z^{- 1}}$ という IIR フィルタの伝達関数を式 $8$ でインパルス応答に変形して、式 $9$ で多項式の因数分解を行っています。式 $8$ の IIR (無限インパルス応答) を適当に打ち切れば FIR フィルタとして計算できます。それだけだと嬉しくありませんが、式 $9$ のように因数分解することで計算量を圧縮するというのがアイデアです。

時間反転を行うときは式 $9$ の $(1 + c^{n} z^{- n})$ を $(c^{n} + z^{- n})$ に置き換えます。

$\begin{matrix} (10) & reverse (H_{c} (z)) = (c + z^{- 1}) (c^{2} + z^{- 2}) (c^{4} + z^{- 4}) (c^{8} + z^{- 8}) \dots \end{matrix}$

以下は式 $9, 10$ で示された 1 次の複素数フィルタの FIR 近似のブロック線図です。

Block diagram of complex 1-pole filter. Time forward and time reversed variants are shown.

この節の式やブロック線図では $8 = 2^{3}$ 次の因数 (factor) までしか計算していませんが、より高い次数を計算することで精度が上がります。以降では $2^{N}$ 次まで計算することを $N - 1$ ステージの計算と呼ぶことにします。

2 次の複素共役フィルタの FIR 近似は 1 次の複素数フィルタの FIR 近似を使って効率よく実装できます。以下は 2 次の複素共役フィルタの伝達関数です。 $c$ は複素数、 $\bar{c}$ は複素数 $c$ の複素共役です。

$\begin{matrix} (11) & H_{c . c .} (z) = \frac{1}{(1 - c z^{- 1}) (1 - \bar{c} z^{- 1})} . \end{matrix}$

$\bar{c} = a - i b$ とすると、適当な入力サンプル $x$ について以下のように伝達関数を変形できます。

$\begin{matrix} (12) & H_{c} x = (1 - a z^{- 1} + i b z^{- 1}) H_{c . c .} x . \end{matrix}$

ここで 1 次の複素数フィルタの出力を $H_{c} x = u + i v$ として実部 $u$ と虚部 $v$ に分解すると以下の式が得られます。

$\begin{matrix} (13) & u = (1 - a z^{- 1}) H_{c . c .} x, v = b z^{- 1} H_{c . c .} x . \end{matrix}$

$v$ の式を $z^{- 1}$ について解き、 $u$ の式に代入すると以下の式が得られます。

$\begin{matrix} (14) & H_{c . c .} x = u + \frac{a}{b} v . \end{matrix}$

上の式は “A New Reverse IIR Filtering Algorithm” の式 11 で示されている 2 次の複素共役フィルタの FIR 近似の効率のいい計算方法です。以下はこの計算のブロック線図です。

Block diagram of complex conjugate 2-pole filter.

以下は C++ と Python 3 による実装へのリンクです。 ComplexIIR クラスの process1PoleForward と process1PoleReversed が 1 次の複素数フィルタの FIR 近似の計算、 process2PoleForward と process2PoleReversed が 2 次の複素共役フィルタの FIR 近似の計算です。

C++ の実装はステージ数をテンプレートパラメータ stage として指定するために、テンプレートの再帰を使っています。理由はよくわかりませんが、再帰を展開した実装よりも、テンプレートの再帰を使った実装のほうが cl.exe (Version 19.36.32532 for x64) では速くなりました。

2 次の複素共役フィルタの FIR 近似を用いて Butterworth フィルタの FIR 近似を行います。この節では「ステップ応答が S 字を描くフィルタ」に掲載した Bessel フィルタを 2 次セクションに分割する手法を流用しています。 Butterworth フィルタも Bessel フィルタと同様に全極フィルタなので同じ手法で 2 次セクションに分割できます。

バイリニア変換した離散系の $N$ 次の Butterworth フィルタは以下の 2 次セクションで表されます。

$\begin{aligned} (15) & H_{Butter . N} (z) & = G \prod_{k = 1}^{N / 2} \frac{1 + 2 z^{- 1} + z^{- 2}}{1 - 2 Re (p_{k}) z^{- 1} + | p_{k} |^{2} z^{- 2}}, \\ (16) & G & = \prod_{k = 1}^{N / 2} \frac{1 - 2 Re (p_{k}) + | p_{k} |^{2}}{4} . \end{aligned}$

$p_{k}$ は 2 次セクションの極の片方、 $G$ はフィルタ全体のゲインです。 Butterworth フィルタの 2 次セクションの極は複素共役となるので、 2 次の複素共役フィルタの伝達関数の分母と一致します。以下は 2 次の複素共役フィルタの伝達関数の再掲です。

$\begin{matrix} (17) & H_{c . c .} (z) = \frac{1}{(1 - c z^{- 1}) (1 - \bar{c} z^{- 1})} . \end{matrix}$

ここで Butterworth の 2 次セクションと 2 次の複素共役フィルタに違いがあります。 Butterworth の 2 次セクションは分子が $1 + 2 z^{- 1} + z^{- 2}$ ですが、 2 次の複素共役フィルタの分子は $1$ です。つまり $1 + 2 z^{- 1} + z^{- 2}$ の部分を追加で計算する必要があります。

Butterworth フィルタの複素共役を除いた極 $p_{k}$ は以下の式で計算できます。 $f_{c}$ はカットオフ周波数、 $f_{s}$ はサンプリング周波数です。

$p_{k} = - 2 π \frac{f_{c}}{f_{s}} \exp (i \frac{π (2 k - N + 1)}{2 N}), k \in [0, N / 2) .$

さらに計算の途中で内部の値が大きくなりすぎることを防ぐために $G$ を以下の式で各セクションのゲイン $g$ に分割します。

$g = G^{\frac{1}{N / 2}} .$

式 $16$ を変形して 2 次セクションごとに適切なゲインを設定することもできますが、保持する変数が増えます。

ここまでをまとめると以下のブロック線図となります。フィルタは線形なので $g$ の位置を Output の直前に動かしても動作します。

Block diagram of FIR Butterworth lowpass and highpass.

以下は C++ による実装へのリンクです。

C++ による Butterworth フィルタの FIR 近似の実装を読む (github.com)

線形位相なのでハイパスは遅延した入力からローパスを減算すれば得られます。つまりハイパスの計算は省略できます。分割と併合についても IIR の Linkwitz-Riley のような補正オールパスは不要です。

Block diagram of FIR Linkwitz-Riley crossover.

以下は C++ による実装へのリンクです。

C++ による FIR Linkwitz-Riley の実装を読む (github.com)

😀 は良、 👹 は悪を表しています。

	素朴な FIR	Linkwitz-Riley IIR	Linkwitz-Riley FIR
線形位相	😀	👹	😀
レイテンシ	👹	😀	👹
実装の容易さ	😀	2 バンド 😀 / 3 バンド以上 👹	👹
クロスオーバーの対称性	👹	😀	😀

低レイテンシが要件なら Linkwitz-Riley IIR が適しています。線形位相が要件かつ手軽に実装したいなら素朴な FIR 、実装に手間をかけてもいいなら Linkwitz-Riley FIR が適しています。

クロスオーバーの対称性はフィルタの振幅特性のスロープやストップバンドでの低減率などがローパスとハイパスで対称かということです。非対称を 👹 としています。

素朴な FIR 、 Linkwitz-Riley IIR 、 FIR Linkwitz-Riley の 3 つについて、 2 バンドの場合の簡単なベンチマークを取りました。

Name	Elapsed [ms]
LinkwitzRileyIIR2Band4n	1.4590
LinkwitzRileyFIR2Band4n	2.4562
WindowedFIR2Band	2.8668

以下はベンチマークに使ったコードへのリンクです。

filter_notes/crossover_filter/code/cpp/benchmark.cpp at master · ryukau/filter_notes · GitHub

以下はパラメータの設定です。

Linkwitz-Riley IIR
- 4 次
Linkwitz-Riley FIR
- 4 次
- 8 ステージ (2^8 サンプルの FIR と等価)
素朴な FIR
- 255 サンプル

結果を見ると IIR の Linkwitz-Riley が速く、 FIR の実装はあまり差がありません。ただし、素朴な FIR は畳み込みの計算について大いに改善の余地があります。また Linkwitz-Riley FIR のステージ数や、素朴な FIR のフィルタ係数の数によって結果が変わることが予想されます。

このベンチマークは実は素朴な FIR で十分ではないか、という疑問を解消するために行ったのでパラメータを細かく変えてトレンドを調べていません。クロスオーバー周波数をリアルタイムで変更する用途であれば、素朴な FIR のフィルタ係数の更新がかなりの負荷になるので、 Linkwitz-Riley FIR の利用価値はあるように見えます。クロスオーバー周波数が固定なら素朴な FIR でも十分な気がします。

2024/10/03
- 「数値計算の安定性を高める」という個所を「計算の途中で内部の値が大きくなりすぎることを防ぐ」というより具体的な表現に変更。
2024/05/17
- 「2 つ直列につないだ Butterworth フィルタ」と「補正オールパスフィルタ」の $Q_{k}$ の式を修正。
2024/05/06
- 線形位相 FIR Linkwitz-Riley の下の図の High と Low が入れ替わっていた誤りを修正。
2024/04/30
- 分子とするべきところで分母と書いていた誤りを修正。
2023/10/30
- 概要の追加。
- 「2 つ直列につないだ Butterworth フィルタ」の式の出典を文中に追加。
- 文章の整理。