Wizard Notes

Pythonを使った音楽信号分析の技術録、作曲活動に関する雑記

マルチバンド(オーディオクロスオーバー)処理を実現するフィルタ(Butterworth, Linkwitz-Riley)のPython実装(2バンド)

低域・高域のような周波数帯域別に信号処理する場合、いい感じに帯域を分割するフィルタが必要になります。

ソフトウェア/プラグインの具体例としてはマルチバンドコンプレッサが挙げられます。

f:id:Kurene:20210923182752p:plain:w500

このような回路の特性としては、音質に影響を与えるような振幅/位相歪みが小さいと嬉しいです。

具体的には,低域/高域通過フィルタの並列接続(加算)した際の伝達関数が以下のような特性を持つように設計します。

  • 振幅特性: 特定の帯域を強調しない(フラットな特性を持つ)
  • 位相特性: 直線位相に近い、遅れが小さい

この記事ではオーディオクロスオーバーを実現するフィルタ(回路)としてバターワースフィルタリンクウィッツライリーフィルタを紹介します。

奇数次バターワースフィルタ

奇数次のバターワースフィルタであれば、フィルタのオーダーが奇数次であれば,低域/高域通過フィルタの並列接続の振幅特性はフラットになります。

1次バターワース

f:id:Kurene:20210923201322p:plain:w500

3次バターワース

f:id:Kurene:20210923201400p:plain:w500

5次バターワース

f:id:Kurene:20210923201427p:plain:w500

偶数時バターワースはなぜダメか?

偶数時のバターワースフィルタの場合はフラットな振幅特性になりません。

2n + 2 次の場合,低域/高域通過フィルタのカットオフ周波数での位相差は180度、つまり逆相となります.

従って,カットオフ周波数において減衰するノッチフィルタ(バンドストップフィルタ)として動きます。

f:id:Kurene:20210923230153p:plain:w500

一方で,2n 次の場合,低域/高域通過フィルタの位相差は0です.

従って,低域/高域通過フィルタの出力は単純に加算されますが,カットオフ周波数での利得は-3dBであるため,並列接続の出力は+3dBとなります.

つまり、カットオフ周波数で最大3dB増幅するピーキングフィルタとして動きます

f:id:Kurene:20210923230212p:plain:w500

Linkwitz-Riley フィルタ

Linkwitz-Riley フィルタはある次数のバターワースフィルタを直接接続することで生成できます.

Linkwitz-Riley フィルタの場合、低域/高域通過フィルタの位相差は0であり、それぞれのフィルタ出力の振幅特性が単純に加算されます。

2つのバターワース直接接続のカットオフ周波数での利得は-6dBであるため,回路全体の振幅特性はフラットになります。

f:id:Kurene:20210923232557p:plain:w500

2次 Linkwitz-Riley

f:id:Kurene:20210923201950p:plain:w500

4次 Linkwitz-Riley

f:id:Kurene:20210923202046p:plain:w500

6次 Linkwitz-Riley

f:id:Kurene:20210923224320p:plain:w500

8次 Linkwitz-Riley

f:id:Kurene:20210923224345p:plain:w500

まとめ

Pythonでscipyのバターワースフィルタを用いて,2つの周波数帯域に分割するオーディオクロスオーバー処理を実現するフィルタ処理(Butterworth, Linkwitz-Riley)を実装しました。

余談ですが、参考Webページに書かれている通り,これらのオーディオクロスオーバーフィルタバンク(回路全体)は,振幅特性は変わらず位相特性だけ変化しているのでオールパスフィルタと言えます。

参考文献

プロット用コード

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


# dB変換
to_dB = lambda h : 20 * np.log10(np.abs(h))

# 振幅特性プロット用
def plot_mag(w, h, label):
    if type(h) != list:
        sum_h = h
    else:
        sum_h = 0.0
        for tmp_h in h:
            sum_h = sum_h + tmp_h
    plt.plot(w, to_dB(sum_h), label=label, alpha=0.5)

# 位相特性プロット用
def plot_phase(w, h, label):
    if type(h) != list:
        sum_h = h
    else:
        sum_h = 0.0
        for tmp_h in h:
            sum_h = sum_h + tmp_h
    plt.plot(w, np.angle(sum_h) * 180 / np.pi, label=label, alpha=0.5)
    
# パラメタ
nfc = 0.5 # 正規化カットオフ周波数
order = int(sys.argv[1]) # フィルタ次数

# バターワースフィルタのフィルタ係数算出
b, a         = scipy.signal.butter(order, nfc, btype='low')
print(b, a)
w, h_bw_low  = scipy.signal.freqz(b, a)
b, a         = scipy.signal.butter(order, nfc, btype='high')
w, h_bw_high = scipy.signal.freqz(b, a)


# バターワースフィルタの応答プロット
plt.subplot(2,1,1)
plot_mag(w, h_bw_low, f"{order}-butter low")
plot_mag(w, h_bw_high, f"{order}-butter high")
plot_mag(w, [h_bw_low, h_bw_high], f"{order}-butter sum")

plt.ylabel('Magnitude [dB]')
plt.xlabel('Frequency [rad/sample]')
plt.ylim(-25, 5)
plt.grid()
plt.legend(loc='lower left')

plt.subplot(2,1,2)
plot_phase(w, h_bw_low, f"{order}-butter low")
plot_phase(w, h_bw_high, f"{order}-butter high")
plot_phase(w, [h_bw_low, h_bw_high], f"{order}-butter sum")

plt.ylabel('Phase (deg)')
plt.xlabel('Frequency [rad/sample]')
plt.grid()
plt.legend(loc='lower left')
plt.ylim(-180, 180)

plt.tight_layout()
plt.show()


# Linkwitz-Riley フィルタの応答算出
h_lr_low  = h_bw_low**2
h_lr_high = h_bw_high**2
# 2, 6, ... 次Linkwitz-Riley フィルタの場合、片方の極性を反転
h_lr_high = -h_lr_high if (order*2+2) % 4 == 0 else h_lr_high


# Linkwitz-Riley フィルタの応答プロット
plt.subplot(2,1,1)
plot_mag(w, h_lr_low, f"{order*2}-Linkwitz-Riley low")
plot_mag(w, h_lr_high, f"{order*2}-Linkwitz-Riley high")
plot_mag(w, [h_lr_low, h_lr_high], f"{order*2}-Linkwitz-Riley sum")

plt.ylabel('Magnitude [dB]')
plt.xlabel('Frequency [rad/sample]')
plt.ylim(-25, 5)
plt.grid()
plt.legend(loc='lower left')

plt.subplot(2,1,2)
plot_phase(w, h_lr_low, f"{order*2}-Linkwitz-Riley low")
plot_phase(w, h_lr_high, f"{order*2}-Linkwitz-Riley high")
plot_phase(w, [h_lr_low, h_lr_high], f"{order*2}-Linkwitz-Riley sum")

plt.ylabel('Phase (rad)')
plt.xlabel('Frequency [rad/sample]')
plt.ylim(-180, 180)
plt.grid()
plt.legend(loc='lower left')

plt.tight_layout()
plt.show()