高速フーリエ変換

FFT

離散フーリエ変換 はO(n2)O(n^2)かかってしまって遅いのでO(nlogn)O(nlogn)で解けるアルゴリズム用いてフーリエ変換しよう

いくつか種類がある
クーリー–テューキー
任意の基数で分割されたグループで DFT を解いてそれらを合成することで最終結果を得るという分割統治法的な考え方
なのでデータ長は基数の累乗である必要がある、基数2なら2の累乗
Rader's FFT Algorithm
データ長が素数の場合に有効

基数2のクーリー–テューキー型 FFT をやるぞ

まずこれは DFT です、 WNW_N を回転因子呼んだりする
F(k)=∑n=0N−1f(n)WNkn,WN=e−2πi1NF(k) = \sum_{n=0}^{N-1}f(n)W_N^{kn}, W_N = e^{-2\pi i\frac{1}{N}}

基数2で分けるということは偶数と奇数で分けることになるのでこうなる
F(k)=∑n=0N/2−1f(2n)WNk(2n)+∑n=0N/2−1f(2n+1)WNk(2n+1)F(k) = \sum_{n=0}^{N/2-1}f(2n)W_N^{k(2n)} + \sum_{n=0}^{N/2-1}f(2n+1)W_N^{k(2n+1)}

さらに変形してこうなる
F(k)=∑n=0N/2−1f(2n)WN2kn+WNk∑n=0N/2−1f(2n+1)WN2knF(k) = \sum_{n=0}^{N/2-1}f(2n)W_N^{2kn} + W_N^k\sum_{n=0}^{N/2-1}f(2n+1)W_N^{2kn}

回転因子をこんな感じで変形すると
WN2kn=e−2πi1N2kn=e−4πikn1N=e−2πikn1N/2=WN/2knW_N^{2kn} = e^{-2\pi i\frac{1}{N}2kn} = e^{-4\pi ikn\frac{1}{N}} = e^{-2\pi ikn\frac{1}{N/2}}=W_{N/2}^{kn}

元のデータ長の半分の偶数、奇数 dft で求まることが分かる
F(k)=∑n=0N/2−1f(2n)WN/2kn+WNk∑n=0N/2−1f(2n+1)WN/2knF(k) = \sum_{n=0}^{N/2-1}f(2n)W_{N/2}^{kn} + W_N^k\sum_{n=0}^{N/2-1}f(2n+1)W_{N/2}^{kn}

Ek=∑n=0N/2−1f(2n)WN/2kn,Ok=∑n=0N/2−1f(2n+1)WN/2knE_k = \sum_{n=0}^{N/2-1}f(2n)W_{N/2}^{kn} , O_k = \sum_{n=0}^{N/2-1}f(2n+1)W_{N/2}^{kn} とすると

F(k)=Ek+WNkOkF(k) = E_k + W_N^k O_kとなり、 Ek,OkE_k, O_kは長さ N/2N/2 の DFT なので周期もN/2N/2なので
Ek=Ek+N/2,Ok=Ok+N/2E_k = E_{k+N/2}, O_k = O_{k+N/2} で同じだが、

回転因子は
WNk+N/2=WNkWNN/2=WNke−πiW_N^{k+N/2} = W_N^k W_N^{N/2} = W_N^k e^{-\pi i}
WNN/2=e−2πi1NN2=e−2πi12=e−πiW_N^{N/2} = e^{-2\pi i\frac{1}{N}\frac{N}{2}} = e^{-2\pi i\frac{1}{2}} = e^{-\pi i}

WNke−πiW_N^k e^{-\pi i} はWNkW_N^kくるっと半周だけ回すので対称性により
WNk+N/2=−WNkW_N^{k+N/2} = -W_N^k

なので
F(k)=Ek+WNkOkF(k) = E_k + W_N^k O_k
F(k+N/2)=Ek−WNkOkF(k+N/2) = E_k - W_N^k O_k
k=(0,1...,N/2−1)k = (0,1...,N/2-1)
ということが言える（バタフライ演算と言うらしい）

つまり
データ長NNのDFTの前半kk番目はそのデータの偶数だけ集めたDFTのkk番目とそのデータの奇数だけを集めたDFTのkk番目をWNkW_N^kだけ回したものを足したものと一致し、
データ長NNのDFTの後半kk番目はそのデータの偶数だけ集めたDFTのkk番目とそのデータの奇数だけを集めたDFTのkk番目をWNkW_N^kだけ回したものを引いたものと一致

この分割の操作をデータ長が1になるまで分割して出来たDFTを合成していく感じ
データ長8ならまずx[0], x[2], x[4], x[6]とx[1], x[3], x[5], x[7]に分ける
それをさらに分けてx[0], x[4] , x[2], x[6] , x[1], x[5] , x[3] x[7]
それをさらに分けてx[0] , x[4] , x[2] , x[6] , x[1] , x[5] , x[3] , x[7]
この分割の操作はビット反転で上手いこと出来るらしいです

それを愚直にコードで解くとこうなる
https://www.mycompiler.io/view/97kyyi4Fg3F

それを for とかで関数にパッケージしたものがこれ (Gemini 製)
https://www.mycompiler.io/view/18ZhASeC4Bw

逆FFT
DFTにバタフライ演算を逆回転でして出来たデータ配列の要素をデータ長で割るだけ
逆DFTと一緒でゲースー

Gemini 曰く基数4のクーリー–テューキーが実践向けでよく実装されてるらしいがほんとかな？
まあ基数4のバタフライ演算の導出はいつかやりますよ ToDo

Backlinks

離散コサイン変換