雷德演算法

雷德演算法（英語：Rader's algorithm）是一種於1968年由麻省理工學院林肯實驗室之查爾斯·M·雷德（Charles M. Rader）提出的快速傅立葉轉換演算法^[1]。當訊號的資料點數量為質數時，此演算法可藉由將離散傅立葉轉換重新表示為圓周摺積，快速計算出該訊號之離散傅立葉轉換結果。另一種稱為布魯斯坦演算法的作法也是透過類似的方式將離散傅立葉轉換改寫為摺積完成轉換，且同樣限制訊號長度需為質數。

由於雷德演算法之運作原理只依賴於具有週期性的離散傅立葉轉換核，故它也可直接套用於其它具有類似特性的轉換（惟訊號長度需為質數），例如數論轉換或離散哈特利轉換。

當所欲轉換的訊號資料皆為實數時，可透過重新索引或排列將原轉換拆成兩個長度為一半的實數循環摺積，如此稍微修改後的雷德演算法可進一步使運算時間減半^[2]；另一種快速計算實數訊號之離散傅立葉轉換的方法則是使用離散哈特利轉換^[3]。

薩繆爾·威諾格拉德（英语：Shmuel Winograd）進一步延伸了雷德演算法，使其也可用於長度為質數之次方數 ${\displaystyle p^{m}}$的訊號之離散傅立葉轉換^[4][5]；也因此，雷德演算法亦可被視為威諾格拉德快速傅立葉轉換演算法（亦稱為乘法性傅立葉轉換演算法^[6]）的特例之一，其中後者可用的訊號長度範圍較廣。然而，當訊號長度為合數（如質數之次方數）時，使用庫利－圖基快速傅立葉轉換演算法更加簡單且實作上也較容易，故雷德演算法一般只用於庫利－圖基演算法之遞迴拆解下較大質數之基本情況^[3]。

演算法

將雷德演算法之離散傅里葉轉換矩陣視覺化的結果。圖中上色之時鐘圖示代表了大小為11的矩陣中各元素之相位。除了第一行與第一列外，在使用11之原根（即2）重新排列矩陣後，原始之離散傅里葉轉換矩陣即形成一循环矩阵。將訊號與一循環矩陣相乘即等同於圓周摺積。

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

若${\displaystyle N}$是一個質數，則${\displaystyle n=1,\ldots ,N-1}$之非零索引數集合即形成一整数模N乘法群。運用数论的一個結論，可知此類的群中存在一個生成元（有時亦稱為原根，可由穷举搜索或其它較有效率之演算法快速找到^[7]）——一個整數 ${\displaystyle g,}$使得對於任意的非零索引數 ${\displaystyle n\in \{1,2,\cdots ,N-1\},}$都存在唯一的 ${\displaystyle q\in \{0,1,\cdots ,N-2\}}$令 ${\displaystyle n=g^{q}{\pmod {N}}}$，即形成一${\displaystyle q}$至非零${\displaystyle n}$的對射；同樣地，對於任意的非零索引數 ${\displaystyle k\in \{1,2,\cdots ,N-1\},}$都存在唯一的 ${\displaystyle p\in \{0,1,\cdots ,N-2\}}$令 ${\displaystyle k=g^{-p}{\pmod {N}}}$，其中指數的負號代表的是 ${\displaystyle g^{p}}$之模反元素。這代表所求之離散傅立葉轉換可用新的索引數 ${\displaystyle p}$及 ${\displaystyle q}$改寫如下：

${\displaystyle X_{0}=\sum _{n=0}^{N-1}x_{n},}$

${\displaystyle X_{g^{-p}}=x_{0}+\sum _{q=0}^{N-2}x_{g^{q}}e^{-{\frac {2\pi i}{N}}g^{-(p-q)}}\qquad p=0,\dots ,N-2}$

其中${\displaystyle x_{n}}$和${\displaystyle X_{k}}$皆對${\displaystyle N}$隱含週期性，而${\displaystyle e^{2\pi i}=1}$。因此，所有索引數以及指數皆可依群算術之要求取模${\displaystyle N}$之結果。

上式中最後的加總即為長度${\displaystyle N-1}$（${\displaystyle q=0,\ldots ,N-2}$）之兩數列${\displaystyle a_{q}}$和${\displaystyle b_{q}}$的圓周摺積：

${\displaystyle a_{q}=x_{g^{q}}}$

${\displaystyle b_{q}=e^{-{\frac {2\pi i}{N}}g^{-q}}}$

摺積計算

由於${\displaystyle N-1}$必為合數，上述之圓周摺積可直接由摺積定理以及其它常用之快速傅立葉轉換演算法求得。然而，若 ${\displaystyle N-1}$本身具有較大之質因數，則此作法須遞迴使用雷德演算法，而較不具效率。替代方法之一乃是將原長度為 ${\displaystyle N-1}$之數列補零至長度大於${\displaystyle 2(N-1)-1}$，甚或是2的次方數，便可透過快速演算法在${\displaystyle O(N\log {N})}$的時間內求得而不須遞迴使用雷德演算法。

如此一來，此演算法則需要${\displaystyle O(N)}$之加法運算以及${\displaystyle O(N\log {N})}$之摺積運算。實作上，${\displaystyle O(N)}$之加法常可被包含至後續的摺積運算中：若摺積是由一對快速傅立葉轉換求得，則 ${\displaystyle x_{n}}$的加總即是 ${\displaystyle x_{0}}$以及 ${\displaystyle a_{q}}$離散傅立葉轉換的第0項輸出（即DC項）之和。同時，可在逆轉換前先將 ${\displaystyle x_{0}}$加至摺積之DC項，這樣可使最後的輸出皆已包含${\displaystyle x_{0}}$。不過，此演算所需的運算步驟仍然比其它接近之合數長度的快速傅立葉轉換來得多，實務上耗時是其 3 至 10 倍。

若雷德演算法未使用如上的補零法，而直接透過長度${\displaystyle N-1}$之快速傅立葉轉換計算，則其效率將取決於${\displaystyle N}$值以及雷德演算法本身遞迴呼叫的次數。最壞情況乃是當${\displaystyle N-1=2N_{2}}$且${\displaystyle N_{2}}$為質數，${\displaystyle N_{2}-1=2N_{3}}$且${\displaystyle N_{3}}$為質數，依此類推；在此情況下，若上述的質數鍊一直延伸至某界值，則遞迴雷德演算法之複雜度將為${\displaystyle O(N^{2})}$。諸如此類的${\displaystyle N_{j}}$稱為索菲·熱爾曼質數，而它們所形成的質數數列則稱為第一類康寧漢鍊（英语：Cunningham chain）。然而，康寧漢鍊之長度增長的速度一般而言遠較${\displaystyle \log _{2}{N}}$慢，故如此使用雷德演算法之複雜度應不為${\displaystyle O(N^{2})}$，但在最壞情況下複雜度仍應比${\displaystyle O(N\log {N})}$高。不過，使用上述的補零法，便可達到 ${\displaystyle O(N\log {N})}$之複雜度。

參考資料

1. C. M. Rader, "Discrete Fourier transforms when the number of data samples is prime," Proc. IEEE 56, 1107–1108 (1968).
2. S. Chu and C. Burrus, "A prime factor FTT [sic] algorithm using distributed arithmetic," IEEE Transactions on Acoustics, Speech, and Signal Processing 30 (2), 217–227 (1982).
3. Matteo Frigo and Steven G. Johnson, "The Design and Implementation of FFTW3 （页面存档备份，存于互联网档案馆）," Proceedings of the IEEE 93 (2), 216–231 (2005).
4. S. Winograd, "On Computing the Discrete Fourier Transform", Proc. National Academy of Sciences USA, 73(4), 1005–1006 (1976).
5. S. Winograd, "On Computing the Discrete Fourier Transform", Mathematics of Computation, 32(141), 175–199 (1978).
6. R. Tolimieri, M. An, and C.Lu, Algorithms for Discrete Fourier Transform and Convolution, Springer-Verlag, 2nd ed., 1997.
7. Donald E. Knuth, The Art of Computer Programming, vol. 2: Seminumerical Algorithms, 3rd edition, section 4.5.4, p. 391 (Addison–Wesley, 1998).