2021/1/3: alphaの算出方法がsin(w0)/Qとなっていたため、sin(w0)/(2Q)に修正

導入

ラウドネスメーターで使われるLKFSやLUFSの算出における信号処理の流れは、以下のブロック図のように、

K周波数重み付け（K-weighting）
各チャネルでエネルギーを計算し集約
常用対数の領域でゲート処理

となっています。

f:id:Kurene:20210102121022p:plain — "Simplified block diagram of multichannel loudness algorithm", Rec. ITU-R BS.1770-4 3

最初の処理である K周波数重み付け（K-weighting）ですが、2つの2次IIRフィルタで信号をフィルタリングします。

それぞれ、下記のような振幅特性となっています。

f:id:Kurene:20210102121647p:plain — "Response of stage 1 of the pre-filter used to account for the acoustic effects of the head", Rec. ITU-R BS.1770-4 3

f:id:Kurene:20210102121811p:plain — "Second stage weighting curve", Rec. ITU-R BS.1770-4 3

課題

Rec. ITU-R BS.1770 に掲載されているK-weighting 用の2次IIRフィルタの係数は、サンプリング周波数が48kHzの場合のものです。

そのため、他のサンプリング周波数（例えば、44.1kHz ）に対応するには工夫が必要になります。

解決方法

1. 入力信号を 48 kHz にリサンプリング

入力信号をリサンプリングすることで、48 kHzのフィルタ係数を適用できるようにします。

既存のフィルタ係数をそのまま使え、リサンプリングの関数があればすぐ実装できる方法ですが、一般的にリサンプリングの計算コストは低くはないため、あまり良い方法ではありません。

2. 一般双2次フィルタのパラメタと係数の連立方程式を解く

Frequency-weighted filter - MATLAB - MathWorks 日本で紹介されている以下の論文の方法です。

AES E-Library » Implementation and Evaluation of Autonomous Multi-Track Fader Control

サンプリング周波数が48kHzの時のディジタルフィルタ係数から、アナログ領域での一般双2次フィルタのパラメタを推定する方法です。

この方法でやってみようとしましたが、Rec. ITU-R BS.1770 のフィルタ係数は a0=1.0 となっており、論文の方法では面倒な連立方程式を解く必要がでてきたのでやめました。

3. 双2次ハイパス／ハイシェルフフィルタのパラメタを推定

Rec. ITU-R BS.1770 を読むと、2つの双2次IIRフィルタはそれぞれ、

1st stage: ハイシェルフフィルタ
2nd stage: ハイパスフィルタ

で良さそうなことが分かります。

従って、それぞれに適当である双2次フィルタの形式に当てはめて、少数のパラメタを人力もしくは数理最適化で調整すれば、所望のフィルタ係数を得られそうです。

BiQuadフィルタの料理法 | g200kg Music & Softwares

以上を参考に、それぞれのフィルタで調整すべきパラメタをまとめると、

1st stage: ハイシェルフフィルタ
- カットオフ周波数
- ゲイン
- Q値（バンド幅）
2nd stage: ハイパスフィルタ
- カットオフ周波数
- Q値

となります。

このうち、カットオフ周波数とゲインに関しては、以下の情報を参考にしました。

How to calculate LUFS? - #4 by ImJimmi - JUCE

Q値については、分からなかったため、総当たりでRec. ITU-R BS.1770に掲載されているディジタルフィルタ係数との誤差が小さくなるQ値を探しました。ただし、Q値があまり大きくなりすぎないように、それっぽい値を選びました。

K-Weighting フィルタのパラメタ推定結果

最終的に、パラメタは以下のようにしました。

1st stage: ハイシェルフフィルタ
- カットオフ周波数 1500Hz
- ゲイン 4dB
- Q値（バンド幅）0.707
2nd stage: ハイパスフィルタ
- カットオフ周波数 37.5Hz
- Q値 0.86

以上の条件での双2次フィルタと、Rec. ITU-R BS.1770の周波数特性をプロットしました。

f:id:Kurene:20210102125201p:plain

目検ですが、ほぼ一致していると思います。また、フィルタ係数の値も、以下のようにほぼ一致していることを確認しました。

f:id:Kurene:20210103213248p:plain

実装

クリックで展開

def compute_k_weighting_highshelf(
        sr,
        fc=1500,
        q=1.414,
        gain_db=4.0
    ):

    omega = 2 * np.pi * fc / sr
    amp   = np.sqrt(10.0 ** (gain_db/20.0))
    alpha = np.sin(omega) / q * 0.5

    a, b = np.zeros(3), np.zeros(3)

    a[0] = (amp + 1.0) - (amp - 1.0) * np.cos(omega)\
         + 2.0 * np.sqrt(amp) * alpha
    a[1] = 2.0 * ((amp - 1.0) - (amp + 1.0) * np.cos(omega))
    a[2] = (amp + 1.0) - (amp - 1.0) * np.cos(omega)\
         - 2.0 * np.sqrt(amp) * alpha

    b[0] = amp * ((amp + 1.0) + (amp - 1.0) * np.cos(omega)\
         + 2.0 * np.sqrt(amp) * alpha)
    b[1] = -2.0 * amp * ((amp - 1.0) + (amp + 1.0) * np.cos(omega))
    b[2] = amp * ((amp + 1.0) + (amp - 1.0) * np.cos(omega)\
         - 2.0 * np.sqrt(amp) * alpha)

    b /= a[0]
    a /= a[0]

    return b, a



def compute_k_weighting_highpass(
        sr,
        fc=37.5,
        q=1.0
    ):

    omega = 2 * np.pi * fc / sr
    alpha = np.sin(omega) / q * 0.5

    a, b = np.zeros(3), np.zeros(3)

    a[0] = 1.0 + alpha
    a[1] = -2.0 * np.cos(omega)
    a[2] = 1.0 - alpha

    b[0] = (1.0 + np.cos(omega)) * 0.5
    b[1] = (1.0 + np.cos(omega)) * (-1.0)
    b[2] = (1.0 + np.cos(omega)) * 0.5

    b /= a[0]
    a /= a[0]

    return b, a