查看“︁雷德演算法”︁的源代码

{{noteTA
|G1=Signals and Systems
}}
'''雷德演算法'''（{{lang-en|Rader's algorithm}}）是一種於1968年由[[麻省理工學院]][[林肯實驗室]]之查爾斯·M·雷德（Charles M. Rader）提出的[[快速傅里叶变换|快速傅立葉轉換]]演算法<ref>C. M. Rader, "Discrete Fourier transforms when the number of data samples is prime," ''Proc. IEEE'' 56, 1107–1108 (1968).</ref>。當訊號的資料點數量為[[質數]]時，此演算法可藉由將[[離散傅立葉變換|離散傅立葉轉換]]重新表示為[[圓周摺積]]，快速計算出該訊號之離散傅立葉轉換結果。另一種稱為[[Chirp-Z轉換#布魯斯坦演算法|布魯斯坦演算法]]的作法也是透過類似的方式將離散傅立葉轉換改寫為摺積完成轉換，且同樣限制訊號長度需為質數。

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

當所欲轉換的訊號資料皆為實數時，可透過重新索引或排列將原轉換拆成兩個長度為一半的實數循環摺積，如此稍微修改後的雷德演算法可進一步使運算時間減半<ref>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&#x2013;227 (1982).</ref>；另一種快速計算實數訊號之離散傅立葉轉換的方法則是使用[[離散哈特利轉換]]<ref name="Frigo05">Matteo Frigo and Steven G. Johnson, "[http://fftw.org/fftw-paper-ieee.pdf The Design and Implementation of FFTW3] {{Wayback|url=http://fftw.org/fftw-paper-ieee.pdf |date=20190716053010 }}," ''Proceedings of the IEEE'' '''93''' (2), 216–231 (2005).</ref>。

{{Tsl|en|Shmuel Winograd|薩繆爾·威諾格拉德}}進一步延伸了雷德演算法，使其也可用於長度為質數之次方數 <math>p^m</math>的訊號之離散傅立葉轉換<ref>S. Winograd, "On Computing the Discrete Fourier Transform", ''Proc. National Academy of Sciences USA'', '''73'''(4), 1005&ndash;1006 (1976).</ref><ref>S. Winograd, "On Computing the Discrete Fourier Transform", ''Mathematics of Computation'', '''32'''(141), 175&ndash;199 (1978).</ref>；也因此，雷德演算法亦可被視為[[威諾格拉德快速傅立葉變換演算法|威諾格拉德快速傅立葉轉換演算法]]（亦稱為乘法性傅立葉轉換演算法<ref>R. Tolimieri, M. An, and C.Lu, ''Algorithms for Discrete Fourier Transform and Convolution'', Springer-Verlag, 2nd ed., 1997.</ref>）的特例之一，其中後者可用的訊號長度範圍較廣。然而，當訊號長度為[[合數]]（如質數之次方數）時，使用[[库利－图基快速傅里叶变换算法|庫利－圖基快速傅立葉轉換演算法]]更加簡單且實作上也較容易，故雷德演算法一般只用於庫利－圖基演算法之[[遞迴]]拆解下較大質數之基本情況<ref name=Frigo05/>。

==演算法==
[[File:FFT visual Rader 11.jpg|thumb|將雷德演算法之[[離散傅里葉變換矩陣|離散傅里葉轉換矩陣]]視覺化的結果。圖中上色之時鐘圖示代表了大小為11的矩陣中各元素之相位。除了第一行與第一列外，在使用11之[[原根]]（即2）重新排列矩陣後，原始之離散傅里葉轉換矩陣即形成一[[循环矩阵]]。將訊號與一循環矩陣相乘即等同於[[圓周摺積]]。]]
:<math> X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N} nk }
\qquad
k = 0,\dots,N-1 </math>

若<math> N </math>是一個[[質數]]，則<math> n = 1,\ldots,N-1 </math>之非零索引數集合即形成一[[整数模n乘法群|整数模''N''乘法群]]。運用[[数论]]的一個結論，可知此類的[[群]]中存在一個[[群的生成集合|生成元]]（有時亦稱為[[原根]]，可由[[穷举搜索]]或其它較有效率之演算法快速找到<ref>Donald E. Knuth, ''The Art of Computer Programming, vol. 2: Seminumerical Algorithms'', 3rd edition, section 4.5.4, p. 391 (Addison–Wesley, 1998).</ref>）——一個整數 <math>g,</math>使得對於任意的非零索引數 <math>n \in \{1, 2, \cdots, N-1\},</math>都存在唯一的 <math>q \in \{0, 1, \cdots, N-2\}</math>令 <math>n=g^q\pmod N</math>，即形成一<math>q</math>至非零<math>n</math>的[[對射]]；同樣地，對於任意的非零索引數 <math>k \in \{1, 2, \cdots, N-1\},</math>都存在唯一的 <math>p \in \{0, 1, \cdots, N-2\}</math>令 <math>k=g^{-p}\pmod N</math>，其中指數的負號代表的是 <math>g^p</math>之[[模反元素]]。這代表所求之離散傅立葉轉換可用新的索引數 <math>p</math>及 <math>q</math>改寫如下：

:<math> X_0 =  \sum_{n=0}^{N-1} x_n,</math>

:<math> 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 </math>

其中<math> x_n </math>和<math> X_k </math>皆對<math> N </math>隱含週期性，而<math> e^{2\pi i}=1 </math>。因此，所有索引數以及指數皆可依群算術之要求取模<math> N </math>之結果。

上式中最後的加總即為長度<math> N-1 </math>（<math> q=0,\ldots,N-2 </math>）之兩數列<math> a_q </math>和<math> b_q </math>的[[圓周摺積]]：

:<math>a_q = x_{g^q}</math>
:<math>b_q = e^{-\frac{2\pi i}{N} g^{-q} }</math>

===摺積計算===

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

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

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

==參考資料==
{{reflist}}

{{Authority control}}
[[Category:快速傅立葉轉換演算法]]