高速フーリエ変換

高速フーリエ変換（こうそくフーリエへんかん、英: fast Fourier transform, FFT）は、離散フーリエ変換（英: discrete Fourier transform, DFT）を計算機上で高速に計算するアルゴリズムである。高速フーリエ変換の逆変換を逆高速フーリエ変換（英: inverse fast Fourier transform, IFFT）と呼ぶ。

概要

複素関数 f(x) の離散フーリエ変換である複素関数 F(t) は以下で定義される。

${\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(N²) である。

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

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

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

逆変換

逆変換は正変換と同じと考えて良いが、指数の符号が逆であり、係数 1/N が掛かる。

高速フーリエ変換のプログラム中、どの符号が逆転するかを一々分岐させると、分岐の判定に時間がかかり、パフォーマンスが落ちる。一方、正変換のプログラムと、逆変換のプログラムを両方用意しておくことも考えられるが、共通部分が多いため、無駄が多くなる。このため、複素共役を使った次のような方法が考えられる。

離散フーリエ変換を

${\displaystyle F(t)=\sum _{x=0}^{N-1}f(x)\exp \left(-i{\frac {2\pi tx}{N}}\right)}$

で定義したとき、逆変換は

${\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アルゴリズム

→詳細は「クーリー–テューキー型FFTアルゴリズム（英語版）」を参照

クーリー–テューキー型アルゴリズムは、代表的な高速フーリエ変換 (FFT) アルゴリズムである。

分割統治法を使ったアルゴリズムで、N = N₁ N₂ のサイズの変換を、より小さいサイズである N₁, N₂ のサイズの変換に分割していくことで高速化を図っている。

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

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

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

原理の簡単な説明

データ数1２の離散フーリエ変換の模式図。時計を模した図形は1の12乗根の一つを表している。時計の針の向きと色は1の12乗根の偏角を表す。この図で表される行列をデータ列にかけることで離散フーリエ変換が得られる。上図で表されるような列の並べ替えを行うことで、元の行列のパターンはデータ数6の離散フーリエ変換のパターンに分解できる。この繰り返しにより最終的にはデータ数３のフーリエ変換に帰着される。

データ数100の離散フーリエ変換の模式図。色は1の100乗根の偏角を表す。バタフライ演算により元の行列のパターンは最終的にデータ数5の離散フーリエ変換のパターンに分解される。

FFTのバタフライ演算

離散フーリエ係数は、1の原始 N 乗根の1つ W_N = e^−2πi/N を使うと、次のように表せる。

${\displaystyle F(t)=\sum _{x=0}^{N-1}f(x)W_{N}^{tx}.}$

例えば、N = 4 のとき、${\displaystyle F(t)=X_{t}}$、${\displaystyle f(k)=x_{k}}$とすれば、離散フーリエ係数は行列を用いて表現すると（W = W₄ と略記）

${\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}}}$

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

${\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}}}$

(${\displaystyle \because W^{k+N}=W^{k}}$)

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

${\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の中には、このバタフライ演算のプログラムを容易にするため、ビット反転アドレッシングを備えているものがある。

原理の説明

FFTの原理は、N = PQ として、N 次離散フーリエ変換を、P 次離散フーリエ変換とQ 次離散フーリエ変換に分解することである^[2]。

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 である。

${\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}}}$

すると

${\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}}}$

ここで、

${\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)}$

と置くと、

${\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 について

${\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次の離散フーリエ変換

${\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 について

${\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 次離散フーリエ変換の問題　⇒　Q 次離散フーリエ変換の実施　⇒　N/Q 次離散フーリエ変換の問題に帰着

特に、Q が2または4の場合は、Q次離散フーリエ変換は非常に簡単な計算になる。

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

Q = 2 かQ = 4 の場合のこの部分のQ次離散フーリエ変換のことを、バタフライ演算と言う。

また､N = Q^k (P = Q^{k − 1}) の場合には、上を繰り返すことができ、Q 次の離散フーリエ変換と回転因子の掛け算を繰り返すことだけで次数を下げられ、最終的に1次離散フーリエ変換（何もしないことと同じ）にまで下げると、F(t)を求めることができる。
このため、2の累乗あるいは4の累乗次の離散フーリエ変換は、両方の性質を利用でき、非常に簡単に計算できる。

また、Q = 2かQ = 4の場合において、計算を終了するまでに何回の「掛け算」が必要かを考える。符号の逆転、実部虚部の交換は「掛け算」として数えなければ、回転因子の掛け算のみが「掛け算」である。N = Q^kの次数を1落とすためにN回の「掛け算」が必要であり、次数をkから0に落とすにはそれをk回繰り返す必要があるため、「掛け算」の数は Nk = N log_Q N となる。高速フーリエ変換の計算において時間がかかるのは「掛け算」の部分であるため、これが「高速フーリエ変換では計算速度は O(N log N) になる」ことの根拠になっている。

ビットの反転

上記の説明で、${\displaystyle N=Q^{k}(P=Q^{k-1})}$ の場合、N = Q^k 個のデータ${\displaystyle f(qQ^{k-1}+p)}$から、N = Q^k 個の計算結果

${\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 を利用し、計算結果 ${\displaystyle f_{1}(p,r)}$ を元データ${\displaystyle f(rQ^{k-1}+p)}$ のあった場所に格納することが多い。これが次の次数 Q^{k − 1} でも繰り返されるため、${\displaystyle p=q_{2}Q^{k-2}+p_{2}}$とすると、次の次数の計算結果${\displaystyle f_{2}(p_{2},q_{2},q)}$は${\displaystyle f(qQ^{k-1}+q_{2}Q^{k-2}+p_{2})}$のあった場所に格納される。繰り返せば、${\displaystyle t=q_{1}Q^{k-1}+q_{2}Q^{k-2}+\cdots +q_{k}}$とすると、計算結果${\displaystyle f_{k}(p_{k},q_{k},q_{k-1},\dots ,q_{2},q_{1})}$は${\displaystyle f(q_{1}Q^{k-1}+q_{2}Q^{k-2}+\cdots +q_{k-1}Q+p_{k})}$のあった場所に格納される。

一方、

${\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 を変数とした Q^{k − 1} 次離散フーリエ変換と見なして、${\displaystyle s=s_{2}Q+r_{2}}$とすると、

${\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)}$

となる。繰り替えせば、

${\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})}$

となるが、左辺について

${\displaystyle s_{k}Q^{k}+r_{k}Q^{k-1}+\cdots +r_{2}Q+r_{1}<Q^{k}}$

より s_k = 0, また右辺について

${\displaystyle Q^{k-k}-1=0}$

より p_k = 0。このため、

${\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}).}$

これは ${\displaystyle f(r_{1}Q^{k-1}+r_{2}Q^{k-2}+\cdots +r_{k-1}Q+r_{k})}$ のあった場所に格納されている。

このように、求める解 ${\displaystyle F(r_{k}Q^{k-1}+\cdots +r_{2}Q+r_{1})}$ が ${\displaystyle f(r_{1}Q^{k-1}+r_{2}Q^{k-2}+\cdots +r_{k-1}Q+r_{k})}$ のあった場所に格納されていることを、ビット反転と言う。これは、Q 進法で表示した場合、${\displaystyle r_{k}Q^{k-1}+\cdots +r_{2}Q+r_{1}}$ は ${\displaystyle (r_{k}r_{k-1}\dots r_{2}r_{1})_{Q}}$となるのに対し、${\displaystyle r_{1}Q^{k-1}+r_{2}Q^{k-2}+\cdots +r_{k-1}+r_{k}}$は逆から読んだ${\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 ^ (deg - i)
 For j0 = 0 To 4 ^ (i - 1) - 1
  For j1 = 0 To Pdeg - 1
   j = j1 + j0 * Pdeg * 4
   'バタフライ演算(Q=4)
   w1 = CR(j) + CR(j + Pdeg) + CR(j + 2 * Pdeg) + CR(j + 3 * Pdeg)
   w2 = CI(j) + CI(j + Pdeg) + CI(j + 2 * Pdeg) + CI(j + 3 * Pdeg)
   w3 = CR(j) + CI(j + Pdeg) - CR(j + 2 * Pdeg) - CI(j + 3 * Pdeg)
   w4 = CI(j) - CR(j + Pdeg) - CI(j + 2 * Pdeg) + CR(j + 3 * Pdeg)
   w5 = CR(j) - CR(j + Pdeg) + CR(j + 2 * Pdeg) - CR(j + 3 * Pdeg)
   w6 = CI(j) - CI(j + Pdeg) + CI(j + 2 * Pdeg) - CI(j + 3 * Pdeg)
   w7 = CR(j) - CI(j + Pdeg) - CR(j + 2 * Pdeg) + CI(j + 3 * Pdeg)
   w8 = CI(j) + CR(j + Pdeg) - CI(j + 2 * Pdeg) - CR(j + 3 * Pdeg)
   CR(j) = w1
   CI(j) = w2
   CR(j + Pdeg) = w3
   CI(j + Pdeg) = w4
   CR(j + 2 * Pdeg) = w5
   CI(j + 2 * Pdeg) = w6
   CR(j + 3 * Pdeg) = w7
   CI(j + 3 * Pdeg) = w8
   '回転因子
   For k = 0 To 3
    w1 = Cos(2 * pi * j * k / Pdeg / 4)
    w2 = -Sin(2 * pi * j * k / Pdeg / 4)
    w3 = CR(j + k * Pdeg) * w1 - CI(j + k * Pdeg) * w2
    w4 = CR(j + k * Pdeg) * w2 + CI(j + k * Pdeg) * w1
    CR(j + k * Pdeg) = w3
    CI(j + k * Pdeg) = w4
   Next k
  Next j1
 Next j0
Next i
'ビット反転
For i = 0 To Ndeg - 1
 k = i
 k1 = 0
 For j = 1 To deg
  k1 = k1 + (k - Int(k / 4) * 4) * 4 ^ (deg - j)
  k = Int(k / 4)
 Next j
 FR(i) = CR(k1)
 FI(i)=CI(k1)
Next i

この例では、最深部 (For k、Next k の間の部分)の繰り返し回数が Ndeg log₄ Ndeg となっている。

その他のアルゴリズム

• Prime Factor Algorithm（英語版） (PFA)
• Bruun's FFT algorithm（英語版）
• レーダーのFFTアルゴリズム
• Bluestein's FFT algorithm（英語版） (see "Chirp Z-transform")　任意長のデータ列に対する変換が高速に可能である。
• オドリツコ・ショーンハーゲ法（英語版） - アンドリュー・オドリツコ（英語版）、アーノルド・ショーンハーゲ（英語版）。
• Fast Walsh–Hadamard transform（英語版）

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

多くの応用において、FFTに対する入力データは実数の列(実入力)であり、このとき変換された出力の列は次の対称性を満たす（    は複素共役）：

${\displaystyle F(-t)={\overline {F(t)}}.}$

そこで、多くの効率的なFFTアルゴリズム^[3] は入力データが実数であることを前提に設計されている。

入力データが実数の場合の効率化の手段としては、次のようなものがある。

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

かつては実数の入力データに対するフーリエ係数を求めるのには、実数計算だけで行える離散ハートリー変換（英語版） (discrete Hartley transform, DHT)を用いると効率的であろうと思われていた。しかしその後に、最適化された離散フーリエ変換 (discrete Fourier transform, DFT) アルゴリズムの方が、離散ハートリー変換アルゴリズムに比べて必要な演算回数が少ないということが判明した。また当初は、実数入力に対してブルーン (Bruun) FFT アルゴリズムは有利であると云われていたが、その後そうではないことが判った。

また、偶奇の対称性を持つ実入力の場合には、DFTはDCTやDST（英語版）となるので、演算と記憶に関してほぼ2倍の効率化が得られる。よって、そのような場合にはDFTのアルゴリズムをそのまま適用するよりも、DCTやDSTを適用してフーリエ係数を求める方が効率的である。

応用

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

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

• フーリエ変換NMR

  核磁気共鳴 (NMR) スペクトルを得るために使用される。

• コンピュータ断層撮影 (CT)、核磁気共鳴画像法 (MRI) 等

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

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

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

• 電波天文学

  FX型デジタル分光相関器等を使用して星間分子のスペクトルを解析する^[5][6]。

歴史

高速フーリエ変換といえば一般的には1965年、ジェイムズ・クーリー（英語版） (J. W. Cooley) とジョン・テューキー (J. W. Tukey) が発見した^[1] とされているクーリー–テューキー型FFTアルゴリズム（英語版）のことをさす^[7]。同時期に高橋秀俊がクーリーとテューキーとは全く独立にフーリエ変換を高速で行うためのアルゴリズムを考案していた^[8]。しかし、1805年頃に既にガウスが同様のアルゴリズムを独自に発見していた。それは彼の没後に刊行された全集に収録されている^[9]（本ページの外部リンク先に同じ文章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'.

ライブラリ

特定のデバイスに限定していない汎用の実装

• FFTPACK
• FFTW

ハードウェアベンダーによる、特定のデバイス向けの実装

• Accelerate vDSP - Apple のデバイス用^[10]
• AOCL-FFTW - AMD CPU 用^[11]
• Arm Performance Libraries - Arm64 用^[12]
• cuFFT - NVIDIA GPU 用^[13]
• Intel oneAPI Math Kernel Library - Intel CPU, GPU 用
• MathKeisan - NEC SX 用^[14]
• rocFFT - AMD GPU 用^[15]