統計解析のための言語(R)で音をどうこうしようというのがそもそも間違いなのだ。ということで数値解析ソフトウェアであるGNU Octaveを使ってみることに。
インストール
家ではMacを使用しているのでMac版を入れました。High Performance Computing for Mac OS Xからバイナリがダウンロードできるようなので、これをそのまま落としました。Octaveは単体だとグラフの描画などができないので、GnuPlotも一緒にインストールします。こちらはMacPortsからインストールしました。
sudo port install gnuplot
使ってみた
組み込み関数のwavreadでwavファイルを読み込めるようです。wavreadで実数配列を取得してから、Rと同様にfft関数でFFTを実行します。
y = abs(fft(wavread("hoge.wav"))) % パワー
t = linspace(0, 44100, length(y))
plot(t, log10(y))
axis([0, 44100/2])
ここまではRと同じ。次にフィルタなどかけてみたかったのですが、いまいちfilter関数の使い方がわからないので、試しに特定の周波数以下の値を切り捨てるような処理を書いてみました。
src = fft(wavread("hoge.wav"))
y = src
t = linspace(0, 44100, length(y))
y(1:find(t > 5000)(0)) = 0
plot(t, log10(abs(y)))
axis([0, 44100/2])
Octaveのベクトル(配列)は [] ではなく () で要素を参照するようです。このとき、 (x:y) とすると x〜y までの値を参照できます。 y(1:100) = 0 とすると 1〜100 までの要素がすべて 0 になります。R もそうですが Octave の配列のインデックスは 1 始まりなので注意。
Octave (in Japanese)によると Octave の for 文は遅いらしく、実際体感的にもかなり遅かったので、上記のような書き方を覚えておいた方が良さそうです。
find 関数は引数の条件に合致するインデックスをすべて返します。 find(t > 5000) とすると t の値が 5000 より大きくなっているインデックスの一覧が返ります。tはソートされているので find(t > 5000)(0) で5000を超えた最初の位置がわかります。
y(1:find(t > 5000)(0)) = 0 とすると、周波数(t)が5000以下の要素をすべて0にするということになります。
で、今度はこれを逆FFTして時間軸波形に戻します。
dst = real(ifft(y)) plot(dst)
ちなみに元の波形は以下。
wavwrite関数を使うとwavファイルを書き出すことができるようなので、保存してみます。
wavwrite(dst, 44100, "after.wav") % 44100 = サンプリング周波数
再生してみると音が変わっていることが確認できますが、まぁ音が小さくなるくらいですね。ノイズのある音声データからノイズを取り除けるようなことができると面白いのでしょうけれど、疲れたのでこの辺で。
