Tables for
Volume B
Reciprocal space
Edited by U. Shmueli

International Tables for Crystallography (2010). Vol. B, ch. 1.3, pp. 58-59   | 1 | 2 |

Section Multidimensional Cooley–Tukey factorization

G. Bricognea

aGlobal Phasing Ltd, Sheraton House, Suites 14–16, Castle Park, Cambridge CB3 0AX, England, and LURE, Bâtiment 209D, Université Paris-Sud, 91405 Orsay, France Multidimensional Cooley–Tukey factorization

| top | pdf |

Let us now assume that this decimation can be factored into d successive decimations, i.e. that[{\bf N} = {\bf N}_{1} {\bf N}_{2} \ldots {\bf N}_{d-1} {\bf N}_{d}]and hence[{\bf N}^{T} = {\bf N}_{d}^{T} {\bf N}_{d - 1}^{T} \ldots {\bf N}_{2}^{T} {\bf N}_{1}^{T}.]Then the coset decomposition formulae corresponding to these successive decimations (Section[link]) can be combined as follows:[\eqalign{{\bb Z}^{n} &= \bigcup_{{\bf k}_{1}}\, ({\bf k}_{1} + {\bf N}_{1} {\bb Z}^{n}) \cr &= \bigcup_{{\bf k}_{1}}\, \left\{{\bf k}_{1} + {\bf N}_{1} \left[\bigcup_{{\bf k}_{2}}\, ({\bf k}_{2} + {\bf N}_{2} {\bb Z}^{n})\right]\right\} \cr &= \ldots \cr &= \bigcup_{{\bf k}_{1}} \ldots \bigcup_{{\bf k}_{d}}\, ({\bf k}_{1} + {\bf N}_{1} {\bf k}_{2} + \ldots + {\bf N}_{1} {\bf N}_{2} \times \ldots \times {\bf N}_{d - 1} {\bf k}_{d} + {\bf N} {\bb Z}^{n})}]with [{\bf k}_{j} \in {\bb Z}^{n} / {\bf N}_{j} {\bb Z}^{n}]. Therefore, any [{\bf k} \in {\bb Z} / {\bf N} {\bb Z}^{n}] may be written uniquely as[{\bf k} = {\bf k}_{1} + {\bf N}_{1} {\bf k}_{2} + \ldots + {\bf N}_{1} {\bf N}_{2} \times \ldots \times {\bf N}_{d - 1} {\bf k}_{d}.]Similarly:[\eqalign{{\bb Z}^{n} &= \bigcup_{{\bf k}_{d}^{*}}\, ({\bf k}_{d}^{*} + {\bf N}_{d}^{T} {\bb Z}^{n}) \cr &= \ldots \cr &= \bigcup_{{\bf k}_{d}^{*}} \ldots \bigcup_{{\bf k}_{1}^{*}} \,({\bf k}_{d}^{*} + {\bf N}_{d}^{T} {\bf k}_{d - 1}^{*} + \ldots + {\bf N}_{d}^{T} \times \ldots \times {\bf N}_{2}^{T} {\bf k}_{1}^{*} \cr &\quad + {\bf N}^{T} {\bb Z}^{n})}]so that any [{\bf k}^{*} \in {\bb Z}^{n} / {\bf N}^{T} {\bb Z}^{n}] may be written uniquely as[{\bf k}^{*} = {\bf k}_{d}^{*} + {\bf N}_{d}^{T} {\bf k}_{d - 1}^{*} + \ldots + {\bf N}_{d}^{T} \times \ldots \times {\bf N}_{2}^{T} {\bf k}_{1}^{*}]with [{\bf k}_{j}^{*} \in {\bb Z}^{n} / {\bf N}_{j}^{T} {\bb Z}^{n}]. These decompositions are the vector analogues of the multi-radix number representation systems used in the Cooley–Tukey factorization.

We may then write the definition of [\bar{F} ({\bf N})] with [d = 2] factors as[\eqalign{X^{*} ({\bf k}_{2}^{*} + {\bf N}_{2}^{T} {\bf k}_{1}^{*}) &= {\textstyle\sum\limits_{{\bf k}_{1}}} {\textstyle\sum\limits_{{\bf k}_{2}}}\, X ({\bf k}_{1} + {\bf N}_{1} {\bf k}_{2}) \cr &\quad \times e[({\bf k}_{2}^{*T} + {\bf k}_{1}^{*T}{\bf N}_2) {\bf N}_{2}^{-1} {\bf N}_{1}^{-1} ({\bf k}_{1} + {\bf N}_{1} {\bf k}_{2})].}]The argument of e(–) may be expanded as[{\bf k}_{2}^{*} \cdot ({\bf N}^{-1} {\bf k}_{1}) + {\bf k}_{1}^{*} \cdot ({\bf N}_{1}^{-1} {\bf k}_{1}) + {\bf k}_{2}^{*} \cdot ({\bf N}_{2}^{-1} {\bf k}_{2}) + {\bf k}_{1}^{*} \cdot {\bf k}_{2}.]The first summand may be recognized as a twiddle factor, the second and third as the kernels of [\bar{F} ({\bf N}_{1})] and [\bar{F} ({\bf N}_{2})], respectively, while the fourth is an integer which may be dropped. We are thus led to a `vector-radix' version of the Cooley–Tukey algorithm, in which the successive decimations may be introduced in all n dimensions simultaneously by general integer matrices. The computation may be decomposed into five stages analogous to those of the one-dimensional algorithm of Section[link]:

  • (i) form the [|{\bf N}_{1}|] vectors [{\bf Y}_{{\bf k}_{1}}] of shape [{\bf N}_{2}] by[Y_{{\bf k}_{1}} ({\bf k}_{2}) = X ({\bf k}_{1} + {\bf N}_{1} {\bf k}_{2}),\quad {\bf k}_{1} \in {\bb Z}^{n} / {\bf N}_{1} {\bb Z}^{n},\quad {\bf k}_{2} \in {\bb Z}^{n} / {\bf N}_{2} {\bb Z}^{n}\hbox{\semi}]

  • (ii) calculate the [|{\bf N}_{1}|] transforms [{\bf Y}_{{\bf k}_{1}}^{*}] on [|{\bf N}_{2}|] points:[Y_{{\bf k}_{1}}^{*} ({\bf k}_{2}^{*}) = {\textstyle\sum\limits_{{\bf k}_{2}}}\, e[{\bf k}_{2}^{*} \cdot ({\bf N}_{2}^{-1} {\bf k}_{2})] Y_{{\bf k}_{1}} ({\bf k}_{2}),\quad {\bf k}_{1} \in {\bb Z}^{n} / {\bf N}_{1} {\bb Z}^{n}\hbox{\semi}]

  • (iii) form the [|{\bf N}_{2}|] vectors [{\bf Z}_{{\bf k}_{2}^{*}}] of shape [{\bf N}_{1}] by[\displaylines{Z_{{\bf k}_{2}^{*}} ({\bf k}_{1}) = e[{\bf k}_{2}^{*} \cdot ({\bf N}^{-1} {\bf k}_{1})] Y_{{\bf k}_{1}}^{*} ({\bf k}_{2}^{*}),\quad {\bf k}_{1} \in {\bb Z}^{n} / {\bf N}_{1} {\bb Z}^{n},\cr {\bf k}_{2}^{*} \in {\bb Z}^{n} / {\bf N}_{2}^{T} {\bb Z}^{n}\hbox{\semi}}]

  • (iv) calculate the [|{\bf N}_{2}|] transforms [{\bf Z}_{{\bf k}_{2}^{*}}^{*}] on [|{\bf N}_{1}|] points:[Z_{{\bf k}_{2}^{*}}^{*} ({\bf k}_{1}^{*}) = {\textstyle\sum\limits_{{\bf k}_{1}}}\, e[{\bf k}_{1}^{*} \cdot ({\bf N}_{1}^{-1} {\bf k}_{1})] Z_{{\bf k}_{2}^{*}} ({\bf k}_{1}),\quad {\bf k}_{2}^{*} \in {\bb Z}^{n} / {\bf N}_{2}^{T} {\bb Z}^{n}\hbox{\semi}]

  • (v) collect [X^{*} ({\bf k}_{2}^{*} + {\bf N}_{2}^{T} {\bf k}_{1}^{*})] as [Z_{{\bf k}_{2}^{*}}^{*} ({\bf k}_{1}^{*})].

The initial [|{\bf N}|]-point transform [\bar{F} ({\bf N})] can thus be performed as [|{\bf N}_{1}|] transforms [\bar{F} ({\bf N}_{2})] on [|{\bf N}_{2}|] points, followed by [|{\bf N}_{2}|] transforms [\bar{F} ({\bf N}_{1})] on [|{\bf N}_{1}|] points. This process can be applied successively to all d factors. The same decomposition applies to [F ({\bf N})], up to the complex conjugation of twiddle factors, the normalization factor [1 / |{\bf N}|] being obtained as the product of the factors [1 / |{\bf N}_{j}|] in the successive partial transforms [F ({\bf N}_{j})].

The geometric interpretation of this factorization in terms of partial transforms on translates of sublattices applies in full to this n-dimensional setting; in particular, the twiddle factors are seen to be related to the residual translations which place the sublattices in register within the big lattice. If the intermediate transforms are performed in place, then the quantity[X^{*} ({\bf k}_{d}^{*} + {\bf N}_{d}^{T} {\bf k}_{d - 1}^{*} + \ldots + {\bf N}_{d}^{T} {\bf N}_{d - 1}^{T} \times \ldots \times {\bf N}_{2}^{T} {\bf k}_{1}^{*})]will eventually be found at location[{\bf k}_{1}^{*} + {\bf N}_{1} {\bf k}_{2}^{*} + \ldots + {\bf N}_{1} {\bf N}_{2} \times \ldots \times {\bf N}_{d - 1} {\bf k}_{d}^{*},]so that the final results will have to be unscrambled by a process which may be called `coset reversal', the vector equivalent of digit reversal.

Factoring by 2 in all n dimensions simultaneously, i.e. taking [{\bf N} = 2{\bf M}], leads to `n-dimensional butterflies'. Decimation in time corresponds to the choice [{\bf N}_{1} = 2{\bf I}, {\bf N}_{2} = {\bf M}], so that [{\bf k}_{1} \in {\bb Z}^{n} / 2{\bb Z}^{n}] is an n-dimensional parity class; the calculation then proceeds by[\displaylines{Y_{{\bf k}_{1}} ({\bf k}_{2}) = X ({\bf k}_{1} + 2{\bf k}_{2}),\quad{\bf k}_{1} \in {\bb Z}^{n} / 2{\bb Z}^{n},\quad {\bf k}_{2} \in {\bb Z}^{n} / {\bf M}{\bb Z}^{n}, \cr Y_{{\bf k}_{1}}^{*} = \bar{F} ({\bf M}) [{\bf Y}_{{\bf k}_{1}}],\quad{\bf k}_{1} \in {\bb Z}^{n} / 2{\bb Z}^{n}\hbox{\semi} \cr \eqalign{X^{*} ({\bf k}_{2}^{*} + {\bf M}^{T} {\bf k}_{1}^{*}) &= {\textstyle\sum\limits_{{\bf k}_{1} \in {\bb Z}^{n} / 2{\bb Z}^{n}}} (-1)^{{\bf k}_{1}^{*} \cdot {\bf k}_{1}} \cr &\quad \times e[{\bf k}_{2}^{*} \cdot ({\bf N}^{-1} {\bf k}_{1})] Y_{{\bf k}_{1}}^{*} ({\bf k}_{2}^{*}).}\cr}]Decimation in frequency corresponds to the choice [{\bf N}_{1} = {\bf M}], [{\bf N}_{2} = 2{\bf I}], so that [{\bf k}_{2} \in {\bb Z}^{n} / 2{\bb Z}^{n}] labels `octant' blocks of shape M; the calculation then proceeds through the following steps:[\eqalign{Z_{{\bf k}_{2}^{*}} ({\bf k}_{1}) &= \left[{\textstyle\sum\limits_{{\bf k}_{2} \in {\bb Z}^{n} / 2{\bb Z}^{n}}} (-1)^{{\bf k}_{2}^{*} \cdot {\bf k}_{2}} X ({\bf k}_{1} + {\bf M}{\bf k}_{2})\right] \cr &\quad \times e[{\bf k}_{2}^{*} \cdot ({\bf N}^{-1} {\bf k}_{1})], \cr {\bf Z}_{{\bf k}_{2}^{*}}^{*} &= \bar{F} ({\bf M}) [{\bf Z}_{{\bf k}_{2}^{*}}], \cr X^{*} ({\bf k}_{2}^{*} + 2{\bf k}_{1}^{*}) &= Z_{{\bf k}_{2}^{*}}^{*} ({\bf k}_{1}^{*}),}]i.e. the [2^{n}] parity classes of results, corresponding to the different [{\bf k}_{2}^{*} \in {\bb Z}^{n} / 2{\bb Z}^{n}], are obtained separately. When the dimension n is 2 and the decimating matrix is diagonal, this analysis reduces to the `vector radix FFT' algorithms proposed by Rivard (1977)[link] and Harris et al. (1977)[link]. These lead to substantial reductions in the number M of multiplications compared to the row–column method: M is reduced to [3M/4] by simultaneous [2 \times 2] factoring, and to [15M/32] by simultaneous [4 \times 4] factoring.

The use of a nondiagonal decimating matrix may bring savings in computing time if the spectrum of the band-limited function under study is of such a shape as to pack more compactly in a nonrectangular than in a rectangular lattice (Mersereau, 1979[link]). If, for instance, the support K of the spectrum Φ is contained in a sphere, then a decimation matrix producing a close packing of these spheres will yield an aliasing-free DFT algorithm with fewer sample points than the standard algorithm using a rectangular lattice.


Harris, D. B., McClellan, J. H., Chan, D. S. K. & Schuessler, H. W. (1977). Vector radix fast Fourier transform. Rec. 1977 IEEE Internat. Conf. Acoust. Speech Signal Process. pp. 548–551.
Mersereau, R. M. (1979). The processing of hexagonally sampled two-dimensional signals. Proc. IEEE, 67, 930–949.
Rivard, G. E. (1977). Direct fast Fourier transform of bivariate functions. IEEE Trans. Acoust. Speech Signal Process. 25, 250–252.

to end of page
to top of page