Cooley–Tukey FFTアルゴリズム

Cooley –Tukeyアルゴリズムは、JW CooleyとJohn Tukeyにちなんで名付けられた、最も一般的な高速フーリエ変換（FFT）アルゴリズムです。任意の合成サイズの離散フーリエ変換（DFT）を、 N _{1 個}のより小さなサイズN ₂のDFTで再帰的に表現することで、高度に合成されたN（滑らかな数）の場合、計算時間をO( N log N )に短縮します。このアルゴリズムの重要性から、特定の変種や実装スタイルは、以下に示すように、それぞれ独自の名前で知られています。 ${\displaystyle N=N_{1}N_{2}}$

Cooley-TukeyアルゴリズムはDFTをより小さなDFTに分割するため、他のDFTアルゴリズムと任意に組み合わせることができます。例えば、RaderアルゴリズムやBluesteinアルゴリズムは、Cooley-Tukeyでは分解できない大きな素因数を扱うために使用できます。また、素因数アルゴリズムは、互いに素な因数を分離する際に効率を高めるために利用できます。

このアルゴリズムとその再帰的応用は、カール・フリードリヒ・ガウスによって発明されました。160年後、クーリーとテューキーが独立して再発見し、普及させました。

歴史

このアルゴリズムは、その再帰的応用を含め、1805年頃にカール・フリードリヒ・ガウスによって発明されました。ガウスはこれを小惑星パラスとジュノの軌道の補間に使用しましたが、彼の研究は広く認知されませんでした（ガウスの死後に新ラテン語で出版されたため）。^{[ 1 ] [ 2 ]}しかし、ガウスは漸近的な計算時間を分析しませんでした。19世紀から20世紀初頭にかけて、様々な限定形式が何度か再発見されました。^{[ 2 ]} FFTは、IBMのジェームズ・クーリーとプリンストンのジョン・テューキーが1965年にこのアルゴリズムを再発明し、^{[ 2 ]}コンピュータ上で簡単に実行する方法を説明する論文を発表したことで人気を博しました。^{[ 3 ]}

伝えられるところによると、テューキーはケネディ大統領の科学諮問委員会の会議中にこのアイデアを思いついた。この会議では、国外に設置した地震計を使ってソ連での核兵器実験を検知する方法について議論していた。これらのセンサーは地震の時系列を生成する。しかし、このデータの分析には、センサーの数と時間の長さから、DFTを計算するための高速アルゴリズムが必要になる。この作業は、ソ連の施設を訪問することなく違反を検知できるように、提案された核実験禁止の批准に非常に重要だった。^{[ 4 ] [ 5 ]} 会議のもう一人の参加者であるIBMのリチャード・ガーウィンはこの方法の可能性を認識し、テューキーをクーリーに紹介した。しかし、ガーウィンはクーリーに本来の目的を知らせないようにした。その代わりに、クーリーには、これはヘリウム3の3次元結晶におけるスピン配向の周期性を決定するために必要であると伝えた。その後、Cooley 氏と Tukey 氏は共同論文を発表し、最大 300  kHzのレートでサンプリングできるアナログ - デジタル コンバーターが同時に開発されたため、急速に広く採用されるようになりました。

Gauss が同じアルゴリズム (ただし漸近的コストは分析していない) を記述していたという事実は、Cooley と Tukey の 1965 年の論文の数年後まで認識されませんでした。^{[ 2 ]}彼らの論文では、現在では素因数 FFT アルゴリズム(PFA)と呼ばれているIJ Good の研究のみがインスピレーションとして引用されていました。 ^{[ 3 ]} Good のアルゴリズムは当初 Cooley–Tukey アルゴリズムと同等であると考えられていましたが、PFA はまったく異なるアルゴリズムであることがすぐに認識されました ( Cooley–Tukey が任意の合成サイズをサポートしているのとは異なり、互いに素な因数を持つサイズに対してのみ機能し、中国剰余定理に依存しています)。^{[ 6 ]}

基数2のDITケース

基数2の時間間引き（DIT）FFTは、Cooley-Tukeyアルゴリズムの最も単純かつ最も一般的な形式ですが、高度に最適化されたCooley-Tukey実装では、通常、以下に説明するように、他の形式のアルゴリズムが使用されます。基数2のDITは、サイズNのDFTを、各再帰ステージでサイズN /2の2つのインターリーブDFT（「基数2」と呼ばれる）に分割します。

離散フーリエ変換 (DFT) は次の式で定義されます。

${\displaystyle X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{\frac {2\pi i}{N}}nk},}$

ここでは0 から までの整数です。 ${\displaystyle n}$ ${\displaystyle N-1}$

基数2のDITは、まず偶数インデックスの入力 と奇数インデックスの入力のDFTを計算し、次にこれら2つの結果を組み合わせてシーケンス全体のDFTを生成します。この考え方は再帰的に実行することで、全体の実行時間をO( N log N )に短縮できます。この簡略化された形式では、Nは2のべき乗であると仮定しています。サンプルポイントの数Nは通常、アプリケーションによって自由に選択できるため（例えば、サンプルレートやウィンドウの変更、ゼロパディングなど）、これは多くの場合重要な制約ではありません。 ${\displaystyle (x_{2m}=x_{0},x_{2},\ldots,x_{N-2})}$ ${\displaystyle (x_{2m+1}=x_{1},x_{3},\ldots ,x_{N-1})}$

基数2のDITアルゴリズムは、関数のDFTを偶数インデックスの合計と奇数インデックスの合計の2つの部分に再配置します。 ${\displaystyle x_{n}}$ ${\displaystyle n={2m}}$ ${\displaystyle n={2m+1}}$

${\displaystyle X_{k}=\sum _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N}}(2m)k}+\sum _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N}}(2m+1)k}}$

下の式に示すように、2番目の和から共通乗数を因数分解することができます。すると、2つの和は関数の偶数インデックス部分のDFTと奇数インデックス部分のDFTであることが明らかです。E番目の偶数インデックス入力のDFTを、 O番目の奇数インデックス入力のDFTを と表すと、以下の式が得られます 。${\displaystyle e^{-{\frac {2\pi i}{N}}k}}$ ${\displaystyle x_{2m}}$ ${\displaystyle x_{2m+1}}$ ${\displaystyle x_{n}}$ ${\displaystyle x_{2m}}$ ${\displaystyle E_{k}}$ ${\displaystyle x_{2m+1}}$ ${\displaystyle O_{k}}$

${\displaystyle X_{k}=\underbrace {\sum \limits _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N/2}}mk}} _{{\text{x_{n}の偶数インデックス部分のDFT}{}+e^{-{\frac {2\pi i}{N}}k}\underbrace {\sum \limits _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N/2}}mk}} _{{\text{x_{n}の奇数インデックス部分のDFT=E_{k}+e^{-{\frac {2\pi i}{N}}k}O_{k}\qquad {\text{k=0,\dots ,{\frac {N}{2}}-1.}$

等式は に対して成立することに注意してください。しかし、肝心なのは、 とは に対してのみこのように計算されるということです。複素指数関数 の周期性により、はおよびからも得られます。 ${\displaystyle k=0,\dots ,N-1}$ ${\displaystyle E_{k}}$ ${\displaystyle O_{k}}$ ${\displaystyle k=0,\dots ,{\frac {N}{2}}-1}$ ${\displaystyle X_{k+{\frac {N}{2}}}}$ ${\displaystyle E_{k}}$ ${\displaystyle O_{k}}$

${\displaystyle {\begin{aligned}X_{k+{\frac {N}{2}}}&=\sum \limits _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N/2}}m(k+{\frac {N}{2}})}+e^{-{\frac {2\pi i}{N}}(k+{\frac {N}{2}})}\sum _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N/2}}m(k+{\frac {N}{2}})}\\&=\sum _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N/2}}mk}e^{-2\pi mi}+e^{-{\frac {2\pi i}{N}}k}e^{-\pi i}\sum _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N/2}}mk}e^{-2\pi mi}\\&=\sum _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N/2}}mk}-e^{-{\frac {2\pi i}{N}}k}\sum _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N/2}}mk}\\&=E_{k}-e^{-{\frac {2\pi i}{N}}k}O_{k}\end{aligned}}}$

と を次のように書き直すことができます。 ${\displaystyle X_{k}}$ ${\displaystyle X_{k+{\frac {N}{2}}}}$

${\displaystyle {\begin{aligned}X_{k}&=E_{k}+e^{-{\frac {2\pi i}{N}}{k}}O_{k}\\X_{k+{\frac {N}{2}}}&=E_{k}-e^{-{\frac {2\pi i}{N}}{k}}O_{k}\end{aligned}}}$

この結果は、長さNのDFTをサイズN /2の2つのDFTで再帰的に表現したもので、基数2のDIT高速フーリエ変換の中核を成しています。このアルゴリズムは、中間計算の結果を再利用して複数のDFT出力を計算することで高速化を図っています。最終出力は、とを+/-で組み合わせることで得られることに注意してください。これは単にサイズ2のDFT（この文脈ではバタフライと呼ばれることもあります）です。これを以下のより大きな基数に一般化すると、サイズ2のDFTはより大きなDFT（それ自体はFFTで評価できます）に置き換えられます。 ${\displaystyle E_{k}}$ ${\displaystyle O_{k}\exp(-2\pi ik/N)}$

N =8のデータ フロー ダイアグラム: 時間デシメーションの基数 2 の FFT は、長さN のDFT を 2 つの長さN /2 の DFT に分割し、その後に「バタフライ」操作と呼ばれるN /2 のサイズ 2 の DFTで構成される結合ステージが続きます(データ フロー ダイアグラムの形状からこのように呼ばれています)。

このプロセスは、分割統治アルゴリズムの一般的な手法の例です。ただし、従来の多くの実装では、明示的な再帰は回避され、代わりに幅優先方式で計算ツリーを走査します。

上記のサイズN のDFT を 2 つのサイズN /2 の DFTとして再表現したものは、ダニエルソン–ランチョスの補題と呼ばれることがあります。これは、この恒等式が 1942 年にこの 2 人の著者によって指摘されたためです^{[ 7 ]} (ルンゲの1903 年の研究^{[ 2 ]}の影響を受けて)。彼らは補題を「後方」再帰的に適用し、変換スペクトルが収束するまで DFT サイズを繰り返し2 倍にしました(ただし、彼らは達成した線型的な(つまり、次数N  log  N ) 漸近的複雑性を認識していなかったようです)。ダニエルソン – ランチョスの作業は、機械式または電子式のコンピュータが広く普及する前のことであり、手計算(おそらく加算機などの機械的な補助を使用)が必要でした。彼らは、3 ～ 5 桁の有効桁数の実数入力で動作するサイズ 64 の DFT の計算時間が 140 分であると報告しました。 CooleyとTukeyの1965年の論文では、IBM 7094（おそらく36ビット単精度、約8桁）でサイズ2048の複素DFTを実行した場合、実行時間は0.02分であると報告されています。 ^{[ 3 ]} 時間を演算回数で再スケールすると、これは約80万倍の高速化に相当します。（手計算にかかる時間を比較すると、サイズ64で140分は、浮動小数点演算あたり平均最大16秒に相当し、そのうち約20%は乗算です。）

擬似コード

擬似コードでは、以下の手順を書くことができます。^{[ 8 ]}

X 0,..., N −1 ← ditfft2 ( x , N , s ): (x 0 , x s , x 2 s , ..., x ( N -1) s ) の DFT : N = 1 の場合、X 0 ← x 0自明なサイズ 1 DFT 基本ケース、それ以外の場合、X 0,..., N /2−1 ← ditfft2 ( x , N /2, 2 s ) (x 0 , x 2 s , x 4 s , ..., x ( N -2) s ) の DFT X N /2,..., N −1 ← ditfft2 ( x +s, N /2, 2 s ) k = 0 から (N/2)-1までの (x s , x s +2 s , x s +4 s , ..., x ( N -1) s ) の DFT は、2 つの半分の DFT を結合します: p ← X k q ← exp(−2π i / N k ) X k + N /2 X k ← p + q X k + N /2 ← p − q 終了 終了if

ここで、ditfft2( x , N ,1) は、基数 2 の DIT FFT によってX =DFT( x )をアウトオブプレースで計算します。ここで、 Nは 2 の整数乗であり、s =1 は入力x配列のストライドです。 x + s は、 x _sで始まる配列を示します。

(結果はX内で正しい順序になっており、それ以上のビット反転の置換は必要ありません。よく言われる個別のビット反転ステージの必要性は、以下で説明する特定のインプレース アルゴリズムの場合にのみ発生します。)

高性能FFTの実装では、この単純な疑似コードに比べて、アルゴリズムの実装に多くの変更が加えられています。例えば、再帰のオーバーヘッドを償却するためにN = 1よりも大きな基本ケースを使用したり、回転因子を 事前に計算したり、キャッシュ上の理由からより大きな基数を使用したりできます。これらとその他の最適化を組み合わせることで、パフォーマンスを1桁以上向上させることができます。^{[ 8 ]} （多くの教科書的な実装では、深さ優先再帰は非再帰的な幅優先アプローチに置き換えられていますが、深さ優先再帰の方がメモリ局所性が高いと主張されています。^{[ 8 ] [ 9 ]}）これらのアイデアのいくつかについては、以下でさらに詳しく説明します。 ${\displaystyle \exp[-2\pi ik/N]}$

アイデア

一般的な因数分解におけるCooley-Tukey FFTの基本ステップは、1次元DFTを2次元DFTのようなものとして再解釈することと見ることができます。長さN = N ₁ N _{2の1次元入力配列は}、列優先順序で格納されたN ₁ × N ₂の2次元行列として再解釈されます。まず、N ₂方向（非連続方向）に沿って小さな1次元DFTを実行し、次に位相係数（回転係数）を乗算し、最後にN ₁方向に沿って1次元DFTを実行します。転置ステップは、ここに示すように途中で実行することも、最初または最後に実行することもできます。これは、小さな変換に対して再帰的に実行されます。

より一般的には、クーリー・テューキーアルゴリズムは、複合サイズN = N ₁ N ₂のDFTを再帰的に次のように再表現する：^{[ 10 ]}

1. サイズN _2のN ₁ DFTを実行します。
2. 複素数根(回転因子と呼ばれることが多い) を掛けます。
3. サイズN _1のN ₂ DFTを実行します。

通常、N ₁またはN ₂は小さな因数（必ずしも素数ではない）で、基数（再帰のステージ間で異なる場合があります）と呼ばれます。N ₁が基数の場合、時間デシメーション（DIT）アルゴリズムと呼ばれます。一方、N ₂が基数の場合、周波数デシメーション（DIF、Sande–Tukey アルゴリズムとも呼ばれます）と呼ばれます。上記のバージョンは基数 2 の DIT アルゴリズムでした。最終的な式では、奇数変換を乗算する位相が回転因子であり、偶数と奇数の変換の +/- 組み合わせ（バタフライ）がサイズ 2 の DFT です。（基数の小さな DFT はバタフライと呼ばれることもあります。これは、基数 2 の場合の データフロー ダイアグラムの形状に由来しています。）

バリエーション

Cooley-Tukey アルゴリズムには他にも多くのバリエーションがあります。

• 混合基数実装では、2 に加えてさまざまな (通常は小さな) 因数を持つ複合サイズを処理し、通常は (常にではありませんが) 再帰の素数ベースのケースに O( N ^{2 ) アルゴリズムを採用します (} RaderやBluesteinのアルゴリズムなど、素数ベースのケースにN  log  Nアルゴリズムを採用することもできます)。
• 分割基数は、基数2の最初の変換では回転係数を必要としないという事実を利用して、基数2と基数4を結合し、2のべき乗のサイズで長い間知られている最も低い算術演算回数を達成しました。 ^{[ 10 ]}ただし、最近のバリエーションではさらに低い回数を達成しています。^{[ 11 ] [ 12 ]}（現在のコンピュータでは、パフォーマンスは厳密な演算回数よりもキャッシュとCPUパイプラインの考慮によって決まります。十分に最適化されたFFT実装では、より大きな基数や、かなりのサイズのハードコードされた基本ケース変換が使用されることがよくあります。）^{[ 13 ]}

Cooley–Tukey アルゴリズムを別の角度から見ると、サイズN の1 次元 DFT をN ₁ x N ₂の2 次元 DFT (および回転子) として再表現したもので、出力行列は転置されます。基数 2 アルゴリズムの場合、これらすべての転置の最終結果は、入力 (DIF) または出力 (DIT) インデックスのビット反転に相当します。小さな基数を使用する代わりに、およそ √ N の基数と明示的な入出力行列転置を使用する場合、それは4 段階 FFTアルゴリズム (または転置の回数によっては6 段階) と呼ばれ、当初はメモリの局所性を向上させるために提案されました^{[ 14 ] [ 15 ]}たとえば、キャッシュ最適化やアウトオブコア操作のために提案され、後にキャッシュを無視する最適なアルゴリズムであることが示されました^{[ 16 ]}

一般的なクーリー・テューキー分解は、添え字kとnをそれぞれとに書き換える。ここで添え字k _aとn _aは0からN _a -1（aは1または2）までの範囲である。つまり、入力（n）と出力（k）をそれぞれ列優先と行優先のN ₁ × N ₂の2次元配列に再添え字付けする。これらの添え字付けの違いは、前述のように転置である。この再添え字付けをnkのDFT式に代入すると、交差項は消滅し（その指数は1）、残りの項は次のようになる 。${\displaystyle k=N_{2}k_{1}+k_{2}}$ ${\displaystyle n=N_{1}n_{2}+n_{1}}$ ${\displaystyle N_{1}n_{2}N_{2}k_{1}}$

${\displaystyle X_{N_{2}k_{1}+k_{2}}=\sum _{n_{1}=0}^{N_{1}-1}\sum _{n_{2}=0}^{N_{2}-1}x_{N_{1}n_{2}+n_{1}}e^{-{\frac {2\pi i}{N_{1}N_{2}}}\cdot (N_{1}n_{2}+n_{1})\cdot (N_{2}k_{1}+k_{2})}}$

${\displaystyle =\sum _{n_{1}=0}^{N_{1}-1}\left[e^{-{\frac {2\pi i}{N_{1}N_{2}}}n_{1}k_{2}}\right]\left(\sum _{n_{2}=0}^{N_{2}-1}x_{N_{1}n_{2}+n_{1}}e^{-{\frac {2\pi i}{N_{2}}}n_{2}k_{2}}\right)e^{-{\frac {2\pi i}{N_{1}}}n_{1}k_{1}}}$

${\displaystyle =\sum _{n_{1}=0}^{N_{1}-1}\left(\sum _{n_{2}=0}^{N_{2}-1}x_{N_{1}n_{2}+n_{1}}e^{-{\frac {2\pi i}{N_{2}}}n_{2}k_{2}}\right)e^{-{\frac {2\pi i}{N_{1}N_{2}}}n_{1}(N_{2}k_{1}+k_{2})}}$ 。

ここで、各内側の合計はサイズN ₂の DFT であり、各外側の合計はサイズN ₁の DFT であり、[...] 括弧内の項は回転因子です。

^{Cooley と Tukey [ 3 ]}および Gauss (基数 3 と基数 6 の手順の例を示した)の両方によって示されているように、任意の基数r (および混合基数) を使用できます。 ^{[ 2 ]} Cooley と Tukey は当初、基数バタフライには O( r ² ) の作業が必要であると想定し、基数rの複雑性をO( r ² N / r  log _r N ) = O( N  log ₂ ( N )  r /log ₂ r ) と計算しました。 rの値が2 から 12 までの整数の場合のr /log ₂ rの値を計算すると、最適な基数は 3 ( eに最も近い整数で、 r /log ₂ rを最小化する)であることがわかります。^{[ 3 ] [ 17 ]} しかし、この分析は誤りであった。基数バタフライもDFTであり、FFTアルゴリズムを使用してO( r log r )の演算で実行できるため、基数rは実際には計算量O( r  log( r )  N / r  log _r N )で打ち消され、最適なrはより複雑な考慮によって決定される。実際には、例えば現代のプロセッサの多数のプロセッサレジスタを効果的に活用するためには、かなり大きなr （32または64）が重要であり、^{[ 13 ]}また、無制限の基数r = √ NでもO( N  log  N )の計算量を達成し、前述のように大きなNに対して理論的および実用的な利点がある。 ^{[ 14 ] [ 15 ] [ 16 ]} 

データの並べ替え、ビット反転、インプレースアルゴリズム

上記のDFTの抽象的なCooley-Tukey分解は、何らかの形でアルゴリズムのすべての実装に適用されますが、FFTの各段階におけるデータの順序付けとアクセスの手法には、はるかに多様なものがあります。特に興味深いのは、O(1)の補助記憶容量のみを使用して、入力データを出力データで上書きする インプレースアルゴリズムを考案するという問題です。

最もよく知られている並べ替え手法では、インプレース radix-2 アルゴリズムで明示的なビット反転が使用されます。ビット 反転とは、インデックスnのデータが、数字b ₄ b ₃ b ₂ b ₁ b ₀ (たとえば、N =32 入力の場合は 5 桁) で 2 進数で記述され、数字が反転されたb ₀ b ₁ b ₂ b ₃ b ₄でインデックスに転送される置換です。上記に示したような、出力が入力にインプレースで上書きされる radix-2 DIT アルゴリズムの最終段階について考えます。と がサイズ 2 の DFT と組み合わされると、これら 2 つの値は出力によって上書きされます。ただし、2 つの出力値は、最上位ビットb ₄ ( N =32 の場合)に対応する出力配列の前半と後半に入る必要があります。一方、2 つの入力 とは、最下位ビットb ₀に対応する偶数要素と奇数要素にインターリーブされます。したがって、出力を正しい場所にするためには、b _{0 を}b ₄の代わりに使用し、インデックスはb ₀ b ₄ b ₃ b ₂ b ₁になります。次の再帰ステージでは、これらの最下位 4 ビットはb ₁ b ₄ b ₃ b ₂になります。基数 2 DIT アルゴリズムの再帰ステージをすべて含める場合、すべてのビットを逆にする必要があり、したがって、順序どおりの出力を得るためには、入力をビット反転で前処理 (または出力を後処理) する必要があります (サイズN /2 の各サブ変換が連続するデータに対して動作する場合、DIT入力はビット反転によって前処理されます)。同様に、すべての手順を逆の順序で実行すると、後処理 (または前処理) でビット反転を行う基数 2 DIF アルゴリズムが得られます。 ${\displaystyle E_{k}}$ ${\displaystyle O_{k}}$ ${\displaystyle E_{k}}$ ${\displaystyle O_{k}}$

このアルゴリズムで使用される対数(log) は、底が 2 の対数です。

以下はビット反転置換法を用いて実装された反復基数2FFTアルゴリズムの疑似コードである。^{[ 18 ]}

アルゴリズムiterative-fftの入力: n 個の複素数値の配列a (n は 2 の累乗)。 出力: a の DFT である配列A。 ビット逆コピー(a, A) n ← a .length s = 1からlog( n )まで、m ← 2 s ω m ← exp(−2π i / m ) k = 0からn -1まで、mによって、 ω ← 1 j = 0からm / 2 – 1まで、t ← ω A [ k + j + m /2] u ← A [ k + j ] A [ k + j ] ← u + t A [ k + j + m /2] ← u – t ω ← ω ω mAを返す

ビット逆コピー手順は次のように実装できます。

アルゴリズムbit-reverse-copy( a , A )の入力: n 個の複素数値の配列a (n は 2 の累乗)。 出力:サイズnの配列A。n ← a .length 、k = 0からn – 1 までA [rev(k)] := a [k] を実行する

あるいは、一部のアプリケーション (畳み込みなど) はビット反転されたデータでも同様に機能するため、ビット反転を行わずに順方向変換、処理、逆変換をすべて実行して、自然な順序で最終結果を生成できます。

しかし、多くの FFT ユーザーは自然順序出力を好み、ビット反転を明示的に行う段階は計算時間に無視できない影響を与える可能性があります。^{[ 13 ]}ビット反転は O( N ) 時間で実行でき、多くの研究の対象となってきましたが。^{[ 19 ] [ 20 ] [ 21 ]}また、基数 2 の場合は置換がビット反転ですが、混合基数の場合はより一般的には任意の（混合基数）桁反転となり、置換アルゴリズムの実装がより複雑になります。 さらに、多くのハードウェア アーキテクチャでは、FFT アルゴリズムの中間段階を並べ替えて、連続する（または少なくともより局所的な）データ要素に対して動作するようにすることが望ましいです。 これらの目的のために、Cooley–Tukey アルゴリズム用に、中間段階での別個のビット反転や追加の置換を必要とせずに済む代替実装方式がいくつか考案されています。

アウトオブプレースの場合、つまり出力配列が入力配列と異なる場合、あるいは等価的に同じサイズの補助配列が利用できる場合、問題は大幅に単純化されます。ストックハム自動ソートアルゴリズム^{[ 22 ] [ 23 ]}は、FFT の各ステージをアウトオブプレースで実行します。通常は 2 つの配列間で書き込みを行い、ステージごとにインデックスの 1 つの「桁」を転置します。このアルゴリズムは特にSIMDアーキテクチャで人気があります。 ^{[ 23 ] [ 24 ]}さらに大きな SIMD の利点（連続アクセスの増加）がピースアルゴリズム^{[ 25 ]} に提案されています。このアルゴリズムもステージごとにアウトオブプレースで並べ替えますが、この方法では別個のビット/桁反転と O(NlogN深さ優先を使用して、クーリー–テューキー分解の定義を直接適用することもできます。これは、別個の置換ステップなしで自然な順序のアウトオブプレース出力を生成し（上記の疑似コードのように）、階層型メモリキャッシュを無視できる局所性の利点。 ^{[ 9 ] [ 13 ] [ 26 ]}

A typical strategy for in-place algorithms without auxiliary storage and without separate digit-reversal passes involves small matrix transpositions (which swap individual pairs of digits) at intermediate stages, which can be combined with the radix butterflies to reduce the number of passes over the data.^{[13][27][28][29][30]}

References

1. Gauss, Carl Friedrich (1866). "Theoria interpolationis methodo nova tractata" [Theory regarding a new method of interpolation]. Nachlass (Unpublished manuscript). Werke (in Latin and German). Vol. 3. Göttingen, Germany: Königlichen Gesellschaft der Wissenschaften zu Göttingen. pp. 265–303.
2. Heideman, M. T., D. H. Johnson, and C. S. Burrus, "Gauss and the history of the fast Fourier transform," IEEE ASSP Magazine, 1, (4), 14–21 (1984)
3. Cooley, James W.; Tukey, John W. (1965). "An algorithm for the machine calculation of complex Fourier series". Math. Comput.19 (90): 297–301. doi:10.2307/2003354. JSTOR 2003354.
4. Cooley, James W.; Lewis, Peter A. W.; Welch, Peter D. (1967). "Historical notes on the fast Fourier transform"(PDF). IEEE Transactions on Audio and Electroacoustics. 15 (2): 76–79. CiteSeerX 10.1.1.467.7209. doi:10.1109/tau.1967.1161903.
5. Rockmore, Daniel N., Comput. Sci. Eng.2 (1), 60 (2000). The FFT — an algorithm the whole family can use Special issue on "top ten algorithms of the century "Barry A. Cipra. "The Best of the 20th Century: Editors Name Top 10 Algorithms"(PDF). SIAM News. 33 (4). Archived from the original(PDF) on 2009-04-07. Retrieved 2009-03-31.
6. James W. Cooley, Peter A. W. Lewis, and Peter W. Welch, "Historical notes on the fast Fourier transform," Proc. IEEE, vol. 55 (no. 10), p. 1675–1677 (1967).
7. Danielson, GC、C. Lanczos、「実用的なフーリエ解析の改良と液体からのX線散乱への応用」、 J. Franklin Inst. 233、365–380、435–452 (1942)。
8. S. G. Johnson と M. Frigo、「FFT の実際の実装」、『高速フーリエ変換』（CS Burrus 編）、第 11 章、ライス大学、ヒューストン、テキサス州：Connexions、2008 年 9 月。
9. シングルトン、リチャード・C. (1967). 「高速フーリエ変換の計算について」 . Commun. ACM . 10 (10): 647– 654. doi : 10.1145/363717.363771 . S2CID 6287781 . 
10. Duhamel, P., M. Vetterli, "高速フーリエ変換：チュートリアルレビューと最新技術" Signal Processing 19 , 259–299 (1990)
11. ^kの実数FFTと畳み込みへの新しい行列アプローチ」、 Computing 80、23–45（2007）。
12. Johnson, SG、M. Frigo、「より少ない算術演算を備えた修正分割基数FFT」、 IEEE Trans. Signal Process. 55 (1)、111–119 (2007)。
13. Frigo, M.; Johnson, SG (2005). 「FFTW3の設計と実装」(PDF) . Proceedings of the IEEE . 93 (2): 216– 231. Bibcode : 2005IEEEP..93..216F . CiteSeerX 10.1.1.66.3097 . doi : 10.1109/JPROC.2004.840301 . S2CID 6644892 .  
14. Gentleman WM、G. Sande、「高速フーリエ変換 - 楽しみと利益のために」、Proc. AFIPS 29、563–578 (1966)。
15. Bailey, David H.、「外部メモリまたは階層メモリにおけるFFT」、J. Supercomputing 4 (1), 23–35 (1990)
16. M. Frigo, CE Leiserson, H. Prokop, S. Ramachandran. キャッシュ忘却アルゴリズム. Proceedings of the 40th IEEE Symposium on Foundations of Computer Science (FOCS 99), p.285-297. 1999. IEEE の拡張概要、Citeseer に掲載。
17. Cooley, JW, P. Lewis, P. Welch, "高速フーリエ変換とその応用", IEEE Trans on Education 12 , 1, 28–34 (1969)
18. コーメン、トーマス H.チャールズ・ライザーソン。ロナルド・リベスト。スタイン、クリフォード (2009)。アルゴリズム入門(第 3 版)。マサチューセッツ州ケンブリッジ：MIT Press。ページ 915–918。ISBN​​ 978-0-262-03384-8。
19. Karp, Alan H. (1996). 「ユニプロセッサにおけるビット反転」. SIAM Review . 38 (1): 1– 26. CiteSeerX 10.1.1.24.2913 . doi : 10.1137/1038001 . JSTOR 2132972 .  
20. カーター, ラリー; ガトリン, カン・スー (1998). 「最適なビット反転順列プログラムに向けて」. Proceedings 39th Annual Symposium on Foundations of Computer Science (Cat. No.98CB36280) . pp.  544– 553. CiteSeerX 10.1.1.46.9319 . doi : 10.1109/SFCS.1998.743505 . ISBN  978-0-8186-9172-0. S2CID  14307262 .
21. Rubio, M.; Gómez, P.; Drouiche, K. (2002). 「新しい超高速ビット反転アルゴリズム」. International Journal of Adaptive Control and Signal Processing . 16 (10): 703– 707. doi : 10.1002/acs.718 . S2CID 62201722 . 
22. もともとはWT Cochranら著「 What is the fast Fourier transform?」、 Proc. IEEE vol. 55, 1664–1674 (1967)でStockhamに帰属していた
23. P. N. Swarztrauber、「ベクトルコンピュータのためのFFTアルゴリズム」、並列コンピューティングvol.1、45–63（1984）。
24. Swarztrauber, PN (1982). 「FFTのベクトル化」 . Rodrigue, G. (編).並列計算. ニューヨーク: Academic Press. pp.  51–83 . ISBN 978-0-12-592101-5。
25. Pease, MC (1968). 「並列処理のための高速フーリエ変換の適応」 . J. ACM . 15 (2): 252– 264. doi : 10.1145/321450.321457 . S2CID 14610645 . 
26. フリーゴ、マッテオ;ジョンソン、スティーブン G. 「FFTW」。Cooley-Tukeyアルゴリズムを使用して、任意のサイズの1次元以上の離散フーリエ変換を計算するための無料（GPL ）Cライブラリ
27. Johnson, HW; Burrus, CS (1984). 「インプレース・インオーダー基数2 FFT」. Proc. ICASSP : 28A.2.1–28A.2.4.
28. Temperton, C. (1991). 「自己ソート型インプレース高速フーリエ変換」. SIAM Journal on Scientific and Statistical Computing . 12 (4): 808– 823. doi : 10.1137/0912043 .
29. Qian, Z.; Lu, C.; An, M.; Tolimieri, R. (1994). 「最小作業スペースを備えた自己ソート型インプレースFFTアルゴリズム」. IEEE Transactions on Signal Processing . 52 (10): 2835– 2836. Bibcode : 1994ITSP...42.2835Q . doi : 10.1109/78.324749 .
30. Hegland, M. (1994). 「ベクトル処理と並列処理に適した自己ソート型インプレース高速フーリエ変換アルゴリズム」. Numerische Mathematik . 68 (4): 507– 547. CiteSeerX 10.1.1.54.5659 . doi : 10.1007/s002110050074 . S2CID 121258187 .  

外部リンク

• 「高速フーリエ変換 - FFT」。Cooley -Tukey法。記事。10. C++によるシンプルで教育的な基数2のアルゴリズム
• 「KISSFFT」 . GitHub . 2022年2月11日. C言語によるシンプルな混合基数Cooley-Tukey実装
• GitHubのDsplib
• 「時間FFTアルゴリズムにおける基数2のデシメーション」 。2017年10月31日時点のオリジナルよりアーカイブ。「Алг​​оритм БПФ по основанию два с прореживанием по времени」（ロシア語）。
• 「周波数FFTアルゴリズムにおける基数2のデシメーション」 。2017年11月14日時点のオリジナルよりアーカイブ。「Алг​​оритм БПФ по основанию два с прореживанием по частоте」（ロシア語）。