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 code 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));

真的，覺得還是紙上方便。