互质因子算法(Prime-factor FFT algorithm, PFA),又称为Good-Thomas算法[1]
[2],是一种快速傅立叶变换(FFT),把N = N1N2大小的离散傅立叶变换重新表示为N1 * N2大小的二维离散傅立叶变换,其中N1与N2需互质。变成N1和N2大小的傅立叶变换后,可以继续递回使用PFA,或用其他快速傅立叶变换算法来计算。
较流行的Cooley-Tukey算法经由mixed-radix一般化后,也是把N = N1N2大小的离散傅立叶变换分割为N1和N2大小的转换,但和互质因子算法 (PFA)作法并不相同,不应混淆。Cooley-Tukey算法的N1与N2不需互质,可以是任何整数。然而有个缺点是比PFA多出一些乘法,和单位根twiddle factors相乘。相对的,PFA的缺点则是N1与N2需互质 (例如N 是2次方就不适用),而且要借由中国剩余定理来进行较复杂的re-indexing。互质因子算法 (PFA)可以和mixed-radix Cooley-Tukey算法相结合,前者将N 分解为互质的因数,后者则用在重复质因数上。
PFA也与nested Winograd FFT算法密切相关,后者使用更为精巧的二维折积技巧分解成N1 * N2的转换。因而一些较古老的论文把Winograd算法称为PFA FFT。
尽管PFA和Cooley-Tukey算法并不相同,但有趣的是Cooley和Tukey在他们1965年发表的有名的论文中,没有发觉到高斯和其他人更早的研究,只引用Good在1958年发表的PFA作为前人的FFT结果。刚开始的时候人们对这两种作法是否不同有点困惑。
离散傅立叶变换(DFT)的定义如下:
![{\displaystyle X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{\frac {2\pi i}{N}}nk}\qquad k=0,\dots ,N-1}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f4cf021c3cf4d0d691655ea39d75aa8bfe4408a3)
PFA将输入和输出re-indexing,代入DFT公式后转换成二维DFT。
设N = N1N2,N1与N2两者互质,然后把输入n 和输出k 一一对应到
![{\displaystyle n=n_{1}N_{2}+n_{2}N_{1}\mod N}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a36f2a9d6b866e9f99a0da93c1e4d0ef18d6679f)
因N1与N2 互质,故根据最大公因数表现定理,对每个n 都存在满足上式的整数n1与n2,且在同余N 之下n1可以调整至0~N1 –1之间,n2可以调整至0~N2 –1之间。并根据同余理论易知满足上式且在以上范围内的整数n1与n2是唯一的。这称为Ruritanian 映射 (或Good's 映射),
![{\displaystyle k=k_{1}\mod N_{1}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7fafda94575b1e0ff9bd6d53bfab3718aaeaf168)
![{\displaystyle k=k_{2}\mod N_{2}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9e60696c06586ae0fd288380924dc4ed5907c338)
举例来说:
如果
对于任一
都可以对应到
因N1与N2 互质,故根据中国剩余定理,对于每组 ( k1 , k2 ) (其中k1在0~N1 – 1之间, k2在0~N2 – 1之间),都有存在且唯一的k 在0~N - 1之间且满足上两式。这称为 CRT 映射。
CRT 映射的另一种表示法如下
![{\displaystyle k=k_{1}N_{2}^{-1}N_{2}+k_{2}N_{1}^{-1}N_{1}\mod N}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a1fdde1bcb3196e049c2660431b1602edf657d8d)
其中N1-1表示N1在模N2之下的反元素,N2-1反之。
( 也可以改成对输入n 用 CRT 映射以及对输出k 用Ruritanian 映射)
对于有效re-indexing (理想上是达到原地)的方法有许多研究[3],以减少耗费时间的模运算。
表示方法一:
将以上的re-indexing代入DFT公式里指数部分的nk 之中,
![{\displaystyle 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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3c029599f7949ad7849a2d7197c7924845441491)
( 因为e2πi = 1,所以两个指数的k 部分可以分别模N1与N2 )。剩下的部分变成
![{\displaystyle 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}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fa7272a47064bad08270a8b0ad5926cd7734fccc)
则内部和外部的总和分别转换成大小为N2与N1的DFT。
表示方法二:
如果令
令
,
相当于取
的余数,
,
对于每一个
都要做一个
点的
,而因为
有
个,所以需要
个
点
,
对于每一组
都要做一个
点的
,而因为
为常数,
有
个,所以需要
个
点
,
因此如果要计算复杂度,可以乘法器的数量当作考量,
假设
点的
需要
个乘法器,
假设
点的
需要
个乘法器,
则总共需要
个乘法器。
以N = 6为例,有两种可能,N1 = 2, N2 = 3或N1 = 3, N2 = 2。
N1 = 2, N2 = 3
N1 = 3, N2 = 2
第一种情形所产生的流程图如左图所示。先做2次3点DFT后再做3次2点DFT。
第二种情形所产生的流程图如右图所示。先做3次2点DFT后再做2次3点DFT。
其中2点DFT的部分因构造单纯,皆以交错的蝴蝶图来显示。
可以看出即使在这个简单的例子中,输入和输出的index也都经过有点复杂的重新排列。
如首段所述,Cooley-Tukey算法和互质因子算法 (PFA)曾被误认为很类似。两者皆有各自优点可适用于不同状况,因此分辨它们的不同是很重要的。在1965年著名的论文中发表的Cooley-Tukey算法,是在DFT的定义
![{\displaystyle X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{\frac {2\pi i}{N}}nk}\qquad k=0,\dots ,N-1}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f4cf021c3cf4d0d691655ea39d75aa8bfe4408a3)
中代入n = n1 + n2N1 , k = k1N2 + k2,则
![{\displaystyle 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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/916e079e01c17f7574f24a2baf9f2b035e375874)
![{\displaystyle 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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7e3a83b757b3ae692102fc82e7cc20c5fe380896)
比PFA多了一些要乘的因子
(称为twiddle factors ),但index较为简单,且适用于任何N1、N2。在J. Cooley稍后发表的关于FFT历史探讨的论文[4]中使用N = 24点FFT为例,显示两种作法在index结构上的不同。
- ^ I. J. Good, The interaction algorithm and practical Fourier analysis, J. R. Statist. Soc. B, 1958, 20(2): 361–372
- ^ L. H. Thomas, Using a computer to solve problems in physics, Applications of Digital Computers, 1963
- ^ S. C. Chan and K. L. Ho, On indexing the prime-factor fast Fourier transform algorithm, IEEE Trans. Circuits and Systems, 1991, 38(8): 951–953 .
- ^ J. Cooley, P. Lewis, and P. Welch, Historical notes on the fast Fourier transform, IEEE Transactions on Audio and Electroacoustics, 1967, 15(2): 76–79
- P. Duhamel and M. Vetterli, Fast Fourier transforms: a tutorial review and a state of the art, Signal Processing, 1990, 19: 259–299