以前の記事では，低域・高域を別々に処理するようなマルチバンド処理（オーディオクロスオーバー）を実現するフィルタバンクの例として、2バンドの場合の実装方法を紹介しました．
ただ，より実用的なプラグイン／ソフトウェアを作ることを考えると，3バンド以上に分割できることが望ましいです．
しかし，3バンド以上のフィルタバンクは2バンドの構成方法のみで作ろうとすると非常に計算効率・使い勝手が悪いです．
そこで，二分木構造の回路によるマルチバンドフィルタバンクの設計方法を紹介します．
元論文：Complementary N-Band IIR Filterbank Based on 2-Band Complementary Filters, Favrot, Alexis and Faller, Christof, IWAENC 2010 Proceedings.
ソース：N-band Linkwitz-Riley crossovers - Page 2 - DSP and Plug-in Development Forum - KVR Audio
- 非実用的な実装：2バンド補完フィルタの直列接続
- 2バンド補間フィルタ x 2種による3バンドフィルタバンク設計
- 4バンド (AVL木)
- 5バンド
- まとめ
- サンプルコード
- パラメタ，プロット用スクリプト
- 参考文献
非実用的な実装：2バンド補完フィルタの直列接続
マルチバンドフィルタバンクの最も単純な実現方法は，2バンドの補完フィルタ並列回路の直列接続です．
以前の記事にも書きましたが，2バンド補完フィルタ (H_LPn, H_HPn) において H_LPn + H_HPnはオールパスフィルタとなるため，補完フィルタの並列回路の出力の振幅は変わりません。
しかし，バンド数だけ直列接続しているため計算効率はとても悪いです．
そのため，全ての補完フィルタを直列接続しない方法を考えていきます．
2バンド補間フィルタ x 2種による3バンドフィルタバンク設計
ダメな例
まず，単純に2バンド補間フィルタのどちらかの出力に2バンド補間フィルタを追加することを考えてみます．
この木構造のような回路にしていけば，計算効率は良くなりそうです．
しかし，上記回路全体の伝達関数を計算すると，実用的な振幅が変化しない回路は実現が難しいことが分かります*1．
3次IIR補完フィルタ（バターワース）を使った回路全体の周波数応答を計算してみると，振幅特性が一定でないことが確認できます*2。
オールパスフィルタによる回路全体のオールパスフィルタ化
論文で紹介されているトリックです．
2つ目の補完フィルタを追加したパスとは逆のパスに，オールパスフィルタ H_AP1 を追加します．
この H_AP1 は H_AP1=H_LP1 + H_HP1 を満たす、すなわち補完フィルタの並列回路と同じ位相特性を持つオールパスフィルタです．
すなわち，追加した2バンド補完フィルタの並列接続の伝達関数と位相特性が等しいオールパスフィルタを他のパスに挿入することによって，回路全体をオールパスフィルタ化できます．
実際に回路全体の周波数応答を計算してみると，振幅特性が一定となっています．
4バンド (AVL木)
3バンドから、さらに分割を増やしてみます．
先ほどと同様に，新たに追加した2バンド補完フィルタに対応するオールパスフィルタ H_AP2 も追加します．
周波数応答は以下の通りです．
5バンド
さらに分割を増やす（木構造を深くする）場合も，基本的にはこれまでと同様の方針でオールパスフィルタを追加していきます．
注意としては，4バンドの時は H_AP2 を1つ追加するだけでしたが，この5バンドの回路の場合は H_AP3 を2つ挿入しています。
木が深くなると，分割を増やした時に挿入するオールパスフィルタの数が増えることに注意が必要です．
周波数応答は以下の通りです．
まとめ
複数の補完フィルタを使ったマルチバンド処理のPython実装（3, 4, 5バンドの場合）を紹介しました．
木構造となっているため，そもそも単純な直列接続よりも効率が良く，また，並列化が効きやすくなると考えられます．
ただ，どうしても分割帯域数が多くなったり帯域幅が狭くなる場合は，周波数変換ベースの処理とどちらがよいのか計算コストや音質の観点で比較・検討するべきだと思います．
サンプルコード
補完フィルタのフィルタ係数・伝達関数計算
def compute_complementary_filters(fc, sr, order=3, return_coef=False): if order % 2 == 0: sys.exit() b, a = {}, {} # フィルタ係数・伝達関数算出(バターワースフィルタ) b['low'], a['low'] = scipy.signal.butter(order, fc/(sr/2), btype='low') w, h_lp = scipy.signal.freqz(b['low'], a['low']) b['high'], a['high'] = scipy.signal.butter(order, fc/(sr/2), btype='high') w, h_hp = scipy.signal.freqz(b['high'], a['high']) # 低域／高域用ハーフフィルタの並列接続時の伝達関数オールパスフィルタ h_ap = h_lp + h_hp if return_coef: return w, h_lp, h_hp, h_ap, b, a else: return w, h_lp, h_hp, h_ap
パラメタ，プロット用スクリプト
sr = 44100 #fc = [sr//8, sr//16] # (3-band) #fc = [sr//8, sr//16, sr//4] # (4-band) fc = [sr//8, sr//16, sr//4, sr//6] # (5-band) n_fc = len(fc) h_lp, h_hp, h_ap = [None] * n_fc, [None] * n_fc, [None] * n_fc for k, tmp_fc in enumerate(fc): w, h_lp[k], h_hp[k], h_ap[k] = compute_complementary_filters(tmp_fc, sr) # 回路全体の伝達関数 #h = h_lp[0] * (h_lp[1] + h_hp[1]) + h_hp[0] #h = h_lp[0] * (h_lp[1] + h_hp[1]) + h_hp[0] * h_ap[1] h = h_lp[0] * h_ap[2] * (h_lp[1] + h_hp[1]) + h_hp[0] * h_ap[1] * (h_lp[2] + h_hp[2]) h = h_lp[0] * h_ap[2] * h_ap[3] * (h_lp[1] + h_hp[1])\ + h_hp[0] * h_ap[1] * h_lp[2] * (h_lp[3] + h_hp[3])\ + h_hp[0] * h_ap[1] * h_hp[2] * h_ap[3] h_all = h_lp + h_hp + [h] labels = [f"LowPass[{k}]" for k in range(n_fc) ]\ + [f"HighPass[{k}]" for k in range(n_fc) ]\ + ["All"] # plot bode_plot( w, h_all, labels=labels, sr=sr, fc=fc, )
