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

1、一、對FFT的介紹1. FFT （ Fast Fourier Transformation ），即為快速傅里葉變換，是離散 傅里葉變換的快速算法，它是根據(jù)離散傅里葉變換的奇、偶、虛、實(shí)等 特性，對離散傅里葉變換的算法進(jìn)行改進(jìn)獲得的。2. FFT算法的基本原理FFT算法是把長序列的DFT逐次分解為較短序列的 DFT。按照抽取方式的不同可分為 DIT-FFT（按時間抽取）和DIF-FFT（按 頻率抽?。┧惴?。按蝶形運(yùn)算的構(gòu)成不同可分為基 2，基4，基8，以及 任意因子的類型。3. 迭代關(guān)系冷 -1e2 r ? ! NfjJ =0y / 2 -1N J 2 - 1二 Ze2刃比(2 J)/Nfl j

2、 +y長此（2十：八0J = o.V 2 - 1N2-1二 2e(2 j ) yf 2 j +W KY總j = Q/ = 0二 F /+W疋 jza4、本次程序的基本過程我們這次所研究的是數(shù)字信號處理中的FFT算法，我們這次所用的數(shù)字信號是復(fù)數(shù)類型的。(1) 所以首先，我們先定義了一個復(fù)數(shù)結(jié)構(gòu)體，因?yàn)槭沁M(jìn)行復(fù)數(shù)的運(yùn)算，我們又相繼定義復(fù)數(shù)的加減乘運(yùn)算的函數(shù)。(2) 緊接著，我們定義了進(jìn)行 FFT計算的fft()快速傅里葉變換函數(shù)ini tW()初始化變換核函數(shù)即旋轉(zhuǎn)因子的計算，cha nge()變址函數(shù)，output()輸出傅里葉變換的結(jié)果的函數(shù)。(3) 定義主函數(shù)，并調(diào)用定義好的相關(guān)子函數(shù)，利

3、用fft()中的蝶形運(yùn)算以及change()函數(shù)來完成從時間域上選取的DIT-FFT。二、FFT中碼位倒置排序1、碼位倒置的實(shí)現(xiàn)方法：(1) 簡單的利用按位與、或循環(huán)實(shí)現(xiàn)(2) 利用公式推導(dǎo)的迭代方法2、為什么要進(jìn)行碼位倒置因?yàn)橛捎贔FT的計算特性，如果按照正常順序輸入，而沒有進(jìn)行碼位倒置的話，就會以亂序輸出，就不便于我們后續(xù)對信號的相關(guān)性質(zhì)進(jìn)行研究了，所以DIT-FFT算法就是在進(jìn)行FFT計算之前，進(jìn)行分奇偶后的碼位倒置運(yùn)算，即二進(jìn)制數(shù)的倒位。3、倒位序由奇偶分組造成，以 N=8為例，說明如下：Op心樂15卩三、蝶形運(yùn)算由心（力）：心（疋除示x （k助運(yùn)算是一種特殊運(yùn)算1x (疋)=兀 1(

4、£ )+ W x2(k )前一半x (Ar ) = Xj (k )- fF x2(k )后一半(k =(k = 1按照上述公式的規(guī)律進(jìn)行逐級分解，直到2點(diǎn)DFT，如下是N=8時的蝶形算法分析圖：四、FFT算法中蝶形算法的基本思想分析（1）我們知道N點(diǎn)FFT運(yùn)算可以分成Iog2 （ N ）級，每一級都有N/2 個碟形，F(xiàn)FT的基本思想是用3層循環(huán)完成全部運(yùn)算（N點(diǎn)FFT）。（2）第一層循環(huán)：由于 N=2m 需要m級計算，第一層循環(huán)對運(yùn)算的 級數(shù)進(jìn)行控制。（stages）（3）第二層循環(huán)：由于第L級有2YL-1）個蝶形因子（乘數(shù)），第二層循環(huán)根據(jù)乘數(shù)進(jìn)行控制，保證對于每一個蝶形因子第三層

5、循環(huán)要執(zhí)行 一次，這樣，第三層循環(huán)在第二層循環(huán)控制下，每一級要進(jìn)行2A（L-1）次循環(huán)計算.（選擇W）（4）第三層循環(huán)：由于第L級共有N/2AL即2A （n-L）個群，并且同一級內(nèi)不同群的乘數(shù)分布相同，當(dāng)?shù)诙友h(huán)確定某一乘數(shù)后，第三層循環(huán)要將本級中每個群中具有這一乘數(shù)的蝶形計算一次，即第三層循環(huán)每執(zhí)行完一次要進(jìn)行 N/2AL個碟形計算。（執(zhí)行不同group中具有 相同W的蝶形運(yùn)算）(5 )可以得出結(jié)論：在每一級中，第三層循環(huán)完成N/2U 個碟形計算;第二層循環(huán)使第三層循環(huán)進(jìn)行2YL-1)次，因此，第二層循環(huán)完成時，個碟形計算。實(shí)質(zhì)是：第二、第三層循共進(jìn)行 2A(L-1) *N/2AL=N/2

6、環(huán)完成了第L級的計算。五、用c語言實(shí)現(xiàn)的FFT算法如下:<span style="font-size:18px;">#include <stdio.h>#in clude <math.h>#i nclude <stdlib.h>#defi ne N 1000/*定義復(fù)數(shù)類型*/typedef structdouble real;double img;complex;complex xN, *W; /*輸入序列，變換核*/int size_x=0; /*輸入序列的大小，在本程序中僅限2的次幕*/double PI;/*圓周率*/

7、void fft(); /*快速傅里葉變換*/void ini tW(); /*初始化變換核*/void cha nge(); /*變址*/void add(complex ,complex ,complex *); /*復(fù)數(shù)加法*/void mul(complex ,complex ,complex *); /*復(fù)數(shù)乘法*/void sub(complex ,complex ,complex *); /*復(fù)數(shù)減法*/void output();/*輸出快速傅里葉變換的結(jié)果*/int mai n()int i;system("cls");PI=ata n*4; printf

8、("果n");/*輸出結(jié)果*/輸岀DIT方法實(shí)現(xiàn)的FFT結(jié)輸入序列的大小輸入序列的實(shí)部和虛部prin tf("Please in put the size of x: n");/ sea nf("%d", &size_x);printf("Please in put the data in xN:n");/ for(i=0;i<size_x;i+)printf(" 請輸入第%d個序列：",i); sea nf("%lf%lf", &xi.real, &a

9、mp;xi.img);printf(“輸岀倒序后的序列n");ini tW();/ 調(diào)用變換核fft();/調(diào)用快速傅里葉變換printf("輸岀FFT后的結(jié)果n");output();/調(diào)用輸岀傅里葉變換結(jié)果函數(shù)return 0;/*快速傅里葉變換*/void fft()int i=O,j=O,k=O,l=O;complex up,dow n,product;eha nge(); / 調(diào)用變址函數(shù) for(i=0;i< log(size_x)/log(2) ;i+)/*l=1<<i;for(j=0;j<size_x;j+= 2*l )/*

10、p的蝶形因子乘數(shù)不同*/for(k=0;k<l;k+)/*蝶形運(yùn)算*/一級蝶形運(yùn)算stage */一組蝶形運(yùn)算group,每個grou一個蝶形運(yùn)算每個group內(nèi)的mul(xj+k+l,Wsize_x*k/2/l,&product); add(xj+k,product，&up);sub(xj+k,product，& dow n);xj+k=up; xj+k+l=dow n;/*初始化變換核，定義一個變換核，相當(dāng)于旋轉(zhuǎn)因子WAP*/void ini tW() int i;生成變換核用歐拉公式計算旋轉(zhuǎn)因子W=(complex *)malloc(sizeof(compl

11、ex) * size_x); / for(i=0;i<size_x;i+)Wi.real=cos(2*PI/size_x*i); / Wi.img=-1*si n(2*PI/size_x*i);/*變址計算，將 x(n)碼位倒置*/void cha nge()complex temp;un sig ned short i=O,j=O,k=O;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)

12、;k=k>>1;if(j>i) /將x(n)的碼位互換temp=xi;xi=xjxj=temp;output();/*輸出傅里葉變換的結(jié)果*/void output()int i;printf("The result are as follows: n");for(i=0;i<size_x;i+)prin tf("%.4f',xi.real);if(xi.img>=0.0001)pri ntf("+%.4fjn",xi.img);else if(fabs(xi.img)<0.0001)pr in tf

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