インパルス応答の分析（残響時間を求める）

ここに書いた方法はISO 3382に記載されたものとは異なり、本当の残響時間とは異なる値が出てしまいます。正しい計算をするプログラムは「残響時間を求めるRのプログラム - 丸井綜研」にあります。

Matlabで部屋のインパルス応答の分析をしています。分析するのには多種多様な方法があるのですが、まずは残響時間の計算ができて、直接音と間接音の比が求まれば、とりあえずの分析材料にはなります。残響時間を求めるために、以前作成した残響曲線の計算をするMatlab用のプログラムからご紹介。

function y = reverb_envelope(x)
% REVERB_ENVELOPE   Calculates reverberation envelope
%    y=reverb_envelope(x) calculates reverberation envelope y from
%    an impulse response x.  This function can be used to see how
%    the reverberation decays in the corresponding room.
%
%    Reference: Ooga, Yamazaki, and Kaneda
%              "Acoustic Systems and the Digital Processing" (p.163)
%
% 2007-05-28 MARUI Atsushi

if ndims(x)>2
  return;
end

if size(x,1) < size(x,2)
  x = x';
end

y = flipud(cumsum(flipud(x .^ 2)));

これをもとにグラフを書いて、マイナス60dBになる部分を探せばOKです。具体的には

re = reverb_envelope(x);
mm = max(10*log10(re));
t = [0:length(re)-1]/fs;
z = 10*log10(re)-mm;
ind = sum(z >= -60);

のようにやれば、indに何サンプル目がマイナス60dBか、というのが入ります。tにはファイル先頭からの時刻が秒で入っているので、t(ind)が残響時間です。ただ、残響曲線の形によっては、-10dBの点と-30dBの点などを直線で外挿してやる、などの処理が必要になる場合もあります。

直接音と間接音の比は明日。

2008年11月18日追記

残響時間の計算部分がイマイチ正確ではなかったので、2008年11月18日のエントリーで修正しました。