カバー画像：合成信号とその周波数成分を表す、AI生成の概念的なイラスト。当時の計算機や実験データを描いたものではない。
信号の中にある成分を聴き分ける
音楽の和音を録音すると想像してみよう。マイクが返すのは振幅の列だが、知りたいのはその中にある音である。フーリエ解析は、この二つの表し方を数学的につなぐ。振動や画像の模様を調べるときにも、似た問いが生じる。
離散フーリエ変換（DFT）は、有限の列を離散的な周波数成分で表す。高速フーリエ変換（FFT）は、その同じ変換を効率よく計算する方法である。新しい種類のスペクトルを作るわけではない。NumPyの技術文書は、入力を時間領域の信号、出力をその周波数領域での表現として説明している。DFTの文書
実用上の難しさは、計算の繰り返しにある。直接計算では、入力の各点を出力の各周波数と組み合わせる。点数を2倍にすると、主要な演算量は4倍になる。科学上の問いを扱うために数千点が必要なら、数学的には簡単な式でも、計算上の障害になりうる。
クーリーとテューキーは、一般的な因数分解と、長さが2のべき乗の場合に便利な実装を示した。入力配列の記憶領域を再利用する方法も含まれていた。論文は抽象的な提案だけでなく、IBM 7094向けのプログラムと計算時間を報告している。編集部が原稿を受け取ったのは1964年8月17日で、Mathematics of Computationの1965年4月号に掲載された。正確な発行日は、ここでは確認できていない。原論文、大学所蔵の写し
長い前史を持つ転機
この考えは、1965年に突然生まれたわけではない。歴史研究者のマイケル・ハイデマン、ドン・ジョンソン、シドニー・バラスは過去の文献を調べ、ガウスの研究に同等の分解を見いだした。彼らは執筆時期を1805年ごろと推定したが、原稿自体に明示された日付はなく、死後の1866年に公刊された。この研究では、ダニエルソンとランチョスによる1942年の研究や、グッドの素因数を用いた方法など、その後の手法も論じている。グッドの方法の制約は、一般的なクーリー–テューキーの因数分解とは異なる。歴史的な調査と閲覧可能な著者版
1993年の回想で、クーリーとテューキーは、研究を結びつけ、開発を後押ししたリチャード・ガーウィンの役割を述べている。以前の発見も認めたうえで、大規模な電子計算がなぜこの方法の価値を高めたのかを説明している。これは当事者による回顧的な記述であり、記憶に基づく説明の限界がある。クーリーとテューキーの回想
長く影響を残した転機は、再利用できる計算法が実務に広まったことだった。IEEEのマイルストーンを紹介するプリンストン大学の記事は、それを科学計算、信号処理、医用画像、データ伝送と結びつけている。これらの応用には、計測機器、物理モデル、ほかのアルゴリズムも必要であり、FFTは強力な計算上の構成要素として貢献した。後年の顕彰は歴史的重要性を裏づけるものであって、最近の記念行事の日付が発見の日になるわけではない。プリンストン大学とIEEEのマイルストーン
4点の例：一般式の前に再利用の仕組みを見る
無次元の列[1, 2, 3, 4]を使う。これは意図的に小さくした教育用の計算であり、1965年の論文のデータではない。iを虚数単位とし、i² = −1とする。尺度係数を掛けず、指数の符号を負とするDFTは、次の4個の出力を与える。
| 周波数の添字 | 出力 |
|---|---|
| 0 | 10 |
| 1 | −2 + 2i |
| 2 | −2 |
| 3 | −2 − 2i |
最初の出力は、四つの入力すべての和である。残りの出力は、正、負、虚数の重みを組み合わせる。複素数の出力は二つの成分を記録し、その大きさと位相が一つの周波数成分をともに表す。
ここで、添字が偶数の入力[1, 3]と、奇数の入力[2, 4]に分ける。それぞれの組に必要なのは、和と差だけである。
- 偶数側の組：E = [4, −2]。
- 奇数側の組：O = [6, −2]。
これらの短い結果を組み合わせる。最初のバタフライ演算は4 + 6 = 10と4 − 6 = −2を作る。二つ目では、まず奇数側の結果に−iを掛けて回転させ、(−i)(−2) = 2iを得る。その二つの出力は−2 + 2iと−2 − 2iとなる。
バタフライとは、このように和と差を対にして求める演算の名称である。図の交差する線は計算の整理を表し、波の間の物理的な結合ではない。二つの出力は同じ中間結果を再利用する。規模が大きくなれば、短い変換をさらに分割できる。再利用は一度だけの近道ではなく、階層になる。
具体例から計算量の増え方へ
入力値がN個あるとき、本記事の規約を次のように定める。
Xk=n=0∑N−1xne−2πink/N.
nは入力、kは出力の添字であり、どちらも0からN − 1まで動く。各因子は、大きさが1の複素数による回転である。本記事の順変換には正規化係数を付けず、逆変換には正の指数と1/Nの係数を使う。ほかの規約もある。原論文は正の指数を持つフーリエ和から出発しているため、符号と正規化を明示的に比較する必要がある。現代の規約
Nが偶数なら、和を偶数添字と奇数添字の入力に分ける。それぞれの長さN/2の変換をEₖとOₖと書き、W_N = exp(−2πi/N)と置く。元の変換は、0 ≤ k < N/2に対して次のように再構成される。
Xk=Ek+WNkOk,Xk+N/2=Ek−WNkOk.
二つ目の等式は、円を半周すると回転因子の符号が反転することを使っている。したがって、二つの短い変換を計算し、Nに比例する回数の組み合わせを行えばよい。N = 2ᵐでこの分割を繰り返すと、m = log₂N段になる。
各段の作業量は列の長さに比例する。そのため、全体はN log₂Nに従って増える。O(N log N)という記法はこの増え方を表し、プロセッサの命令数や秒数を正確に数えたものではない。ここで説明した基数2の構成は2のべき乗に限られるが、一般的なクーリー–テューキーの因数分解では、ほかの合成数の長さも扱える。
実機での速度は、データの配置、メモリを通る転送、実装にも依存する。演算回数が少ないというだけでは、手元の機器で実際の経過時間がどれだけ短くなるかは確定しない。
高速化では解決できないこと
等間隔の測定の時間間隔がΔt秒なら、標本化周波数は1/Δtヘルツ、周波数ビンの間隔は1/(NΔt)ヘルツである。大きい側の出力添字を解釈するには、正と負の周波数にどの並び順を採用したかを考慮する必要がある。周波数の規約
有限の標本化による推論の限界は残る。ここでtは秒単位の時間、fはヘルツ単位の信号の周波数、fₛ = 1/Δtはヘルツ単位の標本化周波数である。たとえば、t = n/fₛで取ったcos(2πft)の標本は、fをf + fₛに置き換えても変わらない。増えた位相が2πnだからである。この標本だけからは、どれほど高速なアルゴリズムでも二つの信号を区別できない。雑音や不十分な測定モデルも、依然として測定の問題である。
数値演算には別の限界もある。複素数の回転因子には一般に数値的な近似が必要であり、演算の順序を変えると丸めも変わりうる。FFTWの精度に関する解説は、回転因子の精度の重要性を強調し、実装や環境間の違いを記録している。数学的に等価でも、浮動小数点数のビット列まで同一になる保証はない。精度に関する解説
教育用の計算を確かめる
本記事の検査にはCPython 3.12.10と標準ライブラリだけを使った。4点の根[1, −i, −1, i]は厳密に表現できるため、この小さな整数入力には完全一致という判定が適切である。この基準を、長い浮動小数点変換へ機械的に適用してはならない。
記録した検査では、次の中心的な計算を実行した。
ROOTS = (1, -1j, -1, 1j)
def direct(x, roots=ROOTS):
return [
sum(x[n] * roots[(n*k) % 4] for n in range(4))
for k in range(4)
]
def butterfly(x):
even = (x[0] + x[2], x[0] - x[2])
odd = (x[1] + x[3], x[1] - x[3])
return [
even[0] + odd[0], even[1] - 1j*odd[1],
even[0] - odd[0], even[1] + 1j*odd[1],
]両方の関数を、[1, 2, 3, 4]、インパルス**[1, 0, 0, 0]、一定の列[1, 1, 1, 1]**、ゼロの列、[1+i, −2, 3−i, 2i]で比較する。五つすべてが完全に一致した。最初のケースは上の表とも一致した。根の虚数部分の符号を反転すると、非対称なケースのスペクトルが変わり、検査はそれを異なる規約として不合格にした。
これで検証されたのは、固定された教育用のケースである。あらゆる実装の正しさを証明したわけでも、歴史的な性能を独立に再現したわけでもない。原論文には、そのような再現に必要なプログラムの全体、ベンチマーク入力、完全な時間測定条件がそろっていない。
この記事は、引用した資料に基づいてAIシステムが執筆・翻訳し、Newsの編集手法に従って自動検査を行った。教育用の計算は記載のとおり実行されたが、これは査読、人間による科学的な検証、歴史的なプログラムの再現を意味しない。