数(ryフーリエ変換2

#前回の~~ラブライブ！~~まとめ

離散フーリエ変換は

F[k] = \sum^{N-1}_{m=0}f[m]e^{-i\frac{2\pi}{N}mk} \\
f[n] = \dfrac{1}{N}\sum^{N-1}_{k=0}F[k]e^{i\frac{2\pi}{N}nk}

のようになります。

[導出] https://qiita.com/DivineJK/items/f42f89beb939cba3a3ee

#今回の目標

高速フーリエ変換を実装するところまでやります。

#離散フーリエ変換を素直に実装

まずは、式にしたがって、離散フーリエ変換のコードを書いてみます。

simple_dft.py

import math

# 標準入力の受け取り
N = int(input())
f = list(map(int, input().split()))

F = [0 for _ in range(N)] # ここに変換後の値を代入

omega = 2 * math.pi / N
grow = 1
fund = complex(math.cos(theta), -math.sin(theta)) # 1のN乗根

for k in range(N):
    for m in range(N):
        F[k] += f[m] * grow
        grow *= fund

print(F) # 出力

$%$これはさらっと見ると2重ループがあるので少なくとも$O(N^2)$です。高速フーリエ変換では、これを改善します。

#高速フーリエ変換

$N$が2の累乗であるとします。つまり、ある正の整数$M$を用いて$N=2^M$とあらわせるとします。もし、与えられたデータの個数が2の累乗でないならば、データ数が2の累乗になるまで0埋めします。$N=1$も2の0乗とみなせますが、その離散フーリエ変換は自分自身になるので扱いません。

$%$まずは、離散フーリエ変換の式に現れる複素指数関数をみてみます。
$m$が偶数、つまり$m=2l(l\in \mathbb{Z})$のとき、

e^{-i\frac{2\pi}{N}mk} = e^{-i\frac{2\pi}{N}2lk} = e^{-i\frac{2\pi}{N/2}lk}

$m$が奇数、つまり$m=2l+1$のとき、

e^{-i\frac{2\pi}{N}mk} = e^{-i\frac{2\pi}{N}(2l+1)k} = e^{-i\frac{2\pi}{N}k}e^{-i\frac{2\pi}{N/2}lk}

どちらも$e^{-i\frac{2\pi}{N/2}lk}$が現れていることがわかると思います。このことから、離散フーリエ変換の式は

\begin{align}
F[k] = \sum^{N/2-1}_{l=0}f[2l]e^{-i\frac{2\pi}{N/2}lk} + e^{-i\frac{2\pi}{N}k}\sum^{N/2-1}_{l=0}f[2l+1]e^{-i\frac{2\pi}{N/2}lk}
\end{align}

$%$という、二つの和に分けることができます。これらはそれぞれデータ数$N/2$の離散フーリエ変換（に、適当な係数がかかったもの）になっています。
次に、$0\le k\le N/2-1$として、$F[k+N/2]$を考えます。これは、

\begin{align}
F[k+N/2] &= \sum^{N/2-1}_{l=0}f[2l]e^{-i\frac{2\pi}{N/2}l(k+N/2)} + e^{-i\frac{2\pi}{N}(k+N/2)}\sum^{N/2-1}_{l=0}f[2l+1]e^{-i\frac{2\pi}{N/2}l(k+N/2)} \\
&= \sum^{N/2-1}_{l=0}f[2l]e^{-i\frac{2\pi}{N/2}lk} - e^{-i\frac{2\pi}{N}k}\sum^{N/2-1}_{l=0}f[2l+1]e^{-i\frac{2\pi}{N/2}lk}
\end{align}

$%$となります。$F[k]$とほぼ同じ形の式が出てきました。この結果は、データの配列$f$の偶数項目と奇数項目を取り出してできた長さ$N/2$の2つの配列について、離散フーリエ変換を求めることを$N/2$回行っていけば良いということになります。

これをもとに擬似コードを書いてみます。

def simply_fft(f, N):
    if N == 1:
        return f # 要素数1の場合、自分自身
    else:
        fe = [f[2*i] for i in range(N//2)]   # fの偶数項目コレクション
        fo = [f[2*i+1] for i in range(N//2)] # fの奇数項目コレクション
        fe = simply_fft(fe, N//2)
        fo = simply_fft(fo, N//2)
        theta = 2 * math.pi / N
        grow = 1
        fund = complex(math.cos(theta), -math.sin(theta)) # 1のN乗根
        FFT = [0 for i in range(N)] # このリストは作らず直接fに代入してもok
        for i in range(N//2):
            FFT[i] = fe[i] + grow * fo[i]
            FFT[i+N//2] = fe[i] - grow * fo[i]
            grow *= fund
        return FFT

このアルゴリズムの実行時間を$T_N$で表すと、

\begin{align}
T_N = \begin{cases}
O(1) & (N = 1) \\
T_{N/2} + O(N) & (\mathrm{otherwise})
\end{cases}
\end{align}

となり、これを解くと$T_N = O(N\log N)$です。

#参考