Ruby: Simple Fast Fourier Transform

 1#!/usr/bin/env ruby
 2require 'complex'
 3
 4class Array
 5  # DFT and inverse.
 6  # 
 7  # Algorithm from 
 8  # 'Meyberg, Vachenauer: Hoehere Mathematik II, Springer Berlin, 1991, page 332'
 9  #
10  # See http://blog.mro.name/2011/04/simple-ruby-fast-fourier-transform/ 
11  #
12  def fft doinverse = false
13    src = self
14    # raise ArgumentError.new "Expected array input but was '#{src.class}'" unless src.kind_of? Array
15    n = src.length
16    nlog2 = Math.log( n ) / Math.log( 2 )
17    raise ArgumentError.new "Input array size must be a power of two but was '#{n}'" unless nlog2.floor - nlog2 == 0.0
18    n2 = n / 2
19    phi = Math::PI / n2
20    if doinverse
21      phi = -phi
22    else
23      src.collect!{|c| c /= n.to_f}
24    end
25
26    # build a sine/cosine table
27    wt = Array.new n2
28    wt.each_index { |i| wt[i] = Complex.new Math.cos(i * phi), Math.sin(i * phi) }
29
30    # bit reordering
31    n1 = n - 1
32    j = 0
33    1.upto(n1) do |i|
34      m = n2
35      while j >= m
36        j -= m
37        m /= 2
38      end
39      j += m
40      src[i],src[j] = src[j],src[i] if j > i
41    end
42
43    # 1d(N) Ueberlagerungen
44    mm = 1
45    begin
46      m = mm
47      mm *= 2
48      0.upto(m - 1) do |k|
49        w = wt[ k * n2 ]
50        k.step(n1, mm) do |i|
51          j = i + m
52          src[j] = src[i] - (temp = w * src[j])
53          src[i] += temp
54        end
55      end
56      n2 /= 2
57    end while mm != n
58    src
59  end
60end
61
62class String
63  # parse Complex.new.to_s
64  def to_c
65    m = @@PATTERN.match self
66    return nil if m.nil?
67    Complex.new m[1].to_f, m[2].to_f
68  end
69private
70  # float_pat = /(-?\d+(?:\.\d+)?(?:e[+-]?\d+)?)/
71  @@PATTERN = /^[ \t\r\n]*(-?\d+(?:\.\d+)?(?:e[+-]?\d+)?)?\s*((?:\s+|[+-])\d+(?:\.\d+)?(?:e[+-]?\d+)?i)?[ \t\r\n]*$/
72end
73
74values = []
75$stdin.each_line do |l|
76  c = l.to_c
77  if c.nil?
78    $stderr.puts "unmatched '#{l}'"
79  else
80    values << c
81  end
82end
83INVERSE = ARGV[0] == 'inverse'
84$stderr.puts INVERSE ? 'inverse' : 'forward'
85
86values.fft(INVERSE).each {|i| puts "#{i.real} #{i.image}i"}