DSP課程設計——FFT的DSP實現(xiàn)

1、 FFT的DSP實現(xiàn)簡介：快速傅里葉變換是一種高效實現(xiàn)離散傅里葉變換的的快速算法 ，是數(shù)字信號處理中最為重要的工具之一，它在聲學、語音、電信和信號處理等領域有著廣泛的應用。一設計目的:1.加深對DFT算法原理和基本性質的理解；2.熟悉FFT的算法原理和FFT子程序的算法流程和應用；3.學習用FFT對連續(xù)信號和時域信號進行頻譜分析的方法；4.學習DSP中FFT的設計和編程思想；5.學習使用CCS的波形觀察窗口觀察信號波形和頻譜情況。二設計內容:用DSP匯編語言及C語言進行編程，實現(xiàn)FFT運算，對輸入信號進行頻譜分析。三設計原理:1 離散傅里葉變換DFT:對于長度為N的有限長序列x(n)，它的離散

2、傅里葉變換（DFT）為X(k)= N-nk ,k=0,1,2N-1 （1）式中，WN=e-j*2/N ，稱為旋轉因子或蝶形因子。從DFT的定義可以看出，在x(n)為復數(shù)序列的情況下，對某個k值，直接按（1）式計算X(k) 只需要N次復數(shù)乘法和（N-1）次復數(shù)加法。因此，對所有N個k值，共需要N2次復數(shù)乘法和N(N-1)次復數(shù)加法。對于一些相當大有N值（如1024點）來說，直接計算它的DFT所需要的計算量是很大的，因此DFT運算的應用受到了很大的限制。2快速傅里葉變換FFT旋轉因子WN 有如下的特性。對稱性： WNk+N/2=-WNk周期性：WNn(N-k)=WNk(N-n)=WN-nk利用這些

3、特性，既可以使DFT中有些項合并，減少了乘法積項，又可以將長序列的DFT分解成幾個短序列的DFT。FFT就是利用了旋轉因子的對稱性和周期性來減少運算量的。FFT的算法是將長序列的DFT分解成短序列的DFT。例如：N為偶數(shù)時，先將N點的DFT分解為兩個N/2點的DFT，使復數(shù)乘法減少一半：再將每個N/2點的DFT分解成N/4點的DFT，使復數(shù)乘又減少一半，繼續(xù)進行分解可以大大減少計算量。最小變換的點數(shù)稱為基數(shù)，對于基數(shù)為2的FFT算法，它的最小變換是2點DFT。一般而言，F(xiàn)FT算法分為按時間抽取的FFT（DIT）和按頻率抽取的（DIF FFT）兩大類。IF FFT算法是在時域內將每一級輸入序列依

4、次按奇偶分成個短序列進行計算。而DIF FFT算法是在頻域內將每一級輸入序列依次奇偶分成個短序列進行計算。兩者的區(qū)別是旋轉因子出現(xiàn)的位置不同，得算法是一樣的。在DIF FFT算法中，旋轉因子出現(xiàn)在輸入端，而在DIF FFT算法中它出現(xiàn)在輸入端。假定序列x(n)的點數(shù)N是2的冪，按照DIF FFT算法可將其分為偶序列和奇序列。偶序列：x(2r)=x1(r)奇序列：x(2r+1)=x2(r)其中：r=0,1,2,N/2-1 則x(n)的DFT表示為式中，x1(k)和x2(k)分別為x1(r)和x2(r)的N/2的DFT式中，x1(k)和x2(k)分別為x1(r)和x2(r)的N/2的DFT。由于對

5、稱性，WNk+N/2=-WNk。因此，N點DFT可分為兩部分：前半部分：x(k)=x1(k)+WkNx2(k) k=0,1,N/2-1 （4）后半部分：x(N/2+k)=x1(k)-WkNx2(k) 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，重復抽取過程，就可以使N點的DFT用上組2點的 DFT來計算，這樣就可以大減少運算量?；? DIF FFT的蝶形運算如圖(a)所示。設蝶形輸入為X1(K)和X2(K)，輸出為x(k)和x(N/2+K)，

6、則有x(k)=x1(k)+WkNx2(k) （6）x(N/2+k)=x1(k)-WkNx2(k) （7）在基數(shù)為2的FFT中，設N=2M，共有M級運算，每級有N/2個2點FFT蝶形運算，因此，N點FFT總共有MN/2個蝶形運算。圖(a) 基2 DIF FFT的蝶形運算例如：基數(shù)為2的FFT，當N=8時，共需要3級，12個基2 DIT FFT的蝶形運算。其信號流程如圖(b)所示。x(0) x(0)1. WN0x(4) x(1) -1 WN0x(2) x(2) -1 WN0 WN2x(6) x(3) -1 -1 WN0x(1) x(4) -1 WN0 WN1x(5) x(5) -1 -1 WN0

7、WN2x(3) x(6) -1 -1 WN0 WN2 WN3x(7) x(7) -1 -1 -1 圖(b) 8點基2 DIF FFT蝶形運算從圖(b)可以看出，輸入是經過比特反轉的倒位序列，稱為位碼倒置，其排列順序為x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)。輸出的是按自然順序排列，其順序為x(0),x(1),x(2),x(3),x(4),x(5),x(6),x(7).四FFT算法的DSP實現(xiàn)過程：DSP芯片的出現(xiàn)使FFT的實現(xiàn)方法變得更為方便。由于大多數(shù)DSP芯片都具有在單指令周期內完成乘法累加操作，并且提供了專門的FFT指令，使得FFT算法在DSP芯片實

8、現(xiàn)的速度更快。FFT算法可以分為按時間抽取FFT和按頻率抽取FFT兩大類，輸入也有實數(shù)和復數(shù)之分，一般情況下，都假定輸入序列為復數(shù)。（一）FFT運算序列的存儲分配FFT運算時間是衡量DSP芯片性能的一個重要指標，因此提高FFT的運算速度是非常重要的。在用DSP芯片實現(xiàn)FFT算法時，應允許利用DSP芯片所提供的各種軟、硬件資源。如何利用DSP芯片的有限資源，合理地安排好所使用的存儲空間是十分重要的。(二)FFT運算的實現(xiàn)用TMS320C54x的匯編程序實現(xiàn)FFT算法主要分為四步：1.實現(xiàn)輸入數(shù)據(jù)的比特反轉輸入數(shù)據(jù)的比特反轉實際上就是將輸入數(shù)據(jù)進行碼位倒置，以便在整個運算后的輸出序列是一個自然序列

9、。在用匯編指令進行碼位倒置時，使用碼位倒置可以大大提高程序執(zhí)行速度和使用存儲器的效率。在這種尋址方式下，AR0存放的整數(shù)N是FFT點的一半，一個輔助寄存器指向一個數(shù)據(jù)存放的單元。當使用位碼倒置尋址將AR0加到輔助寄存器時，地址將以位碼倒置的方式產生。2.實現(xiàn)N點復數(shù)FFTN點復數(shù)FFT算法的實現(xiàn)可分為三個功能塊，即第一級蝶形運算、第二級蝶形運算、第三級至級蝶形運算。對于任何一個2的整數(shù)冪，總可以通過M次分解最后成為2點的DFT計算。通過這樣的M次分解，可構成M（即）級迭代計算，每級由N/2個蝶形運算組成。3.功率譜的計算用FFT計算想x(n)的頻譜，即計算X（k）=X(k)一般是由實部(k)和

10、虛部(k)組成的復數(shù)，即X（k）=(k)+j(k)因此，計算功率譜時只需將FFT變換好的數(shù)據(jù)，按照實部實部(k)和虛部(k)求它們的平方和，然后對平方和進行開平方運算。但是考慮到編程的難度，對于求FFT變換后數(shù)據(jù)的最大值，不開平方也可以找到最大值，并對功率譜的結果沒有影響，所以在實際的DSP編程中省去了開方運算。4.輸出FFT結果(三)匯編語言程序程序主體由rfft-task、bit-rev、fft和power四個子程序組成。rfft-task：主調用子程序，用來調用其他子程序，實現(xiàn)統(tǒng)一的接口。bit-rev：位碼倒置子程序，用來實現(xiàn)輸入數(shù)據(jù)的比特反轉。fft：FFT算法子程序，用來完成N點F

11、FT運算。在運算過程中，為避免運算結果的溢出，對每個蝶形的運算結果右移一位。fft子程序分為三個功能塊：第一級蝶形運算、第二級蝶形運算、第三級至至級蝶形運算。(四）正弦系數(shù)表和余弦系數(shù)表：正弦系數(shù)表和余弦系數(shù)表可以由數(shù)據(jù)文件coeff.inc給出，主程序通過.copy匯編命令將正弦和余弦系數(shù)表與程序代碼匯編在一起。在本例中，數(shù)據(jù)文件coeff.inc給出1024復數(shù)點FFT的正弦、余弦系數(shù)各512個。利用此系數(shù)表可完成81024點FFT的運算。（五）FFT算法的模擬信號輸入：FFT算法的模擬信號輸入可以采用C語言編程來生成一個文本文件sindata，然后在rfft-task匯編程序中，通過.c

12、opy匯編命令將生成的數(shù)據(jù)文件復制到數(shù)據(jù)存儲器中，作為FFT算法的輸入數(shù)據(jù)參與FFT運算。這種方法的優(yōu)點是程序的可讀性強，缺點是當輸入數(shù)據(jù)修改后，必須重新編譯、匯編和鏈接。五設計步驟：1.啟動CCS，在CCS中建立一個C源文件和一個命令文件，并將這兩個文件添加到工程，再編譯并裝載程序：閱讀Dsp原理及應用中fft 用dsp實現(xiàn)的有關程序。2.雙擊，啟動CCS的仿真平臺的配著選項。選擇C5502 Simulator。3.啟動ccs2后建立工程文件FFT.pjt4.建立源文件FFT.c與鏈接文件FFT.cmd5.將這兩個文件加到FFT.pjt這個工程中。6.創(chuàng)建out文件7.加載out文件8.加載

13、數(shù)據(jù)9.觀察輸入輸出波形輸入波形（時域）輸出圖形（頻域）10.改變信號的頻率可以再做次實驗。也可作512點或更多點的FFT.六實驗程序：.title "rfft_task.asm".mmregs.copy "coeff.inc".def rfft_tasksine: .usect "sine", 512cosine: .usect "cosine", 512fft_data: .usect "fft_data", 2048d_input: .usect "fft_data",

14、 2048fft_out: .usect "fft_out", 1024;d_input: .copy sindataSTACK: .usect "STACK", 10K_DATA_IDX_1 .set 2K_DATA_IDX_2 .set 4K_DATA_IDX_3 .set 8K_FLY_COUNT_3 .set 4K_TWID_TBL_SIZE .set 512K_TWID_IDX_3 .set 128K_FFT_SIZE .set 32K_LOGN .set 5.bss d_twid_idx, 1.bss d_data_idx, 1.bss d_

15、grps_cnt, 1.sect "rfft_prg"rfft_rask:SSBX FRCTSTM #STACK+10, SPSTM #sine, AR1RPT #K_TWID_TBL_SIZE-1MVPD sine1, *AR1+STM #cosine, AR1RPT #K_TWID_TBL_SIZE-1MVPD cosine1, *AR1+CALL bit_revCALL fftCALL powerRET* 位碼倒置程序 bit_rev *.asg AR2, REORDERED.asg AR3,ORIGINAL_INPUT.asg AR7, DATA_PROC_BUF.

16、sect "rfft_prg"bit_rev:STM #d_input, ORIGINAL_INPUTSTM #fft_data, DATA_PROC_BUFMVMM DATA_PROC_BUF, REORDEREDSTM #K_FFT_SIZE-1, BRCRPTBD bit_rev_end-1STM #K_FFT_SIZE, AR0MVDD *ORIGINAL_INPUT+, *REORDERED+MVDD *ORIGINAL_INPUT-, *REORDERED+MAR *ORIGINAL_INPUT+0B;bit_rev_end:RET;* FFT算法子程序 fft

17、 *;.asg AR1, GROUP_COUNTER.asg AR2, PX.asg AR3, QX.asg AR4, WR.asg AR5, WI.asg AR6, BUTTERFLY_COUNTER.asg AR7, STAGE_COUNTER.sect "rfft_prg"fft:;*第1級蝶形運算 stage1*:STM #0, BKLD #-1, ASMSTM #fft_data, PXLD *PX, 16, ASTM #fft_data+K_DATA_IDX_1,QXSTM #K_FFT_SIZE/2-1, BRCRPTBD stage_end-1STM #K_

18、DATA_IDX_1+1, AR0SUB *QX, 16, A, BADD *QX, 16, ASTH A, ASM, *PX+ST B, *QX+|LD *PX, ASUB *QX, 16, A, BADD *QX, 16, ASTH A, ASM, *PX+0ST B, *QX+0%|LD *PX, A;stage1_end:;*第2級蝶形運算 stage2*;STM #fft_data, PXSTM #fft_data+K_DATA_IDX_2, QXSTM #K_FFT_SIZE/4-1, BRCLD *PX, 16, ARPTBD stage2_end-1STM #K_DATA_ID

19、X_2+1, AR0SUB *QX, 16, A, BADD *QX, 16, ASTH A, ASM, *PX+ST B, *QX+|LD *PX, ASUB *QX, 16, A, BADD *QX, 16, ASTH A, ASM, *PX+STH B, ASM, *QX+MAR *QX+ADD *PX, *QX, ASUB *PX, *QX-, BSTH A, ASM, *PX+SUB *PX, *QX, AST B, *QX|LD *QX+, BST A, *PX|ADD *PX+0%, AST A, *QX+0%|LD *PX, A;stage2_end:;*第3級蝶形運算 sta

20、ge3*;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_FLY_COUNT_3-1, BUTTERFLY_COUNTERST #K_DATA_IDX_3, d_data_idx;stage:STM #fft_data, PXLD d_data_idx, AADD *(PX), ASTLM A, QXMVDK

21、 d_grps_cnt, GROUP_COUNTER;group:MVMD BUTTERFLY_COUNTER, BRCRPTBD butterfly_end-1LD *WR, TMPY *QX+, AMACR *WI+0%, *QX-, AADD *PX, 16, A, BST B, *PX|SUB *PX+, BST B, *QX|MPY *QX+, AMASR *QX, *WR+0%, AADD *PX, 16, A, BST B, *QX+|SUB *PX+, BLD *WR, TST B, *PX+|MPY *QX, A;butterfly_end:PSHM AR0MVDK d_da

22、ta_idx, AR0MAR *PX+0MAR *QX+0BANZD group, *GROUP_COUNTER-POPM AR0MAR *QX-LD d_data_idx, ASUB #1, A, BSTLM B, BUTTERFLY_COUNTERSTL A, 1, d_data_idxLD d_group_cnt, ASTL A, ASM, d_group_cntLD d_twid_idx, ASTL A, ASM, d_twid_idxBANZD stage, *STAGE_COUNTER-MVDK d_twid_idx, AR0;fft_end:RET;*/功率譜計算子程序 power /*;.sect "rfft_prg"power:STM #fft_data, AR2STM #fft_out, AR4STM #K_FFT_SIZE*2-1, BRCRPTB power_end-1SQUR *AR2+, ASQURA *AR2+, ASTH A, *AR4+;power_end:RET.end鏈接命令文件rfft_task.cmd清單：vector.objrfft_task.obj-o rfft_task.obj-m rfft_task.map-e