←戻る

残響時間(T30)の計算 / Computing reverberation time (T30)

ISO 3382-1にしたがって作成したつもりですが、間違いがあるかもしれません。

出力のグラフは、黒がインパルス応答の振幅の2乗波形（デシベル）、緑が残響曲線、赤の破線が(-5,\mathrm{dB})〜(-35,\mathrm{dB})区間の近似直線です（fit.startとfit.endで始まる2行を編集すると近似範囲の指定を変更できます）。近似直線から計算した残響時間はグラフタイトルに書かれています。（低域は信頼性が低いですね）

library(tuneR)
library(signal)

# read impulse response WAV file
wavobj <- readWave(file.path("IR_file.wav"))
fs <- wavobj@samp.rate
nbits <- wavobj@bit
x <- wavobj@left / 2^(nbits-1)

# find direct sound (which starts 20dB below the peak) only for drawing nice x-axis
xx <- 10*log10(x^2)
for (t0 in which.max(xx):1) {
  if (xx[t0] < max(xx)-20) {
  	break;
  }
}

# octave band filter definition (1024-point FIR)
cfreqs <- c(63, 125, 250, 500, 1000, 2000, 4000, 8000)
N <- 1024

for (fc in cfreqs) {
  # apply octave band filter
  fl <- fc / 2^(1/2) / (fs/2)
  fh <- fc * 2^(1/2) / (fs/2)
  k <- fir1(N, c(fl, fh), type="pass")
  y <- fftfilt(k, append(x, rep(0, N)))
  t1 <- t0+N/2
  t <- (-t0:(length(y)-t0-1))/fs

  # Schroeder curve and reverb time (T30)
  sc <- cumsum((y[length(y):1])^2)
  sc <- 10*log10(sc[length(sc):1] / max(sc))
  fit.start <- sum(sc >= -5)
  fit.end <- sum(sc >= -35)
  z <- line(t[fit.start:fit.end], sc[fit.start:fit.end])
  T30 <- -60 / coef(z)[2];

  # plot the result
  yy = 10*log10((y/max(abs(y)))^2)
  plot(t, yy, type="l", ylim=c(-80,0),
       main=sprintf("Reverberation Curve (%dHz, T30=%.2fs)", fc, T30),
	   xlab="Time (s)", ylab="Power (dB)")
  lines(t, sc, col="green")
  abline(coef(z), col="red", lty="dashed")
  grid()

  # save to PDF file
  dev.copy2pdf(file=sprintf("reverb_time_%05d.pdf", fc))
}