International
Tables for Crystallography Volume B Reciprocal space Edited by U. Shmueli © International Union of Crystallography 2006 |
International Tables for Crystallography (2006). Vol. B. ch. 1.3, pp. 53-54
|
The previous two algorithms essentially reduce the calculation of the DFT on N points for N composite to the calculation of smaller DFTs on prime numbers of points, the latter remaining irreducible. However, Rader (1968) showed that the p-point DFT for p an odd prime can itself be factored by invoking some extra arithmetic structure present in .
The ring has the property that its non-zero elements, called units, form a multiplicative group . In particular, all units have a unique multiplicative inverse in , i.e. a unit such that . This endows with the structure of a finite field.
Furthermore, is a cyclic group, i.e. consists of the successive powers of a generator g called a primitive root mod p (such a g may not be unique, but it always exists). For instance, for , is generated by , whose successive powers mod 7 are: [see Apostol (1976), Chapter 10].
The basis of Rader's algorithm is to bring to light a hidden regularity in the matrix by permuting the basis vectors and of as follows: where g is a primitive root mod p.
With respect to these new bases, the matrix representing will have the following elements: Thus the `core' of matrix , of size , formed by the elements with two non-zero indices, has a so-called skew-circulant structure because element depends only on . Simplification may now occur because multiplication by is closely related to a cyclic convolution. Introducing the notation we may write the relation in the permuted bases as where Z is defined by , .
Thus may be obtained by cyclic convolution of C and Z, which may for instance be calculated by where × denotes the component-wise multiplication of vectors. Since p is odd, is always divisible by 2 and may even be highly composite. In that case, factoring by means of the Cooley–Tukey or Good methods leads to an algorithm of complexity p log p rather than for . An added bonus is that, because , the elements of can be shown to be either purely real or purely imaginary, which halves the number of real multiplications involved.
This idea was extended by Winograd (1976, 1978) to the treatment of prime powers , using the cyclic structure of the multiplicative group of units . The latter consists of all those elements of which are not divisible by p, and thus has elements. It is cyclic, and there exist primitive roots g modulo such that The elements divisible by p, which are divisors of zero, have to be treated separately just as 0 had to be treated separately for .
When , then with . The results are p-decimated, hence can be obtained via the -point DFT of the -periodized data Y: with
When , then we may write where contains the contributions from and those from . By a converse of the previous calculation, arises from p-decimated data Z, hence is the -periodization of the -point DFT of these data: with (the -periodicity follows implicity from the fact that the transform on the right-hand side is independent of ).
Finally, the contribution from all may be calculated by reindexing by the powers of a primitive root g modulo , i.e. by writing then carrying out the multiplication by the skew-circulant matrix core as a convolution.
Thus the DFT of size may be reduced to two DFTs of size (dealing, respectively, with p-decimated results and p-decimated data) and a convolution of size . The latter may be `diagonalized' into a multiplication by purely real or purely imaginary numbers (because ) by two DFTs, whose factoring in turn leads to DFTs of size and . This method, applied recursively, allows the complete decomposition of the DFT on points into arbitrarily small DFTs.
When , the same method can be applied, except for a slight modification in the calculation of . There is no primitive root modulo for : the group is the direct product of two cyclic groups, the first (of order 2) generated by −1, the second (of order ) generated by 3 or 5. One then uses a representation and the reindexed core matrix gives rise to a two-dimensional convolution. The latter may be carried out by means of two 2D DFTs on points.
References
Apostol, T. M. (1976). Introduction to analytic number theory. New York: Springer-Verlag.Google ScholarRader, C. M. (1968). Discrete Fourier transforms when the number of data samples is prime. Proc. IEEE, 56, 1107–1108.Google Scholar
Winograd, S. (1976). On computing the discrete Fourier transform. Proc. Natl Acad. Sci. USA, 73, 1005–1006.Google Scholar
Winograd, S. (1978). On computing the discrete Fourier transform. Math. Comput. 32, 175–199.Google Scholar