DSP常见算法的实现

1、3.6常见的算法实现在实际应用中虽然信号处理的方式多种多样，但其算法的基本要素却大多相同，在本节 中介绍几种较为典型的算法实现，希望通过对这些例子(单精度，16bit)的分析，能够让 大家熟悉DSP编程中的一些技巧，在以后的工作中可以借鉴，达到举一反三的效果。1. 函数的产生在高级语言的编程中，如果要使用诸如正弦、余弦、对数等数学函数，都可以直接调用 运行库中的函数来实现，而在DSP编程中操作就不会这样简单了。虽然TI公司提供的实时 运行库中有一些数学函数，但它们所耗费的时间大多太长，而且对于大多数定点程序使用双 精度浮点数的返回结果有点“大材小用”的感觉，因此需要编程人员根据自身的要求“定制

2、” 数学函数。实现数学函数的方法主要有查表法、迭代法和级数逼近法等，它们各有特点，适 合于不同的应用。查表法是最直接的一种方法，程序员可以根据运算的需要预先计算好所有可能出现的函 数值，将这些结果编排成数据表，在使用时只需要根据输入查出表中对应的函数值即可。它 的特点是速度快，但需要占用大量的存储空间，且灵活度低。当然，可以对上述查表法作些 变通，仅仅将一些关键的函数值放置在表中，对任意一个输入，可根据和它最接近的数据采 用插值方法来求得。这样占用的存储空间有所节约，但数值的准确度有所下降。迭代法是一种非常有用的方法，在自适应信号处理中发挥着重要的作用。作为函数产生 的一种方法，它利用了自变量

3、取值临近的函数值之间存在的关系，如时间序列分析中的AR、 MA、ARMA等模型，刻画出了信号内部的特征。因为它只需要存储信号模型的参量和相关 的状态变量，所以所占用的存储空间相对较少，运算时间也较短。但它存在一个致命的弱点， 由于新的数值的产生利用了之前的函数值，所以它容易产生误差累积，适合精度要求不高的 场合。级数逼近法是用级数的方法在某一自变量取值范闱内去追近数学函数，而将自变量取值 在此范围外的函数值利用一些数学关系，用该范围内的数值来表示。这种方法最大的优点是 灵活度高，且不存在误差累积，数值精度由程序员完全控制。该方法的关键在于选择一个合 适的自变量取值区间和寻找相应的系数。下而通过

4、正弦函数的实现，具体对上述三种方法作比较。查表法较简单，只需要自制一张数据表，也可以利用C5400 DSP ROM内的正弦函数表。迭代法的关键是寻找函数值间的递推关系。假设函数采样时间间隔为T,正弦函数的角 频率为e,那么可以如下推导：令 siiMe + tyT)： asin(°)+ /7sin(-<y7')等式的左边展开为left_ side = sin (pcoscoT + cossin coT等式的右边展开为right _ side = s in(p(a + cos coT )-/7 cos 夕 s in coT对比系数，可以得到c = 2cosa；/7 = -l

5、。令(p = nT,便可以得到如下的递推式：sii = 2coscoT-1-sn - 2令sl = O,s-2 = _4sin，逐一迭代就能够获得采样间隔为T的正弦序列了。从迭 代公式可以更清楚地看出，这种方法存在误差累积。再来看级数逼近法，首先需要寻找一个合适的自变量取值区间和寻找相应的系数。从正弦函 数的对称性可知，只需要计算取值在0,笈内的函数就可以推断出所有取值范围内的函数。接下来寻找系数，对于sin(x)可以作如下变换sin(x) = sin(金)= sin_new(y),那么y的取值范围在0,1内，而sin_new( sin()同构，所以在下面的分析中将sin_ncw(潜代sin(

6、),0x0000OxFFFF提到的正弦函数即指sin_new()。若汇编编程时采用的数据为Q15格式，那么取值与实际的 弧度的对应关系如下图所示。0X7FFFn0x8000 -图3-算法取值与弧度的对应关系在0,1内，正弦函数的修正级数(五次)展开如下式:sin_new(x)= 3.140625x+ 0.02026367c2 -5.325196x3 +0.5446778x4 +1.800293x5根据上式，可以写出正弦函数的生成程序。:compute polynomialstlmA.T；T=xstm#SinePolyCoeff, AR2Id*AR2+, 16, A;AH=c5Id*AR2+,

7、16,B;BH=c4poiy*AR2+;A=c5*x+c4psy*AR2+;A=c5*xA2+c4*x+c3psy*AR2+;A=c5*xA3+c4*xA2+c3*x+c2poly*AR2+;A=c5*xA4+c4*xA3+c3*xA2+c2*x+clmpyaAsftaA, 3;adjust AH to Q15SinePolyCoeff:.word Oxlcce.word 0x08b7;in QI2 format;1.800293 (coef for xA5 = c5);0.5446778 (coef for xA4 = c4).word.word.wordOxaacc;-5.325196(c

8、oef for xA3 = c3)0x0053;0.02026367(coef for xA2 = c2)0x3240;3.140625(coef for xAl = cl)在编程过程中，使用到了 POLY语句，它能够使多项式的计算简洁快速地完成。该函数的结 果可以通过实验X来验证。2. FIR滤波器的实现FIR滤波器由于具有线性相位而且延迟能够确定，因而在信号处理中广泛应用。FIR的 基本模型如下图所示.图3- FIR模型NT其数学表达式为»巾N一,I，根据该表达式可以给出一种最为直接的实现形 i=0式。直接形式中采用线性地址来存放数据，如图3-所示。xnhN-1 xn-1 h N

9、-2 i，x n - ( N - 2)h Ux n - ( N -1)h 0图3-直接形式程序中可以采用MACD来实现程序如下：Id"(Input), AstlA, *(x_n)stm#x n Nm1, AR2mpy *AR2-,#h_Nm1, Brpt#N-2macd “AR2-, h_Nm2, B程序首先将新的数据放置在xn处，然后对状态缓存由下而上计算，在循环语句中每执 行一次状态变量自动向下移动一级，即自动更新。它的计算量基本为N个时钟周期。当然，数据存放还可以采用循环缓存。另外，由于FIR的系数存在对称性，其结构如 图3-所示。那么可以利用这个特点，实现更为快速的计算。st

10、m stm stm stm fir Joop:Idstl"(Input), AA. *AR2:filtering addrptzfirssth*AR2+0%, *AR3+0%. AB, #SIZE-1*AR2+0%, *AR3+0%, HR.CocffB, *(Output);store outputmarmarmvdd b为方便理解,*+AR2(2)%*AR3+%*AR2, *AR3+0% ;update bufferfir_loop在图3-中,将状态变量更新的过程标明,左边的是Buffer.new,右边的是BuffejokL开始时，指针AR2和AR3分别指向Buffejnew和B

11、uffer.old中的x与x-19, 将最“新”的数据存进Buffejnew （步骤1 ）。利用FIRS实现FIR,结果放在BH中。计算 结束后，AR2和AR3分别指向x-l和x-18（步骤2）o然后调整AR2,让它指向Buffer.new 中最“老”的数据x-9（步骤3）；调整AR3,让它指向Buffejold中最“老”的数据x-19图3-对称结构的FIR为计算方便将状态变量分别存放在两个缓存区间内，一块命名为Buffejnew,存放图3-上 半部分的数据，另一块命名为Buffer_old,存放图3-下半部分的数据。它们都当循环缓存 使用，大小为FIR阶数的一半，即N/2 （常数SIZE）。

12、根据上述原理编写的汇编程序如下：#Buffer_new, AR2#Buffer_old. AR3#SIZE. BK#-l, ARO;read input（步骤4）。接下来进行数据更新，将Buffejnew中最“老”的数据放进Buffejold中，成 为最“新”的数据。最后AR2指向x-9的位置，这个位置将在下次计算时放置新的输入： AR3指向即Buffejold中最“老”的数据，便于下次计算（步骤5）。x0x-19©x-9X-1Ox-2x-17®x-1x-18图3-对称结构的FIR实现中状态变量的更新利用对称结构的实现在计算量上有了减少，大约为N/2个时钟周期。当然，利用上

13、述结构必 须注意安排好数据的位置，以保证能进行正常的循环寻址。3. HR滤波器的实现HR滤波器也是数字信号处理的主要工具之一，但由于它不具备线性相位，而且无法准 确知道其延迟，所以使用较FIR少。下而，来观察一下HR的结构，其数学表达式如下:MNyW = 一 日 + 一 a /=0/=1从其数学公式可以看出，我们可以仿照在FIR处理来直接实现。定义两段缓存分别对 应于x和y,然后采用类似于FIR的计算方式,采用MACD语句,便能很快地完成 IIR滤波。在直接实现中同时需要保存x和y,当其阶数较大时，会占用比较大的数据空间。 为此，我们来考察HR的另一种结构。如图3-,在这种结构中存储的状态变量

14、为d, 所以存储空间大大地减少了。图3-IIR的另一种结构通过对FIR和HR算法的分析，一方面向读者介绍这些基本处理中的设计技巧，另一方 而也提醒读者在算法实现过程中应充分考察算法本身的特点，以求利用它们获得高效的 运算和存储空间的节省。4. FFT的实现在信号处理中，为突出信号的特征，经常在时域和频域之间作变换，而傅立叶变换是这 当中最为典型的。基2的快速傅立叶变换有比特翻转排序和蝶形运算组成，前者在3.x 节已经作了介绍，这里着重介绍蝶形运算的实现。蝶形运算是快速傅立叶变换中设计得极为精巧的部分，它充分揭示了傅立叶变换系数间 内在的关系，而且还能实现同址冲算，如图所示。完成一次蝶形运算需要

15、一次复数的乘 法和两次复数的加法.图3-蝶形运算的示意图对于时间抽取的FFT而言，在其第一级是乘法的系数为1匕§ （也就是1）,那么这个复 数的乘法就名存实亡了，因而在计算FFT时，第一级可以直接用复数的加法实现。第 一级的程序如下：;* stage 1 *stm#0. BKIdstm#FFT_Data, PXstm#FFT Data+K DATAIDX 1,QXId*PX, 16,A:AH=PX.xstm#K_FFT_SIZE/2-l, BRCrptbdstage 1 end-1stm#K_DATAJDX_1+LAROsub add sth stII Id sub add sth

16、stII Id stage lend: 在代码中,*QX,16,A.B *QX, 16,A A, ASM, *PX+ B. *QX+ *PX,A *QX,16,A.B *QX, 16,A A. ASM, *PX+O B. *QX+O% *PX.A:BH=PX,X-QX.x ;AH=PX.x+QX.x:AH=PX.y:BH=PX,y-QX.y:AH=PX.y+QX.y为方便与图3-对应,使用.asg伪指令将寄存器的名字替换成示意图中的表达方式，PX和QX分别指向蝶形运算的数据的地址。可见所有的操作完全是由加减实 现的。在第二级中乘法的系数为卬1和1明/4 （即+i和小），那么乘法操作只是影响参数

17、的符 号，在复数的加减运算时只要弄清与-j相乘造成的结果，所有的操作仍然可以只用加减 法来实现。第二级的实现代码如下：产* Stage 2 *stm#FFT_Data, PXstm#FFT Data+K DATA IDX 2. QXId*PX, 16,A:AH=PX.xstm#K_FFT_SIZE/4-l, BRCrptbdstage2end-lstm#K_DATA_IDX_2+ L ARO;1st butterfly:BH=PX.X-QX.x;AH=PX.x+QX.x:AH=PX.y:BH=PX,y-QX.y:AH=PX.y+QX.ysub*QX, 16,A, Badd*QX, 16,Ast

18、hA, ASM, *PX+stB. *QX+II Id*PX, Asub*QX, 16,A, Badd*QX, 16,AsthA, ASM. *PX+sthB.ASM. *QX+;2nd butterflymar*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.*QXIIId*QX+. B:BH=QX.xstA, *PXIIadd*PX+0%, A;AH=PX.y+QX.xstA, *QX+0%IIId*PX,Astage2end:第二

19、级与第一级相同，计算中不包含乘法。从第三级开始，乘法系数的特征就没有第一、 第二级那样好了，所以只能直接采用图3-的方法来计算，代码如下。产* Stage 3 through Stage logN *stm#K_TWID_TBL_SIZE. BKst#K_TWID_IDX_3, *(djwid_idx)stm#K_TWID_IDX_3. AROstm#Cos, WRstm#Sin, WIstm#K_LOGN-2-k STAGE.COUNTERst#K_FFT_SIZE/8-l, *(d_grps_cnt)stm#K_FLY_COUNT_3-1, BUTTERFLY.COUNTERst#K_DA

20、TA_IDX_3, *(d_datajdx)stage:stm Id add stlm mvdk group:mvmd rptbd#FFT_Data, PX*(d_datajdx), A*(PX), AA,QX*(d_grps_cnt), GROUP.COUNTERBUTTERFLY.COUNTER, BRC butterflyend-1Idmpy*WR.T*QX+, A；T := WR;A := QR*WR II QX->QImacr*WI+O%, *QX-, A；A := QR*WR+QI*WI：ll QX->QRadd*PX, 16.A,B;B := (QR*WR+QI*WI

21、)+PRstB,*PX:PR':=(QR*WR+QI*WI)+PR)/2IIsub*PX+, B;B=PR-(QR*WR+QI*WI);II PX->PIstB,*QX;QR':= (PR-(QR*WR+QI*WI)/2IImpy*QX+, A；A := QR*WI T=WI;ll QX->QImasr*QX, *WR+0%. A；A := QR*WI-QI*WRadd*PX, 16,A.B；B :=(QR*WI-QI*WR)+PIstB, *QX+;QF:=(QR*WI-QI*WR)+PI)/2II sub*PX,B;B=PI-(QR*WI-QOWR)Id*WR.T；T := WRstB. *PX+:pr：= (PL(QR*WI-QI*WR)/2II mpy*QX+. A；ll PX->PR；A := QR*WR II QX->QIbutterflyend:; Update pointers for next grouppshmAROmvdk*(d_datajdx), AROmar*PX+0mar*QX+0banzdgroup, *GROUP_COUNTER-popinAROmar*QX-; Update counters and