非线性运算的定点快速实现

Aguilera

非线性运算的定点快速实现 [复制链接]

在数值运算中，除基本的加减乘除运算外，还有其它许多非线性运算，如，对数运算，开方运算，指数运算，三角函数运算等，实现这些非线性运算的方法一般有：(1)调用DSP编译系统的库函数；(2)查表法；(3)混合法。下面我们分别介绍这三种方法。
       1.调用DSP编译系统的库函数
       TMS320C2X/C5X的C编译器提供了比较丰富的运行支持库函数。在这些库函数中，包含了诸如对数、开方、三角函数、指数等常用的非线性函数。在C程序中(也可在汇编程序中)只要采用与库函数相同的变量定义，就可以直接调用。例如，在库函数中，定义了以10为底的常用对数log10()：
#include＜math.h＞
double，log10(double x)；
复制代码

       在C程序中按如下方式调用：
float x，y;
X=10.0;
y=log10(x)；
复制代码

       从上例可以看出，库函数中的常用对数log10()要求的输入值为浮点数，返回值也为浮点数，运算的精度完全可以保证。直接调用库函数非常方便，但由于运算量大，很难在实时DSP中得到应用。
       2.查表法
       在实时DSP应用中实现非线性运算，一般都采取适当降低运算精度来提高程序的运算速度。查表法是快速实现非线性运算最常用的方法。采用这种方法必须根据自变量的范围和精度要求制作一张表格。显然输人的范围越大，精度要求越高，则所需的表格就越大，即存储量也越大。查表法求值所需的计算就是根据输入值确定表的地址，根据地址就可得到相应的值，因而运算量较小。查表法比较适合于非线性函数是周期函数或已知非线性函数输入值范围这两种情况、例1.12和例1. 13分别说明这两种情况。
       例1.12  已知正弦函数y=cos(x)，制作一个512点表格，并说明查表方法。由于正弦函数是周期函数，函数值在-1至+1之间，用查表法比较合适。由于Q15的表示范围为1-至32767/32768之间，原则上讲-1至＋1的范围必须用Q14表示。但一般从方便和总体精度考虑，类似情况仍用Q15表示，此时+1用32767来表示。
       (1)产生5l2点值的C语言程序如下所示。
#define N 512
#define pi 3.14l59
int sin_tab[5l2]；
void  main()
{
int i;
for(i=0；i＜N；i++)sin_tab=(int)(32767*sin(2*pi*i/N))；
复制代码

       (2)查表
       查表实际上就是根据输人值确定表的地址。设输入x在0~2π之间，则x对应于512点表的地址为：index=(int)(512*x/2π)，则y=sin(x)=sin_tab[index]如果x用Q12定点数表示，将512/2π用Q8表示为20861，则计算正弦表的地址的公式为。
index=(x*20861L)>>20；
复制代码

       例1.12用查表法求以2为底的对数，已知自变量值范围为0.5-1，要求将自变量范围均匀划分为10等分。试制作这个表格并说明查表方法。
       (1)作表：
       y=log2(x)，由于x在0.5到1之间，因此y在-1到0之间，x和y均可用Q15表示。由于对x均匀划分为10段，因此，10段对应于输入x的范围如表3.2所示。若每一段的对数值都取第一点的对数值，则表中第一段的对数值为y0(Q15)=(int)(log(O.5)*32768)，第二段的对数值为y1(Q15)=(int)(log2(0.55)*32768)，依次类推，如表3.2所示。
       (2)查表：
       查表时，先根据输人值计算表的地址，计算方法为：
index=((x-16384)*20)>>15；
复制代码

       式中， index就是查表用的地址。例如，已知输人x=26869，则index=6，因此，y= -10549。
表1.2  logtab0 10点对数表
地址    输入值    对数值(Q15)
0    0.50-0.55    -32768
1    0.55-0.60    -28262
2    0.60-0.65    -24149
3    0.65-0.70    -20365
4    0.70-0.75    -16862
5    0.75-0.80    -13600
6    0.80-0.85    -10549
7    0.85-0.90    -7683
8    0.90-0.95    -4981
9    0.95-1.00    -2425
复制代码

       3.混合法
       (1)提高查表法的精度
       上述方法查表所得结果的精度随表的大小而变化，表越大，则精度越高，但存储量也越大。当系统的存储量有限而精度要求也较高时，查表法就不太适合。那么能否在适当增加运算量的情况下提高非线性运算的精度呢？下面介绍一种查表结合少量运算来计算非线性函数的混合法，这种方法适用于在输入变量的范围内函数呈单调变化的情形。混合法是在查表的基础上来用计算的方法以提高当输入值处于表格两点之间时的精度。提高精度的一个简便方法是采用折线近似法，如图1.1所示。
       图1.1提高精度的折线近似法”
       仍以求以2为底的对数为例(例1.12)。设输入值为x，则精确的对数值为y，在表格值的两点之间作一直线，用y'作为y的近似值，则有：
y'=y0＋△y
复制代码

       其中y0由查表求得。现在只需在查表求得y0的基础上增加△y既可。△y的计算方法如下：                   △y=(△x/△x0)△y=△x(△y0/△x0)
       其中△y0/△x0对每一段来说是一个恒定值，可作一个表格直接查得。此外计算此时需用到每段横坐标的起始值，这个值也可作一个表格。这佯共有三个大小均为10的表格，分别为存储每段起点对数值的表logtab0、存储每段△y0/△x0值的表logtab1和存储每段输入起始值x0的表logtab2，表logtab1和表logtab2可用下列两个数组表示。
int  logtab1[10]={22529，20567，18920，17517，16308，
15255，14330，13511，12780，12124}；/*△y0/△x0：Q13*/
int  logtab2[10]={16384，18022，19660，21299，22938，
24576，26214，27853，29491，31130}；/*x0：Q15*/
复制代码

       综上所述，采用混合法计算对数值的方法可归纳为：
       (1)根据输人值，计算查表地址：index=((x-16384)*20)>>15；
       (2)查表得y0=logtab0[index]；
       (3)计算△x=x-logtab2[index]；
       (4)计算△y=(△x*logtab1[index])>>13；
       (5)计算得结果y=y0＋△y。
       例1.13已知x=0.54，求log2(x)。
       0.54的精确对数值为y=log2(0.54)=-0.889。
       混合法求对数值的过程为：
       (1)定标Q15，定标值x=0.54*32768=17694；
       (2)表地址index=((x-16384)*20)>>15=0；
       (3)查表得y0=logtab0[0]=-32768；
       (4)计算△x=x-logtab2[0]=17694-16384=1310；
       (5)计算△y=(△xlogtab1[0]>>13=(13l0*22529L)>>13=3602
       (6)计算结果y=y0＋△y=-32768＋3602=-29166。
       结果y为Q15定标，析算成浮点数为-29166/32768=-0.89，可见精度较高。
       (2)扩大自变量范围
       如上所述，查表法比较适用于周期函数或自变量的动态范围不是太大的情形。对于像对数这样的非线性函数，输入值和函数值的变化范围都很大。如果输入值的变化范围很大，则作表就比较困难。那么能否比较好地解决这个问题，即不便表格太大，又能得到比较高的精度呢？下面我们来讨论一种切实可行的方法。
       设x是一个大于0.5的数，则x可以表示为下列形式：
x=m*2e
式中，0.5<=m<=1.0，e为整数。则求x的对数可以表示为：
    log2(x)=log2(m*2e)=log2(m)＋log2(2e)=e＋log2(m)
复制代码

       也就是说，求x的对数实际上只要求m的对数就可以了，而由于m的数值在0.5和1.0之间，用上面介绍的方法是完全可以实现的。例如：
    log2(10000)=log2(0.61035*214)=log2(0.61035)+14 =13.2877
复制代码

       可见，如果一个数可以用比较简便的方法表示为上面的形式，则求任意大小数的对数也比较方便的。TMS320C2X/C5X指令集提供了一条用于对ACC中的数进行规格化的指令NORM，该指令的作用就是使累加器中的数左移，直至数的最高位被移至累加器的第30位。例如，对数值10000进行规格化的TMS320C25程序为。
LAC #10000
SACL TEMP
ZALH TEMP
LAR AR1，#0FH
RPT 14
NORM * -
复制代码

       上述程序执行后，AR1=#0eH，ACCH=2000(10进制)。对一个16位整数x进行上述程序处理实际上就是作这样一个等效变换：
    x=[(x*2e)/32768]*215-Q
复制代码

       其中，寄存器AR1包含的值为15-Q累加器ACC高16位包含的值为x.2Q，其数值在16384至32768之间。
       例1.14实现以2为底的对数的C定点模拟程序。

int logtab0[10]={-32768，-28262，-24149，-20365，-16862，
-13600)，-1O549，-7683，-4981，-2425};/*Q15*/
int logtab1[10]={122529，20567，18920，175l7，16308，
15255，14330，13511，12780，12124};/*Q13*/
int logtab2[10]={16384，l8022，19660，21299，22938，
24576，26214，27853，29491，31130};/*Q15*/
int log2_fast(int Am)
{
int point，point1；
int index，x0，dx，dy，y;
point=0；
while(Am＜16384){point＋＋；Am=Am＜＜1；}/*对Am进行规格化*/
point1=(15-point-4)*512；/*输入为Q4，输出为Q9*/
index=((Am-16384)*20L)＞＞15；/*求查表地址*/
dx=Am-logtab2[index]；
dy=((long)dx*logtab1[index])＞＞13；
y=(dy＋longtab0[index])＞＞6；/*Q9*/
y=point1＋y;
return(y)；
}
复制代码

       上述程序中，输入值Am采用Q4表示，输出采用Q9表示，如果输入输出的Q值与上面程序中的不同，则应作相应的修改。
       以上讨论了DSP芯片进行定点运算所涉及的一些基本问题，这些问题包括：数的定标，DSP程序的定点模拟，DSP芯片的足点运算以及定点实现非线性函数的快速实现方法等。充分理解这些问题对于用定点芯片实现DSP算法具有非常重要的作用。

lzz8031

快速实现，可以在程序中实现吗？