快速傅裡葉變換FFT的C語言演算法徹底研究

快速傅裡葉變換FFT的C語言演算法徹底研究

LED音樂頻譜顯示的核心演算法就是快速傅裡葉變換，FFT的理解和編程還是比較難的，特地撰寫此文分享一下研究成果。

 一、徹底理解傅裡葉變換

快速傅裡葉變換（Fast Fourier Transform)是離散傅裡葉變換的一種快速演算法，簡稱FFT，通過FFT可以將一個訊號從時域變換到頻域。

類比訊號經過A/D轉換變為數字訊號的過程稱為採樣。為保證採樣後訊號的頻譜形狀不失真，採樣頻率必須大於訊號中最高頻率成分的2倍，這稱之為採樣定理。

假設採樣頻率為fs，採樣點數為N，那麼FFT結果就是一個N點的複數，每一個點就對應著一個頻率點，某一點n(n從1開始)表示的頻率為：fn=(n-1)*fs/N。

舉例說明：用1kHz的採樣頻率採樣128點，則FFT結果的128個資料即對應的頻率點分別是0，1k/128，2k/128，3k/128，…，127k/128 Hz。

這個頻率點的幅值為：該點複數的模值除以N/2（n=1時是直流分量，其幅值是該點的模值除以N）。

二、傅裡葉變換的C語言編程

1、對於快速傅裡葉變換FFT，第一個要解決的問題就是碼位倒序。

    假設一個N點的輸入序列，那麼它的序號位元位元就是t=log2N.

    碼位倒序要解決兩個問題：①將t位位元倒序；②將倒序後的兩個儲存單元進行交換。

①將t=3位位元倒序

②將倒序後的兩個儲存單元進行交換

 如果輸入序列的自然順序號i用位元表示，例如若最大序號為15，即用4位就可表示n3n2n1n0，則其倒序後j對應的位元就是n0n1n2n3，那麼怎樣才能實現倒序呢。利用C語言的移位功能。

程式如下，我不多說，看不懂者智商一定在180以下。

複數類型定義及其運算

#define N 64      //64點

#define log2N 6   //log2N=6

/*複數類型*/

typedef struct

{

 float real;

 float img;

}complex;

complex xdata x[N]; //輸入序列

/*複數加法*/

complex add(complex a,complex b)

{

 complex c;

 c.real=a.real+b.real;

 c.img=a.img+b.img;

 return c;

}

/*複數減法*/

complex sub(complex a,complex b)

{

 complex c;

 c.real=a.real-b.real;

 c.img=a.img-b.img;

 return c;

}

/*複數乘法*/

complex mul(complex a,complex b)

{

 complex c;

 c.real=a.real*b.real - a.img*b.img;

 c.img=a.real*b.img + a.img*b.real;

 return c;

}

/***碼位倒序函數***/

void Reverse(void)

{

 unsigned int i,j,k;

 unsigned int t;

 complex temp;//臨時交換變數

 for(i=0;i<N;i++)//從第0個序號到第N-1個序號

 {

  k=i;//當前第i個序號

  j=0;//儲存倒序後的序號，先初始化為0

  for(t=0;t<log2N;t++)//共移位t次，其中log2N是事先宏定義算好的

  {

   j<<=1;

   j|=(k&1);//j左移一位然後加上k的最低位

   k>>=1;//k右移一位，次低位變為最低位

  }

  if(j>i)//如果倒序後大於原序數，就將兩個儲存單元進行交換（判斷j>i是為了防止重複交換）

  {

   temp=x;

   x=x[j];

   x[j]=temp;

  }

 }

}

2、第二個要解決的問題就是蝶形運算

①第1級（第1列）每個蝶形的兩節點“距離”為1，第2級每個蝶形的兩節點“距離”為2，第3級每個蝶形的兩節點“距離”為4，第4級每個蝶形的兩節點“距離”為8。由此推得，

 第m級蝶形運算，每個蝶形的兩節點“距離”L=2m-1。

②對於16點的FFT，第1級有16組蝶形，每組有1個蝶形；第2級有4組蝶形，每組有2個蝶形；第3級有2組蝶形，每組有4個蝶形；第4級有1組蝶形，每組有8個蝶形。由此可推出，

對於N點的FFT，第m級有N/2L組蝶形，每組有L=2m-1個蝶形。

③旋轉因子的確定

以16點FFT為例，第m級第k個旋轉因子為，其中k=0～2m-1-1，即第m級共有2m-1個旋轉因子，根據旋轉因子的可約性，，所以第m級第k個旋轉因子為，其中k=0～2m-1-1。

為提高FFT的運算速度，我們可以事先建立一個旋轉因子數組，然後通過查表法來實現。

complex code WN[N]=//旋轉因子數組
{ //為節省CPU計算時間，旋轉因子採用查表處理
  //★根據實際FFT的點數N，該表資料需自行修改
  //以下結果通過Excel自動產生
  //  WN[k].real=cos(2*PI/N*k);
  //  WN[k].img=-sin(2*PI/N*k);

  {1.00000,0.00000},{0.99518,-0.09802},{0.98079,-0.19509},{0.95694,-0.29028},
  {0.92388,-0.38268},{0.88192,-0.47140},{0.83147,-0.55557},{0.77301,-0.63439},
  {0.70711,-0.70711},{0.63439,-0.77301},{0.55557,-0.83147},{0.47140,-0.88192},
  {0.38268,-0.92388},{0.29028,-0.95694},{0.19509,-0.98079},{0.09802,-0.99518},
  {0.00000,-1.00000},{-0.09802,-0.99518},{-0.19509,-0.98079},{-0.29028,-0.95694},
  {-0.38268,-0.92388},{-0.47140,-0.88192},{-0.55557,-0.83147},{-0.63439,-0.77301},
  {-0.70711,-0.70711},{-0.77301,-0.63439},{-0.83147,-0.55557},{-0.88192,-0.47140},
  {-0.92388,-0.38268},{-0.95694,-0.29028},{-0.98079,-0.19509},{-0.99518,-0.09802},
  {-1.00000,0.00000},{-0.99518,0.09802},{-0.98079,0.19509},{-0.95694,0.29028},
  {-0.92388,0.38268},{-0.88192,0.47140},{-0.83147,0.55557},{-0.77301,0.63439},
  {-0.70711,0.70711},{-0.63439,0.77301},{-0.55557,0.83147},{-0.47140,0.88192},
  {-0.38268,0.92388},{-0