一个用于机器计算复数傅里叶变化的算法

一个用于机器计算复数傅里叶变化的算法
By James W. Cooley and John W. Tukey
一个有效地计算一个2m析因试验的交互算法被Yates发现，也因此与他的的名字广为流传。

更一般3m的算法被Box等人[1]给出。

Good[2]总结了这些方法，并给出了优雅的实用算法来处理一类傅里叶变化。

在这些算法最普适的形式下，Good的方法可以实用地解决一定的问题：在使用过程中，人们总是必须将一个N阶向量乘以N * N的矩阵，而且这个矩阵必须分解成m的稀疏矩阵（其中，m 与logN成正比）。

这个结果的算术复杂度为O（NlogN）而优于O（N2）。

这些方法的本质将被引用在这篇文章中来进行复数傅里叶变化（在这些方法中，当数据点的个数是，或能够被选择成一个高度合成数时，他们是有用的），但是我们在这里提出的算法将会衍生出一个相当不一样的形式。

我们需要注意的是对于N的选择--因为我们可以发现，如果我们使用二进制计算机来计算N=2m的傅里叶变化时，这个算法会有很大的好处。

此外，我们还可以得到，仅仅使用这N 元数组的存储空间，就可以进行完整的计算。

我们考虑如下的复数傅里叶变化
在这个变化中，傅里叶系数是复数形式，而且W是首位的第N个单位根；
若直接利用（1）式进行计算将会需要N2次运算（在这篇文章中，“运算”意味着一次复乘和一次复加）。

而我们在这里提出来的算法基于复数傅里叶变化幅值矩阵的迭代，而且仅需要少于2Nlog2N次运算，且不需要多于数组A的数据储存空间。

为了导出这个算法，我们假设N是一个合数，比如说N=r1*r2。

然后我们让（1）式中的指数表示成
然后，我们变形为
因为
这时，对于k1而言的内和，仅仅取决于j0与k0，而且可以被定义成新的数组
结果则可以表示成
此时有N个元素在数组A1中，每一个要求r1次运算，所以得到A1需要一共计算Nr1次。

相似的，我们需要计算Nr2次从A1中再得到X。

因此，这两步的算法--（6）式与（7）式给出的--一共需要
我们可以清楚的发现，这个算法存在着很好的后续性的应用，从(6)式开始，给出了一个m步的算法，要求
次运算，其中
我们可以得到，如果r j=s j t j且s j,t j>1,那么s j+ t j<r j(当s j=t j=2时，s j+t j=r j)。

总的说来，如果我们尽可能多地分解N，我们就可以为（9）式的T找到一个最小值，而如果使用2作为因子，则没有任何增加/减少。

如果我们能够选择一个N 成为高度合成数，我们可以获得真正的便利。

假设如果所有的r j等于r，那么从（10）我们可以得到
与所有的操作数为
如果N=r m s n t p…，那么我们可以发现
所以
将会是一个带有权重的平均量对于下列各式
而它们的值分别是
使用r j=3是形式上最有效的，但是得到的便利也仅仅之比使用2或者4高了6%，而2或4还有另外层面上的优势。

如若必须，取r j的值高达10的时候，会将运算的数量增加接近于50%。

因此，我们可以发现，对于N来说，高度合成数的值仅仅只占据很小的比例在任意给定的大数中。

无论何时应用，使用r=2或4带入N=r m会为二进制的计算带来显著地优势，不仅在寻址上还是在乘法的经济程度上。

带入r=2的算法可以表示成
其中j v和k v都等于0或1，且就是在二进制表示j和k时，各个位地址上的值。

现在所有的数组都可以表示成关于各位地址上的值的函数。

于是我们可以将（1）写成
其中各累加和有关于各个k v=0，1的计算。

又因为
式（15）中最大的内和（也就是k m-1上的计算），仅仅取决于j0，k m-2，...，k0，于是该式可以被写成
继续运算至下一个最大的内和（k m-2上的计算），然后不断循环运算下去，不断使用
我们可以得到一个后续的数组，
其中l=1，2，...，m。

将累加和写成下述这种形式
根据进制表示的规定，这个和将会保存在二进制编码为下述的地址上。

我们可以发现，在式（20）中，仅仅两个在2m-l的位地址上用二进制编码的存储地址被调用来处理运算。

又因为对于式（20）中表示的运算，计算机可以同时执行各个位地址j0，...，j l-2和k0，...，k m-l-1的值的计算，所以并行计算是可行
的。

在一些应用中*，这种表示可以方便的用A l-2来表达A l，假定使用等价于r=4的算法。

最后一个计算得到的矩阵就是我们所要的傅里叶变化的和，
在这样的规则顺序下，X的二进制编码就必须倒序来产生正确的数组A m的编码。

在一些应用中，傅里叶变化的累加和需要被应用两次，此时上述的过程就不需要再进行倒序的过程。

例如，考虑一个差分方程的解，
我们的方法可以首先用来计算傅里叶变化中F(j)的幅值，其中，F(j)将带入下式：
于是，解的幅值就为
其中，B(k)和A(k)矩阵都应该成倒序，但是根据式（20）的修正，A(k)就可以用来表示正确的解。

在IBM7409上运行的一个计算机程序用来计算应用上述方法的3阶傅里叶变化和。

用来运算2a*2b*2c的数组大小的数据所花的计算时间如下表所示：
*在ler和S.Winogard (IBM Watson Research Center)设计的基于此算法多步循环的程序中，发现带入r=4是最实用的。