残響時間を求めるRのプログラム

今日はR言語がどのていど信号処理っぽい目的に使えるか試してみました。とは言ってもガリガリと信号処理をするわけではなく、簡単にフィルタがかけられて、結果がグラフにできればいいので、室内のインパルス応答から残響時間を求めるプログラムを作成してみました。以前にMatlabで作成したことがあるので、それを思い出しつつ、解説系のウェブサイトは参考にせずRのhelpコマンドのみ利用するルールで進めました。

tuneRパッケージでWAVファイルの読み書きができるのは知っていましたし、signalパッケージで簡単な信号処理ができるのはわかっていたのですが、思ったより時間がかかってしまいました。

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 is 20dB below the peak)
xx <- 10*log10(x^2)
for (t0 in which.max(xx):1) {
  if (xx[t0] < max(xx)-20) {
  	break;
  }
}

# apply octave band filter
cfreqs <- c(63, 125, 250, 500, 1000, 2000, 4000, 8000)
N <- 1024
for (fc in cfreqs) {
  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(abs(y[length(y):1]))
  sc <- 20*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))
}

以前に作成したプログラムとほぼ同じ残響時間が計算できています。おそらくR言語の流儀からはだいぶ外れた、無駄の多いプログラムでしょうが、とりあえず目的達成。めでたしめでたし。

下図、黒い線はインパルス応答のエネルギー、緑の実線はシュレーダー曲線（残響曲線）、赤の破線は-5dBから-35dBでの回帰直線です。