July 14, 2005

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です。

19:30 fenrir が投稿 : 固定リンク | | このエントリーを含むはてなブックマーク | この記事をdel.icio.usでブックマーク | トラックバック
このエントリーのトラックバックURL: http://fenrir.naruoka.org/mt/mt-tb.cgi/408
コメント

お疲れ様.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
コメントする









名前、アドレスを登録しますか?
(次回以降コメント入力が楽になります)
  • 匿名でのコメントは受け付けておりません。
  • 名前(ハンドル名可)とメールアドレスは必ず入力してください。
  • メールアドレスを表示されたくないときはURLも必ず記入してください。
  • コメント欄でHTMLタグは使用できません。
  • コメント本文に日本語(全角文字)がある程度多く含まれている必要があります。
  • コメント欄内のURLと思われる文字列は自動的にリンクに変換されます。