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（

2、m）,即N點(diǎn)DFT變換大約就需要N2次運(yùn)算。當(dāng)N=1024點(diǎn)甚至更多的時候，需要N2=1048576次運(yùn)算，在FFT中，利用WN的周期性和對稱性，把一個N項(xiàng)序列（設(shè)N=2k,k為正整數(shù)），分為兩個N/2項(xiàng)的子序列，每個N/2點(diǎn)DFT變換需要（N/2）2次運(yùn)算，再用N次運(yùn)算把兩個N/2點(diǎn)的DFT變換組合成一個N點(diǎn)的DFT變換。這樣變換以后，總的運(yùn)算次數(shù)就變成N+2（N/2）2=N+N2/2。繼續(xù)上面的例子，N=1024時，總的運(yùn)算次數(shù)就變成了525312次，節(jié)省了大約50%的運(yùn)算量。而如果我們將這種“一分為二”的思想不斷進(jìn)行下去，直到分成兩兩一組的DFT運(yùn)算單元，那么N點(diǎn)的DFT變換就只需要Nl

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

4、分成兩組：則x(n)的DFT: 其中 再利用周期性求X(k)的后半部分： 復(fù)數(shù)乘法： 當(dāng)N = 2L時，共有L級蝶形，每級N / 2個蝶形，每個蝶形有1次復(fù)數(shù)乘法2次復(fù)數(shù)加法。 2）、運(yùn)算量 復(fù)數(shù)乘法：2、運(yùn)算量 復(fù)數(shù)加法：比較DFT 3）、算法特點(diǎn)原位計(jì)算 倒位序 蝶形運(yùn)算 對N = 2L點(diǎn)FFT，輸入倒位序，輸出自然序，第m級運(yùn)算每個蝶形的兩節(jié)點(diǎn)距離為 2m1第m級運(yùn)算： 蝶形運(yùn)算兩節(jié)點(diǎn)的第一個節(jié)點(diǎn)為k值，表示成L位二進(jìn)制數(shù)，左移L m位，把右邊空出的位置補(bǔ)零，結(jié)果為r的二進(jìn)制數(shù)。 蝶形運(yùn)算兩節(jié)點(diǎn)的第一個節(jié)點(diǎn)為k值，表示成L位二進(jìn)制數(shù)，左移L m位，把右邊空出的位置補(bǔ)零，結(jié)果為r的二進(jìn)制數(shù)

5、。 存儲單元 輸入序列x(n) : N個存儲單元 系數(shù) ：N / 2個存儲單元 3，快速傅立葉變換的C語言實(shí)現(xiàn)方法我們要衡量一個處理器有沒有足夠的能力來運(yùn)行FFT算法，根據(jù)以上的簡單介紹可以得出以下兩點(diǎn)：1. 處理器要在一個指令周期能完成乘和累加的工作，因?yàn)閺?fù)數(shù)運(yùn)算要多次查表相乘才能實(shí)現(xiàn)。 2. 間接尋址，可以實(shí)現(xiàn)增/減1個變址量，方便各種查表方法。FFT要對原始序列進(jìn)行反序排列，處理器要有反序間接尋址的能力。 #include <stdio.h> #include <math.h> #include <stdlib.h> #define N 1000 ty

6、pedef struct double real; double img; complex; void fft(); void ifft(); void initW(); /*初始化變化核*/ void change(); void add(complex ,complex ,complex *); void mul(complex ,complex ,complex *); void sub(complex ,complex ,complex *); void divi(complex ,complex ,complex *); void output(); complex xN, *W;/

7、*輸出序列的值*/ int size_x=0;/*輸入序列的長度，只限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"); /*輸

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

9、x up,down,product; change(); for(i=0;i< log(size_x)/log(2) ;i+) /*一級蝶形運(yùn)算*/ l=1<<i; for(j=0;j<size_x;j+= 2*l ) /*一組蝶形運(yùn)算*/ for(k=0;k<l;k+) /*一個蝶形運(yùn)算*/ 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

10、,j=0,k=0,l=size_x; complex up,down; for(i=0;i< (int)( log(size_x)/log(2) );i+) /*一級蝶形運(yùn)算*/ l/=2; for(j=0;j<size_x;j+= 2*l ) /*一組蝶形運(yùn)算*/ for(k=0;k<l;k+) /*一個蝶形運(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,Wsize_x*k/2/l,&

11、;down); xj+k=up; xj+k+l=down; change(); /*初始化變化核*/ void initW() 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*sin(2*PI/size_x*i); /*變址計(jì)算，將x(n)碼位倒置*/ void change() complex temp; unsigned short i=0,j=0,k=0; double t; for(i=0;i<si

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

13、intf("+%.4fjn",xi.img); else if(fabs(xi.img)<0.0001) printf("n"); else printf("%.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(complex a,complex b,complex *c) c->real=a.real-b.