快速傅利葉變換的C語言實現

快速傅利葉變換(fft)利用了旋轉因子e−j

2πne^}

e−jn2π

的週期性、共軛對稱性以及可約性極大簡化了dft的計算量，具體可以查閱清華大學出版社數字訊號處理程佩青第四版第四章。

理論推導一大堆，最重要的是對於蝶形圖的理解。書上有8點的蝶形圖，這裡給出16點的蝶形圖。

這張圖其實有點小問題，所有的橫線都畫漏了。

下面討論c語言實現，此處實現dit的基2 fft演算法。分為兩步，第一步是對於輸入序列的重新排序，以確保輸出的頻域值是順序的；第二步是進行複數的乘法和加法，每乙個蝶形需要兩次複數加法（其實是一次加一次減）以及一次複數乘法。輸入的點個數必須是2的冪次，如果不是，需要補0讓輸入點數滿足2的冪次。

第一步採用rader演算法，就是將乙個數轉化成倒位序。以8點為例：

原序列fft序列

原序列二進位制（i）

fft序列二進位制（j）00

00000014

00110022

01001036

01111041

10000155

10110163

11001177

111111

從表中可以看出，要將輸入序列換成需要的fft序列，則只需要將原序列的二進位制表示倒過來就可以了。(當初我發現這一點時驚呆了！)

設輸入點數為n，則蝶形的級數為t

=log⁡2

nt=\log_2n

t=log2n

，表示序列的2進製位數也是t位。

rader演算法實現原理是：i，j都從0開始，若已知某個倒位序j，要求下乙個倒位序數，則應先判斷j的最高位是否為0，將j與k=n/2相比較，如果k>j，則j的最高位為0，只要把該位變為1（j與k=n/2相加即可），就得到下乙個倒位序數；如果k<=j，則j的最高位為1，可將最高位變為0（j與k=n/2相減即可）。然後還需判斷次高位，這可與k=n/4相比較，若次高位為0，則需將它變為1（加n/4即可）其他位不變，既得到下乙個倒位序數；若次高位是1，則需將它也變為0。然後再判斷下一位，以此類推。

c語言實現：

nv2=fft_n/2;                  //變址運算，即把自然順序變成倒位序，採用雷德演算法 
nm1=fft_n-1; 
for(i=0;i再給出一種實現方式，利用按位與來實現。
void change() 
if(j>i) //將x(n)的碼位互換，其實只需要迴圈size_x/2次就行了（待考證），因為i>size_x/2時，倒序必定比i小。 
} output(); 
}

接下來是蝶形運算的處理

這裡以8點fft為例子，放上一張：

再放上fft核心c**：

void fft()  
{ int i=0,j=0,k=0,l=0; 
complex up,down,product; 
change(); //呼叫變址函式 
for(i=0;i< log(size_x)/log(2) ;i++) /*一級蝶形運算 stage */ 
{ 
l=1⁡2n
\log_2n
log2n
級，我們這裡的8個點就是3級。我們計算fft結果的過程其實就是將x[0]->x[7]完成**蝶形運算的過程，通過三重迴圈實現，每一重迴圈完成一級的運算。例如在圖上已經標註出來的，第一級運算之後x[0]變成了x[0]+x[1]*w[0]，把新x[0]再帶入下一級的運算中。**當然了，在第一級運算開始之前，你需要change函式完成倒序；還需要表示出旋轉因子wn
i_n^i
ni。
我還手寫了乙份上述c**迴圈執行過程，幫助理解：
小結一下，dft是dtft的離散形式，fft是dft的計算方法，fft只是乙個算係數的過程。

快速傅利葉變換的C語言實現

離散傅利葉變換快速傅利葉變換C 實現

c語言做快速傅利葉變換和快速逆傅利葉變換

快速傅利葉變換C 遞迴演算法實現

快速傅利葉變換的C語言實現

離散傅利葉變換 快速傅利葉變換C 實現

c語言做快速傅利葉變換和快速逆傅利葉變換

快速傅利葉變換C 遞迴演算法實現

相關推薦

離散傅利葉變換快速傅利葉變換C 實現