FFT最詳細(xì)的源代碼和解釋.docx

我自己的一些詳細(xì)標(biāo)注，有利于深入了解FFT，后面附加幾位網(wǎng)友對(duì)FFT的理解及源代碼，讓廣大朋友更迅速的掌握FFT#include DSP281x_Device.h / DSP281x Headerfile Include File，添加所有頭文件#include DSP281x_Examples.h / DSP281x Examples Include File，條件編譯而已#include f2812a.h /一些變量的宏定義而已#includemath.h#define PI 3.1415926 /前變后常#define SAMPLENUMBER 128/#define SAMPLENUMBER 512void InitForFFT();void MakeWave();/void FFT(float dataRSAMPLENUMBER,float dataISAMPLENUMBER);int INPUTSAMPLENUMBER,DATASAMPLENUMBER;float fWaveRSAMPLENUMBER,fWaveISAMPLENUMBER,wSAMPLENUMBER;float sin_tabSAMPLENUMBER,cos_tabSAMPLENUMBER;/逐級(jí)計(jì)算FFT，一級(jí)一級(jí)遞推void FFT(float dataRSAMPLENUMBER,float dataISAMPLENUMBER)int x0,x1,x2,x3,x4,x5,x6,xx;int i,j,k,b,p,L;float TR,TI,temp;/* following code invert sequence */倒序for ( i=0;iSAMPLENUMBER;i+ )/就是碼位倒置嘛，二進(jìn)制各個(gè)位獨(dú)立出來再反向 /128七位二進(jìn)制表示，/號(hào)代表右移嘛x0=x1=x2=x3=x4=x5=x6=0;x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01;x4=(i/16)&0x01; x5=(i/32)&0x01; x6=(i/64)&0x01;xx=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;/最低位，最高位反過來dataIxx=dataRi;for ( i=0;iSAMPLENUMBER;i+ )dataRi=dataIi; dataIi=0; /對(duì)應(yīng)過來/* following code FFT */for ( L=1;L0 ) b=b*2; i-; /* b= 2(L-1) */for ( j=0;j0 ) /* p=pow(2,7-L)*j; */p=p*2; i-;p=p*j;for ( k=j;k128;k=k+2*b ) /* for (3) */TR=dataRk; TI=dataIk; temp=dataRk+b;dataRk=dataRk+dataRk+b*cos_tabp+dataIk+b*sin_tabp;dataIk=dataIk-dataRk+b*sin_tabp+dataIk+b*cos_tabp;dataRk+b=TR-dataRk+b*cos_tabp-dataIk+b*sin_tabp;dataIk+b=TI+temp*sin_tabp-dataIk+b*cos_tabp; /遞推嘛，防止立馬調(diào)用結(jié)果， /* END for (3) */ /引入一個(gè)中間變量存原始值， /* END for (2) */ /防止上一步對(duì)下一步的影響 /* END for (1) */for ( i=0;iSAMPLENUMBER/2;i+ ) /對(duì)稱性，前半部分即可 wi=sqrt(dataRi*dataRi+dataIi*dataIi); /* END FFT */main()int i;InitForFFT();MakeWave();for ( i=0;iSAMPLENUMBER;i+ )fWaveRi=INPUTi;fWaveIi=0.0f;wi=0.0f;FFT(fWaveR,fWaveI);/輸入波形進(jìn)行FFT變換,此處引入起始實(shí)參即可遞推下去for ( i=0;iSAMPLENUMBER;i+ )DATAi=wi;/變換后的波形轉(zhuǎn)換到輸出接口while ( 1 );/ break point/旋轉(zhuǎn)因子事先初始化好，方便調(diào)用void InitForFFT()int i;for ( i=0;iSAMPLENUMBER;i+ )sin_tabi=sin(PI*2*i/SAMPLENUMBER);cos_tabi=cos(PI*2*i/SAMPLENUMBER);/旋轉(zhuǎn)因子事先初始化好，方便調(diào)用 /利用這個(gè)，確實(shí)能產(chǎn)生各種各樣的諧波void MakeWave() /利用這個(gè)，確實(shí)能產(chǎn)生各種各樣的諧波int i;for ( i=0;iSAMPLENUMBER;i+ )/1024是相應(yīng)的幅值嘛，只要弄出個(gè)基波，以之為標(biāo)準(zhǔn)即可INPUTi=sin(PI*2*i/SAMPLENUMBER)*1024 +sin(PI*2*i/SAMPLENUMBER*3)*1024/3 +sin(PI*2*i/SAMPLENUMBER*5)*1024/5 +sin(PI*2*i/SAMPLENUMBER*7)*1024/7 +sin(PI*2*i/SAMPLENUMBER*9)*1024/9; /INPUTi=sin(PI*2*i/SAMPLENUMBER*3)*1024;FFT原理及實(shí)現(xiàn)(Radix-2)哈!經(jīng)過連續(xù)幾個(gè)晚上的奮戰(zhàn),終于弄懂了FFT推導(dǎo)過程及實(shí)現(xiàn)!HappyJ基2 FFT總的思想是將輸入信號(hào)對(duì)半分割,再對(duì)半分割,再再對(duì)半分割(以下省略10000個(gè)再再.J)直至分割到2點(diǎn).兩點(diǎn)DFT簡化假設(shè)輸入為x0,x1;輸出為X0,X1.偽代碼如下:/ -#defineN2#definePI3.1415926/ -inti, jfor(i=0, Xi=0.0; iN; i+)for(j=0; jN; j+)Xi += xj * ( cos(2*PI*i*j/N) - sin(2*PI*i*j/N) );注意到(我想Audio編解碼很多時(shí)候都是對(duì)cos,sin進(jìn)行優(yōu)化!)j=0j=1i=0cos1sin0tw1cos1sin0tw1i=1cos1Sin0tw1cos-1sin0tw-1X0=x0*(1-0) + x1*(1-0)=x0 + 1*x1;X1=x0*(1-0) + x1*(-1-0) =x0 - 1*x1;這就是單個(gè)2點(diǎn)蝶形算法.FFT實(shí)現(xiàn)流程圖分析(N=8,以8點(diǎn)信號(hào)為例)FFT implementation of an 8-point DFT as two 4-point DFTs and four 2-point DFTs8點(diǎn)FFT流程圖(Layer表示層, gr表示當(dāng)前層的顆粒)下面以LayerI為例.LayerI部分,具有4個(gè)顆粒,每個(gè)顆粒2個(gè)輸入(注意2個(gè)輸入的來源,由時(shí)域信號(hào)友情提供,感謝感謝J)我們將輸入xk分為兩部分x_rk, x_ik.具有實(shí)部和虛部,時(shí)域信號(hào)本沒有虛部的,因此可以讓x_ik為0.那么為什么還要畫蛇添足分為實(shí)部和虛部呢?這是因?yàn)長ayerII, LayerIII的輸入是復(fù)數(shù),為了編碼統(tǒng)一而強(qiáng)行分的.當(dāng)然你編碼時(shí)可以判斷當(dāng)前層是否為1來決定是否分.但是我想每個(gè)人最后都會(huì)傾向分的.旋轉(zhuǎn)因子tw = cos(2*PI*k/N)-j*sin(2*PI*k/N);也可以分為實(shí)部和虛部,令其為tw_r, tw_i;則tw = tw_r - j*tw_i;Xk = (x_rk + j*x_ik) + (tw_rj*tw_i) * (x_rk+N/2+j*x_ik+N/2)則X_Rk = x_rk + tw_r*x_rk+N/2 + tw_i*x_ik+N/2;X_Ik = x_ik - tw_i*x_rk+N/2 + tw_r*x_ik+N/2;LayerII部分,具有2個(gè)顆粒,每個(gè)顆粒4個(gè)輸入(注意4個(gè)輸入的來源,由LayerI友情提供,感謝感謝J)LayerIII部分,具有1個(gè)顆粒,每個(gè)顆粒8個(gè)輸入(注意8個(gè)輸入的來源,由LayerII友情提供,感謝感謝J)LayerI, LayerII, LayerIII從左往右,蝶形信號(hào)運(yùn)算流非常明顯!假令輸入為xk, xk+N/2,輸出為Xk, Xk+N/2. xk分解為x_rk, x_ik部分則該蝶形運(yùn)算為Xk= (x_rk-j*x_ik) + (x_rk+N/2-j*x_ik+N/2)*(cos(2*PI*k/N)-j*sin(2*PI*k/N);再令cos(2*PI*k/N)為tw1, sin(2*PI*k/N)為tw2則Xk = (x_rk-j*x_ik) + (x_rk+N/2-j*x_ik+N/2)*(tw1-j*tw2);X_Rk = x_rk + x_rk+N/2*tw1 -x_ik+N/2*tw2;X_IK = x_ikx_rk = x_rk + x_rk+b*tw1 + x_ik+b*tw2;x_ik = x_ik - x_rk+b*tw2 + x_ik+b*tw1;譬如8點(diǎn)輸入x81.先分割成2部分:x0, x2, x4, x6和x1, x3, x5, x72.信號(hào)x0, x2, x4, x6再分割成x0, x4和x2, x6信號(hào)x1, x3, x5, x7再分割成x1, x5和x3, x73.無法分割了,已經(jīng)分割成2點(diǎn)了J.如上圖:在LayerI的時(shí)候,我們是對(duì)2點(diǎn)進(jìn)行DFT.(一共4次DFT )輸入為x0&x4;x2&x6;x1&x5;x3&x7輸出為y0,y1;Y2,y3;Y4,y5;Y6,y7;流程:I.希望將輸入直接轉(zhuǎn)換為x0, x4, x2, x6, x1, x5, x3, x7的順序II.對(duì)轉(zhuǎn)換順序后的信號(hào)進(jìn)行4次DFT步驟I代碼實(shí)現(xiàn)/*反轉(zhuǎn)算法.這個(gè)算法效率比較低!先用起來在說,之后需要進(jìn)行優(yōu)化.*/staticvoidbitrev(void)intp=1, q, i;intbit_rev N ;floatxx_r N ;bit_rev 0 = 0;while( p N )for(q=0; qp; q+)bit_rev q = bit_rev q * 2;bit_rev q + p = bit_rev q + 1;p *= 2;for(i=0; iN; i+)xx_r i = x_r i ;for(i=0; iN; i+)x_ri = xx_r bit_revi ;/ -此刻序列x重排完畢-步驟II代碼實(shí)現(xiàn)int j;float TR;/臨時(shí)變量float tw1;/旋轉(zhuǎn)因子/*兩點(diǎn)DFT */for(k=0; kN; k+=2)/兩點(diǎn)DFT簡化告訴我們tw1=1TR = x_rk;/ TR就是A, x_rk+b就是B.x_rk= TR + tw1*x_rk+b;x_rk+b = TR - tw1*x_rk+b;在LayerII的時(shí)候,我們希望得到z,就需要對(duì)y進(jìn)行DFT.y0,y2;y1,y3;y4,y6;y5,y7;z0,z1;z2,z3;z4,z5;z6,z7;在LayerIII的時(shí)候,我們希望得到v,就需要對(duì)z進(jìn)行DFT.z0,z4;z1,z5;z2,z6;z3,z7;v0,v1;v2,v3;v4,v5;v6,v7;準(zhǔn)備令輸入為xs, xs+N/2,輸出為ys, ys+N/2這個(gè)N絕對(duì)不是上面的8,這個(gè)N是當(dāng)前顆粒的輸入樣本總量對(duì)于LayerI而言N是2;對(duì)于LayerII而言N是4;對(duì)于LayerIII而言N是8復(fù)數(shù)乘法:(a+j*b) * (c+j*d)實(shí)部= a*c bd;虛部= ad + bc;旋轉(zhuǎn)因子:實(shí)現(xiàn)(C描述)#include#include#include/#include complex.h/ -#defineN8/64#defineM3/6/2m=N#definePI3.1415926/ -floattwiddleN/2 = 1.0, 0.707, 0.0, -0.707;floatx_rN = 1, 1, 1, 1, 0, 0, 0, 0;floatx_iN;/N=8/*floattwiddleN/2 = 1,0.9951, 0.9808, 0.9570, 0.9239, 0.8820, 0.8317, 0.7733,0.7075, 0.6349,0.5561,0.4721,0.3835,0.2912,0.1961,0.0991,0.0000,-0.0991,-0.1961,-0.2912,-0.3835,-0.4721,-0.5561,-0.6349,-0.7075,-0.7733, 0.8317,-0.8820,-0.9239,-0.9570,-0.9808,-0.9951;/N=64floatx_rN=1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,;floatx_iN;*/FILE *fp;/ - func -/*初始化輸出虛部*/staticvoidfft_init(void)inti;for(i=0; iN; i+)x_ii = 0.0;/*反轉(zhuǎn)算法.將時(shí)域信號(hào)重新排序.*這個(gè)算法有改進(jìn)的空間*/staticvoidbitrev(void)intp=1, q, i;intbit_rev N ;/floatxx_r N ;/bit_rev 0 = 0;while( p N )for(q=0; qp; q+)bit_rev q = bit_rev q * 2;bit_rev q + p = bit_rev q + 1;p *= 2;for(i=0; iN; i+)xx_r i = x_r i ;for(i=0; iN; i+)x_ri = xx_r bit_revi ;/* - add by sshc625 - */staticvoidbitrev2(void)return;/*/voiddisplay(void)printf(/n/n);inti;for(i=0; iN; i+)printf(%f/t%f/n, x_ri, x_ii);/*/voidfft1(void)fp = fopen(log1.txt,a+);intL, i, b, j, p, k, tx1, tx2;floatTR, TI, temp;/臨時(shí)變量floattw1, tw2;/*深M.對(duì)層進(jìn)行循環(huán). L為當(dāng)前層,總層數(shù)為M. */for(L=1; L 0)b *= 2;i-;/ -是否外層對(duì)顆粒循環(huán),內(nèi)層對(duì)樣本點(diǎn)循環(huán)邏輯性更強(qiáng)一些呢! -/* outter對(duì)參與DFT的樣本點(diǎn)進(jìn)行循環(huán)* L=1,循環(huán)了1次(4個(gè)顆粒,每個(gè)顆粒2個(gè)樣本點(diǎn))* L=2,循環(huán)了2次(2個(gè)顆粒,每個(gè)顆粒4個(gè)樣本點(diǎn))* L=3,循環(huán)了4次(1個(gè)顆粒,每個(gè)顆粒8個(gè)樣本點(diǎn))*/for(j=0; j 0)p = p*2;i-;p= p * j;tx1 = p % N;tx2 = tx1 + 3*N/4;tx2 = tx2 % N;/ tw1是cos部分,實(shí)部; tw2是sin部分,虛數(shù)部分.tw1 = ( tx1=N/2)? -twiddletx1-N/2: twiddle tx1 ;tw2 = ( tx2=N/2)? -twiddletx2-(N/2) : twiddletx2;/* inner對(duì)顆粒進(jìn)行循環(huán)* L=1,循環(huán)了4次(4個(gè)顆粒,每個(gè)顆粒2個(gè)輸入)* L=2,循環(huán)了2次(2個(gè)顆粒,每個(gè)顆粒4個(gè)輸入)* L=3,循環(huán)了1次(1個(gè)顆粒,每個(gè)顆粒8個(gè)輸入)*/for(k=j; kN; k=k+2*b)TR = x_rk;/ TR就是A, x_rk+b就是B.TI = x_ik;temp = x_rk+b;/*如果復(fù)習(xí)一下(a+j*b)(c+j*d)兩個(gè)復(fù)數(shù)相乘后的實(shí)部虛部分別是什么*就能理解為什么會(huì)如下運(yùn)算了,只有在L=1時(shí)候輸入才是實(shí)數(shù),之后層的*輸入都是復(fù)數(shù),為了讓所有的層的輸入都是復(fù)數(shù),我們只好讓L=1時(shí)候的*輸入虛部為0* x_ik+b*tw2是兩個(gè)虛數(shù)相乘*/fprintf(fp,tw1=%f, tw2=%f/n, tw1, tw2);x_rk= TR + x_rk+b*tw1 + x_ik+b*tw2;x_ik= TI - x_rk+b*tw2 + x_ik+b*tw1;x_rk+b = TR - x_rk+b*tw1 - x_ik+b*tw2;x_ik+b = TI + temp*tw2- x_ik+b*tw1;fprintf(fp,k=%d, x_rk=%f, x_ik=%f/n, k, x_rk, x_ik);fprintf(fp,k=%d, x_rk=%f, x_ik=%f/n, k+b, x_rk+b, x_ik+b);/* - add by sshc625 -*該實(shí)現(xiàn)的流程為* for( Layer )*for( Granule )*for( Sample )*/voidfft2(void)fp = fopen(log2.txt,a+);intcur_layer, gr_num, i, k, p;floattmp_real, tmp_imag, temp;/臨時(shí)變量,記錄實(shí)部floattw1, tw2;/旋轉(zhuǎn)因子,tw1為旋轉(zhuǎn)因子的實(shí)部cos部分, tw2為旋轉(zhuǎn)因子的虛部sin部分.intstep;/步進(jìn)intsample_num;/顆粒的樣本總數(shù)(各層不同,因?yàn)楦鲗宇w粒的輸入不同)/*對(duì)層循環(huán)*/for(cur_layer=1; cur_layer 0)i-;gr_num *= 2;/*每個(gè)顆粒的輸入樣本數(shù)N */sample_num= (int)pow(2, cur_layer);/*步進(jìn).步進(jìn)是N/2 */step= sample_num/2;/*/k = 0;/*對(duì)顆粒進(jìn)行循環(huán)*/for(i=0; igr_num; i+)/*對(duì)樣本點(diǎn)進(jìn)行循環(huán),注意上限和步進(jìn)*/for(p=0; psample_num/2; p+)/旋轉(zhuǎn)因子,需要優(yōu)化.tw1 = cos(2*PI*p/pow(2, cur_layer);tw2 = -sin(2*PI*p/pow(2, cur_layer);tmp_real = x_rk+p;tmp_imag = x_ik+p;temp = x_rk+p+step;/*(tw1+jtw2)(x_rk+jx_ik)* real : tw1*x_rk - tw2*x_ik* imag : tw1*x_ik + tw2*x_rk*我想不抽象出一個(gè)* typedef struct *double real;/實(shí)部*double imag;/虛部* complex;以及針對(duì)complex的操作*來簡化復(fù)數(shù)運(yùn)算是否是因?yàn)樾噬系目紤]!*/*蝶形算法*/x_rk+p= tmp_real + ( tw1*x_rk+p+step - tw2*x_ik+p+step );x_ik+p= tmp_imag + ( tw2*x_rk+p+step + tw1*x_ik+p+step );/* Xk = A(k)+WB(k)* Xk+N/2 = A(k)-WB(k)的性質(zhì)可以優(yōu)化這里*/旋轉(zhuǎn)因子,需要優(yōu)化.tw1 = cos(2*PI*(p+step)/pow(2, cur_layer);tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer);x_rk+p+step= tmp_real + ( tw1*temp - tw2*x_ik+p+step );x_ik+p+step= tmp_imag + ( tw2*temp + tw1*x_ik+p+step );printf(k=%d, x_rk=%f, x_ik=%f/n, k+p, x_rk+p, x_ik+p);printf(k=%d, x_rk=%f, x_ik=%f/n, k+p+step, x_rk+p+step, x_ik+p+step);/*開跳!:) */k += 2*step;/*后記:*究竟是顆粒在外層循環(huán)還是樣本輸入在外層,好象也差不多,復(fù)雜度完全一樣.*但以我資質(zhì)愚鈍花費(fèi)了不少時(shí)間才弄明白這數(shù)十行代碼.*從中我發(fā)現(xiàn)一個(gè)于我非常有幫助的教訓(xùn),很久以前我寫過一部分算法,其中絕大多數(shù)都是遞歸.*將數(shù)據(jù)量減少,減少再減少,用歸納的方式來找出數(shù)據(jù)量加大代碼的規(guī)律*比如FFT* 1.先寫死LayerI的代碼;然后再把LayerI的輸出作為LayerII的輸入,又寫死代碼; .*大約3層就可以統(tǒng)計(jì)出規(guī)律來.這和遞歸也是一樣,先寫死一兩層,自然就出來了!* 2.有的功能可以寫偽代碼,不急于求出結(jié)果,降低復(fù)雜性,把邏輯結(jié)果定出來后再添加.*比如旋轉(zhuǎn)因子就可以寫死,就寫1.0.流程出來后再寫旋轉(zhuǎn)因子.*寥寥數(shù)語,我可真是流了不少汗! Happy!*/voiddft(void)inti, n, k, tx1, tx2;floattw1,tw2;floatxx_rN,xx_iN;/* clear any data in Real and Imaginary result arrays prior to DFT*/for(k=0; k=N-1; k+)xx_rk = xx_ik = x_ik = 0.0;/ caculate the DFTfor(k=0; k=(N-1); k+)for(n=0; n= (N/2)tw1 = -twiddletx1-(N/2);elsetw1 = twiddletx1;if(tx2 = (N/2)tw2 = -twiddletx2-(N/2);elsetw2 = twiddletx2;xx_rk = xx_rk+x_rn*tw1;xx_ik = xx_ik+x_rn*tw2;xx_ik = -xx_ik;/ displayfor(i=0; iN; i+)printf(%f/t%f/n, xx_ri, xx_ii);/ -intmain(void)fft_init( );bitrev( );/ bitrev2( );/fft1( );fft2( );display( );system(pause);/ dft();return1;#include #include /* 快速福利葉變換C函數(shù)函數(shù)簡介：此函數(shù)是通用的快速傅里葉變換C語言函數(shù)，移植性強(qiáng)，以下部分不依 賴硬件。此函數(shù)采用聯(lián)合體的形式表示一個(gè)復(fù)數(shù)，輸入為自然順序的復(fù) 數(shù)（輸入實(shí)數(shù)是可令復(fù)數(shù)虛部為0），輸出為經(jīng)過FFT變換的自然順序的 復(fù)數(shù)使用說明：使用此函數(shù)只需更改宏定義FFT_N的值即可實(shí)現(xiàn)點(diǎn)數(shù)的改變，F(xiàn)FT_N的 應(yīng)該為2的N次方，不滿足此條件時(shí)應(yīng)在后面補(bǔ)0函數(shù)調(diào)用：FFT(s);時(shí) 間：2010-2-20版 本：Ver1.0參考文獻(xiàn)： */#include#define PI 3.1415926535897932384626433832795028841971 /定義圓周率值#define FFT_N 128 /定義福利葉變換的點(diǎn)數(shù)struct compx float real,imag; /定義一個(gè)復(fù)數(shù)結(jié)構(gòu)struct compx sFFT_N; /FFT輸入和輸出：從S1開始存放，根據(jù)大小自己定義/*函數(shù)原型：struct compx EE(struct compx b1,struct compx b2) 函數(shù)功能：對(duì)兩個(gè)復(fù)數(shù)進(jìn)行乘法運(yùn)算輸入?yún)?shù)：兩個(gè)以聯(lián)合體定義的復(fù)數(shù)a,b輸出參數(shù)：a和b的乘積，以聯(lián)合體的形式輸出*/struct compx EE(struct compx a,struct compx b) struct compx c; c.real=a.real*b.real-a.imag*b.imag; c.imag=a.real*b.imag+a.imag*b.real; return(c);/*函數(shù)原型：void FFT(struct compx *xin,int N)函數(shù)功能：對(duì)輸入的復(fù)數(shù)組進(jìn)行快速傅里葉變