Editing Fast Fourier transform (section)

==Multidimensional FFTs==
<!-- This section is linked from [[Discrete Fourier transform]] -->

As defined in the [[Discrete Fourier transform#Multidimensional DFT|multidimensional DFT]] article, the multidimensional DFT

:<math>X_\mathbf{k} = \sum_{\mathbf{n}=0}^{\mathbf{N}-1} e^{-2\pi i \mathbf{k} \cdot (\mathbf{n} / \mathbf{N})} x_\mathbf{n}</math>

transforms an array {{math|''x''{{sub|'''n'''}}}} with a {{mvar|d}}-dimensional [[coordinate vector|vector]] of indices <math display="inline">\mathbf{n} = \left(n_1, \ldots, n_d\right)</math> by a set of {{mvar|d}} nested summations (over <math display="inline">n_j = 0 \ldots N_j - 1</math> for each {{mvar|j}}), where the division <math display="inline">\mathbf{n} / \mathbf{N} = \left(n_1/N_1, \ldots, n_d/N_d\right)</math> is performed element-wise.  Equivalently, it is the composition of a sequence of ''d'' sets of one-dimensional DFTs, performed along one dimension at a time (in any order).

This compositional viewpoint immediately provides the simplest and most common multidimensional DFT algorithm, known as the '''row-column''' algorithm (after the two-dimensional case, below). That is, one simply performs a sequence of {{mvar|d}} one-dimensional FFTs (by any of the above algorithms): first you transform along the {{math|''n''{{sub|1}}}} dimension, then along the {{math|''n''{{sub|2}}}} dimension, and so on (actually, any ordering works). This method is easily shown to have the usual <math display="inline">O(n \log n)</math> complexity, where <math display="inline">n = n_1 \cdot n_2 \cdots n_d</math> is the total number of data points transformed. In particular, there are {{math|''n''/''n''{{sub|1}}}} transforms of size {{math|''n''{{sub|1}}}}, etc., so the complexity of the sequence of FFTs is:

:<math>\begin{align}
      & \frac{n}{n_1} O(n_1 \log n_1) + \cdots + \frac{n}{n_d} O(n_d \log n_d) \\[6pt]
  ={} & O\left(n \left[\log n_1 + \cdots + \log n_d\right]\right) = O(n \log n).
\end{align}</math>

In two dimensions, the ''x''<sub>'''k'''</sub> can be viewed as an <math>n_1 \times n_2</math> [[matrix (mathematics)|matrix]], and this algorithm corresponds to first performing the FFT of all the rows (resp. columns), grouping the resulting transformed rows (resp. columns) together as another <math>n_1 \times n_2</math> matrix, and then performing the FFT on each of the columns (resp. rows) of this second matrix, and similarly grouping the results into the final result matrix.

In more than two dimensions, it is often advantageous for [[Cache (computing)|cache]] locality to group the dimensions recursively.  For example, a three-dimensional FFT might first perform two-dimensional FFTs of each planar slice for each fixed ''n''<sub>1</sub>, and then perform the one-dimensional FFTs along the ''n''<sub>1</sub> direction. More generally, an [[asymptotically optimal]] [[cache-oblivious algorithm]] consists of recursively dividing the dimensions into two groups <math display="inline">(n_1, \ldots, n_{d/2})</math> and <math display="inline">(n_{d/2+1}, \ldots, n_d)</math> that are transformed recursively (rounding if {{mvar|d}} is not even) (see Frigo and Johnson, 2005).<ref name="Frigo_Johnson_2005"/> Still, this remains a straightforward variation of the row-column algorithm that ultimately requires only a one-dimensional FFT algorithm as the base case, and still has <math>O(n \log n)</math> complexity.  Yet another variation is to perform matrix [[transpose|transpositions]] in between transforming subsequent dimensions, so that the transforms operate on contiguous data; this is especially important for [[out-of-core]] and [[distributed memory]] situations where accessing non-contiguous data is extremely time-consuming.

There are other multidimensional FFT algorithms that are distinct from the row-column algorithm, although all of them have <math display="inline">O(n \log n)</math> complexity. Perhaps the simplest non-row-column FFT is the [[vector-radix FFT algorithm]], which is a generalization of the ordinary Cooley–Tukey algorithm where one divides the transform dimensions by a vector <math display="inline">\mathbf{r} = \left(r_1, r_2, \ldots, r_d\right)</math> of radices at each step. (This may also have cache benefits.) The simplest case of vector-radix is where all of the radices are equal (e.g. vector-radix-2 divides ''all'' of the dimensions by two), but this is not necessary.  Vector radix with only a single non-unit radix at a time, i.e. <math display="inline">\mathbf{r} = \left(1, \ldots, 1, r, 1, \ldots, 1\right)</math>, is essentially a row-column algorithm. Other, more complicated, methods include polynomial transform algorithms due to Nussbaumer (1977),<ref name="Nussbaumer_1977"/> which view the transform in terms of convolutions and polynomial products. See Duhamel and Vetterli (1990)<ref name="Duhamel_Vetterli_1990"/> for more information and references.