C語言編寫FFT程序

1、sinclude 短序列的DFT,從而減少其運算量。FFT算法分類：時間抽選法DIT: Decimation-In-Time；頻率抽選法DIF: Dec imat ionTn - Frequency按時間抽選的基-2FFT算法1、算法原理設序列點數(shù)N二2L, L為整數(shù)。若不滿足，貝IJ補零。N為2的整數(shù)幕的FFT算法稱基-2FFT算法。將序列x(n) 按n的奇偶分成兩t(2r) = x1(r).N/ =U, !,，1x(2r + l) = x2(r)2則x(n)的)=餐(”)哼=%()“)(w汙+叫12(廠)(w汗r=()r=0二叱2+說藝2(廠)吒2?=()r=()= X) +睨 X2(k)

2、(以=0,.善1)其中 41卜1X“)=$(2門嘴r-0r-()卜1卜1x& ) = 5)唱=1(2)呼伙=0 丄上-1)r=0*r=02再利用周期性求X(k)的后半部分： X伙),蜀伙)是以孕為周期的X+分 X) X2”+號卜 X2又比戶二昨吆=-觀X(k) = X) + W；X2(R) N.X(k + ) = Xi(k)-WX2(k)w乙X2(A)X1(*) + wv xeX心)-環(huán) X2(k)圖41時間抽選法蝶形運算流圖符號畀為偶數(shù)打為奇數(shù)2sinclude std:oxi(O)-x(O) A71(1 )=(2) xi(2)-v(4) xi(3)=x(6) X2(0)-X(l) 心1)=

3、5V X2(2)=x(5) .v2(3)=v(7)乎點DFTDFTkVi(O)r x（o）Yi(l)Vi (2)M2)PVi(3)X(3)X(l)A(4)W)0(1)M2)Xz(3)-圖42按時間抽選，將一個”點DFT分解為兩個必點DFT|分解后的運算量:復數(shù)乘法復數(shù)加法一個 N/2 點 DFT(N/2)2N/2(N/2_1)兩個N/2點DFT2/2N(N/2-1)一個蝶形12N/2個蝶形N/2N總計N2/2+N/2 N2HN(N/2-l) + NN2!2運算量滅少了近一半!2運算量當N二2L時，共有L級蝶形，每級N/ 2個蝶形，每NN復數(shù)乘加二厶二log。N 法：22復數(shù)加 u 二 NL =

4、 NojN 法比較DFTmF（DFT） N?2Nm/FFT匕隅屮=3).算法特原位計算X 腫= XM)+ XM)W(Xm(j) = X 心伙)xX心)=X,”.“)+Xi(/)0 s%心)一也昭圖4.7按時間抽選蝶形運算結構倒位x(no 件(呼/J0 0 04蝶形運sinclude stdo倒位序000100仗）+心 i （j）爐左移L - m010110Yjm./0 圖47按時間抽選蝶形運算結構%（/）*】（/）%001V；的確定4100心伙）自然眾）+ XM + 2MX,”伙+ 2心）=X心-Xg伙+ 2心）叫7101兩彳活的第為k值，表示成L位二進制數(shù)，左移LOIL m位，魏右邊空：6的

5、位置減製）結果冷r的二進制數(shù)。111 存旳單 7111對N二2L點FFT,輸入倒位序，輸出自然序,2m - 1系數(shù) ：N/2個存儲單元3, 快速傅立葉變換的C語言實現(xiàn)方法我們要衡量一個處理器有沒有足夠的能力來運行FFT算法，根據(jù)以上的簡單介紹可 以得出以下兩點：1. 處理器要在一個指令周期能完成乘和累加的工作，因為復數(shù)運算要多 次査表相乘才能實現(xiàn)。2. 間接尋址，可以實現(xiàn)增/減1個變址量，方便各種査表方法。FFT要對 原始序列進行反序排列，處理器要有反序間接尋址的能力。#include#include#includeSdefine N 1000typedef structidoublereal

6、;doubleimg;complex;voidvoidifftO;voidinitWO; /*初始化變化核*/voidchange 0;voidadd(complext complex,complex*);voidmul(complex,complex,complex*);voidsub(complex,complex,complex*);voiddivi(complex,complex,complex*);voidoutput 0;complex x N,*W; /*輸出序列的值*/intsize_x=0;/*輸入序列的長度，只限2的N次方*/double PI;int mainOsincl

7、ude stdioint i,method;system(cls);PI=atan(1)*4;/*pi等于4乘以10的正切值*/printf (/Please input the size of x:rT);/*輸入序列的長度*/scanf&size_x);printf (Please input the data in xN: (such as:5 6)n);/*輸入序列對應的值*/for(i=0;isize_x;i+)scanf (“%If %lfz &xi real, &xi img);initW();/*選擇FFT或逆FFT運算*/printf CUse FFT(O) or IFFT(

8、l)?n);scanf f fonethod);if (method=0)fft 0 ；elseifftO;output 0;return 0;/*進行基-2 FFT運算*/void fft0int i=0, j=0, k=0, 1=0;complex up, down, product;change 0;for (i=0; i log (size_x)/log (2); i+) /*一級蝶形運算*/l=li;for (j二0; jsize_x; j+=2*1) /*一組蝶形運算*/for(k=0;kl;k+) /*一個蝶形運算*/mul (x L j +k+1 , W1 size_x*k4/

9、2/1 , product);add(xLj+k, product, &up);sub (xlj+k, product, &down); xj+k二up;xj+k+l=down;void ifft ()inti=0, j二0, k=0, 1 二size_x;complex up, down;for(i=0; i (int) ( log(size_x)/log(2) ) ; i+) /*一級蝶形運算*/ l/=2;for (j二0; jsize_x; j+=2*1) /*一組蝶形運算*/for(k=0;kl;k+) /*一個蝶形運算*/add(xj+k, xj+k+l, &up);up.real

10、/=2;up img/=2;sub (xLj+k, xj+k+l, &down);down real/=2;down img/=2;divi(down, Wsize_x*k/2/l, &down);xj+k二up;xj+k+1=down;change 0;/*初始化變化核*/void initWOint i;W=(complex*)malloc(sizeof(complex)* size_x);for(i=0;isize_x;i+)Wi real=cos(2*PI/size_x*i);Wi. img=-l*sin(2*PI/size_x*i); /*變址計算，將x(n)碼位倒置*/void c

11、hange 0complex temp;unsigned short i=0,j=0,k=0;double t;for(i=0;i0 )J=J1；j =(k&1);k=kl;if(ji)temp二xi;xi=xj;xj二temp;void output 0/*輸岀結果*/int i;followsn);printf (Z/The result are asfor(i=0;i=0. 0001)printf (+% 4f jn, xi img);else if(fabs(xij. img)real=a real+b real;c-img=a img+b img;void mul(complex 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 afcomplex b, complex *c)c-real=a real-b real;c-img=a img-