課程設(shè)計報告FFT的DSP實現(xiàn)

1、DSP原理及應(yīng)用課程設(shè)計報告FFT 的 DSP實現(xiàn)一、設(shè)計目的1、加深對 DFT算法原理和基本性質(zhì)的理解；2、了解并學(xué)習(xí)使用 FFT算法，以及其在 TMS320C54X上的運用；3、學(xué)習(xí) DSP中 FFT的設(shè)計和編程思想；4、練習(xí)使用 CCS 的探針和圖形工具來觀察器觀察波形和頻譜情況。二、設(shè)計內(nèi)容用 C 語言及匯編語言進行編程，實現(xiàn) FFT 運算，對于 C 語言，實現(xiàn) 8 點和 16 點的 FFT運算，對于匯編語言，需調(diào)試出 8點的 FFT運算結(jié)果。三、設(shè)計原理快速傅里葉變換 （FFT）是一種高效實現(xiàn)離散傅里葉變換 （DFT）的快速算法， 是數(shù)字信號處理中最為重要的工具之一， 它在聲學(xué)， 語

2、音，電信和信號處理等領(lǐng) 域有著廣泛的應(yīng)用。1、離散傅里葉變換 DFT對于長度為 N的有限長序列 x（n） ，它的離散傅里葉變換（ DFT）為 n1nkX（k）x（n）WNnk ,k 0,1, N 1（1）n0j 2 / N式中， WN e ，稱為旋轉(zhuǎn)因子或蝶形因子。從 DFT的定義可以看出， 在 x（n） 為復(fù)數(shù)序列的情況下， 對某個 k 值，直 接按（ 1）式計算 X（k） 只需要 N次復(fù)數(shù)乘法和（ N-1）次復(fù)數(shù)加法。因此，對所 有 N個 k 值，共需要 N2次復(fù)數(shù)乘法和 N（N-1） 次復(fù)數(shù)加法。對于一些相當大有 N 值（如 1024 點）來說，直接計算它的 DFT所需要的計算量是很大的

3、，因此 DFT 運算的應(yīng)用受到了很大的限制。2、快速傅里葉變換 FFT 旋轉(zhuǎn)因子 WN 有如下的特性 : 對稱性： WNkWNk N /2周期性： WNk WNk N 利用這些特性，既可以使 DFT中有些項合并，減少了乘法積項，又可以將長 序列的 DFT分解成幾個短序列的 DFT。FFT 就是利用了旋轉(zhuǎn)因子的對稱性和周期 性來減少運算量的。FFT的算法是將長序列的 DFT分解成短序列的 DFT。例如：N為偶數(shù)時， 先將 N點的 DFT分解為兩個 N/2 點的 DFT，使復(fù)數(shù)乘法減少一半：再將每個 N/2 點的 DFT分解成 N/4 點的 DFT，使復(fù)數(shù)乘又減少一半，繼續(xù)進行分解可以大大減少計

4、算量。最小變換的點數(shù)稱為基數(shù)，對于基數(shù)為 2的 FFT算法，它的最小變換是 2 點 DFT。一般而言， FFT算法分為按時間抽取的 FFT( DIT FFT)和按頻率抽取的 FFT (DIF FFT)兩大類。 DIF FFT算法是在時域內(nèi)將每一級輸入序列依次按奇 /偶分 成 2 個短序列進行計算，而 DIF FFT 算法是在頻域內(nèi)將每一級輸入序列依次奇 / 偶分成 2個短序列進行計算。 兩者的區(qū)別是旋轉(zhuǎn)因子出現(xiàn)的位置不同， 但算法是 一樣的。在 DIF FFT算法中，旋轉(zhuǎn)因子 W Nk 出現(xiàn)在輸入端，而在 DIF FFT 算法中 它出現(xiàn)在輸入端。假定序列 x(n) 的點數(shù) N是2的冪，按照 D

5、IF FFT 算法可將其分為偶序列和 奇序列，記偶序列為 x(0), x(2), x(4), x(N - 2),即x1 x(2r),r 0,1, N/2 1 ，奇序 列為 x(1),x(3),x(5), x(N-1),即x2 x(2r 1),r 0,1, N/2 1，則 x(n) 的 FFT表示為 N 1N 1X (k)x(n)W Nnkx(n)WNnkn 0n 0n 為偶數(shù)n為奇數(shù)N / 2 1N / 211)W N( 2r 1)k2 rkx(2r )W N2rkx(2rr0r0N / 2 1N /2 1(r )WN2rk2rkx1 (r)WN2rkW Nkx2(2)r0r0由于WN2e j

6、(2 /N) 2 e j2 /(N/2)WN/2 ，則(3)式可表示為N /2 1N / 2 1X (k)x1 ( r )WNrk/ 2 WNkx2 (r)WNrk/2r 0 r 0X1(k) WNkX2(k)k 0,1, N /2 1 (3)式中， X1(k)和X2(k)分別為 x1(n)和x2 (n)的N/2 的DFT。由于對稱性，WNk N/2WNK , 則 X (k N /2) X1(k) WNkX 2(k) 。因此， N 點X(k) 可分為兩部分：前半部分： X(k) X1(k) WNkX2(k)k 0,1, N /2 1(4)后半部分： X(k N/2) X1(k) WNkX2(k

7、)k 0,1, N /2 1(5)從式(4)和式(5)可以看出，只要求出 0N/2-1 區(qū)間 X1(k)和 X2(k)的值， 就可求出 0N-1區(qū)間X(k)的 N點值。以同樣的方式進行抽取，可以求得 N/4 點的 DFT，重復(fù)抽取過程，就可以使 N點的 DFT用上組 2點的 DFT來計算，這樣就可以大減少運算量。在基數(shù)為 2 的FFT中，設(shè) N=2M，共有 M級運算，每級有 N/2 個2點 FFT蝶形 運算，因此， N點 FFT總共有 MN/2個蝶形運算。蝶形運算如圖 1 所示。圖 1 蝶形運算設(shè)蝶形輸入為 xm1(p)和xm 1(q) ，輸出為 xm(p)和xm(q)，則有k xm(p) x

8、m 1(p) xm 1(q)WN(6)k xm (q) xm 1(p) xm 1(q)WN (7) 在基數(shù)為 2 的FFT中，設(shè) N=2M，共有 M級運算，每級有 N/2 個2點 FFT蝶形 運算，因此， N點 FFT總共有 (N/2)log2 N個蝶形運算。例如：基數(shù)為 2的FFT，當 N=8時，共需要 3級，12個基 2 DIT FFT的蝶形圖 28 點基 2 DIF FFT蝶形運算從圖 2可以看出，輸入是經(jīng)過比特反轉(zhuǎn)的倒位序列， 稱為位碼倒置， 其排列 順序為 x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7) 。輸出是按自然順序排列，其順序 為 x(0),x

9、(1), ,x(6),x(7)。3、FFT運算的實現(xiàn)(1) 實現(xiàn)輸入數(shù)據(jù)的比特反轉(zhuǎn)輸入數(shù)據(jù)的比特反轉(zhuǎn)實際上就是將輸入數(shù)據(jù)進行碼位倒置， 以便在整個運算 后輸出序列是一個自然序列。 在用匯編指令進行碼位倒置時， 使用位碼倒置尋址 可以大大提高程序執(zhí)行速度和使用存儲器的效率。在這種尋址方式下，AR0存放的整數(shù) N的 FFT點的一半，一個輔助寄存器指向一個數(shù)據(jù)存放單元。 當使用位碼 倒置尋址將 AR0加到輔助寄存器時，地址將以位碼倒置的方式產(chǎn)生。(2) 實現(xiàn) N 點復(fù)數(shù) FFTN點復(fù)數(shù) FFT算法的實現(xiàn)可分為三個功能塊，即第一級蝶形運算、第二級蝶 形運算、第三級至 log2N 級蝶形運算。隊以任何一

10、個 2 的整數(shù)冪 N=2M，總可以 通過 M 次分解后最后成為二點的 DFT運算。通過這樣的 M次分解，可以構(gòu)成 M級迭代計算，每級由 N/2 個蝶形運算完成。 （4）輸出 FFT 結(jié)果 FFT算法的程序流程圖如下圖所示。圖 3 FFT 運算開的始程序流程圖四、基于 C 語言的 FFT算法1、源程序 FFT.c #include myapp.h #include ICETEK-VC5509-EDU.h #include scancode.h #include #define PI 3.1415926 #define SAMPLENUMBER 16 void MakeWave(); int IN

11、PUTSAMPLENUMBER; struct compxfloat rea struct compx EE(struc struct compx xinSAM struct compx fWaveSAMPLENUMBER; float ySAMPLENUMBER,dataSAMPL struct compx EE(struc l,imag;t compx,struct com PLENUMBER;設(shè)定輸入序列循環(huán) mm=1 到N3 級蝶形運算求出蝶形運算級數(shù) m=3struct compx b3;BER;Yt compx bY px求）; 該級旋 轉(zhuǎn)因 子下標 Nm循環(huán) 該組 1 到3-mm

12、個蝶形運算2mm-1組蝶形運算2real-b1.imag*b2.imag;Yb3.real=b1.real*b2b2.imag+b1.imag*b2.real;b3.imag=b1.realreturn(b3);main()int i;MakeWave();for ( i=0;iSAMP LENUMBER;i+ )fWavei.rea fWavei.imag=0.0f; yi=0.0f;l=INPUTi;FFT(fWave);for ( i=0;iSAMPLENUMBER;i+ ) datai=yi;while ( 1 ); / break point計算一個蝶形運算單元序列倒序后繪圖結(jié)束vo

13、id FFT(struct compx *xin)int f,m,nv2,nm1,i,k,j=0,l;struct compx v,w,t;nv2=SAMPLENUMBER/2;f=SAMPLENUMBER;for(m=1;(f=f/2)!=1;m+);nm1=SAMPLENUMBER-1;for(i=0;inm1;i+)if(ij)t=xinj;xinj=xini;xini=t;k=nv2;while(k=j)j=j-k;k=k/2;j=j+k;int le,lei,ip;for(l=1;l=m;l+)le=2(l-1);lei=le/2;v.real=1.0;v.imag =0.0;w.r

14、eal =cos(PI/lei);w.imag =-sin(PI/lei); for(j=0;j=lei-1;j+)for(i=j;i=SAMPLENUMBER;i=i+le)ip=i+lei;t=EE(xinip,v);xinip.real=xini.real-t.real; xinip.imag =xini.imag-t.imag; xini.real=xini.real+t.real; xini.imag=xini.imag+t.imag;v=EE(v,w); for ( i=0;iSAMPLENUMBER;i+ ) yi=sqrt(xini.real*xini.real+xini.im

15、ag*xini.imag); void MakeWave() int i;for ( i=0;i DARAM.vectors: VECT .trcinit: DARAM .gblinit: DARAM frt: DARAM .cinit: DARAM .pinit: DARAM .sysinit: DARAM .bss: DARAM2.far: DARAM2.const: DARAM2 .switch: DARAM2.sysmem: DARAM2.cio: DARAM2.MEM$obj: DARAM2.sysheap: DARAM2.sysstack DARAM2.stack: DARAM2

16、3、FFT程序的使用方法（1）根據(jù) N 值，修改 FFT.c 中的中的常數(shù)， 如 N=8，將#define SAMPLENUMBER 16語句中的“ 16”修改為 8。（ 2）編譯、匯編、鏈接，得到 .out 文件，加載。（3）將 data 加入觀察窗，可看到 FFT 運算輸出結(jié)果。4、運行結(jié)果8點的 FFT運算，且輸入為 1、2、3、 8 時，運算結(jié)果如圖 4所示， 16 點的 FFT運算，且輸入為 1、2、3、 16時，運算結(jié)果如圖 5 所示。圖 48 點 FFT 運算結(jié)果圖 5 16 點 FFT 運算結(jié)果五、基于匯編語言的 FFT算法1、匯編源程序 fft.asm.titlefft.as

17、m.mmregs.include .include .defcoeff.incin.incstartsine:.usectsine,512sine1:.usectsine1,512cosine:.usectcosine,512cosine1:.usectcosine1,512fft_data:.usectfft_data,1024fft_out:.usectfft_out,512STACK.usectSTACK,10K_DATA_IDX_1.set 2K_DATA_IDX_2.set 4K_DATA_IDX_3.set 8K_TWID_TBL_SIZE .set 512K_TWID_IDX_3

18、.set 128K_FLY_COUNT_3.set 4K_FFT_SIZE.set 8K_LOGN.set 3PA0.set 0.bssd_twid_idx,1.bssd_data_idx,1.bssd_grps_cnt,1.sectfft_prg位碼倒置程序 *.asgAR2,REORDERED.asgAR3,ORIGINAL_INPUT.asgAR7,DATA_PROC_BUFstart:SSBXFRCTSTM#STACK+10,SPSTM#sine,AR1RPT#K_TWID_TBL_SIZE-1MVPDSTMRPT#sine1,*AR1+cosine,AR1#K_TWID_TBL_SI

19、ZE-1STM#d_input,ORIGINAL_INPUTSTM#fft_data,DA TA_PROC_BUFMVMMDATA_PROC_BUF,REORDEREDSTM#K_FFT_SIZE-1,BRCRPTBDbit_rev_end-1STM#K_FFT_SIZE,AR0MVDD*ORIGINAL_INPUT+,*REORDERED+MVDD*ORIGINAL_INPUT-,*REORDERED+MVPD#cosine1,*AR1+MAR*ORIGINAL_INPUT+0Bbit_rev_end:* FFT CODE *.asgAR1,GROUP_COUNTER.asgAR2,PX.a

20、sgAR3,QX.asgAR4,WR.asgAR5,WI.asgAR6,BUTTERFL Y_COUNTER.asgAR7,STAGE_COUNTER第一級蝶形運算 stage1* STMLDSTMSTMSTMLDRPTBDSTMSUBADDSTHST|LDSUBADDSTHST |LD#0,BK#-1,ASM#fft_data,PX#fft_data+K_DA TA_IDX_1,QXK_FFT_SIZE/2-1,BRC*PX,16,Astage1end-1#K_DATA_IDX_1+1,AR0*QX,16,A,B*QX,16,AA,ASM,*PX+B,*QX+*PX,A*QX,16,A,B*

21、QX,16,AA,ASM,*PX+0%B,*QX+0%*PX,Astage1end:第二級蝶形運算 stage2* STMSTMSTMLDRPTBDSTM#fft_data,PX#fft_data+K_DA TA_IDX_2,QX#K_FFT_SIZE/4-1,BRC*PX,16,Astage2end-1#K_DATA_IDX_2+1,AR0;1st butterflySUB*QX,16,A,BADDSTH*QX,16,AA,ASM,*PX+STB,*QX+ |LD SUB*PX,A*QX,16,A,BADDSTHSTH*QX,16,AA,ASM,*PX+B,ASM,*QX+;2nd butt

22、erflyMAR*QX+ADD*PX,*QX,ASUB*PX,*QX-,BSTHA,ASM,*PX+SUB*PX,*QX,ASTB,*QX|LD*QX+,BSTA,*PX|ADD*PX+0%,ASTA,*QX+0%|LD*PX,Astage2end:第三級至最后一級蝶形運算 *STM#K_TWID_TBL_SIZE,BKST#K_TWID_IDX_3,d_twid_idxSTM#K_TWID_IDX_3,AR0STM#cosine,WRSTM#sine,WISTM#K_LOGN-2-1,STAGE_COUNTERST#K_FFT_SIZE/8-1,d_grps_cntSTM#K_FL Y_CO

23、UNT_3-1,BUTTERFLST#K_DATA_IDX_3,d_data_idxstage:STM#fft_data,PXLDd_data_idx,AADD*(PX),ASTLMA,QXMVDKd_grps_cnt,GROUP_COUNTERgroup:MVMDBUTTERFL Y_COUNTER,BRCRPTBDbutterflyend-1LD*WR,TMPY*QX+,AMAC*WI+0%,*QX-,AADD*PX,16,A,BSTB,*PX|SUB*PX+,BSTB,*QX|MPY*QX+,AMAS*QX,*WR+0%,ACOUNTERADD ST |SUB LD ST |MPY bu

24、tterflyend:PSHM MVDK MAR MAR BANZD POPM MAR LD SUB STLM STL LD STL LD STL BANZD MVDK*PX,16,A,BB,*QX+*PX,B*WR,TB,*PX+*QX+,AAR0d_data_idx,AR0*PX+0*QX+0group,*GROUP_COUNTER-AR0*QX- d_data_idx,A#1,A,BB,BUTTERFL Y_COUNTERA,1,d_data_idx d_grps_cnt,A A,ASM,d_grps_cnt d_twid_idx,A A,ASM,d_twid_idx stage,*STAGE_COUNTER- d_twid_idx,AR0fft_end:STM; STMSTMSTMRPTBSQURSQURASTHpow