※ITU-R BS, 1770によるラウドネスレベルの算出は以下の記事をご覧ください。

はじめに

ラウドネスの算出・実装方法の勉強のため、心理音響評価モジュール／開発フレームMoSQIToの非定常信号向けラウドネス（時変ラウドネス）算出モジュールのソースコードを分析／コメント翻訳してみました。

ソースコードの分析において、以下の文献を参考にしています。

MoSQIToや、MoSQIToを使ったラウドネス算出については、以下の記事をご参照ください。

内容に誤り等ありましたら、適宜修正したいと思います。

はじめに
MoSQIToを使ったラウドネス算出のソースコード
各信号処理ブロックの説明
まとめ（＋疑問）
補足：参考文献のラウドネス算出ブロック図／説明
- Zwicker Loudness Measurement with R&S® UPV Audio Analyzer Application Note
- MathWorks - acousticloudness

`MoSQITo`を使ったラウドネス算出のソースコード

ラウドネス算出処理の流れを追いやすくするため、MoSQIToのラウドネス算出チュートリアルで使われているコードをベースに、 mosqito.functions.loudness_zwicker.comp_loudnessおよびmosqito/functions/loudness_zwicker/loudness_zwicker_timeを展開したコードとなっています。

また、可読性向上のため一部の変数名を変えております。

github.com

クリックで展開

import numpy as np
import matplotlib.pyplot as plt

from mosqito.functions.shared.load import load
from mosqito.functions.loudness_zwicker.comp_loudness import comp_loudness
from mosqito.functions.oct3filter.comp_third_spectrum import comp_third_spec
from mosqito.functions.loudness_zwicker.loudness_zwicker_shared import calc_main_loudness
from mosqito.functions.loudness_zwicker.loudness_zwicker_nonlinear_decay import calc_nl_loudness
from mosqito.functions.loudness_zwicker.loudness_zwicker_shared import calc_slopes
from mosqito.functions.loudness_zwicker.loudness_zwicker_temporal_weighting import loudness_zwicker_temporal_weighting


dirpath  = "./mosqito_master/mosqito/validations/loudness_zwicker/data/ISO_532-1/Annex B.5/"
filepath ="Test signal 24 (woodpecker).wav"

# パラメタ
is_stationary = False # 非定常信号ならばFalse, 定常信号ならば True
calib = 2 * 2**0.5    # 振幅値補正
field_type = 'free'

# 分析対象である信号をロード（サンプリング周波数は48kHz）
signal, fs = load(is_stationary, dirpath+filepath, calib=calib)
n_samples = signal.shape[0]

# 1. 1/3オクターブバンド分析
# 25, 31.5, ..., 10000, 12500, 計28バンド
# third_octave_levels.shape == (28, n_sample_level)
# n_sample_level = n_samples // 24 + 1 
third_spec = comp_third_spec(is_stationary, signal, fs)
third_octave_levels = third_spec['values']  
n_sample_level = third_octave_levels.shape[1]

# 2. レベル補正、コアラウドネス算出
core_loudness = np.zeros((21, n_sample_level))
for k in np.arange(n_sample_level - 1):
    core_loudness[:, k] = calc_main_loudness(third_octave_levels[:, k], field_type)

# 3. 時間マスキング
core_loudness = calc_nl_loudness(core_loudness)

# 4. バーク尺度へのマッピング，周波数マスキング
n_bark = int(24 / 0.1)
spec_loudness = np.zeros((n_bark, n_sample_level)) # ラウドネス密度
loudness = np.zeros(n_sample_level) # ラウドネス
for k in np.arange(n_sample_level):
    loudness[k], spec_loudness[:, k] = calc_slopes(core_loudness[:, k])

plt.subplot(2,1,1)
plt.plot(np.mean(core_loudness, axis=1))
plt.xlabel("Critical band")
plt.title("Core loudness")
plt.subplot(2,1,2)
plt.plot(np.mean(spec_loudness, axis=1))
plt.xlabel("Bark scale")
plt.title("Specific loudness")

plt.tight_layout()
plt.show()

# 5. ラウドネスへの時間重み付け
filt_loudness = loudness_zwicker_temporal_weighting(loudness)

# デシメーションによって、ラウドネス、ラウドネス密度を 0.5ms から 2msごとのデータに変換
dec_factor = 4
N      = filt_loudness[::dec_factor]
N_spec = spec_loudness[:, ::dec_factor]

bark_axis = np.linspace(0.1, 24, n_bark)  # バーク尺度
time = np.linspace(0, 0.002*(N.size - 1), N.size) # 2msごとに算出されるため

# プロット
plt.plot(time, N)
plt.xlabel("Time [s]")
plt.ylabel("Loudness, [sone]")
plt.title(filepath)
plt.show()

各信号処理ブロックの説明

MoSQIToでのラウドネス算出の信号処理の流れは以下のようになっています。

各信号処理では、周波数重み付け、周波数／時間マスキングなどといった、人間の聴覚特性を考慮した処理であることにに注意してください。

f:id:Kurene:20210118084814p:plain

以下では、各処理の詳細を追っていきます。

1. 1/3オクターブバンド分析

mosqito.functions.oct3filter.comp_third_spectrum.comp_third_spec関数によって、1/3オクターブバンド分析が行われます。

非定常音の場合はmosqito.functions.oct3filter.calc_third_octave_levels.calc_third_octave_levels関数、定常音の場合はmosqito.functions.oct3filter.oct3spec.oct3spec関数が使われます。

calc_third_octave_levels関数では、28バンドごとに2次のバンドパスフィルタを用いて出力信号を算出しています。注意すべき点として、単なる1/3オクターブバンド分析ではなく、周波数重み付け、平滑化、SPLへの変換、デシメーションなどの細かい処理があります。

MoSQITo/calc_third_octave_levels.py at master · Eomys/MoSQITo · GitHub

# コアの部分のみ抽出
for i_bands in range(n_level_band):
    # 2次バンドパスフィルタ (ISO 532-1 section 6.3 and A.2参照)
    # filter_gainで周波数重み付け
    coeff = third_octave_filter_ref - third_octave_filter[i_bands,:,:]
    sig_filt = filter_gain[i_bands] * signal.sosfilt(coeff, sig)
    # 各周波数バンドの中心周波数
    center_freq = 10**((i_bands-16) / 10) * 1000
    # 信号のパワー算出および平滑化
    # 周波数バンドの中心周波数に応じたローパスフィルタ (ISO 532-1 section 6.3 参照)
    sig_filt = square_and_smooth(sig_filt, center_freq, 48000)
    # SPLへの変換およびデシメーション (sr=2kHz相当)
    third_octave_level[i_bands,:] = 10 * np.log10((sig_filt[::dec_factor]
            + tiny_value) / i_ref)

なお、周波数重み付けは低周波での耳の感度が低いことを考慮したもので、周波数重み付け用の係数配列filter_gainは以下のようになっています。
f:id:Kurene:20210117181459p:plain

2. 臨界帯域レベル（コアラウドネス）の算出

1/3オクターブバンド分析で算出した信号レベルを入力として、コアラウドネスを算出しています。

この処理は、デシメーションされた1サンプルごとに行われます。時間換算すると、24/48000 = 0.5msec ごとにコアラウドネスを算出していることになります。

なお、この処理を担うmosqito.functions.loudness_zwicker.loudness_zwicker_shared.calc_main_loudness関数は定常信号向けのラウドネス算出でも使われています。

MoSQITo/loudness_zwicker_shared.py at master · Eomys/MoSQITo · GitHub

def calc_main_loudness(spec_third, field_type):
    #...
    ti = np.zeros((dll.shape[1], 1))
    for i in np.arange(dll.shape[1]):
        j = 0
        while spec_third[i] > rap[j] - dll[j, i] and j < dll.shape[0]:
            j += 1
        xp = spec_third[i] + dll[j, i]
        ti[i] = np.power(10, (xp / 10))
    
    # Determination of levels LCB(1), LCB(2) and LCB(3) within the
    # first three critical bands
    gi = np.zeros(3)
    gi[0] = ti[0:6].sum()
    gi[1] = ti[6:9].sum()
    gi[2] = ti[9:11].sum()
    lcb = np.zeros(3)
    lcb[gi > 0] = 10 * np.log10(gi[gi > 0])
    
    # Calculation of main loudness
    # 低周波数帯域の補正
    nm = np.zeros(20)
    le = spec_third[8:]  #160Hz以上の1/3オクターブバンドレベル
    le = le.reshape((20))
    le[0:3] = lcb    　　# 160Hz, 200Hz, 250Hzを lcbに置き換え 
 
   # 耳の伝達特性 a0 で補正
    le = le-a0               

    # 拡散音場の場合、ddfで補正
    if field_type == "diffuse":
        le += ddf           
    
    # コアラウドネス算出
    # 閾値処理 ltq以上のleのみ
    s = 0.25
    i = le > ltq
    le[i] -= dcb[i] #臨界帯域補正
    mp1 = 0.0635 * np.power(10,0.025 * ltq[i])
    mp2 = np.power(1 - s + s * np.power(10,0.1 * (le[i] - ltq[i])),0.25 ) -1 
    nm[i] = mp1 * mp2
    nm[nm < 0] = 0
    nm = np.append(nm, 0)
    
    # 低域の補正
    korry = 0.4 + 0.32 * nm[0] ** 0.2
    if korry <= 1:
        nm[0] *= korry
    return nm

細かい計算はさておき、低周波数帯域のレベル補正、臨界帯域スケールに変換するための補正、絶対閾値処理などのようなレベル補正を行い、コアラウドネスを算出しているようです。

なお、配列の次元数は21となっているが、この時点では21番目の値は0.0が入っているため、他の文献と同様に出力は20帯域であると見なせます。

3. 時間マスキング (Temporal_masking)

時変ラウドネス計算では、経時マスキング (時間マスキング)の影響を加味しています。

聴神経の発火・減衰特性に基づく時間的な聴覚特性であり、詳細は以下のページを参考にしてください。

"小野測器 - 音質評価とは 5.6 時間マスキング"

具体的な処理は複雑なため割愛しますが、diode-capacitor-resistor network と呼ばれる回路をシミュレートしているようです。

MoSQITo/loudness_zwicker_nonlinear_decay.py at master · Eomys/MoSQITo · GitHub

4. バーク尺度へのマッピング，周波数マスキング

mosqito.functions.loudness_zwicker.loudness_zwicker_shared.calc_slopes関数を使い、コアラウドネスから部分ラウドネスおよびラウドネスを算出します。

この処理も、デシメーションされた1サンプルごとに行われます。また、定常信号向けラウドネス算出でも使われています。

MoSQITo/loudness_zwicker_shared.py at master · Eomys/MoSQITo · GitHub

細かい処理は割愛しますが、周波数マスキングの影響を計算した後、20帯域である臨界帯域をバーク尺度へ変換しています。よく心理音響の教科書に書かれている、ラウドネス計算用チャートへの書き込みに対応する処理だと思います。

ある音（"Test signal 24 (woodpecker).wav"）のコアラウドネスと部分ラウドネスの時間平均をそれぞれ上下にプロットしてみました。部分ラウドネスで、特に周波数が高くなる方向にスロープがあることが確認できます。

f:id:Kurene:20210118093616p:plain

5. ラウドネスへの時間重み付け

時変ラウドネスの算出では、部分ラウドネスの総和によりラウドネスを算出した後、ローパスフィルタリングが行われています。

理由としては、loudness_zwicker_temporal_weighting.py のコメントには,

in order to simulate the duration dependent behaviour of loudness perception for short impulses.

（短いインパルスに対するラウドネス知覚の持続時間依存性の挙動を再現するため）

と書かれています。

具体的な処理の流れは以下のようになっています。

時定数の異なる2つのローパスフィルタを並列して適用する
- LP_1: 3.5 msec
- LP_2: 70 msec
2つのローパスフィルタの出力信号を重み付けて加算する
- LP_1: 0.47
- LP_2: 0.53

入力信号はサンプリング周波数＝2kHz相当ですが、ローパスフィルタの計算ではオーバーサンプリング（48kHz）して計算しているようです。

MoSQITo/loudness_zwicker_temporal_weighting.py at master · Eomys/MoSQITo · GitHub

def loudness_zwicker_temporal_weighting(loudness):
    sample_rate = 2000
    tau = 3.5 * 10**-3
    filt_loudness_1 = loudness_zwicker_lowpass_intp(loudness, tau, sample_rate)
    tau = 70 * 10**-3
    filt_loudness_2 = loudness_zwicker_lowpass_intp(loudness, tau, sample_rate)

    loudness = 0.47 * filt_loudness_1 + 0.53 * filt_loudness_2
    
    return loudness

なお、ローパスフィルタリングはmosqito.functions.loudness_zwicker.loudness_zwicker_lowpass_intp.loudness_zwicker_lowpass_intp関数で処理しています。

まとめ（＋疑問）

1点だけ最後まで不明なのが、MoSQIToの時変ラウドネスの算出では、部分ラウドネスに対して時間重み付けが適用されていない点でした。機会があれば、ISO532などを読んで直接確認してみたいと思っています。

補足：参考文献のラウドネス算出ブロック図／説明

Zwicker Loudness Measurement with R&S® UPV Audio Analyzer Application Note

参考までに、参考文献[2]の ”2.2.1 Overview” に書かれているブロック図および各ブロックの説明を引用します。

注意：定常信号向けラウドネス算出のブロック図のため、時間マスキングおよびラウドネスへの時間重み付けは含まれていません。

f:id:Kurene:20210117210802p:plain — Zwicker Loudness Measurement with R&S® UPV Audio Analyzer Application Note: Fig. 2-6:

25 Hzから12.5 kHzまでの中心周波数で1/3オクターブ分析
低周波での耳の感度が低いことを考慮し、低周波の1/3オクターブレベルを重み付け
低周波（25 ～ 250 Hz）の1/3オクターブ分析のエネルギーを組み合わせ
臨界帯域レベルの算出
- 人の耳の伝達特性に応じた補正
- 拡散音場音の場合の補正
- 1/3オクターブ分析から得られた周波数帯域の実際の臨界帯域幅からの乖離に対する臨界帯域レベル補正
臨界帯域のレベルをコアラウドネスに変換
バーク尺度 (critical band rate) へのマッピング
聴覚の周波数分析（選択性）の非対称性を考慮し、より高い臨界帯域への傾きを追加
- スロープの急峻さは、それぞれの臨界帯域のコアラウドネスに依存している
各臨界帯域の傾きとコアラウドネスのうち大きい方を選択し、部分ラウドネス（ラウドネス密度）を算出
部分ラウドネスを臨界帯域（周波数）方向に積分し、ラウドネスを算出

MathWorks - acousticloudness

f:id:Kurene:20210117215135p:plain — https://jp.mathworks.com/help/audio/ref/acousticloudness.html

時間領域の信号レベルを、CalibrationFactorに従って調整
信号を1/3オクターブSPL表現に変換
- フィルタバンクは、25 Hzから12.5 kHzの間の28個のフィルタで構成
- このステージからの出力はdB単位であり、また、基準圧力で正規化済み
バンド依存の、時間方向平滑化フィルター
帯域重み付け、および帯域結合
- 低周波の1/3オクターブ帯は、固定の重み付けテーブルに従って非強調( de-emphasized) する
- 低周波帯域のいくつかは、合計20の臨界帯域を形成するために結合する
レベル補正、コアラウドネス算出
- 臨界帯域レベルは、フィルタの帯域幅と、静寂のしきい値での臨界帯域レベルで補正する
- 音場のタイプを考慮してレベル補正
非線形時間減衰
- diode-capacitor-resistor networkを用いてシュミレート
- 長い信号と比較して、短い信号の後の急峻な知覚低下をモデル化している
コアラウドネスをBark尺度にマッピング
レベルと周波数に依存した傾きの表を用いて周波数の広がりを計算
ラウドネスへの時間重み付け
- ラウドネス知覚の持続時間依存性をシミュレートするため
- 部分ラウドネスを算出
周波数拡散のスロープを考慮し、部分ラウドネスの積分としてラウドネスを計算

Wizard Notes

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

PythonとMoSQIToで学ぶ、ラウドネス (sone) の算出および実装方法 (Zwicker method)

はじめに

`MoSQITo`を使ったラウドネス算出のソースコード

各信号処理ブロックの説明

1. 1/3オクターブバンド分析

2. 臨界帯域レベル（コアラウドネス）の算出

3. 時間マスキング (Temporal_masking)

4. バーク尺度へのマッピング，周波数マスキング

5. ラウドネスへの時間重み付け

まとめ（＋疑問）

補足：参考文献のラウドネス算出ブロック図／説明

Zwicker Loudness Measurement with R&S® UPV Audio Analyzer Application Note

MathWorks - acousticloudness

はじめに

MoSQIToを使ったラウドネス算出のソースコード

各信号処理ブロックの説明

1. 1/3オクターブバンド分析

2. 臨界帯域レベル（コアラウドネス）の算出

3. 時間マスキング (Temporal_masking)

4. バーク尺度へのマッピング，周波数マスキング

5. ラウドネスへの時間重み付け

まとめ（＋疑問）

補足：参考文献のラウドネス算出ブロック図／説明

Zwicker Loudness Measurement with R&S® UPV Audio Analyzer Application Note

MathWorks - acousticloudness

`MoSQITo`を使ったラウドネス算出のソースコード