Ruby script to compute discrete Fourier transform.

discrete_fourier_transform.rb

#! /usr/local/bin/ruby
#*********************************************
# 離散フーリエ変換
#   f(t) = 2 * sin(4 * t) + 3 * cos(2 * t)
#          ( 0 <= t < 2 * pi )
#*********************************************/
#
# 計算クラス
class DiscreteFourierTransform
  N = 100                # 分割数
  CSV_DFT  = "DFT.csv"   # 出力ファイル (DFT)
  CSV_IDFT = "IDFT.csv"  # 出力ファイル (IDFT)

  def initialize
    @src_re, @src_im = Array.new, Array.new
    @dft_re, @dft_im = Array.new, Array.new
  end

  # 元データ作成
  def make_source_data
    N.times do |i|
      @src_re << 2 * Math.sin(4 * (2 * Math::PI / N) * i) \
               + 3 * Math.cos(2 * (2 * Math::PI / N) * i)
      @src_im << 0.0
    end
  rescue => e
    raise
  end

  # 離散フーリエ変換
  def execute_DFT
    # 出力ファイルOPEN
    open(CSV_DFT, "w") do |f|
      # ヘッダ出力 ( k, 角周波数, 元データ(実部), 元データ(虚部), DFT(実部), DFT(虚部) )
      f.puts "k,f,x_re,x_im,X_re,X_im"
      # 計算・結果出力
      N.times do |k|
        dft_re = 0.0
        dft_im = 0.0
        N.times do |n|
          dft_re += @src_re[n] * ( Math.cos((2 * Math::PI / N) * k * n))
                  + @src_im[n] * ( Math.sin((2 * Math::PI / N) * k * n))
          dft_im += @src_re[n] * (-Math.sin((2 * Math::PI / N) * k * n))
                  + @src_im[n] * ( Math.cos((2 * Math::PI / N) * k * n))
        end
        @dft_re << dft_re
        @dft_im << dft_im
        f.puts "%d,%.6f,%.6f,%.6f,%.6f,%.6f" %
          [k, (2 * Math::PI / N) * k, @src_re[k], @src_im[k], dft_re, dft_im]
      end
    end
  rescue => e
    raise
  end

  # 逆離散フーリエ変換
  def execute_IDFT
    # 出力ファイルOPEN
    open(CSV_IDFT, "w") do |f|
      # ヘッダ出力 ( k, 角周波数, DFT(実部), DFT(虚部), IDFT(実部), IDFT(虚部) )
      f.puts "k,f,X_re,X_im,x_re,x_im"
      # 計算・結果出力
      N.times do |n|
        idft_re, idft_im  = 0.0, 0.0
        N.times do |k|
          idft_re += @dft_re[k] * (Math.cos((2 * Math::PI / N) * k * n)) \
                   - @dft_im[k] * (Math.sin((2 * Math::PI / N) * k * n))
          idft_im += @dft_re[k] * (Math.sin((2 * Math::PI / N) * k * n)) \
                   + @dft_im[k] * (Math.cos((2 * Math::PI / N) * k * n))
        end
        idft_re /= N
        idft_im /= N
        f.puts "%d,%.6f,%.6f,%.6f,%.6f,%.6f" %
          [n, (2 * Math::PI / N) * n, @dft_re[n], @dft_im[n], idft_re, idft_im]
      end
    end
  rescue => e
    raise
  end
end

if __FILE__ == $0
  begin
    obj = DiscreteFourierTransform.new  # 計算クラスインスタンス化
    obj.make_source_data                # 元データ作成
    obj.execute_DFT                     # 離散フーリエ変換
    obj.execute_IDFT                    # 逆離散フーリエ変換
  rescue => e
    $stderr.puts "[#{e.class}] #{e.message}"
    e.backtrace.each{ |tr| $stderr.puts "\t#{tr}" }
  end
end