高速フーリエ変換

出典: フリー百科事典『ウィキペディア（Wikipedia）』

移動：案内, 検索

高速フーリエ変換（こうそくフーリエへんかん、Fast Fourier Transform, FFT）とは、離散フーリエ変換 (Discrete Fourier Transform, DFT) を計算機上で高速に計算するアルゴリズム。逆変換をIFFT (Inverse FFT) という。

高速フーリエ変換といえば一般的には1965年、ジェイムズ・クーリー (J. W. Cooley) とジョン・テューキー (J. W. Tukey) が発見した^[1]とされているCooley-Tukey型FFTアルゴリズムを呼ぶ。しかし、1805年ごろにガウスが同様のアルゴリズムを独自に発見していた^[2]。

[編集] 概要

離散フーリエ変換 (DFT) は以下の級数で定義される。

$f_j = \sum_{k=0}^{N-1} x_k \exp\left( - \frac{2 \pi i j k}{N} \right) \qquad j = 0,\dots,N-1$ .

これを直接計算したときの時間計算量は O ( N² ) である（Oはランダウの記号）。

高速フーリエ変換は、この結果を、次数 N が2の累乗のときに O ( N log N ) の計算量で得るアルゴリズムである。より一般的には、次数が $N = \prod n_i$ と素因数分解できるとき、 $O ( N \sum n_i )$ の計算量となる。次数が2の累乗のときが最も高速に計算でき、アルゴリズムも単純になるので、0詰めで次数を調整することもある。

逆変換 (IFFT) は正変換 (FFT) と同じと考えて良いが、指数の符号が逆であり、係数 1 / N が掛かる。高速フーリエ変換のアルゴリズムはこれらの変更点を簡単に適用できるため、逆変換も同じアルゴリズムを使って計算できる。

高速フーリエ変換を使って、畳み込み積分などの計算を高速に求めることができる。これも計算量をO ( N² ) → O ( N log N ) に落とせる。

現在は、初期の手法^[1]をより高速化したアルゴリズムが使用されている。

[編集] Cooley-Tukey型FFTアルゴリズム

Cooley-Tukey型アルゴリズムは、代表的な高速フーリエ変換 (FFT) アルゴリズムである。

分割統治法 (divide and conquer) を使ったアルゴリズムで、再帰的に N = N₁ N₂ のサイズの変換を、ガウス平面における回転因子である1の冪根を N 程度のオーダーかけて、より小さいサイズである N₁, N₂ のサイズの変換にしていくことで高速化を図っている。

最もよく知られたCooley-Tukey型アルゴリズムは、ステップごとに変換のサイズをサイズ N / 2 の2つの変換に分割するので、2の累乗次数に限定される。しかし、一般的には次数は2の累乗にはならないので、素因数が偶数と奇数とで別々のアルゴリズムに分岐する。

伝統的なFFTの処理実装の多くは、再帰的な処理を、系統だった再帰をしないアルゴリズムにより実現している。

Cooley-Tukey型アルゴリズムは変換をより小さい変換に分解していくので、後述のような他の離散フーリエ係数のアルゴリズムと任意に組み合わせることができる。とりわけ、N ≤ 8 あたりまで分解すると、固定次数の高速なアルゴリズムに切り替えることが多い。

[編集] 原理の簡単な説明

FFTのバタフライ演算

離散フーリエ係数は、1の原始 N 乗根の1つ W_N ≡ exp( -2 π i / N ) を使うと、次のように表せる。

$f_j = \sum_{k=0}^{N-1} x_k W_N^{jk} \qquad j = 0,\dots,N-1$ .

例えば、N = 4 のとき、離散フーリエ係数は行列を用いて表現すると（W ≡ W₄ と略記）

$\begin{bmatrix}f_0\\f_1\\f_2\\f_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}$ .

入力列 x_k を添字の偶奇で分けて、以下のように変形する。

$\begin{matrix} \begin{bmatrix}f_0\\f_1\\f_2\\f_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^0W^0&W^0W^0\\ W^0&W^2&W^1W^0&W^1W^2\\ W^0&W^0&W^2W^0&W^2W^0\\ W^0&W^2&W^3W^0&W^3W^2 \end{bmatrix} \begin{bmatrix}x_0\\x_2\\x_1\\x_3\end{bmatrix} \end{matrix}.$

すると、サイズ2のFFTの演算結果を用いて表現でき、サイズの分割ができる。

$\begin{bmatrix}f_0\\f_1\\f_2\\f_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次離散フーリエ変換 $F(t)=\sum_{j=0}^{N-1}f(j)\exp(-2\pi ijt/N)$ を考える。

$0 = < j < = N - 1$ なら、 $j = q P + p$ ただし $0 = < p < = P - 1$ 、 $0 = < q < = Q - 1$ と書ける。

$0 = < t < = N - 1$ なら、 $t = s Q + r$ ただし $0 = < s < = P - 1$ 、 $0 = < r < = Q - 1$ と書ける。

$F(sQ+r)=\sum_{p=0}^{P-1}\sum_{q=0}^{Q-1}f(qP+p)\exp(-2\pi i(qP+p)(sQ+r)/PQ)$ $=\sum_{p=0}^{P-1}\sum_{q=0}^{Q-1}f(qP+p)\exp(-2\pi i(rq/Q+sp/P+pr/PQ))$ $=\sum_{p=0}^{P-1}\exp(-2\pi i(sp/P+pr/PQ))\sum_{q=0}^{Q-1}f(qP+p)\exp(-2\pi irq/Q)$

ここで、 $f_1(p,r)= \exp(-2\pi ipr/PQ)\sum_{q=0}^{Q-1}f(qP+p)\exp(-2\pi irq/Q)$ と置く。

$f 1 (p, r)$ はsの影響を受けない。この具体的計算は、Q次の離散フーリエ変換 $\sum_{q=0}^{Q-1}f(qP+p)\exp(-2\pi irq/Q)$ の実行と、回転因子 $exp( - 2π i p r / P Q)$ の掛け算を全てのp、rの組(PQ通り=N通り)に対して行うことによって行われる。

$F(sQ+r)=\sum_{p=0}^{P-1}\exp(-2\pi isp/P)f_1(p,r)$ ここで、右辺は、rを固定すれば、P次の離散フーリエ変換であるので、 N次＝PQ次の離散フーリエ変換を、 Q次の離散フーリエ変換と回転因子の掛け算の実行により、 Q組（r=0,1,2,,,Q-1）のP次離散フーリエ変換に分解したと見ることができる。

$N = Q k$ の場合（ $P = Q k - 1$ ）には、上を繰り返せば、Q次の離散フーリエ変換と回転因子の掛け算を繰り返すことだけで次数を下げることができ、最終的に1次離散フーリエ変換（＝何もしないことと同じ）にまで下げると、F(t)を求めることができる。特に、Qが2又は4の場合は、Q次の離散フーリエ変換は非常に簡単な計算になる。

Q=2の場合は、 $exp( - 2π i r q / Q)$ は1か-1なので、Q次の離散フーリエ変換は、符号の逆転と足し算だけで計算できる。

Q=4の場合は、 $exp( - 2π i r q / Q)$ は1か-1かiか-iなので、Q次の離散フーリエ変換の計算は、符号の逆転、実部虚部の交換と足し算だけで計算できる。

このため、2の累乗あるいは4の累乗次の離散フーリエ変換は簡単に計算できる。実務的に用いられるのは、Q=2かQ=4の場合のみである。

[編集] ビットの反転

上記の説明で、 $N = Q k$ の場合（ $P = Q k - 1$ の場合）、 $N = Q k$ 個のデータ $f (q Q k - 1 + p)$ から、 $N = Q k$ 個の計算結果 $f_1(p,r)= \exp(-2\pi ipr/Q^k)\sum_{q=0}^{Q-1}f(qQ^{k-1}+p)\exp(-2\pi irq/Q)$ を計算する場合に、メモリーの節約のため、 $0 = < q < = Q - 1$ 、 $0 = < r < = Q - 1$ を利用し、

計算結果 $f 1 (p, r)$ を、元データ $f (r Q k - 1 + p)$ のあった場所に格納することが多い。

これが、次の次数、 $Q k - 1$ でも繰り返されるため、 $p = q 2 Q k - 2 + p 2$ とすると、次の次数の計算結果 $f 2 (p 2, q 2, q)$ は、 $f (q Q k - 1 + q 2 Q k - 2 + p 2)$ のあった場所に格納される。繰り返せば、 $t = q 1 Q k - 1 + q 2 Q k - 2,,, + q k$ とすると、計算結果 $f k (p k, q k, q k - 1,,, q 2, q 1)$ は、 $f (q 1 Q k - 1 + q 2 Q k - 2,,, + q k - 1 Q + p k)$ のあった場所に格納される。

一方、 $F(sQ+r)=\sum_{p=0}^{Q^{k-1}-1}\exp(-2\pi isp/Q^{k-1})f_1(p,r)$ を、rを固定しsを変数とした $Q k - 1$ 次離散フーリエ変換と見なして、 $s = s 2 Q + r 2$ とすると、 $F(s_2Q^2+r_2Q+r)=\sum_{p_2=0}^{Q^{k-2}-1}\exp(-2\pi is_2p_2/Q^{k-2})f_2(p_2,r_2,r)$ となる。繰り替えせば、 $F(s_kQ^k+r_kQ^{k-1}+,,,+r_2Q+r_1)=\sum_{p_k=0}^{Q^{k-k}-1}\exp(-2\pi is_kp_k/Q^{k-k})f_k(p_k,r_k,r_{k-1},,,r_2,r_1)$ となるが、 $s k Q k + r k Q k - 1 + ,,, + r 2 Q + r 1 < Q k$ より $s k = 0$ 、 $\sum_{p_k=0}^{Q^{k-k}-1}$ より $p k = 0$ このため、 $F (r k Q k - 1 + ,,, + r 2 Q + r 1) = f k (0, r k, r k - 1,,, r 2, r 1)$ これは $f (r 1 Q k - 1 + r 2 Q k - 2,,, + r k - 1 Q + r k)$ のあった場所に格納されている。

このように、求める解 $F (r k Q k - 1 + ,,, + r 2 Q + r 1)$ が $f (r 1 Q k - 1 + r 2 Q k - 2,,, + r k - 1 Q + r k)$ のあった場所に格納されているこを、ビット反転と言う。これは、Q進数で表示した場合、 $r k Q k - 1 + ,,, + r 2 Q + r 1$ は $(r k r k - 1,,, r 2 r 1) Q$ となるのに対し、 $r 1 Q k - 1 + r 2 Q k - 2,,, + r k - 1 + r k$ は逆から読んだ $(r 1 r 2,,, r k - 1 r k) Q$ となるためである。

[編集] その他のアルゴリズム

Prime Factor Algorithm (PFA)
Bruun's FFT algorithm
Rader's FFT algorithm
Bluestein's FFT algorithm

この節は執筆中です。加筆、訂正して下さる協力者を求めています。

[編集] 実数および対称的な入力への最適化

多くの実用面において、FFTへの入力は実数のみの列(実入力)であり、このとき出力の列は次の対称性を満たす (*は複素共役):

$f_{ - j} = f_j^* \quad j = 0, \dots ,N - 1$ .

そして、多くの効率的なFFTアルゴリズム（例えば^[3]）はこの条件を前堤に設計されている。

実入力への効率化には、次のような手段がある。

Cooley-Tukey型アルゴリズムなど典型的なアルゴリズムを利用して、時間とメモリーの両方の浪費を低減する。
偶数区分数のフーリエ係数はその半分の長さの複素フーリエ係数として表現することができる（出力の実数/虚数成分は、それぞれ入力の偶関数/奇関数成分に対応する）ことを利用する。

かつては実入力に対するフーリエ係数は離散ハートリー変換 (DHT : Discrete Hartley transform) でより効率的に計算できると思われていた。しかし、その後最適化された離散フーリエ変換(DFT：Discrete Fourier Transform)アルゴリズムが、離散ハートリー変換アルゴリズムよりも必要とする演算回数を下まわることがわかった。また、Bruun's FFTアルゴリズムは、はじめこそ実入力に対して有利といわれていたが、今ではそういわれてはいない。

さらに偶奇の対称性を持つ実入力の場合にはさらに最適化ができ、DFTはDCTやDSTになり、時間とメモリーの（ほとんど）2つに関して最適化が得られる。この場合には直接DFTのアルゴリズムを応用しないでも、DCTやDSTの計算によってフーリエ係数を求めることができる。

[編集] 応用

スペクトラムアナライザ
OFDM変復調器

日本および欧州で地上デジタルテレビジョン放送に用いられる変調方式であるOFDMの実装は、LSI化されたFFTおよびIFFT（逆変換）をそれぞれ復調器および変調器を用いて行われている。

CT、MRI等

受像素子を360度回転させながら連続撮影した映像をフーリエ変換する事により、回転面の透過画像を合成する。

多倍精度の乗除算
自動列車停止装置 (例: JR西日本の最新型車両。地上子が発振する周波数に高速フーリエ変換を用いている)
FFTアナライザ

周波数の分布を調べるために使用される。近年はデジタルオシロスコープにFFTの機能を内蔵している物もある。また、以前はハードウェアで信号を処理していたが近年はCPUの性能が向上した為ソフトウェアで処理している。また、ノートPCとUSBで接続して使用するものもある。

[編集] 参考文献

[ヘルプ]

^ ^a ^b J. W. Cooley and J. W. Tukey: Math. of Comput. 19 (1965) 297.
^ Carl Friedrich Gauss, "Nachlass: Theoria interpolationis methodo nova tractata," Werke band 3, 265?327 (Konigliche Gesellschaft der Wissenschaften, Gottingen, 1866). See also M. T. Heideman, D. H. Johnson, and C. S. Burrus, "Gauss and the history of the fast Fourier transform," IEEE ASSP Magazine 1 (4), 14?21 (1984).
^ H. V. Sorensen, D. L. Jones, M. T. Heideman, and C. S. Burrus, "Real-valued fast Fourier transform algorithms," IEEE Trans. Acoust. Speech Sig. Processing ASSP-35, 849–863 (1987).

[編集] 関連記事

フーリエ変換
離散フーリエ変換 (DFT)

[Math-0] J. W. Cooley and J. W. Tukey: Math. of Comput. 19 (1965) 297.

[1] Carl Friedrich Gauss, "Nachlass: Theoria interpolationis methodo nova tractata," Werke band 3, 265?327 (Konigliche Gesellschaft der Wissenschaften, Gottingen, 1866). See also M. T. Heideman, D. H. Johnson, and C. S. Burrus, "Gauss and the history of the fast Fourier transform," IEEE ASSP Magazine 1 (4), 14?21 (1984).

[2] H. V. Sorensen, D. L. Jones, M. T. Heideman, and C. S. Burrus, "Real-valued fast Fourier transform algorithms," IEEE Trans. Acoust. Speech Sig. Processing ASSP-35, 849–863 (1987).

[1]

[2]

[3]