(快速傅里葉變換)C語(yǔ)言程序

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

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

3、 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)函數(shù)功能：對(duì)輸入的復(fù)數(shù)組進(jìn)行快速傅里葉變換（ 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;法

4、變址運(yùn)算，即把自然順序變成倒位序，采用雷德算nm1=FFT_N-1;for(i=0;i<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;如果i<j,即進(jìn)行變址求j的下一個(gè)倒位序如果k<=j,表示j的最高位為1把最高位變成0/k/2，比較次高位，依次類推，逐個(gè)比較，直到某個(gè)位為把0改為1int le,lei,ip;f=FFT_N;for(l=1;(f=f/2)!=1;l+);/FFT運(yùn)算核，使用蝶形運(yùn)算完成FFT運(yùn)算計(jì)算l的值，即計(jì)算蝶形級(jí)數(shù)for(m=1;m<=

5、l;m+)l=log Nle=2<<(m-1);相距l(xiāng)e點(diǎn)lei=le/2;u.real=1.0;u.imag=0.0;w.real=cos(PI/lei);w.imag=-sin(PI/lei); for(j=0;j<=lei-1;j+)/控制蝶形結(jié)級(jí)數(shù)/m表示第m級(jí)蝶形，l為蝶形級(jí)總數(shù)/le蝶形結(jié)距離，即第 m級(jí)蝶形的蝶形結(jié)同一蝶形結(jié)中參加運(yùn)算的兩點(diǎn)的距離/u為蝶形結(jié)運(yùn)算系數(shù)，初始值為1/w為系數(shù)商，即當(dāng)前系數(shù)與前一個(gè)系數(shù)的商控制計(jì)算不同種蝶形結(jié)，即計(jì)算系數(shù)不同的蝶形結(jié)for(i=j;i<=FFT_N-1;i=i+le)結(jié)ip=i+lei;t=EE(xinip,u)

6、;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é)運(yùn)算，即計(jì)算系數(shù)相同蝶形/i, ip分別表示參加蝶形運(yùn)算的兩個(gè)節(jié)點(diǎn) 蝶形運(yùn)算，詳見(jiàn)公式改變系數(shù)，進(jìn)行下一個(gè)蝶形運(yùn)算/*函數(shù)原型函數(shù)功能輸入?yún)?shù)輸出參數(shù)void main()測(cè)試FFT變換，演示函數(shù)使用方法無(wú)無(wú)*/void main()int i;for(i=0;i<FFT_N;i+)si.real=sin(2*3/p>

7、9793*i/FFT_N); / si.imag=0;/給結(jié)構(gòu)體賦值實(shí)部為正弦波FFT_N點(diǎn)采樣，賦值為1/虛部為0FFT(s)；/進(jìn)行快速傅立葉變換for(i=0;i<FFT_N;i+)部部分/求變換后結(jié)果的模值，存入復(fù)數(shù)的實(shí)si.real=sqrt(si.real*si.real+si.imag*si.imag);while(1);#include <iom128.h>#include <intrinsics.h>/*快速傅立葉變換C程序包函數(shù)簡(jiǎn)介：此程序包是通用的快速傅里葉變換C語(yǔ)言函數(shù)，移植性強(qiáng)，以下部分不依賴硬件。此程序包采用聯(lián)合體的形式表示一個(gè)復(fù)數(shù)，輸

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

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

10、以聯(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ù)功能輸入?yún)?shù)輸出參數(shù)創(chuàng)建一個(gè)正弦采樣表，采樣點(diǎn)數(shù)與傅立葉變換點(diǎn)數(shù)相同*sin_t存放正弦表的數(shù)組指針無(wú)*/void create_sin_tab(float *sin_t)int i;for(i=0;i<

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

12、_N/2;return a;/*函數(shù)原型函數(shù)功能輸入?yún)?shù)輸出參數(shù)void cos_tab(float pi)采用查表的方法計(jì)算一個(gè)數(shù)的余弦值pi所要計(jì)算余弦值弧度值，范圍 0-2*PI ,不滿足時(shí)需要轉(zhuǎn)換 輸入值 pi的余弦值*/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ù)功能：對(duì)輸入的復(fù)數(shù)組進(jìn)行快速傅里葉變換(FFT)輸入?yún)?shù)：*xin復(fù)數(shù)結(jié)構(gòu)體組的首地址指針，s

13、truct型輸出參數(shù)：無(wú)*/void FFT(struct compx *xin)int f,m,nv2,nm1,i,k,l,j=0;struct compx u,w,t;nv2=FFT_N/2;變址運(yùn)算，即把自然順序變成倒位序，采用雷德算法nm1=FFT_N-1;for(i=0;i<nm1;i+)if(i<j)如果i<j,即進(jìn)行變址t=xinj;xinj=xini;xini=t;k=nv2;求j的下一個(gè)倒位序while(k<=j)如果k<=j,表示j的最高位為1j=j-k;把最高位變成0k=k/2;/k/2，比較次高位，依次類推，逐個(gè)比較，直到某個(gè)位為 0j=j

14、+k;把0改為1/FFT運(yùn)算核，使用蝶形運(yùn)算完成FFT運(yùn)算計(jì)算l的值，即計(jì)算蝶形級(jí)數(shù)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);/控制蝶形結(jié)級(jí)數(shù)/m表示第m級(jí)蝶形，l為蝶形級(jí)總數(shù)l=log (2) N/le蝶形結(jié)距離，即第 m級(jí)蝶形的蝶形結(jié)

15、相距l(xiāng)e點(diǎn)/同一蝶形結(jié)中參加運(yùn)算的兩點(diǎn)的距離/u為蝶形結(jié)運(yùn)算系數(shù)，初始值為1/不適用查表法計(jì)算sin值和cos值/w為系數(shù)商，即當(dāng)前系數(shù)與前一個(gè)系數(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);/控制計(jì)算不同種蝶形結(jié)，即計(jì)算系數(shù)不同的蝶形結(jié)控制同一蝶形

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

17、的模值，存入復(fù)數(shù)的實(shí)部部分si.real=sqrt(si.real*si.real+si.imag*si.imag);while(1);#include <iom128.h>#include <intrinsics.h>/*快速傅立葉變換C程序包函數(shù)簡(jiǎn)介：此程序包是通用的快速傅里葉變換C語(yǔ)言函數(shù)，移植性強(qiáng)，以下部分不依賴硬件。此程序包采用聯(lián)合體的形式表示一個(gè)復(fù)數(shù)，輸入為自然順序的復(fù)數(shù)(輸入實(shí)數(shù)是可令復(fù)數(shù)虛部為0),輸出為經(jīng)過(guò)FFT變換的自然順序的復(fù)數(shù).此程序包可在初始化時(shí)調(diào)用create_sin_tab()函數(shù)創(chuàng)建正弦函數(shù)表，以后的可采用查表法計(jì)算耗時(shí)較多的sin和co

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

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

20、,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ù)原型函數(shù)功能輸入?yún)?shù)輸出參數(shù)/*void create_sin_tab(float *sin_t)創(chuàng)建一個(gè)正弦采樣表，采樣點(diǎn)數(shù)與傅立葉變換點(diǎn)數(shù)相同*sin_t存放正弦表的數(shù)組指針無(wú)*/ 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ù)原型

21、函數(shù)功能輸入?yún)?shù)輸出參數(shù)/*void sin_tab(float pi)采用查表的方法計(jì)算一個(gè)數(shù)的正弦值pi所要計(jì)算正弦值弧度值，范圍 0-2*PI ,不滿足時(shí)需要轉(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

22、&&n<3*FFT_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;/*函數(shù)原型函數(shù)功能輸入?yún)?shù)輸出參數(shù)void cos_tab(float pi)采用查表的方法計(jì)算一個(gè)數(shù)的余弦值pi所要計(jì)算余弦值弧度值，范圍 0-2*PI ,不滿足時(shí)需要轉(zhuǎn)換 輸入值 pi的余弦值 */ float cos_tab(float pi) float a,pi2;pi2=pi+PI/2;if(pi2>2*PI)pi2-=2*P

23、I;a=sin_tab(pi2); return a;/* 函數(shù)原型：void FFT(struct compx *xin,int N)函數(shù)功能：對(duì)輸入的復(fù)數(shù)組進(jìn)行快速傅里葉變換( FFT) 輸入?yún)?shù)：*xin復(fù)數(shù)結(jié)構(gòu)體組的首地址指針，struct型輸出參數(shù)：無(wú)*/ 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)t=xinj;xinj=xini;xini=t;k=nv2;while(k<=j)j=j-k;k=k/2;j=j+k;變址運(yùn)算，即把自然順序變成倒位序，采用雷德算/如果i<j,即進(jìn)行變址求j的下一個(gè)倒位序/如果k<=j,表示j的最高位為1把最高位變成0/k/2，比較次高位，依次類推，逐個(gè)比較，直到某個(gè)位為把0改為1int le,lei,ip;f=FFT_N;for(l=1;(f=f/2)!=1;l+)；for(m=1;m<=l;m+)l=log Nle=2<<(m-1);相距l(xiāng)e點(diǎn)lei=le/2;u.real=1.0;u.imag=0.0;w.real=cos(PI/lei);/ w.imag=-sin(PI/