July 14, 2005
RubyでFFT
レポートでフーリエ変換をどうたらしなければならないので、FFTをRubyで実装してみることにしました。といってもRubyという時点で金持ちプログラミング決定なので効率無視です。
Cooley-Tukey 型 FFTで実装してみました。これよりも効率のよいPrime Factor 型 FFTというのもあるのですが(詳しくはFFTの概略と設計法参照のこと)、因数分解とかがぱっと実装思いつかなかったので、簡単なCooley-Tukey 型にしました。
コードは続きをどうぞ。改善の余地がかなりあると思います。
FFT.rb
# FFT
class FFT
def FFT.fft(data)
#まず何回2で割り切れるか数える
power2 = 0
size = data.size
while (size & 0b1) == 0
power2 += 1
size >>= 1
end
# FFTをする前に並べ替え
border = data.size
sorted = Array::new([data])
(0...power2).each{|i|
border >>= 1
temp_plus = []
temp_minus = []
sorted.each{|element|
sorted_plus = []
sorted_minus = []
(0...border).each{|j|
sorted_plus << (element[j] + element[j + border])
theta = 2 * j * Math::PI / element.size
sorted_minus << (element[j] - element[j + border]) * Complex(Math::cos(theta), -Math::sin(theta))
}
temp_plus << sorted_plus
temp_minus << sorted_minus
}
sorted = temp_plus + temp_minus
#p sorted
}
# FFTする
result = Array::new()
sorted.each_with_index{|element, i|
(0...element.size).each{|j|
k = j * (2**power2)
result[k + i] = ft(element, k)
#p k + i
}
}
return result
end
# X_{k} = \sum_{n=0}^{N-1} x_{n} \exp^{\frac{2 \pi n k}{N}}
def FFT.ft(data, k)
result = 0
data.each_with_index{|v, n|
theta = 2 * Math::PI * n * k / data.size
result += v * Complex::new(Math::cos(theta), -Math::sin(theta))
}
return result
end
end
ドライブするには
open(DATA_FILE){|src|
src.each{|line|
data << line.chop.to_f
}
}
FFT::fft(data).each{|c|
puts "#{c.real}\t#{c.image}\t#{c.abs2}"
}
とかすればOKです。
お疲れ様．RubyでFFTしようとするとGSLを利用するとか各種ライブラリを利用するとか結構メンドいんですよね．
かく言うわたしもハマりつつ丁度FFTのコードを書こうとしてたところだったり．もらってこうかな(をゐ)
実はですね、ショートカットしようとしてライブラリを導入しようとしたんですよ。しかしながらldがundefined Symbol吐きまくって調べるのやになっちゃったんで、ええい書いてくれるわ、ってことで書いたら1時間で実装できてしまいましたとさ。
ご自由にお使いください。無保証ですが(笑)
＞調べるのやになっちゃったんで、
あるある(苦笑)
＞ご自由にお使いください。無保証ですが(笑)
どうもありがとうございます．いつもより気力多目(当社比20%増しくらい)なので拡張ライブラリも書いてみようかな．