DSP常见算法的实现.doc

3.6 常见的算法实现在实际应用中虽然信号处理的方式多种多样，但其算法的基本要素却大多相同，在本节中介绍几种较为典型的算法实现，希望通过对这些例子（单精度，16bit）的分析，能够让大家熟悉DSP编程中的一些技巧，在以后的工作中可以借鉴，达到举一反三的效果。1. 函数的产生在高级语言的编程中，如果要使用诸如正弦、余弦、对数等数学函数，都可以直接调用运行库中的函数来实现，而在DSP编程中操作就不会这样简单了。虽然TI公司提供的实时运行库中有一些数学函数，但它们所耗费的时间大多太长，而且对于大多数定点程序使用双精度浮点数的返回结果有点“大材小用”的感觉，因此需要编程人员根据自身的要求“定制”数学函数。实现数学函数的方法主要有查表法、迭代法和级数逼近法等，它们各有特点，适合于不同的应用。查表法是最直接的一种方法，程序员可以根据运算的需要预先计算好所有可能出现的函数值，将这些结果编排成数据表，在使用时只需要根据输入查出表中对应的函数值即可。它的特点是速度快，但需要占用大量的存储空间，且灵活度低。当然，可以对上述查表法作些变通，仅仅将一些关键的函数值放置在表中，对任意一个输入，可根据和它最接近的数据采用插值方法来求得。这样占用的存储空间有所节约，但数值的准确度有所下降。迭代法是一种非常有用的方法，在自适应信号处理中发挥着重要的作用。作为函数产生的一种方法，它利用了自变量取值临近的函数值之间存在的关系，如时间序列分析中的AR、MA、ARMA等模型，刻画出了信号内部的特征。因为它只需要存储信号模型的参量和相关的状态变量，所以所占用的存储空间相对较少，运算时间也较短。但它存在一个致命的弱点，由于新的数值的产生利用了之前的函数值，所以它容易产生误差累积，适合精度要求不高的场合。级数逼近法是用级数的方法在某一自变量取值范围内去逼近数学函数，而将自变量取值在此范围外的函数值利用一些数学关系，用该范围内的数值来表示。这种方法最大的优点是灵活度高，且不存在误差累积，数值精度由程序员完全控制。该方法的关键在于选择一个合适的自变量取值区间和寻找相应的系数。下面通过正弦函数的实现，具体对上述三种方法作比较。查表法较简单，只需要自制一张数据表，也可以利用C5400 DSP ROM内的正弦函数表。迭代法的关键是寻找函数值间的递推关系。假设函数采样时间间隔为T，正弦函数的角频率为，那么可以如下推导：令等式的左边展开为等式的右边展开为对比系数，可以得到。令，便可以得到如下的递推式：令，逐一迭代就能够获得采样间隔为T的正弦序列了。从迭代公式可以更清楚地看出，这种方法存在误差累积。再来看级数逼近法，首先需要寻找一个合适的自变量取值区间和寻找相应的系数。从正弦函数的对称性可知，只需要计算取值在内的函数就可以推断出所有取值范围内的函数。接下来寻找系数，对于可以作如下变换，那么y的取值范围在内，而sin_new( )与sin( )同构，所以在下面的分析中将sin_new( )替代sin( )，提到的正弦函数即指sin_new( )。若汇编编程时采用的数据为Q15格式，那么取值与实际的弧度的对应关系如下图所示。图3- 算法取值与弧度的对应关系在内，正弦函数的修正级数（五次）展开如下式：根据上式，可以写出正弦函数的生成程序。;compute polynomialstlmA, T;T=xstm#SinePolyCoeff, AR2ld*AR2+, 16, A;AH=c5ld*AR2+, 16, B;BH=c4poly*AR2+;A=c5*x+c4poly*AR2+;A=c5*x2+c4*x+c3poly*AR2+;A=c5*x3+c4*x2+c3*x+c2poly*AR2+;A=c5*x4+c4*x3+c3*x2+c2*x+c1mpyaAsftaA, 3;adjust AH to Q15SinePolyCoeff:;in Q12 format.word0x1cce; 1.800293(coef for x5 = c5).word0x08b7; 0.5446778(coef for x4 = c4).word0xaacc; -5.325196(coef for x3 = c3).word0x0053; 0.02026367(coef for x2 = c2).word0x3240; 3.140625(coef for x1 = c1)在编程过程中，使用到了POLY语句，它能够使多项式的计算简洁快速地完成。该函数的结果可以通过实验X来验证。2. FIR滤波器的实现FIR滤波器由于具有线性相位而且延迟能够确定，因而在信号处理中广泛应用。FIR的基本模型如下图所示。图3- FIR模型其数学表达式为，根据该表达式可以给出一种最为直接的实现形式。直接形式中采用线性地址来存放数据，如图3- 所示。图3- 直接形式程序中可以采用MACD来实现程序如下：ld*(Input), AstlA, *(x_n)stm#x_n_Nm1, AR2mpy*AR2-, #h_Nm1, Brpt#N-2macd*AR2-, h_Nm2, B程序首先将新的数据放置在xn处，然后对状态缓存由下而上计算，在循环语句中每执行一次状态变量自动向下移动一级，即自动更新。它的计算量基本为N个时钟周期。当然，数据存放还可以采用循环缓存。另外，由于FIR的系数存在对称性，其结构如图3- 所示。那么可以利用这个特点，实现更为快速的计算。图3- 对称结构的FIR为计算方便将状态变量分别存放在两个缓存区间内，一块命名为Buffer_new，存放图3- 上半部分的数据，另一块命名为Buffer_old，存放图3- 下半部分的数据。它们都当循环缓存使用，大小为FIR阶数的一半，即N/2（常数SIZE）。根据上述原理编写的汇编程序如下：stm#Buffer_new, AR2stm#Buffer_old, AR3stm#SIZE, BKstm# -1, AR0fir_loop:;read inputld*(Input), AstlA, *AR2;filteringadd*AR2+0%, *AR3+0%, ArptzB, #SIZE-1firs*AR2+0%, *AR3+0%, FIR_CoeffsthB, *(Output);store outputmar*+AR2(2)%mar*AR3+%mvdd*AR2, *AR3+0%;update bufferbfir_loop为方便理解，在图3- 中，将状态变量更新的过程标明，左边的是Buffer_new，右边的是Buffer_old。开始时，指针AR2和AR3分别指向Buffer_new和Buffer_old中的x0与x-19，将最“新”的数据存进Buffer_new（步骤1）。利用FIRS实现FIR，结果放在BH中。计算结束后，AR2和AR3分别指向x-1和x-18（步骤2）。然后调整AR2，让它指向Buffer_new中最“老”的数据x-9 （步骤3）；调整AR3，让它指向Buffer_old中最“老”的数据x-19（步骤4）。接下来进行数据更新，将Buffer_new中最“老”的数据放进Buffer_old中，成为最“新”的数据。最后AR2指向x-9的位置，这个位置将在下次计算时放置新的输入；AR3指向x-18，即Buffer_old中最“老”的数据，便于下次计算（步骤5）。图3- 对称结构的FIR实现中状态变量的更新利用对称结构的实现在计算量上有了减少，大约为N/2个时钟周期。当然，利用上述结构必须注意安排好数据的位置，以保证能进行正常的循环寻址。3. IIR滤波器的实现IIR滤波器也是数字信号处理的主要工具之一，但由于它不具备线性相位，而且无法准确知道其延迟，所以使用较FIR少。下面，来观察一下IIR的结构，其数学表达式如下：从其数学公式可以看出，我们可以仿照在FIR处理来直接实现。定义两段缓存分别对应于x 和y ，然后采用类似于FIR的计算方式，采用MACD语句，便能很快地完成IIR滤波。在直接实现中同时需要保存x 和y ，当其阶数较大时，会占用比较大的数据空间。为此，我们来考察IIR的另一种结构。如图3- ，在这种结构中存储的状态变量为d ，所以存储空间大大地减少了。图3- IIR的另一种结构通过对FIR和IIR算法的分析，一方面向读者介绍这些基本处理中的设计技巧，另一方面也提醒读者在算法实现过程中应充分考察算法本身的特点，以求利用它们获得高效的运算和存储空间的节省。4. FFT的实现在信号处理中，为突出信号的特征，经常在时域和频域之间作变换，而傅立叶变换是这当中最为典型的。基2的快速傅立叶变换有比特翻转排序和蝶形运算组成，前者在3.x节已经作了介绍，这里着重介绍蝶形运算的实现。蝶形运算是快速傅立叶变换中设计得极为精巧的部分，它充分揭示了傅立叶变换系数间内在的关系，而且还能实现同址计算，如图所示。完成一次蝶形运算需要一次复数的乘法和两次复数的加法。图3- 蝶形运算的示意图对于时间抽取的FFT而言，在其第一级是乘法的系数为（也就是1），那么这个复数的乘法就名存实亡了，因而在计算FFT时，第一级可以直接用复数的加法实现。第一级的程序如下：;* stage 1 * stm#0, BK ld#-1, ASM stm#FFT_Data, PX stm#FFT_Data+K_DATA_IDX_1, QX ld*PX, 16, A;AH=PX.x stm#K_FFT_SIZE/2-1, BRC rptbdstage1end-1 stm #K_DATA_IDX_1+1, AR0 sub*QX, 16, A, B;BH=PX.x-QX.x add*QX, 16, A;AH=PX.x+QX.x sthA, ASM, *PX+ stB, *QX+|ld*PX, A;AH=PX.y sub*QX, 16, A, B;BH=PX.y-QX.y add*QX, 16, A;AH=PX.y+QX.y sthA, ASM, *PX+0 stB, *QX+0%|ld*PX, Astage1end:在代码中，为方便与图3- 对应，使用.asg伪指令将寄存器的名字替换成示意图中的表达方式，PX和QX分别指向蝶形运算的数据的地址。可见所有的操作完全是由加减实现的。在第二级中乘法的系数为和（即+1和-j），那么乘法操作只是影响参数的符号，在复数的加减运算时只要弄清与-j相乘造成的结果，所有的操作仍然可以只用加减法来实现。第二级的实现代码如下：;* Stage 2 * stm#FFT_Data, PX stm#FFT_Data+K_DATA_IDX_2, QX ld*PX, 16, A;AH=PX.x stm#K_FFT_SIZE/4-1, BRC rptbdstage2end-1 stm#K_DATA_IDX_2+1, AR0; 1st butterfly sub*QX, 16, A, B;BH=PX.x-QX.x add*QX, 16, A;AH=PX.x+QX.x sthA, ASM, *PX+ stB, *QX+|ld*PX, A;AH=PX.y sub*QX, 16, A, B;BH=PX.y-QX.y add*QX, 16, A;AH=PX.y+QX.y sthA, ASM, *PX+ sthB, ASM, *QX+; 2nd butterfly mar*QX+ add*PX, *QX, A;AH=PX.x+QX.y sub*PX, *QX-, B;BH=PX.x-QX.y sthA, ASM, *PX+ sub*PX, *QX, A;AH=PX.y-QX.x stB, *QX|ld*QX+, B;BH=QX.x stA, *PX|add*PX+0%, A;AH=PX.y+QX.x stA, *QX+0%|ld*PX, Astage2end:第二级与第一级相同，计算中不包含乘法。从第三级开始，乘法系数的特征就没有第一、第二级那样好了，所以只能直接采用图3- 的方法来计算，代码如下。;* Stage 3 through Stage logN * stm#K_TWID_TBL_SIZE, BK st#K_TWID_IDX_3, *(d_twid_idx) stm#K_TWID_IDX_3, AR0 stm#Cos, WR stm#Sin, WI stm#K_LOGN-2-1, STAGE_COUNTER st#K_FFT_SIZE/8-1, *(d_grps_cnt) stm#K_FLY_COUNT_3-1, BUTTERFLY_COUNTER st#K_DATA_IDX_3, *(d_data_idx)stage: stm#FFT_Data, PX ld*(d_data_idx), A add*(PX), A stlmA, QX mvdk*(d_grps_cnt), GROUP_COUNTERgroup: mvmdBUTTERFLY_COUNTER, BRC rptbdbutterflyend-1 ld*WR, T; T := WR mpy*QX+, A ; A := QR*WR | QX-QI macr*WI+0%, *QX-, A ; A := QR*WR+QI*WI ; | QX-QR add*PX, 16, A, B ; B := (QR*WR+QI*WI)+PR stB, *PX ; PR:=(QR*WR+QI*WI)+PR)/2|sub*PX+, B; B=PR-(QR*WR+QI*WI); | PX-PI stB, *QX ; QR:= (PR-(QR*WR+QI*WI)/2|mpy*QX+, A ; A := QR*WI T=WI ; | QX-QI masr*QX, *WR+0%, A ; A := QR*WI-QI*WR add*PX, 16, A, B ; B := (QR*WI-QI*WR)+PI stB, *QX+ ; QI:=(QR*WI-QI*WR)+PI)/2|sub*PX, B; B=PI-(QR*WI-QI*WR) ld*WR, T ; T := WR stB, *PX+ ; PI:= (PI-(QR*WI-QI*WR)/2|mpy*QX+, A ; | PX-PR ; A := QR*WR | QX-QIbutterflyend:; Update pointers for next group pshmAR0 mvdk*(d_data_idx), AR0 mar*PX+0 mar*QX+0 banzdgroup, *GROUP_COUNTER- popmAR0 mar*QX-; Update counters and indices for next stage ld*(d_data_idx), A sub#1, A