高速フーリエ変換
高速フーリエ変換(こうそくフーリエへんかん、fast Fourier transform, FFT)は、離散フーリエ変換(discrete Fourier transform, DFT)を計算機上で高速に計算するアルゴリズムである。高速フーリエ変換の逆変換を逆高速フーリエ変換(inverse fast Fourier transform, IFFT)と呼ぶ。
複素関数 f(x) の離散フーリエ変換である複素関数 F(t) は以下で定義される。
F(t)= \sum_{x=0}^{N-1} f(x) \exp\left(-i\frac{2 \pi t x}{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) = \sum_{x=0}^{N-1} f(x) \exp\left(-i\frac{2 \pi t x}{N}\right)
で定義したとき、逆変換は
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 t x}{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) = \sum_{x=0}^{N-1} f(x) W_N^{tx} .
例えば、N = 4 のとき、F(t)=X_t、f(k
) = x_kとすれば、離散フーリエ係数は行列を用いて表現すると(W = W4 と略記)
\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 を添字の偶奇で分けて、以下のように変形する。
\begin{align}
\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^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{align}
(\because W^{k+N}=W^k)
すると、サイズ 2 のFFTの演算結果を用いて表現でき、サイズの分割ができる。
\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\\
\end{bmatrix}
また、この分割手順を図にすると蝶のような図になることから、バタフライ演算とも呼ばれる。
バタフライ演算は、計算機上ではビット反転で実現される。DSPの中には、このバタフライ演算のプログラムを容易にするため、ビット反転アドレッシングを備えているものがある。
N = PQ とする。N 次離散フーリエ変換を、以下のようにP 次離散フーリエ変換とQ 次離散フーリエ変換に分解する。
N 次離散フーリエ変換:
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 である。
\begin{align}
n &= sQ+r \qquad 0\le s\le P-1,\,0\le r\le Q-1 \\
k &= qP+p \qquad 0\le p\le P-1,\;0\le q\le Q-1
\end{align}
\begin{align}
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{align}
ここで、
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(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\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次の離散フーリエ変換
\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(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}) の場合、N = Qk 個のデータf(qQ^{k-1}+p)から、N = Qk 個の計算結果
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) を元データf(rQ^{k-1}+p) のあった場所に格納することが多い。これが次の次数  Qk − 1 でも繰り返されるため、p=q_2Q^{k-2}+p_2とすると、次の次数の計算結果f_2(p_2,q_2,q)はf(qQ^{k-1}+q_2Q^{k-2}+p_2)のあった場所に格納される。繰り返せば、t=q_1Q^{k-1}+q_2Q^{k-2}+\cdots+q_kとすると、計算結果f_k(p_k,q_k,q_{k-1},\dots,q_2,q_1)はf(q_1Q^{k-1}+q_2Q^{k-2}+\cdots+q_{k-1}Q+p_k)のあった場所に格納される。
一方、
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_2Q+r_2とすると、
F(s_2Q^2+r_2Q+r) = \sum_{p_2=0}^{Q^{k-2}-1}\exp\left(-2\pi i\frac{s_2p_2}{Q^{k-2}\right)f_2(p_2,r_2,r)
となる。繰り替えせば、
F(s_kQ^k+r_kQ^{k-1}+\cdots+r_2Q+r_1)
= \sum_{p_k=0}^{Q^{k-k}-1}\exp\left(-2\pi i\frac{s_kp_k}{Q^{k-k}\right)f_k(p_k,r_k,r_{k-1},\dots,r_2,r_1)
となるが、左辺について
s_kQ^k+r_kQ^{k-1}+\cdots+r_2Q+r_1<Q^k
より sk = 0, また右辺について
Q^{k-k}-1=0
より pk = 0。このため、
F(r_kQ^{k-1}+\cdots+r_2Q+r_1)=f_k(0,r_k,r_{k-1},\dots,r_2,r_1).
これは f(r_1Q^{k-1}+r_2Q^{k-2}+\cdots+r_{k-1}Q+r_k) のあった場所に格納されている。
このように、求める解 F(r_kQ^{k-1}+\cdots+r_2Q+r_1) が f(r_1Q^{k-1}+r_2Q^{k-2}+\cdots+r_{k-1}Q+r_k) のあった場所に格納されていることを、ビット反転と言う。これは、Q 進法で表示した場合、r_kQ^{k-1}+\cdots+r_2Q+r_1 は (r_kr_{k-1}\dots r_2r_1)_Qとなるのに対し、r_1Q^{k-1}+r_2Q^{k-2}+\cdots+r_{k-1}+r_kは逆から読んだ(r_1r_2\dots r_{k-1}r_k)_Qとなるためである。
以下は、高速フーリエ変換のプログラムを Q = 4 の場合にMicrosoft Visual Basicの文法を用いて書いた例である。
この例では、最深部 (For k、Next k の間の部分)の繰り返し回数が Ndeg log4 Ndeg となっている。
Prime Factor Algorithm (PFA)
Bluestein's FFT algorithm (see "Chirp Z-transform") 任意長のデータ列に対する変換が高速に可能である。
オドリツコ・ショーンハーゲ法 - アンドリュー・オドリツコ、アーノルド・ショーンハーゲ。
多くの応用において、FFTに対する入力データは実数の列(実入力)であり、このとき変換された出力の列は次の対称性を満たす(    は複素共役):
F(-t) = \overline{F(t)} .
そこで、多くの効率的なFFTアルゴリズム は入力データが実数であることを前提に設計されている。
入力データが実数の場合の効率化の手段としては、次のようなものがある。
クーリー-テューキー型アルゴリズムなど典型的なアルゴリズムを利用して、時間とメモリーの両方のコストを低減する。
入力データが偶数の長さのフーリエ係数はその半分の長さの複素フーリエ係数として表現できる(出力の実数/虚数成分は、それぞれ入力の偶関数/奇関数成分に対応する)ことを利用する。
かつては実数の入力データに対するフーリエ係数を求めるのには、実数計算だけで行える離散ハートリー変換 (discrete Hartley transform, DHT)を用いると効率的であろうと思われていた。しかしその後に、最適化された離散フーリエ変換 (discrete Fourier transform, DFT) アルゴリズムの方が、離散ハートリー変換アルゴリズムに比べて必要な演算回数が少ないということが判明した。また当初は、実数入力に対してブルーン (Bruun) FFT アルゴリズムは有利であると云われていたが、その後そうではないことが判った。
また、偶奇の対称性を持つ実入力の場合には、DFTはDCTや離散サイン変換となるので、演算と記憶に関してほぼ2倍の効率化が得られる。よって、そのような場合にはDFTのアルゴリズムをそのまま適用するよりも、DCTやDSTを適用してフーリエ係数を求める方が効率的である。
OFDM(日本および欧州で地上デジタルテレビジョン放送やADSL等に用いられる変調方式)の実装は、LSI化されたFFTおよびIFFT(逆変換)をそれぞれ復調器および変調器を用いて行われている。
核磁気共鳴 (NMR) スペクトルを得るために使用される。
コンピュータ断層撮影 (CT)、核磁気共鳴画像法 (MRI) 等
受像素子を360度回転させながら連続撮影した映像をフーリエ変換する事により、回転面の透過画像を合成する。
自動列車停止装置 (例: JR西日本の最新型車両。地上子が発振する周波数の検出に、高速フーリエ変換が用いられている)
周波数の分布を調べるために使用される。以前はハードウェアで信号を処理していたが、近年はCPUの性能が向上した為ソフトウェアで処理される。ノートパソコンとUSBで接続して使用するもの や、近年はデジタルオシロスコープにFFTの機能を内蔵している物もある。
FX型デジタル分光相関器等を使用して星間分子のスペクトルを解析する。
高速フーリエ変換といえば一般的には1965年、ジェイムズ・クーリー (J. W. Cooley) とジョン・テューキー (J. W. Tukey) が発見した とされているクーリー–テューキー型FFTアルゴリズムを呼ぶ。同時期に高橋秀俊がクーリーとテューキーとは全く独立にフーリエ変換を高速で行うためのアルゴリズムを考案していた。しかし、1805年頃に既にガウスが同様のアルゴリズムを独自に発見していた(本ページの外部リンク先に同じ文章PDFへのリンクがある)。ガウスの論文以降、地球物理学や気候や潮位解析などの分野などで測定値に対する調和解析は行われていたので、計算上の工夫を必要とする応用分野で受け継がれていたようである(たとえば、Robart L. Nowack: "Development of the FFT and Applications in Geophysics", in Proceedings of the Cornelius Lanczos International Centenary Conference,SIAM, ISBN 978-0898713398 (1994), pp.395--397、の中では Danielson and Lanczos(1942年)などの先行例をあげている。和書でも沼倉三郎:「測定値計算法」、森北出版、(1956年)には,一般の合成数Nに対してではないが人が計算を行う場合にある程度の大きさの合成数Nに対してどのように計算すればよいかについての説明をみることができる)。
以下の書籍にも、天体観測の軌道の補間のためにガウスが高速フーリエ変換を利用したことが書かれている。
Elena Prestini:"The Evolution of Applied Harmonic Analysis", Springer, ISBN 978-0-8176-4125-2 (2004)のSec.3.10 'Gauss and the asteroids: history of the FFT'.
離散フーリエ変換 (DFT)
ショーンハーゲ・ストラッセン法(乗算アルゴリズム)
Sparse Fourier transform ※ (高速)スパース・フーリエ変換(SFT)
Henri J. Nussbaumer: "Fast Fourier Transform and Convolution Algorithms",2nd Ed.,Springer-Verlag, ISBN 978-3-540-11825-1 (1982年).
E.Oran Brigham:「高速フーリエ変換」、科学技術出版社 (1985年).
Henri J. Nussbaumer:「高速フーリエ変換のアルゴリズム」、科学技術出版社、ISBN 978-4876530069 (1989年).
William L. Briggs and Van Emden Henson: "The DFT: An Owners' Manual for the Discrete Fourier Transform", SIAM, ISBN 978-0-898713-42-8 (1995年).
Gerlind Plonka, Daniel Potts, Gabriele Steidl and Manfred Tasche: "Numerical Fourier Analysis", Birkhaeuser, ISBN 978-3030043056 (2019年2月).
谷萩隆嗣:「高速アルゴリズムと並列信号処理」、コロナ社、ISBN 4-339-01124-X(2000年7月26日)。