インパルス応答の分析（直間比）

昨日の続き。直接音と間接音のレベル比を計算します。ISO 3382という規格にclarityという物理指標が出てきます。インパルス応答の前半部分（直接音に対応）と後半部分（間接音に対応）のレベル比をdBで求めるもので、例えば以下のような式で表されます。

$C_{80} = 10\log_{10}\left(\frac{\int_{0\,\mathrm{ms}}^{80}p^2(\tau)d\tau}{\int_{80\,\mathrm{ms}}^{\infty}p^2(\tau)d\tau}\right)$

この式は直接音から80ミリ秒を境目として、その前後の音圧比を求めます。80ミリ秒はだいたい初期反射音と残響音との境目だということです。C80ではなくC50やC30を使えば、室の大きさによって異なる反射音と残響音の境目に対応することが出来ます。また、C1とかC0.1などのように、小さい値を用いれば直接音と間接音の境目にすることもできるでしょう。

ただし、この式を使うためには、インパルス応答の直接音に相当する部分を自動的に見つけなければいけません。こんなかんじでやってみました。

[z, zind] = max(log10(abs(diff(x(1:fs)))) > -2);   % roughly look for the first peak
z = diff(x(1:fs));
for n=zind:length(x)
  if z(n-1) * z(n) < 0, break;   % find zero crossing of the difference of the impulse response
  end
end
zind = n;   % zero crossing is where the direct sound is

信号波形の微分値が初めてゼロを上から下に挟む場所（波形の傾きがプラスからマイナスに変化する極大点）が直接音のピークなのですが、単純に波形の先頭から計算してしまうとノイズにだまされてしまいます。なので、まず粗くピークを見つけておいて、そこから細かく見ていくわけです。

あとは上記数式どおりに計算をすればOK。

sepa = floor(80 * fs / 1000);
Cnumer = sum(x(zind:zind+sepa) .^ 2);
Cdenom = sum(x(zind+sepa+1:end) .^ 2);
C = 10 * log10(Cnumer / Cdenom);