TI庫實(shí)現(xiàn)2048點(diǎn)或更大點(diǎn)數(shù)FFT..

1、在DSP運(yùn)算中，經(jīng)常需要把輸入時(shí)域信號(hào)在頻域進(jìn)行處理之后，再還原為時(shí)域信號(hào)，這樣就需要進(jìn)行FFT和IFFT運(yùn)算：x(n)-FFT-X(f)-頻域處理-Y(f)-IFFT-y(n)而一般的DSP芯片只支持整數(shù)運(yùn)算，也就是說只能進(jìn)行定點(diǎn)小數(shù)計(jì)算。N點(diǎn)FFT計(jì)算出0N-1，N個(gè)復(fù)數(shù)：0,A，N/2,A*，A為(N/2-1)個(gè)復(fù)數(shù)，A*為A的共軛復(fù)數(shù)。FFT的公式為：NX(k)=sumx(n)*exp(-j*2*pi*(k-1)*(n-1)/N)、1=k=N.n=1IFFT的公式為：Nx(n)=(1/N)sumX(k)*exp(j*2*pi*(k-1)*(n-1)/N)、1=n1024FFT-Q5-1

2、024IFFT-Q5這樣，經(jīng)過FFT和IFFT變換之后，數(shù)據(jù)從Q15變成了Q5，丟失了10bit的信息。為了使得定點(diǎn)FFT和IFFT之后能夠還原為原來的數(shù)據(jù)，必須使用32bit定點(diǎn)FFT，輸入數(shù)據(jù)雖然是16bit，在進(jìn)行FFT的時(shí)候?qū)⑺D(zhuǎn)換為32bit運(yùn)算，這樣輸入的數(shù)據(jù)為Q31,1024點(diǎn)FFT的輸出數(shù)據(jù)為Q21，滿足16bit的精度要求。這樣整個(gè)計(jì)算流程變?yōu)椋篞15-Q31-1024FFT32bit-Q21-1024IFFT32bit-Q21-Q15在TI的DSPLIB中的FFT和IFFT都提供了兩種選擇：SCALE和NOSCALE。仍以1024點(diǎn)計(jì)算為例,它們的意思分別為：Q15-NO

3、SCALEFFT-Q15（結(jié)果可能溢出）Q15-SCALEFFT-Q5（結(jié)果不會(huì)溢出，但是精度降低）Q15-NOSCALEIFFT-Q15（結(jié)果可能溢出）Q15-SCALEIFFT-Q5（結(jié)果不會(huì)溢出，但是精度降低）32位定點(diǎn)FFT,IFFT與此類似，不再重復(fù)。因此，用DSPLIB進(jìn)行FFT和IFFT計(jì)算時(shí)，注意必須采用32bit精度，而且FFT變換時(shí)采用SCALE，而IFFT變換時(shí)采用NOSCALEo關(guān)于TIFFT庫的使用2010-01-2512:53星期一這篇文章是應(yīng)一些找我討論DSP的同學(xué)所寫，貼在這里大家一起學(xué)習(xí)。曾有不少論壇上的同學(xué)（包括DSP算法討論群里的一些同學(xué)）問過我關(guān)于TI的

4、FFT庫的使用，這里我將我使用過的一些經(jīng)驗(yàn)說一下。TI的這個(gè)FFT庫在計(jì)算速度、計(jì)算精度以及數(shù)據(jù)存儲(chǔ)等方面是做了不少優(yōu)化的，比如數(shù)據(jù)存儲(chǔ)，若作N點(diǎn)的FFT,供查表用的旋轉(zhuǎn)因子必須有N/2點(diǎn)的正弦值與N/2點(diǎn)的余弦值，這個(gè)庫將其壓縮成3N/4點(diǎn)的正弦值，因此就節(jié)省了N/4點(diǎn)的存儲(chǔ)空間；另外計(jì)算N點(diǎn)實(shí)數(shù)FFT時(shí)，一般簡(jiǎn)單的做法是將N點(diǎn)實(shí)數(shù)的虛部全化為0來處理，而這個(gè)庫則把N點(diǎn)實(shí)數(shù)數(shù)據(jù)打包成N/2點(diǎn)復(fù)數(shù)數(shù)據(jù)來處理，在計(jì)算速度和存儲(chǔ)空間都有很大改進(jìn)。之前我瀏覽helloDSP論壇的帖子，有很多人發(fā)出疑問：為什么我計(jì)算出的mag值全為零？這樣的帖子真不少見；TI的官方工程師論壇不少老外也發(fā)問：Whydi

5、dIgetallzeros?事實(shí)上我想主要原因是輸入數(shù)據(jù)格式不對(duì)。我認(rèn)為使用這個(gè)庫主要注意一下兩點(diǎn)：1. 數(shù)據(jù)輸入輸出的Q格式；2. 存儲(chǔ)空間分配。下面以32位實(shí)數(shù)FFT為例來說明。注意到文檔的40頁有如下說明：1. 在函數(shù)voidcalc（RFFT32_handle）有如下一句：NotethattheinputandoutputdataareinQ31format.；2. 在函數(shù)voidmag（FFT128R_handle）有如下一句：NotethatthemagnitudeoutputisstoredinQ30format.因?yàn)?8x系列DSP是定點(diǎn)處理器，而FFT計(jì)算涉及到不少浮點(diǎn)計(jì)算，

6、TI使用Q格式來解決這個(gè)問題（Q格式說明可參考sprcO87_IQmath）。事實(shí)上輸入數(shù)據(jù)采用Q31格式能在避免計(jì)算溢出前提下獲得最好的計(jì)算精度。對(duì)AD采樣的數(shù)據(jù)進(jìn)行FFT計(jì)算，定義計(jì)算緩沖區(qū)數(shù)組：longipcbN+2;因?yàn)锳D結(jié)果寄存器是12位的，在數(shù)據(jù)左對(duì)齊的情況下直接左移15位即可：ipcbConversionCount二（unsignedlong）AdcRegs.ADCRESULT0）15;另外我記得有份文檔提到在某些存儲(chǔ)器下是右對(duì)齊，此時(shí)則需左移19位，大家在使用時(shí)注意這個(gè)問題。當(dāng)采樣完成后，按照FFT庫的文檔上的說明或者仿照文檔附帶的例程，進(jìn)行相應(yīng)函數(shù)調(diào)用來實(shí)現(xiàn)自己的FFT計(jì)算

7、，比如按計(jì)算點(diǎn)數(shù)來定義各個(gè)變量數(shù)組，是否加窗，是否求解幅值平方值等。在進(jìn)行存儲(chǔ)器分配時(shí)，文檔上要求（128點(diǎn)實(shí)數(shù)FFT為例）：FFTipcbALIGN（256）:L0L1RAMPAGE1FFTmagL0L1RAMPAGE1FFTtfNVMEMPAGE0/*Nonvolatilememory*/.econstNVMEMPAGE0/*Nonvolatilememory*/注意兩點(diǎn)：FFT計(jì)算緩沖區(qū)FFTipcb需在page1上連續(xù)分配2N個(gè)位置（以ALIGN來指定），F(xiàn)FTtf（旋轉(zhuǎn)因子）需放在Nonvolatilememory的pageO內(nèi)（事實(shí)上如何才能為Nonvolatile我也不清楚）。F

8、FTtf位置這點(diǎn)我之前在調(diào)試時(shí)對(duì)計(jì)算FFT影響很大，因?yàn)閠widdlefactor若因存儲(chǔ)沖突肯定會(huì)造成查表值不準(zhǔn)確，那計(jì)算FFT時(shí)肯定就不對(duì)了。有個(gè)論壇帖子作者說一定要放在origin=0x008000開始位置，其實(shí)也不對(duì)，大家可自己去試驗(yàn)；后面我也會(huì)給出我的存儲(chǔ)配置文件（即.cmd文件內(nèi)容）。無圖無真相，下面給出計(jì)算實(shí)例。假設(shè)有一信號(hào)包含兩個(gè)諧波頻率值，分別為413.0Hz（幅值設(shè)為1.00V）和287.0Hz（幅值設(shè)為0.400V），利用函數(shù)發(fā)生器產(chǎn)生這兩路信號(hào)再混合，加上偏置后送入AD采樣。設(shè)采樣頻率1024Hz，共采樣2048點(diǎn)。圖1的采樣點(diǎn)均是右對(duì)齊的12位采樣結(jié)果值。AD采樣得到

9、的采樣點(diǎn)圖利用TI的FFT庫進(jìn)行計(jì)算，查看mag數(shù)組，得到圖2.以1024Hz采樣2048點(diǎn)，采樣時(shí)間2s,對(duì)兩個(gè)諧波頻率可采樣到整數(shù)倍周期；從另一個(gè)角度理解，此時(shí)最小頻率分辨率為0.5Hz,413.0Hz與287.0Hz均是其整數(shù)倍數(shù)，故不會(huì)發(fā)生頻譜展寬或混疊情況，計(jì)算得到的頻譜圖應(yīng)該為兩根尖峰線。從圖2結(jié)果也能看出這一點(diǎn)。查看mag數(shù)組，可知第一根尖峰線下標(biāo)574,第二根尖峰線下標(biāo)826,故真實(shí)頻率值分別為：574*0.5=287Hz,826*0.5=413Hz若要計(jì)算幅值，按照輸出的Q30格式除相應(yīng)系數(shù)即可。注意最好另外定義浮點(diǎn)數(shù)組來做除運(yùn)算，因?yàn)檎蛿?shù)據(jù)做除運(yùn)算（或者右移位操作）會(huì)丟失

10、小數(shù)位數(shù)據(jù)。若要驗(yàn)證DSP的FFT計(jì)算結(jié)果，可將AD數(shù)據(jù)從CCS導(dǎo)入到MATLAB中做對(duì)比計(jì)算。有些同學(xué)不清楚如何導(dǎo)入導(dǎo)出，下面說一下步驟。1. 點(diǎn)菜單欄file-data-save，選擇保存類型Integer（若是其他進(jìn)制還需在MATLAB中轉(zhuǎn)換），點(diǎn)“確定”后，在address”欄填入要保存數(shù)據(jù)的起始地址（填變量數(shù)組名或真實(shí)存儲(chǔ)器地址皆可），在“l(fā)ength”欄內(nèi)填入數(shù)據(jù)長度，“page”肯定選“data”頁了；全部設(shè)好后點(diǎn)“0K”2. 找到保存的數(shù)據(jù)文件，將后綴改成.txt,再用記事本打開（也可以不改后綴直接用記事本打開），刪去第一行數(shù)據(jù)；3. 打開MATLAB，點(diǎn)開“l(fā)oaddataf

11、ile”，選中剛才的數(shù)據(jù)文件，然后按照提示一步步往下導(dǎo)入即可。最后不妨用變量temp來保存這些數(shù)據(jù)。在MATLAB內(nèi)運(yùn)行如下代碼：%f_sample=1024;%采樣頻率N=2048;%采樣點(diǎn)Ny=temp;%temp即為導(dǎo)入數(shù)據(jù)的變量名n=0:N-1;t二n/f_sample;%做采樣點(diǎn)圖plot(t,y)figure;stem(t,y,.);figure;%fft變換并作圖fft_result二fft(y);mag=abs(fft_result)/(N/2);mag(1)=0;%為觀察諧波分量，此處特意將直流分量置為0%求解真實(shí)頻率值f=(0:length(fft_result)-l)*f

12、_sample/length(fft_result);%作頻譜圖stem(f(l:N/2),mag(l:N/2),b.);grid;%可將MATLAB計(jì)算結(jié)果與CCS內(nèi)的結(jié)果做些對(duì)比。附：存儲(chǔ)器分配(.cmd文件內(nèi)容)，我用的是2808的板子，2812或其他的稍作改動(dòng)即可。MEMORYPAGE0:BEGIN:origin=0x000000,length=0x000002RAMM0:origin=0x000002,length=0x0003FEPRAMH0:origin=0x3FA000,length=0x002000RESET:origin=0x3FFFC0,length=0x000002BO

13、OTROM:origin=0x3FF000,length=0x000FC0TESARAM:origin=0x008000,length=0x002000PAGE1:BOOT_RSVD:origin=0x000400,length=0x000080RAMM1:origin=0x00A000,length=0x001000L0L1RAM:origin=0x00B000,length=0x001000HL0SARAM:origin=0x3F8000,length=0x001010DRAMH0:origin=0x3F9010,length=0x000ff0SECTIONS/*Setupforboott

14、oSARAMmode:Thecodestartsection(foundinDSP28_CodeStartBranch.asm)re-directsexecutiontothestartofusercode.*/codestart:BEGIN,PAGE=0ramfuncs:RAMM0PAGE=0.text:PRAMH0,PAGE=0.cinit:RAMM0,PAGE=0.pinit:RAMM0,PAGE=0.switch:RAMM0,PAGE=0.reset:RESET,PAGE=0,TYPE=DSECT/*notused,*/FFTipcbALIGN(4096):HL0SARAM,PAGE=

15、1FFTtf:DRAMH0,PAGE=1FFTmag:L0L1RAM,PAGE=1.const:DRAMH0,PAGE=1.bss:.stack:DRAMH0,PAGE=1RAMM1,PAGE=1.sysmem:RAMM1,PAGE=1.ebss:DRAMH0,PAGE=1.econst:PRAMH0,PAGE=0.esysmem:L0L1RAM,PAGE=1IQmath:PRAMH0,PAGE=0IQmathTables:BOOTROM,type=NOLOAD,PAGE=0另外，有些同學(xué)一開始運(yùn)行FFT附帶的例程，會(huì)遇上兩個(gè)問題：一是缺少函數(shù)文件，這個(gè)到TI官網(wǎng)下那個(gè)sprc083_SGEN

16、包就好了，或找我我用郵箱發(fā)給大家也行；二是提示編譯不成功，找到上述配置文件的這一行.reset:RESET,PAGE=0,TYPE=DSECT/*notused,*/在后面加上TYPE=DSECT/*notused,*/。這些是我憑印象寫的，因好久沒做這個(gè)FFT，可能會(huì)寫出錯(cuò)誤，歡迎大家提出來：-)有問題就解決問題之TI庫實(shí)現(xiàn)2048點(diǎn)或更大點(diǎn)數(shù)FFT2010-05-1813:20星期二TI提供的FFT庫(sprc081)使用起來很方便，只要直接調(diào)用即可；更重要的是其中做了好些優(yōu)化，而且是用匯編寫成，運(yùn)算速度肯定要比自己寫的C代碼更快好多。這個(gè)庫最大能運(yùn)算1024點(diǎn)的實(shí)數(shù)FFT，但實(shí)際中我們可

17、能還需要更大點(diǎn)數(shù)的FFT，那如何去修改呢？其實(shí)修改的地方并不多，以前我曾花了一段時(shí)間瀏覽了一下那些頭疼的匯編代碼，發(fā)現(xiàn)模塊化寫得很好啊(佩服TI的工程師！)，這樣修改起來就快多了。1. 2048點(diǎn)FFT直接在工程中加入原來的庫文件fft.lib，并在代碼中作如下修改：第一步：在有關(guān)FFT代碼的頭文件fft.h中增添如下一段代碼：#defineRFFT32_2048P_DEFAULTS(long*)NULL,(long*)NULL,1024,10,(long*)NULL,(long*)NULL,0,0,1,(void(*)(void*)CFFT32_init,(void(*)(void*)CFF

18、T32_calc,(void(*)(void*)RFFT32_split,(void(*)(void*)RFFT32_mag,(void(*)(void*)RFFT32_win結(jié)構(gòu)體RFFT32關(guān)鍵說明：1. 第一個(gè)(long*)NULL,，是計(jì)算緩沖區(qū)(ipcb)的指針；2. 第二個(gè)(long*)NULL,，是查表旋轉(zhuǎn)因子(ft)的指針；3. 1024指的是FFT的大小，注意這里是2048點(diǎn)實(shí)數(shù)FFT被打包成1024點(diǎn)的復(fù)數(shù)，這是這個(gè)庫的實(shí)數(shù)FFT運(yùn)算的特點(diǎn)，這樣節(jié)省一半運(yùn)算時(shí)間和存儲(chǔ)空間；4. 10即1024是2的10次方，因?yàn)檫@里采用的是基2時(shí)間抽取算法；5. 下面的兩個(gè)(long*)N

19、ULL,分別是運(yùn)算結(jié)果中存放幅值(mag)以及加窗表格(win)的指針；6. 下面的兩個(gè)0分別是計(jì)算結(jié)果中的峰值幅值與峰值頻率，這里初始化0；7. 下面的1指的是查表的步進(jìn)值，因在這個(gè)sprc081文檔的庫已經(jīng)有Q31格式(以及Q30格式)的1024點(diǎn)的正弦因子表，因此對(duì)于1024點(diǎn)復(fù)數(shù)FFT查表步進(jìn)值設(shè)為1即可。8. 下面其他的就是要調(diào)用的函數(shù)指針了，不再贅述。明白上述結(jié)構(gòu)體的定義后，改寫其他點(diǎn)數(shù)就很容易了，比如這個(gè)庫最少能做到128點(diǎn)的實(shí)數(shù)FFT，我們還可以相應(yīng)改成64點(diǎn)或者更低點(diǎn)數(shù)的實(shí)數(shù)FFT，比如64點(diǎn)的結(jié)構(gòu)體就是如下的一段：#defineRFFT32_64P_DEFAULTS(lon

20、g*)NULL,(long*)NULL,32,5,(long*)NULL,(long*)NULL,0,0,32,(void(*)(void*)CFFT32_init,(void(*)(void*)CFFT32_calc,(void(*)(void*)RFFT32_split,(void(*)(void*)RFFT32_mag,(void(*)(void*)RFFT32_win第二步：修改源程序代碼里的變量定義：#includefft.h#defineN2048/定義N的點(diǎn)數(shù)，2048點(diǎn)/分別在數(shù)據(jù)存儲(chǔ)區(qū)為ipcb和mag分配足夠的地址，具體見下面的連接器文件#pragmaDATA_SECTIO

21、N(ipcb,FFTipcb);#pragmaDATA_SECTION(mag,FFTmag);/定義變量RFFT32fft=RFFT32_2048P_DEFAULTS;longipcbN+2;longmagN/2+l;/其他部分函數(shù)調(diào)用都是一樣的，這里不再列出第三步：在連接器(.cmd)文件里分配存儲(chǔ)地址：這里給出我在2808上的分配方法：MEMORYPAGE1:L0L1RAM:origin=OxOOBOOO,length=0x001000HL0SARAM:origin=0x3F8000,length=0x001010DRAMH0:origin=0x3F9010,length=0x000ff

22、0SECTIONSFFTipcbALIGN(4096):HL0SARAM,PAGE=1FFTtf:DRAMH0,PAGE=1FFTmag:L0L1RAM,PAGE=1好，到此為止就修改好了，可以使用這個(gè)FFT庫實(shí)現(xiàn)2048點(diǎn)實(shí)數(shù)FFT了。2.4096點(diǎn)或更大點(diǎn)數(shù)的FFT第一步：修改匯編代碼并制成庫文件其實(shí)要實(shí)現(xiàn)4096點(diǎn)的FFT需要改動(dòng)的地方也不多，主要因?yàn)樵瓉淼倪@個(gè)庫只提供了1024點(diǎn)的正弦值，我們需要自己算出2048點(diǎn)的Q31(或Q30)格式的一周期正弦值，并將其代替原來匯編代碼里的1024點(diǎn)的表格。匯編代碼的其他地方不需要改動(dòng)，不然可能產(chǎn)生不可預(yù)知的后果。在匯編代碼文件sel_q.asm

23、里查看設(shè)置：;SelectTwiddlefactorQformatTF_QFMAT.setQ30若是TF表設(shè)為Q30,則只需將計(jì)算出的正弦值移位成Q30格式即可。另外需注意的是,2048點(diǎn)的FFT需要2048點(diǎn)的TF表，但這個(gè)庫將其優(yōu)化成只需2048X3/4=1536點(diǎn)，即1.5K大小的表格(見sprc081文檔第10頁，注意到正弦的第二部分與余弦的第一部分是相重合的，這樣便節(jié)省了1/4周期的點(diǎn)數(shù))。計(jì)算好正弦值并加入到匯編代碼里，然后將所有的匯編代碼導(dǎo)進(jìn)CCS里做成一個(gè)新的庫，到時(shí)直接添加到工程里調(diào)用。為了方便這里我已將其做成一個(gè)新的庫命名為FFT_new.lib，大家直接用即可。第二步：修改頭文件fft.h結(jié)構(gòu)體定義增添如下一段代碼：#defineRFFT32_4096P_DEFAULTS(long*)NULL,(long*)NULL,2048,11,(long*)NULL,(long*)NULL,0,0,1,(void(*)(void*)CFFT32_init,(void(*)(void*)CFFT32_calc,(void(*)(void*)RFFT32_spl