C語言編寫FFT程序

1、DSP課程作業(yè)用C語言編寫FFT程序1快速傅里葉變換FFT簡介快速傅氏變換（FFT）,是離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的奇、 偶、虛、實等特性，對離散傅立葉變換的算法進(jìn)行改進(jìn)獲得的。它對傅氏變換的理論并 沒有新的發(fā)現(xiàn)，但是對于在計算機(jī)系統(tǒng)或者說數(shù)字系統(tǒng)中應(yīng)用離散傅立葉變換，可以說 是進(jìn)了 一大步。我們假設(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ù)

2、序列的X（ m ,即N點DFT變換大約就 需要NT次運算。當(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+（N/2）2=N+N2/2。 繼續(xù)上面的例子，N=1024時，總的運算次數(shù)就變成了525312次，節(jié)省了大約50%勺運算量。而如果我們將這種“一分為二”的思想不斷進(jìn)行下去，直到分成兩兩一組的DFT運算

3、單元，那么N點的DFT變換就只需要 Nlog2N次的運算，N在1024點時，運算量僅有 10240次，是先前的直接算法的1%點數(shù)越多，運算量的節(jié)約就越大，這就是 FFT的優(yōu)越性。2, FFT算法的基本原理FFT 算法的基本思想：利用DFT系數(shù)的特性，合并DFT運算中的某些項， 吧長序列的 DFT-短序列的DFT,從而減少其運算量。FFT 算法分類：時 間 抽選法 DIT: Decimation-ln-Time;頻率抽選法 DIF:Decimati on-I n-F reque ncy 按時間抽選的基-2FFT算法1、算法原理設(shè)序列點數(shù)N = 2L , L為整數(shù)。若不滿足，則補(bǔ)零。N為2的整數(shù)幕

4、的FFT算法稱基-2FFT算法。將序列 x（n）按n的奇偶分成兩組x2rrr N d0,1,., 1x2r1x2 r2則 x(n)的 DFT:N1Nnk1nkN 1nkX kx nWNx n WNx n WN2 ix1(r)W；kr 02N I2 Ir 0x(2r)叫x(2r)W2(k 0,1,.； 1)Xi (A) -V2(A)NN彳2I2 1x 2r wN2rkx 2r2r 1 k1 Wn丁1°r 0N 122 rkk2 rkX1 rWnWNX2r Wnr 0r 0N 1N2 1X1 r WN2WNkX2r W,/2r 0r 0kX1 k WnX2k(r,k0,1,.N 1)2其

5、中X!(k)Xi(k)再利用周期性求X(k)的后半部分：QXi k汛2 k是以N為周期的X1 k NX1 k X2 k N X2 k2 2k NN又Wn 刁 Wn2W,w,kX(k) Xi(k) WnX2(GNkX(k ) Xi(k) WzX2(k)n為偶數(shù)圖4時聞抽選法蝶形運算流圖符號n為奇數(shù)V(0)點V(<2)DI TXiO)4»v<1)x；(! )=a (3)點DFTa-(3)=a<7)ffl4 2按時間抽選，將仆V點DFT分鵬為 兩個22總DFTVz( 2) <分解后的運算量:妲數(shù)乘法復(fù)數(shù)加法一個N/2點 DFT(V/2)2N/2(N/2-l)兩個N/

6、2點DF1N2/2/V(V/2 1)個蝶形12 '和2個蝶形/V/2N總計/VJ/2 + A/2*Af2/2N(N/2-)-Nf/2運JTSW少了近半!2 ）、運算量當(dāng)N = 2L時，共有L級蝶形，每級N / 2個蝶形, 每比較DFTmF (DFT )N 22 Nm f ( FFT )2 log 2 Nlog 2 N3 )、算法特點原位計算蝶形運算兩節(jié)點的第一個節(jié)點為 把右邊空出的位置補(bǔ)零，結(jié)果為k值，表示成L位二進(jìn)制數(shù)，左移 Lr的二進(jìn)制數(shù)。m位,Xm(k) Xmi(k) Xmi(j)WNXm(j) Xmi(k) Xmi(j)WNA i (/) W §Xi® 一計

7、一,匕(/)=匕戀)_K泊呎你-I圖47按時間抽選蝶形運算結(jié)構(gòu)倒位序蝶形運算對N = 2L點FFT，輸入倒位序，輸出自然 序，第m級運算每個蝶形的兩節(jié)點距離為2m -倒位序自然序00000000100410010102201011063()1100114100101551010113611011177111m 1 rXm(k) Xmi(k) Xmi(k 2 )WnXm(k 2m1) Xmi(k) Xmi(k 2m1)wN.VnUAU因一兀*朋叭V(A)-Vm.|(A)+.；ll0,Yju_?(/) *»ffl 4-7按吋間抽選蝶形運算結(jié)構(gòu)wN的確定蝶形運算兩節(jié)點的第一個節(jié)點為k值，表

8、示成L位二進(jìn)制數(shù)，左移L -m位，把右邊空出的位置補(bǔ)零，結(jié)果為 r的二進(jìn)制數(shù)。存儲單元輸入序列x(n) : N個存儲單元系數(shù)WN : N / 2個存儲單元3，快速傅立葉變換的C語言實現(xiàn)方法我們要衡量一個處理器有沒有足夠的能力來運行FFT算法，根據(jù)以上的簡單介紹可以得出以下兩點:1. 處理器要在一個指令周期能完成乘和累加的工作，因為復(fù)數(shù)運算要多次查表相乘 才能實現(xiàn)。2. 間接尋址，可以實現(xiàn)增/減1個變址量，方便各種查表方法。FFT要對原始序列進(jìn) 行反序排列，處理器要有反序間接尋址的能力。#i nclude <stdio.h>#in clude <math.h>#in cl

9、ude <stdlib.h>#defi ne N 1000typedef structdouble real; double img; complex;voidfft();voidifft();voidinitW(); /*初始化變化核 */voidchange();voidadd(complex,complex,complexvoidmul(complex,complex,complexvoidsub(complex,complex,complexvoiddivi(complex,complex,complexvoidoutput();*);*);*);*);complex xN

10、, *W;/* 輸出序列的值 */int size_x=O;/* 輸入序列的長度，只限 2的N次方*/ double PI;int main()int i,method;system("cls");PI=atan(1)*4;/*pi等于4乘以 1.0 的正切值 */printf("Please input the size of x:n");/* 輸入序列的長度 */scanf("%d",&size_x);printf("Please input the data in xN:(such as:5 6)n"

11、); /* 輸入序列對應(yīng)的值 */for(i=0;i<size_x;i+)scanf("%lf %lf",&xi.real,&xi.img);initW();/*選擇FFT或逆FFT運算*/printf("Use FFT(0) or IFFT(1)?n"); scanf("%d",&method);if(method=0)fft();elseifft(); output();return 0;/* 進(jìn)行基 -2 FFT 運算 */void fft()int i=0,j=0,k=0,l=0;complex

12、up,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(xj+k+l,Wsize_x*k/2/l,&product);add(xj+k,product,&up);sub(xj+k,product,&down);xj+k=up;xj+k+l=down;void ifft()int i=0,j=0,k

13、=0,l=size_x;一級蝶形運算 */complex up,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(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,Wsize_x*k/2/l,&down);xj

14、+k=up;xj+k+l=down;change();/* 初始化變化核 */void initW()int i;size_x);W=(complex *)malloc(sizeof(complex)for(i=0;i<size_x;i+)Wi.real=cos(2*PI/size_x*i);Wi.img=-1*sin(2*PI/size_x*i);/* 變址計算，將 x(n) 碼位倒置 */void change()complex temp;unsigned short i=0,j=0,k=0;double t;for(i=0;i<size_x;i+)k=i;j=0;t=(log

15、(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;printf("The result are as followsn"); for(i=0;i<size_x;i+)printf("%.4f",xi.real);if(xi.img>=0.0001) printf("+%.4fjn",xi.img);

16、else if(fabs(xi.img)<0.0001) printf("n");elseprintf("%.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(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(com