インパルス応答から残響時間を計算する(python)

  • by

今回はインパルス応答から残響時間を計算する方法を紹介します。

言語は相変わらずpythonでやっています。
音の入出力という意味でも、信号処理という意味でも、
圧倒的にMATLABの方が便利ですし、私個人の所属機関の関係で
無料で使える環境にはあるのですが、
やっぱりこれからもタダで使えてオープンな言語のほうが
いいだろうということでpythonでの公開になります。

Contents

  1. 残響時間について
  2. 求め方
  3. 実装
  4. まとめ

1. 残響時間について

コンサートホールなどの建築物について詳しい人や、
かつての高度経済成長時代にオーケストラなどに興味があった人なら
「残響時間2秒がいい」という謎の(?)フレーズをご存知かと思います。

残響時間とは、とっても平たく言うと物理的な「残響の量」です。
音が発せられて、停止した後に残っている音のエネルギーの量を
時間で表したものです。

もっと具体的に言えば、
音が鳴り終わったあと、一定の量減衰するまでにかかる時間になります。

より正しい定義では、
残響時間とは減衰曲線が60dB低下するまでの時間を指します (Fig. 1)。

Fig. 1 残響時間の定義

数字の目安として、一般的な居室では0.3秒、大型ホールであれば0.7~2秒、
大きな教会や残響室では10秒くらいだといわれています。
特殊な例だと建設されて間もなくまだ開通していないトンネルでの
測定で30秒くらい(うろ覚え)という調査もあります。

前回の記事では明瞭度という心理量との対応のある指標を紹介しましたが、
残響時間という物理量は心理量である「残響感」と対応があるわけでは
ありません。残響感と対応があるのは初期残響時間です。
残響時間はあくまで物理量であるという点で注意が必要です。

表記の仕方は「T」や「RT」「RT60」のような書き方がメジャーだと思います。

2. 残響時間の求め方

残響時間を求める方法は主に二つあります。
ノイズを放出して停止した後の減衰を直接見るノイズ断続法と、
インパルス応答から色々やって求める方法があります。
今回は後者です。

残響時間を求める流れを以下に示します。

  1. インパルス応答を得る(必要に合わせてオクターブバンドフィルタをかける)
  2. シュレーダー積分式を使って減衰曲線を得る
  3. 一定の区間で回帰分析をして傾きを求める
  4. 残響時間を計算する

1のインパルス応答は特に問題ないと思います。
過去の記事も参考になりそうです。

2で突然、シュレーダー積分式とかいうものが出てきます。
数式での証明はしませんが、これの何がうれしいかというと、
ノイズ断続法では何度も減衰曲線を求めて平均を出した(同期加算した)上で
傾きを見るのですが、これを使うと平均された状態の減衰曲線が出てきます。

$$ \langle S^{2}(t) \rangle =\int_{0}^{\infty}h^{2}(t) {\rm dt} – \int_{0}^{t}h^{2}(t){\rm dt}$$

$$ \langle S^{2}(t) \rangle : 減衰波形の集合平均, h^{2}(t): インパルス応答の二乗波形 $$

3、4についてですが、60dBの減衰を考える場合には、
必然的に60dB以上のSN比が必要となります。
しかし、実際そこまで綺麗にとれる場合がないので、
30dBの減衰にかかった時間を2倍したり( \(RT_{30}\) )、
20dBの減衰にかかった時間を3倍したり ( \(RT_{20}\) )、します。
いずれの場合も傾きを求めることと同じなので、回帰分析でできちゃいます。

そしてもうひとつ、残響時間という概念は拡散音場というおおざっぱで
(割と)非現実的な仮定のもと成り立っています。
詳細は省きますが、直接音や偏りのある音場の影響が減衰開始直後に
現れるため、回帰する範囲は-5dB~-35dBのように-5dBスタートである
ことがほとんどです。

3. 実装(プログラム)

またもや今回も室内音響指標ベンチマークのデータ(A-11)を拝借。
500 [Hz]帯での残響時間 ( \(T_{30}\) ) はベンチマーク上だと0.49~0.50秒です。

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
decay_lebel_dB = 30
times = 60 / decay_lebel_dB

# 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_bpf = sosfilt(bpf_coef, ir)

# decay curve
ir_bpf_square = ir_bpf ** 2.0
ir_bpf_square_sum = np.sum(ir_bpf_square)
temp = ir_bpf_square_sum
curve=[]
for i in range(len(ir)):
    temp = temp - ir_bpf_square[i]
    curve.append(temp)
curve_dB = 10.0 * np.log10(curve)
curve_offset = max(curve_dB)
decay_curve = curve_dB - curve_offset

# find regression target
i = 0
while decay_curve[i] > -5.0:
    i += 1
start_sample = 1
while decay_curve[i] > -1.0*decay_lebel_dB - 5.0:
    i += 1
end_sample = i
regression_target = decay_curve[start_sample:end_sample]

# linear regression for T
x = np.linspace(start_sample, end_sample, end_sample-start_sample)
a, b = np.polyfit(x, regression_target, 1)
rt_sec = (-1.0*decay_lebel_dB/a)*times/fs
print(rt_sec)

実行するとおよそ0.49と出力されるので良さそうです。

4. まとめ

残響時間の概念や求め方についての説明、
インパルス応答から残響時間を計算するプログラムの公開、
そしてベンチマークとの比較をしました。


Rosa foetida

ロサ・フェティダ

黄薔薇さまで通じる方が果たして
どれほどいらっしゃるのだろうか

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です