※ISO 532 (Zwicker method) による推定ラウドネスおよびラウドネスレベルの算出は以下の記事をご覧ください。

先日、下記の記事で統合ラウドネス値（Integrated loudness） LKFS/LUFS のPython実装 pyloudnorm をご紹介しました。

pyloudnormで平均ラウドネス値 (LUFS) 算出／ラウドネス正規化 - Wizard Notes

こちらのライブラリのコードは、Rec. ITU-R BS.1770-4 に沿って実装されているため、勧告の内容・数式と具体的な算出／実装方法の勉強に役立ちます。

勉強のため、pyloudnorm](https://github.com/csteinmetz1/pyloudnorm)での計算処理の流れを日本語コメントで解説したコードを作成しましたので、具体的な算出／実装方法が気になる方はご覧ください。

LUFSの算出のための処理を整理すると、以下の４つの処理に分けることができます。それぞれの処理の要点を以下にまとめました。

K-weighting filterの適用
- 頭部伝達関数や聴覚特性を考慮？
各チャネルのMean Squareの算出
- 処理ブロックの時間長は 400ms
- オーバーラップは75%, つまり100msごとに計算
ゲートなしラウドネスレベル算出
- 全チャネルのMean Squareを集約（重み付けて総和）した値
- チャネルごとにGで重み付けする
平均ラウドネス値 (Integrated loudness) 算出のためのゲート処理
- 絶対値ゲート
  - ゲートなしラウドネスレベルの値が-70未満の処理は除去
  - 無音や極端に小さいブロックを除去するため
- 相対値ゲート
  - 絶対値ゲートを通った信号の(3)を見て、(3)の平均値と比べて-10LU (dB) 低いブロックは除去

コード

クリックで展開

import numpy as np
import scipy.signal
from future.utils import iteritems

"""
Rec. ITU-R BS.1770-4 Integrated Loudness 
ゲート処理適応後のラウドネス値計算の解説 by Kurene
Original: https://github.com/csteinmetz1/pyloudnorm/blob/master/pyloudnorm/meter.py
License: MIT
"""

def integrated_loudness(x, meter):
    if x.ndim == 1:
        x = np.reshape(x, (x.shape[0], 1))
    # x.shape は (numSamples, numChannels)    
    numSamples, numChannels = x.shape[0], x.shape[1]
    
    rate       = meter.rate
    block_size = meter.block_size

    # K-weighting フィルタリング
    y = x.copy()
    for (filter_class, filter_stage) in iteritems(meter._filters):
        for ch in range(numChannels):
            #x[:,ch] = scipy.signal.lfilter(filter_stage.b, filter_stage.a, x[:,ch])
            y[:,ch] = filter_stage.apply_filter(y[:,ch])
    
    # LUFS算出のための変数
    G = [1.0, 1.0, 1.0, 1.41, 1.41] # 各チャンネルへの重み付け (L, C, R, Ls, Rs)
    T_g = block_size # 処理ブロックの時間長（秒数）
    Gamma_a = -70.0       # -70 LKFS = ゲート処理におけるスレッショルド
    overlap = 0.75        # ゲート処理の時間フレームのオーバーラップ
    step = 1.0 - overlap  # ブロックサイズの何%ホップするか
    
    T = numSamples / rate # 入力信号の秒数
    numBlocks = int(np.round(((T - T_g) / (T_g * step)))+1) # ゲート処理のブロック数
    j_range = np.arange(0, numBlocks) # ゲート処理ブロックのインデックスを格納した配列
    z = np.zeros(shape=(numChannels,numBlocks)) # Momentary Loudness

    # 各チャネルごとに、Mean Square の算出
    for i in range(numChannels): 
        for j in j_range: # 処理ブロック単位で処理
            l = int(T_g * (j * step    ) * rate)
            u = int(T_g * (j * step + 1) * rate)
            # 式 (1)
            z[i,j] = (1.0 / (T_g * rate)) * np.sum(np.square(y[l:u, i]))

    # ゲートなしラウドネスレベル の算出
    # 各チャネルのMean Squareを重み付け総和
    l_all = [0.0 for _ in j_range]
    for j in j_range:
        tmp = 0.0
        for i in range(numChannels):
            tmp += G[i] * z[i,j]
            
        with np.errstate(divide='ignore'):
            # 式 (4)
            l_all[j] = -0.691 + 10.0 * np.log10(tmp) 
        
    # ゲート処理
    # まず、絶対値ゲート処理
    #l_allが閾値gamma_a 以上の処理ブロックのインデックス配列 J_g を作る
    # 
    J_g = [j for j,l_j in enumerate(l_all) if l_j >= Gamma_a]
    # 各チャネルごとに平均を取る。従って、z_avg_gated は長さがチャネル数 (numChannels) の配列
    # 式(5)
    z_avg_gated = [ np.mean([z[i,j] for j in J_g]) for i in range(numChannels)]
    # 相対値ゲートの閾値計算
    # z_avg_gated を Gでチャンネル重み付けをして総和をとり、閾値 Gamma_r を計算
    # 式 (6)
    Gamma_r = -0.691 + 10.0 * np.log10(np.sum([G[i] * z_avg_gated[i] for i in range(numChannels)])) - 10.0

    # 絶対値ゲートおよび相対値ゲート処理
    # l_allが閾値gamma_a 以上の処理ブロックのインデックス配列 J_g を作る
    J_g = [j for j,l_j in enumerate(l_all) if (l_j > Gamma_r and l_j > Gamma_a)]
    # 各チャネルごとに平均を取る
    z_avg_gated = np.nan_to_num(np.array([np.mean([z[i,j] for j in J_g]) for i in range(numChannels)]))
    # ゲート処理されたラウドネスの計算
    # G でチャンネル重み付けをして総和
    # 式(7) 
    with np.errstate(divide='ignore'):
        LUFS = -0.691 + 10.0 * np.log10(np.sum([G[i] * z_avg_gated[i] for i in range(numChannels)]))

    return LUFS