基于DSP的FFT實(shí)現(xiàn)設(shè)計(jì)報(bào)告丁濤

1、基于DSP實(shí)驗(yàn)平臺(tái)的FFT算法實(shí)現(xiàn) 組員： 丁 濤, 賈 冰 , 葛云峰, 程自奇, 日期：2015.01.15一、實(shí)驗(yàn)?zāi)康?. 加深對(duì)DFT算法原理和基本性質(zhì)的理解；2. 熟悉FFT的算法原理和FFT子程序的算法流程和應(yīng)用；3. 學(xué)習(xí)用FFT對(duì)連續(xù)信號(hào)和時(shí)域信號(hào)進(jìn)行頻譜分析的方法；4. 學(xué)習(xí)DSP中FFT的設(shè)計(jì)和編程思想；5. 學(xué)習(xí)使用CCS的波形觀察器觀察波形和頻譜情況；二、實(shí)驗(yàn)內(nèi)容 用DSP匯編語(yǔ)言及C語(yǔ)言進(jìn)行編程，實(shí)現(xiàn)FFT運(yùn)算、對(duì)輸入信號(hào)進(jìn)行頻譜分析。三、實(shí)驗(yàn)原理快速傅里葉變換FFT旋轉(zhuǎn)因子WN 有如下的特性。對(duì)稱性：WNk+N/2=-WNk （2）周期性：WNn(N-k)=WNk(

2、N-n)=WN-nk （3）利用這些特性，既可以使DFT中有些項(xiàng)合并，減少了乘法積項(xiàng)，又可以將長(zhǎng)序列的DFT分解成幾個(gè)短序列的DFT。FFT就是利用了旋轉(zhuǎn)因子的對(duì)稱性和周期性來(lái)減少運(yùn)算量的。FFT的算法是將長(zhǎng)序列的DFT分解成短序列的DFT。例如：N為偶數(shù)時(shí)，先將N點(diǎn)的DFT分解為兩個(gè)N/2點(diǎn)的DFT，使復(fù)數(shù)乘法減少一半：再將每個(gè)N/2點(diǎn)的DFT分解成N/4點(diǎn)的DFT，使復(fù)數(shù)乘又減少一半，繼續(xù)進(jìn)行分解可以大大減少計(jì)算量。最小變換的點(diǎn)數(shù)稱為基數(shù)，對(duì)于基數(shù)為2的FFT算法，它的最小變換是2點(diǎn)DFT。一般而言，F(xiàn)FT算法分為按時(shí)間抽取的FFT（DIT FFT）和按頻率抽取的FFT（DIF FFT）兩

3、大類。DIF FFT算法是在時(shí)域內(nèi)將每一級(jí)輸入序列依次按奇偶分成2個(gè)短序列進(jìn)行計(jì)算。而DIF FFT算法是在頻域內(nèi)將每一級(jí)輸入序列依次奇偶分成2個(gè)短序列進(jìn)行計(jì)算。兩者的區(qū)別是旋轉(zhuǎn)因子出現(xiàn)的位置不同，得算法是一樣的。在DIF FFT算法中，旋轉(zhuǎn)因子WN出現(xiàn)在輸入端，而在DIF FFT算法中它出現(xiàn)在輸入端。假定序列x(n)的點(diǎ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。 由于對(duì)稱性，

4、WNk+N/2=-WNk。因此，N點(diǎn)DFT可分為兩部分：前半部分：x(k)=x1(k)+WkNx2(k) （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點(diǎn)值。以同樣的方式進(jìn)行抽取，可以求得N/4點(diǎn)的DFT，重復(fù)抽取過(guò)程，就可以使N點(diǎn)的DFT用上組2點(diǎn)的DFT來(lái)計(jì)算，這樣就可以大減少運(yùn)算量。基2 DIF FFT的蝶形運(yùn)算如圖3.1所示。設(shè)蝶形輸入為x1(k)和x2(k)，輸出為x(k)和x(N/2+K)，則有 x(k)=x1(k)+Wk

5、Nx2(k) （6）x(N/2+k)=x1(k)-WkNx2(k) （7）在基數(shù)為2的FFT中，設(shè)N=2M，共有M級(jí)運(yùn)算，每級(jí)有N/2個(gè)2點(diǎn)FFT蝶形運(yùn)算，因此，N點(diǎn)FFT總共有MN/2個(gè)蝶形運(yùn)算。 圖3.1 基2 DIF FFT的蝶形運(yùn)算例如：基數(shù)為2的FFT，當(dāng)N=8時(shí)，共需要3級(jí)，12個(gè)基2 DIT FFT的蝶形運(yùn)算。其信號(hào)流程如圖3.2所示。x(0) x(0) 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 WN2x(3) x(6) -1

6、-1 WN0 WN2 WN3x(7) x(7) -1 -1 -1圖3.2 8點(diǎn)基2 DIF FFT蝶形運(yùn)算從圖(b)可以看出，輸入是經(jīng)過(guò)比特反轉(zhuǎn)的倒位序列，稱為位碼倒置，其排列順序?yàn)閤(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)，輸出是按自然順序排列，其順序?yàn)閤(0),x(1),x(2),x(3),x(4),x(5),x(6),x(7).程序設(shè)計(jì)順序DSP初始化串口設(shè)置AD設(shè)置功率譜計(jì)算觀看轉(zhuǎn)換結(jié)果，保存數(shù)據(jù)串口發(fā)送轉(zhuǎn)換結(jié)果FFT運(yùn)算位碼倒置串口接收，AD采樣設(shè)置信號(hào)源類型、頻率幅值、和采樣點(diǎn)數(shù) 四、FFT算法的DSP實(shí)現(xiàn)過(guò)程：DSP芯片的出現(xiàn)使FFT的實(shí)現(xiàn)方法

7、變得更為方便。由于大多數(shù)DSP芯片都具有在單指令周期內(nèi)完成乘法累加操作，并且提供了專門的FFT指令，使得FFT算法在DSP芯片實(shí)現(xiàn)的速度更快。FFT算法可以分為按時(shí)間抽取FFT和按頻率抽取FFT兩大類，輸入也有實(shí)數(shù)和復(fù)數(shù)之分，一般情況下，都假定輸入序列為復(fù)數(shù)。（一）FFT運(yùn)算序列的存儲(chǔ)分配FFT運(yùn)算時(shí)間是衡量DSP芯片性能的一個(gè)重要指標(biāo)，因此提高FFT的運(yùn)算速度是非常重要的。在用DSP芯片實(shí)現(xiàn)FFT算法時(shí)，應(yīng)允許利用DSP芯片所提供的各種軟、硬件資源。如何利用DSP芯片的有限資源，合理地安排好所使用的存儲(chǔ)空間是十分重要的。(二)FFT運(yùn)算的實(shí)現(xiàn)用TMS320C54x的匯編程序?qū)崿F(xiàn)FFT算法主要

8、分為四步：1.實(shí)現(xiàn)輸入數(shù)據(jù)的比特反轉(zhuǎn)輸入數(shù)據(jù)的比特反轉(zhuǎn)實(shí)際上就是將輸入數(shù)據(jù)進(jìn)行碼位倒置，以便在整個(gè)運(yùn)算后的輸出序列是一個(gè)自然序列。在用匯編指令進(jìn)行碼位倒置時(shí)，使用碼位倒置可以大大提高程序執(zhí)行速度和使用存儲(chǔ)器的效率。在這種尋址方式下，AR0存放的整數(shù)N是FFT點(diǎn)的一半，一個(gè)輔助寄存器指向一個(gè)數(shù)據(jù)存放的單元。當(dāng)使用位碼倒置尋址將AR0加到輔助寄存器時(shí)，地址將以位碼倒置的方式產(chǎn)生。2.實(shí)現(xiàn)N點(diǎn)復(fù)數(shù)FFTN點(diǎn)復(fù)數(shù)FFT算法的實(shí)現(xiàn)可分為三個(gè)功能塊，即第一級(jí)蝶形運(yùn)算、第二級(jí)蝶形運(yùn)算、第三級(jí)至log2N級(jí)蝶形運(yùn)算。對(duì)于任何一個(gè)2的整數(shù)冪，總可以通過(guò)M次分解最后成為2點(diǎn)的DFT計(jì)算。通過(guò)這樣的M次分解，可構(gòu)成

9、M（即log2N）級(jí)迭代計(jì)算，每級(jí)由N/2個(gè)蝶形運(yùn)算組成。3.功率譜的計(jì)算用FFT計(jì)算想x(n)的頻譜，即計(jì)算X（k）=X(k)一般是由實(shí)部XR(k)和虛部Xl(k)組成的復(fù)數(shù)，即X（k）=XR(k)+jXl(k)因此，計(jì)算功率譜時(shí)只需將FFT變換好的數(shù)據(jù)，按照實(shí)部實(shí)部XR(k)和虛部Xl(k)求它們的平方和，然后對(duì)平方和進(jìn)行開平方運(yùn)算。但是考慮到編程的難度，對(duì)于求FFT變換后數(shù)據(jù)的最大值，不開平方也可以找到最大值，并對(duì)功率譜的結(jié)果沒有影響，所以在實(shí)際的DSP編程中省去了開方運(yùn)算。4.輸出FFT結(jié)果(三)匯編語(yǔ)言程序程序主體由rfft-task、bit-rev、fft和power四個(gè)子程序組成

10、。rfft-task：主調(diào)用子程序，用來(lái)調(diào)用其他子程序，實(shí)現(xiàn)統(tǒng)一的接口。bit-rev：位碼倒置子程序，用來(lái)實(shí)現(xiàn)輸入數(shù)據(jù)的比特反轉(zhuǎn)。fft：FFT算法子程序，用來(lái)完成N點(diǎn)FFT運(yùn)算。在運(yùn)算過(guò)程中，為避免運(yùn)算結(jié)果的溢出，對(duì)每個(gè)蝶形的運(yùn)算結(jié)果右移一位。fft子程序分為三個(gè)功能塊：第一級(jí)蝶形運(yùn)算、第二級(jí)蝶形運(yùn)算、第三級(jí)至至log2N級(jí)蝶形運(yùn)算。(四）正弦系數(shù)表和余弦系數(shù)表：正弦系數(shù)表和余弦系數(shù)表可以由數(shù)據(jù)文件coeff.inc給出，主程序通過(guò).copy匯編命令將正弦和余弦系數(shù)表與程序代碼匯編在一起。在本例中，數(shù)據(jù)文件coeff.inc給出1024復(fù)數(shù)點(diǎn)FFT的正弦、余弦系數(shù)各512個(gè)。利用此系數(shù)表可

11、完成81024點(diǎn)FFT的運(yùn)算。（五）FFT算法的模擬信號(hào)輸入：FFT算法的模擬信號(hào)輸入可以采用C語(yǔ)言編程來(lái)生成一個(gè)文本文件sindata，然后在rfft-task匯編程序中，通過(guò).copy匯編命令將生成的數(shù)據(jù)文件復(fù)制到數(shù)據(jù)存儲(chǔ)器中，作為FFT算法的輸入數(shù)據(jù)參與FFT運(yùn)算。這種方法的優(yōu)點(diǎn)是程序的可讀性強(qiáng)，缺點(diǎn)是當(dāng)輸入數(shù)據(jù)修改后，必須重新編譯、匯編和鏈接。五、設(shè)計(jì)步驟：1.啟動(dòng)CCS，在CCS中建立一個(gè)C源文件和一個(gè)命令文件，并將這兩個(gè)文件添加到工程，再編譯并裝載程序：閱讀Dsp原理及應(yīng)用中fft 用dsp實(shí)現(xiàn)的有關(guān)程序。2.雙擊，啟動(dòng)CCS的仿真平臺(tái)的配著選項(xiàng)。選擇C5502 Simulator

12、。3.啟動(dòng)ccs2后建立工程文件FFT.pjt4.建立源文件FFT.c與鏈接文件FFT.cmd5.將這兩個(gè)文件加到FFT.pjt這個(gè)工程中。6.創(chuàng)建out文件7.加載out文件 六、實(shí)驗(yàn)過(guò)程1編譯并下載程序。2打開觀察窗口：*選擇菜單View-Graph-Time/Frequency進(jìn)行如下圖所示設(shè)置。圖1圖2圖33清除顯示：在以上打開的窗口中單擊鼠標(biāo)右鍵，選擇彈出式菜單中“Clear Display”功能。4設(shè)置斷點(diǎn)：在程序FFT.c 中有注釋“break point”的語(yǔ)句上設(shè)置軟件斷點(diǎn)。圖45運(yùn)行并觀察結(jié)果。 選擇“Debug”菜單的“Animate”項(xiàng)，或按Alt+F5 鍵運(yùn)行程序。

13、觀察“Test Wave”窗口中時(shí)域圖形；圖5 在“Test Wave”窗口中點(diǎn)擊右鍵，選擇屬性，更改圖形顯示為FFT。觀察頻域圖形。圖6 觀察“FFT”窗口中的由CCS 計(jì)算出的正弦波的FFT。圖7七、實(shí)驗(yàn)結(jié)果通過(guò)觀察頻域和時(shí)域圖，程序計(jì)算出了測(cè)試波形的功率譜，與CCS 計(jì)算的FFT 結(jié)果相近。八、心得體會(huì)本來(lái)打算想自己弄的,但一方面時(shí)間有限,另一方面確實(shí)能力有限,對(duì)硬件不太懂,所以只能做到這種程度了。附 錄Cmd源文件代碼：-f 0-w-stack 500-sysstack 500-1 rts55.libMEMORY DARAN: o=0x100, 1=0x7f00 VECT:o=0x80

14、00, 1=0x100 DARAM2:o=0x8100,1=0x7f00 SARAM:o=0x10000,1=0x30000 SDRAM:o=0x40000,1=0x3e0000 SECTIONS .text:DARAM .vectors:VECT .trcinit:DARAM .gblinit:DARAM .frt:DARAM .cinit:DARAM .pinit:DARAM .sysinit:DARAM2 .far:DARAM2 .const:DARAM2 .switch:DARAM2 .sysmem:DARAM2 .cio:DARAM2 .MEM$obj:DARAM2 .sysheap

15、:DARAM2 .sysstack:DARAM2 .stack:DARAM2 .input:DARAM2 .fftcode:DARAM2C文件源碼：#include math.h#define sample_1 256#define signal_1_f 60#define signal_2_f 200#define signal_sample_f 512#define pi 3.int inputsample_1;float fwaversample_1,fwaveisample_1,wsample_1;float sin_tabsample_1;float cos_tabsample_1;

16、void init_fft_tab();void fft(float datarsample_1,float dataisample_1);void main() int i; init_fft_tab(); input_data(); for(i=0;isample_1;i+) fwaveri=inputi; fwaveii=0.0f; wi=0.0f; fft(fwaver,fwavei); while(1);void init_fft_tab() float wt1; float wt2; int i; for(i=0;isample_1;i+) wt1=2*pi*i*signal_1_

17、f; wt1=wt1/signal_sample_f; wt2=2*pi*i*signal_2_f; wt2=wt2/signal_sample_f; inputi=(cos(wt1)+cos(wt2)/2*32768; void input_data() int i; for(i=0;isample_1;i+) sin_tabi=sin(2*pi*i/sample_1); cos_tabi=cos(2*pi*i/sample_1); void fft(float datarsample_1,float dataisample_1) int x0,x1,x2,x3,x4,x5,x6,x7,xx

18、; int i,j,k,b,p,L; float TR,TI,temp; for(i=0;isample_1;i+) 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;x7=(i/128)&0x01; xx=x0*128+x1*64+x2*32+x3*16+x4*8+x5*4+x6*2+x7; dataixx=datari; for(i=0;isample_1;i+) datari=dataii;dataii=0; for(L=1;L0) b=b*2;i-; for(j=0;j0) p=p*2;i-; p=p*j; for(k=j;k256;k=k+2*b) TR=da