專題八之三-通用數(shù)字信號處理方法的DSP實現(xiàn)v2013

專題八

通用數(shù)字信號處理方法的DSP實現(xiàn)DSP常見的幾種信號處理算法：DSP基本算術(shù)運算指令除法運算平方根運算級數(shù)展開產(chǎn)生正弦波FIR濾波器的實現(xiàn)FFT的實現(xiàn)自適應(yīng)濾波（LMS）的實現(xiàn)實現(xiàn)整數(shù)（定點）加法指令C5400，C5500中提供了多條用于加法的指令，如ADD。加法指令中不同的情況，如用于無符號數(shù)的加法運算，用于帶進(jìn)位的加法運算（如32位擴(kuò)展精度加法），以及用于立即數(shù)的加法。它們的具體加法指令有所不同。要實現(xiàn)飽和運算加法，需要設(shè)置標(biāo)志位SATD。C5500下還提供了同時完成加/減運算的指令A(yù)DDSUB等，以及條件執(zhí)行加減，輔助寄存器加法等增強(qiáng)加法指令。實現(xiàn)C6000加法指令C6000也提供了ADD加法指令。無符號加法（ADDU），飽和運算（SADD）加法，立即數(shù)加法（ADDK）等都有專用指令完成，等等也提供同時完成加/減運算的ADDSUB、ADDSUB2指令。也提供ADD2，ADD4的SIMD操作。C674x以及C6700浮點加法指令使用ADDSP，ADDDP完成循環(huán)尋址的加法C5500的循環(huán)尋址通過ARx輔助寄存器在指令中實現(xiàn)，可在任意間接尋址模式中使用。使用時可設(shè)置對應(yīng)的ARnLC比特位（ST2_55）。也可在指令后加.CRC6000的循環(huán)尋址可以通過指令A(yù)DDAB，ADDAD，ADDAH，ADDAW，并使用B4-B7/A4-A7來實現(xiàn)實現(xiàn)減法操作減法運算指令與加法類似，畢竟加一個負(fù)數(shù)就等效于做一個減法操作。減法指令中，SUBC為移位減，DSP中的除法就是用該指令來實現(xiàn)的。C5400，C5500，C6000都提供SUBC指令。參考教材的5.1.1小節(jié)（Page247），如何實現(xiàn)整數(shù)除法以及求模運算。教材Page252程序示例C6X整數(shù)除法C程序。實現(xiàn)16位定點整數(shù)除法在C54X中沒有提供專門的除法指令，一般有兩種方法來完成除法。一種是用乘法來代替，除以某個數(shù)相當(dāng)于乘以其倒數(shù)，所以先求出其倒數(shù)，然后相乘。這種方法對于除以常數(shù)特別適用。另一種方法是使用SUBC指令，重復(fù)16次減法完成除法運算。temp1/temp2為例，說明如何使用SUBC指令實現(xiàn)整數(shù)除法:

ldtemp1,T

；將被除數(shù)裝入T寄存器

mpytemp2,A

；除數(shù)與被除數(shù)相乘，結(jié)果放入A寄存器

ldtemp2,B

；將除數(shù)temp2裝入B寄存器的低16位

absB

；求絕對值

stlB,temp2

；將B寄存器的低16位存回temp2

ldtemp1,B

；將被除數(shù)temp1裝入B寄存器的低16位

absB

；求絕對值

rpt#15

；重復(fù)SUBC指令16次

subc temp2,b

；使用SUBC指令完成除法運算

bcd div_end,agt

；延時跳轉(zhuǎn)，先執(zhí)行下面兩條指令，

;然后判斷A，若A>0，則跳轉(zhuǎn)到標(biāo)號

;div_end，結(jié)束除法運算stlB,quot_i ;將商（B寄存器的低16位）存入變量quot_isthB,remain_i

；將余數(shù)（B寄存器的高16位）存入變量remain_ixorB

；若兩數(shù)相乘的結(jié)果為負(fù)，則商也應(yīng)為負(fù)。Subquot_i,B

；將商反號stlB,quot_i

；存回變量quot_i中div_end:

實現(xiàn)乘法運算乘法指令在C5500，C6000中都有大量的指令來實現(xiàn)。在C5000中提供了大量的乘法運算指令，16位相乘，其結(jié)果都是32位，放在累加器中。乘數(shù)在C5000的乘法指令很靈活，可以是T寄存器、立即數(shù)、存貯單元和累加器的高16位。有符號與無符號數(shù)乘法是不同的指令。實現(xiàn)乘法運算-C6000C6000也有MPY乘法指令，但其乘法的乘數(shù)有16位，32位，其結(jié)果也有32位，64位的不同指令。也提供有符號與無符號的區(qū)別指令。C674x以及C67+提供MPYSP/MPYDP浮點數(shù)乘法指令。C6000的乘法指令需要2個或更多的時鐘周期，如MPYSP（4），MPYDP（10），等。具體請參考指令手冊。小數(shù)乘法小數(shù)乘法與整數(shù)乘法主要差別在乘積多出一個符號比特，需要左移一位來消除C5000下使用狀態(tài)寄存器的FRCT位來控制，當(dāng)FRCT=1時進(jìn)行小數(shù)乘法。C6000使用獨立指令來實現(xiàn)乘積左移，實現(xiàn)小數(shù)乘法，如SMPY，SMPY32，SMPY2，SMPYHL，等待16位定點整數(shù)乘法例子在C54X中，小數(shù)的乘法與整數(shù)乘法基本一致，只是由于兩個有符號的小數(shù)相乘，其結(jié)果的小數(shù)點的位置在次高的后面，所以必須左移一位，才能得到正確的結(jié)果。C54X中提供了一個狀態(tài)位FRCT，將其設(shè)置為1時，系統(tǒng)自動將乘積結(jié)果左移移位。rsbx FRCT

；清FRCT標(biāo)志，準(zhǔn)備整數(shù)乘ld temp1,T

；將變量temp1裝入T寄存器mpy temp2,a

；完成temp2*temp1，結(jié)果放；入A寄存器（32位）16位定點小數(shù)乘法例子0ccdH（十進(jìn)制的0.1）x0599aH(十進(jìn)制的0.7），兩數(shù)相乘后B寄存器的內(nèi)容為08f5f0a4H（十進(jìn)制的0.07000549323857）。可以使用RND或使用MPYR指令對低16位做四舍五入的處理。SsbxFRCT

；FRCT=1，準(zhǔn)備小數(shù)乘法ld temp1,16,a

；將變量temp1裝入寄存器A的高16位mpya temp2

；完成temp2乘寄存器A的高16位，結(jié)；果在B中，同時將temp2裝入T寄存器sth b,mpy_f

；將乘積結(jié)果的高16位存入變量mpy_f實現(xiàn)16定點小數(shù)除法在C54X中實現(xiàn)16位的小數(shù)除法與前面的整數(shù)除法基本一致，也是使用循環(huán)的SUBC指令來完成。但有兩點需要注意：第一，小數(shù)除法的結(jié)果一定是小數(shù)（小于1），所以被除數(shù)一定小于除數(shù)。這與整數(shù)除法正好相反。所以在執(zhí)行SUBC指令前，應(yīng)將被除數(shù)裝入A或B寄存器的高16位，而不是低16位。其結(jié)果的格式與整數(shù)除法一樣，A或B寄存器的高16位為余數(shù)，低16位為商。第二，與小數(shù)乘法一樣，應(yīng)考慮符號位對結(jié)果小數(shù)點的影響。所以應(yīng)對商右移一位。乘累加運算-C5500MAC/MAS是DSP的核心運算，C5500提供直接的乘累加或乘累減指令。C5500在完成MAC/MAS同時，還可以直接完成數(shù)據(jù)的搬移，如MACMZC5500提供MAC/MPY/MOV的并行執(zhí)行能力，如MAC::MAC，MAC::MAS，MAC::MPY，MACMZ::MOV，等等乘累加運算-C6000C62x沒有直接的MAC指令，主要通過功能單元的并行能力來實現(xiàn)這個DSP核心操作！但從C64開始，增加了DOTP2，DOTPN2等乘累加操作指令。從C64+開始，增加了DDOTP2等多次乘累加操作指令（需要pipeline配合）。在C66x，有DOTP4H等4個乘累加操作。特殊乘法指令C5500的FIRSADD/FIRSSUB，LMS，SQA，SQDST，SQR等等C64+的復(fù)數(shù)乘法：CMPYC64+的Galoisfield乘法（伽羅瓦域）：GMPY，XORMPY,等?？捎糜诩用芴幚淼亩囗検竭\算C6700的平方根指令：RSQRSP，等C6700的倒數(shù)指令：RCPSP，RCPDP，等除法的實現(xiàn)5.1.4小節(jié)（Page257）除法的實現(xiàn)DSP中沒有直接除法指令5.1.1小節(jié)和5.1.2小節(jié)可以使用SUBC指令實現(xiàn)整數(shù)和小數(shù)除法更通常的做法是將除法改為乘以倒數(shù)。故除法問題轉(zhuǎn)為倒數(shù)計算。5.1.4小節(jié)給出了使用牛頓迭代法計算倒數(shù)。Newton-RhapsonC6700提供倒數(shù)計算指令，但精度有限。牛頓迭代法計算倒數(shù)先定義函數(shù)，a為輸入?yún)?shù)則f(x)的根就是a的倒數(shù)。利用牛頓迭代法求去f(x)的根。假設(shè)f’(x)為f(x)的一階導(dǎo)數(shù)，則如果x[n]為f(x)的根，則迭代方程為：Matlab代碼%計算輸入數(shù)據(jù)的倒數(shù)，其迭代的初始值應(yīng)該保證x2始終為正！a=input('Plsinputthenumber=');x2=0.3;error=1;number=0;whileerror>0.00001x1=x2;x2=x1*(2-a*x1);number=number+1;error=abs(x2-x1);disp([num2str(number),'',num2str(x2)]);enddisp(['Resultis',num2str(x2)]);C6700的除法實現(xiàn)教材的Page259提供了單精度浮點除法運算的匯編程序?qū)嵗a。同時給出了擴(kuò)展精度的倒數(shù)計算方法。教材Page260給出了C代碼函數(shù)實現(xiàn)除法例子。如果是其它定點DSP，可以參考前面的定點數(shù)運算，用牛頓迭代法計算倒數(shù)，從而實現(xiàn)定點除法計算！平方根的實現(xiàn)5.2小節(jié)（Page267）平方根計算求平方根也是數(shù)字信號處理中常見的算術(shù)運算之一。利用上一節(jié)介紹的Newton-Rhapson迭代算法，可以很方便地完成平方根或者平方根倒數(shù)的運計算。將平方根運算作如下變換：由上式可以看出，求平方根和求平方根倒數(shù)的運算實際上是等價的。平方根計算因此，定義函數(shù)只要求得方程g(y)=0的根，進(jìn)行x=ay的運算可得到a的平方根。為了使用Newton-Rhapson迭代算法，計算g(y)的一階導(dǎo)數(shù)并帶入迭代方程：平方根計算-C6700若使用C67xx，可將y[0]=RSQRSP(a)為初始值，即使用C67xx的求倒數(shù)的指令計算值作為初值。X0=_rsqrsp(arg1) ;getinitialTLUvalueto8bitsV=0.5*arg1 ;generatetermjustonceX1=X0*(1.5-V*X02) ;1stiterationto16-bitsX2=X1*(1.5–V*X12) ;2nditerationto32-bitsAns=arg1*X2 ;generatesquareroot*_rsqrsp為C6700的平方根計算指令，這里用作初值，提高精度。級數(shù)求和的實現(xiàn)5.3小節(jié)（PLOY指令的應(yīng)用）利用POLY指令計算泰勒級數(shù)三角函數(shù)、對數(shù)、指數(shù)等超越函數(shù)都可以用級數(shù)展開，如常用的泰勒級數(shù)。不同的級數(shù)展開有不同的收斂速度。以指數(shù)的展開為例：（實際一般取9項）利用POLY指令計算泰勒級數(shù)利用POLY指令計算泰勒級數(shù)POLY指令的含義:polySmemRound(A(32-16))xT+B->ASmem<16->B指數(shù)展開

x->TSmem->1/n!taylor:STMa9,AR3 ;AR3pointstocoefficientin ;Taylor’sequ.LDX,T ;setuprunningenvironmentusingLD*AR3+,16,A;powerfulpolyinstructionon54xDSPLD*AR3+,16,B

RPT#7 ;loop8timesenoughforaudioapp.POLY*AR3+ ;AH=fractionalpartinQ15 ;format利用POLY指令計算泰勒級數(shù)1,7,46,273,1365,5461,16384,32767,0,0C5500計算級數(shù)C5400提供了POLY指令實現(xiàn)級數(shù)迭代計算，但C5500，C6000都沒有類似指令。但C5500可以通過并行指令來實現(xiàn)：上面的并行指令完全可以替代POLY的功能C6000則通過功能單元并行來實現(xiàn)級數(shù)運算MACAC1,T0,AC2,AC1::MOV*AR1+<<#16,AC2正弦波的產(chǎn)生5.4小節(jié)（Page278）查表法產(chǎn)生正弦波可以構(gòu)造一個查找表，直接給出正弦函數(shù)的值。需要建立一個±/2范圍輸入所對應(yīng)的正弦(余弦)函數(shù)的值。不妨將表的大小初設(shè)為257項,表中的第一個值對應(yīng)于0°，最后一個值對應(yīng)于180°，或者說。這樣，表中相鄰兩點之間的間隔為:180/256=0.7031250°正弦函數(shù)和余弦函數(shù)之間只有一個90°的相移。查表法產(chǎn)生正弦波表格中的第1項是sin(0°)的值，第2項是sin(0.7031250°)的值，第3項是sin(1.406250°)的值，依此類推，最后一項是sin(180°)的值。利用該表所能構(gòu)造的波形的相位步進(jìn)，通常是0.703125°的整數(shù)倍。記錄波形的初始相位和當(dāng)前相位，就可以根據(jù)相位值去查找正弦表中相應(yīng)的位置，得到當(dāng)前輸出點的幅度值。查表法產(chǎn)生正弦波例如，在語音處理中，如果需要產(chǎn)生的正弦波的頻率為10kHz，采樣頻率為44.1kHz，對應(yīng)的相位步進(jìn)為：根據(jù)輸出信號的初始相位確定第一個點，則下一個輸出點的相位為81.6°，再下一個輸出點的相位為163.2°，依此類推。得到這些相位值后，就可以從查找表中的對應(yīng)位置，讀取當(dāng)前輸出點的幅度。需要產(chǎn)生余弦波時，將相位值加上/2（90°）數(shù)字振蕩器產(chǎn)生正弦數(shù)字振蕩器的本質(zhì)是，使用一個IIR（InfiniteImpulseResponse）濾波器，將一對共軛極點放在單位圓上，來產(chǎn)生正弦或余弦函數(shù)。利用正弦波sinx的指數(shù)形式（即系統(tǒng)的沖擊相應(yīng)）：可以得到正弦序列x[n]的z變換其中，A=2cos(T)，B=-1，C=sin(T)數(shù)字振蕩器產(chǎn)生正弦轉(zhuǎn)換為二階差分方程該系統(tǒng)的結(jié)構(gòu)框圖：通過遞推可以逐點計算系統(tǒng)輸出，即正弦信號數(shù)字振蕩器產(chǎn)生正弦例如，設(shè)該振蕩器的頻率為

，采

樣頻率為

，則

系數(shù)計算為：教材的5.4.2小節(jié)給出DSP實現(xiàn)的代碼例子

FIR濾波器

教材6章（6.3節(jié)）

FIR濾波器的結(jié)構(gòu)y(n)=h(0)x(n)+h(1)x(n-1）

+...+h(N-1)x[n-(N-1)]DSP實現(xiàn)FIR濾波器可以使用數(shù)字濾波器輔助設(shè)計軟件包或自行計算FIR濾波器的系數(shù)。本實驗例子中使用的是一個34階的對稱結(jié)構(gòu)的FIR低通濾波器，其采樣頻率Fs為25KHZ，通帶截止頻率

1.5KHZ，阻帶截止頻率為3KHZ，阻帶衰減為-40dB。(系數(shù)使用DFDP4自動生成）FIR濾波器的數(shù)據(jù)存儲方式FIR濾波器的DSP實現(xiàn)為了能正確使用循環(huán)尋址，數(shù)據(jù)緩沖區(qū)和系數(shù)的開始地址必須正確指定。濾波系數(shù)指針初始化時指向h(N-1)，經(jīng)過一次FIR濾波計算后，在循環(huán)尋址的作用下，仍然指向h(N-1)。而數(shù)據(jù)緩沖區(qū)指針指向的是需要更新的數(shù)據(jù)，如x(n)。在寫入新數(shù)據(jù)并完成FIR運算后，該指針指向x(n-(N-1))。使用循環(huán)尋址可以方便地完成濾波窗口數(shù)據(jù)的自動更新.FIR濾波器的DSP實現(xiàn)使用帶MAC指令的循環(huán)尋址模式實現(xiàn)FIR濾波器，程序片段如下：(輸入數(shù)據(jù)在AL中,濾波結(jié)果在AH中)STM #1,AR0 ;AR0=1STM #N,BK ;BK=N,循環(huán)尋址BUFFER大小為NSTL A,*FIR_DATA_P+%;更新濾波窗口中的采樣數(shù)據(jù)RPTZA,#(N-1)

；重復(fù)MAC指令N次,先將A清零MAC*FIR_DATA_P+0%,*FIR_COFF_P+0%,A

;完成濾波計算。注意FIR濾波系數(shù)存

；放在數(shù)據(jù)存儲區(qū)

.title"exampleforFIRfilter!".mmregs.globalmainstartOFF_INTR_3.set0Chk_win_size.set34;34tapsFIRLPFilterwin_data.set2000h;FIRdatawindow.datafilter_coff.word8ch;low_pass,1.5kHz(3kHz),34h(N-1).word59h;fs=25kHz,fc=1.5kHz.word0FF78h

.word0Fe56h.word0FF78h.word59hend_coff.word8ch;h(0)

.text_c_int00:

ssbxintm;disableallinterrupt!

stm#VECTOR,ar0;vector'startaddress->ar0

st#INSTR_B,*ar0(#OFF_INTR_3)

st#fir,*ar0(#OFF_INTR_3+1) ;initA/Dintvector!

ld#temp,dp;thefollowingcodesforMAC

st#win_data,t_ar3;ar3->2000hdatawindows

st#filter_coff,t_ar2;ar2->filtercoff

xora,a;cleara

xorb,b;clearb

rsbxintm;enableallint!g:idle1

bg;;Thesecodesmaybecalledbyserialrevint!allregistersdon't;change!;Whenenterthissubroutine,theADdatahasbeenputinA!;fir:

ld#temp,dp

stla,temp;a->ADdata!

calllow_pass_mac

lda,-16,b

ld#0,dp

stlb,2,TDXR;sentresulttoDA!

rete;********************************************************;LOWPASSFILTER(useMAC);Inputdata->A,outputdata->A;usedAR2->coff,AR3->data_buffer!;********************************************************low_pass_mac:

pshmst1

pshmst0

pshmar0

pshmbk

mvdm#t_ar2,ar2;restorear2

mvdm#t_ar3,ar3;restorear3

stm#1,ar0

stm#N,bk;setcircularaddressingsize

stla,*ar3+%

rptza,#(N-1);0->a,thenrepeat34times

mac*ar2+0%,*ar3+0%,a;doneFIRfilter,resultina

mvmdar3,#t_ar3;savear3

mvmdar2,#t_ar2;savear2

popmbk

popmar0

popmst0

popmst1

ret程序、數(shù)據(jù)的存儲器安排程序功能框圖相關(guān)外設(shè)的準(zhǔn)備：DSP，AD/DA，……相關(guān)外設(shè)的軟件設(shè)置：McBSP串口初始化、AC01的初始化、……硬件電路設(shè)計、調(diào)試軟件設(shè)計、調(diào)試FIR濾波器的工程實現(xiàn)初始化串口初始化AC01

等待新數(shù)據(jù)？調(diào)用濾波程序串口發(fā)送中斷服務(wù)程序串口接受中斷服務(wù)程序是llCCS圖形工具顯示FIR濾波效果圖形顯示FIR輸入/輸出頻譜FIRS指令來實現(xiàn)FIR濾波器一種有限單位沖激響應(yīng)呈現(xiàn)對中心點對稱的FIR濾波器。長度為N的線性相位FIR濾波器的輸出表達(dá)式為：FIRS指令來實現(xiàn)FIR濾波器FIRSXmem,Ymem,pmad含義：FIRS指令來實現(xiàn)FIR濾波器

B+（A(32-16))xPmad->B(Xmem+Ymem)<<16->APAR++Pmad尋址FIR濾波器系數(shù)，Xmem和Ymem分別指向窗口數(shù)據(jù)的上下兩部分。16點FIRS濾波數(shù)據(jù)存放*ar2*ar3new*ar2++16點FIRS濾波數(shù)據(jù)存放FIR*AR2+0%,*AR3+0%#FIR_COFAR2--,AR3-=2*ar2*ar3new*ar2++使用FIRS指令完成濾波利用FIRS指令,需要將輸入數(shù)據(jù)緩沖分成兩個,大小為N/2。初始狀態(tài)將AR2指到緩沖區(qū)1的頂部，AR3指到緩沖區(qū)2的底部。每次濾波之前,應(yīng)先將緩沖區(qū)1頂部的數(shù)據(jù)移到緩沖區(qū)2的底部,新來的一個樣本存儲到緩沖區(qū)1中時,并對緩沖區(qū)1指針AR2加1(使用循環(huán)尋址)。使用FIRS指令完成濾波處理器然后使用FIRS指令進(jìn)行乘加運算。當(dāng)然,在使用FIRS指令前,需要預(yù)先計算一次求和,以初始化A。在RPTZ重復(fù)指令和循環(huán)尋址的配合下,完成FIR濾波。濾波完成后，需要對兩個數(shù)據(jù)緩沖的指針進(jìn)行修正，以便對下一個點進(jìn)行處理。將Buffer1的指針減1和Buffer2的指針減2，使他們指向各自緩沖的數(shù)據(jù)隊列的最后。STM #1,AR0

；AR0=1STM #（N/2）,BK

；BK=N/2,循環(huán)尋址BUFFER大小為NMVDD*ar2,*ar3

；更新Buffer2STL A,*ar2+%

；更新濾波窗口中的采樣數(shù)據(jù)ADD*ar2+0%,*ar3+0%

；初始化ARPTZB,#(N/2-1)

；重復(fù)FIRS指令N/2次,先將B清零FIRS*ar2+0%,*ar3+0%,filter_coff+N/2

；完成濾波計算。注意FIR濾波系數(shù)存放；在程序存貯區(qū)，filter_coff為系數(shù)起始地址MAR*ar2-%

；修改Buffer1指針MAR*+ar3(-2)%

；修改Buffer2指針使用帶FIRS指令的循環(huán)尋址模式實現(xiàn)FIR濾波器，程序片段如下：(輸入數(shù)據(jù)在AL中,濾波結(jié)果在B中)FIRSADDAcx,ACy,Cmem,Xmem,YmemFIRSADD指令-C5500

ACy=ACy+(ACx*Cmem)::ACx=(Xmem<<16)+(Ymem<<16)FIRADD.CR可使得后面的Cmem，Xmem，Ymem間接尋址都是用循環(huán)尋址。在滿足C5500的并行條件時，可以與其它指令并行執(zhí)行。C6000實現(xiàn)濾波器-乘累加loop:

ldh.d1 *A8++,A2

||ldh.d2 *B9++,B3

nop 4 mpy.m1x A2,B3,A4

nop

add.l1 A4,A6,A6 sub.l2 B0,1,B0

[b0]b.s1 loop

nop 5Whatcanyou

putinparallel?LoadInstructionsMPY2x=+a1a0A0LDW.D1*A4++,A0x1x0B0||LDW.D2*B4++,B0A2A3a1*x1a0*x0MPY2A0,B0,A3:A2||MPYH.M2A0,B0,B5+a1x1+a3x3...a0x0+a2x2...ADD.L1A2,A6,A6||ADD.L2A3,B6,B6A6B6finalsumADD.L1A6,B6,A4A4+DOTP2withLDDW=+a2a0A1:A0LDDW.D1*A4++,A1:A0||LDDW.D2*B4++,B1:B0A2B2a3*x3+a2*x2a1*x1+a0*x0DOTP2A0,B0,A2||DOTP2A1,B1,B2+intermediatesumADDA2,A3,A3a1a3:A5x2x0B1:B0x1x3:finalsumADDA3,B3,A4A4+||ADDB2,B3,B3intermediatesumA3B3InCh8,we'llgetalltheseinstructionsworkinginparalleld0*c0d1*c1d2*c2d3*c3+yd4*c4d5*c5d6*c6d7*c7coefdataBlockRealFIRfor(i=0;I<ndata;i++){ sum=0; for(j=0;j<ncoef;j++){ sum=sum+(d[i+j]*c[j]); }

y[i]=sum;}loopIteration[i,j][0,0][0,1]d0c0d1c1d1c0d2c2d2c1d3c3d3c2...BlockRealFIRExample(DDTOPL2)for(i=0;I<ndata;i++){ sum=0; for(j=0;j<ncoef;j++){ sum=sum+(d[i+j]*c[j]); }

y[i]=sum;}loopIteration[i,j][0,0][0,1]d0c0

+d1c1d1c0+d2c1d2c2d3c3d3c2...Four16x16multipliesIneach.Muniteverycycle

addsupto8MACs/cycle,or 8000MMACSBottomLine:TwoloopiterationsforthepriceofoneDDOTPL2d3d2:d1d0,

c1c0,

sum1:sum0loopIteration[i,j][0,0][0,1][0,2][0,3][0,4][0,5]d0c0

+d1c1d1c0+

d2c1d2c2d2c0d3c3d3c2d3c1d3c0d4c4d4c3d4c2+d5c3d4c1d4c0d5c5d5c4d5c2+d6c3d5c1d6c0d6c6d6c5d6c4d6c2d6c1d7c7d7c6d7c5d7c4d7c3d6c2d8c7d8c6d8c5d8c4d6c3

DDOTPL2.M1d3d2:d1d0,c1c0,sum1:sum0||DDOTPL2.M2d7d6:d5d4,c3c2,sum3:sum2ParallelDDOTPL2’sRealBlockFIR[A_j]SPLOOPD4||MVC.S2B_i0,ILC;setILC||ADD.L2B_i0,1,B_i0;T/4||ADDAH.D2B_DLYaddr,nCoefs-1+4,B_DLYOUTaddr||MVK.S1nCoefs/4-3,A_T||ADDAB.D1A_DLYaddr,8,A_DLYaddr*-stageA*LDDW.D1T2*++A_DLYaddr,B_d7d6:B_d5d4;[1,1]||LDDW.D2T1*++B_DLYaddr,A_d3d2:A_d1d0;[1,1]LDDW.D2T2*++B_COEFaddr,B_c3c2:B_c1c0;[2,1]||LDDW.D1T1*++A_COEFaddr,A_c3c2:A_c1c0;[2,1]SPMASK||LDDW.D1T1*A_DLYaddr[1],A_dbda:A_d9d8;[3,1]||^LDDW.D2T2*B_INaddr++,B_TEMP1:B_TEMP0;ld1stinput||^MVC.S2B_i0,RILC;setRILCSPMASK||^LDDW.D2T1*B_INaddr++,A_TEMP1:A_TEMP0;ld2ndinput||^ZERO.L1A_st;clearstflag||^ADDAB.D1DP,outputs+8,A_OUTaddr||^MVK.S2nCoefs/4-1,B_TC||^MVK.S1nCoefs/4-1,A_TC*-stageB*SPMASK||^SUB.L2XA_OUTaddr,8,B_OUTaddrNOP1

DMV.S2XB_d5d4,A_d3d2,B_d5d4_:B_d3d2_;[7,1]||DDOTPL2.M1A_d3d2:A_d1d0,A_c1c0,A_pb:A_pf;[7,1]||DDOTPL2.M2B_d7d6:B_d5d4,B_c3c2,B_p2:B_p6;[7,1]SPMASK||DMV.S1XA_d9d8,B_d7d6,A_d9d8_:A_d7d6_;[8,1]||DDOTPL2.M2B_d5d4_:B_d3d2_,B_c3c2,B_pa:B_pe;[8,1]||DDOTPL2.M1A_dbda:A_d9d8,A_c3c2,A_p0:A_p4;[8,1]||^STNDW.D2T2B_TEMP1:B_TEMP0,*B_DLYOUTaddr++;st1stinputNewdouble-throughput16-bitmpyinstructionsABCABCABCDDDABCD4ABCABCABCDDDABCDO1O2O3O4O5e1e2,p1e3,p2e4,p3p4DDOTP4DDOTP4(.unit)src1,src2,dst_o:dst_eBlockRealFIRC64xImplementationDOTP2

d1d0,c1c0,s0;d[1]*c[1]+d[0]*c[0]

DOTP2

d2d1,c1c0,s1;d[2]*c[1]+d[1]*c[0]C64xPlusImplementationDDOTPL2d3d2:d1d0,c1c0,s1:s0Reduces.MrequirementbyhalfC64x: 194cycles,624bytes(N=40,T=16)

C64x+: 126cycles,496bytes(N=40,T=16)N=Lengthofexampleblock,T=datasizeCMPY...ComplexMultiply(CMPY)A0r1i1xxA1r2i2==CMPYA0,A1,A3:A2r1*r2-i1*i2:i1*r2+r1*i2

32-bits

32-bits

Four16x16multipliesper.MunitUsingtwoCMPYs,atotalofeight16x16multipliespercycleFloating-pointversion(CMPYSP)uses:64-bitinputs(registerpair)128-bitpackedproducts(registerquad)Youthenneedtoadd/subtracttheproductstogetthefinalresultsingle.MunitComplexMPYCMPY(.M)src1,src2,dst_o:dst_eFour16-bitinputs(real0,imag0,real1,imag1)Two32-bitoutputs(real,imag)CMPYR(.M)src1,src2,dstFour16-bitinputs(real0,imag0,real1,imag1)Onepacked32-bitoutput(16bitreal,16bitimag)Roundedbyadding2**15CMPYR1(.M)src1,src2,dstFour16-bitinputs(real0,imag0,real1,imag1)Onepacked32-bitoutput(real,imag)Roundedbyadding2**14plusleftshiftBlockComplexFIRC64xImplementationDOTP2 dre_dim,cim_cre,pim

DOTPN2 dre_dim,cre_cim,preC64xPlusImplementationCMPY dre_dim,cre_cim,pre:pimReduces.MrequirementbyhalfC64x: 674cycles,572bytes(N=40,T=16)

C64x+: 344cycles,460bytes(N=40,T=16)N=Lengthofexampleblock,T=datasizeDSPLIBOptimizedDSPFunctionLibraryforCprogrammersusingC62x/C67xandC64xdevicesTheseroutinesaretypicallyusedincomputationallyintensivereal-timeapplicationswhereoptimalexecutionspeediscritical.Byusingtheseroutines,youcanachieveexecutionspeedsconsiderablyfasterthanequivalentcodewritteninstandardANSIClanguage.Andtheseready-to-usefunctionscansignificantlyshortenyourdevelopmenttime.TheDSPlibraryfeatures:C-callableHand-codedassembly-optimizedTestedagainstCmodelandexistingrun-time-supportfunctionsAdaptive

filteringMathDSP_firlms2DSP_dotp_sqrCorrelationDSP_dotprodDSP_autocorDSP_maxvalFFTDSP_maxidxDSP_bitrev_cplxDSP_minvalDSP_radix2DSP_mul32DSP_r4fftDSP_neg32DSP_fftDSP_recip16DSP_fft16x16rDSP_vecsumsqDSP_fft16x16tDSP_w_vecDSP_fft16x32MatrixDSP_fft32x32DSP_mat_mulDSP_fft32x32sDSP_mat_transDSP_ifft16x32MiscellaneousDSP_ifft32x32DSP_bexpFilters&convolutionDSP_blk_eswap16DSP_fir_cplxDSP_blk_eswap32DSP_fir_genDSP_blk_eswap64DSP_fir_r4DSP_blk_moveDSP_fir_r8DSP_fltoq15DSP_fir_symDSP_minerrorDSP_iirDSP_q15toflTechnicalTrainingOrganizationT

TOFFT實現(xiàn)教材7章(7.3小節(jié)及7.4小節(jié))FFT是數(shù)字信號處理中重要的工具FFT是一種高效實現(xiàn)離散付氏變換的算法。(8.2.2小節(jié)介紹了Goertzel算法）

離散付氏變換的目的是把信號由時域變換到頻域，從而可以在頻域分析處理信息，得到的結(jié)果再由付氏逆變換到時域。DFT的定義為：DFT的定義可以方便的把它改寫為如下形式：WN（旋轉(zhuǎn)因子）的周期性是DFT的關(guān)鍵性質(zhì)之一。為了強(qiáng)調(diào)起見，常用表達(dá)式WN取代W以便明確其周期是N。FFT是DFT的快速算法由DFT的定義可以看出，在x[n]為復(fù)數(shù)序列的情況下，完全直接運算N點DFT需要（N-1）2次復(fù)數(shù)乘法和N（N-1）次加法。FFT的基本思想在于，將原有的N點序列序列分成兩個較短的序列，這些序列的DFT可以很簡單的組合起來得到原序列的DFT。FFT是DFT的快速算法例如，若N為偶數(shù)，將原有的N點序列分成兩個（N/2）點序列，那么計算N點DFT將只需要約[(N/2)2·2]=N2/2次復(fù)數(shù)乘法。即比直接計算少作一半乘法。該處理方法可以反復(fù)使用，即（N/2）點的DFT計算也可以化成兩個（N/4）點的DFT（假定N/2為偶數(shù)），從而又少作一半的乘法。這樣一級一級的劃分下去一直到最后就劃分成兩點的FFT運算的情況。FFT是DFT的快速算法比如，一個N=8點的FFT運算按照這種方法來計算FFT可以用下面的流程圖來表示：實數(shù)FFT運算對于離散傅立葉變換（DFT）的數(shù)字計算，F(xiàn)FT是一種有效的方法。一般假定輸入序列是復(fù)數(shù)。當(dāng)實際輸入是實數(shù)時，利用對稱性質(zhì)可以使計算DFT非常有效。一個優(yōu)化的實數(shù)FFT算法是一個組合以后的算法。原始的2N個點的實輸入序列組合成一個N點的復(fù)序列，之后對復(fù)序列進(jìn)行N點的FFT運算，最后再由N點的復(fù)數(shù)輸出拆散成2N點的復(fù)數(shù)序列。實數(shù)FFT運算使用這種方法，在組合輸入和拆散輸出的操作中，F(xiàn)FT運算量減半。這樣利用實數(shù)FFT算法來計算實輸入序列的DFT的速度幾乎是一般復(fù)FFT算法的兩倍。本實驗就用這種方法實現(xiàn)了一個256點實數(shù)FFT（2N=256）運算。實數(shù)FFT運算序列的存儲分配基二實數(shù)FFT運算的算法第一步，輸入數(shù)據(jù)的組合和位倒序：（1）把輸入序列作位倒序，是為了在整個運算最后的輸出中得到的序列是自然順序。（2）原始的輸入的2N=256個點的實數(shù)序列復(fù)制放到標(biāo)記有“d_input_addr”的相鄰單元，當(dāng)成N=128點的復(fù)數(shù)序列d[n]。其中，奇數(shù)地址是d[n]的實部，偶數(shù)地址是d[n]的虛部。復(fù)數(shù)序列經(jīng)過位倒序，存儲在數(shù)據(jù)處理緩沖器中“fft-data”。bit_rev:

STM#d_input_addr,ar3

；在AR3中放入輸入地址

STM#fft_data,ar7

；在AR7中放入處理后輸出的地址

MVMMDATA_PROC_BUF,ar2;AR2中裝入第一個位倒序

;數(shù)據(jù)指針

STM#K_FFT_SIZE-1,BRC

STM#K_FFT_SIZE,ar0;AR0=輸入數(shù)據(jù)數(shù)目的一半(128)

RPTBbit_rev_end

MVDD*ar3+,*ar7+ ;將原始輸入緩沖中的數(shù)據(jù)放入到

;位倒序緩沖中去之后輸入緩沖

;(AR3)指針加1,位倒序緩沖(AR2)指

;針也加1

MVDD*ar3-,*ar7+ ;將原始輸入緩沖中的數(shù)據(jù)放入到

;位倒序緩沖中去之后輸入緩沖(AR3);指針減1以保證位倒序?qū)ぶ氛_

MAR*ar3+0B

；按位倒序?qū)ぶ贩绞叫薷腁R3bit_rev_end:基二實數(shù)FFT運算的算法第二步，N點復(fù)數(shù)FFT：（1）在數(shù)據(jù)處理緩沖器里進(jìn)行N點復(fù)數(shù)FFT運算。由于在FFT運算中要用到旋轉(zhuǎn)因子WN，它是一個復(fù)數(shù)。我們把它分為正弦和余弦部分，用Q15格式將它們存儲在兩個分離的表中。每個表中有128項，對應(yīng)從0度到180度。因為采用循環(huán)尋址來對表尋址，128=27<28，因此每張表排隊的開始地址就必須是8個LSB位為0的地址。基二實數(shù)FFT運算的算法我們把128點的復(fù)數(shù)FFT分為七級來算，第一級是計算兩點的FFT蝶形結(jié)，第二級是計算四點的FFT蝶形結(jié)，然后是八點、十六點、三十二點六十四點、一百二十八點的蝶形結(jié)計算。最后所得的結(jié)果表示為：其中，R[k]、I[k]分別是D[k]的實部和虛部。D[k]=F{d[n]}=R[k]+jI[k]這一步中，實現(xiàn)FFT計算的具體程序如下：;stage1：計算FFT的第一步，兩點的FFT.asgAR1,GROUP_COUNTER ；定義FFT計算的組指針

.asgAR2,PX ；定義AR2為指向參加蝶形運算第一個數(shù)據(jù)的指針

.asgAR3,QX ；定義AR3為指向參加蝶形運算第二個數(shù)據(jù)的指針

.asgAR4,WR ；定義AR4為指向余弦表的指針

.asgAR5,WI ；定義AR5為指向正弦表的指針

.asgAR6,BUTTERFLY_COUNTER

；定義AR6為指向蝶形結(jié)的指針

.asgAR7,DATA_PROC_BUF

；定義在第一步中的數(shù)據(jù)處理緩沖指針

.asgAR7,STAGE_COUNTER

；定義剩下幾步中的數(shù)據(jù)處理緩沖指針

pshmst0

pshmar0

pshmbk

；保存環(huán)境變量

SSBXSXM

；開啟符號擴(kuò)展模式

STM#K_ZERO_BK,BK

；讓BK=0使*ARn+0%=*ARn+0

LD#-1,ASM

；為避免溢出在每一步輸出時右移一位

MVMMDATA_PROC_BUF,PX

；PX指向參加蝶形結(jié)運算的

;第一個數(shù)的實部（PR）

LD*PX,16,A

；AH:=PR

STM#fft_data+K_DATA_IDX_1,QX

；QX指向參加蝶形

;結(jié)運算的第二個數(shù)的實部（QR）

STM#K_FFT_SIZE/2-1,BRC

；設(shè)置塊循環(huán)計數(shù)器

RPTBDstage1end-1

；語句重復(fù)執(zhí)行的范圍到地址

;stage1end-1處

STM#K_DATA_IDX_1+1,AR0

；延遲執(zhí)行的兩字節(jié)的指令

;（該指令不重復(fù)執(zhí)行）

SUB*QX,16,A,B

；BH:=PR-QR

ADD*QX,16,A

；AH:=PR+QR

STHA,ASM,*PX+

；PR':=(PR+QR)/2

STB,*QX+

；QR':=(PR-QR)/2

||LD*PX,A

；AH:=PI

SUB*QX,16,A,B

；BH:=PI-QI

ADD*QX,16,A

；AH:=PI+QI

STHA,ASM,*PX+0%

；PI':=(PI+QI)/2

STB,*QX+0%

；QI':=(PI-QI)/2

||LD*PX,A

；AH:=nextPRstage1end:(R[0]+jI[0])+(R[1]+jI[1])*W0

;Stage2：計算FFT的第二步，四點的FFT

MVMMDATA_PROC_BUF,PX

；PX指向參加蝶形結(jié)

;運算第一個數(shù)據(jù)的實部（PR）

STM#fft_data+K_DATA_IDX_2,QX

；QX指向參加蝶形結(jié)

;運算第二個數(shù)據(jù)的實部（QR）

STM#K_FFT_SIZE/4-1,BRC

；設(shè)置塊循環(huán)計數(shù)器

LD*PX,16,A

；AH:=PR

RPTBDstage2end-1

；語句重復(fù)執(zhí)行的范圍到地址

;stage1end-1處

STM#K_DATA_IDX_2+1,AR0；初始化AR0以被循環(huán)尋址;以下是第二步運算的第一個蝶形結(jié)運算過程

SUB*QX,16,A,B

；BH:=PR-QR

ADD*QX,16,A

；AH:=PR+QR

STHA,ASM,*PX+

；PR':=(PR+QR)/2

STB,*QX+

；QR':=(PR-QR)/2

||LD*PX,A

；AH:=PI

SUB*QX,16,A,B

；BH:=PI-QI

ADD*QX,16,A

；AH:=PI+QI

STHA,ASM,*PX+

；PI':=(PI+QI)/2

STHB,ASM,*QX+

；QI':=(PI-QI)/2;以下是第二步運算的第二個蝶形結(jié)運算過程

MAR*QX+

；QX中的地址加一

ADD*PX,*QX,A

；AH:=PR+QI

SUB*PX,*QX-,B

；BH:=PR-QI

STHA,ASM,*PX+

；PR':=(PR+QI)/2

SUB*PX,*QX,A

；AH:=PI-QR

STB,*QX

；QR':=(PR-QI)/2

||LD*QX+,B

；BH:=QR

STA,*PX ;PI':=(PI-QR)/2

||ADD*PX+0%,A

；AH:=PI+QR

STA,*QX+0%

；QI':=(PI+QR)/2

||LD*PX,A

；AH:=PRstage2end:;Stage3thruStagelogN-1：從第三步到第六步的過程如下(略）基二實數(shù)FFT運算的算法第三步，分離復(fù)數(shù)FFT的輸出為奇部分和偶部分：分離FFT輸出為相關(guān)的四個序列：RP、RM、IP和IM，即偶實數(shù)，奇實數(shù)、偶虛數(shù)和奇虛數(shù)四部分，以便第四步形成最終結(jié)果。利用信號分析的理論我們把D[k]通過下面的公式分為偶實數(shù)RP[k]、奇實數(shù)RM[k]、偶虛數(shù)IP[k]和奇虛數(shù)IM[k]：RP[k]=RP[N-k]=0.5*(R[k]+R[N-k])RM[k]=-RM[N-k]=0.5*(R[k]-R[N-k])IP[k]=IP[N-k]=0.5*(I[k]+I[N-k])IM[k]=-IM[N-k]=0.5*(I[k]-I[N-k])RP[0]=R[0]IP[0]=I[0]RM[0]=IM[0]=RM[N/2]=IM[N/2]=0RP[N/2]=R[N/2]IP[N/2]=I[N/2]

這一過程的程序代碼(略）基二實數(shù)FFT運算的算法第四步，產(chǎn)生最后的N=256點的復(fù)數(shù)FFT結(jié)果:產(chǎn)生2N=256個點的復(fù)數(shù)輸出，它與原始的256個點的實輸入序列的DFT一致。輸出駐留在數(shù)據(jù)緩沖器中。通過下面的公式由RP[k]、RM[n]、IP[n]和IM[n]四個序列可以計算出a[n]的DFT：AR[k]=AR[2N-k]=RP[k]+cos(kπ/N)*IP[k]-sin(kπ/N)*RM[k]AI[k]=-AI[2N-k]=IM[k]-cos(kπ/N)*RM[k]-sin(kπ/N)*IP[k]AR[0]=RP[0]+IP[0]AI[0]=IM[0]-RM[0]AR[N]=R[0]-I[0]AI[N]=0其中：

A[k]=A[2N-k]=AR[k]+jAI[k]=F{a(n)}

這一過程的程序代碼(略）基二實數(shù)FFT運算的算法計算所求信號的功率：由于最后所得的FFT數(shù)據(jù)是一個復(fù)數(shù)，為了能夠方便的在虛擬頻譜儀上觀察該信號的特征，我們通常對所得的FFT數(shù)據(jù)進(jìn)行處理——取其實部和虛部的平方和，即求得該信號的功率。CCS中的圖形工具顯示FFT結(jié)果DSPLIBOptimizedDSPFunctionLibraryforCprogrammersusingC62x/C67xandC64xdevicesTheseroutinesaretypicallyusedincomputationallyintensivereal-timeapplicationswhereoptimalexecutionspeediscritical.Byusingtheseroutines,youcanachieveexecutionspeedsconsiderablyfasterthanequivalentcodewritteninstandardANSIClanguage.Andtheseready-to-usefunctionscansignificantlyshortenyourdevelopmenttime.TheDSPlibraryfeatures:C-callableHand-codedassembly-optimizedTestedagainstCmodelandexistingrun-time-supportfunctionsAdaptive

filteringMathDSP_firlms2DSP_dotp_sqrCorrelationDSP_dotprodDSP_autocorDSP_maxvalFFTDSP_maxidxDSP_bitrev_cplxDSP_minvalDSP_radix2DSP_mul32DSP_r4fftDSP_neg32DSP_fftDSP_recip16DSP_fft16x16rDSP_