快速傅立葉變換(FFT)C語言函式
阿新 • • 發佈:2019-02-18
/********************************************************************* 快速福利葉變換C函式 函式簡介:此函式是通用的快速傅立葉變換C語言函式,移植性強,以下部分不依 賴硬體。此函式採用聯合體的形式表示一個複數,輸入為自然順序的復 數(輸入實數是可令複數虛部為0),輸出為經過FFT變換的自然順序的 複數 使用說明:使用此函式只需更改巨集定義FFT_N的值即可實現點數的改變,FFT_N的 應該為2的N次方,不滿足此條件時應在後面補0 函式呼叫:FFT(s); 作 者:吉帥虎 時 間:2010-2-20 版 本:Ver1.0 參考文獻: **********************************************************************/ //#include <reg52.h> //AT89C52 //#include <iom128.h> //atmeg128 //#include <intrinsics.h> #include <math.h> #define PI 3.1415926535897932384626433832795028841971 //定義圓周率值 #define FFT_N 128 //定義福利葉變換的點數 struct compx {float real,imag;}; //定義一個複數結構 struct compx s[FFT_N]; //FFT輸入和輸出:從S[1]開始存放,根據大小自己定義 /******************************************************************* 函式原型:struct compx EE(struct compx b1,struct compx b2) 函式功能:對兩個複數進行乘法運算 輸入引數:兩個以聯合體定義的複數a,b 輸出引數:a和b的乘積,以聯合體的形式輸出 *******************************************************************/ struct compx EE(struct compx a,struct compx b) { struct compx c; c.real=a.real*b.real-a.imag*b.imag; c.imag=a.real*b.imag+a.imag*b.real; return(c); } /***************************************************************** 函式原型:void FFT(struct compx *xin,int N) 函式功能:對輸入的複數組進行快速傅立葉變換(FFT) 輸入引數:*xin複數結構體組的首地址指標,struct型 *****************************************************************/ void FFT(struct compx *xin) { int f , m, nv2, nm1, i, k, l, j = 0; struct compx u,w,t; nv2 = FFT_N / 2; //變址運算,即把自然順序變成倒位序,採用雷德演算法 nm1 = FFT_N - 1; for(i = 0; i < nm1; i++) { if(i < j) //如果i<j,即進行變址 { t = xin[j]; xin[j] = xin[i]; xin[i] = t; } k = nv2; //求j的下一個倒位序 while( k <= j) //如果k<=j,表示j的最高位為1 { j = j - k; //把最高位變成0 k = k / 2; //k/2,比較次高位,依次類推,逐個比較,直到某個位為0 } j = j + k; //把0改為1 } { int le , lei, ip; //FFT運算核,使用蝶形運算完成FFT運算 f = FFT_N; for(l = 1; (f=f/2)!=1; l++) //計算l的值,即計算蝶形級數 ; for( m = 1; m <= l; m++) // 控制蝶形結級數 { //m表示第m級蝶形,l為蝶形級總數l=log(2)N le = 2 << (m - 1); //le蝶形結距離,即第m級蝶形的蝶形結相距le點 lei = le / 2; //同一蝶形結中參加運算的兩點的距離 u.real = 1.0; //u為蝶形結運算係數,初始值為1 u.imag = 0.0; w.real = cos(PI / lei); //w為係數商,即當前係數與前一個係數的商 w.imag = -sin(PI / lei); for(j = 0;j <= lei - 1; j++) //控制計算不同種蝶形結,即計算係數不同的蝶形結 { for(i = j; i <= FFT_N - 1; i = i + le) //控制同一蝶形結運算,即計算係數相同蝶形結 { ip = i + lei; //i,ip分別表示參加蝶形運算的兩個節點 t = EE(xin[ip], u); //蝶形運算,詳見公式 xin[ip].real = xin[i].real - t.real; xin[ip].imag = xin[i].imag - t.imag; xin[i].real = xin[i].real + t.real; xin[i].imag = xin[i].imag + t.imag; } u = EE(u, w); //改變係數,進行下一個蝶形運算 } } } } /************************************************************ 函式原型:void main() 函式功能:測試FFT變換,演示函式使用方法 輸入引數:無 輸出引數:無 ************************************************************/ void main() { int i; for(i = 0; i < FFT_N; i++) //給結構體賦值 { s[i].real = 1 + 2*sin(2*3.141592653589793*i / FFT_N); //實部為正弦波FFT_N點取樣,賦值為1 s[i].imag = 0; //虛部為0 } FFT(s); //進行快速福利葉變換 for(i = 0; i < FFT_N; i++) //求變換後結果的模值,存入複數的實部部分 s[i].real = sqrt(s[i].real*s[i].real + s[i].imag * s[i].imag); while(1); }