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です。
コメント
実はですね、ショートカットしようとしてライブラリを導入しようとしたんですよ。しかしながらldがundefined Symbol吐きまくって調べるのやになっちゃったんで、ええい書いてくれるわ、ってことで書いたら1時間で実装できてしまいましたとさ。
ご自由にお使いください。無保証ですが(笑)
＞調べるのやになっちゃったんで、
あるある(苦笑)
＞ご自由にお使いください。無保証ですが(笑)
どうもありがとうございます．いつもより気力多目(当社比20%増しくらい)なので拡張ライブラリも書いてみようかな．
コメントする
- 匿名でのコメントは受け付けておりません。
- お名前（ハンドル名可）とメールアドレスは必ず入力してください。
- メールアドレスを表示されたくないときはURLも必ず記入してください。
- コメント欄でHTMLタグは使用できません。
- コメント本文に日本語(全角文字)がある程度多く含まれている必要があります。
- コメント欄内のURLと思われる文字列は自動的にリンクに変換されます。
- 投稿ボタンを押してエラーがでなければ、投稿は成功しています。反映されるまでには少し時間がかかります。
お疲れ様．RubyでFFTしようとするとGSLを利用するとか各種ライブラリを利用するとか結構メンドいんですよね．
Posted by: あおき : July 14, 2005 10:59 PMかく言うわたしもハマりつつ丁度FFTのコードを書こうとしてたところだったり．もらってこうかな(をゐ)