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 , 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:

${\displaystyle 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 ${\displaystyle N=N_{1}N_{2}}$, where ${\displaystyle N_{1}}$ and ${\displaystyle N_{2}}$ are relatively prime, i.e. ${\displaystyle \gcd(N_{1},N_{2})=1}$. Then the re-indexing is performed using two bijective mappings between ${\displaystyle \mathbb {Z} _{N_{1}}\times \mathbb {Z} _{N_{2}}}$ and ${\displaystyle \mathbb {Z} _{N}}$.

The first map

${\displaystyle 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 ${\displaystyle \mathbb {Z} _{N_{1}}\times \mathbb {Z} _{N_{2}}\to \mathbb {Z} _{N}}$, because ${\displaystyle \forall n_{1},m_{1}\in \mathbb {Z} _{N_{1}}}$ and ${\displaystyle \forall n_{2},m_{2}\in \mathbb {Z} _{N_{2}}}$:

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

The second map

${\displaystyle k=\kappa (k_{1},k_{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 ${\displaystyle \kappa :\mathbb {Z} _{N_{1}}\times \mathbb {Z} _{N_{2}}\to \mathbb {Z} _{N}}$, in which ${\displaystyle N'_{1}}$ and ${\displaystyle N'_{2}}$ are any solution to the linear diophantine equation equation ${\displaystyle 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 ${\displaystyle n_{1}\in \mathbb {Z} _{N_{1}}}$ and ${\displaystyle n_{2}\in \mathbb {Z} _{N_{2}}}$ to distinct values of ${\displaystyle n\in \mathbb {Z} _{N}}$ and also pairs ${\displaystyle k_{1}\in \mathbb {Z} _{N_{1}}}$ and ${\displaystyle k_{2}\in \mathbb {Z} _{N_{2}}}$ to ${\displaystyle 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 ${\displaystyle nk}$ in the exponent. Because ${\displaystyle \forall t\in \mathbb {Z} :e^{2\pi it}=1}$, this exponent is evaluated modulo ${\displaystyle N}$. Similarly, ${\displaystyle X_{k}}$ and ${\displaystyle x_{n}}$ are implicitly periodic in ${\displaystyle N}$, so their subscripts are evaluated modulo ${\displaystyle N}$.

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

${\displaystyle 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 ${\displaystyle k}$ to produce

${\displaystyle {\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}}}$

Likewise, substitution of the CRT mapping in place of ${\displaystyle n}$ and Ruritanian mapping in place of ${\displaystyle k}$ yields

${\displaystyle X_{k_{1}N_{2}+k_{2}N_{1}}=\sum _{n_{1}=0}^{N_{1}-1}e^{-{\frac {2\pi ik_{1}n_{1}}{N_{1}}}}\sum _{n_{2}=0}^{N_{2}-1}x_{n_{1}N_{2}N'_{2}+n_{2}N_{1}N'_{1}}e^{-{\frac {2\pi ik_{2}n_{2}}{N_{2}}}}}$.

In both cases the inner and outer sums are simply DFTs of size ${\displaystyle N_{2}}$ and ${\displaystyle 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