Ruby: Simple Fast Fourier Transform

by far not as powerful as the Fastest Fourier Transform in the West but maybe sometimes useful for a quick data analysis or de-noising. Reads stdin and writes to stdout.

Algorithm taken from Meyberg, Vachenauer: Höhere Mathematik II and ported to plain ruby myself.

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
#!/usr/bin/env ruby
require 'complex'
 
class Array
  # DFT and inverse.
  # 
  # Algorithm from 
  # 'Meyberg, Vachenauer: Hoehere Mathematik II, Springer Berlin, 1991, page 332'
  #
  # See http://blog.mro.name/2011/04/simple-ruby-fast-fourier-transform/ 
  #
  def fft doinverse = false
    src = self
    # raise ArgumentError.new "Expected array input but was '#{src.class}'" unless src.kind_of? Array
    n = src.length
    nlog2 = Math.log( n ) / Math.log( 2 )
    raise ArgumentError.new "Input array size must be a power of two but was '#{n}'" unless nlog2.floor - nlog2 == 0.0
    n2 = n / 2
    phi = Math::PI / n2
    if doinverse
      phi = -phi
    else
      src.collect!{|c| c /= n.to_f}
    end
 
    # build a sine/cosine table
    wt = Array.new n2
    wt.each_index { |i| wt[i] = Complex.new Math.cos(i * phi), Math.sin(i * phi) }
 
    # bit reordering
    n1 = n - 1
    j = 0
    1.upto(n1) do |i|
      m = n2
      while j >= m
        j -= m
        m /= 2
      end
      j += m
      src[i],src[j] = src[j],src[i] if j > i
    end
 
    # 1d(N) Ueberlagerungen
    mm = 1
    begin
      m = mm
      mm *= 2
      0.upto(m - 1) do |k|
        w = wt[ k * n2 ]
        k.step(n1, mm) do |i|
          j = i + m
          src[j] = src[i] - (temp = w * src[j])
          src[i] += temp
        end
      end
      n2 /= 2
    end while mm != n
    src
  end
end
 
class String
  # parse Complex.new.to_s
  def to_c
    m = @@PATTERN.match self
    return nil if m.nil?
    Complex.new m[1].to_f, m[2].to_f
  end
private
  # float_pat = /(-?\d+(?:\.\d+)?(?:e[+-]?\d+)?)/
  @@PATTERN = /^[ \t\r\n]*(-?\d+(?:\.\d+)?(?:e[+-]?\d+)?)?\s*((?:\s+|[+-])\d+(?:\.\d+)?(?:e[+-]?\d+)?i)?[ \t\r\n]*$/
end
 
values = []
$stdin.each_line do |l|
  c = l.to_c
  if c.nil?
    $stderr.puts "unmatched '#{l}'"
  else
    values << c
  end
end
INVERSE = ARGV[0] == 'inverse'
$stderr.puts INVERSE ? 'inverse' : 'forward'
 
values.fft(INVERSE).each {|i| puts "#{i.real} #{i.image}i"}

WordPress messes up the angle brackets as usual, but there’s a gist for that.

Flattr this!

Post a Comment

Your email is never published nor shared. Required fields are marked *