Created
April 4, 2011 23:46
-
-
Save mro/902724 to your computer and use it in GitHub Desktop.
Discrete Fast Fourier Transform from 'Meyberg, Vachenauer: Hoehere Mathematik II, Springer Berlin, 1991, page 332'
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/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 | |
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 | |
raise ArgumentError.new "Couldn't convert to complex number: '#{self.strip}'" 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| | |
begin | |
values << l.to_c | |
rescue ArgumentError => e | |
$stderr.puts e | |
end | |
end | |
INVERSE = ARGV[0] == 'inverse' | |
$stderr.puts INVERSE ? 'inverse' : 'forward' | |
values.fft!(INVERSE).each {|i| puts "#{i.real} #{i.image}i"} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment