fft算法代碼注釋及流程框圖

1、精選優(yōu)質(zhì)文檔-傾情為你奉上基2的DIT蝶形算法源代碼及注釋如下： /*FFT*/ /整個程序輸入和輸出利用同一個空間xN，節(jié)約空間 #include <stdio.h> #include <math.h> #include <stdlib.h> #define N 1000 /定義輸入或者輸出空間的最大長度 typedef struct double real; double img; complex; /定義復(fù)數(shù)型變量的結(jié)構(gòu)體 void fft(); /快速傅里葉變換函數(shù)聲明 void initW(); /計算W(0)W(size_x-1)的值函數(shù)聲明 v

2、oid change(); /碼元位置倒置函數(shù)函數(shù)聲明 void add(complex,complex,complex *); /*復(fù)數(shù)加法*/ void mul(complex,complex,complex *); /*復(fù)數(shù)乘法*/ void sub(complex,complex,complex *); /*復(fù)數(shù)減法*/ void divi(complex,complex,complex *); /*復(fù)數(shù)除法*/ void output(); /*輸出結(jié)果*/ complex xN,*W; /*輸出序列的值*/ int size_x=0; /*輸入序列的長度，只限2的N次方*/ dou

3、ble PI; /pi的值 int main() int i; system("cls"); PI=atan(1)*4; printf("Please input the size of x:n"); /*輸入序列的長度*/ scanf("%d",&size_x); printf("Please input the data in xN:(such as:5 6)n"); /*輸入序列對應(yīng)的值*/ for(i=0;i<size_x;i+) scanf("%lf %lf",&

4、xi.real,&xi.img); initW(); /計算W(0)W(size_x-1)的值 fft(); /利用fft快速算法進行DFT變化 output(); /順序輸出size_x個fft的結(jié)果 return 0; /*進行基-2 FFT運算,蝶形算法。這個算法的思路就是， 先把計算過程分為log(size_x)/log(2)-1級（用i控制級數(shù)）； 然后把每一級蝶形單元分組（用j控制組的第一個元素起始下標）； 最后算出某一級某一組每一個蝶形單元（用k控制個數(shù)，共l個）。*/ void fft() int i=0,j=0,k=0,l=0; complex up,down,pro

5、duct; change(); /實現(xiàn)對碼位的倒置 for(i=0;i<log(size_x)/log(2);i+) /循環(huán)算出fft的結(jié)果 l=1<<i; for(j=0;j<size_x;j+=2*l) /算出第m=i級的結(jié)果【i從0到（log(size_x)/log(2)）-1】 for(k=0;k<l;k+) /算出第i級內(nèi)j組蝶形單元的結(jié)果 /算出j組中第k個蝶形單元 mul(xj+k+l,W(size_x/2/l)*k,&product); /*size/2/l是該級W的相鄰上標差，l是該級該組取的W總個數(shù)*/ add(xj+k,product

6、,&up); sub(xj+k,product,&down); xj+k=up; /up為蝶形單元右上方的值 xj+k+l=down; /down為蝶形單元右下方的值 void initW() /計算W的實現(xiàn)函數(shù) int i; W=(complex *)malloc(sizeof(complex) * size_x); /*申請size_x個復(fù)數(shù)W的空間(這部申請的空間有點多，實際上只要申請size_x/2個即可)*/ for(i=0;i<(size_x/2);i+) /*預(yù)先計算出size_x/2個W的值，存放,由于蝶形算法只需要前size_x/2個值即可*/ Wi.r

7、eal=cos(2*PI/size_x*i); /計算W的實部 Wi.img=-1*sin(2*PI/size_x*i); /計算W的虛部 void change() /輸入的碼組碼位倒置實現(xiàn)函數(shù) 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(size_x)/log(2); while(t-)>0) j=j<<1; j|=(k & 1); k=k>>1; if(j>i) temp=xi; xi=xj; xj=te

8、mp; void output() /輸出結(jié)果實現(xiàn)函數(shù) int i; printf("The result are as followsn"); for(i=0;i<size_x;i+) printf("%.4f",xi.real); /輸出實部 if(xi.img>=0.0001) /如果虛部的值大于0.0001，輸出+jx.img的形式 printf("+j%.4fn",xi.img); else if(fabs(xi.img)<0.0001) /如果虛部的絕對值小于0.0001，無需輸出 printf(&quo

9、t;n"); else printf("-j%.4fn",fabs(xi.img); /如果虛部的值小于-0.0001，輸出-jx.img的形式 void add(complex a,complex b,complex *c) /復(fù)數(shù)加法實現(xiàn)函數(shù) c->real = a.real + b.real; /復(fù)數(shù)實部相加 c->img = a.img + b.img; /復(fù)數(shù)虛部相加 void mul(complex a,complex b,complex *c) /復(fù)數(shù)乘法實現(xiàn)函數(shù) c->real = a.real*b.real - a.img*b.

10、img; /獲取相乘結(jié)果的實部 c->img = a.real*b.img + a.img*b.real; /獲取相乘結(jié)果的虛部 void sub(complex a,complex b,complex *c) /復(fù)數(shù)減法實現(xiàn)函數(shù) c->real = a.real - b.real; /復(fù)數(shù)實部相減 c->img = a.img - b.img; /復(fù)數(shù)虛部相減 void divi(complex a,complex b,complex *c) /復(fù)數(shù)除法實現(xiàn)函數(shù) c->real = (a.real*b.real + a.img*b.img) / (b.real*b.real+b.img*b.img); /獲取相除結(jié)果的實部 c->img = (a.img*b.real - a.real*b.img) / (b.real*b.real+b.img*b.img); /獲取相除結(jié)果的虛部 由于輸入序列和輸出序列存儲（包括中間結(jié)果）在同一個空間數(shù)組x【size_x】，故每進行一個蝶形算法，相應(yīng)的存儲單元有兩個值被替換。該程序框圖如下：開始計算出前size_x/2個exp(-j*2*k/size_x)個值，即W的值輸入序列對應(yīng)值（例如5+j3，輸入5 3）輸入序