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

------------------------------------------------------------------------C語言實(shí)現(xiàn)FFT(快速傅里葉變換)#include<iom128.h>#include<intrinsics.h>/*********************************************************************快速福利葉變換C函數(shù)函數(shù)簡(jiǎn)介：此函數(shù)是通用的快速傅里葉變換C語言函數(shù)，移植性強(qiáng)，以下部分不依賴硬件。此函數(shù)采用聯(lián)合體的形式表示一個(gè)復(fù)數(shù)，輸入為自然順序的復(fù)數(shù)（輸入實(shí)數(shù)是可令復(fù)數(shù)虛部為0），輸出為經(jīng)過FFT變換的自然順序的復(fù)數(shù)使用說明：使用此函數(shù)只需更改宏定義FFT_N的值即可實(shí)現(xiàn)點(diǎn)數(shù)的改變，F(xiàn)FT_N的應(yīng)該為2的N次方，不滿足此條件時(shí)應(yīng)在后面補(bǔ)0函數(shù)調(diào)用：FFT(s);時(shí)間：2010-2-20版本：Ver1.0參考文獻(xiàn)：**********************************************************************/#include<math.h>#definePI3.1415926535897932384626433832795028841971//定義圓周率值#defineFFT_N128//定義福利葉變換的點(diǎn)數(shù)structcompx{floatreal,imag;};//定義一個(gè)復(fù)數(shù)結(jié)構(gòu)structcompxs[FFT_N];//FFT輸入和輸出：從S[1]開始存放，根據(jù)大小自己定義/*******************************************************************函數(shù)原型：structcompxEE(structcompxb1,structcompxb2)函數(shù)功能：對(duì)兩個(gè)復(fù)數(shù)進(jìn)行乘法運(yùn)算輸入?yún)?shù)：兩個(gè)以聯(lián)合體定義的復(fù)數(shù)a,b輸出參數(shù)：a和b的乘積，以聯(lián)合體的形式輸出*******************************************************************/structcompxEE(structcompxa,structcompxb){structcompxc;c.real=a.real*b.real-a.imag*b.imag;c.imag=a.real*b.imag+a.imag*b.real;return(c);}/*****************************************************************函數(shù)原型：voidFFT(structcompx*xin,intN)函數(shù)功能：對(duì)輸入的復(fù)數(shù)組進(jìn)行快速傅里葉變換（FFT）輸入?yún)?shù)：*xin復(fù)數(shù)結(jié)構(gòu)體組的首地址指針，struct型*****************************************************************/voidFFT(structcompx*xin){intf,m,nv2,nm1,i,k,l,j=0;structcompxu,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=xin[j];xin[j]=xin[i];xin[i]=t;}k=nv2;//求j的下一個(gè)倒位序while(k<=j)//如果k<=j,表示j的最高位為1{j=j-k;//把最高位變成0k=k/2;//k/2，比較次高位，依次類推，逐個(gè)比較，直到某個(gè)位為0}j=j+k;//把0改為1}{intle,lei,ip;//FFT運(yùn)算核，使用蝶形運(yùn)算完成FFT運(yùn)算f=FFT_N;for(l=1;(f=f/2)!=1;l++)//計(jì)算l的值，即計(jì)算蝶形級(jí)數(shù);for(m=1;m<=l;m++)//控制蝶形結(jié)級(jí)數(shù){//m表示第m級(jí)蝶形，l為蝶形級(jí)總數(shù)l=log（2）Nle=2<<(m-1);//le蝶形結(jié)距離，即第m級(jí)蝶形的蝶形結(jié)相距l(xiāng)e點(diǎn)lei=le/2;//同一蝶形結(jié)中參加運(yùn)算的兩點(diǎn)的距離u.real=1.0;//u為蝶形結(jié)運(yùn)算系數(shù)，初始值為1u.imag=0.0;w.real=cos(PI/lei);//w為系數(shù)商，即當(dāng)前系數(shù)與前一個(gè)系數(shù)的商w.imag=-sin(PI/lei);for(j=0;j<=lei-1;j++)//控制計(jì)算不同種蝶形結(jié)，即計(jì)算系數(shù)不同的蝶形結(jié){for(i=j;i<=FFT_N-1;i=i+le)//控制同一蝶形結(jié)運(yùn)算，即計(jì)算系數(shù)相同蝶形結(jié){ip=i+lei;//i，ip分別表示參加蝶形運(yùn)算的兩個(gè)節(jié)點(diǎn)t=EE(xin[ip],u);//蝶形運(yùn)算，詳見公式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);//改變系數(shù)，進(jìn)行下一個(gè)蝶形運(yùn)算}}}}/************************************************************函數(shù)原型：voidmain()函數(shù)功能：測(cè)試FFT變換，演示函數(shù)使用方法輸入?yún)?shù)：無輸出參數(shù)：無************************************************************/voidmain(){inti;for(i=0;i<FFT_N;i++)//給結(jié)構(gòu)體賦值{s[i].real=sin(2*3.141592653589793*i/FFT_N);//實(shí)部為正弦波FFT_N點(diǎn)采樣，賦值為1s[i].imag=0;//虛部為0}FFT(s);//進(jìn)行快速福利葉變換for(i=0;i<FFT_N;i++)//求變換后結(jié)果的模值，存入復(fù)數(shù)的實(shí)部部分s[i].real=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);while(1);}

#include<iom128.h>#include<intrinsics.h>/*********************************************************************快速福利葉變換C程序包函數(shù)簡(jiǎn)介：此程序包是通用的快速傅里葉變換C語言函數(shù)，移植性強(qiáng)，以下部分不依賴硬件。此程序包采用聯(lián)合體的形式表示一個(gè)復(fù)數(shù)，輸入為自然順序的復(fù)數(shù)（輸入實(shí)數(shù)是可令復(fù)數(shù)虛部為0），輸出為經(jīng)過FFT變換的自然順序的復(fù)數(shù).此程序包可在初始化時(shí)調(diào)用create_sin_tab()函數(shù)創(chuàng)建正弦函數(shù)表，以后的可采用查表法計(jì)算耗時(shí)較多的sin和cos運(yùn)算，加快可計(jì)算速度使用說明：使用此函數(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>#defineFFT_N128//定義福利葉變換的點(diǎn)數(shù)#definePI3.1415926535897932384626433832795028841971//定義圓周率值structcompx{floatreal,imag;};//定義一個(gè)復(fù)數(shù)結(jié)構(gòu)structcompxs[FFT_N];//FFT輸入和輸出：從S[0]開始存放，根據(jù)大小自己定義floatSIN_TAB[FFT_N/2];//定義正弦表的存放空間/*******************************************************************函數(shù)原型：structcompxEE(structcompxb1,structcompxb2)函數(shù)功能：對(duì)兩個(gè)復(fù)數(shù)進(jìn)行乘法運(yùn)算輸入?yún)?shù)：兩個(gè)以聯(lián)合體定義的復(fù)數(shù)a,b輸出參數(shù)：a和b的乘積，以聯(lián)合體的形式輸出*******************************************************************/structcompxEE(structcompxa,structcompxb){structcompxc;c.real=a.real*b.real-a.imag*b.imag;c.imag=a.real*b.imag+a.imag*b.real;return(c);}/******************************************************************函數(shù)原型：voidcreate_sin_tab(float*sin_t)函數(shù)功能：創(chuàng)建一個(gè)正弦采樣表，采樣點(diǎn)數(shù)與福利葉變換點(diǎn)數(shù)相同輸入?yún)?shù)：*sin_t存放正弦表的數(shù)組指針輸出參數(shù)：無******************************************************************/voidcreate_sin_tab(float*sin_t){inti;for(i=0;i<FFT_N/2;i++)sin_t[i]=sin(2*PI*i/FFT_N);}/******************************************************************函數(shù)原型：voidsin_tab(floatpi)函數(shù)功能：采用查表的方法計(jì)算一個(gè)數(shù)的正弦值輸入?yún)?shù)：pi所要計(jì)算正弦值弧度值，范圍0--2*PI，不滿足時(shí)需要轉(zhuǎn)換輸出參數(shù)：輸入值pi的正弦值******************************************************************/floatsin_tab(floatpi){intn;floata;n=(int)(pi*FFT_N/2/PI);if(n>=0&&n<FFT_N/2)a=SIN_TAB[n];elseif(n>=FFT_N/2&&n<FFT_N)a=-SIN_TAB[n-FFT_N/2];returna;}/******************************************************************函數(shù)原型：voidcos_tab(floatpi)函數(shù)功能：采用查表的方法計(jì)算一個(gè)數(shù)的余弦值輸入?yún)?shù)：pi所要計(jì)算余弦值弧度值，范圍0--2*PI，不滿足時(shí)需要轉(zhuǎn)換輸出參數(shù)：輸入值pi的余弦值******************************************************************/floatcos_tab(floatpi){floata,pi2;pi2=pi+PI/2;if(pi2>2*PI)pi2-=2*PI;a=sin_tab(pi2);returna;}/*****************************************************************函數(shù)原型：voidFFT(structcompx*xin,intN)函數(shù)功能：對(duì)輸入的復(fù)數(shù)組進(jìn)行快速傅里葉變換（FFT）輸入?yún)?shù)：*xin復(fù)數(shù)結(jié)構(gòu)體組的首地址指針，struct型輸出參數(shù)：無*****************************************************************/voidFFT(structcompx*xin){intf,m,nv2,nm1,i,k,l,j=0;structcompxu,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=xin[j];xin[j]=xin[i];xin[i]=t;}k=nv2;//求j的下一個(gè)倒位序while(k<=j)//如果k<=j,表示j的最高位為1{j=j-k;//把最高位變成0k=k/2;//k/2，比較次高位，依次類推，逐個(gè)比較，直到某個(gè)位為0}j=j+k;//把0改為1}{intle,lei,ip;//FFT運(yùn)算核，使用蝶形運(yùn)算完成FFT運(yùn)算f=FFT_N;for(l=1;(f=f/2)!=1;l++)//計(jì)算l的值，即計(jì)算蝶形級(jí)數(shù);for(m=1;m<=l;m++)//控制蝶形結(jié)級(jí)數(shù){//m表示第m級(jí)蝶形，l為蝶形級(jí)總數(shù)l=log（2）Nle=2<<(m-1);//le蝶形結(jié)距離，即第m級(jí)蝶形的蝶形結(jié)相距l(xiāng)e點(diǎn)lei=le/2;//同一蝶形結(jié)中參加運(yùn)算的兩點(diǎn)的距離u.real=1.0;//u為蝶形結(jié)運(yùn)算系數(shù)，初始值為1u.imag=0.0;//w.real=cos(PI/lei);//不適用查表法計(jì)算sin值和cos值//w.imag=-sin(PI/lei);w.real=cos_tab(PI/lei);//w為系數(shù)商，即當(dāng)前系數(shù)與前一個(gè)系數(shù)的商w.imag=-sin_tab(PI/lei);for(j=0;j<=lei-1;j++)//控制計(jì)算不同種蝶形結(jié)，即計(jì)算系數(shù)不同的蝶形結(jié){for(i=j;i<=FFT_N-1;i=i+le)//控制同一蝶形結(jié)運(yùn)算，即計(jì)算系數(shù)相同蝶形結(jié){ip=i+lei;//i，ip分別表示參加蝶形運(yùn)算的兩個(gè)節(jié)點(diǎn)t=EE(xin[ip],u);//蝶形運(yùn)算，詳見公式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);//改變系數(shù)，進(jìn)行下一個(gè)蝶形運(yùn)算}}}}/************************************************************函數(shù)原型：voidmain()函數(shù)功能：測(cè)試FFT變換，演示函數(shù)使用方法輸入?yún)?shù)：無輸出參數(shù)：無************************************************************/voidmain(){inti;create_sin_tab(SIN_TAB);for(i=0;i<FFT_N;i++)//給結(jié)構(gòu)體賦值{s[i].real=sin(2*3.141592653589793*i/FFT_N);//實(shí)部為正弦波FFT_N點(diǎn)采樣，賦值為1s[i].imag=0;//虛部為0}FFT(s);//進(jìn)行快速福利葉變換for(i=0;i<FFT_N;i++)//求變換后結(jié)果的模值，存入復(fù)數(shù)的實(shí)部部分s[i].real=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);while(1);}

#include<iom128.h>#include<intrinsics.h>/*********************************************************************快速福利葉變換C程序包函數(shù)簡(jiǎn)介：此程序包是通用的快速傅里葉變換C語言函數(shù)，移植性強(qiáng)，以下部分不依賴硬件。此程序包采用聯(lián)合體的形式表示一個(gè)復(fù)數(shù)，輸入為自然順序的復(fù)數(shù)（輸入實(shí)數(shù)是可令復(fù)數(shù)虛部為0），輸出為經(jīng)過FFT變換的自然順序的復(fù)數(shù).此程序包可在初始化時(shí)調(diào)用create_sin_tab()函數(shù)創(chuàng)建正弦函數(shù)表，以后的可采用查表法計(jì)算耗時(shí)較多的sin和cos運(yùn)算，加快可計(jì)算速度.與Ver1.1版相比較，Ver1.2版在創(chuàng)建正弦表時(shí)只建立了1/4個(gè)正弦波的采樣值，相比之下節(jié)省了FFT_N/4個(gè)存儲(chǔ)空間使用說明：使用此函數(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.2參考文獻(xiàn)：**********************************************************************/#include<math.h>#defineFFT_N128//定義福利葉變換的點(diǎn)數(shù)#definePI3.1415926535897932384626433832795028841971 //定義圓周率值structcompx{floatreal,imag;}; //定義一個(gè)復(fù)數(shù)結(jié)構(gòu)structcompxs[FFT_N]; //FFT輸入和輸出：從S[0]開始存放，根據(jù)大小自己定義floatSIN_TAB[FFT_N/4+1]; //定義正弦表的存放空間/*******************************************************************函數(shù)原型：structcompxEE(structcompxb1,structcompxb2)函數(shù)功能：對(duì)兩個(gè)復(fù)數(shù)進(jìn)行乘法運(yùn)算輸入?yún)?shù)：兩個(gè)以聯(lián)合體定義的復(fù)數(shù)a,b輸出參數(shù)：a和b的乘積，以聯(lián)合體的形式輸出*******************************************************************/structcompxEE(structcompxa,structcompxb){structcompxc;c.real=a.real*b.real-a.imag*b.imag;c.imag=a.real*b.imag+a.imag*b.real;return(c);}/******************************************************************函數(shù)原型：voidcreate_sin_tab(float*sin_t)函數(shù)功能：創(chuàng)建一個(gè)正弦采樣表，采樣點(diǎn)數(shù)與福利葉變換點(diǎn)數(shù)相同輸入?yún)?shù)：*sin_t存放正弦表的數(shù)組指針輸出參數(shù)：無******************************************************************/voidcreate_sin_tab(float*sin_t){inti;for(i=0;i<=FFT_N/4;i++)sin_t[i]=sin(2*PI*i/FFT_N);}/******************************************************************函數(shù)原型：voidsin_tab(floatpi)函數(shù)功能：采用查表的方法計(jì)算一個(gè)數(shù)的正弦值輸入?yún)?shù)：pi所要計(jì)算正弦值弧度值，范圍0--2*PI，不滿足時(shí)需要轉(zhuǎn)換輸出參數(shù)：輸入值pi的正弦值******************************************************************/floatsin_tab(floatpi){intn;floata;n=(int)(pi*FFT_N/2/PI);if(n>=0&&n<=FFT_N/4)a=SIN_TAB[n];elseif(n>FFT_N/4&&n<FFT_N/2){n-=FFT_N/4;a=SIN_TAB[FFT_N/4-n];}elseif(n>=FFT_N/2&&n<3*FFT_N/4){n-=FFT_N/2;a=-SIN_TAB[n];}elseif(n>=3*FFT_N/4&&n<3*FFT_N){n=FFT_N-n;a=-SIN_TAB[n];}returna;}/******************************************************************函數(shù)原型：voidcos_tab(floatpi)函數(shù)功能：采用查表的方法計(jì)算一個(gè)數(shù)的余弦值輸入?yún)?shù)：pi所要計(jì)算余弦值弧度值，范圍0--2*PI，不滿足時(shí)需要轉(zhuǎn)換輸出參數(shù)：輸入值pi的余弦值******************************************************************/floatcos_tab(floatpi){floata,pi2;pi2=pi+PI/2;if(pi2>2*PI)pi2-=2*PI;a=sin_tab(pi2);returna;}/*****************************************************************函數(shù)原型：voidFFT(structcompx*xin,intN)函數(shù)功能：對(duì)輸入的復(fù)數(shù)組進(jìn)行快速傅里葉變換（FFT）輸入?yún)?shù)：*xin復(fù)數(shù)結(jié)構(gòu)體組的首地址指針，struct型輸出參數(shù)：無*****************************************************************/voidFFT(structcompx*xin){intf,m,nv2,nm1,i,k,l,j=0;structcompxu,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=xin[j];xin[j]=xin[i];xin[i]=t;}k=nv2;//求j的下一個(gè)倒位序while(k<=j)//如果k<=j,表示j的最高位為1{j=j-k;//把最高位變成0k=k/2;//k/2，比較次高位，依次類推，逐個(gè)比較，直到某個(gè)位為0}j=j+k;//把0改為1}{intle,lei,ip;//FFT運(yùn)算核，使用蝶形運(yùn)算完成FFT運(yùn)算f=FFT_N;for(l=1;(f=f/2)!=1;l++)//計(jì)算l的值，即計(jì)算蝶形級(jí)數(shù);for(m=1;m<=l;m++)//控制蝶形結(jié)級(jí)數(shù){//m表示第m級(jí)蝶形，l為蝶形級(jí)總數(shù)l=log（2）Nle=2<<(m-1);//le蝶形結(jié)距離，即第m級(jí)蝶形的蝶形結(jié)相距l(xiāng)e點(diǎn)lei=le/2;