DSP常見算法的實現(xiàn)

3.6常見的算法實現(xiàn)在實際應(yīng)用中雖然信號處理的方式多種多樣，但其算法的基本要素卻大多相同，在本節(jié)中介紹幾種較為典型的算法實現(xiàn)，希望通過對這些例子（單精度，16bit）的分析，能夠讓大家熟悉DSP編程中的一些技巧，在以后的工作中可以借鑒，達(dá)到舉一反三的效果。1.函數(shù)的產(chǎn)生在高級語言的編程中，如果要使用諸如正弦、余弦、對數(shù)等數(shù)學(xué)函數(shù)，都可以直接調(diào)用運(yùn)行庫中的函數(shù)來實現(xiàn)，而在DSP編程中操作就不會這樣簡單了。雖然 TI公司提供的實時運(yùn)行庫中有一些數(shù)學(xué)函數(shù)，但它們所耗費(fèi)的時間大多太長，而且對于大多數(shù)定點程序使用雙精度浮點數(shù)的返回結(jié)果有點“大材小用”的感覺，因此需要編程人員根據(jù)自身的要求“定制”數(shù)學(xué)函數(shù)。實現(xiàn)數(shù)學(xué)函數(shù)的方法主要有查表法、迭代法和級數(shù)逼近法等，它們各有特點，適合于不同的應(yīng)用。查表法是最直接的一種方法，程序員可以根據(jù)運(yùn)算的需要預(yù)先計算好所有可能出現(xiàn)的函數(shù)值，將這些結(jié)果編排成數(shù)據(jù)表，在使用時只需要根據(jù)輸入查出表中對應(yīng)的函數(shù)值即可。 它的特點是速度快，但需要占用大量的存儲空間，且靈活度低。 當(dāng)然，可以對上述查表法作些變通，僅僅將一些關(guān)鍵的函數(shù)值放置在表中，對任意一個輸入，可根據(jù)和它最接近的數(shù)據(jù)采用插值方法來求得。這樣占用的存儲空間有所節(jié)約，但數(shù)值的準(zhǔn)確度有所下降。迭代法是一種非常有用的方法，在自適應(yīng)信號處理中發(fā)揮著重要的作用。作為函數(shù)產(chǎn)生的一種方法，它利用了自變量取值臨近的函數(shù)值之間存在的關(guān)系，如時間序列分析中的 AR、MAARMA等模型，刻畫出了信號內(nèi)部的特征。因為它只需要存儲信號模型的參量和相關(guān)的狀態(tài)變量，所以所占用的存儲空間相對較少，運(yùn)算時間也較短。但它存在一個致命的弱點，由于新的數(shù)值的產(chǎn)生利用了之前的函數(shù)值，所以它容易產(chǎn)生誤差累積，適合精度要求不高的場合。級數(shù)逼近法是用級數(shù)的方法在某一自變量取值范圍內(nèi)去逼近數(shù)學(xué)函數(shù)， 而將自變量取值在此范圍外的函數(shù)值利用一些數(shù)學(xué)關(guān)系，用該范圍內(nèi)的數(shù)值來表示。這種方法最大的優(yōu)點是靈活度高，且不存在誤差累積，數(shù)值精度由程序員完全控制。該方法的關(guān)鍵在于選擇一個合適的自變量取值區(qū)間和尋找相應(yīng)的系數(shù)。下而通過正弦函數(shù)的實現(xiàn)，具體對上述三種方法作比較。查表法較簡單，只需要自制一張數(shù)據(jù)表，也可以利用 C5400DSPROM內(nèi)的正弦函數(shù)表。迭代法的關(guān)鍵是尋找函數(shù)值間的遞推關(guān)系。假設(shè)函數(shù)采樣時間間隔為 T,正弦函數(shù)的角頻率為，那么可以如下推導(dǎo)：令sinTsin sinT等式的左邊展開為left_sidesincosTcossinT等式的右邊展開為right_sidesincosTcossinT對比系數(shù)，可以得到2cosT, 1。令 nT,便可以得到如下的遞推式:sn2cosTsn1sn2

令s1 0,s2AsinT,逐一迭代就能夠獲得采樣間隔為 t的正弦序列了。從迭代公式可以更清楚地看岀，這種方法存在誤差累積。再來看級數(shù)逼近法，首先需要尋找一個合適的自變量取值區(qū)間和尋找相應(yīng)的系數(shù)。 從正弦函數(shù)的對稱性可知，只需要計算取值在 ［0,］內(nèi)的函數(shù)就可以推斷出所有取值范圍內(nèi)的函數(shù)。接下來尋找系數(shù)，對于sinx可以作如下變換sinxsiny值范圍在［0,1］內(nèi)，而sin.new()與sin()同構(gòu)，所以在下面的分析中將sin(),提到的正弦函數(shù)即指 sin_new()o若匯編編程時采用的數(shù)據(jù)為sin_new(y),那么y的取sin_new()sin_new(y),那么y的取sin_new()替代Q15格式,那么取2圖3-算法取值與弧度的對應(yīng)關(guān)系在［0,1］內(nèi)，正弦函數(shù)的修正級數(shù)(五次)展開如下式：sin_new(x)3.140625x0.02026367k2-5.325196x3 0.5446778xl 1.800293x5根據(jù)上式，可以寫岀正弦函數(shù)的生成程序。SinePolyCoeff:?wordOxlcceSinePolyCoeff:?wordOxlcce;inQ12formatA;1.800293(coefforx5=c5)stlmA,T；T=xstm#SinePolyCoeff,AR2Id*AR2+,16,A;AH=c5Id*AR2+,16,B;BH二c4poly*AR2+;A=c5*x+c4poly*AR2+;A=c5*xA2+c4*x+c3poly*AR2+;A=c5*xA3+c4*xA2+c3*x+c2poly*AR2+;A=c5*xA4+c4*xA3+c3*xA2+c2*x+clmpyaAsftaA,3;adjustAHtoQ15

?wordAOxaacc ;-5?325196(coefforx3=c3)?word0x0053 ;0.02026367(coefforxA2=c2)?word0x3240 ；3.140625(coefforxAl=cl)在編稈過稈中，使用到了 FOLY語句，它能夠使多項式的計算簡潔快速地完成。該函數(shù)的結(jié)果可以通過實驗X來驗證。2.FIR濾波器的實現(xiàn)FIR的FIR濾波器由于具有線性相位而且延遲能夠確定，因而在信號處理中廣泛應(yīng)用。FIR的x(n)h(0)z1z1h(1)?dy(n)x(n)h(0)z1z1h(1)?dy(n)基木模型如下圖所示。圖3-FIR基木模型如下圖所示。圖3-FIR模型N1其數(shù)學(xué)表達(dá)式為yn hii0ni,根據(jù)該表達(dá)式可以給出一種最為直接的實現(xiàn)形式。直接形式中采用線性地址來存放數(shù)據(jù)，如圖3-所示。式。直接形式中采用線性地址來存放數(shù)據(jù)，如圖3-所示。程序中可以采用MACD來實現(xiàn)程序如下:Id*(Input),AstlA,*(x_n)stm#x_n_Nml,AR2mpy*AR2~,#h_Nml,Brpt#N—2

macc*AR2~,h_Nm2,B程序首先將新的數(shù)據(jù)放置在xM處，然后對狀態(tài)緩存由下而上計算，在循環(huán)語句中每執(zhí)行一次狀態(tài)變量自動向下移動一級，即自動更新。它的計算量基本為 N個時鐘周期。當(dāng)然，數(shù)據(jù)存放還可以采用循環(huán)緩存。 另外，由于FIR的系數(shù)存在對稱性， 其結(jié)構(gòu)如圖3-所示。那么可以利用這個特點，實現(xiàn)更為快速的計算。x[n]h[0]yin]h[0]yin]圖3-對稱結(jié)構(gòu)的FIR為計算方便將狀態(tài)變量分別存放在兩個緩存區(qū)間內(nèi)，一塊命名為 Buffer.new,存放圖3-上半部分的數(shù)據(jù)，另一塊命名為Buffer_old ,存放圖3-下半部分的數(shù)據(jù)。它們都當(dāng)循環(huán)緩存使用，大小為FIR階數(shù)的一半，即N/2（常數(shù)SIZE） 。根據(jù)上述原理編寫的匯編程序如下:stm#Buffer_new,AR2stm#Buffer_o1d,AR3stm#SIZE,BKstm#-1,AROfir_loop::readinputId*(lnput),AstlA,*AR2;filteriadd*AR2+0%,*AR3+0%,ArptzB,#SIZE-1firs*AR2+0%>*AR3+0%,FIR_CoeffsthB,*(Output);storeoutputmar*+AR2(2)%mar*AR3+%mvdc*AR2,*AR3+0%:updatebufferbfir_loop為方便理解，在圖 3-中，將狀態(tài)變量更新的過程標(biāo)明，左邊的是 Buffer_new,右邊的是Buffer_old。開始時,指針AR2和AR3分別指向Buffer_newBuffer_old中的x[0]與x[-19],將最"新"的數(shù)據(jù)存進(jìn) Buffer_new（步驟1）。利用FIRS實現(xiàn)FIR,結(jié)果放在BH中。計算結(jié)束后，AR2和AR3分別指向x[-1]和x[-18]（步驟2）。然后調(diào)整AR2讓它指向Buffer_new41最"老"的數(shù)據(jù)x[-9[（步驟3）；調(diào)整AR3,讓它指向Buffer_old中最"老"的數(shù)據(jù)x[-19]（步驟4）。接下來進(jìn)行數(shù)據(jù)更新，將 Buffer_new中最〃老”的數(shù)據(jù)放進(jìn)Buffer_old中，成為最“新”的數(shù)據(jù)。最后AR2指向x[-9]的位置，這個位置將在下次計算時放置新的輸入;AR3扌旨向AR3扌旨向x[-18],即Buffer_old中最“老”的數(shù)據(jù)，便于下次計算（步驟5）ox[0]x[-19]驟5）ox[0]x[-19]x[?9]x[-10]X[-2]X[-1]X[-2]X[-1]x[-17]x[-18]圖3-對稱結(jié)構(gòu)的FIR實現(xiàn)中狀態(tài)變量的更新利用對稱結(jié)構(gòu)的實現(xiàn)在計算量上有了減少,大約為利用對稱結(jié)構(gòu)的實現(xiàn)在計算量上有了減少,大約為N/2個時鐘周期。當(dāng)然，利用上述結(jié)構(gòu)必須注意安排好數(shù)據(jù)的位置，以保證能進(jìn)行正常的循環(huán)尋址。3.IIR濾波器的實現(xiàn)IIR濾波器也是數(shù)字信號處理的主要工具之一，但由于它不具備線性相位， 而且無法準(zhǔn)確知道其延遲，所以使用較FIR少。下面，來觀察一下IIR的結(jié)構(gòu)，其數(shù)學(xué)表達(dá)式如下：TOC\o"1-5"\h\zM Ny[n]bix[ni]可y[ni]\o"CurrentDocument"i0 i1從其數(shù)學(xué)公式可以看出，我們可以仿照在FIR處理來直接實現(xiàn)。定義兩段緩存分別對應(yīng)于x[]和y[],然后采用類似于FIR的計算方式，采用MACD語句，便能很快地完成IIR濾波。在直接實現(xiàn)中同時需要保存 譏]和y[],當(dāng)其階數(shù)較大時，會占用比較大的數(shù)據(jù)空間。為此，我們來考察IIR的另一種結(jié)構(gòu)。如圖3-，在這種結(jié)構(gòu)中存儲的狀態(tài)變量為 d[],所以存儲空間大大地減少了。

x[n]d[n]b0圖3-IIR的另一種結(jié)構(gòu)x[n]d[n]b0圖3-IIR的另一種結(jié)構(gòu)y[n]另一方通過對FIR和IIR算法的分析，一方面向讀者介紹這些基木處理中的設(shè)計技巧,另一方以求利用它們獲得高效的面也提醒讀者在算法實現(xiàn)過程中應(yīng)充分考察算法本身的特點,運(yùn)算和存儲空間的節(jié)省。以求利用它們獲得高效的4.FFT的實現(xiàn)在信號處理中，為突出信號的特征，經(jīng)常在時域和頻域之間作變換， 而傅立葉變換是這當(dāng)中最為典型的?；?的快速傅立葉變換有比特翻轉(zhuǎn)排序和蝶形運(yùn)算組成， 前者在3.x節(jié)己經(jīng)作了介紹，這里著重介紹蝶形運(yùn)算的實現(xiàn)。蝶形運(yùn)算是快速傅立葉變換中設(shè)計得極為精巧的部分， 它充分揭示了傅立葉變換系數(shù)間內(nèi)在的關(guān)系，而且還能實現(xiàn)同址計算，如圖所示。完成一次蝶形運(yùn)算需要一次復(fù)數(shù)的乘法和兩次復(fù)數(shù)的加法。Xm(P) * 1Xm+i(p)xm(q) * xm+1(q)Wn 1圖3-蝶形運(yùn)算的示意圖

對于時間抽取的FFT而言，在其第一級是乘法的系數(shù)為WN（也就是1）,那么這個復(fù)數(shù)的乘法就名存實亡了，因而在計算FFT時，第一級可以直接用復(fù)數(shù)的加法實現(xiàn)。第-級的程序如下:；***stage1***stm#0,BKIdti-1,ASMstm#FFT_Data,PXstm#FFT_Data+K_DATA_IDX_l,QX

Idstm*PX,16,A;AH二PX?xrptbdstm#K_FFT_SIZE/2T,BRCstagelend-1#KDATAIDX1+1,ARO*QX,16,A,B;BH二PX.x-QX.x*QX,16,A;AH二PX.x+QX.Idstm*PX,16,A;AH二PX?xrptbdstm#K_FFT_SIZE/2T,BRCstagelend-1#KDATAIDX1+1,ARO*QX,16,A,B;BH二PX.x-QX.x*QX,16,A;AH二PX.x+QX.xA,ASM,*PX+B,*QX+*PX,A;AH二PX.y*QX,16,A,B;BH二PX.v-QX.ysubadd*QX,16,A;AH二PX.y+QX.ysthst IdASM,*PX+0*QX+O%subadd*PX,Asthst Idstagelend:在代碼中，為方便與圖3-對應(yīng)，使用偽指令將寄存器的名字替換成示意圖中的表達(dá)方式，PX和QX分別指向蝶形運(yùn)算的數(shù)據(jù)的地址。 可見所有的操作完全是由加減實現(xiàn)的。在第二級中乘法的系數(shù)為wN和wN'1（即+1和-J）,那么乘法操作只是影響參數(shù)的符號，在復(fù)數(shù)的加減運(yùn)算時只要弄清與-J相乘造成的結(jié)果，所有的操作仍然可以只用加減法來實現(xiàn)。第二級的實現(xiàn)代碼如下：stmstm#FFUata,PX#FFT_Data+K_DATA_IDX_2,QX*PX,16,A;AH=PX.xIdstm#K_FFT_SIZE/4-l,BRCstage2end-lrptbdstm#K_DATA_IDX_2+1,ARO;1stbutterfly*QX,16,A,B;BH二PX.x-QX.x;AH二PX.x+QX.xsubadd*QX,16,AASM,*QX+*PX+sthst1IId*PX,A;AH二PX.y*QX,16,A,B;BH二PX.y-QX.ysub*QX,16,A;AH二PX.y+QX.yaddsthA,ASM,*PX+sthB,ASM,*QX+;2ndbutterfly

mar＞：＜QX+add*PX,＞：＜QX,A;AH二PX.x+QX.ysub*PX,＞：＜QX-,B;BH二PX.x-QX.ysthA,ASM,*PX+sub*PX,＞：＜QX,A;AH二PX.y-QX.xstB,＞：＜QX11Id＞：＜QX+,B;BH二QX.xStA,*PXadd*PX+0%,A;AH二PX.y+QX.xStA,＞：＜QX+O%11Id*PX,Astage2end:第二級與第一級相同，計算中不包含乘法。從第三級開始，乘法系數(shù)的特征就沒有第第二級那樣好了，所以只能直接采用圖3-的方法來計算，代碼如下。;***Stage3throughStagelogN***stm#K_TWID_TBL_SIZE,BKst#K_TWID_IDX_3,*(d_twid_idx)stmffK_TWID_IDX_3,AROstm#Cos,WRstm#Sin,WIstm#K_LOGN-2-l,STAGE_COUNTERst#K_FFT_SIZE/8-1,*(d_grps_cnt)stm#K_FLY_COUNT_3-1,BUTTERFLY.COUNTERst#K_DATA_IDX_3,*(d_data_idx)stage:stm#FFT_Data,PXId*(d_data_idx),Aadd*(PX),AstlmA,QXmvdk*(d_grps_cnt),GR0UP_C0UNTERgroup:mvmdBUTTERFLY.COUNTER,BRCrptbdIdmpybutterflyend-1*WR,T*QX+,A;T:二WR;A:二QR＞：＜WR QX-＞QImacr*WI+0%,*QX-,A;A:二QR＞：＜WR+Q"WI9QX-＞QRadd*PX,16,A,B;B:二(QR＞：＜WR+QI*WI)+PRstB,*PX;PR'：二((QR*WR+QI*WI)+PR)/2sub*PX+,B;B二PR-(QR*WR+QI*WI)；:PX-＞PIstB,*QX;QR'：二(PR-(QR*WR+QI*WI))/2mpy*QX+,A:A：=QR*WI[T=WI]

QX->QIIIaddstsub*PX,16,A,B ;B:=(QR*WI-QI*WR)+PI;QI1:=((QR*WI-QI*WR)+PI)/2B,水QX+*PX,B;B=PI-(QR*WI-QI*WR)Id水WR,T;T：=WRstB,*PX+;Pf：=(PI-(QR*WI-QI*WR))/2IImpy水QX+,A:1IPX->PR;A：=QR*WR QX->QIbutterflyend：;UpdatepointersfornextgrouppshmAROmvdk*(d_data_idx),AROmar*PX+Omar*QX+Obanzdgroup,*GROUP_COUNTER-popmAROmar*QX-;Updatecountersandindicesf