Discrete Fourier Transform

这次的数字图像处理作业是用C/C++ 实现离散Fourier 变换和FFT。我偷了个大懒，花了两小时用Ruby 实现了一个。

首先是O(n²)的离散Fourier 变换。其实都很简单。

require 'complex' a = [1, 0, -1, 0] def fourier_transform(ary) n = ary.size a = Array.new(n) ang = -2*Math::PI/n unit = Complex(Math.cos(ang), Math.sin(ang)) (0 ... n).each do |i| a[i] = 0 (0 ... n).each do |j| a[i] += ary[j] * unit**(i*j) end end a end puts fourier_transform(a) # => 0i # => 2.0+1.22460635382238e-016i # => -2.44921270764475e-016i # => 2.0+3.67381906146713e-016i

然后找了本台湾教授李家同写的《算法设计与分析导论》，伪码写得很详细。

require 'complex' a = [1, 0, -1, 0] def fast_fourier_transform(ary) n = ary.size a = Array.new(n) ang = 2*Math::PI/n unit = Complex(Math.cos(ang), Math.sin(ang)) if n == 2 a[0] = ary[0] + ary[1] a[1] = ary[0] - ary[1] else even, odd = [], [] ary.each_with_index do |x,i| if i %2 == 0 then even << x else odd << x end end b = fast_fourier_transform even c = fast_fourier_transform odd (0 ... n/2).each do |j| a[j] = b[j] + unit**j * c[j] a[j + n/2] = b[j] + unit**(j + n/2) * c[j] end end a end puts fast_fourier_transform(a) # => 0i # => 2.0+0.0i # => 0.0i # => 2.0+0.0i

接着用benchmark 比较一下性能：

require 'benchmark' a = [1, 0, -1, 0] * 2**4 Benchmark.bm(5) do |timer| timer.report('DFT') {fourier_transform(a)} timer.report('FFT') {fast_fourier_transform(a)} end # => user system total real # =>DFT 1.281000 0.000000 1.281000 ( 1.281000) # =>FFT 0.032000 0.000000 0.032000 ( 0.032000)

这是我run 了近一个小时才得到的珍贵数据：

还有一个小问题，这个程序只能处理2 的幂次个元素的序列。如果你传入的参数不是2 的幂次个，就会得到一堆稀奇古怪的错误和Ruby 不适合数值计算的结论。

P.S.:

1. 李家同是个文采学问都做得很好的人，极力推荐他的一篇散文生活思索: 樂盡所能 祇取所需
2. Ruby最吸引我的一点就是它真正让我专注于我需要的东西（在这里是算法。如果你看过伪码就会发现Ruby 代码真的很“干净”）。