何かあれば GitHub のリポジトリに issue を作るか ryukau@gmail.com までお気軽にどうぞ。
Update: 2024-08-06
Martin Vicanek さんによる Matched Second Order Digital Filters と Matched One-Pole Digital Shelving Filters で紹介されていた整合フィルタ (matched filter) を実装します。整合フィルタとはアナログの周波数特性に一致する離散フィルタのことです。
Audio EQ
Cookbook
などで紹介されている二次フィルタはバイリニア変換を使っています。バイリニア変換は
Python 3 でフィルタ係数
ユーザから入力されるパラメータの一覧です。
cutoffRadian
: カットオフ周波数。単位はラジアン (rad)
。Q
: Q 値
(quality factor) 。レゾナンス。G
:
ピークゲイン。例えばゲインが 10 倍なら
G=10
。ピーキングとシェルビングフィルタのみで利用。各フィルタの実装の後に掲載しているプロットは
以下は C++ による実装です。 MatchedBiquad
がフィルタのクラスです。
“Matched Second Order Digital Filters” の式 13 から式 16
に掲載されている、 Orfanidis
によるピーキングフィルタ のフィルタ係数の計算式を実装します。
def orfanidisPeaking(cutoffRadian, Q, G):
0 = cutoffRadian
ω= 1
G0 = 1
G1 = ω0 / Q
Δω = np.sqrt(G)
GB
0 = np.tan(0.5 * ω0)
Ω
= G * G
G_2 = GB * GB
GB_2 = G0 * G0
G0_2 = G1 * G1
G1_2
= np.sqrt((G_2 - G1_2) / (G_2 - G0_2)) * Ω0 * Ω0
W_2 = (1 + np.sqrt((GB_2 - G0_2) / (GB_2 - G1_2)) * W_2) * np.tan(0.5 * Δω)
ΔΩ
= (ΔΩ * ΔΩ) * np.abs(GB_2 - G1_2) - 2 * W_2 * (
C abs(GB_2 - G0 * G1) - np.sqrt((GB_2 - G0_2) * (GB_2 - G1_2))
np.
)= 2 * W_2 * (np.abs(G_2 - G0 * G1) - np.sqrt((G_2 - G0_2) * (G_2 - G1_2)))
D
= np.sqrt((C + D) / np.abs(G_2 - GB_2))
A = np.sqrt((G_2 * C + GB_2 * D) / np.abs(G_2 - GB_2))
B
= 1 + W_2 + A
a0 = -2 * (1 - W_2) / a0
a1 = (1 + W_2 - A) / a0
a2 = (G1 + G0 * W_2 + B) / a0
b0 = -2 * (G1 - G0 * W_2) / a0
b1 = (G1 + G0 * W_2 - B) / a0
b2
return [[b0, b1, b2, 1, a1, a2]] # scipy.signal の sos 形式。
“Matched Second Order Digital Filters” の式 20 から式 23 に掲載されている、 Massberg による二次整合ローパスフィルタのフィルタ係数の計算式を実装します。
def massbergLowpass(cutoffRadian, Q):
0 = cutoffRadian
ω
= Q * Q
Q2 = np.pi * np.pi / (ω0 * ω0)
t1 = 1 - t1
t2 = t1 / Q2
t3 = 1 / (np.sqrt(t2 * t2 + t3 * t3))
g1
= 0
Ωs if Q > np.sqrt(0.5):
= 2 * Q2 / np.sqrt(4 * Q2 - 1)
gr = ω0 * np.sqrt(1 - 1 / (2 * Q2))
ωr = np.tan(ωr / 2) * np.power((gr * gr - g1 * g1) / (gr * gr - 1), 1 / 4)
Ωs else:
= ω0 * np.sqrt(
ωm 1 - 1 / 2 / Q2 + np.sqrt((1 - 4 * Q2) / (4 * Q2 * Q2) + 1 / g1)
)= min(0.5 * ω0 * np.power(1 - g1 * g1, 1 / 4), np.tan(ωm / 2))
Ωs
= 2 * np.arctan(Ωs / np.sqrt(g1))
ωz = ωz * ωz / (ω0 * ω0)
z_tmp1 = 1 - z_tmp1
z_tmp2 = 1 / (z_tmp2 * z_tmp2 + z_tmp1 / Q2)
gz
= 2 * np.arctan(Ωs)
ωp = ωp * ωp / (ω0 * ω0)
p_tmp1 = 1 - p_tmp1
p_tmp2 = 1 / (p_tmp2 * p_tmp2 + p_tmp1 / Q2)
gp
= gz * gz
gz_2 = gp * gp
gp_2 = g1 - 1
β = np.sqrt(g1 * g1 * (gp_2 - gz_2) / (gz_2 * (g1 + gp_2) * β * β))
Qz = np.sqrt(g1 * (gp_2 - gz_2) / ((g1 + gz_2) * β * β))
Qp
= Ωs * Ωs
Ωs_2 = np.sqrt(g1)
sqrt_g1 0 = Ωs_2 + sqrt_g1 * Ωs / Qz + g1
β1 = 2 * (Ωs_2 - g1)
β2 = Ωs_2 - sqrt_g1 * Ωs / Qz + g1
β= Ωs_2 + Ωs / Qp + 1
γ 1 = 2 * (Ωs_2 - 1)
α2 = Ωs_2 - Ωs / Qp + 1
α
= β0 / γ
b0 = β1 / γ
b1 = β2 / γ
b2 = α1 / γ
a1 = α2 / γ
a2
return [[b0, b1, b2, 1, a1, a2]] # scipy.signal の sos 形式。
Vicanek の整合フィルタは伝達関数の分母に関する計算が共通しています。以下は関連する “Matched Second Order Digital Filters” の式の一覧です。
以下は共通部の実装です。
def solveDenominator(ω0, Q):
= 0.5 / Q
q = -2 * np.exp(-q * ω0)
a1 if q <= 1:
*= np.cos(np.sqrt(1 - q * q) * ω0)
a1 else:
*= np.cosh(np.sqrt(q * q - 1) * ω0)
a1 = np.exp(-2 * q * ω0)
a2
= np.sin(ω0 / 2)
sn 0 = 1 - sn * sn
φ1 = sn * sn
φ2 = 4 * φ0 * φ1
φ
= (1 + a1 + a2) ** 2
A0 = (1 - a1 + a2) ** 2
A1 = -4 * a2
A2
return (a1, a2, φ0, φ1, φ2, A0, A1, A2)
“Matched Second Order Digital Filters” の式 31 から式 34
の実装です。式 30 に続く文で
def matchedLowpass(cutoffRadian, Q):
0 = cutoffRadian
ω
0, φ1, φ2, A0, A1, A2 = solveDenominator(ω0, Q)
a1, a2, φ
= 1 + a1 + a2
sqrt_B0 = A0
B0
= Q * Q * (A0 * φ0 + A1 * φ1 + A2 * φ2)
R1 = (R1 - B0 * φ0) / φ1
B1
= 0.5 * (sqrt_B0 + np.sqrt(B1))
b0 = sqrt_B0 - b0
b1
return [[b0, b1, 0, 1, a1, a2]] # scipy.signal の sos 形式。
“Matched Second Order Digital Filters” の式 36 の実装です。
def matchedHighpass(cutoffRadian, Q):
0 = cutoffRadian
ω
0, φ1, φ2, A0, A1, A2 = solveDenominator(ω0, Q)
a1, a2, φ
= Q * np.sqrt(A0 * φ0 + A1 * φ1 + A2 * φ2) / (4 * φ1)
b0 = -2 * b0
b1 = b0
b2
return [[b0, b1, b2, 1, a1, a2]] # scipy.signal の sos 形式。
“Matched Second Order Digital Filters” の式 37 から式 41 の実装です。
def matchedBandpass(cutoffRadian, Q):
0 = cutoffRadian
ω
0, φ1, φ2, A0, A1, A2 = solveDenominator(ω0, Q)
a1, a2, φ
= A0 * φ0 + A1 * φ1 + A2 * φ2
R1 = -A0 + A1 + 4 * (φ0 - φ1) * A2
R2
= (R1 - R2 * φ1) / (4 * φ1 * φ1)
B2 = R2 - 4 * (φ0 - φ1) * B2
B1
= -0.5 * np.sqrt(B1)
b1 = 0.5 * (np.sqrt(B2 + b1 * b1) - b1)
b0 = -b0 - b1
b2
return [[b0, b1, b2, 1, a1, a2]] # scipy.signal の sos 形式。
“Matched Second Order Digital Filters” の式 44 から式 45
の実装です。
def matchedPeaking(cutoffRadian, Q, gain):
0 = cutoffRadian
ω= gain
G
0, φ1, φ2, A0, A1, A2 = solveDenominator(ω0, Q)
a1, a2, φ
= G * G * (A0 * φ0 + A1 * φ1 + A2 * φ2)
R1 = G * G * (-A0 + A1 + 4 * (φ0 - φ1) * A2)
R2
= A0
B0 = (R1 - R2 * φ1 - B0) / (4 * φ1 * φ1)
B2 = R2 + B0 - 4 * (φ0 - φ1) * B2
B1
= 1 + a1 + a2
sqrt_B0 = np.sqrt(B1)
sqrt_B1
= 0.5 * (sqrt_B0 + sqrt_B1)
W = 0.5 * (W + np.sqrt(W * W + B2))
b0 = 0.5 * (sqrt_B0 - sqrt_B1)
b1 = -B2 / (4 * b0)
b2
return [[b0, b1, b2, 1, a1, a2]] # scipy.signal の sos 形式。
Vicanek による、フィルタ係数の計算を簡略化した整合フィルタを実装します。このフィルタは計算が軽いだけでなく、より安定しているので、オーディオレートでフィルタのパラメータを変調するときに適しているそうです。アナログ特性との整合性はバイリニア変換よりはいいですが、前節の整合フィルタに比べるとやや劣ります。
solveDenominator()
を使いまわしています。
a1, a2
だけを使うので本格的な実装では計算を簡略化できます。
“Matched Second Order Digital Filters” の式 46 から式 47 の実装です。
def simpleMatchedLowpass(cutoffRadian, Q):
0 = cutoffRadian
ω0_2 = ω0 * ω0
ω
= solveDenominator(ω0, Q)
a1, a2, _, _, _, _, _, _
= 1 + a1 + a2
r0 = (1 - a1 + a2) * ω0_2 / np.sqrt((1 - ω0_2) ** 2 + ω0_2 / Q / Q)
r1
= 0.5 * (r0 + r1)
b0 = r0 - b0
b1
return [[b0, b1, 0, 1, a1, a2]]
“Matched Second Order Digital Filters” の式 48 から式 49 の実装です。
def simpleMatchedHighpass(cutoffRadian, Q):
0 = cutoffRadian
ω0_2 = ω0 * ω0
ω
= solveDenominator(ω0, Q)
a1, a2, _, _, _, _, _, _
= (1 - a1 + a2) / np.sqrt((1 - ω0_2) ** 2 + ω0_2 / Q / Q)
r1
= 0.25 * r1
b0 = -2 * b0
b1 = b0
b2
return [[b0, b1, b2, 1, a1, a2]]
“Matched Second Order Digital Filters” の式 50 から式 51 の実装です。
def simpleMatchedBandpass(cutoffRadian, Q):
0 = cutoffRadian
ω0_2 = ω0 * ω0
ω
= solveDenominator(ω0, Q)
a1, a2, _, _, _, _, _, _
= (1 + a1 + a2) / (ω0 * Q)
r0 = (1 - a1 + a2) * ω0 / Q / np.sqrt((1 - ω0_2) ** 2 + ω0_2 / Q / Q)
r1
= 0.5 * r0 + 0.25 * r1
b0 = -0.5 * r1
b1 = -b0 - b1
b2
return [[b0, b1, b2, 1, a1, a2]] # scipy.signal の sos 形式。
“Matched One-Pole Digital Shelving Filters” で紹介されていた整合一次ハイシェルビングフィルタのフィルタ係数の計算を実装します。対応する式は 4, 5, 10, 11, 12 です。
一次なのでレゾナンス
ローシェルビングフィルタはハイシェルビングの出力にゲインの逆数を乗算すれば作ることができます。ローシェルビング出力を
def matchedShelvingOnePole(cutoffRadian, gain):
= cutoffRadian / np.pi
fc = gain
G
= 0.9
fm = 1 - np.cos(np.pi * fm)
φm
def alphabeta(V):
return 2 / (np.pi * np.pi) * (1 / (fm * fm) + V) - 1 / φm
= alphabeta(1 / (G * fc * fc))
α = alphabeta(G / (fc * fc))
β
= -α / (1 + α + np.sqrt(1 + 2 * α))
a1 = -β / (1 + β + np.sqrt(1 + 2 * β))
b = (1 + a1) / (1 + b)
b0 = b * b0
b1
return [[b0, b1, 0, 1, a1, 0]] # scipy.signal の sos 形式。