RubyでFFT

レポートでフーリエ変換をどうたらしなければならないので、FFTをRubyで実装してみることにしました。といってもRubyという時点で金持ちプログラミング決定なので効率無視です。

Cooley-Tukey 型 FFTで実装してみました。これよりも効率のよいPrime Factor 型 FFTというのもあるのですが(詳しくはFFTの概略と設計法参照のこと)、因数分解とかがぱっと実装思いつかなかったので、簡単なCooley-Tukey 型にしました。

コードは続きをどうぞ。改善の余地がかなりあると思います。

FFT.rb

require 'complex'

# 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

ドライブするには

data = []
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です。

July 14, 2005 19:30 fenrir が投稿 : 固定リンク | | このエントリーを含むはてなブックマーク

コメント

お疲れ様.RubyでFFTしようとするとGSLを利用するとか各種ライブラリを利用するとか結構メンドいんですよね.
かく言うわたしもハマりつつ丁度FFTのコードを書こうとしてたところだったり.もらってこうかな(をゐ)

Posted by: あおき : July 14, 2005 10:59 PM

実はですね、ショートカットしようとしてライブラリを導入しようとしたんですよ。しかしながらldがundefined Symbol吐きまくって調べるのやになっちゃったんで、ええい書いてくれるわ、ってことで書いたら1時間で実装できてしまいましたとさ。
ご自由にお使いください。無保証ですが(笑)

Posted by: fenrir : July 14, 2005 11:20 PM

>調べるのやになっちゃったんで、
あるある(苦笑)

>ご自由にお使いください。無保証ですが(笑)
どうもありがとうございます.いつもより気力多目(当社比20%増しくらい)なので拡張ライブラリも書いてみようかな.

Posted by: あおき : July 15, 2005 05:01 PM

コメントする