算法設(shè)計(jì)與實(shí)現(xiàn)

1、FFT算法研究報(bào)告 1、 程序設(shè)計(jì)背景(FFT算法理解)FFT(fast fourier transformation),快速傅里葉變換。是對DFT算法的改進(jìn)，其利用了WNnk的周期性、共軛對稱性和可約性，使得DFT中有些項(xiàng)可以合并，大大減小了計(jì)算量。按輸入序列在時(shí)間上的次序是屬于偶數(shù)還是奇數(shù)來分解稱為“按時(shí)間抽取法”(DIT)。另一種是把輸出序列X(k)按順序的奇偶分解為越來越短的序列，稱為按頻率抽樣的FFT算法（DIF）。DIT算法是先作復(fù)乘后作加減，而DIF的復(fù)乘只出現(xiàn)在減法之后。本次程序采用DIT算法實(shí)現(xiàn)FFT。用c語言實(shí)現(xiàn)FFT的難點(diǎn)在于數(shù)據(jù)倒位序的處理，以及各級蝶形運(yùn)算的實(shí)現(xiàn)。倒位

2、序的實(shí)現(xiàn)可以使用“反向進(jìn)位加法”，即倒位序二進(jìn)制數(shù)的下面一個數(shù)是上面一個數(shù)在最高位加一并由高位向低位進(jìn)位而得到的。對于點(diǎn)數(shù)為N = 2L的FFT運(yùn)算，可以分解為L階蝶形圖級聯(lián)，第M階蝶形圖內(nèi)又分為2(L-M)個蝶形組，每個蝶形組內(nèi)包含2(M-1)個蝶形。而且旋轉(zhuǎn)因子與蝶形階數(shù)和蝶形分組內(nèi)的蝶形個數(shù)存在關(guān)聯(lián)。因此我們就可以構(gòu)造循環(huán)來實(shí)現(xiàn)蝶形運(yùn)算。2、 FFT算法流程圖輸入數(shù)組指針a，長度n倒位序流程圖：n=0? Yi=0Ni<n?N i%2=0?Y NTempi=ai/2temp(n+i)/2=ai Yi=i+1a=tempn=n/2結(jié)束Y3、 FFT程序代碼#include <st

3、dio.h> #include <stdlib.h> #include <math.h> #include <time.h>#define PI 3.141592653 #define SIZE 16#define MAX 10 /定義復(fù)數(shù)結(jié)構(gòu)體 typedef struct Complex double real; double imag; comp; /定義復(fù)數(shù)乘法，加法，減法 void Add_(comp * a1,comp *a2,comp *b) b->imag=a1->imag+a2->imag; b->real=a

4、1->real+a2->real; void Sub_(comp * a1,comp *a2,comp *b) b->imag=a1->imag-a2->imag; b->real=a1->real-a2->real; void Mul_(comp * a1,comp *a2,comp *b) double r1=0.0,r2=0.0; double i1=0.0,i2=0.0; r1=a1->real; r2=a2->real; i1=a1->imag; i2=a2->imag; b->imag=r1*i2+r2*

5、i1; b->real=r1*r2-i1*i2; /利用srand()、rand()隨機(jī)生成一個輸入并顯示數(shù)據(jù) void Input(double *data,int n) printf("輸入信號為:n"); srand(int)time(0); for(int i=0;i<SIZE;i+) datai=rand()%MAX; printf("%lfn",datai); /定義WN的n次方項(xiàng)void WN(double n,double size_n,comp * b) double x=2.0*PI*n/size_n; b->ima

6、g=-sin(x); b->real=cos(x); /處理FFT的數(shù)據(jù)位置，輸入倒位序，輸出自然順序（正序），用遞歸結(jié)構(gòu)實(shí)現(xiàn) int reverse(double * a,int n) if(n=1) return 0; double * temp=(double *)malloc(sizeof(double)*n); for(int i=0;i<n;i+) if(i%2=0) tempi/2=ai; else temp(n+i)/2=ai; for(int i=0;i<n;i+) ai=tempi; free(temp); reverse(a, n/2); reverse

7、(a+n/2, n/2); return 1; /定義FFT函數(shù) void FFT(double * a,comp * b,int n) reverse(a, n); /對輸入信號進(jìn)行倒位序處理 int k=n; int m=0; while (k/=2) /記錄需要蝶形運(yùn)算的級數(shù) m+; k=m; comp * a_comp=(comp*)malloc(sizeof(comp)*n); for(int i=0;i<n;i+) a_compi.real=ai; a_compi.imag=0; /將輸入數(shù)據(jù)賦值給 a_comp for(int i=0;i<k;i+) /依次計(jì)算各級蝶

8、形運(yùn)算 z=0; for(int j=0;j<n;j+) if(j/(2(i-1)%2=1) comp wn; WN(z, n, &wn); Mul_(&a_compj, &wn,&a_compj); z+=2(k-i-2); comp temp; int company=j-(2(i-1); temp.real=a_compcompany.real; temp.imag=a_compcompany.imag; Add_(&temp, &a_compj, &a_compcompany ); Sub_(&temp, &

9、a_compj, &a_compj); else m=0; printf("nnFFT處理結(jié)果:n"); for(int i=0;i<n;i+) if(a_compi.imag>=0.0) printf("%lf+%lfjn",a_compi.real,a_compi.imag); else printf("%lf%lfjn",a_compi.real,a_compi.imag); for(int i=0;i<n;i+) bi.imag=a_compi.imag; bi.real=a_compi.real;

10、int main() double arraySIZE; comp bSIZE; Input(array,SIZE); /隨機(jī)生成16個輸入數(shù)據(jù) FFT(array, b, SIZE); 4、 與MATLAB自帶FFT函數(shù)運(yùn)行結(jié)果的比較自編FFT結(jié)果：將輸入信號在matlab中進(jìn)行FFT運(yùn)算，結(jié)果如下：由以上截圖可知，函數(shù)成功地完成了對于序列xn的FFT計(jì)算，但是可以看出matlab的準(zhǔn)確度更高，應(yīng)該是matlab計(jì)算使用的pi值更為精確。為比較二者的時(shí)間，將數(shù)據(jù)長度重新定義為215Matlab FFT程序： 自編FFT程序： 由于計(jì)算matlab用時(shí)使用的數(shù)據(jù)與自編FFT程序所用數(shù)據(jù)不同，可能會有不可預(yù)料的誤差