あけましておめでとうございます。R Advent Calendar 2013の32日目の記事です。
僕の専門は音響と心理との学際領域で、ふだんは音に関する主観評価実験の準備・分析にMatlabやRを使用しています。今回は室内音響分野で使われる計算をRでやってみます（その分野だと、みなさんC/C++やMatlab、最近だとPythonとか使っているみたいです）。Advent Calendar経由で来た方にとっては耳慣れない言葉が出てくるかもしれません。また、室内音響については単純化した説明をしているので専門の方にはアラが気になるかもしれません。ご容赦を。
室内音響指標
さてここから本題。年末年始でコンサートに出かけた方もいるかもしれません。コンサートホールは、ステージで音が発せられて、それを鑑賞するという、たったひとつの目的のための建物なのに、室内の形状は様々です。これは、それぞれのホールが目的とする音楽の様式だったり、演奏だけじゃなく講演会にも使いたいという狙いだったり、収容人数や土地面積や建築方法や予算制約など、いろいろなものに影響を受けて形状が違ってくるのですね。
形状が違うことにより、同じ曲の演奏でもコンサートホールごとに響き方が違ってきます。「音が良いと評判のホールはなぜ音が良いのか」や「あのホールとこのホールはどこが違うのか」を説明できる指標（数値）があるといいなぁということで、いろいろと室内音響指標が提案されています。特にメジャーなのは、ホールの物理測定データ（インパルス応答）から室内音響指標を計算するもので、測定方法や計算方法が「ISO 3382-1:2009 Acoustics - Measurement of room acoustic parameters, Part 1: Performance spaces」という国際規格にまとめられています。
インパルス応答
インパルス応答というのは、ホールの中で時間極短・エネルギー極大な音を出したときに、ホール壁面などから戻ってくる響きを録音したものだと考えて差し支えないと思います。ホールの場合はステージ上で音を出し、それが客席にどのように届くのかを測定することが多いです。
ホールのようにある程度大きな空間の場合は下図のように、まずステージからまっすぐに届いた音（直接音）、続いて床や壁で反射した音（初期反射音）、それから反射音どうしの区別がつかなくなった残響音の3つの部分に分けて考えます。直接音は1つだけ。初期反射音と残響音の境目はホールの大きさ・形状などによって違うものの、80〜100ミリ秒が一般的です。
C値
さて今回は、人間が感じる「音の明瞭さ」に対応すると言われるC値を計算してみます。C値はインパルス応答の前半と後半のエネルギー比をデシベルで表したもので、以下のような式で計算されます。
ここでp(t)は時刻tにおける音圧、teは前半と後半を分割する時刻です*1。前半というのは直接音と初期反射音を含む部分、後半は残響音の部分とするので、上記の80〜100ミリ秒などを使用します。log中の分数では、分子が前半、分母が後半のエネルギーを計算しています。デシベルに変換するために、常用対数をとって10をかけています。（実際にはp(t)をオクターブ帯域ごとに分割して、それぞれの帯域でC値が計算し、500〜1000Hzの2帯域の算術平均値を使ったりします。）
ISO 3382-1によると、te=50ミリ秒だとスピーチ、te=80ミリ秒だと音楽の「明瞭さ」によく対応するそうで、だいたいC値は±5dBの範囲におさまることが多いとのこと。C値が正の値だと直接音＋初期反射音が強いため、残響音にかき消されずにステージからの音がはっきり聞こえるだろう、というのがこの式の意味するところです。*2
Rでの実装
インパルス応答はモノラルのWAVファイルとして保存されているとします。tuneRというパッケージを使うとWAVファイルの読み書きができるので、使わせてもらいます。
# パッケージの読み込み
library(tuneR)
# WAVファイルとして保存してあるインパルス応答を読み込む
wavobj <- readWave("IR.wav")
fs <- wavobj@samp.rate
nbits <- wavobj@bit
x <- wavobj@left / (2^(nbits - 1))
t <- seq(from=0, to=length(x)-1) / fsWAVファイルを読み込みますが、16ビットの場合はデータ範囲が-32768〜+32767となっているので、それを-1〜+1の範囲に変換しています。24ビットの場合もあるかもしれないので、とりあえず決めうちはしないようにしています。それぞれのサンプルの時刻を計算して、tというベクトルも作りました。
# 直接音はピークから20dB下がったところから開始するとして、そのポイントを探す
x.dB <- 10 * log10(x * x)
peak <- which.max(x.dB)
for (n in peak:1) {
if (x.dB[n] < -20) {
peak.start <- n
break;
}
}WAVファイルの先頭が必ず直接音であるとは限らないので、直接音がどの時刻にあるかを探さなければなりません。インパルス応答は直接音が最大値であることが多いので、今回もそれを仮定して瞬間エネルギーが最大になる時刻を探します。ただし、そのピークを計算開始点とすると、エネルギーがピーク前半と後半のふたつに分かれてしまいますので、ここでは「最大値から前方に20dB下がったところから直接音の立ち上がりが開始する」として探索しています。
# C値の計算（今回は80ミリ秒（＝0.080秒）を使用） t.value <- 0.080 x.early <- x[1:(peak.start+round(t.value*fs))] x.late <- x[(peak.start+1+round(t.value*fs)):length(x)]
そして、今回は80ミリ秒のところで、前半と後半に分割しました。前半・後半を色分けしてプロットすると以下のようになります。前半には大きな直接音の後に、数発の初期反射音を見ることができます。後半は突出した音はなく、なだらかに減衰しています。
ここから、オクターブバンドごとに分割して、それぞれの帯域のC80値を求めます。帯域分割をするのにバターワース・フィルタを使用しますが、signalパッケージに入っているものを使います。
library(signal)
cfreqs <- c(125, 250, 500, 1000, 2000, 4000, 8000)
for (cf in cfreqs) {
bf <- butter(4, c(cf/2^(1/2), cf*2^(1/2))/(fs/2))
y.early <- filter(bf, x.early)
y.late <- filter(bf, x.late)
Ct <- 10*log10(sum(y.early^2) / sum(y.late^2))
cat(sprintf("%5dHz: %+.2fdB\n", cf, C80))
}中心周波数は125Hz〜8000Hzとしました。1オクターブ（＝周波数2倍）の幅にしたいので、中心周波数を21/2で除算／乗算してあげます。たとえば1000Hzを中心とした1オクターブバンドは、1000÷21/2=707.1Hzから1000×21/2=1414.2Hzです。707Hzと1414Hzとがちょうと周波数2倍になります。フィルタ次数は適当に4次としておきました。
インパルス応答の前半・後半それぞれにフィルタをかけ、C値の数式どおりに計算すると、こんな感じにC80の計算結果が表示されました。
125Hz: -0.20dB
250Hz: -3.24dB
500Hz: -0.88dB
1000Hz: -0.50dB
2000Hz: -2.03dB
4000Hz: -1.43dB
8000Hz: +2.35dB
代表的な値ひとつだけを報告したいときは500Hz帯の-0.88と1000Hz帯の-0.50を算術平均した-0.69dBを使うことになります。
まとめ
僕はふだんの信号処理にはMatlabを使っているのですが、Matlabは学生さんにお勧めできる値段じゃないし*3、かといってGNU Octaveはインストールが煩雑でやっぱりお勧めしにくい*4。しかも実験準備の段階では信号処理、実験後の分析には統計処理、論文執筆時にはグラフ作成をしますが、こういった用途ごとに異なる言語を覚えたくない、という欲張りをかなえてくれる環境はないかと探した結果、Rにたどり着きました。
なので、まずはRを統計分析で使っていましたが、このところは信号処理などにも使い始めています。今回紹介したC値の計算も、そうしたRで信号処理をやってみる一環です。Rで本格的な信号処理をするには処理速度が遅くて涙が出ますが*5、ちょっとした波形の切り貼りだったりグラフ化など、実験の準備から分析までがひとつの流れの中で実行できるのは有り難いです。
*1:te：earlyはここまで、を表す時刻t、ということみたい。
*2:C値の式：現在、C値の妥当性・有効性について議論が続いています。確かに、この式は音の到来方向を無視しているし、エネルギー比だけを計算するので位相ズレが起きていても値が変化しないとか、いろいろと問題はあります。
*3:Matlab学生版がありますが、それでも1万円程度しますし、卒業してからも使い続けたいとなると大人版購入は相当な出費です。
*4:数年前まではプリコンパイルされたバイナリが配布されていたんですが、このところはHomebrewなどのパッケージマネージャが推奨されているようです。しかも頻繁にパッチが変更されてインストールできない状況になり、結構苦労します。
