c - FFT 和 FFT 的逆

Question

我是从事电信项目的计算机程序员。
在我们的项目中，我必须将一系列复数更改为它们的傅立叶变换。因此我需要一个高效的 FFT 代码来满足 C89 标准。
我正在使用以下代码，它运行良好:

    short FFT(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(true);
}

但是这段代码，只支持大小为2^m的数组。如CLRS书本代码。
我们应该转换的数组不在这个范围内，添加零会很昂贵，所以我正在寻找另一种解决方案来帮助我输入任何大小。
类似于 IT++ 和 matlab 所做的事情。但正如我们希望在纯 C 中使用它们是不可能的。此外，IT++ 代码在我检查时被阻止

Answer 1

如果您正在使用任何大众市场计算平台(带有 Windows 或 OS X、iOS 等的英特尔)，那么供应商或制造商会提供高性能 FFT 实现。

否则，你应该评估FFTW .

为除 2 的幂以外的大小编写高性能 FFT 是一项复杂的任务。

如果您要使用自己的实现，那么，关于二次方大小:

您展示的实现在 FFT 期间计算 sqrt。大多数高性能 FFT 实现提前计算常量并将它们存储在表中。

缩放包含除法运算，在 x[i]/= n 和 y[i]/= n 中。编译器很可能将这些实现为除法指令。除法在普通处理器上通常是一条缓慢的指令。计算一次 scale = 1./n 并乘以 scale 而不是除以 n 会更好。

更好的办法是完全省略比例。转换通常在没有比例的情况下很有用，或者可以从单个转换中省略比例并仅作为聚合比例应用一次。 (例如，与其执行两个缩放操作，一个在正向 FFT 中，一个在反向 FFT 中，将缩放操作从 FFT 例程中移除，并在正向 FFT 和反向 FFT 之后只执行一次。)

如果可以接受按位反转顺序排列频域数据，则您可以省略位反转排列。

如果保留位反转排列，则可以对其进行优化。执行此操作的技术取决于平台。某些平台有反转整数位的指令(例如，ARM 有 rbit)。如果您的平台没有，您可能希望将位反转索引保存在表中或研究比当前代码更快地计算它们的方法。

如果同时保留位反转排列和缩放，则应考虑同时进行。排列使用大量内存运动，缩放使用处理器的算术单元。大多数现代处理器可以同时执行这两项操作，因此您可以从重叠操作中获益。

您当前的代码使用 radix-2 butterfly。 Radix-4 通常更好，因为它从以下事实中获得了一些优势:只需更改使用的数据片段并将一些加法更改为减法，反之亦然即可完成与 i 的乘法。

如果您的数组长度接近处理器上一级内存缓存的大小，FFT 实现的某些部分将扰动缓存并显着降低速度，除非您设计适当的代码来处理这个问题(通常通过将数组的部分复制到一个临时缓冲区)。

如果您的目标处理器具有 SIMD 功能，您绝对应该在 FFT 中使用这些功能；它们极大地加快了 FFT 性能。

以上内容应该告诉您，编写高效的 FFT 是一项复杂的任务。除非您想为此花费大量精力，否则最好使用 FFTW 或其他现有实现。