Rader's FFT algorithm: Difference between revisions

Revision as of 00:55, 4 October 2019

Latest revision as of 21:35, 10 December 2024

(6 intermediate revisions by 5 users not shown)

Line 1:

{{Short description|Discrete Fourier transform for prime sizes}}

'''Rader's algorithm''' (1968),<ref>C. M. Rader, "Discrete Fourier transforms when the number of data samples is prime," ''Proc. IEEE'' 56, 1107–1108 (1968).</ref> named for Charles M. Rader of [[MIT Lincoln Laboratory]], is a [[fast Fourier transform]] (FFT) algorithm that computes the [[discrete Fourier transform]] (DFT) of [[prime number|prime]] sizes by re-expressing the DFT as a cyclic [[convolution]] (the other algorithm for FFTs of prime sizes, [[Bluestein's FFT algorithm|Bluestein's algorithm]], also works by rewriting the DFT as a convolution).

Line 8 ⟶ 9:

==Algorithm==

[[File:FFT visual Rader 11.jpg|thumb|Visual representation of a [[DFT matrix]] in Rader's FFT algorithm,. The array consists of colored clocks representrepresenting a DFT matrix of size 11. By permuting rows and columns using(except athe sequencefirst generatedof byeach) theaccording primitiveto rootsequences ofgenerated 11by (whichthe ispowers 2) exceptof the 1stprimitive rowroot andof column11, the original DFT matrix becomes a [[circulant matrix]]. Multiplying a circulantdata matrixsequence towith a datacirculant sequencematrix is equivalent to the [[cyclic convolution]] with the matrix's row vector. This relation is an example of the fact that the [[Multiplicative group | multiplicative group]] is cyclic: <math>(\mathbb Z/p\mathbb Z)^\times \cong C_{p-1}</math>.]]

Begin with the definition of the discrete Fourier transform:

:<math> X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2\pi i}{N} nk }

\qquad

k = 0,\dots,N-1. </math>

If ''N'' is a prime number, then the set of non-zero indices ''<math>n'' =\in{} \{1,...\dots,''N''&ndash;-1\}</math> forms a [[group (mathematics)|group]] under multiplication [[modular arithmetic|modulo]] ''N''. One consequence of the [[number theory]] of such groups is that there exists a [[generating set of a group|generator]] of the group (sometimes called a [[Primitive root modulo n|primitive root]], which can be found quickly by exhaustive search or slightly better algorithms<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>),. This generator is an integer ''g'' such that ''<math>n'' = ''g''<sup>''^q'' \pmod N</supmath> (mod ''N'') for any non-zero index ''n'' and for a unique ''<math>q'' \in{} \{0,...\dots,''N''&ndash;-2\}</math> (forming a [[bijection]] from ''q'' to non-zero ''n''). Similarly, ''<math>k'' = ''g''<sup>&ndash;''^{-p''} \pmod N</supmath> (mod ''N'') for any non-zero index ''k'' and for a unique ''<math>p'' \in{} \{0,...\dots,''N''&ndash;-2\}</math>, where the negative exponent denotes the [[modular multiplicative inverse|multiplicative inverse]] of ''g''<supmath>''g^p'' \mod N</supmath> modulo ''N''. That means that we can rewrite the DFT using these new indices ''p'' and ''q'' as:

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

Line 21 ⟶ 24:

p = 0,\dots,N-2. </math>

(Recall that ''x''_''n'' and ''X''_''k'' are implicitly periodic in ''N'', and also that ''e''<supmath>2πi e^{2\pi i}=1 </supmath>=1. ([[Euler's identity]]). Thus, all indices and exponents are taken modulo ''N'' as required by the group arithmetic.)

The final summation, above, is precisely a cyclic convolution of the two sequences ''a''_''q'' and ''b''_''q'' (of length ''N''&ndash;1, (''because <math>q'' =\in{} \{0,...\dots,''N''&ndash;-2\}</math>) defined by:

:<math>a_q = x_{g^q}</math>

Line 34 ⟶ 37:

This algorithm, then, requires O(''N'') additions plus O(''N'' log ''N'') time for the convolution. In practice, the O(''N'') additions can often be performed by absorbing the additions into the convolution: if the convolution is performed by a pair of FFTs, then the sum of ''x''_''n'' is given by the DC (0th) output of the FFT of ''a''_''q'' plus ''x''₀, and ''x''₀ can be added to all the outputs by adding it to the DC term of the convolution prior to the inverse FFT. Still, this algorithm requires intrinsically more operations than FFTs of nearby composite sizes, and typically takes 3&ndash;10 times as long in practice.

If Rader's algorithm is performed by using FFTs of size ''N''&ndash;1 to compute the convolution, rather than by zero padding as mentioned above, the efficiency depends strongly upon ''N'' and the number of times that Rader's algorithm must be applied recursively. The worst case would be if ''N''&ndash;1 were 2''N''₂ where ''N''₂ is prime, with ''N''₂&ndash;1 = 2''N''₃ where ''N''₃ is prime, and so on. In such cases, supposing that the chain of primes extended all the way down to some bounded value, the recursive application of Rader's algorithm would actually require O(''N''²) time. Such ''N''_j are called [[Sophie Germain prime]]s, and such a sequence of them is called a [[Cunningham chain]] of the first kind. The lengths of Cunningham chainsHowever, however,the arealternative observedof tozero growpadding morecan slowlyalways thanbe log₂(''N''),employed soif Rader's algorithm applied in this way is probably not [[Big O notation|O]](''N''²),&ndash;1 though it is possibly worse than O(''N'' log ''N'') for the worst cases. Fortunately,has a guaranteelarge ofprime O(''N'' log ''N'') complexity can be achieved by zero paddingfactor.

==References==