mk-mode BLOG

このブログは自作の自宅サーバに構築した Debian GNU/Linux で運用しています。
PC・サーバ構築等の話題を中心に公開しております。(クローンサイト: GitHub Pages

ブログ開設日2009-01-05
サーバ連続稼働時間
Reading...
Page View 合計
Reading...
今日
Reading...
昨日
Reading...

Ruby - (離散)フーリエ変換!

[ プログラミング, 数学 ] [ Ruby ]

こんばんは。

前回、「離散フーリエ変換」の C++ での実装に関する記事を紹介しました。

今回は、同じアルゴリズムを Ruby で実装してみました。
実際、ほとんど同じです。

0. 前提条件

  • Linux Mint 14 Nadia (64bit) での作業を想定。
  • Ruby 2.0.0-p195

1. (離散)フーリエ変換について

前回の記事を参照ください。

2. Ruby スクリプト作成

実際に実装してみた。

  • 変換元の周期関数は、 とした。
  • t の範囲は に限定している。
  • 分割数は N の値を変更して対応する。(以下の例では N=100 としている)
  • 「元データ作成」、「元データを離散フーリエ変換」、「離散フーリエ変換されたデータを逆離散フーリエ変換」としている。
discrete_fourier_transform.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
#*********************************************
# 離散フーリエ変換
#   f(t) = 2 * sin(4 * t) + 3 * cos(2 * t)
#          ( 0 <= t < 2 * pi )
#*********************************************/
#
N = 100                # 分割数
CSV_DFT  = "DFT.csv"   # 出力ファイル (DFT)
CSV_IDFT = "IDFT.csv"  # 出力ファイル (IDFT)

# 計算クラス
class DiscreteFourierTransform
  def initialize
    @src_re  = Array.new
    @src_im  = Array.new
    @dft_re  = Array.new
    @dft_im  = Array.new
  end

  # 元データ作成
  def make_source_data
    0.upto(N - 1) 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
  end

  # 離散フーリエ変換
  def execute_DFT
    # 出力ファイルOPEN
    out_file = open(CSV_DFT, "w")

    # ヘッダ出力 ( k, 角周波数, 元データ(実部), 元データ(虚部), DFT(実部), DFT(虚部) )
    out_file.puts "k,f,x_re,x_im,X_re,X_im"

    # 計算・結果出力
    0.upto(N - 1) do |k|
      dft_re = 0.0
      dft_im = 0.0
      0.upto(N - 1) 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
      out_file.printf("%d,%.6f,%.6f,%.6f,%.6f,%.6f\n",
        k, (2 * Math::PI / N) * k, @src_re[k], @src_im[k], dft_re, dft_im)
    end

    # 出力ファイルCLOSE
    out_file.close
  end

  # 逆離散フーリエ変換
  def execute_IDFT
    # 出力ファイルOPEN
    out_file = open(CSV_IDFT, "w")

    # ヘッダ出力 ( k, 角周波数, DFT(実部), DFT(虚部), IDFT(実部), IDFT(虚部) )
    out_file.puts "k,f,X_re,X_im,x_re,x_im"

    # 計算・結果出力
    0.upto(N - 1) do |n|
      idft_re = 0.0
      idft_im = 0.0
      0.upto(N - 1) 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
      out_file.printf("%d,%.6f,%.6f,%.6f,%.6f,%.6f\n",
        n, (2 * Math::PI / N) * n, @dft_re[n], @dft_im[n], idft_re, idft_im)
    end

    # 出力ファイルCLOSE
    out_file.close
  end
end

# メイン処理
begin
  # 計算クラスインスタンス化
  obj_calc = DiscreteFourierTransform.new

  # 元データ作成
  obj_calc.make_source_data

  # 離散フーリエ変換
  obj_calc.execute_DFT

  # 逆離散フーリエ変換
  obj_calc.execute_IDFT
rescue => e
  # エラーメッセージ
  puts "[例外発生] #{e}"
end

離散フーリエ変換と逆離散フーリエ変換の処理は、実質符号を変えるだけなので、1つにまとめても良いかもしれない。

3. 実行

実際に実行して検証してみる。

1
$ ruby discrete_fourier_transform

コンソールには特に何も表示しない。
アプリと同じディレクトリに “DFT.csv”, “IDFT.csv” という CSV ファイルが作成される。
内容は以下のようになっているはず。

DFT.csv
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
k,f,x_re,x_im,X_re,X_im
0,0.000000,3.000000,0.000000,0.000000,0.000000
1,0.062832,3.473724,0.000000,0.000000,-0.000000
2,0.125664,3.869257,0.000000,150.000000,-0.000000
3,0.188496,4.158424,0.000000,0.000000,-0.000000
4,0.251327,4.317576,0.000000,0.000000,-100.000000
5,0.314159,4.329164,0.000000,-0.000000,0.000000
6,0.376991,4.182959,0.000000,0.000000,0.000000
7,0.439823,3.876846,0.000000,0.000000,0.000000
8,0.502655,3.417134,0.000000,0.000000,0.000000
9,0.565487,2.818364,0.000000,0.000000,-0.000000
         :
====< 途中省略 >====
         :
90,5.654867,-0.248520,0.000000,-0.000000,-0.000000
91,5.717699,-0.263689,0.000000,-0.000000,0.000000
92,5.780530,-0.202174,0.000000,-0.000000,-0.000000
93,5.843362,-0.052303,0.000000,-0.000000,0.000000
94,5.906194,0.190852,0.000000,0.000000,0.000000
95,5.969026,0.524938,0.000000,-0.000000,-0.000000
96,6.031858,0.940264,0.000000,0.000000,100.000000
97,6.094690,1.420235,0.000000,-0.000000,-0.000000
98,6.157522,1.942242,0.000000,150.000000,-0.000000
99,6.220353,2.478964,0.000000,0.000000,-0.000000
IDFT.csv
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
k,f,X_re,X_im,x_re,x_im
0,0.000000,0.000000,0.000000,3.000000,-0.000000
1,0.062832,0.000000,-0.000000,3.473724,-0.000000
2,0.125664,150.000000,-0.000000,3.869257,-0.000000
3,0.188496,0.000000,-0.000000,4.158424,-0.000000
4,0.251327,0.000000,-100.000000,4.317576,-0.000000
5,0.314159,-0.000000,0.000000,4.329164,-0.000000
6,0.376991,0.000000,0.000000,4.182959,-0.000000
7,0.439823,0.000000,0.000000,3.876846,-0.000000
8,0.502655,0.000000,0.000000,3.417134,-0.000000
9,0.565487,0.000000,-0.000000,2.818364,-0.000000
         :
====< 途中省略 >====
         :
90,5.654867,-0.000000,-0.000000,-0.248520,0.000000
91,5.717699,-0.000000,0.000000,-0.263689,-0.000000
92,5.780530,-0.000000,-0.000000,-0.202174,-0.000000
93,5.843362,-0.000000,0.000000,-0.052303,0.000000
94,5.906194,0.000000,0.000000,0.190852,-0.000000
95,5.969026,-0.000000,-0.000000,0.524938,-0.000000
96,6.031858,0.000000,100.000000,0.940264,0.000000
97,6.094690,-0.000000,-0.000000,1.420235,0.000000
98,6.157522,150.000000,-0.000000,1.942242,0.000000
99,6.220353,0.000000,-0.000000,2.478964,-0.000000

離散フーリエ変換後のデータを逆離散フーリエ変換して、元のデータに戻っていることが分かる。
また、C++ 版と同じ結果になっていることも分かる。(結果が同じなのでグラフ出力は省略)


電気工学、音響学、振動、光学等でよく使用する重要な概念です。応用範囲は広いので他にも利用できるかと思います。

近いうちに、離散フーリエ変換を高速に行う「高速フーリエ変換」について考えてみたいとも思っています。

以上。

Comments