はじめに

音声認識において、学習用データを水増し(Augumentation)する際に
任意の残響時間のリバーブを付加したくなったので、
任意の残響時間から残響曲線を求めてリバーブフィルタを作成しました。

インパルス応答から残響曲線を求める方法

インパルス応答から残響曲線を求めるとき、
通常の場合、インパルス応答2乗積分法[1]を用います。

インパルス応答2乗積分法(Integrated impulse response method)は、インパルス応答から下式を用いて残響エネルギーの減衰曲線を求める方法です。この積分は開発者の名前からシュレーダ積分法(Schroeder integration method)ともいいます。

E(t) = \int^\infty_t p^2 (\tau)dt　\tag{1}

(1)式の$p(t)$にインパルス応答を代入すると、残響曲線$E(t)$が求められます。
(インパルス応答は音圧の次元、残響曲線はパワーの次元です。)

シュレーダ積分を解く

(1)式を$p(t)$について解くと、任意の残響時間を模擬したインパルス応答が求められ（る気がし）ます。

E(t) = \frac{1}{3} \bigr[p^3\bigr]^\infty _t \\
E(t) = \frac{1}{3} p^3(\infty) - \frac{1}{3} p^3(t)\\
\frac{1}{3}p^3(t) = p^3(\infty) - E(t)

$p(\infty) = 0$とすると、

p^3(t) = -E(t) \\
p(t) = -E^{\frac{1}{3}}(t)

残響曲線を作る

任意の残響時間を模擬した残響曲線を作りましょう。
残響時間とは、以下で定義されます[2]。

残響時間は、音源が発音を止めてから、残響音が60デシベル減衰するまでの時間をいう。

すなわち、t=0のとき0dB、t=[残響時間]のとき-60dBとなるような減衰を擬似的に作れば良さそうです。
簡単のため、残響曲線が直線を描くとすると、
残響時間がRTのとき、擬似残響曲線$E'(t)$は、

E'(t) = -\frac{60}{RT} t

で表せます。

Pythonで実装

必要な材料が整ったので、Pythonで実装してみます。

# import modules
import numpy as np

# define const.
rt = 0.1
fs = 48000

t = np.linspace(0, ir_len, int(ir_len*fs))
E = np.power(10, -3 * t / rt)

reverb = np.power(E, 2/3)
rand = np.random.randint(0, 2, reverb.size)
rand = np.where(rand==1, 1, -1)
reverb *= rand

例えば、RT=0.3sのとき、以下のようなインパルス応答が得られました。

まとめ

簡易的に残響曲線から残響付加インパルス応答を算出する方法について紹介しました。
あくまで簡易的なので、初期反射音・後期反射音の反射音密度の違いや、
残響時間の周波数による違いまでは再現できていませんが、
簡単な再現でも良い場合は、使えるかもしれません。

参考文献

[1] インパルス応答から求める室内音響指標、羽入敏樹、日本音響学会誌69巻6号(2013)
[2] 残響時間 - Wikipedia