高速フーリエ変換

高速フーリエ変換（こうそくフーリエへんかん、英: fast Fourier transform, FFT）は、離散フーリエ変換（英: discrete Fourier transform, DFT）を計算機上で高速に計算するアルゴリズムである。高速フーリエ変換の逆変換を逆高速フーリエ変換（英: inverse fast Fourier transform, IFFT）と呼ぶ。 複素関数 f(x) の離散フーリエ変換である複素関数 F(t) は以下で定義される。 F ( t ) = ∑ x = 0 N − 1 f ( x ) exp ⁡ ( − i 2 π t x N ) . {\displaystyle F(t)=\sum _{x=0}^{N-1}f(x)\exp \left(-i{\frac {2\pi tx}{N}}\right).} このとき、{x = 0, 1, 2, ..., N − 1} を標本点と言う。 これを直接計算したときの時間計算量は、ランダウの記号を用いて表現すると O(N2) である。 高速フーリエ変換は、この結果を、次数Nが2の累乗のときに O(N log N) の計算量で得るアルゴリズムである。より一般的には、次数が N = ∏ ni と素因数分解できるとき、O(N∑ni) の計算量となる。次数が 2 の累乗のときが最も高速に計算でき、アルゴリズムも単純になるので、0 詰めで次数を調整することもある。 高速フーリエ変換を使って、畳み込み積分などの計算を高速に求めることができる。これも計算量を O(N2) から O(N log N) まで落とせる。 現在は、初期の手法 をより高速化したアルゴリズムが使用されている。 逆変換は正変換と同じと考えて良いが、指数の符号が逆であり、係数 1/N が掛かる。 高速フーリエ変換のプログラム中、どの符号が逆転するかを一々分岐させると、分岐の判定に時間がかかり、パフォーマンスが落ちる。一方、正変換のプログラムと、逆変換のプログラムを両方用意しておくことも考えられるが、共通部分が多いため、無駄が多くなる。このため、複素共役を使った次のような方法が考えられる。 離散フーリエ変換を F ( t ) = ∑ x = 0 N − 1 f ( x ) exp ⁡ ( − i 2 π t x N ) {\displaystyle F(t)=\sum _{x=0}^{N-1}f(x)\exp \left(-i{\frac {2\pi tx}{N}}\right)} で定義したとき、逆変換は f ( x ) = 1 N ∑ t = 0 N − 1 F ( t ) exp ⁡ ( i 2 π t x N ) = 1 N ∑ t = 0 N − 1 F ( t ) ¯ exp ⁡ ( − i 2 π t x N ) ¯ {\displaystyle f(x)={\frac {1}{N}}\sum _{t=0}^{N-1}F(t)\exp \left(i{\frac {2\pi {tx}}{N}}\right)={\frac {1}{N}}{\overline {\sum _{t=0}^{N-1}{\overline {F(t)}}\exp \left(-i{\frac {2\pi tx}{N}}\right)}}} となる。 このため、F(t) の離散フーリエ逆変換を求めるには、 (1) 複素共役を取り、F(t) を求める、 (2) F(t) の正変換の離散フーリエ変換を高速フーリエ変換で行う、 (3) その結果の複素共役を取り、N で割る とすれば良く、正変換の高速フーリエ変換のプログラムがあれば、逆変換は容易に作ることができる。 クーリー–テューキー型アルゴリズムは、代表的な高速フーリエ変換 (FFT) アルゴリズムである。 分割統治法を使ったアルゴリズムで、N = N1 N2 のサイズの変換を、より小さいサイズである N1, N2 のサイズの変換に分割していくことで高速化を図っている。 最もよく知られたクーリー–テューキー型アルゴリズムは、ステップごとに変換のサイズをサイズ N/2 の2つの変換に分割するので、2 の累乗次数に限定される。しかし、一般的には次数は 2 の累乗にはならないので、素因数が偶数と奇数とで別々のアルゴリズムに分岐する。 伝統的なFFTの処理実装の多くは、再帰的な処理を、系統だった再帰をしないアルゴリズムにより実現している。 クーリー–テューキー型アルゴリズムは変換をより小さい変換に分解していくので、後述のような他の離散フーリエ係数のアルゴリズムと任意に組み合わせることができる。とりわけ、N ≤ 8 あたりまで分解すると、固定次数の高速なアルゴリズムに切り替えることが多い。 離散フーリエ係数は、1の原始 N 乗根の1つ WN = e−2πi/N を使うと、次のように表せる。 F ( t ) = ∑ x = 0 N − 1 f ( x ) W N t x . {\displaystyle F(t)=\sum _{x=0}^{N-1}f(x)W_{N}^{tx}.} 例えば、N = 4 のとき、 F ( t ) = X t {\displaystyle F(t)=X_{t}} 、 f ( k ) = x k {\displaystyle f(k)=x_{k}} とすれば、離散フーリエ係数は行列を用いて表現すると（W = W4 と略記） [ X 0 X 1 X 2 X 3 ] = [ W 0 W 0 W 0 W 0 W 0 W 1 W 2 W 3 W 0 W 2 W 4 W 6 W 0 W 3 W 6 W 9 ] [ x 0 x 1 x 2 x 3 ] {\displaystyle {\begin{bmatrix}X_{0}\\X_{1}\\X_{2}\\X_{3}\end{bmatrix}}={\begin{bmatrix}W^{0}&W^{0}&W^{0}&W^{0}\\W^{0}&W^{1}&W^{2}&W^{3}\\W^{0}&W^{2}&W^{4}&W^{6}\\W^{0}&W^{3}&W^{6}&W^{9}\end{bmatrix}}{\begin{bmatrix}x_{0}\\x_{1}\\x_{2}\\x_{3}\end{bmatrix}}} となる。入力列 xk を添字の偶奇で分けて、以下のように変形する。 [ X 0 X 1 X 2 X 3 ] = [ W 0 W 0 W 0 W 0 W 0 W 2 W 1 W 3 W 0 W 4 W 2 W 6 W 0 W 6 W 3 W 9 ] [ x 0 x 2 x 1 x 3 ] = [ W 0 W 0 W 0 W 0 W 0 W 0 W 0 W 2 W 1 W 0 W 1 W 2 W 0 W 0 W 2 W 0 W 2 W 0 W 0 W 2 W 3 W 0 W 3 W 2 ] [ x 0 x 2 x 1 x 3 ] {\displaystyle {\begin{aligned}{\begin{bmatrix}X_{0}\\X_{1}\\X_{2}\\X_{3}\end{bmatrix}}&={\begin{bmatrix}W^{0}&W^{0}&W^{0}&W^{0}\\W^{0}&W^{2}&W^{1}&W^{3}\\W^{0}&W^{4}&W^{2}&W^{6}\\W^{0}&W^{6}&W^{3}&W^{9}\end{bmatrix}}{\begin{bmatrix}x_{0}\\x_{2}\\x_{1}\\x_{3}\end{bmatrix}}\\&={\begin{bmatrix}W^{0}&W^{0}&W^{0}W^{0}&W^{0}W^{0}\\W^{0}&W^{2}&W^{1}W^{0}&W^{1}W^{2}\\W^{0}&W^{0}&W^{2}W^{0}&W^{2}W^{0}\\W^{0}&W^{2}&W^{3}W^{0}&W^{3}W^{2}\end{bmatrix}}{\begin{bmatrix}x_{0}\\x_{2}\\x_{1}\\x_{3}\end{bmatrix}}\end{aligned}}} ( ∵ W k + N = W k {\displaystyle \because W^{k+N}=W^{k}} ) すると、サイズ 2 のFFTの演算結果を用いて表現でき、サイズの分割ができる。 [ X 0 X 1 X 2 X 3 ] = [ 1 0 W 0 0 0 1 0 W 1 1 0 W 2 0 0 1 0 W 3 ] [ W 2 0 W 2 0 0 0 W 2 0 W 2 1 0 0 0 0 W 2 0 W 2 0 0 0 W 2 0 W 2 1 ] [ x 0 x 2 x 1 x 3 ] {\displaystyle {\begin{bmatrix}X_{0}\\X_{1}\\X_{2}\\X_{3}\end{bmatrix}}={\begin{bmatrix}1&0&W^{0}&0\\0&1&0&W^{1}\\1&0&W^{2}&0\\0&1&0&W^{3}\end{bmatrix}}\,{\begin{bmatrix}W_{2}^{0}&W_{2}^{0}&0&0\\W_{2}^{0}&W_{2}^{1}&0&0\\0&0&W_{2}^{0}&W_{2}^{0}\\0&0&W_{2}^{0}&W_{2}^{1}\end{bmatrix}}\,{\begin{bmatrix}x_{0}\\x_{2}\\x_{1}\\x_{3}\end{bmatrix}}} また、この分割手順を図にすると蝶のような図になることから、バタフライ演算とも呼ばれる。 バタフライ演算は、計算機上ではビット反転で実現される。DSPの中には、このバタフライ演算のプログラムを容易にするため、ビット反転アドレッシングを備えているものがある。 N = PQ とする。N 次離散フーリエ変換を、以下のようにP 次離散フーリエ変換とQ 次離散フーリエ変換に分解する。 N 次離散フーリエ変換： F ( n ) = ∑ k = 0 N − 1 f ( k ) exp ⁡ ( − 2 π i n k N ) {\displaystyle F(n)=\sum _{k=0}^{N-1}f(k)\exp \left(-2\pi i{\frac {nk}{N}}\right)} を、n = 0, 1, ..., N − 1 について計算することを考える。n, k を次のように書き換える。ただし 0 ≤ n ≤ N − 1 また 0 ≤ k ≤ N − 1 である。 n = s Q + r 0 ≤ s ≤ P − 1 , 0 ≤ r ≤ Q − 1 k = q P + p 0 ≤ p ≤ P − 1 , 0 ≤ q ≤ Q − 1 {\displaystyle {\begin{aligned}n&=sQ+r\qquad 0\leq s\leq P-1,\,0\leq r\leq Q-1\\k&=qP+p\qquad 0\leq p\leq P-1,\;0\leq q\leq Q-1\end{aligned}}} すると F ( n ) = F ( s Q + r ) = ∑ p = 0 P − 1 ∑ q = 0 Q − 1 f ( q P + p ) exp ⁡ ( − 2 π i ( q P + p ) ( s Q + r ) P Q ) = ∑ p = 0 P − 1 ∑ q = 0 Q − 1 f ( q P + p ) exp ⁡ ( − 2 π i ( r q Q + s p P + p r P Q ) ) = ∑ p = 0 P − 1 [ exp ⁡ ( − 2 π i ( s p P + p r P Q ) ) ∑ q = 0 Q − 1 f ( q P + p ) exp ⁡ ( − 2 π i r q Q ) ] {\displaystyle {\begin{aligned}F(n)&=F(sQ+r)=\sum _{p=0}^{P-1}\sum _{q=0}^{Q-1}f(qP+p)\exp \left(-2\pi i{\frac {(qP+p)(sQ+r)}{PQ}}\right)\\&=\sum _{p=0}^{P-1}\sum _{q=0}^{Q-1}f(qP+p)\exp \left(-2\pi i\left({\frac {rq}{Q}}+{\frac {sp}{P}}+{\frac {pr}{PQ}}\right)\right)\\&=\sum _{p=0}^{P-1}\left[\exp \left(-2\pi i\left({\frac {sp}{P}}+{\frac {pr}{PQ}}\right)\right)\sum _{q=0}^{Q-1}f(qP+p)\exp \left(-2\pi i{\frac {rq}{Q}}\right)\right]\end{aligned}}} ここで、 f 1 ( p , r ) = exp ⁡ ( − 2 π i p r P Q ) ∑ q = 0 Q − 1 f ( q P + p ) exp ⁡ ( − 2 π i r q Q ) {\displaystyle f_{1}(p,r)=\exp \left(-2\pi i{\frac {pr}{PQ}}\right)\sum _{q=0}^{Q-1}f(qP+p)\exp \left(-2\pi i{\frac {rq}{Q}}\right)} と置くと、 F ( n ) = F ( s Q + r ) = ∑ p = 0 P − 1 exp ⁡ ( − 2 π i s p P ) f 1 ( p , r ) {\displaystyle F(n)=F(sQ+r)=\sum _{p=0}^{P-1}\exp \left(-2\pi i{\frac {sp}{P}}\right)f_{1}(p,r)} となる。即ち、F(n) = F(sQ + r) の計算は、次の2ステップになる。 ステップ1 p = 0, 1, ..., P − 1 と r = 0, 1, ..., Q − 1 について f 1 ( p , r ) = exp ⁡ ( − 2 π i p r P Q ) ∑ q = 0 Q − 1 f ( q P + p ) exp ⁡ ( − 2 π i r q Q ) {\displaystyle f_{1}(p,r)=\exp \left(-2\pi i{\frac {pr}{PQ}}\right)\sum _{q=0}^{Q-1}f(qP+p)\exp \left(-2\pi i{\frac {rq}{Q}}\right)} を計算する。これは、Q次の離散フーリエ変換 ∑ q = 0 Q − 1 f ( q P + p ) exp ⁡ ( − 2 π i r q Q ) {\displaystyle \sum _{q=0}^{Q-1}f(qP+p)\exp \left(-2\pi i{\frac {rq}{Q}}\right)} の実行と、回転因子 exp(−2πipr/PQ) の掛け算を、全ての p, r の組（PQ = N 通り）に対して行うことと見ることができる。 ステップ2 s = 0, 1, ..., P − 1 と r = 0, 1, ..., Q − 1 について F ( s Q + r ) = ∑ p = 0 P − 1 exp ⁡ ( − 2 π i s p P ) f 1 ( p , r ) {\displaystyle F(sQ+r)=\sum _{p=0}^{P-1}\exp \left(-2\pi i{\frac {sp}{P}}\right)f_{1}(p,r)} を計算する。ここで、右辺は r を固定すれば、P 次の離散フーリエ変換である。 ステップ1、2は、N = PQ 次の離散フーリエ変換を、Q 次の離散フーリエ変換と回転因子の掛け算の実行により、Q 組 (r = 0, 1, ..., Q − 1) の P 次離散フーリエ変換に分解したと見ることができる。 N = Qk (P = Qk − 1) の場合には、上を繰り返せば、Q 次の離散フーリエ変換と回転因子の掛け算を繰り返すことだけで次数を下げることができ、最終的に1次離散フーリエ変換（何もしないことと同じ）にまで下げると、F(t)を求めることができる。特に、Q が2または4の場合は、Q次の離散フーリエ変換は非常に簡単な計算になる。 Q = 2 の場合は、exp(−2πirq/Q) は 1 か −1 なので、Q 次の離散フーリエ変換は符号の逆転と足し算だけで計算できる。 Q = 4 の場合は、exp(−2πirq/Q) は 1, −1, i, −i のいずれかなので、Q 次の離散フーリエ変換の計算は、符号の逆転、実部虚部の交換と足し算だけで計算できる。 このため、2の累乗あるいは4の累乗次の離散フーリエ変換は簡単に計算できる。実務的に用いられるのは、Q = 2 か Q = 4 の場合のみである。なお、Q = 2 かQ = 4 の場合のこの部分のQ次の離散フーリエ変換のことを、バタフライ演算と言う。 また、Q = 2かQ = 4の場合において、計算を終了するまでに何回の「掛け算」が必要かを考える。符号の逆転、実部虚部の交換は「掛け算」として数えなければ、回転因子の掛け算のみが「掛け算」である。N = Qkの次数を1落とすためにN回の「掛け算」が必要であり、次数をkから0に落とすにはそれをk回繰り返す必要があるため、「掛け算」の数は Nk = N logQ N となる。高速フーリエ変換の計算において時間がかかるのは「掛け算」の部分であるため、これが「高速フーリエ変換では計算速度は O(N' log N) になる」ことの根拠になっている。 上記の説明で、 N = Q k ( P = Q k − 1 ) {\displaystyle N=Q^{k}(P=Q^{k-1})} の場合、N = Qk 個のデータ f ( q Q k − 1 + p ) {\displaystyle f(qQ^{k-1}+p)} から、N = Qk 個の計算結果 f 1 ( p , r ) = exp ⁡ ( − 2 π i p r Q k ) ∑ q = 0 Q − 1 f ( q Q k − 1 + p ) exp ⁡ ( − 2 π i r q Q ) {\displaystyle f_{1}(p,r)=\exp \left(-2\pi i{\frac {pr}{Q^{k}}}\right)\sum _{q=0}^{Q-1}f(qQ^{k-1}+p)\exp \left(-2\pi i{\frac {rq}{Q}}\right)} を計算する場合に、メモリの節約のため、0 ≤ q ≤ Q − 1 と 0 ≤ r ≤ Q − 1 を利用し、計算結果 f 1 ( p , r ) {\displaystyle f_{1}(p,r)} を元データ f ( r Q k − 1 + p ) {\displaystyle f(rQ^{k-1}+p)} のあった場所に格納することが多い。これが次の次数 Qk − 1 でも繰り返されるため、 p = q 2 Q k − 2 + p 2 {\displaystyle p=q_{2}Q^{k-2}+p_{2}} とすると、次の次数の計算結果 f 2 ( p 2 , q 2 , q ) {\displaystyle f_{2}(p_{2},q_{2},q)} は f ( q Q k − 1 + q 2 Q k − 2 + p 2 ) {\displaystyle f(qQ^{k-1}+q_{2}Q^{k-2}+p_{2})} のあった場所に格納される。繰り返せば、 t = q 1 Q k − 1 + q 2 Q k − 2 + ⋯ + q k {\displaystyle t=q_{1}Q^{k-1}+q_{2}Q^{k-2}+\cdots +q_{k}} とすると、計算結果 f k ( p k , q k , q k − 1 , … , q 2 , q 1 ) {\displaystyle f_{k}(p_{k},q_{k},q_{k-1},\dots ,q_{2},q_{1})} は f ( q 1 Q k − 1 + q 2 Q k − 2 + ⋯ + q k − 1 Q + p k ) {\displaystyle f(q_{1}Q^{k-1}+q_{2}Q^{k-2}+\cdots +q_{k-1}Q+p_{k})} のあった場所に格納される。 一方、 F ( s Q + r ) = ∑ p = 0 Q k − 1 − 1 exp ⁡ ( − 2 π i s p Q k − 1 ) f 1 ( p , r ) {\displaystyle F(sQ+r)=\sum _{p=0}^{Q^{k-1}-1}\exp \left(-2\pi i{\frac {sp}{Q^{k-1}}}\right)f_{1}(p,r)} を、r を固定し s を変数とした Qk − 1 次離散フーリエ変換と見なして、 s = s 2 Q + r 2 {\displaystyle s=s_{2}Q+r_{2}} とすると、 F ( s 2 Q 2 + r 2 Q + r ) = ∑ p 2 = 0 Q k − 2 − 1 exp ⁡ ( − 2 π i s 2 p 2 Q k − 2 ) f 2 ( p 2 , r 2 , r ) {\displaystyle F(s_{2}Q^{2}+r_{2}Q+r)=\sum _{p_{2}=0}^{Q^{k-2}-1}\exp \left(-2\pi i{\frac {s_{2}p_{2}}{Q^{k-2}}}\right)f_{2}(p_{2},r_{2},r)} となる。繰り替えせば、 F ( s k Q k + r k Q k − 1 + ⋯ + r 2 Q + r 1 ) = ∑ p k = 0 Q k − k − 1 exp ⁡ ( − 2 π i s k p k Q k − k ) f k ( p k , r k , r k − 1 , … , r 2 , r 1 ) {\displaystyle F(s_{k}Q^{k}+r_{k}Q^{k-1}+\cdots +r_{2}Q+r_{1})=\sum _{p_{k}=0}^{Q^{k-k}-1}\exp \left(-2\pi i{\frac {s_{k}p_{k}}{Q^{k-k}}}\right)f_{k}(p_{k},r_{k},r_{k-1},\dots ,r_{2},r_{1})} となるが、左辺について s k Q k + r k Q k − 1 + ⋯ + r 2 Q + r 1 < Q k {\displaystyle s_{k}Q^{k}+r_{k}Q^{k-1}+\cdots +r_{2}Q+r_{1}<Q^{k}} より sk = 0, また右辺について Q k − k − 1 = 0 {\displaystyle Q^{k-k}-1=0} より pk = 0。このため、 F ( r k Q k − 1 + ⋯ + r 2 Q + r 1 ) = f k ( 0 , r k , r k − 1 , … , r 2 , r 1 ) . {\displaystyle F(r_{k}Q^{k-1}+\cdots +r_{2}Q+r_{1})=f_{k}(0,r_{k},r_{k-1},\dots ,r_{2},r_{1}).} これは f ( r 1 Q k − 1 + r 2 Q k − 2 + ⋯ + r k − 1 Q + r k ) {\displaystyle f(r_{1}Q^{k-1}+r_{2}Q^{k-2}+\cdots +r_{k-1}Q+r_{k})} のあった場所に格納されている。 このように、求める解 F ( r k Q k − 1 + ⋯ + r 2 Q + r 1 ) {\displaystyle F(r_{k}Q^{k-1}+\cdots +r_{2}Q+r_{1})} が f ( r 1 Q k − 1 + r 2 Q k − 2 + ⋯ + r k − 1 Q + r k ) {\displaystyle f(r_{1}Q^{k-1}+r_{2}Q^{k-2}+\cdots +r_{k-1}Q+r_{k})} のあった場所に格納されていることを、ビット反転と言う。これは、Q 進法で表示した場合、 r k Q k − 1 + ⋯ + r 2 Q + r 1 {\displaystyle r_{k}Q^{k-1}+\cdots +r_{2}Q+r_{1}} は ( r k r k − 1 … r 2 r 1 ) Q {\displaystyle (r_{k}r_{k-1}\dots r_{2}r_{1})_{Q}} となるのに対し、 r 1 Q k − 1 + r 2 Q k − 2 + ⋯ + r k − 1 + r k {\displaystyle r_{1}Q^{k-1}+r_{2}Q^{k-2}+\cdots +r_{k-1}+r_{k}} は逆から読んだ ( r 1 r 2 … r k − 1 r k ) Q {\displaystyle (r_{1}r_{2}\dots r_{k-1}r_{k})_{Q}} となるためである。 以下は、高速フーリエ変換のプログラムを Q = 4 の場合にMicrosoft Visual Basicの文法を用いて書いた例である。 Const pi As Double = 3.14159265358979 '円周率 Dim Ndeg As Long '4^deg Dim Pdeg As Long '4^(deg-i) Dim CR() As Double '入力実数部 Dim CI() As Double '入力虚数部 Dim FR() As Double '出力実数部 Dim FI() As Double '出力虚数部 deg=5 '任意に設定。5ならN=4^5=1024で計算 Ndeg=4^deg ReDim CR(Ndeg - 1) As Double '入力実数部 ReDim CI(Ndeg - 1) As Double '入力虚数部 ReDim FR(Ndeg - 1) As Double '出力実数部 ReDim FI(Ndeg - 1) As Double '出力虚数部 'ここで、変換される関数の実部をCR(0)からCR(Ndeg-1)に、虚部をCI(0)からCI(Ndeg-1)に入力しておくこと 'フーリエ変換 For i = 1 To deg Pdeg = 4