残響時間を求める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))
}