DSP常見算法的實現

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

令s1 0,s2AsinT,逐一迭代就能夠獲得采樣間隔為 t的正弦序列了。從迭代公式可以更清楚地看岀，這種方法存在誤差累積。再來看級數逼近法，首先需要尋找一個合適的自變量取值區(qū)間和尋找相應的系數。 從正弦函數的對稱性可知，只需要計算取值在 ［0,］內的函數就可以推斷出所有取值范圍內的函數。接下來尋找系數，對于sinx可以作如下變換sinxsiny值范圍在［0,1］內，而sin.new()與sin()同構，所以在下面的分析中將sin(),提到的正弦函數即指 sin_new()o若匯編編程時采用的數據為sin_new(y),那么y的取sin_new()sin_new(y),那么y的取sin_new()替代Q15格式,那么取2圖3-算法取值與弧度的對應關系在［0,1］內，正弦函數的修正級數(五次)展開如下式：sin_new(x)3.140625x0.02026367k2-5.325196x3 0.5446778xl 1.800293x5根據上式，可以寫岀正弦函數的生成程序。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語句，它能夠使多項式的計算簡潔快速地完成。該函數的結果可以通過實驗X來驗證。2.FIR濾波器的實現FIR的FIR濾波器由于具有線性相位而且延遲能夠確定，因而在信號處理中廣泛應用。FIR的x(n)h(0)z1z1h(1)?dy(n)x(n)h(0)z1z1h(1)?dy(n)基木模型如下圖所示。圖3-FIR基木模型如下圖所示。圖3-FIR模型N1其數學表達式為yn hii0ni,根據該表達式可以給出一種最為直接的實現形式。直接形式中采用線性地址來存放數據，如圖3-所示。式。直接形式中采用線性地址來存放數據，如圖3-所示。程序中可以采用MACD來實現程序如下:Id*(Input),AstlA,*(x_n)stm#x_n_Nml,AR2mpy*AR2~,#h_Nml,Brpt#N—2

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

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

對于時間抽取的FFT而言，在其第一級是乘法的系數為WN（也就是1）,那么這個復數的乘法就名存實亡了，因而在計算FFT時，第一級可以直接用復數的加法實現。第-級的程序如下:；***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-對應，使用偽指令將寄存器的名字替換成示意圖中的表達方式，PX和QX分別指向蝶形運算的數據的地址。 可見所有的操作完全是由加減實現的。在第二級中乘法的系數為wN和wN'1（即+1和-J）,那么乘法操作只是影響參數的符號，在復數的加減運算時只要弄清與-J相乘造成的結果，所有的操作仍然可以只用加減法來實現。第二級的實現代碼如下：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:第二級與第一級相同，計算中不包含乘法。從第三級開始，乘法系數的特征就沒有第第二級那樣好了，所以只能直接采用圖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