Prime-factor FFT algorithm

The prime-factor algorithm (PFA), also called the Good–Thomas algorithm (1958/1963), is a fast Fourier transform (FFT) algorithm that re-expresses the discrete Fourier transform (DFT) of a size N = N₁N₂ as a two-dimensional N₁×N₂ DFT, but only for the case where N₁ and N₂ are relatively prime. These smaller transforms of size N₁ and N₂ can then be evaluated by applying PFA recursively or by using some other FFT algorithm.

PFA should not be confused with the mixed-radix generalization of the popular Cooley–Tukey algorithm, which also subdivides a DFT of size N = N₁N₂ into smaller transforms of size N₁ and N₂. The latter algorithm can use any factors (not necessarily relatively prime), but it has the disadvantage that it also requires extra multiplications by roots of unity called twiddle factors, in addition to the smaller transforms. On the other hand, PFA has the disadvantages that it only works for relatively prime factors (e.g. it is useless for power-of-two sizes) and that it requires a more complicated re-indexing of the data based on the Chinese remainder theorem (CRT). Note, however, that PFA can be combined with mixed-radix Cooley–Tukey, with the former factorizing N into relatively prime components and the latter handling repeated factors.

PFA is also closely related to the nested Winograd FFT algorithm, where the latter performs the decomposed N₁ by N₂ transform via more sophisticated two-dimensional convolution techniques. Some older papers therefore also call Winograd's algorithm a PFA FFT.

(Although the PFA is distinct from the Cooley–Tukey algorithm, Good's 1958 work on the PFA was cited as inspiration by Cooley and Tukey in their 1965 paper, and there was initially some confusion about whether the two algorithms were different. In fact, it was the only prior FFT work cited by them, as they were not then aware of the earlier research by Gauss and others.)

Algorithm

Recall that the DFT is defined by the formula:

$X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{\frac {2\pi i}{N}}nk}\quad k=0,\,\dots ,\,N-1.$

The PFA involves a re-indexing of the input and output arrays, which when substituted into the DFT formula transforms it into two nested DFTs (a two-dimensional DFT).

Re-indexing

Suppose that $N=N_{1}N_{2}$, where $N_{1}$ and $N_{2}$ are relatively prime, i.e. $\gcd(N_{1},N_{2})=1$. Then the re-indexing is performed using two bijective mappings between $\mathbb {Z} _{N_{1}}\times \mathbb {Z} _{N_{2}}$ and $\mathbb {Z} _{N}$.

The first map

$n=\eta (n_{1},n_{2})=(n_{1}N_{2}+n_{2}N_{1}){\bmod {N}}$

is a bijection called the Ruritanian mapping (also Good's mapping).

Indeed it is a homomorphism $\mathbb {Z} _{N_{1}}\times \mathbb {Z} _{N_{2}}\to \mathbb {Z} _{N}$, because $\forall n_{1},m_{1}\in \mathbb {Z} _{N_{1}}$ and $\forall n_{2},m_{2}\in \mathbb {Z} _{N_{2}}$:

${\begin{aligned}\eta \left((n_{1}+m_{1}){\bmod {N}}_{1},(n_{2}+m_{2}){\bmod {N}}_{2}\right)&=\\&=\left(\left((n_{1}+m_{1}){\bmod {N}}_{1}\right)N_{2}+\left((n_{2}+m_{2}){\bmod {N}}_{2}\right)N_{1}\right){\bmod {N}}\\&=\left(\left((n_{1}+m_{1})N_{2}\right){\bmod {(}}N_{1}N_{2})+\left((n_{2}+m_{2})N_{1}\right){\bmod {(}}N_{1}N_{2})\right){\bmod {N}}\\&=\left((n_{1}+m_{1})N_{2}+(n_{2}+m_{2})N_{1}\right){\bmod {N}}\\&=\left(n_{1}N_{2}+m_{1}N_{2}+n_{2}N_{1}+m_{2}N_{1}\right){\bmod {N}}\\&=\left((n_{1}N_{2}+n_{2}N_{1}){\bmod {N}}+(m_{1}N_{2}+m_{2}N_{1}){\bmod {N}}\right){\bmod {N}}\\&=\left(\eta (n_{1},n_{2})+\eta (m_{1},m_{2})\right){\bmod {N}}.\end{aligned}}$

Therefore, according to the first isomorphism theorem, $\eta$ is an injection to the quotient group $\mathbb {Z} _{N}/\ker(\eta )$. Here the kernel $\ker(\eta )=\{0\}$, because otherwise there would exist a pair $n_{1}\in \mathbb {Z} _{N_{1}}$ and $n_{2}\in \mathbb {Z} _{N_{2}}$, which are not simultaneously zero, such that $n_{1}N_{2}+n_{2}N_{1}=tN_{1}N_{2}$ for some nonzero $t\in \mathbb {N}$. Since $n_{1}N_{2}+n_{2}N_{1}\leq 2N_{1}N_{2}-(N_{1}+N_{2})<2$, the only remaining value of $t$ is 1. In this case $n_{1}=N_{1}-{\frac {n_{2}N_{1}}{N_{2}}}$ would be an integer from $\mathbb {Z} _{N_{1}}$ which is impossible, because for the fraction ${\frac {n_{2}N_{1}}{N_{2}}}$ to yield an integral value, $n_{2}$ must be a multiple of $N_{2}$ (since $\gcd(N_{1},N_{2})=1$). But this would contradict $n_{2}\in \mathbb {Z} _{N_{2}}$. Thus, $\mathbb {Z} _{N}/\ker(\eta )\cong \mathbb {Z} _{N}$, that is $\eta$ injects to $\mathbb {Z} _{N}$. Now since $\operatorname {card} \left(\mathbb {Z} _{N_{1}}\times \mathbb {Z} _{N_{2}}\right)=\operatorname {card} \left(\mathbb {Z} _{N}\right)=N_{1}N_{2}$, $\eta$ is indeed bijective, i.e. for distinct values of the pair $n_{1},n_{2}$ it produces distinct values of $n$ throughout the whole set $\mathbb {Z} _{N}$.

The second map

$k=\kappa (n_{1},n_{2})=(k_{1}N'_{2}N_{2}+k_{2}N'_{1}N_{1}){\bmod {N}}$

is called the CRT mapping. The name refers to the Chinese remainder theorem which provides the bijective mapping $\kappa :\mathbb {Z} _{N_{1}}\times \mathbb {Z} _{N_{2}}\to \mathbb {Z} _{N}$, in which $N'_{1}$ and $N'_{2}$ are any solution to the equation $N'_{1}N_{1}+N'_{2}N_{2}=\gcd(N_{1},N_{2})=1$ (see Bézout's identity).

In order to perform the DFT, one needs to map different pairs of $n_{1}\in \mathbb {Z} _{N_{1}}$ and $n_{2}\in \mathbb {Z} _{N_{2}}$ to distinct values of $n\in \mathbb {Z} _{N}$ and also pairs $k_{1}\in \mathbb {Z} _{N_{1}}$ and $k_{2}\in \mathbb {Z} _{N_{2}}$ to $k\in \mathbb {Z} _{N}$. To do it one can use the Ruritanian mapping to produce indices in the input vector and the CRT mapping to evaluate indices of the output vector or use the two mappings the opposite way.

A great deal of research has been devoted to schemes for evaluating this re-indexing efficiently, ideally in-place, while minimizing the number of costly modulo (remainder) operations (Chan, 1991, and references).

DFT re-expression

The above re-indexing is then substituted into the formula for the DFT, and in particular into the product $nk$ in the exponent. Because $\forall t\in \mathbb {Z} :e^{2\pi it}=1$, this exponent is evaluated modulo $N$. Similarly, $X_{k}$ and $x_{n}$ are implicitly periodic in $N$, so their subscripts are evaluated modulo $N$.

First, substitute the Ruritanian mapping into the formula for DFT:

$X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{\frac {2\pi i}{N}}nk}=\sum _{n=0}^{N_{1}N_{2}-1}x_{n}e^{-{\frac {2\pi i}{N_{1}N_{2}}}nk}=\sum _{n_{1}=0}^{N_{1}-1}\sum _{n_{2}=0}^{N_{2}-1}x_{n_{1}N_{2}+n_{2}N_{1}}e^{-2\pi ik{\frac {n_{1}N_{2}+n_{2}N_{1}}{N_{1}N_{2}}}}=\sum _{n_{1}=0}^{N_{1}-1}e^{-{\frac {2\pi i}{N_{1}}}n_{1}k}\sum _{n_{2}=0}^{N_{2}-1}x_{n_{1}N_{2}+n_{2}N_{1}}e^{-{\frac {2\pi i}{N_{2}}}n_{2}k}$.

Now substitute the CRT mapping in place of $k$ to produce

${\begin{aligned}X_{k_{1}N_{2}N'_{2}+k_{2}N_{1}N'_{1}}&=\sum _{n_{1}=0}^{N_{1}-1}e^{-2\pi in_{1}{\frac {k_{1}N_{2}N'_{2}+k_{2}N_{1}N'_{1}}{N_{1}}}}\sum _{n_{2}=0}^{N_{2}-1}x_{n_{1}N_{2}+n_{2}N_{1}}e^{-2\pi in_{2}{\frac {k_{1}N_{2}N'_{2}+k_{2}N_{1}N'_{1}}{N_{2}}}}\\&=\sum _{n_{1}=0}^{N_{1}-1}e^{-2\pi in_{1}k_{1}{\frac {N_{2}N'_{2}}{N_{1}}}}\sum _{n_{2}=0}^{N_{2}-1}x_{n_{1}N_{2}+n_{2}N_{1}}e^{-2\pi in_{2}k_{2}{\frac {N_{1}N'_{1}}{N_{2}}}}\\&=\sum _{n_{1}=0}^{N_{1}-1}e^{-2\pi in_{1}k_{1}{\frac {1-N_{1}N'_{1}}{N_{1}}}}\sum _{n_{2}=0}^{N_{2}-1}x_{n_{1}N_{2}+n_{2}N_{1}}e^{-2\pi in_{2}k_{2}{\frac {1-N_{2}N'_{2}}{N_{2}}}}\\&=\sum _{n_{1}=0}^{N_{1}-1}e^{-{\frac {2\pi in_{1}k_{1}}{N_{1}}}}\sum _{n_{2}=0}^{N_{2}-1}x_{n_{1}N_{2}+n_{2}N_{1}}e^{-{\frac {2\pi in_{2}k_{2}}{N_{2}}}}.\end{aligned}}$

The inner and outer sums are simply DFTs of size $N_{2}$ and $N_{1}$, respectively.

References

• Good, I. J. (1958). "The interaction algorithm and practical Fourier analysis". Journal of the Royal Statistical Society, Series B. 20 (2): 361–372. JSTOR 2983896. Addendum, ibid. 22 (2), 373-375 (1960) JSTOR 2984108.
• Thomas, L. H. (1963). "Using a computer to solve problems in physics". Applications of Digital Computers. Boston: Ginn.
• Duhamel, P.; Vetterli, M. (1990). "Fast Fourier transforms: a tutorial review and a state of the art". Signal Processing. 19 (4): 259–299. doi:10.1016/0165-1684(90)90158-U.
• Chan, S. C.; Ho, K. L. (1991). "On indexing the prime-factor fast Fourier transform algorithm". IEEE Trans. Circuits and Systems. 38 (8): 951–953. doi:10.1109/31.85638.

See also

• Rader's FFT algorithm
• Bluestein's FFT algorithm