室内音響指標算出をpythonで
インパルス応答は対象とする系の情報を持っています。
中でも室に対するインパルス応答は音の反射の情報を持っており、
反射音の集合である残響は音の聴こえに少なくない影響を与えます。
室内で測定したインパルス応答から聴こえを評価する指標は
いくつもあり、明瞭度はそのひとつです。
Contents
- 明瞭度の定義
- 算出フロー
- 実装
- 室内音響指標ベンチマークとの比較
- まとめ
1. 明瞭度の定義
明瞭度、C値、Clarity、もしくはC50やC80と呼びます。
明瞭度はインパルス応答を前半と後半に分けて(Fig. 1)、
それぞれのエネルギーの和の比で定義されています。
前後半に分ける時間は50[ms]と80[ms]が主に用いられています。
式で表すと以下の通り。
\[C_{t_{e}}=10{\rm log}{10}\frac{\int_{0}^{t_{e}} p^{2}(t)dt}{\int_{t_{e}}^{\infty} p^{2}(t)dt}\] \[C_{t_{e}}: 明瞭度[{\rm dB}]，{t_{e}}: 前後半の境界(50[{\rm ms}] か 80 [{\rm ms}])，p(t): インパルス応答波形\]50 [ms]は音声の明瞭度、80 [ms]は音楽の明瞭度に関係があるとされています。
明瞭度の提案者によるオーケストラ演奏における心理実験 (Reichardt, 1974) の
結果から、曲によって異なるがC80において0[dB]～-3[dB]が
適切であるとしています。
また居室における会話に着目した調査(佐藤，2010)では、
聴取においてC50が3.1[dB]以上、発話において12.6[dB]以上を、
設計目標の値として提案しています。
一般的に同一の室だと±5dBに収まるとかなんとか…
音源と受音点の位置でだいぶ変わったりするので、
室の用途や、解析したい内容に合わせて適宜場所を変えたり、
複数点でとって平均したりします。
2. 明瞭度の算出フロー
- 直接音到来時間を推定
- インパルス応答を直接音到来時間を起点（＝0）として前後半に分割
- 分割して短くなった分を0埋め
- オクターブバンドパスフィルターをかける
- 2乗して和を求める
- 比を求める
明瞭度は1/1オクターブバンドフィルターをかけて表すことがほとんどです。
各中心周波数ごとにプロットしたり、ひとつの値として出す場合は中心周波数
500[Hz]と1000[Hz]で明瞭度を求めて算術平均の結果をとったりします。
詳細は室内音響指標ベンチマークのコラムなども参考にすると良いです。
3. 実装（pythonプログラム）
入力：
ir インパルス応答 (array like)
必要なデータ：
fc 分析対象オクターブバンドの中心周波数 (int)
fs インパルス応答のサンプリング周波数 (int)
rate 直接音到来時間推定におけるしきい値 (int)
bound_ms 前後半を分ける時間(50[ms] or 80[ms]) (int)
出力：
clarity 明瞭度 (int)
分析対象は室内音響指標ベンチマークの基本問題A-11です。
import numpy as np
from scipy.io.wavfile import read
from math import log10
from scipy.signal import sosfilt, butter
# read data
fs, ir = read("a11.wav")
# analysis configure
fc = 500
bound_ms = 80
rate = 20
# 1. estimate direct sound arrival time (sample)
ir_square = ir ** 2.0
ir_square_max = max(ir_square)
i=0
arrival_sample = 0
ir_square_rate = ir_square[i]/ir_square_max
for i in range(len(ir)):
if ir_square[i] != 0.0:
arrival_sample = i
ir_square_rate = ir_square[i] / ir_square_max
if 10.0*log10(ir_square_rate) > -1.0*rate:
break
else:
continue
# 2. split signal
bound_sample = int(bound_ms*0.001*fs)+arrival_sample
ir_first = ir[arrival_sample:bound_sample]
ir_latter = ir[bound_sample:]
# 3. zero padding
analysis_length = len(ir[arrival_sample:])
first_pad = np.zeros(len(ir_latter))
ir_first_pad = np.hstack((ir_first, first_pad))
latter_pad = np.zeros(len(ir_first))
ir_latter_pad = np.hstack((latter_pad, ir_latter))
# 4. band pass filter
nyquist = 0.5*fs
low_cutoff = fc/(2**(1/2)) / nyquist
high_cutoff = fc*(2**(1/2)) / nyquist
bpf_coef = butter(4, [low_cutoff, high_cutoff], btype='bandpass', output='sos')
ir_first_pad_bpf = sosfilt(bpf_coef, ir_first_pad)
ir_latter_pad_bpf = sosfilt(bpf_coef, ir_latter_pad)
# 5. calcurate square sum
c_numer = np.sum(ir_first_pad_bpf**2)
c_denom = np.sum(ir_latter_pad_bpf**2)
# 6. obtain clarity from log
clarity = 10.0 * np.log10(c_numer/c_denom)
print(clarity)
コメントの数字は算出フローでのナンバリングと対応しています。
1の直接音の到来では、インパルス応答2乗波形の最大値を基準として、
-20[dB]以下、かつそれまでで最大の点を直接音到来時刻としています。
これはISOの規格によって定められていますが、実際はインパルス応答自体の
SN比によって調整する必要があります。
2では直接音到来～50or80[ms]とそれ以降で切り分けています。
3のゼロ埋めはバンドパスフィルタによる遅延した分のエネルギーを、
計算に含めるために追加しています。
4のバンドパスフィルタには、バタワースフィルタを用いました。
バタワースフィルタは通過帯域内の振幅がフラットになるという特性があり、
python上ではカットオフ周波数からのフィルタ係数の生成(butter)と、
係数にもとづく2次セクション型IIRのフィルタリングをしています(sosfilt)。
5でそれぞれlog10の真数部分にあたる分子分母を計算し、
6で比として出るようにlog10の計算をしています。
4. 室内音響指標ベンチマークでの評価
まだ公開はしていませんが、上で挙げたようなプログラムを使いやすいように
pythonモジュールとしてまとめています。
その開発中のもの（今回紹介したプログラムと同じアルゴリズム）を、
使って室内音響指標ベンチマークと比較してみました。
室内音響指標を計算するプログラム作った時、その精度を確認するための
ベンチマーク問題が日本建築学会から「室内音響指標ベンチマーク」として、
公開されているので、これを利用しました。
125Hz付近が危うい。
全体を通して低域が低めに出て、中心周波数ごとの差が少し
減っているような感じがしますが、おおむね一致しています。
5. まとめ
今回はインパルス応答から明瞭度を計算する方法を、具体的な
プログラムと一緒に紹介しました。
こうやって計算する、みたいな感じで計算式や望ましい信号処理の仕方などは
色々見つけられますが、具体的なコードを交えて議論できる場が
あまりないので、個人的にはどんどん情報共有されるといいなと思います。
P.S. お知らせ
7/3からさくらレンタルサーバのPHPのバージョンアップに伴い、
プログラムを掲載しているページが正しく表示されない不具合が
発生しております。
シンタックスハイライトを表示できるプラグインである
「Crayon Syntax Highlighter」が原因であり、
適宜変更していく予定です。（本記事は大丈夫です。）
菜の花
葛西臨海公園にて撮影
菜の花は花としての可愛さもあるけど、おひたしが美味しくて好き