用c語(yǔ)言實(shí)現(xiàn)的FFT

1、一、對(duì) FFT 的介紹1. FFT（ Fast Fourier Transformation）,即為快速傅里葉變換，是離散傅里葉變換的快速算法，它是根據(jù)離散傅里葉變換的奇、偶、虛、實(shí)等特性，對(duì)離散傅里葉變換的算法進(jìn)行改進(jìn)獲得的。2.FFT 算法的基本原理FFT 算法是把長(zhǎng)序列的DFT 逐次分解為較短序列的DFT 。按照抽取方式的不同可分為 DIT-FFT( 按時(shí)間抽取 )和 DIF-FFT( 按頻率抽取 )算法。按蝶形運(yùn)算的構(gòu)成不同可分為基 2，基 4，基 8，以及任意因子的類(lèi)型。3.迭代關(guān)系4、本次程序的基本過(guò)程.我們這次所研究的是數(shù)字信號(hào)處理中的FFT 算法，我們這次所用的數(shù)字信號(hào)是復(fù)數(shù)類(lèi)

2、型的。(1) 所以首先，我們先定義了一個(gè)復(fù)數(shù)結(jié)構(gòu)體，因?yàn)槭沁M(jìn)行復(fù)數(shù)的運(yùn)算，我們又相繼定義復(fù)數(shù)的加減乘運(yùn)算的函數(shù)。(2) 緊接著，我們定義了進(jìn)行 FFT 計(jì)算的 fft() 快速傅里葉變換函數(shù)initW()初始化變換核函數(shù)即旋轉(zhuǎn)因子的計(jì)算，change()變址函數(shù)，output()輸出傅里葉變換的結(jié)果的函數(shù)。(3) 定義主函數(shù)， 并調(diào)用定義好的相關(guān)子函數(shù)， 利用 fft() 中的蝶形運(yùn)算以及 change() 函數(shù)來(lái)完成從時(shí)間域上選取的 DIT-FFT 。二、 FFT 中碼位倒置排序1、碼位倒置的實(shí)現(xiàn)方法：(1) 簡(jiǎn)單的利用按位與、或循環(huán)實(shí)現(xiàn)(2) 利用公式推導(dǎo)的迭代方法2、為什么要進(jìn)行碼位倒置

3、因?yàn)橛捎?FFT 的計(jì)算特性，如果按照正常順序輸入，而沒(méi)有進(jìn)行碼位倒置的話，就會(huì)以亂序輸出， 就不便于我們后續(xù)對(duì)信號(hào)的相關(guān)性質(zhì)進(jìn)行研究了，所以 DIT-FFT 算法就是在進(jìn)行 FFT 計(jì)算之前，進(jìn)行分奇偶后的碼位倒置運(yùn)算，即二進(jìn)制數(shù)的倒位。3、倒位序由奇偶分組造成，以N=8 為例，說(shuō)明如下：.三、蝶形運(yùn)算由.按照上述公式的規(guī)律進(jìn)行逐級(jí)分解，直到2 點(diǎn) DFT ，如下是 N=8 時(shí)的蝶形算法分析圖：四、 FFT 算法中蝶形算法的基本思想分析（1）我們知道 N 點(diǎn) FFT 運(yùn)算可以分成 log2 （N ）級(jí)，每一級(jí)都有 N/2 個(gè)碟形， FFT 的基本思想是用 3 層循環(huán)完成全部運(yùn)算 (N 點(diǎn) F

4、FT)。.（ 2）第一層循環(huán)：由于 N=2m 需要 m 級(jí)計(jì)算 ,第一層循環(huán)對(duì)運(yùn)算的級(jí)數(shù)進(jìn)行控制。（ stages ）（ 3）第二層循環(huán)：由于第 L 級(jí)有 2(L-1) 個(gè)蝶形因子（乘數(shù)），第二層循環(huán)根據(jù)乘數(shù)進(jìn)行控制， 保證對(duì)于每一個(gè)蝶形因子第三層循環(huán)要執(zhí)行一次，這樣，第三層循環(huán)在第二層循環(huán)控制下，每一級(jí)要進(jìn)行2(L-1)次循環(huán)計(jì)算 .(選擇 W)（4）第三層循環(huán)：由于第L 級(jí)共有 N/2L即 2 （ n-L ）個(gè)群，并且同一級(jí)內(nèi)不同群的乘數(shù)分布相同，當(dāng)?shù)诙友h(huán)確定某一乘數(shù)后，第三層循環(huán)要將本級(jí)中每個(gè)群中具有這一乘數(shù)的蝶形計(jì)算一次， 即第三層循環(huán)每執(zhí)行完一次要進(jìn)行 N/2L 個(gè)碟形計(jì)算。 （

5、執(zhí)行不同 group 中具有相同 W 的蝶形運(yùn)算）（ 5）可以得出結(jié)論：在每一級(jí)中，第三層循環(huán)完成 N/2L 個(gè)碟形計(jì)算；第二層循環(huán)使第三層循環(huán)進(jìn)行 2(L-1) 次，因此，第二層循環(huán)完成時(shí)，共進(jìn)行 2(L-1) *N/2L=N/2個(gè)碟形計(jì)算。實(shí)質(zhì)是：第二、第三層循.環(huán)完成了第 L 級(jí)的計(jì)算。五、用 c 語(yǔ)言實(shí)現(xiàn)的 FFT 算法如下：<span style="font-size:18px;">#include <stdio.h>#include <math.h>#include <stdlib.h>#define N 1000

6、/* 定義復(fù)數(shù)類(lèi)型 */typedef structdouble real;double img;complex;complex xN, *W; /*輸入序列 , 變換核 */int size_x=0;/*輸入序列的大小，在本程序中僅限2 的次冪 */double PI;/*圓周率 */void fft(); /*快速傅里葉變換 */void initW();/*初始化變換核 */void change(); /*變址 */void add(complex ,complex ,complex *); /*復(fù)數(shù)加法 */void mul(complex ,complex ,complex *);

7、 /*復(fù)數(shù)乘法 */void sub(complex ,complex ,complex *); /*復(fù)數(shù)減法 */void output();/*輸出快速傅里葉變換的結(jié)果*/int main()int i;/*輸出結(jié)果 */system("cls");PI=atan(1)*4;printf("輸出 DIT 方法實(shí)現(xiàn)的FFT 結(jié)果n");printf("Please input the size of x:n");/輸入序列的大小scanf("%d",&size_x);printf("Please

8、 input the data in xN:n");/輸入序列的實(shí)部和虛部for(i=0;i<size_x;i+)printf("請(qǐng)輸入第 %d 個(gè)序列： ",i);scanf("%lf%lf",&xi.real,&xi.img);printf("輸出倒序后的序列n");initW();/調(diào)用變換核fft();/調(diào)用快速傅里葉變換printf("輸出 FFT 后的結(jié)果 n");.output();/調(diào)用輸出傅里葉變換結(jié)果函數(shù)return 0;/* 快速傅里葉變換*/void fft(

9、)int i=0,j=0,k=0,l=0;complex up,down,product;change(); /調(diào)用變址函數(shù)for(i=0;i< log(size_x)/log(2) ;i+)/*一級(jí)蝶形運(yùn)算stage */l=1<<i;for(j=0;j<size_x;j+= 2*l )/*一組蝶形運(yùn)算group,每個(gè) group 的蝶形因子乘數(shù)不同*/for(k=0;k<l;k+)/*一個(gè)蝶形運(yùn)算每個(gè) group內(nèi)的蝶形運(yùn)算 */mul(xj+k+l,Wsize_x*k/2/l,&product);add(xj+k,product,&up);s

10、ub(xj+k,product,&down);xj+k=up;xj+k+l=down;/* 初始化變換核，定義一個(gè)變換核，相當(dāng)于旋轉(zhuǎn)因子WAP*/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);/用歐拉公式計(jì)算旋轉(zhuǎn)因子Wi.img=-1*sin(2*PI/size_x*i);/* 變址計(jì)算，將x(n)碼位倒置 */void change()complex temp;unsigned short

11、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 )/利用按位與以及循環(huán)實(shí)現(xiàn)碼位顛倒j=j<<1;j|=(k & 1);k=k>>1;if(j>i)/將 x(n)的碼位互換temp=xi;xi=xj;xj=temp;output();/* 輸出傅里葉變換的結(jié)果*/void output()int i;printf("The result are as follows： n");for(i=0;i<s

12、ize_x;i+)printf("%.4f",xi.real);if(xi.img>=0.0001)printf("+%.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) /復(fù)數(shù)加法的定義c->real=a.real+b.real;c->img=a.img+b.img;void mul(compl

13、ex a,complex b,complex *c) /復(fù)數(shù)乘法的定義void sub(complex a,complex b,complex *c) /復(fù)數(shù)減法的定義.c->real=a.real-b.real;c->img=a.img-b.img;</span>六、 FFT 原理的理解和程序設(shè)計(jì)中遇到的相關(guān)問(wèn)題及解決方法1、遇到的相關(guān)問(wèn)題：(1) 首先一開(kāi)始不知道什么是 FFT，以及 FFT 原理是什么(2) 不理解 FFT 中迭代關(guān)系的推導(dǎo)以及緣由(3) 不理解變址計(jì)算的原理(4) 對(duì)蝶形運(yùn)算的推導(dǎo)與原理的理解不透徹(5) 編程過(guò)程中對(duì)變址計(jì)算即對(duì)按位與的變換形式的不理解2、解決的方法：(1) 到圖書(shū)館借相關(guān)的書(shū)籍理解相關(guān)的原理和過(guò)程(2) 到百度收索相關(guān)的資料促進(jìn)理解相應(yīng)的原理過(guò)程(3) 向