定点DSP中高精度除法的实现方法

引言

各种集成化单片数字信号处理器(DSP)以其功能强、集成度高、应用灵活、性价比高等优点，在信号处理和系统控制中的主导性地位日益明显。许多信号处理和控制需要运用除法运算。一般的数字信号处理器中没有现成的除法指令。十多年前诞生的浮点DSP，由于其用硬件完成浮点数的运算，在数据处理和运算能力上大大超出定点DSP，处理除法运算也比定点DSP更为简单。但是定点DSP每器件产品的价格更低，这对大规模的大众市场应用而言是相当重要的优势，也是定点器件至今仍是业界主流的主要原因。所以，讨论定点DSP中除法的实现仍不失其意义。在定点DSP中虽然已经有人给出除法的算法，但是由于其运算的复杂和精度难以如愿，致使一些带有除法的好的算法在信号处理中难以得到应用。

为了提高运算结果的精度，本文在已有除法算法的基础上进行了一些改进，最大限度地保证了结果的精度。最后在TI公司的TMS320C5416芯片里具体实现并验证了这一高精度除法。

1 经典算法

DSP中没有现成的除法指令，除法是靠被除数与除数之间的移位相减来实现的。在C54X系列里利用减法指令SUBC和循环指令RPT实现2个16位数的相除。下面以C54X为例来具体实现经典的除法：

C54X提供的SuBC指令仅对无符号数进行操作，所以在移位相减开始之前必须先将被除数和除数取绝对值，仅考虑2个正数的除法。此时除法运算有两种情况：

当|被除数|<|除数|时，将|被除数|存放在累加器的高16位，然后用SUBC完成15次移位相减，相减之后在累加器A的低16位中存放商的绝对值。根据运算前被除数和除数的符号是否相同来决定是否要改变所得结果的符号。

当|被除数|≥|除数|时，将|被除数|存放在累加器的低16位，然后用SUBC完成16次移位相减，相减之后在累加器A的低16位中存放商的绝对值。根据运算前被除数和除数的符号是否相同来决定是否要改变所得结果的符号。

从实现的过程分析，当|被除数|<|除数|时，移位相减开始时|被除数|和|除数|的小数点位置正好相差一位。第一次相减后在累加器A的O位最低位存进的数值正是商的最高位，该位为商的小数点后第一位。在15次移位相减之后，累加器A低16位所得的结果为Q值为15的小数。当|被除数|≥|除数|时，在第l6次相减时，|被除数|位于A的高16位(30～15位)上，小数点位在A的15位后，和|除数|的小数点位正好对齐，则此次相减后在A的 0位加上的值正好是商的最低有效整数位，相当于十进制数中的个位。所以在16次移位相减之后，累加器A低16位所得的结果为Q值为0的整数。以此分析，当商的精确值不是整数，或者超出Q值15所表示的范围时，此算法所得结果就达不到16位数据所能表达的精确度。

表1中任取几组数据来说明。


2 高精度算法

由于经典算法在结果上存在精度不高的问题，所以应考虑精度尽可能高的除法算法。

和经典算法一致，本算法也采用移位相减实现除法。但不同的是，在进行移位相减之前算法对除数与被除数进行了一些处理，并且在移位相减时不再分两种情况。下面具体介绍算法的流程。

首先给出被除数y和除数x。由于定点DSP中所有的数据都以整数的形式存在，所以在这里就把y和x看成是整数，不考虑它们本身的Q值。在图1所示的流程图中，A和B是C5416的2个累加器。算法第一步是将x和y相乘，结果存放在B中，以便在最后从B中提取商的符号(商的符号与B的符号相同)。第二步分别对x和y取绝对值以进行无符号数相除。第三步，分别判断x和y在最高有效位之前共有多少个O，分别记为n和m，即在x和y的分别有(16一n)和 (16一m)个有效位。然后将x和y分别左移(n一1)位和(m一2)位，即x=x·2n-1，y=y·2m-2。第四步，将前面所得x和y再进行移位相减，把y存放在A的高位，即AH=y，利用SUBC指令进行1 5次移位相减。在所得的结果当中，A的低位即为商的绝对值，B的符号即为商的符号，两者结合得到quot值。由于在移位相减之前分别对除数和被除数有一个左移放大的操作，所以在结果中就必须考虑商的Q值。



假设给定的除数和被除数是x0和y0，则进行移位相减的两数分别为x=x0·2n-1，y=y0·2m-2。根据计算的过程可以得到商的Q值为 (14+m-n)，所以最后的结果应该包括quot值和Q值(14+m—n)两部分。即商=quot·2-Q。

3 实验结果及讨论

前一节中详细介绍了高精度算法的具体流程，可以发现其与经典算法相比，不同之处在于进行移位相减之前把除数和被除数进行了不同位数的左移，保证进行移位相减的除数的数据位最高位是l，即14位为1，而被除数的第13位为1。正是这个左移使结果的精度得到了提高。

在所引用的经典算法中，当1被除数l<l除数I时所得的商为小数，而当|被除数|≥|除数|时所得的商为整数。实际上商为小数时其Q值为 15，商为整数时其Q值为O。根据后面给出的算法，结果的Q值为(14+m一n)。当|被除数|<|除数|时，可以得到n≤m，则商的Q值≥14；当|被除数|≥|除数|时，可以得到n≥m，则商的Q值≤14。从对比中可以发现，当|被除数|≥|除数|时结果的Q值提高了(14+m—n)。

利用高精度算法在TMS320C5416中计算上面的4组数据，所得结果如表2所列。



上面的实验说明，该算法使计算结果的精度有了很大的提高。取表中第三组数据进行详细分析，y=479、x=240，利用经典算法结果为整数1。利用改进后的算法，得到的结果有两部分：quot值=3FDDH，Q值=13。由这两部分可以算得精确的商值1．995 727 539 062 5，而479／240=1．995 83，计算结果与实际值之间的误差为O．000105 794 270 83。当Q值=13时，定点数据表示形式允许的误差为2-13=O．000 122 070 312 5，可见计算结果误差在允许的误差范围之内，也就是说计算结果满足所推导出的精度的范围。对比前后两种算法，可以发现这里的计算结果精度由原来的1提高到了2-13。

通过分析可知，该算法使计算结果在精度上有了非常大的提高，但是在效率上则有所降低。如表3所列：



表3中的时钟周期指的是DSP的机器时钟周期。由表3可以发现高精度算法对于不同的数据所需的周期数是不一样的，而且算法中计算被除数和除数的高位无效位数占用了大部分的周期数。

经典算法在效率上比高精度算法高，主要是因为在高精度算法中高位无效位数的计算效率过低。用C5416计算|被除数|和|除数|的高位无效位数，效率低是由于在C5416中没有现成的指令来探测数据中哪一位为最高有效位。表3给出的时钟周期数是在笔者的程序中所得到的结果，由于时间匆忙和优化方面能力的局限性，所以并不能保证在程序优化方面做得很好，所以在C5416中对该部分程序再进行更好的优化以减少执行的时间是很可能的。

在TI公司C6000系列的定点DSP中，有专门的指令LMBD来快速得到累加器中从符号位开始无效0位的个数，只需几个周期就能完成C5416里高位无效位数的计算，所以在C6000中完成整个算法的效率与经典算法相比肯定能够提高。并且在C6000系列中对不同的数据进行计算所需的时钟周期数是一样的，因此在C6O0O系列定点DSP中该算法有着更好的实用性。