高速フーリエ変換の実装

2009-05-16
#ruby #math

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

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

$$ \begin{align*} 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 \tag{4.53} \end{align*} $$

とすると、

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

となる。

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

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

ただし、

$$ \begin{cases} 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{cases} \tag{4.57} $$

とする。

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

$$ \begin{align*} \begin{cases} f(\omega_N^k) &= p(\omega_{N/2}^k) + \omega_N^k q(\omega_{N/2}^k), \\ f(\omega_N^{N/2+k}) &= p(\omega_{N/2}^k) - \omega_N^k q(\omega_{N/2}^k) \end{cases} \\ (k = 0, 1, …, N/2-1) \end{align*} \tag{4.61} $$

と計算することができる。\(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

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

関連記事

新着記事