まず、分析対象となる音ファイルを読み込みます。例で使ったのはギターの音です。
library(tuneR)
# ファイルの読み込み
wavobj <- readWave("guitar.wav")
x <- wavobj@left # 左チャンネルだけを使う（wavobj@.Dataを使う場合も）
fs <- wavobj@samp.rate # 標本化周波数
nbits <- wavobj@bit # 量子化ビット数
x <- x[1:(fs*5)] # 先頭の5秒間だけを分析対象として抜き出す
試しに波形を見てみます。
t <- seq(from=0, to=length(x)-1) / fs # 各サンプルの時刻を計算
plot(t, x, type="l", xlab="Time (s)", ylab="Amplitude", ylim=c(-(2^(nbits-1)), (2^(nbits-1))-1))
grid()
振幅は量子化ビット数が16なので、(-(2^{16-1})) 〜 (2^{16-1}-1)（つまり-32768〜32767）の範囲を縦軸に設定しています。
音信号の周波数特性を見たい時にはフーリエ変換を用います。Rでもfft()を用いるとフーリエ変換を行うことができます。
# フーリエ変換
y <- fft(x)
RでFFTを行うと、入力（時間波形x）と出力（周波数応答y）のデータ長は同じになります。しかし、yは以下のようにナイキスト周波数での応答を対称にして周波数応答とその複素共役の形になっています。
試しに絶対値をプロットしてみると左右対称になっていることが分かります。
y.tmp <- Mod(y) # 複素数の絶対値を計算
plot(y.tmp, type="h")
詳しくは述べませんが、これはDC成分とナイキスト成分以外のエネルギーが左右に分かれている状態なので、それを考慮して以下のように振幅スペクトルを計算します。
# 振幅スペクトル（amplitude spectrum）の計算
y.tmp <- Mod(y) # 複素数の絶対値を計算
y.ampspec <- y.tmp[1:(length(y)/2+1)] # DC〜ナイキストの範囲
# 対応する周波数の計算（DC〜ナイキストの範囲をビン数で均等割）
f <- seq(from=0, to=fs/2, length=length(y)/2+1)
# グラフ作成
plot(f, y.ampspec, type="h", xlab="Frequency (Hz)", ylab="Amplitude Spectrum", xlim=c(0, 5000))
grid()
今回は楽器音を分析しましたから、基音と倍音成分がきれいに並んでいる様子を見ることができます。
上記は横軸を周波数、縦軸を振幅でプロットしましたが、両対数軸にしたほうが分かりやすいこともあります。以下では、横軸はlog="x"で対数にし、縦軸は振幅値を2乗してパワーにしてからデシベルに変換し、実質両対数軸にしています。
y.dB <- 10*log10(y.ampspec ^ 2)
plot(f, y.dB - max(y.dB), log="x", type="l", xlab="Frequency (Hz)", ylab="Power (dB)", xlim=c(20, 20000), ylim=c(-100, 0))
## Warning in xy.coords(x, y, xlabel, ylabel, log): 1 x value <= 0 omitted from logarithmic plot
grid()
周波数fはDC（0 Hz）から始まりますが、「0の対数は定義できない」という警告が表示されています。グラフの表示範囲はxlimで20〜20000 Hzと指定していますので、警告は無視します。
# 逆フーリエ変換で時間波形に戻す
x2 <- Re(fft(y, inverse=TRUE) / length(x))
# 原信号との差はほとんどないことを確認
plot(x - x2, type="l") # 誤差をプロット
10*log10(sum((x - x2)^2)) # 誤差の総和をデシベルで計算
## [1] -178.395
計算誤差程度の差があることが確認できますが、-178 dB程度なので大きな問題になることはないでしょう。
ある周波数を境に低域のエネルギと高域のエネルギがつりあうとき、その周波数のことをスペクトル中心（spectral centroid）と呼びます。スペクトル中心は、音の明るさと対応のよい音響指標値として用いられています。
一般には下式のようにFFTの各ビンの中心周波数(f)とそれに対応するスペクトル振幅値(p)を用いて計算が行われます（(n)はビン番号）。
[ f_\mathrm{centroid} = \frac{\sum_{n=1}^{N}f(n)\cdot p(n)}{\sum_{n=1}^{N}p(n)} ]
周波数を線型にとるか対数でとるか、スペクトルを振幅でとるかパワーでとるかで計算結果が大きく異なってきます。いろいろな考え方があるとは思いますが、ここでは、対数周波数とパワースペクトルを用いてスペクトル中心を計算します。（これが感覚的な音の明るさに直感として合う気がしています。）
library(tuneR)
wavobj <- readWave("guitar.wav")
x <- wavobj@left
fs <- wavobj@samp.rate
# パワースペクトルの計算
y <- fft(x)
y.tmp <- Mod(y) # 複素数の絶対値を計算
y.ampspec <- y.tmp[1:(length(y)/2+1)]
y.ampspec[2:(length(y)/2)] <- y.ampspec[2:(length(y)/2)] * 2
y.powspec <- y.ampspec ^ 2
# 対応する周波数の計算
f <- seq(from=0, to=fs/2, length=length(y)/2+1)
# DCとナイキストを除去
y.powspec <- y.powspec[ 2 : (length(y.powspec)-1) ]
f <- f[ 2 : (length(f)-1) ]
# スペクトル中心の計算（周波数を対数でとるためにlogと^を使用）
(speccent <- 2^(sum(log2(f) * y.powspec) / sum(y.powspec)))
## [1] 429.8408
というように、この音のスペクトル中心が429.8 Hzだということが計算できました。
なお音の明るさの音響指標値としては、Grey & Gordon (1978)ではbarkとsoneを用いてより聴感に則した値の計算を試みていますし、von Bismarck (1974)によるSharpnessもあります。