串行FFT遞歸算法蝶式遞歸計(jì)算原理求傅里葉變換

1、串行FFT遞歸算法（蝶式遞歸計(jì)算原理）求傅里葉變換摘要 FFT，即為快速傅氏變換，是離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的奇、偶、虛、實(shí)等特性，對(duì)離散傅立葉變換的算法進(jìn)行改進(jìn)獲得的。它對(duì)傅氏變換的理論并沒(méi)有新的發(fā)現(xiàn)，但是對(duì)于在計(jì)算機(jī)系統(tǒng)或者說(shuō)數(shù)字系統(tǒng)中應(yīng)用離散傅立葉變換，可以說(shuō)是進(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ù)加法）

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

3、N點(diǎn)的DFT變換就只需要Nlog(2)(N)次的運(yùn)算，N在1024點(diǎn)時(shí)，運(yùn)算量?jī)H有10240次，是先前的直接算法的1%，點(diǎn)數(shù)越多，運(yùn)算量的節(jié)約就越大，這就是FFT的優(yōu)越性。關(guān)鍵字：FFT 蝶式計(jì)算 傅里葉變換16 / 18目錄一題目與要求11.1題目1二設(shè)計(jì)算法、算法原理12.1算法原理與設(shè)計(jì)12.2設(shè)計(jì)步驟2三算法描述、設(shè)計(jì)流程43.1算法描述43.2流程圖5四源程序代碼與運(yùn)行結(jié)果74.1源程序代碼74.2運(yùn)行結(jié)果12五算法分析、優(yōu)缺點(diǎn)135.1算法分析135.2優(yōu)缺點(diǎn)14六總結(jié)15七參考文獻(xiàn)16一題目與要求1.1題目對(duì)給定的，利用串行FFT遞歸算法（蝶式遞歸計(jì)算原理）計(jì)算其傅里葉變換的結(jié)果

4、。1.2要求利用串行遞歸與蝶式遞歸原理，對(duì)給定的向量求解傅里葉變換的結(jié)果。二設(shè)計(jì)算法、算法原理2.1算法原理與設(shè)計(jì)蝶式遞歸計(jì)算原理：令 為n/2次單位元根，則有 ，將b向量的偶數(shù)項(xiàng) 和奇數(shù)項(xiàng) 分別記 為 和 。注意推導(dǎo)中反復(fù)使用： 。圖2.1 公式圖形2.2設(shè)計(jì)步驟對(duì)于以上的分析可畫出如圖 2.2所示的離散傅里葉變換遞歸計(jì)算流圖。圖2.3就是一個(gè)按此遞歸方法計(jì)算的n=8的FFT蝶式計(jì)算圖。FFT的蝶式遞歸計(jì)算圖（有計(jì)算原理推出）：圖2.2 遞歸計(jì)算流圖特別的，的FFT蝶式計(jì)算圖（展開的）：圖2.3 蝶式計(jì)算圖按輸入元素展開，前面將輸出序列之元素 按其偶下標(biāo)（）和（）展開，導(dǎo)出 和遞歸計(jì)算式，按

5、此構(gòu)造出了如圖1所示的FFT遞歸計(jì)算流程圖。事實(shí)上，我們也可以將輸入序列之元素按其偶下標(biāo)（） 和幾下標(biāo)（）展開，則導(dǎo)出另一種形式的FFT遞歸計(jì)算式 。三算法描述、設(shè)計(jì)流程3.1算法描述SISD上的FFT分治遞歸算法：輸入: a=(a0,a1,an-1); 輸出: B=(b0,b1,bn-1) Procedure RFFT(a,b) begin if n=1 then b0=a0 else (1)RFFT(a0,a2,an-2, u0,u1,un/2-1) (2)RFFT(a1,a3,an-1, v0,v1,vn/2-1) (3)z=1 (4)for j=0 to n-1 do (4.1)bj=

6、uj mod n/2+zvj mod n/2 (4.2)z=z endfor endifend 注: (1)算法時(shí)間復(fù)雜度t(n)=2t(n/2)+O(n) t(n)=O(nlogn)n=8的FFT蝶式計(jì)算圖:圖3.1 FFT蝶式計(jì)算圖n=6的FFT遞歸計(jì)算流程圖:圖3.2 FFT遞歸計(jì)算流程圖3.2流程圖開始計(jì)算出前size_x/2個(gè)exp(-j*2*k/size_x)個(gè)值，即W的值輸入序列對(duì)應(yīng)值（例如5+j3，輸入5 3）輸入序列長(zhǎng)度size_x飛級(jí)數(shù)i>=？級(jí)數(shù)i加1是 輸出fft結(jié)果序列結(jié)束 否 該級(jí)該組起始下標(biāo)j>=？計(jì)算出該級(jí)需要的W的個(gè)數(shù)l 是 否該級(jí)該組元素序數(shù)k&

7、gt;=？組起始下標(biāo)加2*lK加1 是Xj+k Xj+klXj+k+l*W(size_x/2/l)*k Xj+k+l -1 否四源程序代碼與運(yùn)行結(jié)果4.1源程序代碼/*FFT*/#include <stdio.h> /整個(gè)程序輸入和輸出利用同一個(gè)空間xN，節(jié)約空間#include <math.h>#include <stdlib.h>#define N 1000 /定義輸入或者輸出空間的最大長(zhǎng)度typedefstructdouble real;doubleimg;complex; /定義復(fù)數(shù)型變量的結(jié)構(gòu)體void fft(); /快速傅里葉變換函數(shù)聲明voi

8、d initW(); /計(jì)算W(0)W(size_x-1)的值函數(shù)聲明void 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; /*輸出序列的值*/intsize_

9、x=0; /*輸入序列的長(zhǎng)度，只限2的N次方*/double PI; /pi的值int main()inti;system("cls");PI=atan(1)*4;printf("Please input the size of x:n"); /*輸入序列的長(zhǎng)度*/scanf("%d",&size_x);printf("Please input the data in xN:(such as:5 6)n");for(i=0;i<size_x;i+) /*輸入序列對(duì)應(yīng)的值*/scanf("%l

10、f %lf",&xi.real,&xi.img);initW(); /計(jì)算W(0)W(size_x-1)的值fft(); /利用fft快速算法進(jìn)行DFT變化output(); /順序輸出size_x個(gè)fft的結(jié)果return 0;/*進(jìn)行基-2 FFT運(yùn)算,蝶形算法。這個(gè)算法的思路就是，先把計(jì)算過(guò)程分為log(size_x)/log(2)-1級(jí)（用i控制級(jí)數(shù)）；然后把每一級(jí)蝶形單元分組（用j控制組的第一個(gè)元素起始下標(biāo)）；最后算出某一級(jí)某一組每一個(gè)蝶形單元（用k控制個(gè)數(shù)，共l個(gè)）。*/voidfft()inti=0,j=0,k=0,l=0;complexup,down,

11、product;change(); /實(shí)現(xiàn)對(duì)碼位的倒置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)for(k=0;k<l;k+) /算出第i級(jí)j組蝶形單元的結(jié)果 /算出j組中第k個(gè)蝶形單元mul(xj+k+l,W(size_x/2/l)*k,&product); /*size/2/l是該級(jí)W的相鄰上標(biāo)差，l是該級(jí)該組取的W總個(gè)數(shù)*/add(xj+k,product,&up);sub(xj+k,product,&down);xj+k=up

12、; /up為蝶形單元右上方的值xj+k+l=down; /down為蝶形單元右下方的值void initW() /計(jì)算W的實(shí)現(xiàn)函數(shù)inti;W=(complex *)malloc(sizeof(complex) * size_x);/*申請(qǐng)size_x個(gè)復(fù)數(shù)W的空間(這部申請(qǐng)的空間有點(diǎn)多，實(shí)際上只要申請(qǐng)size_x/2個(gè)即可)*/for(i=0;i<(size_x/2);i+)/*預(yù)先計(jì)算出size_x/2個(gè)W的值，存放,由于蝶形算法只需要前size_x/2個(gè)值即可*/Wi.real=cos(2*PI/size_x*i); /計(jì)算W的實(shí)部Wi.img=-1*sin(2*PI/size_x

13、*i); /計(jì)算W的虛部void change() /輸入的碼組碼位倒置實(shí)現(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=temp;void output() /輸出結(jié)果實(shí)現(xiàn)函數(shù)inti;printf("The result are as followsn&qu

14、ot;);for(i=0;i<size_x;i+)printf("%.4f",xi.real); /輸出實(shí)部if(xi.img>=0.0001) /如果虛部的值大于0.0001，輸出+jx.img的形式printf("+j%.4fn",xi.img);else if(fabs(xi.img)<0.0001)printf("n");elseprintf("-j%.4fn",fabs(xi.img);/如果虛部的值小于-0.0001，輸出-jx.img的形式void add(complex a,com

15、plexb,complex *c) /復(fù)數(shù)加法實(shí)現(xiàn)函數(shù)c->real = a.real + b.real; /復(fù)數(shù)實(shí)部相加c->img = a.img + b.img; /復(fù)數(shù)虛部相加void mul(complex a,complexb,complex *c) /復(fù)數(shù)乘法實(shí)現(xiàn)函數(shù)c->real = a.real*b.real - a.img*b.img; /獲取相乘結(jié)果的實(shí)部c->img = a.real*b.img + a.img*b.real; /獲取相乘結(jié)果的虛部void sub(complex a,complexb,complex *c) /復(fù)數(shù)減法實(shí)現(xiàn)函數(shù)c

16、->real = a.real - b.real; /復(fù)數(shù)實(shí)部相減c->img = a.img - b.img; /復(fù)數(shù)虛部相減void divi(complex a,complexb,complex *c) /復(fù)數(shù)除法實(shí)現(xiàn)函數(shù)c->real = (a.real*b.real + a.img*b.img) / (b.real*b.real+b.img*b.img);/獲取相除結(jié)果的實(shí)部c->img = (a.img*b.real - a.real*b.img) / (b.real*b.real+b.img*b.img);/獲取相除結(jié)果的虛部4.2運(yùn)行結(jié)果（1）處理器p=

17、8：圖4.1當(dāng)時(shí)串行FFT輸出結(jié)果（2）處理器p=8：當(dāng)時(shí)輸出結(jié)果與計(jì)算結(jié)果相符如圖4.2所示圖4.2運(yùn)行圖五 算法分析、優(yōu)缺點(diǎn)5.1算法分析（1）FFT算法的基本原理是把長(zhǎng)序列的DFT逐次分解為較短序列的DFT。按照抽取方式的不同可分為DIT-FFT（按時(shí)間抽?。┖虳IF-FFT（按頻率抽取）算法。按照蝶形運(yùn)算的構(gòu)成不同可分為基2、基4、基8以與任意因子（2n,n為大于1的整數(shù)），基2、基4算法較為常用。（2）總體結(jié)構(gòu)說(shuō)明 輸入數(shù)據(jù)為串行的數(shù)據(jù)流，故在第一級(jí)蝶形運(yùn)算模塊前加入串并轉(zhuǎn)換模塊，將串行數(shù)據(jù)流轉(zhuǎn)換為并行的兩列數(shù)據(jù)流以適應(yīng)基2蝶形運(yùn)算模塊的輸入信號(hào)要求。 由于每級(jí)蝶

18、形運(yùn)算一次處理的兩個(gè)輸入數(shù)據(jù)不能直接由前一級(jí)蝶形運(yùn)算一次性輸出，故在兩個(gè)蝶形運(yùn)算單元之間插入延時(shí)對(duì)齊模塊，將前一級(jí)蝶形運(yùn)算的結(jié)果（兩列并行的數(shù)據(jù)流）作適當(dāng)?shù)难訒r(shí)并通過(guò)轉(zhuǎn)接器對(duì)齊，形成后一級(jí)蝶形運(yùn)算模塊所需要的2列輸入序列。 在最后一級(jí)蝶形運(yùn)算后加入串并轉(zhuǎn)換模塊，將2列并行的數(shù)據(jù)流合成為1列。最后加入倒序模塊將DIF-FFT得到的倒序輸出序列整理為順序輸出。 旋轉(zhuǎn)因子產(chǎn)生模塊產(chǎn)生各級(jí)基2蝶形運(yùn)算所需的旋轉(zhuǎn)因子。由運(yùn)算流圖可以看出最后一級(jí)的旋轉(zhuǎn)因子其實(shí)是1，故可省略最后一級(jí)蝶形運(yùn)算單元中的旋轉(zhuǎn)因子乘法器。因此用一個(gè)雙口ROM將兩組數(shù)據(jù)分別輸出到第一級(jí)和第二級(jí)的蝶形運(yùn)算單元即可。

19、 基2蝶形運(yùn)算模塊由兩個(gè)復(fù)數(shù)加法器和一個(gè)復(fù)數(shù)乘法器構(gòu)成。旋轉(zhuǎn)因子由ROM產(chǎn)生后，作為復(fù)數(shù)乘法器的輸入之一，與前面復(fù)數(shù)加法器得到的結(jié)果相乘完成一次蝶形運(yùn)算。為提高系統(tǒng)的運(yùn)行速度可在蝶形運(yùn)算單元中插入流水線寄存中間結(jié)果。（3）蝶形運(yùn)算單元如下所示：圖5.1運(yùn)算單元5.2優(yōu)缺點(diǎn)優(yōu)點(diǎn)：結(jié)構(gòu)清晰，可讀性強(qiáng)，而且容易用數(shù)學(xué)歸納法來(lái)證明算法的正確性，因此它為設(shè)計(jì)算法、調(diào)試程序帶來(lái)很大方便。缺點(diǎn)：遞歸算法的運(yùn)行效率較低，無(wú)論是耗費(fèi)的計(jì)算時(shí)間還是占用的存儲(chǔ)空間都比非遞歸算法要多。六總結(jié)經(jīng)過(guò)一周的課程設(shè)計(jì)，使我更加深刻的學(xué)習(xí)課本知識(shí)，又復(fù)習(xí)鞏固了以前學(xué)過(guò)的知識(shí)。雖然這次課程是那么短暫的一周時(shí)間，但是我感覺(jué)