Ruby で FFT (高速フーリエ変換) を書いてみた

Ruby

ref: 【ニコニコ動画】ミクをPCの再生音に合わせて自動で踊らせてみた
↑に触発されて波の処理をしたくなったので、Ruby で FFT (高速フーリエ変換) を書いてみました。
FFT とは、波の形を見て周波数とかを見抜く魔法のことです。数式とか考えたくないので、とにかく Ruby で書いてみました *1 。

def fft(a)
  n = a.size
  return a if n == 1
  w = Complex.polar(1, -2 * Math::PI / n)
  a1 = fft((0 .. n / 2 - 1).map {|i| a[i] + a[i + n / 2] })
  a2 = fft((0 .. n / 2 - 1).map {|i| (a[i] - a[i + n / 2]) * (w ** i) })
  a1.zip(a2).flatten
end

これだけです。短いですよね。
実用上はここから in-place 化とか scramble とか定数オーダの最適化が待っているらしいですが、Ruby で書いてる時点でそんな方向には興味がありません。*2

そんなことよりちょっと動かしてみましょう。まずは入力用の波形を作る。さっきの fft の実装では、サンプル数が 2 の累乗でないとダメです。

N = 64
a = (0...N).map do |n|
  v = Math.sin(2 * 2 * Math::PI * n / N) * 2
  v + Math.cos(5 * 2 * Math::PI * n / N)
end

2 Hz の sin 波と 5 Hz の cos 波を合成した波形。2 Hz の振幅は 2 倍。図示して確認すると

a.map do |v|
  s = [" "] * 20
  min, max = [(-v * 3 + 10).round, 10].sort
  s[min..max] = ["#"] * (max - min)
  s
end.transpose.each do |l|
  puts l.join
end

                                       #
                                      ###
           ##                        #####
          ####                       #####
         ######                     #######
 ###    #######                     #######
################                   #########
################                   #########
################                  ###########                  #
                 #################                  ###########
                 ################                    #########
                 ################                    #########
                  #######    ###                      #######
                  ######                              #######
                   ####                                #####
                    ##                                 #####
                                                        ###
                                                         #

なんかよくわかりませんが、できてるのかな？ これを見て、2 Hz はなんとなくわかりますが、5 Hz が混ざってることはちょっとわからないですよね。さっきの FFT にかけてみます。

fft(a)[0, N/2].each_with_index do |v, n|
  next if n == 0
  puts "%2d Hz: %.3f" % [n, v.abs / (N / 2)]
end

 1 Hz: 0.000
 2 Hz: 2.000
 3 Hz: 0.000
 4 Hz: 0.000
 5 Hz: 1.000
 6 Hz: 0.000
 7 Hz: 0.000
 8 Hz: 0.000
 9 Hz: 0.000
10 Hz: 0.000
11 Hz: 0.000
12 Hz: 0.000
13 Hz: 0.000
14 Hz: 0.000
15 Hz: 0.000
16 Hz: 0.000
17 Hz: 0.000
18 Hz: 0.000
19 Hz: 0.000
20 Hz: 0.000
21 Hz: 0.000
22 Hz: 0.000
23 Hz: 0.000
24 Hz: 0.000
25 Hz: 0.000
26 Hz: 0.000
27 Hz: 0.000
28 Hz: 0.000
29 Hz: 0.000
30 Hz: 0.000
31 Hz: 0.000

ばーん。2 Hz と 5 Hz を見事見破ってくれました。fft の返り値の見かたのポイントは、

  • 意味のある情報は配列の前半だけ
  • 0 番目の値も無視する
  • 値は複素数になっているので大きさ (v.abs) をとる
  • (N / 2) で割った値がその周波数の強さ
  • ちなみに位相の情報は偏角 (v.angle) として入っているらしい (あんま見てない)

という感じです。*3

*1:「FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法」 の「リスト1.2.1-1. 再帰呼び出しによる FFT」の移植です。

*2:まじめに Ruby で FFT とかしたい人は NArray の FFTW3 とか使うといいんだと思います。ぼくはやったことありません。

*3:フーリエ変換の実際 の説明が分かりやすかったです。