FFT（快速傅里叶变换）程序分析

      程序代码来源于：http://local.wasp.uwa.edu.au/~pbourke/miscellaneous/dft
        在该网站中不仅提供了FFT程序的源代码，也进行了些理论的分析。
        感谢程序的作者和信号与系统的研究者。Thanks for all the people contributed to System and SIngal。
        首先帖上程序的源代码：感谢程序源代码作者Peter Cusack！

/*
  Modification of Paul Bourkes FFT co

de by Peter Cusack
  to utilise the Microsoft complex type.

  This computes an in-place complex-to-complex FFT
  x and y are the real and imaginary arrays of 2^m points.
  dir =  1 gives forward transform
  dir = -1 gives reverse transform
*/
void FFT(int dir, long m, complex <double> x[])
{
  long i, i1, i2,j, k, l, l1, l2, n;
  complex <double> tx, t1, u, c;

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

  /* Do the bit reversal */
  i2 = n >> 1;
  j = 0;

/*首先首先位倒序输入*/
  for (i = 0; i < n-1 ; i++)
    {
      if (i < j)
    swap(x[i], x[j]);

      k = i2;

      while (k <= j)
    {
      j -= k;
      k >>= 1;
    }

      j += k;
   }

  /* Compute the FFT */
  c.real(-1.0);
  c.imag(0.0);
  l2 = 1;
  for (l = 0; l < m; l++)
    {
      l1 = l2;
      l2 <<= 1;
      u.real(1.0);
      u.imag(0.0);

/* 使用相同旋转因子的行计算下一列的值 */
      for (j = 0; j < l1; j++)
    {
      for (i = j; i < n; i += l2) {
            i1 = i + l1;
            t1 = u * x[i1];
            x[i1] = x[i] - t1;
            x[i] += t1;
     }

/* 旋转因子乘以旋转因子差值 */
      u = u * c;
    }

     /* 做FFT旋转因子的计算 （后面将进行详细分析）*/
      c.imag(sqrt((1.0 - c.real()) / 2.0));
      if (dir == 1)
    c.imag(-c.imag());
      c.real(sqrt((1.0 + c.real()) / 2.0));
    }

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

基2 FFT算法主要使用的是蝶形运算进行FFT计算。

    图1 : FFT蝶形运算示意图（其中C即是蝶形运行的旋转因子）

  采用蝶形运算，某一列的任何两个节点K,J的变量进行蝶形运算后，得到的仍是下一列蝶形运算的节点K,J。由于蝶形运算采用的是迭代的方式实现，那么旋转因子的计算即成为了蝶形运算的重要步骤。

   

图2：旋转因子的规律（第L级的旋转因子，N=2^M）

  可以看到蝶形运算的旋转因子在同一L时，相邻即J和J+1列旋转因子系数相关2^(M-L)，并且在上一列(L-1)列的相同J(J相同)时旋转因子系数之差存在2倍的关系。我们所以只需要计算L列时旋转因子系数的差值即2^(M-L)，就可以计算得到所有的旋转因子的值。（因为当J=0时，旋转因子为值 1，然后再循环乘以就可以得到所有的旋转因子值）。

  在程序中并没有首先计算所有的旋转因子的值，而是首先计算使用相同旋转因子的行下一列的值。
 

图三：计算过程（比如第二列首先会计算使用的行）

     旋转因子的程序计算方法：由于当前列和上一列旋转因子系数的差值存在2倍关系，即L列旋转因子系数的差值是第L+1列旋转因子系数差值的2倍。所以利用程序迭代和上一次的旋转因子即可以计算当前旋转因子值的（倍差）【乘法关系】。

   旋转因子计算推导：设L-1列旋转因子系数差值为X1, L列旋转因子系数差值为X2
   则旋转因子倍差值，第L-1列倍值为，第L列倍值为，由于系数差值存在倍数关系，
即：X1 = 2 * X2;  设 M = , T = -j * 2 * pi / N;
   = = cos(X1 * T) - j sin(X1 * T) = c1.real + jc1.img;      (1-1)式
   = = cos(X2 * T) - j sin(X2 * T) = c2.real + jc2.img;      (1-2)式
  
   现在的目的即是用c1.real 和 c1.img 表示 c2.real 和 c2.img;

  由1-1可以得到
  cos(X2 * T) * cos(X2 * T) - j * 2 * sin(X2 * T) * cos(X2 * T) = c1.real + j c1.img;  (1-3)
  2cos(X2 * T) * cos(X2 * T) - 1 - sin(X2 * T) * cos(X2 * T) = c1.real + j c1.img;   (1-4)

所以：
2 * c2.real * c2.real = c1.real + 1;         (1-5)
c2.real = sqrt(c1.real + 1) / 2;               (1-6)

由 1-3得 2 * c2.img * c2.real = c1.img;    (1-7)
4 * c2.img * c2.img * c2.real * c2.real = c1.img * c1.img = 1 - c1.real * c1.real;   (1-8)
2 * (c1.real + 1) * c2.img * c2.img = (1 - c1.real) * (1 + c1.real);    (1-9)
2 * c2.img * c2.img = 1 - c1.real;              (1-10)
所以： c2.img = sqrt((1 - c1.real) / 2);        (1-11)

所以: (1-6)和(1-11)即是程序代码表达的公式：
c.imag(sqrt((1.0 - c.real()) / 2.0));
c.real(sqrt((1.0 + c.real()) / 2.0));


真的，觉得还是纸上方便。