C語言實現(xiàn)FFT快速傅里葉變換

1、#include <iom128.h>#include <intrinsics.h>/*快速福利葉變換C函數(shù)函數(shù)簡介：此函數(shù)是通用的快速傅里葉變換C語言函數(shù)，移植性強，以下部分不依賴硬件。此函數(shù)采用聯(lián)合體的形式表示一個復(fù)數(shù)，輸入為自然順序的復(fù) 數(shù)(輸入實數(shù)是可令復(fù)數(shù)虛部為0)，輸出為經(jīng)過FFT變換的自然順序的復(fù)數(shù)使用說明：使用此函數(shù)只需更改宏定義FFT_N的值即可實現(xiàn)點數(shù)的改變，F(xiàn)FT_N的應(yīng)該為2的N次方，不滿足此條件時應(yīng)在后面補0FFT(s); 2010-2-20函數(shù)調(diào)用：時 間:版 本:參考文獻：Ver1.0*/#include<math.h>#def

2、ine PI 3.1415926535897932384626433832795028841971#define FFT_N 128變換的點數(shù)定義圓周率值定義福利葉struct compx float real,imag;struct compx sFFT_N;出：從S1開始存放，根據(jù)大小自己定義定義一個復(fù)數(shù)結(jié)構(gòu)/FFT輸入和輸/*函數(shù)原型函數(shù)功能輸入?yún)?shù)輸出參數(shù)struct compx EE(struct compx b1,struct compx b2)對兩個復(fù)數(shù)進行乘法運算兩個以聯(lián)合體定義的復(fù)數(shù)a,ba和b的乘積，以聯(lián)合體的形式輸出*/struct compx EE(struct com

3、px 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);/*函數(shù)原型：void FFT(struct compx *xin,int N)(int le,lei,ip;f=FFT_N;for(l=1;(f=f/2)!=1;l+);for(m=1;m<=l;m+)(l=log (2) Nle=2<<(m-1);相距l(xiāng)e點lei=le/2;u.real=1.0;u.imag=0.0;w.real=cos(PI

4、/lei);w.imag=-sin(PI/lei);for(j=0;j<=lei-1;j+)函數(shù)功能：對輸入的復(fù)數(shù)組進行快速傅里葉變換(FFT)輸入?yún)?shù)：*xin復(fù)數(shù)結(jié)構(gòu)體組的首地址指針，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=xinj;xinj=xini;xini=t;k=

5、nv2;/求j的下一個倒位序while(k<=j)如果k<=j,表示j的最高位為1(j=j-k;/把最高位變成0k=k/2;/k/2，比較次高位，依次類推，逐個比較，直到某個位為0j=j+k;把0改為1/FFT運算核，使用蝶形運算完成FFT運算計算l的值，即計算蝶形級數(shù)/控制蝶形結(jié)級數(shù)/m表示第m級蝶形，l為蝶形級總數(shù)/le蝶形結(jié)距離，即第 m級蝶形的蝶形結(jié)/同一蝶形結(jié)中參加運算的兩點的距離/u為蝶形結(jié)運算系數(shù)，初始值為1/w為系數(shù)商，即當前系數(shù)與前一個系數(shù)的商/控制計算不同種蝶形結(jié)，即計算系數(shù)不同的蝶形結(jié)(for(i=j;i<=FFT_N-1;i=i+le)結(jié)(ip=i+l

6、ei;t=EE(xinip,u);xinip.real=xini.real-t.real;xinip.imag=xini.imag-t.imag;xini.real=xini.real+t.real;xini.imag=xini.imag+t.imag;u=EE(u,w);控制同一蝶形結(jié)運算，即計算系數(shù)相同蝶形/i, ip分別表示參加蝶形運算的兩個節(jié)點蝶形運算，詳見公式改變系數(shù)，進行下一個蝶形運算/*函數(shù)原型函數(shù)功能輸入?yún)?shù)輸出參數(shù)void main()測試FFT變換，演示函數(shù)使用方法無無*/void main()(int i;for(i=0;i<FFT_N;i+)( si.real=s

7、in(2*3.141592653589793*i/FFT_N); / si.imag=0;/給結(jié)構(gòu)體賦值實部為正弦波FFT_N點采樣，賦值為1虛部為0FFT(s);/進行快速福利葉變換for(i=0;i<FFT_N;i+)部部分/求變換后結(jié)果的模值，存入復(fù)數(shù)的實si.real=sqrt(si.real*si.real+si.imag*si.imag);while(1);#include <iom128.h>#include <intrinsics.h>/*快速福利葉變換C程序包函數(shù)簡介：此程序包是通用的快速傅里葉變換C語言函數(shù)，移植性強，以下部分不依賴硬件。此程序

8、包采用聯(lián)合體的形式表示一個復(fù)數(shù)，輸入為自然順序的復(fù)數(shù)(輸入實數(shù)是可令復(fù)數(shù)虛部為0)，輸出為經(jīng)過FFT變換的自然順序的復(fù)數(shù).此程序包可在初始化時調(diào)用create_sin_tab()函數(shù)創(chuàng)建正弦函數(shù)表，以后的可采用查表法計算耗時較多的sin和cos運算，加快可計算速度使用說明：使用此函數(shù)只需更改宏定義FFT_N的值即可實現(xiàn)點數(shù)的改變， FFT_N的應(yīng)該為2的N次方，不滿足此條件時應(yīng)在后面補0。若使用查表法計算sin值和cos值，應(yīng)在調(diào)用 FFT函數(shù)前調(diào)用create_sin_tab()函數(shù)創(chuàng)建正弦表函數(shù)調(diào)用：FFT(s);時間：2010-2-20版本：Ver1.1參考文獻：*/#include&l

9、t;math.h>#define FFT_N 128定義福利葉變換的點數(shù)#define PI 3.1415926535897932384626433832795028841971定義圓周率值struct compx float real,imag;定義一個復(fù)數(shù)結(jié)構(gòu)struct compx sFFT_N;/FFT輸入和輸出：從 S0開始存放，根據(jù)大小自己定義/*函數(shù)原型函數(shù)功能輸入?yún)?shù)輸出參數(shù)float SIN_TABFFT_N/2;/定義正弦表的存放空間 struct compx EE(struct compx b1,struct compx b2)對兩個復(fù)數(shù)進行乘法運算兩個以聯(lián)合體定義

10、的復(fù)數(shù)a,ba和b的乘積，以聯(lián)合體的形式輸出*/ 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);/*函數(shù)原型：void create_sin_tab(float *sin_t)函數(shù)功能：創(chuàng)建一個正弦采樣表，采樣點數(shù)與福利葉變換點數(shù)相同 輸入?yún)?shù)：*sin_t存放正弦表的數(shù)組指針輸出參數(shù)：無*/ void create_sin_tab(float *sin_t)

11、int i;for(i=0;i<FFT_N/2;i+)sin_ti=sin(2*PI*i/FFT_N);void sin_tab(float pi)采用查表的方法計算一個數(shù)的正弦值pi所要計算正弦值弧度值，范圍 0-2*PI ,不滿足時需要轉(zhuǎn)換 輸入值 pi的正弦值/*函數(shù)原型函數(shù)功能輸入?yún)?shù)輸出參數(shù)*/ float sin_tab(float pi)int n;float a;n=(int)(pi*FFT_N/2/PI);if(n>=0&&n<FFT_N/2)a=SIN_TABn;else if(n>=FFT_N/2&&n<FFT

12、_N)a=-SIN_TABn-FFT_N/2;return a;void cos_tab(float pi)采用查表的方法計算一個數(shù)的余弦值pi所要計算余弦值弧度值，范圍 0-2*PI ,不滿足時需要轉(zhuǎn)換輸入值 pi的余弦值/*函數(shù)原型函數(shù)功能輸入?yún)?shù)輸出參數(shù)*/ float cos_tab(float pi) float a,pi2;pi2=pi+PI/2;if(pi2>2*PI)pi2-=2*PI;a=sin_tab(pi2);return a;*函數(shù)原型：void FFT(struct compx *xin,int N)函數(shù)功能：對輸入的復(fù)數(shù)組進行快速傅里葉變換(FFT)輸入?yún)?shù)：

13、*xin復(fù)數(shù)結(jié)構(gòu)體組的首地址指針，struct型輸出參數(shù)：無 * 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=xinj;xinj=xini;xini=t;k=nv2;/求j的下一個倒位序while(k<=j)如果k<=j,表示j的最高位為1(j=j-k;/把最高位變成0k=k/2;/k/2，比

14、較次高位，依次類推，逐個比較，直到某個位為0j=j+k;把0改為1(int le,lei,ip;f=FFT_N;for(l=1;(f=f/2)!=1;l+)Jfor(m=1;m<=l;m+)(le=2<<(m-1);lei=le/2;u.real=1.0;u.imag=0.0;/w.real=cos(PI/lei);/ w.imag=-sin(PI/lei);w.real=cos_tab(PI/lei);w.imag=-sin_tab(PI/lei);/FFT運算核，使用蝶形運算完成FFT運算/計算l的值，即計算蝶形級數(shù)/控制蝶形結(jié)級數(shù)/m表示第m級蝶形，l為蝶形級總數(shù)l=l

15、og (2) N/le蝶形結(jié)距離，即第 m級蝶形的蝶形結(jié)相距l(xiāng)e點同一蝶形結(jié)中參加運算的兩點的距離/u為蝶形結(jié)運算系數(shù)，初始值為1不適用查表法計算 sin值和cos值/w為系數(shù)商，即當前系數(shù)與前一個系數(shù)的商for(j=0;j<=lei-1;j+)for(i=j;i<=FFT_N-1;i=i+le)ip=i+lei;t=EE(xinip,u);xinip.real=xini.real-t.real;xinip.imag=xini.imag-t.imag;xini.real=xini.real+t.real;xini.imag=xini.imag+t.imag;u=EE(u,w);/控

16、制計算不同種蝶形結(jié)，即計算系數(shù)不同的蝶形結(jié)控制同一蝶形結(jié)運算，即計算系數(shù)相同蝶形結(jié)/i, ip分別表示參加蝶形運算的兩個節(jié)點蝶形運算，詳見公式改變系數(shù)，進行下一個蝶形運算void main()測試FFT變換，演示函數(shù)使用方法無無/*函數(shù)原型函數(shù)功能輸入?yún)?shù)輸出參數(shù)*/ void main() int i;create_sin_tab(SIN_TAB);/給結(jié)構(gòu)體賦值for(i=0;i<FFT_N;i+) si.real=sin(2*3.141592653589793*i/FFT_N); / 實部為正弦波 FFT_N 點采樣，賦值為 si.imag=0;/虛部為 0FFT(s);/進行快速

17、福利葉變換for(i=0;i<FFT_N;i+)求變換后結(jié)果的模值，存入復(fù)數(shù)的實部部分si.real=sqrt(si.real*si.real+si.imag*si.imag);while(1);/*快速福利葉變換C程序包函數(shù)簡介：此程序包是通用的快速傅里葉變換C語言函數(shù)，移植性強，以下部分不依賴硬件。此程序包采用聯(lián)合體的形式表示一個復(fù)數(shù)，輸入為自然順序的復(fù)數(shù)(輸入實數(shù)是可令復(fù)數(shù)虛部為0)，輸出為經(jīng)過FFT變換的自然順序的復(fù)數(shù).此程序包可在初始化時調(diào)用create_sin_tab()函數(shù)創(chuàng)建正弦函數(shù)表，以后的可采用查表法計算耗時較多的sin和cos運算，加快可計算速度.與Ver1.1版相

18、比較，Ver1.2版在創(chuàng)建正弦表時只建立了1/4個正弦波的采樣值，相比之下節(jié)省了FFT_N/4個存儲空間使用說明：使用此函數(shù)只需更改宏定義 FFT_N的值即可實現(xiàn)點數(shù)的改變，F(xiàn)FT_N的應(yīng)該為2的N次方，不滿足此條件時應(yīng)在后面補0。若使用查表法計算sin值和cos值，應(yīng)在調(diào)用FFT函數(shù)前調(diào)用create_sin_tab()函數(shù)創(chuàng)建正弦表函數(shù)調(diào)用：FFT(s);時間：2010-2-20版本：Ver1.2參考文獻：*/#include<math.h>#define FFT_N 128定義福利葉變換的點數(shù)#define PI 3.14159265358979323846264338327

19、95028841971/ 定義圓周率值struct compx float real,imag;定義一個復(fù)數(shù)結(jié)構(gòu)struct compx sFFT_N;/FFT輸入和輸出：從 S0開始存放，根據(jù)大小自己定義float SIN_TABFFT_N/4+1;定義正弦表的存放空間/*函數(shù)原型:struct compx EE(struct compx b1,struct compx b2)函數(shù)功能：對兩個復(fù)數(shù)進行乘法運算輸入?yún)?shù)：兩個以聯(lián)合體定義的復(fù)數(shù)a,b輸出參數(shù)：a和b的乘積，以聯(lián)合體的形式輸出*/struct compx EE(struct compx a,struct compx b)struc

20、t compx c;c.real=a.real*b.real-a.imag*b.imag;c.imag=a.real*b.imag+a.imag*b.real;return(c);/*函數(shù)原型函數(shù)功能輸入?yún)?shù)輸出參數(shù)void create_sin_tab(float *sin_t)創(chuàng)建一個正弦采樣表，采樣點數(shù)與福利葉變換點數(shù)相同*sin_t存放正弦表的數(shù)組指針無*/ void create_sin_tab(float *sin_t)(int i;for(i=0;i<=FFT_N/4;i+)sin_ti=sin(2*PI*i/FFT_N);/*函數(shù)原型函數(shù)功能輸入?yún)?shù)輸出參數(shù)void si

21、n_tab(float pi)采用查表的方法計算一個數(shù)的正弦值pi所要計算正弦值弧度值，范圍 0-2*PI ,不滿足時需要轉(zhuǎn)換 輸入值 pi的正弦值*/ float sin_tab(float pi)(int n;float a;n=(int)(pi*FFT_N/2/PI);if(n>=0&&n<=FFT_N/4) a=SIN_TABn;else if(n>FFT_N/4&&n<FFT_N/2)(n-=FFT_N/4;a=SIN_TABFFT_N/4-n;else if(n>=FFT_N/2&&n<3*FFT_

22、N/4)(n-=FFT_N/2;a=-SIN_TABn;else if(n>=3*FFT_N/4&&n<3*FFT_N)(n=FFT_N-n;a=-SIN_TABn;return a;void cos_tab(float pi)采用查表的方法計算一個數(shù)的余弦值pi所要計算余弦值弧度值，范圍 0-2*PI ,不滿足時需要轉(zhuǎn)換輸入值 pi的余弦值void FFT(struct compx *xin,int N)對輸入的復(fù)數(shù)組進行快速傅里葉變換(FFT)*xin復(fù)數(shù)結(jié)構(gòu)體組的首地址指針，struct型無nv2=FFT_N/2;法nm1=FFT_N-1;for(i=0;i&

23、lt;nm1;i+)(if(i<j)(t=xinj;xinj=xini;xini=t;k=nv2;while(k<=j)(j=j-k;k=k/2;j=j+k;/*函數(shù)原型函數(shù)功能輸入?yún)?shù)輸出參數(shù)*/ float cos_tab(float pi)(float a,pi2;pi2=pi+PI/2;if(pi2>2*PI)pi2-=2*PI;a=sin_tab(pi2);return a;/*函數(shù)原型函數(shù)功能輸入?yún)?shù)輸出參數(shù)*/ void FFT(struct compx *xin) (int f,m,nv2,nm1,i,k,l,j=0;struct compx u,w,t;變址運算，即把自然順序變成倒位序，采用雷德算如果i<j,即進行變址/求j的下一個倒位序如果k<=j,表示j的最高位為1把最高位變成0/k/2 ,比較次高位，依次類推，逐個比較，直到某個位為把0改為1(int le,lei,ip;f=FFT_N;for(l=1;(f=f/2)!=1;l+)；for(m=1;m<=l;m+)(l=log (2) Nle=2<<(m-1);相距l(xiāng)e點lei=le/2;u.real=1.0;u.imag=0.0;/w.real=cos(PI/lei);/ w.imag=-sin(PI/lei);w.