FFT与FFT的逆

Question

我是一个从事电信项目的计算机程序员。

在我们的项目中，我必须将一系列复数转换为它们的傅里叶transform.so，我需要一个有效的C89标准的FFT代码。

我正在使用以下代码，它运行得很好：

    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.like CLRS图书代码大小的数组。

我们的数组，应该转换，不是在这个范围内，加上零将是昂贵的，所以我正在寻找另一个解决方案，帮助我有任何大小的输入。

就像IT++和matlab所做的那样。但是，正如我们希望在纯C中使用的是impossible.Moreover，在我检查时，IT++代码被阻塞了

Answer 1

如果您正在从事任何大众市场计算平台(英特尔与视窗或OS，iOS等)，那么有一个高性能的快速傅立叶变换实现提供的供应商或制造商。

否则，您应该评估FFTW。

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

如果您要使用自己的实现，那么，考虑到两种大小的能力：

您显示的实现在快速傅立叶变换期间计算sqrt。大多数高性能FFT实现提前计算常量并将它们存储在表中。

缩放包含x[i] /= n和y[i] /= n中的除法操作。编译器可能会将这些实现为除法指令。除法通常是通用处理器上的慢速指令。最好是一次计算scale = 1. / n，然后用scale乘以，而不是用n除以。

更好的办法是完全忽略这个规模。转换通常在没有标度的情况下是有用的，或者可以从单个转换中省略该比例，并将其作为聚合规模应用一次。(例如，不执行两次缩放操作，一次在正向FFT中，一次在逆FFT中，将缩放操作排除在FFT例程之外，只在正向FFT和逆FFT之后执行一次。)

如果有位反转顺序的频域数据是可以接受的，您可能可以省略位反转排列。

如果你保持位反转排列，它可以被优化。这样做的技术依赖于平台。有些平台有一个指令来反转整数中的位(例如，ARM有rbit)。如果您的平台没有，您可能希望将位反转索引保存在表中，或者研究比当前代码更快地计算它们的方法。

如果你同时保持位反转排列和缩放，你应该考虑同时做它们。置换使用大量的内存运动，而缩放使用处理器的算术单元。大多数现代处理器可以同时完成这两种操作，因此您可以从重叠操作中获得一些好处。

当前的代码使用基-2蝶形。基数-4通常更好，因为它从这样一个事实中获得了一些优势，即通过改变使用哪一段数据并更改减法的一些加法，就可以实现I的乘法，反之亦然。

如果您的数组长度接近处理器上一级内存缓存的大小，则FFT实现的部分将破坏缓存并显著减慢，除非您设计了适当的代码来处理这个问题(通常是将数组的部分临时复制到缓冲区中)。

如果您的目标处理器具有SIMD特性，那么您应该绝对使用FFT中的那些特性；它们极大地提高了FFT的性能。

以上所述应向您表明，编写高效的FFT是一项复杂的任务。除非您想在这方面花费大量精力，否则使用FFTW或其他现有实现可能会更好。

Answer 2

在您的实现中，我关注以下代码：

     z =  u1 * c1 - u2 * c2;
     u2 = u1 * c2 + u2 * c1;
     u1 = z;

(u1，u2)当l1较大时，会产生大量累积舍入误差。你可能会得到一个不准确的转换。

我会赞同FFTW的建议，但我相信这是非常特定于平台的。(大多数FFT库是。)编辑：不，它实际上是直C89。这正是你说你要找的。

维基百科的FFT页面列出了一组适用于奇怪大小的输入数组的算法。我不知道它们是如何工作的，但我相信一般的想法是使用Rader算法或Bluestein算法作为主要大小的输入，使用库利-图基将组合大小的转换减少为一组素数大小的转换。

Answer 3

要获得FFTW的替代方案，请查看我的混合基FFT，它还处理非2^N FFT长度。C源可以从http://www.corix.dk/Mix-FFT/mix-fft.html获得.