高速フーリエ変換

N個のサンプルの離散フーリエ変換は、N個の項からなる和をN回計算するものです。この和をより小さな和に分割することで、はるかに少ない計算コストで同じ答えが得られます。バタフライダイアグラムで分割の過程を観察するか、小惑星パラスに関する12の数値を用いてガウスが行った方法を見てみましょう。

計算例 — 高速フーリエ変換

何分の1か	5.3

ページを開いたときに計算機が表示する値です。上の入力を変えると再計算されます。

インタラクティブシミュレーションを読み込んでいます...

レッスン

理論 — 高速フーリエ変換

高速フーリエ変換（FFT）とは、丸め誤差を除いて離散フーリエ変換 X[k] = Σ x[n] W^kn （W = e^−2πi/N）を正確に計算するアルゴリズムであり、N² 回の代わりに約 N log N 回の演算で済みます。クーリー・テューキー型は N = N₁N₂ と分割します。まず長さ N₁ の変換を N₂ 個インターリーブし、変換 n₂ の各出力 k₁ に固有の回転因子 W^n₂k₁ を掛け、次にそれらをまたぐように長さ N₂ の変換を N₁ 個実行します。これを N = 2^m に対して繰り返し適用したものが基数2のFFTであり、N/2 個のバタフライからなる m 個のステージになります。

各記号の意味

N: サンプル数、および計算される周波数の数
W: 1の冪根 e^−2πi/N。この累乗が回転因子になる
m: log₂N、基数2のFFTのステージ数

ガウスはクーリーとテューキーの約160年前にFFTの分割を書き留めていた 🖖

ガウス、長さ4の周期×3 を押してください。これらの12の数値はガウスの手稿に基づくもので、小惑星パラスが移動しうる天球上の帯の北縁の赤緯を、赤経30°ごとに分角単位で記したフォン・ツァハの雑誌からの引用です。彼はこれらの値を通る三角級数を求めていました。彼は1つの12項の和を計算する代わりに、3つおきの値（0°、90°、180°、270°、次に 30°、120°、…、次に 60°、…）を取り出し、4つの値からなる各セットを独立して当てはめ、3つの短い結果を結合しました。破線の曲線はそれぞれの当てはめを示していますが、これらは一致しません。4点だけではより速い波打ちを捉えきれないため、セットごとに異なる結果が得られます。この不一致こそが、彼を結合方法の探究へと向かわせた理由でした。結果の下にある表は、彼の計算を再実行したものです。最初の周期からは 779.5、−415.0、−726.0、43.5 という手稿内の数値が得られ、結合された級数は直接求めた12項の級数になります。彼は第2の方法として長さ3の周期を4つ用いて計算し、同じ級数にたどり着きました。1984年にハイデマン、ジョンソン、バラスが示したように、ここでのインデックス演算はクーリーとテューキーが1965年に発表したものと同一です。同時に、この手稿に書かれていないことも重要です。省かれた計算量については全く言及されていません。ガウスは分割した方が速いと述べている箇所で、その証明を実践に委ね、「試す者が実践の中で知るだろう（praxis tentantem docebit）」としています。計算例には 12 = 3 × 4 と 36 = 6 × 6 の2つの因数分解があり、第27項では分割がさらに拡張できるとだけ述べられています。この手稿は彼の生前には発表されず、1866年の著作集でラテン語で掲載されましたが、日付はありませんでした。著作集の編集者はこれを1805年10月に着手されたと思われる草稿の改訂版と呼んでおり、ハイデマン、ジョンソン、バラスは一連の手がかりからこれを1805年の10月か11月と特定しています。その手がかりの1つが、手稿に使われているジュノーの離心率 0.254236 であり、これは1805年5月に出版されたものでした。

削減量は固定ではなく、Nとともに増加する 🖖

1,024点 を押してください。直接和はすべてのサンプルにすべての回転因子を掛けるため、N² = 1,048,576 回の乗算を行います。FFTは、それぞれ1回の乗算を含む N/2 = 512 個のバタフライからなる log₂N = 10 のステージを実行するため、5,120 回になります。これは 204.8 分の1の回数であり、比率は 2N/log₂N となるため、この差は広がり続けます。100万点 を押すと、直接和の 1,099,511,627,776 回の乗算に対して 10,485,760 回となり、比率は 104,857.6 になります。Nを2倍にすると、直接和の計算量は4倍になります。一方FFTの計算量は1,000点の時点で 2.2 倍の増加にとどまり、Nが大きくなるほど2倍に近づいていきます。8点ではわずか 5.3 倍しか削減されないため、ダイアグラム上では複雑な配線の割に利点が少なく見えます。どちらの回数にも W⁰ = 1 による乗算が含まれており、クーリーとテューキー自身の N² も同様です。彼らは乗算と加算のペアを数えて 2N log₂N 未満という結果を得ましたが、いずれにせよ比率は N/log N に比例して増加します。どのようなコンピュータを用いても、回数自体やその比率が変わることはありません。高速なコンピュータは実行時間を変えますが、必ずしも両方が同じ係数で速くなるわけではありません。表の下にあるボタンを使うと、お使いのブラウザ上で両方の時間を計ることができます。測定された比率もNとともに上昇しますが、これらの数値には達しません。これはメモリアクセスや加算に時間がかかる一方、演算回数にはそれが考慮されていないためです。

大きなNでは、高速な手法の方が精度も高くなる 🖖

精度表は、各演算を単精度（有効数字約7桁）に丸めて両方の手法を実行し、倍精度の答えと比較したものです。1,024点では、直接和が答えの大きさに対して 5.75×10⁻⁷ ずれるのに対し、FFTは 1.23×10⁻⁷ となり、約 4.7 分の1に抑えられています。16点から 1,024 点へと 64 倍に増やすと、直接和の誤差は 7.3 倍になり、FFTの誤差は 1.4 倍になりました。直接和の各出力は丸められた積をN回足し合わせたものなので、丸め誤差が 1,024 回の加算を通じて蓄積されます。一方、FFTの各出力はわずか log₂N = 10 個のバタフライしか通過しません。16点では両者にほとんど差はなく、この入力ではFFTの方が 8.74×10⁻⁸ 対 7.91×10⁻⁸ とわずかに劣ることすらあります。どちらの手法も同じ方法で各回転因子を取得し、倍精度で計算して1回だけ丸めているため、この表は2つのアルゴリズムを比較するものであり、サインの取得方法を比較するものではありません。

全プロセスの詳細解説

ガウスのもう一つの例：ジュノー、36の値 4 ステップ
ガウス手稿の第41項では、10°ごとに取られた36の値からジュノーの中心差を補間し、これを 6 × 6 に分割しています。クーリーとテューキーの計算によれば、直接和にはどれだけの演算が必要で、ガウスの分割によってどれだけ削減でき、さらに細かく分割すればより良くなるでしょうか。まず ガウス、長さ4の周期×3 を押してください。このページでは12の値に対しても同じように計算しています。

この問題をツールで開く
1. 直接和、周波数1つ・サンプル1つにつき1回の演算：
  
  $N^{2} = 3 6^{2} = 1, 296$
2. ガウスの分割、長さ6の周期を6つ：クーリー・テューキーの2段階の計算量は、Nと因数の和の積になる。
  
  $T = N (N_{1} + N_{2}) = 36 \times (6 + 6) = 432$
3. 最後まで分割した場合、36 = 2 × 2 × 3 × 3：
  
  $T = 36 \times (2 + 2 + 3 + 3) = 360$
4. クーリーとテューキーの法則：s, t > 1 となる因数 r = st は、s = t = 2 の場合を除き、常に s + t < r となる：
  
  $s + t < s t (s, t > 1, (s, t) \neq = (2, 2)); 6 = 2 \times 3 : 2 + 3 = 5 < 6$
解答

直接和で 1,296 回、ガウスの 6 × 6 の分割で 432 回、最後まで分割すると 360 回の演算になります。 ガウスの分割の時点で計算量は3分の1になっており、最も細かい分割ではそこからさらに6分の1が削減されます。このページの12の値でも同様に、直接和で 144 回、分割で 84 回になります。4を 2 × 2 に分割してもコストは同じなので、最後まで分割しても 12 × (2 + 2 + 3) = 84 回のままです。

$1, 296 \to 432 \to 360; 1 2^{2} = 144 \to 12 \times (4 + 3) = 12 \times (2 + 2 + 3) = 84$

学習の道すじ

フーリエ変換の実際

参考文献 (3)

Gauss's manuscript on trigonometric interpolation, printed after his death. Article 28 is the page's Gauss mode: twelve declinations of the northern limit of Pallas's zodiac at every 30° of right ascension, split into three periods of four and into four periods of three. Article 27 says the split extends to three or more factors and that practice will teach whoever tries it ("praxis tentantem docebit"); article 41 works Juno's equation of the centre from 36 values with the eccentricity 0.254236. The editor's note calls the text a reworking of a draft apparently begun in October 1805: C. F. Gauss, Theoria interpolationis methodo nova tractata (Nachlass). Werke, Band III, Königliche Gesellschaft der Wissenschaften, Göttingen, 1866; second printing 1876, where the treatise begins on p. 265 and arts. 25–28 are on pp. 303–310.
The paper that showed Gauss's split is the Cooley–Tukey FFT: the change of index n = N₂n₁ + n₂, k = k₁ + N₁k₂ written as Gauss's procedure, his examples (N = 12 as 4 × 3 and 3 × 4, N = 36 as 6 × 6), the dating to October–November 1805 from the Juno eccentricity and other clues, and what the manuscript does not contain: "Gauss did not, however, go on to quantify the computational requirements of his method": 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 (October 1984). doi:10.1109/MASSP.1984.1162257
The paper that made the FFT standard: the split N = r₁r₂ with T = N(r₁ + r₂) operations, its repetition to T = N(r₁ + … + r_m), the radix-2 form with bit-reversed order, and the count of N² operations for the direct sum against fewer than 2N log₂N, where an operation is a complex multiplication followed by a complex addition: J. W. Cooley and J. W. Tukey, An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation 19(90), 297–301 (1965). doi:10.1090/S0025-5718-1965-0178586-1

42Math チーム · 公表された資料と照合しています。訂正のご指摘を歓迎します。