提高FFT实现速度

Question

我是编程的初学者，目前正在尝试一个需要快速傅立叶变换实现的项目。

到目前为止，我已经成功地实现了以下内容：

有没有人有任何替代方案和建议来提高程序的速度，而不会失去准确性。

short FFTMethod::FFTcalc(short int dir,long m,double *x,double *y)
{
long n,i,i1,j,k,i2,l,l1,l2;
double c1,c2,tx,ty,t1,t2,u1,u2,z;

/* Calculate the number of points */
n = 1;
for (i=0;i<m;i++) 
    n *= 2;

/* Do the bit reversal */
i2 = n >> 1;
j = 0;
for (i=0;i<n-1;i++) {
  if (i < j) {
     tx = x[i];
     ty = y[i];
     x[i] = x[j];
     y[i] = y[j];
     x[j] = tx;
     y[j] = ty;
  }
  k = i2;
  while (k <= j) {
     j -= k;
     k >>= 1;
  }
  j += k;
}

/* Compute the FFT */
c1 = -1.0; 
c2 = 0.0;
l2 = 1;
for (l=0;l<m;l++) {
   l1 = l2;
   l2 <<= 1;
   u1 = 1.0; 
   u2 = 0.0;
   for (j=0;j<l1;j++) {
     for (i=j;i<n;i+=l2) {
        i1 = i + l1;
        t1 = u1 * x[i1] - u2 * y[i1];
        t2 = u1 * y[i1] + u2 * x[i1];
        x[i1] = x[i] - t1; 
        y[i1] = y[i] - t2;
        x[i] += t1;
        y[i] += t2;
     }
     z =  u1 * c1 - u2 * c2;
     u2 = u1 * c2 + u2 * c1;
     u1 = z;
   }
   c2 = sqrt((1.0 - c1) / 2.0);
   if (dir == 1) 
     c2 = -c2;
     c1 = sqrt((1.0 + c1) / 2.0);
  }

/* Scaling for forward transform */
if (dir == 1) {
   for (i=0;i<n;i++) {
      x[i] /= n;
      y[i] /= n;
   }
 } 


   return(1);
}

Answer 1

我最近在Construction of a high performance FFTs上找到了Eric Postpischil写的这篇很棒的PDF。我自己开发过几个FFT，我知道与商业图书馆竞争有多难。相信我，如果你的FFT只比Intel或FFTW慢4倍，而不是40倍，那你就做得很好！然而，你可以竞争，下面是如何竞争的。

总结这篇文章，作者指出Radix2 FFT简单但效率低下，最有效的结构是radix4 FFT。一种更有效的方法是Radix8，但是这通常不适合CPU上的寄存器，所以首选Radix4。

Radix4可以分阶段构建，因此要计算1024点的FFT，您可以执行10个阶段的FFT (作为2^10 - 1024)，或者执行5个阶段的FFT (4^5 = 1024)。如果你选择的话，你甚至可以在8*4*4*4*2的阶段计算1024点的FFT。更少的级数意味着对内存的读取和写入更少( FFT性能的瓶颈是内存带宽)，因此必须动态选择基数4、8或更高的基数。Radix4阶段是特别有效的，因为所有的权重都是1+0i，0+1i，-1+0i，0-1i和Radix4蝶形代码可以写成完全适合高速缓存。

其次，FFT中的每个阶段都不相同。第一阶段的权重都等于1+0i。计算这个权重甚至乘以它都没有意义，因为它是一个复数乘以1，所以第一阶段可以在没有权重的情况下执行。最后一级也可以被不同地处理，并且可以用于在时间上执行抽取(比特反转)。Eric Postpischil的文档涵盖了所有这些内容。

可以预先计算权重并将其存储在表中。在x86硬件上，Sin/cos计算每次大约需要100-150个周期，因此预计算这些计算可以节省总计算时间的10-20%，因为在这种情况下，内存访问比CPU计算更快。使用快速算法一次性计算sincos特别有好处(请注意cos等于sqrt(1.0 -sine* in )，或者使用查表，cos只是正弦的相移)。

最后，一旦你有了超级流线型的FFT实现，你就可以利用SIMD矢量化在蝶形例程中的每个循环中计算4倍浮点或2倍双浮点操作，以进一步提高100-300%的速度。综合以上所有这些，你会有一个相当圆滑和快速的快速傅立叶变换！

更进一步，您可以通过提供针对特定处理器体系结构的FFT阶段的不同实现来动态执行优化。缓存大小、寄存器计数、SSE/SSE2/3/4指令集等因机器不同而不同，因此选择一种适合所有人的方法通常会被目标例程击败。例如，在FFTW中，许多较小的FFT都是针对特定架构的高度优化的展开(无循环)实现。通过组合这些较小的构造(比如RadixN例程)，您可以为手头的任务选择最快、最好的例程。

Answer 2

虽然我现在不能给你一个性能提示，但我想为你的优化提供一些建议，因为评论太长了：

1. 如果您还没有这样做，那么现在就为您的代码编写一些正确性测试。像“对这个数组进行快速傅立叶变换，看看结果是否与我提供的结果匹配”这样的简单测试就足够了，但在优化代码之前，您需要一个可靠的自动化单元测试来确认优化后的代码是正确的。然后，
2. 会评测您的代码，看看真正的瓶颈在哪里。虽然我怀疑最里面的循环for (i=j;i<n;i+=l2) {，但眼见为实好于相信。--

Answer 3

有几个我可以推荐尝试的方法：

1. 不交换输入元素，而是计算位反转索引。这将为您节省大量的内存读取和writes.
2. Precalculate系数，如果您正在进行许多相同大小的FFT。这将节省一些computations.
3. Use基数-4FFT，而不是基数-2。这将减少内部循环中的迭代次数。

当然，通过分析代码可以找到最终的答案。