C語言編寫FFT程序

1、DSP課程作業(yè)用C語言編寫FFT程序1,快速傅里葉變換FFT簡介快速傅氏變換（FFT）,是離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的奇、偶、虛、實(shí)等特性，對離散傅立葉變換的算法進(jìn)行改進(jìn)獲得的。它對傅氏變換的理論并沒有新的發(fā)現(xiàn)，但是對于在計(jì)算機(jī)系統(tǒng)或者說數(shù)字系統(tǒng)中應(yīng)用離散傅立葉變換，可以說是進(jìn)了一大步。我們假設(shè)x（n）為N項(xiàng)的復(fù)數(shù)序列，由DFT變換，任一X（M）的計(jì)算都需要N次復(fù)數(shù)乘法和N-1次復(fù)數(shù)加法，而一次復(fù)數(shù)乘法等于四次實(shí)數(shù)乘法和兩次實(shí)數(shù)加法，一次復(fù)數(shù)加法等于兩次實(shí)數(shù)加法，即使把一次復(fù)數(shù)乘法和一次復(fù)數(shù)加法定義成一次“運(yùn)算”（四次實(shí)數(shù)乘法和四次實(shí)數(shù)加法），那么求出N項(xiàng)復(fù)數(shù)序列的X（m），

2、即N點(diǎn)DFT變換大約就需要N，次運(yùn)算。當(dāng)N=1024點(diǎn)甚至更多的時(shí)候，需要N2=1048576次運(yùn)算，在FFT中，利用WN勺周期性和對稱性，把一個(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)的DFT變換組合成一個(gè)N點(diǎn)的DFT變換。這樣變換以后，總的運(yùn)算次數(shù)就變成N+2N/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

3、次的運(yùn)算，N在1024點(diǎn)時(shí)，運(yùn)算量僅有10240次，是先前的直接算法的1%點(diǎn)數(shù)越多，運(yùn)算量的節(jié)約就越大，這就是FFT的優(yōu)越性。2,FFT算法的基本原理FFT算法的基本思想：利用DFT系數(shù)的特性，合并DFT運(yùn)算中的某些項(xiàng)，吧長序列的DFT-短序列的DFT,從而減少其運(yùn)算量。FFT算法分類：時(shí)間抽選法DIT:Decimation-In-Time;頻率抽選法DIF:Decimation-In-Frequency按時(shí)間抽選的基-2FFT算法1、算法原理設(shè)序列點(diǎn)數(shù)N=2L,L為整數(shù)。若不滿足，則補(bǔ)零。N為2的整數(shù)哥的的奇偶分成兩組：y2XrAXIx2r1x2rFFT算法稱基-2FFT算法。將序列x(n)

4、按nN,r0,1,.,12則 x(n)的 DFT:X kN1N1N1xnWN1kxnWN1kxnWN1k21x2rWN2rkN-1x2r2r1k1WnNr021x1rr0N121x1rr02rkWN2W；%r0n2kWNkrN2W；rx201X202rkWn2rW%XikkWnX2k(r,k(kN 0,1,2 1)其中X1(k)Xi(k)x(r)Wkr021X1(r)WNrkNI21x(2r)WNrkr0N21x(2r)W；k再利用周期性求X(k)的后半部分:QX1k,X2k是吟為周期的X1kk又WnN2N工XikX2Wn2W：X2kX1(k) WnX2%)kXi(k) WNX2(k)X(k)

5、NX(k2)Vi(A)Xi(A)-V2(A)圖41時(shí)間抽選法蝶形運(yùn)算流圖符號n為偶數(shù)n為奇數(shù)V(0)V(2)DI TM(3)v1)x；( I )=a (3)V2(2)DFTa-(3)=a17)圖4-2按時(shí)間抽選分解后的運(yùn)算量:省數(shù)乘法笈數(shù)加法一個(gè)N/2點(diǎn)DFT(V/2)2N/2(N/2】)兩個(gè)2/2點(diǎn)DFTNq2M(N/21)個(gè)蝶形12-個(gè)蝶形N/2N總計(jì)/VJ/2+A/2力N工，2JV(AT/2-l)+JV運(yùn)算量減少了近半!2）、運(yùn)算量當(dāng)N=2L時(shí)，共有L級蝶形，每級N/2個(gè)蝶形,每比較DFT2mF(DFT)N22NmF(FFT)Nlog2NF7log2Ng23)、算法特點(diǎn)原位計(jì)算蝶形運(yùn)算兩

6、節(jié)點(diǎn)的第一個(gè)節(jié)點(diǎn)為k值，表示成L位二進(jìn)制數(shù)，左移L-m位,把右邊空出的位置補(bǔ)零，結(jié)果為r的二進(jìn)制數(shù)。Xm(k)Xmi(k)Xmi(j)WNXm(j)Xmi(k)Xmi(j)WNXMA)工(燈二Xa圖4-7按時(shí)間抽選蝶形運(yùn)算結(jié)構(gòu)倒位序?qū)=2L點(diǎn)FFT,輸入倒位序，輸出自然序，第m級運(yùn)算每個(gè)蝶形的兩節(jié)點(diǎn)距離為2m1倒位序自然序00000000100410010102201011063()11001141001015510101136110111771112一m1rXm(k)Xmi(k)Xmi(k2必Xm(k2m1)Xmi(k)Xmi(k2m1)W；Xa點(diǎn)）抬尸入入（幻+N”W7一心二居一（知一,

7、5隊(duì)圖4-7按時(shí)間抽選蝶形運(yùn)算結(jié)構(gòu)W；的確定蝶形運(yùn)算兩節(jié)點(diǎn)的第一個(gè)節(jié)點(diǎn)為k值，表示成L位二進(jìn)制數(shù)，左移L-m位，把右邊空出的位置補(bǔ)零，結(jié)果為r的二進(jìn)制數(shù)。存儲(chǔ)單元輸入序列x（n）:N個(gè)存儲(chǔ)單元系數(shù)W；:N/2個(gè)存儲(chǔ)單元3，快速傅立葉變換的C語言實(shí)現(xiàn)方法我們要衡量一個(gè)處理器有沒有足夠的能力來運(yùn)行FFT算法，根據(jù)以上的簡單介紹可以得出以下兩點(diǎn)：1 .處理器要在一個(gè)指令周期能完成乘和累加的工作，因?yàn)閺?fù)數(shù)運(yùn)算要多次查表相乘才能實(shí)現(xiàn)。2 .間接尋址，可以實(shí)現(xiàn)增/減1個(gè)變址量，方便各種查表方法。FFT要對原始序列進(jìn)行反序排列，處理器要有反序間接尋址的能力。#include#include#include#

8、defineN1000typedefstructdoublereal;doubleimg;complex;voidfft();voidifft();voidvoidinitW();/*初始化變化核*/change();voidadd(complex,complex,complex*);voidmul(complex,complex,complex*);voidsub(complex,complex,complex*);voidvoiddivi(complex,complex,complex*);output();complexxN,*W;/*輸出序列值*/intsize_x=0;/*輸入序列的

9、長度，只限2的N次方*/doublePI;intmain()inti,method;system(cls);PI=atan(1)*4;/*pi等于4乘以1.0正切值*/printf(Pleaseinputthesizeofx:n);/*輸入序列長度*/scanf(%d,&size_x);printf(PleaseinputthedatainxN:(suchas:56)n);/*輸入序列對應(yīng)值*/for(i=0;isize_x;i+)scanf(%lf%lf,&xi.real,&xi.img);initW();/*選擇FFT逆FFT運(yùn)算*/printf(UseFFT(0)orIFFT(1)?n)

10、;scanf(%d,&method);if(method=0)fft();elseifft();output();return0;/*進(jìn)行基-2FFT運(yùn)算*/voidfft()inti=0,j=0,k=0,l=0;complexup,down,product;change();for(i=0;ilog(size_x)/log(2);i+)/*一級蝶形運(yùn)算*/l=1i;for(j=0;jsize_x;j+=2*l)/*一組蝶形運(yùn)算*/for(k=0;kl;k+)/*一個(gè)蝶形運(yùn)算*/mul(xj+k+l,Wsize_x*k/2/l,&product);add(xj+k,product,&up);s

11、ub(xj+k,product,&down);xj+k=up;xj+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+)/*一級蝶形運(yùn)算*/l/=2;for(j=0;jsize_x;j+=2*l)/*一組蝶形運(yùn)算*/for(k=0;kl;k+)/*一個(gè)蝶形運(yùn)算*/add(xj+k,xj+k+l,&up);up.real/=2;up.img/=2;sub(xj+k,xj+k+l,&down);down.real/=2;down.img/=2;divi(down

12、,Wsize_x*k/2/l,&down);xj+k=up;xj+k+l=down;change();/*初始化變化核*/voidinitW()inti;W=(complex*)malloc(sizeof(complex)*size_x);for(i=0;isize_x;i+)Wi.real=cos(2*PI/size_x*i);Wi.img=-1*sin(2*PI/size_x*i);/*變址計(jì)算，將x(n)碼位倒置*/voidchange()complextemp;unsignedshorti=0,j=0,k=0;doublet;for(i=0;i0)j=j1;if(ji)temp=xi;xi=xj;xj=temp;voidoutput()/*輸出結(jié)果*/inti;printf(Theresultareasfollowsn);for(i=0;i=0.0001)printf(+%.4fjn,xi.img);elseif(fabs(xi.img)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,comple