快速傅里葉變換FFT的C語言實現(xiàn)及應(yīng)用

快速傅里葉變換FFT的C語言實現(xiàn)及應(yīng)用快速傅里葉變換簡介計算離散傅里葉變換的一種快速算法,簡稱FFT?？焖俑道锶~變換是1965年由J.W.庫利和T.W.圖基提出的。采用這種算法能使計算機計算離散傅里葉變換所需要的乘法次數(shù)大為減少,特別是被變換的抽樣點數(shù)N越多,FFT算法計算量的節(jié)省就越顯著。有限長序列可以通過離散傅里葉變換(DFT)將其頻域也離散化

快速傅里葉變換成有限長序列.但其計算量太大,很難實時地處理問題,因此引出了快速傅里葉變換(FFT).1965年，Cooley和Tukey提出了計算離散傅里葉變換〔DFT〕的快速算法，將DFT的運算量減少了幾個數(shù)量級。從此，對快速傅里葉變換〔FFT〕算法的研究便不斷深入，數(shù)字信號處理這門新興學(xué)科也隨FFT的出現(xiàn)和開展而迅速開展。根據(jù)對序列分解與選取方法的不同而產(chǎn)生了FFT的多種算法，根本算法是基２DIT和基２DIF。FFT在離散傅里葉反變換、線性卷積和線性相關(guān)等方面也有重要應(yīng)用。快速傅氏變換〔FFT〕，是離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的奇、偶、虛、實等特性，對離散傅立葉變換的算法進行改良獲得的。它對傅氏變換的理論并沒有新的發(fā)現(xiàn)，但是對于在計算機系統(tǒng)或者說數(shù)字系統(tǒng)中應(yīng)用離散傅立葉變換，可以說是進了一大步。設(shè)

快速傅里葉變換x(n)為N項的復(fù)數(shù)序列，由DFT變換，任一X〔m〕的計算都需要N次復(fù)數(shù)乘法和N-1次復(fù)數(shù)加法，而一次復(fù)數(shù)乘法等于四次實數(shù)乘法和兩次實數(shù)加法，一次復(fù)數(shù)加法等于兩次實

快速傅里葉變換數(shù)加法，即使把一次復(fù)數(shù)乘法和一次復(fù)數(shù)加法定義成一次“運算”〔四次實數(shù)乘法和四次實數(shù)加法〕，那么求出N項復(fù)數(shù)序列的X〔m〕,即N點DFT變換大約就需要N^2次運算。當(dāng)N=1024點甚至更多的時候，需要N2=1048576次運算，在FFT中，利用WN的周期性和對稱性，把一個N項序列〔設(shè)N=2k,k為正整數(shù)〕，分為兩個N/2項的子序列，每個N/2點DFT變換需要〔N/2〕2次運算，再用N次運算把兩個N/2點的DFT變換組合成一個N點的DFT變換。這樣變換以后，總的運算次數(shù)就變成N+2〔N/2〕2=N+N2/2。繼續(xù)上面的例子，N=1024時，總的運算次數(shù)就變成了525312次，節(jié)省了大約50%的運算量。而如果我們將這種“一分為二”的思想不斷進行下去，直到分成兩兩一組的DFT運算單元，那么N點的DFT變換就只需要Nlog2N次的運算，N在1024點時，運算量僅有10240次，是先前的直接算法的1%，點數(shù)越多，運算量的節(jié)約就越大，這就是FFT的優(yōu)越性。FFT算法的根本原理FFT算法的根本思想：利用DFT系數(shù)的特性，合并DFT運算中的某些項，吧長序列的DFT—>短序列的DFT，從而減少其運算量。FFT算法分類:時間抽選法DIT:Decimation-In-Time；頻率抽選法DIF:Decimation-In-Frequency按時間抽選的基-2FFT算法1、算法原理設(shè)序列點數(shù)N=2L，L為整數(shù)。假設(shè)不滿足，那么補零。N為2的整數(shù)冪的FFT算法稱基-2FFT算法。將序列x(n)按n的奇偶分成兩組：那么x(n)的DFT:其中再利用周期性求X(k)的后半局部：復(fù)數(shù)乘法：復(fù)數(shù)乘法：當(dāng)N=2L時，共有L級蝶形，每級N/2個蝶形，每

個蝶形有1次復(fù)數(shù)乘法2次復(fù)數(shù)加法。2、運算量當(dāng)N=2L時，共有L級蝶形，每級N/2個蝶形，每

個蝶形有1次復(fù)數(shù)乘法2次復(fù)數(shù)加法。復(fù)數(shù)乘法：2、運算量2、運算量復(fù)數(shù)加法：比擬DFT比擬DFT3、算法特點3、算法特點1〕原位計算1〕原位計算2〕倒位序2〕倒位序3〕蝶形運算3〕蝶形運算對N=2L點FFT，輸入倒位序，輸出自然序，

第m級運算每個蝶形的兩節(jié)點距離為2m–1

第m級運算：對N=2L點FFT，輸入倒位序，輸出自然序，

第m級運算每個蝶形的兩節(jié)點距離為2m–1

第m級運算：蝶形運算兩節(jié)點的第一個節(jié)點為k值，表示成L位二進制數(shù)，左移L–m位，把右邊空出的位置補零，結(jié)果為r的二進制數(shù)。蝶形運算兩節(jié)點的第一個節(jié)點為k值，表示成L位二進制數(shù)，左移L–m位，把右邊空出的位置補零，結(jié)果為r的二進制數(shù)。蝶形運算兩節(jié)點的第一個節(jié)點為k值，表示成L位二進制數(shù)，左移L–m位，把右邊空出的位置補零，結(jié)果為r的二進制數(shù)。蝶形運算兩節(jié)點的第一個節(jié)點為k值，表示成L位二進制數(shù)，左移L–m位，把右邊空出的位置補零，結(jié)果為r的二進制數(shù)。4〕存儲單元4〕存儲單元輸入序列x(n):N個存儲單元系數(shù)：N/2個存儲單元快速傅立葉變換的C語言實現(xiàn)方法有了傅立葉變換，我們可以從信號的頻域特征去分析信號。尤其在無線通信系統(tǒng)中，傅里葉變換的重要性就更加明顯了，無論是設(shè)計者還是測試工程師，在工作中都會和傅立葉變換打交道。我們要衡量一個處理器有沒有足夠的能力來運行FFT算法，根據(jù)以上的簡單介紹可以得出以下兩點：處理器要在一個指令周期能完成乘和累加的工作，因為復(fù)數(shù)運算要屢次查表相乘才能實現(xiàn)。間接尋址，可以實現(xiàn)增/減1個變址量，方便各種查表方法。FFT要對原始序列進行反序排列，處理器要有反序間接尋址的能力。下面為一份FFT〔快速傅立葉變換〕的源碼〔基于C〕/************FFT***********/#include<stdio.h>#include<math.h>#include<stdlib.h>#defineN1000typedefstruct{doublereal;/*實部*/doubleimg;/*虛部*/}complex;voidfft();/*快速傅里葉變換*/voidifft();/*快速傅里葉逆變換*/voidinitW();/*初始化變化核*/voidchange();/*變址*/voidadd(complex,complex,complex*);/*復(fù)數(shù)加法*/voidmul(complex,complex,complex*);/*復(fù)數(shù)乘法*/voidsub(complex,complex,complex*);/*復(fù)數(shù)減法*/voiddivi(complex,complex,complex*);/*復(fù)數(shù)除法*/voidoutput();/*輸出結(jié)果*/complexx[N],*W;/*輸出序列的值*/intsize_x=0;/*輸入序列的長度，只限2的N次方*/doublePI;intmain(){inti,method;system("cls");PI=atan(1)*4;/*pi等于4乘以1.0的正切值*/printf("Pleaseinputthesizeofx:\n");/*輸入序列的長度*/scanf("%d",&size_x);printf("Pleaseinputthedatainx[N]:(suchas:56)\n");/*輸入序列對應(yīng)的值*/for(i=0;i<size_x;i++)scanf("%lf%lf",&x[i].real,&x[i].img);initW();/*選擇FFT或逆FFT運算*/printf("UseFFT(0)orIFFT(1)?\n");scanf("%d",&method);if(method==0)fft();elseifft();output();system("pause");return0;}/*進行基-2FFT運算*/voidfft(){inti=0,j=0,k=0,l=0;complexup,down,product;change();for(i=0;i<log(size_x)/log(2);i++)/*一級蝶形運算*/{l=1<<i;for(j=0;j<size_x;j+=2*l)/*一組蝶形運算*/{for(k=0;k<l;k++)/*一個蝶形運算*/{mul(x[j+k+l],W[size_x*k/2/l],&product);add(x[j+k],product,&up);sub(x[j+k],product,&down);x[j+k]=up;x[j+k+l]=down;}}}}voidifft(){inti=0,j=0,k=0,l=size_x;complexup,down;for(i=0;i<(int)(log(size_x)/log(2));i++)/*一級蝶形運算*/{l/=2;for(j=0;j<size_x;j+=2*l)/*一組蝶形運算*/{for(k=0;k<l;k++)/*一個蝶形運算*/{add(x[j+k],x[j+k+l],&up);up.real/=2;up.img/=2;sub(x[j+k],x[j+k+l],&down);down.real/=2;down.img/=2;divi(down,W[size_x*k/2/l],&down);x[j+k]=up;x[j+k+l]=down;}}}change();}/*初始化變化核*/voidinitW(){inti;W=(complex*)malloc(sizeof(complex)*size_x);for(i=0;i<size_x;i++){W[i].real=cos(2*PI/size_x*i);W[i].img=-1*sin(2*PI/size_x*i);}}/*變址計算，將x(n)碼位倒置*/voidchange(){complextemp;unsignedshorti=0,j=0,k=0;doublet;for(i=0;i<size_x;i++){k=i;j=0;t=(log(size_x)/log(2));while((t--)>0){j=j<<1;j|=(k&1);k=k>>1;}if(j>i){temp=x[i];x[i]=x[j];x[j]=temp;}}}voidoutput()/*輸出結(jié)果*/{inti;printf("Theresultareasfollows\n");for(i=0;i<size_x;i++){printf("%.4f",x[i].real);if(x[i].img>=0.0001)printf("+%.4fj\n",x[i].img);elseif(fabs(x[i].img)<0.0001)printf("\n");elseprintf("%.4fj\n",x[i].img);}}voidadd(complexa,complexb,complex*c){c->real=a.real+b.real;c->img=a.img+b.img;}voidmul(complexa,complexb,complex*c){c->real=a.real*b.real-a.img*b.img;c->img=a.real*b.img+a.img*b.real;}voidsub(complexa,complexb,complex*c){c->real=a.real-b.real;c->img=a.img-b.img;}voiddivi(complexa,complexb,complex*c){c->real=(a.real*b.real+a.img*b.img)/(b.real*b.real+b.img*b.img);c->img=(a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img);}快速傅立葉變換的開展前景快速傅立葉變換作為一種數(shù)學(xué)方法，已經(jīng)廣泛地應(yīng)用在幾乎所有領(lǐng)域的頻譜分析中，而且經(jīng)久不衰，因為信號處理方法沒有先進和落后之分，只有經(jīng)典和現(xiàn)代之別，在實際系統(tǒng)中用得最好的方法就是管用的方法。換句話說，信號處理方法與應(yīng)用背景和目的的貼近程度是衡量信號處理方法優(yōu)劣的唯一標(biāo)準(zhǔn)。FFT是快速傅利葉變換(FastFourierTransform簡稱FFT)的英文縮寫,它在當(dāng)今科技世界中的應(yīng)用姻當(dāng)活潑,無論是在時間序列分析領(lǐng)域中,還是在我國剛剛興起的生物頻譜治療的研究與應(yīng)用中,都有著重要的作用。同時,它又是軟件實現(xiàn)數(shù)字濾波器的必備組成局部之一?？焖俑盗⑷~變換的應(yīng)用領(lǐng)域

信號分析，包括濾波、數(shù)據(jù)壓縮、電力系統(tǒng)的監(jiān)控等；

研究偏微分方程，比方求解熱力學(xué)方程的解時，把f(t)展開為三角級數(shù)最為關(guān)鍵

概率與統(tǒng)計，量子力學(xué)等學(xué)科。

我們以快速傅里葉變換在信號分析某一方面為例稍微說明一下應(yīng)用過程。利用快速傅里葉變換FFT將圖像信號從空間轉(zhuǎn)換