音源制作の現場で最も利用されているプラグインの一つとして、ある周波数帯域を増幅／減衰させるイコライザがあります。
その一種であるグラフィックイコライザをPythonで実装してみます*1。
グラフィックイコライザについて
イコライザにおいて重要なパラメタは、
の3つになります。
フィルタの構成方法はいくつがありますが、特に双2次フィルタ（Biquad filter）を使う方法がよく知られています。
www.wizard-notes.com
多くのグラフィックイコライザでは、複数個のフィルタを用意します。そして、3つのパラメタの内の中心周波数と帯域幅を固定し、増幅／減衰させる量のみをユーザが操作できる仕様とします。
iTunesのグラフィックイコライザ
なお、グラフィックイコライザは1950年代から1960年代にかけてポストプロダクションや音声強調で流行し、例えばフェーダー付きであるプログラムイコライザ Langevin Model EQ-251A/EQ-252Aはグラフィックイコライザの先駆けと言われている製品だそうです。
（下記画像の参考元：EQ の概念とその向き合い方。）
引用: http://www.preservationsound.com/
一方で、パラメトリックイコライザの理論や製品は1970年代に出てきているようです（YAMAHA - DEQ7 Digital Equalizeなど）.
オクターブバンド分割（中心周波数と帯域幅の決定）
複数のフィルタを使うグラフィックイコライザでは、それぞれのフィルタで制御する周波数帯域を中心周波数と帯域幅（クロスオーバー周波数）で分けています。
グラフィックイコライザの中心周波数は、1/Nオクターブバンドで与えるのが一般的です。
すなわち、周波数が2倍になったときに対数軸で等間隔になるように周波数をＮ分割しています。
例えば、1000Hz を基準周波数とした1/3オクターブバンドを厳密に計算すると以下のようになります。
>>> import numpy as np
>>> k = np.arange(-12, 12)
>>> 1000 * 2 ** (k/3)
array([ 62.5 , 78.75, 99.21, 125. , 157.49, 198.43,
250. , 314.98, 396.85, 500. , 629.96, 793.7 ,
1000. , 1259.92, 1587.4 , 2000. , 2519.84, 3174.8 ,
4000. , 5039.68, 6349.6 , 8000. , 10079.37, 12699.21])
>>>
また、クロスオーバー周波数は幾何平均で算出することができます。
>>> (62.5*78.75)**0.5
70.1560760020114
ただし、これらの中心周波数・帯域幅は使いにくいため、実用的にはISO 266やJIS C 1514:2002（= IEC 61260:1995）で定められています。
"All about audio equalization: Solutions and frontiers", Table1
1/3オクターブバンドの中心周波数・帯域幅, "All about audio equalization: Solutions and frontiers", Table2
実装
Audio EQ Cookbook に基づいて複数の双二次フィルタを実装します。
注意点として、今回はQ値ではなく帯域幅でalphaの値を算出します。
このとき、帯域幅 bw は上限周波数fupと下限周波数flowから、
bw = np.log2(fup) - np.log2(flow)
として算出することに注意してください。
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import freqz
class GraphicEQ():
def __init__(self, sr, octave="1/3oct"):
self.sr = sr
self.octave = octave
self.set_params()
self.set_coefs()
def dB2amp(self, gain_dB):
return 10.0 ** (gain_dB/40.0)
def f2w(self, f):
return 2.0 * np.pi * f / self.sr
def set_params(self):
if self.octave == "1/1oct":
self.bandfreqs = np.array([
[22, 31.5, 44],
[44, 63, 88],
[88, 125, 177],
[177, 250, 355],
[355, 500, 710],
[710, 1000, 1420],
[1420, 2000, 2840],
[2840, 4000, 5680],
[5680, 8000, 11360],
[11360, 16000, 22720]
])
elif self.octave == "1/3oct":
self.bandfreqs = np.array([
[ 28.2, 31.5, 35.5], [ 35.5, 40, 44.7], [ 44.7, 50, 56.2],
[ 56.2, 63, 70.8], [ 70.8, 80, 89.1], [ 89.1, 100, 112],
[ 112, 125, 141], [ 141, 160, 178], [ 178, 200, 224],
[ 224, 250, 282], [ 282, 315, 355], [ 355, 400, 447],
[ 447, 500, 562], [ 562, 630, 708], [ 708, 800, 891],
[ 891, 1000, 1122], [1122, 1250, 1413], [ 1413, 1600, 1778],
[ 1778, 2000, 2239], [2239, 2500, 2818], [ 2818, 3150, 3548],
[ 3548, 4000, 4467], [4467, 5000, 5623], [ 5623, 6300, 7079],
[ 7079, 8000, 8913], [8913, 10000, 11220], [11220, 12500, 14130],
[14130, 16000, 17780],
])
self.n_band = self.bandfreqs.shape[0]
self.bandwidth = np.zeros(self.n_band)
self.gain_db = np.ones(self.n_band) * 3
self.amp = self.dB2amp(self.gain_db)
for k, (flow, fc, fup) in enumerate(self.bandfreqs):
self.bandwidth[k] = np.log2(fup) - np.log2(flow)
self.b = np.zeros((self.n_band, 3))
self.a = np.zeros((self.n_band, 3))
self.q = np.zeros(self.n_band)
def set_coefs(self):
amp = self.amp
for k in range(self.n_band):
fc = self.bandfreqs[k, 1]
w = self.f2w(fc)
cos_w, sin_w = np.cos(w), np.sin(w)
alpha = sin_w * np.sinh(0.5*np.log(2)*self.bandwidth[k]*w/sin_w)
self.q[k] = 0.5*sin_w/alpha
self.b[k, 0] = 1.0 + alpha * amp[k]
self.b[k, 1] = -2.0 * cos_w
self.b[k, 2] = 1.0 - alpha * amp[k]
self.a[k, 0] = 1.0 + alpha / amp[k]
self.a[k, 1] = -2.0 * cos_w
self.a[k, 2] = 1.0 - alpha / amp[k]
self.b[k] /= self.a[k,0]
self.a[k] /= self.a[k,0]
def freqz(self, worN=4096*10, plot_on=True):
self.h_list = [None for k in range(self.n_band)]
for k in range(self.n_band):
self.h_w, self.h_list[k] = freqz(self.b[k], a=self.a[k], worN=worN)
self.freqs = self.h_w / np.pi * (self.sr//2)
self.h_series = np.ones(worN, dtype=np.complex)
if plot_on:
plt.clf()
plt.subplot(2,1,1)
for k in range(self.n_band):
fc = self.bandfreqs[k, 1]
self.h_series *= self.h_list[k]
if plot_on:
plt.plot(self.freqs, 20*np.log10(np.abs(self.h_list[k])))
if plot_on:
plt.xscale("log")
plt.grid()
plt.xlim([10,self.sr//2])
plt.subplot(2,1,2)
plt.plot(self.freqs, 20*np.log10(np.abs(self.h_series)))
plt.xscale("log")
plt.grid()
plt.xlim([10,self.sr//2])
plt.tight_layout()
plt.show()
if __name__ == "__main__":
eq = GraphicEQ(44100)
eq.freqz()
周波数特性のプロット
1/1オクターブバンドと1/3オクターブバンドの周波数特性をプロットしました。
各フィルタのゲインは、すべて +3dBとしています。
プロットの上段が各フィルタの振幅特性、下段が全フィルタを直列接続した場合の振幅特性です。
周波数特性 (1/1オクターブバンド)
周波数特性 (1/3オクターブバンド)
まとめ、単純なグラフィックイコライザの問題点
Pythonでピーク/ノッチ型双二次フィルタによるグラフィックイコライザを実装しました。
上記のように単純な実装の場合、グラフィックイコライザの全体の周波数応答には問題があります。
具体的には、各フィルタのゲインは+3dBですが、最終的な応答を見ていると 3dB以上のゲインになっています。
これは、各フィルタが他のフィルタが制御する周波数帯域にも影響を与えるためです。
次回は、このような問題を解決する方法として、グラフィックイコライザのフィルタゲイン最適化手法を紹介したいと思います。
参考文献
