普通の(fastじゃない)離散フーリエ変換

2009-05-15
blog

周期の関数を等間隔で分割してそれぞれのサンプル点を

、サンプル点での値をとすると、離散フーリエ変換の係数は

となる。

複素表示を展開すると

となる。

逆変換

サンプル点を復元する逆変換は、

プログラム化

Rubyの配列にベクトル演算を仕込む:

class Array
def vadd(v2)
self.zip(v2).map do |a, b|
a + b
end
end
def vsub(v2)
self.zip(v2).map do |a, b|
a - b
end
end
def vscale(s)
self.map do |x|
x * s
end
end
def vlength()
ll = self.inject(0) do |r, x|
r + x ** 2
end
Math.sqrt(ll)
end
def dot(v2)
self.zip(v2).inject(0) do |r, x|
r + x[0] * x[1]
end
end
end

離散フーリエ変換を行うクラス:


class DFT
def initialize(n)
@n = n
@base_cos = (0...@n).map {|i| make_cos_tbl(8, i)}
@base_sin = (0...@n).map {|i| make_sin_tbl(8, i)}
end

# 変換
def transform(v)
real_parts = @base_cos.map do |b|
v.dot(b) / @n
end
img_parts = @base_sin.map do |b|
-v.dot(b) / @n
end
return real_parts, img_parts
end

# 逆変換
def itransform(real_parts, img_parts)
rc = _extract(real_parts, @base_cos)
is = _extract(img_parts, @base_sin)
return rc.vsub(is)
end

def _extract(coeffs, base)
coeffs.zip(base).inject([0] * @n) do |r, (c, b)|
r.vadd(b.vscale(c))
end
end
private :_extract
end

sin, cosを作るヘルパー関数

def make_cos_tbl(n, k)
(0...n).map do |l|
Math.cos(2 * Math::PI * k * l / n)
end
end

def make_sin_tbl(n, k)
(0...n).map do |l|
Math.sin(2 * Math::PI * k * l / n)
end
end

テスト:

dft8 = DFT.new(8)

#v = make_cos_tbl(8, 1).vscale(2).vadd(make_sin_tbl(8, 2).vscale(3))
v = [1,7,3,2,0,5,0,8]
p v # 元のベクトル

r, i = dft8.transform(v)
p r # 離散フーリエ変換結果の係数:実部
p i # 離散フーリエ変換結果の係数:虚部

tr = dft8.itransform(r, i)
p tr # 係数から復元した値
p tr.vsub(v).vlength # 誤差の大きさ

結果:

[1, 7, 3, 2, 0, 5, 0, 8]  # 元のベクトル
[3.25, 0.8321067811865474, -0.2500000000000002, -0.5821067811865479, -2.25, -0.5821067811865464, -0.2500000000000022, 0.8321067811865511] # 離散フーリエ変換結果の係数:実部
[-0.0, -0.021446609406726047, -0.25, 0.7285533905932745, -1.3471114790620886e-15, -0.7285533905932751, 0.25, 0.02144660940672849] # 離散フーリエ変換結果の係数:虚部
[1.0000000000000018, 7.000000000000002, 3.0000000000000058, 1.9999999999999987, -6.709254972601658e-15, 5.0, -1.3322676295501878e-15, 8.000000000000005] # 係数から復元した値
1.0798297524896547e-14 # 誤差の大きさ
  • 計算量はO(N^2)

リンク