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

1、快速傅里葉變換FFT的C語(yǔ)言實(shí)現(xiàn)及應(yīng)用快速傅里葉變換簡(jiǎn)介計(jì)算離散傅里葉變換的一種快速算法，簡(jiǎn)稱(chēng)FFT?？焖俑道锶~變換是 1965年由J.W.庫(kù)利和T.W.圖基提出的。采用這種算法能使計(jì)算機(jī)計(jì)算離散傅里葉變換所需要的乘法次 數(shù)大為減少，特別是被變換的抽樣點(diǎn)數(shù) N越多,FFT算法計(jì)算量的節(jié)省就越顯著。有限長(zhǎng)序列可以通過(guò)離散傅里葉變換(DFT)將其頻域也離散化K1 墾事快速傅里葉變換成有限長(zhǎng)序列但其計(jì)算量太大，很難實(shí)時(shí)地處理問(wèn)題，因此引出了快速傅里葉變換(FFT). 1965年，Cooley和Tukey提出了計(jì)算離散傅里葉變換(DFT)的快速算法，將DFT的運(yùn)算量減少了幾個(gè)數(shù)量級(jí)。從此，對(duì)快速傅里葉

2、變換(FFT)算法的研究便不斷深入，數(shù)字信號(hào)處理這門(mén)新興學(xué)科也隨FFT的出現(xiàn)和發(fā)展而迅速發(fā)展。根據(jù)對(duì)序列分解與選取方法的不同而產(chǎn)生了FFT的多種算法，基本算法是基2DIT和基2 DIF。FFT在離散傅里葉反變換、線(xiàn)性卷積和線(xiàn)性相關(guān)等方面也有重要應(yīng)用?？焖俑凳献儞Q(FFT),是離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的奇、偶、虛、實(shí)等特性，對(duì)離散傅立葉變換的算法進(jìn)行改進(jìn)獲得的。它對(duì)傅氏變換的理論并沒(méi)有新的發(fā)現(xiàn)，但是對(duì)于在計(jì)算機(jī)系統(tǒng)或者說(shuō)數(shù)字系統(tǒng)中應(yīng)用離散傅立葉變換，可以說(shuō)是進(jìn)了一大步。設(shè)快速傅里葉變換x(n)為N項(xiàng)的復(fù)數(shù)序列，由DFT變換，任一 X ( m)的計(jì)算都需要N次復(fù)數(shù)乘法和N-1次

3、復(fù)數(shù)加法，而一次復(fù)數(shù)乘法等于四次實(shí)數(shù)乘法和兩次實(shí)數(shù)加法，一次復(fù)數(shù)加 法等于兩次實(shí)N21M )(WN=Z x(2r Wrkk 22 rkWN-X2rWNr =0r=0快速傅里葉變換（四次實(shí)數(shù)乘法數(shù)加法，即使把一次復(fù)數(shù)乘法和一次復(fù)數(shù)加法定義成一次“運(yùn)算 和四次實(shí)數(shù)加法），那么求出 N項(xiàng)復(fù)數(shù)序列的 X（ m）,即N點(diǎn)DFT變換大約就需要NA2次運(yùn)算。當(dāng)N=1024點(diǎn)甚至更多的時(shí)候，需要N2=1048576次運(yùn)算，在 FFT中，利用WN的周期性和對(duì)稱(chēng)性，把一個(gè)N項(xiàng)序列（設(shè)N=2k,k為正整數(shù)），分為兩個(gè) N/2項(xiàng)的子序列，每個(gè) N/2點(diǎn)DFT變換需要（N/2）2次運(yùn)算，再用 N次運(yùn)算把兩個(gè)N/2點(diǎn)的D

4、FT變換組合成一個(gè)N點(diǎn)的DFT變換。這樣變換以后，總的運(yùn)算次數(shù)就變成N+2（ N/2）2=N+N2/2。繼續(xù)上面的例子，N=1024時(shí)，總的運(yùn)算次數(shù)就變成了525312次，節(jié)省了大約50%的運(yùn)算量。而如果我們將這種“一分為二”的思想不斷進(jìn)行下去，直到分成兩兩一組的DFT運(yùn)算單元，那么N點(diǎn)的DFT變換就只需要 Nlog2N次的運(yùn)算，N在1024點(diǎn)時(shí)，運(yùn)算量?jī)H有10240次，是先前的直接算法的1%,點(diǎn)數(shù)越多，運(yùn)算量的節(jié)約就越大，這就是FFT的優(yōu)越性。FFT算法的基本原理FFT 算法的基本思想：利用DFT系數(shù)的特性，合并DFT運(yùn)算中的某些項(xiàng)， 吧長(zhǎng)序列的 DFT-短序列的DFT，從而減少其運(yùn)算量。F

5、FT 算法分類(lèi)：時(shí) 間 抽選法 DIT: Decimation-ln-Time ;頻率抽選法 DIF: Decimati on-I n-F reque ncy 按時(shí)間抽選的基-2FFT算法1、算法原理設(shè)序列點(diǎn)數(shù)N = 2L , L為整數(shù)。若不滿(mǎn)足，則補(bǔ)零。N為2的整數(shù)幕的FFT算法稱(chēng)基-2FFT算法。將序列x（n）按n 的奇偶分成兩組：則 x(n)的 DFT:r 二 0,1,.,N -12N -1NN -4X k x n W,kx n W,k' x n W,kn£n=0 n為偶數(shù)n=0 n為奇數(shù)其丿NtVi(2)0點(diǎn)V(7)心叭ik 0辛點(diǎn)門(mén)0.140)0)*7 Y(0)at

6、i( 1 )=x(2) .ri(2)=AX4) OFT2 1 "xnn 期性求iL|XidU的后 2222X, k ,X2 k是以N為周期的X! k N =X, kX2 k N = X2 k2 2k N N 又Wn 2 二 Wn2W,八w.ckX(k) =Xdk) WiX2(k)丄NkX(k ) =Xi(k) -WiX2(k)Xi(A)+»<t X2(A)Vi(A) - W v X2(k)3、算法特點(diǎn)1)原位計(jì)算分解后的運(yùn)算量:個(gè)Nd點(diǎn)DFT復(fù)數(shù)乘法(JV / 2嚴(yán)復(fù)數(shù)加法JV/2(JV/2-1)兩個(gè)砒2點(diǎn)DFTAW2N(N/2 1)一個(gè)蝶形12AV 2個(gè)蝶形N/2N

7、總計(jì)N2/2 + N/2N2/2川(山2-1) + N«JV2/2運(yùn)算量減少了近一半!2、運(yùn)算量當(dāng)N = 2L時(shí)，共有 每復(fù)數(shù)乘法：L級(jí)蝶形，每級(jí) N / 2個(gè)蝶形,zNLF 2N log2 N復(fù)數(shù)加法:aF 二 NL 二 Nlog2N比較DFTlog 2 N2 Nlog 2 NmF (DFT ) m f ( FFT )22Xm(k) = Xm/(k)+Xm(k + 2)WnXm(k +2m)=Xmv(k)- Xm/(k +2m)wN蝶形運(yùn)算兩節(jié)點(diǎn)的第一個(gè)節(jié)點(diǎn)為k值，表示成L位二進(jìn)制數(shù)，左移 L -m位,把右邊空出的位置補(bǔ)零，結(jié)果為r的二進(jìn)制數(shù)。；Xm(k)=Xm(k)+Xm(j)W

8、N.XmUWXm/kLXm/jMNXi 打(/) *圖±1按算結(jié)構(gòu)2Jr00000000100410010102201011063Oil0011410010155101Oil36110111771113)蝶形運(yùn)算對(duì)N = 2L點(diǎn)FFT，輸入倒位序，輸出自然 序，2m -第m級(jí)運(yùn)算每個(gè)蝶形的兩節(jié)點(diǎn)距離為wN的確定蝶形運(yùn)算兩節(jié)點(diǎn)的第一個(gè)節(jié)點(diǎn)為k值，表示成L位二進(jìn)制數(shù)，左移L -m位輸入序右邊空出的位置補(bǔ)個(gè)存儲(chǔ)果元r的二進(jìn)制數(shù)。系數(shù)一：N / 2個(gè)存儲(chǔ)單元4）存儲(chǔ)單元快速傅立葉變換的C語(yǔ)言實(shí)現(xiàn)方法有了傅立葉變換，我們可以從信號(hào)的頻域特征去分析信號(hào)。尤其在無(wú)線(xiàn)通信系統(tǒng) 中，傅里葉變換的重要性

9、就更加明顯了，無(wú)論是設(shè)計(jì)者還是測(cè)試工程師，在工作中都會(huì) 和傅立葉變換打交道。我們要衡量一個(gè)處理器有沒(méi)有足夠的能力來(lái)運(yùn)行FFT算法，根據(jù)以上的簡(jiǎn)單介紹可以得出以下兩點(diǎn):1. 處理器要在一個(gè)指令周期能完成乘和累加的工作，因?yàn)閺?fù)數(shù)運(yùn)算要多次查表相乘 才能實(shí)現(xiàn)。2. 間接尋址，可以實(shí)現(xiàn)增/減1個(gè)變址量，方便各種查表方法。FFT要對(duì)原始序列進(jìn) 行反序排列，處理器要有反序間接尋址的能力。F面為一份FFT （快速傅立葉變換）的源碼（基于C）/* 匚 * /#i nclude <stdio.h>#in clude <math.h>#in clude <stdlib.h>#d

10、efi ne N 1000 typedef struct double real;/* 實(shí)部 */ double img;/* 虛部 */ complex;voidfft(); /*快速傅里葉變換*/voidifft(); /*快速傅里葉逆變換*/voidini tW(); /*初始化變化核*/voidcha nge(); /*變址*/voidadd(complexcomplex ,complex*); /*復(fù)數(shù)加法*/voidmul(complex,complex ,complex*); /*復(fù)數(shù)乘法*/voidsub(complex,complex ,complex*); /*復(fù)數(shù)減法*/

11、voiddivi(complex,complex ,complex*);/*復(fù)數(shù)除法*/voidoutput(); /*輸出結(jié)果*/complex xN, *W;/*輸出序列的值 */int size_x=0;/*輸入序列的長(zhǎng)度，只限 2的N次方*/double PI;int mai n()int i,method;system("cls");Pl=atan(1)*4;/*pi 等于4乘以1.0的正切值*/prin tf("Please in put the size of x:n");/*輸入序列的長(zhǎng)度*/scan f("%d",

12、&size_x);prin tf("Please in put the data in xN:(such as:5 6)n");/*輸入序列對(duì)應(yīng)的值*/for(i=0;i<size_x;i+)scan f("%lf %lf', &xi.real, &xi.img);ini tW();/*選擇FFT或逆FFT運(yùn)算*/prin tf("Use FFT(0) or IFFT(1)?n");scan f("%d",&method);if(method=0)fft();elseifft()

13、;output。；system("pause");return 0;/*進(jìn)行基-2 FFT運(yùn)算*/void fft()int i=O,j=O,k=O,l=O;complex up,dow n, product;cha nge();for(i=0;i< log(size_x)/log(2) ;i+) /*一級(jí)蝶形運(yùn)算 */l=1<<i;for(j=0;j<size_x;j+= 2*l ) /*一組蝶形運(yùn)算 */for(k=0;k<l;k+) /*個(gè)蝶形運(yùn)算 */mul(xj+k+l,Wsize_x*k/2/l, &product);add

14、(xj+k,product，&up);sub(xj+k,product，& dow n);xj+k=up;xj+k+l=dow n;void ifft()int i=0,j=0,k=0,l=size_x;complex up,dow n;for(i=0;i< (int)( log(size_x)/log(2) );i+) /*一級(jí)蝶形運(yùn)算 */l/=2;for(j=0;j<size_x;j+= 2*l ) /*一組蝶形運(yùn)算 */for(k=0;k<l;k+) /*個(gè)蝶形運(yùn)算 */add(xj+k,xj+k+l,&up);up.real/=2;up.im

15、g/=2;sub(xj+k,xj+k+l, &dow n);dow n.real/=2;dow n.img/=2;divi(dow n, Wsize_x*k/2 川，& dow n);xj+k=up;xj+k+l=dow n;cha nge();/*初始化變化核*/void ini tW()int i;W=(complex *)malloc(sizeof(complex) * size_x);for(i=0;i<size_x;i+)Wi.real=cos(2*PI/size_x*i);Wi.img=-1*si n(2*PI/size_x*i);/*變址計(jì)算，將x(n)碼位

16、倒置*/void cha nge()complex temp;un sig ned short i=O,j=O,k=O;double t;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=xi;xi=xj;xj=temp;void output。/*輸出結(jié)果 */int i;prin tf("The result are as follows' n");for(i

17、=0;i<size_x;i+)prin tf("%.4f",xi.real);if(xi.img>=0.0001)prin tf("+%.4fjn",xi.img);else if(fabs(xi.img)<0.0001)prin tf("n");elseprin tf("%.4fjn",xi.img);void add(complex a,complex b,complex *c)c->real=a.real+b.real;c->img=a.img+b.img;void mul(co

18、mplex a,complex b,complex *c) c->real=a.real*b.real - a.img*b.img;c->img=a.real*b.img + a.img*b.real;void sub(complex a,complex b,complex *c)c->real=a.real-b.real;c->img=a.img-b.img;void divi(complex a,complex b,complex *c)c->real=( a.real*b.real+a.img*b.img )/(b.real*b.real+b.img*b.i

19、mg);c->img=( a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img); C:JSSOFTCTuTanbinwteMp. exePlease4PleasecinputinputthetheofdatainX：xNl ： (such as :!j 6>di63?812154UseFFT<0>or1FFT<1>?0Theresultareasfollows:31.3300*29:-5.0000+7.0000j-6.0000-18,W(J(J0jPlease input4Please inputX：A

20、.the data inm TN;such as:S 6>6378124UseFFT<0>op 1FFT<1>?are as followsThe result 7.7500+7.25B0J -1.5000-4.5000j -1.2500+1,7500jB.eeM+i.seeej快速傅立葉變換的發(fā)展前景快速傅立葉變換作為一種數(shù)學(xué)方法，已經(jīng)廣泛地應(yīng)用在幾乎所有領(lǐng)域的頻譜分析中，而且經(jīng)久不衰，因?yàn)樾盘?hào)處理方法沒(méi)有先進(jìn)和落后之分，只有經(jīng)典和現(xiàn)代之別，在實(shí)際 系統(tǒng)中用得最好的方法就是管用的方法。換句話(huà)說(shuō)，信號(hào)處理方法與應(yīng)用背景和目的的 貼近程度是衡量信號(hào)處理方法優(yōu)劣的唯一標(biāo)準(zhǔn)。FFT是快速傅利葉變換(Fast FourierTransform簡(jiǎn)稱(chēng)FFT)的英文縮寫(xiě),它在當(dāng)今科技世界中的應(yīng)用姻當(dāng)活躍,無(wú)論是在時(shí)間序列分析領(lǐng)域中，還是在我國(guó)剛剛興起的生物頻譜治療的研究與應(yīng)用中，都有著重要的作用。同時(shí)，它又是軟件實(shí)現(xiàn)數(shù)字濾波器的必備組成部分之一。快速傅立葉變換的應(yīng)用領(lǐng)域信號(hào)分析，包括濾波、數(shù)據(jù)壓縮、電力系統(tǒng)的監(jiān)控等；研究偏微分方程，比如求解熱力學(xué)方程的解時(shí)，把f(t)展開(kāi)為三角級(jí)數(shù)最為關(guān)鍵概率與統(tǒng)計(jì)，量子力學(xué)等學(xué)科。我們以快速傅里葉變換在信號(hào)分析某一方面為例稍微說(shuō)明一下應(yīng)用過(guò)程。利用快速傅里葉變換FFT將圖像信號(hào)從