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

./串行FFT遞歸算法〔蝶式遞歸計算原理求傅里葉變換摘要FFT,即為快速傅氏變換,是離散傅氏變換的快速算法,它是根據(jù)離散傅氏變換的奇、偶、虛、實等特性,對離散傅立葉變換的算法進(jìn)行改進(jìn)獲得的。它對傅氏變換的理論并沒有新的發(fā)現(xiàn),但是對于在計算機系統(tǒng)或者說數(shù)字系統(tǒng)中應(yīng)用離散傅立葉變換,可以說是進(jìn)了一大步。設(shè)x<n>為N項的復(fù)數(shù)序列,由DFT變換,任一X〔m的計算都需要N次復(fù)數(shù)乘法和N-1次復(fù)數(shù)加法,而一次復(fù)數(shù)乘法等于四次實數(shù)乘法和兩次實數(shù)加法,一次復(fù)數(shù)加法等于兩次實數(shù)加法,即使把一次復(fù)數(shù)乘法和一次復(fù)數(shù)加法定義成一次"運算"〔四次實數(shù)乘法和四次實數(shù)加法,那么求出N項復(fù)數(shù)序列的X〔m,即N點DFT變換大約就需要N^2次運算。當(dāng)N=1024點甚至更多的時候,需要N2=1048576次運算,在FFT中,利用WN的周期性和對稱性,把一個N項序列〔設(shè)N=2k,k為正整數(shù),分為兩個N/2項的子序列,每個N/2點DFT變換需要〔N/2^2次運算,再用N次運算把兩個N/2點的DFT變換組合成一個N點的DFT變換。這樣變換以后,總的運算次數(shù)就變成N+2〔N/2^2=N+N^2/2。繼續(xù)上面的例子,N=1024時,總的運算次數(shù)就變成了525312次,節(jié)省了大約50%的運算量。而如果我們將這種"一分為二"的思想不斷進(jìn)行下去,直到分成兩兩一組的DFT運算單元,那么N點的DFT變換就只需要Nlog<2><N>次的運算,N在1024點時,運算量僅有10240次,是先前的直接算法的1%,點數(shù)越多,運算量的節(jié)約就越大,這就是FFT的優(yōu)越性。關(guān)鍵字：FFT蝶式計算傅里葉變換.目錄一．題目及要求11.1題目1二．設(shè)計算法、算法原理12.1算法原理與設(shè)計12.2設(shè)計步驟2三．算法描述、設(shè)計流程43.1算法描述43.2流程圖5四．源程序代碼及運行結(jié)果74.1源程序代碼74.2運行結(jié)果12五．算法分析、優(yōu)缺點135.1算法分析135.2優(yōu)缺點14六．總結(jié)15七．參考文獻(xiàn)16.一．題目及要求1.1題目對給定的,利用串行FFT遞歸算法〔蝶式遞歸計算原理計算其傅里葉變換的結(jié)果。1.2要求利用串行遞歸與蝶式遞歸原理,對給定的向量求解傅里葉變換的結(jié)果。二．設(shè)計算法、算法原理2.1算法原理與設(shè)計蝶式遞歸計算原理：令為n/2次單位元根,則有,將b向量的偶數(shù)項和奇數(shù)項分別記為和。注意推導(dǎo)中反復(fù)使用：。圖2.1公式圖形2.2設(shè)計步驟對于以上的分析可畫出如圖2.2所示的離散傅里葉變換遞歸計算流圖。圖2.3就是一個按此遞歸方法計算的n=8的FFT蝶式計算圖。FFT的蝶式遞歸計算圖〔有計算原理推出：圖2.2遞歸計算流圖特別的,ｎ＝８的FFT蝶式計算圖〔展開的：圖2.3蝶式計算圖按輸入元素展開,前面將輸出序列之元素按其偶下標(biāo)〔和〔展開,導(dǎo)出和遞歸計算式,按此構(gòu)造出了如圖1所示的FFT遞歸計算流程圖。事實上,我們也可以將輸入序列之元素按其偶下標(biāo)〔和幾下標(biāo)〔展開,則導(dǎo)出另一種形式的FFT遞歸計算式。三．算法描述、設(shè)計流程3.1算法描述SISD上的FFT分治遞歸算法：輸入:a=<a0,a1,…,an-1>;輸出:B=<b0,b1,…,bn-1>ProcedureRFFT<a,b>beginifn=1thenb0=a0else<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>forj=0ton-1do<4.1>bj=ujmodn/2+zvjmodn/2<4.2>z=zωendforendifend注:<1>算法時間復(fù)雜度t<n>=2t<n/2>+O<n>t<n>=O<nlogn>n=8的FFT蝶式計算圖:圖3.1FFT蝶式計算圖n=6的FFT遞歸計算流程圖:圖3.2FFT遞歸計算流程圖3.2流程圖開始開始計算出前size_x/2個exp<-j*2π*k/size_x>個值,即W的值輸入序列對應(yīng)值〔例如5+j3,輸入53輸入序列長度size_x計算出前size_x/2個exp<-j*2π*k/size_x>個值,即W的值輸入序列對應(yīng)值〔例如5+j3,輸入53輸入序列長度size_x級數(shù)i>=QUOTE級數(shù)i>=QUOTE？級數(shù)i加1是輸出fft結(jié)果序列輸出fft結(jié)果序列結(jié)束否結(jié)束該級該組起始下標(biāo)j>=QUOTE該級該組起始下標(biāo)j>=QUOTE？計算出該級需要的W的個數(shù)l否該級該組元素序數(shù)k>=QUOTE該級該組元素序數(shù)k>=QUOTE？組起始下標(biāo)加2*lK加1是K加1X[j+k]X[j+k]lX[j+k+l]*W[<size_x/2/l>*k]X[j+k+l]X[j+k]X[j+k]lX[j+k+l]*W[<size_x/2/l>*k]X[j+k+l]-1四．源程序代碼及運行結(jié)果4.1源程序代碼/************FFT***********/#include<stdio.h>//整個程序輸入和輸出利用同一個空間x[N],節(jié)約空間#include<math.h>#include<stdlib.h>#defineN1000//定義輸入或者輸出空間的最大長度typedefstruct{doublereal;doubleimg;}complex;//定義復(fù)數(shù)型變量的結(jié)構(gòu)體voidfft<>;//快速傅里葉變換函數(shù)聲明voidinitW<>;//計算W<0>~W<size_x-1>的值函數(shù)聲明voidchange<>;//碼元位置倒置函數(shù)函數(shù)聲明voidadd<complex,complex,complex*>;/*復(fù)數(shù)加法*/voidmul<complex,complex,complex*>;/*復(fù)數(shù)乘法*/voidsub<complex,complex,complex*>;/*復(fù)數(shù)減法*/voiddivi<complex,complex,complex*>;/*復(fù)數(shù)除法*/voidoutput<>;/*輸出結(jié)果*/complexx[N],*W;/*輸出序列的值*/intsize_x=0;/*輸入序列的長度,只限2的N次方*/doublePI;//pi的值intmain<>{inti;system<"cls">;PI=atan<1>*4;printf<"Pleaseinputthesizeofx:\n">;/*輸入序列的長度*/scanf<"%d",&size_x>;printf<"Pleaseinputthedatainx[N]:<suchas:56>\n">;for<i=0;i<size_x;i++>/*輸入序列對應(yīng)的值*/scanf<"%lf%lf",&x[i].real,&x[i].img>;initW<>;//計算W<0>~W<size_x-1>的值fft<>;//利用fft快速算法進(jìn)行DFT變化output<>;//順序輸出size_x個fft的結(jié)果return0;}/*進(jìn)行基-2FFT運算,蝶形算法。這個算法的思路就是,先把計算過程分為log<size_x>/log<2>-1級〔用i控制級數(shù)；然后把每一級蝶形單元分組〔用j控制組的第一個元素起始下標(biāo)；最后算出某一級某一組每一個蝶形單元〔用k控制個數(shù),共l個。*/voidfft<>{inti=0,j=0,k=0,l=0;complexup,down,product;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>{for<k=0;k<l;k++>//算出第i級內(nèi)j組蝶形單元的結(jié)果{//算出j組中第k個蝶形單元mul<x[j+k+l],W[<size_x/2/l>*k],&product>;/*size/2/l是該級W的相鄰上標(biāo)差,l是該級該組取的W總個數(shù)*/add<x[j+k],product,&up>;sub<x[j+k],product,&down>;x[j+k]=up;//up為蝶形單元右上方的值x[j+k+l]=down;//down為蝶形單元右下方的值}}}}voidinitW<>//計算W的實現(xiàn)函數(shù){inti;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個值即可*/{W[i].real=cos<2*PI/size_x*i>;//計算W的實部W[i].img=-1*sin<2*PI/size_x*i>;//計算W的虛部}}voidchange<>//輸入的碼組碼位倒置實現(xiàn)函數(shù){complextemp;unsignedshorti=0,j=0,k=0;doublet;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=x[i];x[i]=x[j];x[j]=temp;}}}voidoutput<>//輸出結(jié)果實現(xiàn)函數(shù){inti;printf<"Theresultareasfollows\n">;for<i=0;i<size_x;i++>{printf<"%.4f",x[i].real>;//輸出實部if<x[i].img>=0.0001>//如果虛部的值大于0.0001,輸出+jx.img的形式printf<"+j%.4f\n",x[i].img>;elseif<fabs<x[i].img><0.0001>printf<"\n">;elseprintf<"-j%.4f\n",fabs<x[i].img>>;//如果虛部的值小于-0.0001,輸出-jx.img的形式}}voidadd<complexa,complexb,complex*c>//復(fù)數(shù)加法實現(xiàn)函數(shù){c->real=a.real+b.real;//復(fù)數(shù)實部相加c->img=a.img+b.img;//復(fù)數(shù)虛部相加}voidmul<complexa,complexb,complex*c>//復(fù)數(shù)乘法實現(xiàn)函數(shù){c->real=a.real*b.real-a.img*b.img;//獲取相乘結(jié)果的實部c->img=a.real*b.img+a.img*b.real;//獲取相乘結(jié)果的虛部}voidsub<complexa,complexb,complex*c>//復(fù)數(shù)減法實現(xiàn)函數(shù){c->real=a.real-b.real;//復(fù)數(shù)實部相減c->img=a.img-b.img;//復(fù)數(shù)虛部相減}voiddivi<complexa,complexb,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é)果的虛部}4.2運行結(jié)果〔1處理器p=8：圖4.1當(dāng)時串行FFT輸出結(jié)果〔2處理器p=8：當(dāng)時輸出結(jié)果與計算結(jié)果相符如圖4.2所示圖4.2運行圖五．算法分析、優(yōu)缺點5.1算法分析〔1FFT算法的基本原理是把長序列的DFT逐次分解為較短序列的DFT。按照抽取方式的不同可分為DIT-FFT〔按時間抽取和DIF-FFT〔按頻率抽取算法。按照蝶形運算的構(gòu)成不同可分為基2、基4、基8以及任意因子〔2n,n為大于1的整數(shù),基2、基4算法較為常用?！?總體結(jié)構(gòu)說明輸入數(shù)據(jù)為串行的數(shù)據(jù)流,故在第一級蝶形運算模塊前加入串并轉(zhuǎn)換模塊,將串行數(shù)據(jù)流轉(zhuǎn)換為并行的兩列數(shù)據(jù)流以適應(yīng)基2蝶形運算模塊的輸入信號要求。由于每級蝶形運算一次處理的兩個輸入數(shù)據(jù)不能直接由前一級蝶形運算一次性輸出,故在兩個蝶形運算單元之間插入延時對齊模塊,將前一級蝶形運算的結(jié)果〔兩列并行的數(shù)據(jù)流作適當(dāng)?shù)难訒r并通過轉(zhuǎn)接器對齊,形成后一級蝶形運算模塊所需要的2列輸入序列。在最后一級蝶形運算后加入串并轉(zhuǎn)換模塊,將2列并行的數(shù)據(jù)流合成為1列。最后加入倒序模塊將DIF-FFT得到的倒序輸出序列整理為順序輸出。旋轉(zhuǎn)因子產(chǎn)生模塊產(chǎn)生各級基2蝶形運算所需的旋轉(zhuǎn)因子。由運算流圖可以看出最后一級的旋轉(zhuǎn)因子其實是1,故可省略最后一級蝶形運算單元中的旋轉(zhuǎn)因子乘法器。因此用一個雙口ROM將兩組數(shù)據(jù)分別輸出到第一級和第二級的蝶形運算單元即可?；?蝶形運算模塊由兩個復(fù)數(shù)加法器和一個復(fù)數(shù)乘法器構(gòu)成。旋轉(zhuǎn)因子由ROM產(chǎn)生后,作為復(fù)數(shù)乘法器的輸入之一,與前面復(fù)數(shù)加法器得到的結(jié)果相乘完成一次蝶形運算。為提高系統(tǒng)的運行速度可在蝶形運算單元中插入流水線寄存中間結(jié)果?！?蝶形運算單元如下所示：圖5.1運算單元5.2優(yōu)缺點優(yōu)點：結(jié)構(gòu)清晰,可讀性強,而且容易用數(shù)學(xué)歸納法來證明算