Rで周波数解析

RでFFTなどしてみました。

CRANからライブラリをインストール

wavファイルの読み書きなどを行う、soundライブラリをインストールします。

install.package("sound")

使い方は help で。関数の一覧は library(help="sound") 取得できます。

wavファイルを読み込んでFFT

soundライブラリを使用してwavファイルを読み込みます。loadSampleで読んで、soundで実数配列を得ます。

wav = loadSample("hoge.wav")
w = sound(wav)

プロットするとこんな感じ。

plot(0:(length(w)-1),w,type="l")

FFTは組み込み関数で存在している(fft)ものをそのまま使用できます。

f = fft(w)

パワースペクトルをプロット。

p = (Re(f) ** 2 + Im(f) ** 2) # パワー
x = 0:(length(p)-1)
x = x / (length(x) / rate(wav)) # rateはsoundライブラリの関数
xmax = max(x) / 2
plot(x, p, type="l", xlim=c(0, xmax), log="y")

フィルタとか

biOpsという画像処理のためのライブラリにハイパスフィルタやローパスフィルタがあったので使えるかなーと思ったのですが、画像が対象なのでさすがにそううまくはいかないようです。

使用するには install.package("biOps") でライブラリをインストールしておく必要があります。

wav = loadSample(file)
w = sound(wav)
f = fft(w)
t = imgFFTHighPass(f, 1000)

p = (Re(f) ** 2 + Im(f) ** 2)
p2 = (Re(t) ** 2 + Im(t) ** 2)
x = 0:(length(p)-1)
x = x / (length(x) / rate(wav))

xmax = max(x) / 2
par(bg="grey95", mar=c(5,5,2,1), mfrow=c(2,1))
plot(x, p, type="l", xlim=c(0, xmax), log="y")
par(mar=c(5,5,2,1))
plot(x, p2, type="l", xlim=c(0, xmax), log="y", col=2)