Rader's FFT algorithm: Difference between revisions

Content deleted Content added
m Reverted edits by 203.90.123.194 to last version by Stevenj
changed notation (f,j,k,n) -> (X,k,n,N) for consistency with discrete Fourier transform
Line 7:
Recall that the DFT is defined by the formula
 
:<math> f_jX_k = \sum_{kn=0}^{nN-1} x_kx_n e^{-\frac{2\pi i}{nN} jknk }
\qquad
jk = 0,\dots,nN-1. </math>
 
If ''nN'' is a prime number, then the set of non-zero indices ''kn'' = 1,...,''nN''&ndash;1 forms a [[group (mathematics)|group]] under multiplication [[modular arithmetic|modulo]] ''nN''. One consequence of the [[number theory]] of such groups is that there exists a [[generating set of a group|generator]] of the group, an integer ''g'' such that ''kn'' = ''g''<sup>''q''</sup> (mod ''nN'') for any non-zero index ''kn'' and for a unique ''q'' in 0,...,''nN''&ndash;2 (forming a [[bijection]] from ''q'' to non-zero ''kn''). Similarly ''jk'' = ''g''<sup>&ndash;''p''</sup> (mod ''nN'') for any non-zero index ''jk'' and for a unique ''p'' in 0,...,''nN''&ndash;2, where the negative exponent denotes the multiplicative inverse of ''g''<sup>''p''</sup> modulo ''nN''. That means that we can rewrite the DFT using these new indices ''p'' and ''q'' as:
 
:<math> f_0X_0 = \sum_{kn=0}^{nN-1} x_kx_n,</math>
 
:<math> f_X_{g^{-p}} = x_0 + \sum_{q=0}^{nN-2} x_{g^q} e^{-\frac{2\pi i}{nN} g^{-(p-q)} }
\qquad
p = 0,\dots,nN-2. </math>
 
(Recall that ''x''<sub>''kn''</sub> and ''fX''<sub>''jk''</sub> are implicitly periodic in ''nN'', and also that ''e''<sup>2&pi;i</sup>=1. Thus, all indices and exponents are taken modulo ''nN'' as required by the group arithmetic.)
 
The final summation, above, is precisely a cyclic convolution of the two sequences ''a''<sub>''q''</sub> and ''b''<sub>''q''</sub> of length ''nN''&ndash;1 (''q'' = 0,...,''nN''&ndash;2) defined by:
 
:<math>a_q = x_{g^q}</math>
:<math>b_q = e^{-\frac{2\pi i}{nN} g^{-q} }.</math>
 
===Evaluating the convolution===
 
Since ''nN''&ndash;1 is composite, this convolution can be performed directly via the [[convolution theorem]] and more conventional FFT algorithms. However, that may not be efficient if ''nN''&ndash;1 itself has large prime factors, requiring recursive use of Rader's algorithm. Instead, one can compute a length-(''nN''&ndash;1) cyclic convolution exactly by zero-padding it to a length of at least 2(''nN''&ndash;1)&ndash;1, say to a [[power of two]], which can then be evaluated in O(''nN'' log ''nN'') time without the recursive application of Rader's algorithm.
 
This algorithm, then, requires O(''nN'') additions plus O(''nN'' log ''nN'') time for the convolution. In practice, the O(''nN'') 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''<sub>''kn''</sub> is given by the DC (0th) output of the FFT of ''a''<sub>''q''</sub> plus ''x''<sub>0</sub>, and ''x''<sub>0</sub> 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 ''nN''&ndash;1 to compute the convolution, rather than by zero padding as mentioned above, the efficiency depends strongly upon ''nN'' and the number of times that Rader's algorithm must be applied recursively. The worst case would be if ''nN''&ndash;1 were 2''nN''<sub>2</sub> where ''nN''<sub>2</sub> is prime, with ''nN''<sub>2</sub>&ndash;1 = 2''nN''<sub>3</sub> where ''nN''<sub>3</sub> 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(''nN''<sup>2</sup>) time. Such ''nN''<sub>j</sub> are called [[Sophie Germain prime|Sophie Germain primes]], and such a sequence of them is called a [[Cunningham chain]] of the first kind. The lengths of Cunningham chains, however, are observed to grow more slowly than log<sub>2</sub>(''nN''), so Rader's algorithm applied in this way is probably not [[Big O notation|&Omega;]](''nN''<sup>2</sup>), though it is possibly worse than O(''nN'' log ''nN'') for the worst cases. Fortunately, a guarantee of O(''nN'' log ''nN'') complexity can be achieved by zero padding.
 
==References==