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

2009-05-15
#math #ruby

周期\(２\pi\)の関数を等間隔で\(N\)分割してそれぞれのサンプル点を

$$ \theta_l = \frac{2\pi}{N}, (l = 0, 1, 2 … N-1) $$

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

$$ F_k = \frac{1}{N} \sum_{l=0}^{N-1} f_l e^{-i2\pi kl/N} $$

となる。

複素表示を展開すると

$$ \begin{align*} F_k &= \frac{1}{N} \sum_{l=0}^{N-1} f_l e^{-i2\pi kl/N} \\ &= \frac{1}{N} \sum_{l=0}^{N-1} \left\lbrace f_l \left( \cos(-2\pi kl/N) + i \sin(-2\pi kl/N) \right) \right\rbrace \\ &= \frac{1}{N} \sum_{l=0}^{N-1} \left\lbrace f_l \left( \cos(2\pi kl/N) - i \sin(2\pi kl/N) \right) \right\rbrace \\ &= \frac{1}{N} \sum_{l=0}^{N-1} f_l \cos(2\pi kl/N) - i\frac{1}{N} \sum_{l=0}^{N-1} f_l \sin(2\pi kl/N) \\ \end{align*} $$

となる。

逆変換

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

$$ \begin{align*} f_l &= \sum_{k=0}^{N-1} F_k e^{i2\pi kl/N} \\ &= \sum_{k=0}^{N-1} (Fr_k + i Fi_k) (\cos(2\pi kl/N) + i\sin(2\pi kl/N)) \\ &= \sum_{k=0}^{N-1} (Fr_k \cos(2\pi kl/N) - Fi_k \sin(2\pi kl/N)) + i \sum_{k=0}^{N-1} (Fr_k \sin(2\pi kl/N) + Fi_k \cos(2\pi kl/N)) \\ \end{align*} $$

プログラム化

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)

リンク

関連記事

新着記事