互质因子算法

互质因子算法（Prime-factor FFT algorithm, PFA），又称为Good-Thomas算法 ，是一种快速傅立叶变换（FFT），把N = N₁N₂大小的离散傅立叶变换重新表示为N₁ * N₂大小的二维离散傅立叶变换，其中N₁与N₂需互质。变成N₁和N₂大小的傅立叶变换后，可以继续递回使用PFA，或用其他快速傅立叶变换算法来计算。

较流行的Cooley-Tukey算法经由mixed-radix一般化后，也是把N = N₁N₂大小的离散傅立叶变换分割为N₁和N₂大小的转换，但和互质因子算法 (PFA)作法并不相同，不应混淆。Cooley-Tukey算法的N₁与N₂不需互质，可以是任何整数。然而有个缺点是比PFA多出一些乘法，和单位根twiddle factors相乘。相对的，PFA的缺点则是N₁与N₂需互质 (例如N 是2次方就不适用)，而且要借由中国剩馀定理来进行较复杂的re-indexing。互质因子算法 (PFA)可以和mixed-radix Cooley-Tukey算法相结合，前者将N 分解为互质的因数，后者则用在重复质因数上。

PFA也与nested Winograd FFT算法密切相关，后者使用更为精巧的二维折积技巧分解成N₁ * N₂的转换。因而一些较古老的论文把Winograd算法称为PFA FFT。

尽管PFA和Cooley-Tukey算法并不相同，但有趣的是Cooley和Tukey在他们1965年发表的有名的论文中，没有发觉到高斯和其他人更早的研究，只引用Good在1958年发表的PFA作为前人的FFT结果。刚开始的时候人们对这两种作法是否不同有点困惑。

算法

离散傅立叶变换（DFT）的定义如下:

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

PFA将输入和输出re-indexing，代入DFT公式后转换成二维DFT。

Re-indexing

设N = N₁N₂，N₁与N₂两者互质，然后把输入n 和输出k 一一对应到

$n=n_{1}N_{2}+n_{2}N_{1}\mod N$

因N₁与N₂ 互质，故根据最大公因数表现定理，对每个n 都存在满足上式的整数n₁与n₂，且在同馀N 之下n₁可以调整至0～N₁ –1之间，n₂可以调整至0～N₂ –1之间。并根据同馀理论易知满足上式且在以上范围内的整数n₁与n₂是唯一的。这称为Ruritanian 映射 (或Good's 映射)，

$k=k_{1}\mod N_{1}$

$k=k_{2}\mod N_{2}$

举例来说:

如果$N=15,N_{1}=5,N_{2}=3,n=0,1,2,...,12,13,14,$对于任一 $n$ 都可以对应到

$n=n_{1}N_{2}+n_{2}N_{1}\mod N,n_{1}=0,1,...,N_{1}-1,n_{2}=0,1,...,N_{2}-1$

$0=0\centerdot N_{2}+0\centerdot N_{1}\mod 15$

$1=2\centerdot N_{2}+2\centerdot N_{1}\mod 15$

$2=4\centerdot N_{2}+1\centerdot N_{1}\mod 15$

$3=1\centerdot N_{2}+0\centerdot N_{1}\mod 15$

$4=3\centerdot N_{2}+2\centerdot N_{1}\mod 15$

$5=0\centerdot N_{2}+1\centerdot N_{1}\mod 15$

$6=2\centerdot N_{2}+0\centerdot N_{1}\mod 15$

$7=4\centerdot N_{2}+2\centerdot N_{1}\mod 15$

$8=1\centerdot N_{2}+1\centerdot N_{1}\mod 15$

$9=3\centerdot N_{2}+0\centerdot N_{1}\mod 15$

$10=0\centerdot N_{2}+2\centerdot N_{1}\mod 15$

$11=2\centerdot N_{2}+1\centerdot N_{1}\mod 15$

$12=4\centerdot N_{2}+0\centerdot N_{1}\mod 15$

$13=1\centerdot N_{2}+2\centerdot N_{1}\mod 15$

$14=3\centerdot N_{2}+1\centerdot N_{1}\mod 15$

因N₁与N₂ 互质，故根据中国剩馀定理，对于每组 ( k₁ , k₂ ) (其中k₁在0～N₁ – 1之间, k₂在0～N₂ – 1之间)，都有存在且唯一的k 在0～N - 1之间且满足上两式。这称为 CRT 映射。 CRT 映射的另一种表示法如下

$k=k_{1}N_{2}^{-1}N_{2}+k_{2}N_{1}^{-1}N_{1}\mod N$

其中N₁^-1表示N₁在模N₂之下的反元素，N₂^-1反之。

( 也可以改成对输入n 用 CRT 映射以及对输出k 用Ruritanian 映射)

对于有效re-indexing (理想上是达到原地)的方法有许多研究，以减少耗费时间的模运算。

DFT re-expression

表示方法一:

将以上的re-indexing代入DFT公式里指数部分的nk 之中，

$e^{-{\frac {2\pi i}{N}}nk}=e^{-{\frac {2\pi i}{N}}(n_{1}N_{2}+n_{2}N_{1})k}=e^{-{\frac {2\pi i}{N_{1}}}kn_{1}}e^{-{\frac {2\pi i}{N_{2}}}kn_{2}}=e^{-{\frac {2\pi i}{N_{1}}}k_{1}n_{1}}e^{-{\frac {2\pi i}{N_{2}}}k_{2}n_{2}}$

( 因为e^2πi = 1，所以两个指数的k 部份可以分别模N₁与N₂ )。剩下的部分变成

$X_{k_{1},k_{2}}=\sum _{n_{1}=0}^{N_{1}-1}\left(\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_{2}}\right)e^{-{\frac {2\pi i}{N_{1}}}n_{1}k_{1}}.$

则内部和外部的总和分别转换成大小为N₂与N₁的DFT。

表示方法二:

如果令 $k=k_{1}N_{2}+k_{2}N_{1}\quad for\quad k=0,1,...,N-1,$

令 $n=((n_{1}N_{2}+n_{2}N_{1}))_{N}$，$(\cdot )_{N}$相当于取 $N$的馀数，$n_{1}=0,\dots ,N_{1}-1$, $n_{2}=0,\dots ,N_{2}-1$

$X[((k_{1}N_{2}+k_{2}N_{1}))_{N}]=\sum _{n=0}^{N-1}x[((n_{1}N_{2}+n_{2}N_{1}))_{N}]e^{-j{\frac {2\pi }{N_{2}N_{1}}}(k_{1}N_{2}+k_{2}N_{1})(n_{1}N_{2}+n_{2}N_{1})}$

$=\sum _{n=0}^{N-1}x[((n_{1}N_{2}+n_{2}N_{1}))_{N}]e^{-j{\frac {2\pi }{N_{2}N_{1}}}(k_{1}n_{1}N_{2}N_{2}+k_{2}n_{2}N_{1}N_{1}+k_{1}n_{2}N_{2}N_{1}+k_{2}n_{1}N_{1}N_{2})}$

$=\sum _{n=0}^{N-1}x[((n_{1}N_{2}+n_{2}N_{1}))_{N}]e^{-j{\frac {2\pi }{N_{1}}}(k_{1}n_{1}N_{2})}e^{-j{\frac {2\pi }{N_{2}}}(k_{2}n_{2}N_{1})}$

$=\sum _{n_{2}=0}^{N_{2}-1}\{\sum _{n_{1}=0}^{N_{1}-1}x[((n_{1}N_{2}+n_{2}N_{1}))_{N}]e^{-j{\frac {2\pi }{N_{1}}}(k_{1}n_{1}N_{2})}\}e^{-j{\frac {2\pi }{N_{2}}}(k_{2}n_{2}N_{1})}.$

对于每一个 $n_{2}$ 都要做一个 $N_{1}$ 点的 $DFT$，而因为 $n_{2}=0,\dots ,N_{2}-1$有 $N_{2}$个，所以需要 $N_{2}$ 个 $N_{1}$ 点 $DFT$,

对于每一组$((k_{1}N_{2}))_{N_{1}}$都要做一个 $N_{2}$ 点的 $DFT$，而因为 $N_{2}$为常数，$k_{1}=0,\dots ,N_{1}-1$有 $N_{1}$ 个，所以需要 $N_{1}$ 个 $N_{2}$ 点 $DFT$，

因此如果要计算复杂度，可以乘法器的数量当作考量,

假设$N_{1}$ 点的 $DFT$ 需要 $M_{1}$个乘法器,

假设$N_{2}$ 点的 $DFT$ 需要 $M_{2}$个乘法器,

则总共需要 $N_{2}M_{1}+N_{1}M_{2}$个乘法器。

范例

以N = 6为例，有两种可能，N₁ = 2, N₂ = 3或N₁ = 3, N₂ = 2。

N₁ = 2, N₂ = 3

N₁ = 3, N₂ = 2

第一种情形所产生的流程图如左图所示。先做2次3点DFT后再做3次2点DFT。

第二种情形所产生的流程图如右图所示。先做3次2点DFT后再做2次3点DFT。

其中2点DFT的部份因构造单纯，皆以交错的蝴蝶图来显示。

可以看出即使在这个简单的例子中，输入和输出的index也都经过有点复杂的重新排列。

与Cooley-Tukey算法的比较

如首段所述，Cooley-Tukey算法和互质因子算法 (PFA)曾被误认为很类似。两者皆有各自优点可适用于不同状况，因此分辨它们的不同是很重要的。在1965年著名的论文中发表的Cooley-Tukey算法，是在DFT的定义

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

中代入n = n₁ + n₂N₁ , k = k₁N₂ + k₂，则

$e^{-{\frac {2\pi i}{N}}nk}=e^{-{\frac {2\pi i}{N}}(n_{1}+n_{2}N_{1})(k_{1}N_{2}+k_{2})}=e^{-{\frac {2\pi i}{N_{1}}}n_{1}k_{1}}e^{-{\frac {2\pi i}{N}}n_{1}k_{2}}e^{-{\frac {2\pi i}{N_{2}}}n_{2}k_{2}}$

$X_{k_{1}N_{2}+k_{2}}=\sum _{n_{1}=0}^{N_{1}-1}\left(\sum _{n_{2}=0}^{N_{2}-1}x_{n_{1}+n_{2}N_{1}}e^{-{\frac {2\pi i}{N_{2}}}n_{2}k_{2}}\right)e^{-{\frac {2\pi i}{N}}n_{1}k_{2}}e^{-{\frac {2\pi i}{N_{1}}}n_{1}k_{1}}$

比PFA多了一些要乘的因子$e^{-{\frac {2\pi i}{N}}n_{1}k_{2}}$ (称为twiddle factors )，但index较为简单，且适用于任何N₁、N₂。在J. Cooley稍后发表的关于FFT历史探讨的论文中使用N = 24点FFT为例，显示两种作法在index结构上的不同。

相关条目

• 快速傅立叶变换
• 中国剩馀定理
• Bézout引理

参考文献

1. P. Duhamel and M. Vetterli, Fast Fourier transforms: a tutorial review and a state of the art, Signal Processing, 1990, 19: 259–299 

外部链接

• fft note by Burrus（页面存档备份，存于互联网档案馆）
• cnx（页面存档备份，存于互联网档案馆）