クーリー・テューキーFFTアルゴリズム

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

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

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

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

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

ガウスが同じアルゴリズムを記述していたこと（ただし漸近的なコストは分析していなかった）は、クーリーとテューキーの 1965 年の論文から数年後まで認識されなかった。^{[ 2 ]}彼らの論文では、現在プライムファクター FFT アルゴリズム(PFA)と呼ばれるものに関するIJ グッド の研究のみをインスピレーションとして引用していた。 ^{[ 3 ]}グッドのアルゴリズムは当初クーリー-テューキー アルゴリズムと同等であると考えられていたが、PFA は全く異なるアルゴリズムであることがすぐに認識された (クーリー-テューキーのように任意の合成サイズをサポートするのではなく、相対的に素因数を持つサイズに対してのみ機能し、中国剰余定理に依存している)。^{[ 6 ]}

基数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},}$

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

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

基数2のDITアルゴリズムは関数のDFTを並べ替えます${\displaystyle x_{n}}$2つの部分に分けます。偶数インデックスの合計${\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}}$

共通乗数を因数分解することができる${\displaystyle e^{-{\frac {2\pi i}{N}}k}}$2番目の和から、以下の式に示すように、2つの和は偶数インデックス部分のDFTであることが明らかです。${\displaystyle x_{2m}}$奇数インデックス部分のDFT${\displaystyle x_{2m+1}}$関数の${\displaystyle x_{n}}$偶数インデックスの入力のDFTを表す。${\displaystyle x_{2m}}$による${\displaystyle E_{k}}$そして、 Oddインデックス入力のDFT${\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{ for }}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{整列}}}$

書き換えることができる${\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 出力を計算することで高速化を実現しています。最終出力は、次の +/− の組み合わせによって得られることに注意してください。${\displaystyle E_{k}}$そして${\displaystyle O_{k}\exp(-2\pi ik/N)}$これは単にサイズ2のDFT（この文脈ではバタフライと呼ばれることもあります）です。これを以下でより大きな基数に一般化すると、サイズ2のDFTはより大きなDFTに置き換えられます（これはFFTで評価できます）。

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

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

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 ) (x s , x s +2 s , x s +4 s , ..., x ( N -1) s ) のDFT、 k = 0 から (N/2)-1 まで、 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 end for end if

より一般的には、Cooley–Tukey アルゴリズムは、合成サイズN = N ₁ N ₂の DFT を再帰的に次のように再表現します。^{[ 10 ]}

1. サイズN2_のN1_個のDFTを実行します。
2. 複素数の1の根（しばしば「回転因子」と呼ばれる）を掛けます。
3. サイズN1の_N2DFTを実行します。

通常、N ₁またはN _2は、基数と呼ばれる小さな因子（必ずしも素数である必要はない）であり、再帰の各段階で異なる場合があります。N ₁が基数の場合、これは時間間引き（DIT）アルゴリズムと呼ばれ、 N ₂が基数の場合、これは周波数間引き（DIF、サンデ・テューキーアルゴリズムとも呼ばれる）アルゴリズムと呼ばれます。上記で示したバージョンは基数 2 の DIT アルゴリズムでした。最終式では、奇数変換に乗算される位相はツイドル因子であり、偶数変換と奇数変換の +/- 組み合わせ（バタフライ）はサイズ 2 の DFT です。（基数の小さな DFT は、基数 2 の場合のデータフロー図の形状から、バタフライと呼ばれることもあります。）

クーリー・テューキーアルゴリズムには、他にも多くのバリエーションが存在する。

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

Cooley–Tukey アルゴリズムを別の視点から見ると、これはサイズN の1 次元 DFT を、出力行列が転置されたN ₁ × N ₂の 2 次元 DFT (およびツイドル) として再表現するものです。基数 2 のアルゴリズムの場合、これらの転置の最終的な結果は、入力 (DIF) または出力 (DIT) インデックスのビット反転に対応します。小さな基数を使用する代わりに、おおよそ√ Nの基数と明示的な入力/出力行列の転置を使用する場合、それは4 ステップ FFTアルゴリズム (転置の数に応じて6 ステップ) と呼ばれ、最初はメモリ局所性を改善するために提案され、^{例えばキャッシュ}最適化やアウトオブコア操作に^使用され、後に最適なキャッシュ非依存アルゴリズムであることが示されました。^{[ 14 ] [ 15 ]}

一般的なクーリー・テューキー因数分解では、インデックスkとn を次のように書き換えます。${\displaystyle k=N_{2}k_{1}+k_{2}}$そして${\displaystyle n=N_{1}n_{2}+n_{1}}$それぞれ、インデックスk _aとn _aは 0.. N _a -1 ( aが 1 または 2 の場合) からなります。つまり、入力 ( n ) と出力 ( k ) を、それぞれ列優先と行優先のN ₁ × N ₂の2 次元配列として再インデックスします。これらのインデックスの違いは、前述のように転置です。この再インデックスをnkの DFT 式に代入すると、${\displaystyle N_{1}n_{2}N_{2}k_{1}}$交差項は消滅し（その指数は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 であり、[...]括弧内の項はツイドル因子です。

任意の基数r (および混合基数) を使用できることは、Cooley と Tukey ^{[ 3 ]}および Gauss (基数 3 および基数 6 のステップの例を示した) によって示されました。^{[ 2 ]} Cooley と Tukey は当初、基数バタフライには O( r ^{2 ) の作業が必要であると想定し、基数}rの複雑さをO( r ² N / r log _r N ) = O( N log ₂ ( N ) r /log ₂ r ) と計算しました。r の整数値 2 から 12 までのr / log ₂ rの値を計算すると、最適な基数は 3 ( r /log ₂ rを最小化するeに最も近い整数) であることがわかります。^{[ 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)の補助記憶領域のみを使用して、入力データを出力データで上書きするインプレースアルゴリズムを考案する問題です。

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

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

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

アルゴリズムiterative-fftの入力: n 個の複素数値の配列a (n は 2 のべき乗)。 出力:配列A (a の DFT)。 ビット反転コピー(a, A) n ← a.length s = 1から log( n )まで繰り返すm ← 2s ω 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 for k = 0 to n – 1 do A [rev(k)] := a [ k]

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

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

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

補助ストレージや個別の桁反転パスを使用しないインプレースアルゴリズムの典型的な戦略は、中間段階で小さな行列転置（個々の桁のペアを交換する）を行うことであり、これは基数バタフライと組み合わせることでデータに対するパスの数を減らすことができる。^{[ 13 ] [ 27 ] [ 28 ] [ 29 ] [ 30 ]}

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