高速フーリエ変換の実装

2009-05-16
#ruby #math

これなら分かる 応用数学教室4.7 「高速フーリエ変換」

 １の原始N乗根による表現による離散フーリエ変換に引き続き。

$$ f(x) = a_0 + a_1 x + a_2 x^2 + … + a_{N-1} x^{N-1} = \sum_{l=0}^{N-1} a_l x^l $$

とすると、

$$ b_k = \sum_{l=0}^{N-1} a_l \omega_N^{kl} = f(\omega_N^k) $$

となる。

\(f(x)\)を書き直すと

$$ \begin{align*} f(x) &= (a_0 + a_2 x^2 + … + a_{N-2} x^{N-2}) + x (a_1 + a_3 x^2 + … + a_{N-1} x^{N-2}) \\ &= p(x^2) + xq(x^2) \\ \end{align*} $$

ただし、

$$ \begin{align*} p(x) &= (a_0 + a_2 x + … + a_{N-2} x^{N/2-1}) \\ q(x) &= (a_1 + a_3 x + … + a_{N-1} x^{N/2-1}) \\ \end{align*} $$

とする。

\(\omega_N = e^{i 2\pi/N}\) から \(\omega_N^{2k} = \omega_{N/2}^k\)、\(\omega_N^{N/2+k} = -\omega_N^k\) を使うと

$$ b_k = \Bigg\lbrace \begin{align*} f(\omega_N^j) &= p(\omega_{N/2}^j) + \omega_N^j q(\omega_{N/2}^j), & \\ f(\omega_N^{N/2+j}) &= p(\omega_{N/2}^j) - \omega_N^j q(\omega_{N/2}^j) & \\ & (j = 0, 1, …, N/2-1) \end{align*} $$

と計算することができる。\(p\)と\(q\)は\(f\)と同じ形なので、この分割を再帰的に行って計算量を減らすことができる。

プログラム化

require 'complex'

# 原始N乗根
def primitiveNthRoot(n, k=1)
  t = 2 * Math::PI * k / n
  Complex(Math.cos(t), Math.sin(t))
end

# 奇数？
def odd(x)
  (x & 1) != 0
end

# バタフライ
def f(a)
  n = a.length
  if n == 1
    a
  elsif odd(n)
    raise "not even"
  else
    n2 = n / 2
    ae = (0...n2).map {|l| a[l*2]}  # 偶数成分
    ao = (0...n2).map {|l| a[l*2 + 1]}  # 奇数成分

    p = f(ae)
    q = f(ao)

    b = Array.new(n)
    (0...n2).each do |k|
      wkN = primitiveNthRoot(n, k)
      b[k]      = p[k] + wkN * q[k]
      b[k + n2] = p[k] - wkN * q[k]
    end
    b
  end
end

# 高速フーリエ変換
def fft(fs)
  n = fs.length
  b = f(fs)
  b.map {|be| be.conjugate / n}
end

テスト：

if $0 == __FILE__
  v = nil
  if ARGV.empty?
    # 固定の信号
    v = [1,7,3,2,0,5,0,8]
  else
    # 指定の長さの信号をランダムで生成
    n = ARGV.shift.to_i
    v = (0...n).map {|i| rand * 2 - 1}
  end

  # 高速フーリエ変換
  r = fft(v)
  p r

  #=== 検算

  # 逆フーリエ変換
  def ifft(cs)
    b = f(cs)
    b.map {|c| c.real}
  end
  s = ifft(r)

  # 元の値との誤差
  p v.vsub(s).vlength
end

うまく説明できないので理解したい方は本を参照してください（＾＾；；

関連記事

新着記事