与afreez一起学习DSP中浮点转定点运算

http://blog.csdn.net/ganxingming/article/details/1449526
------------------------------------------------------------------------------------------------
这段时间在学习DSP浮点运算转定点运算过程中,在网上看到了一位网友给出的一段程序.如下所示: 　　---------------------------------------------------------
　　声明：
　　 此文为原创，欢迎转载，转载请保留如下信息
　　 作者：聂飞（afreez） 北京-中关村
　　 联系方式：afreez@sina.com （欢迎与作者交流）
　　 初次发布时间：2006-11-28
　　 不经本人同意，不得用语商业或赢利性质目的，否则，作者有权追究相关责任！
　　---------------------------------------------------------
　　/***************************************************************
　　// 将32位浮点数fval转换为32位整数并存储在ival中
　　// 小数部分将被裁剪掉
　　//
　　void TruncToInt32 (int &ival, float fval)
　　{
　　ival = *(int *)&fval;
　　// 提取尾数
　　// 注意实际的尾数前面还有一个被省略掉的1
　　int mantissa = (ival & 0x07fffff) | 0x800000;
　　// 提取指数
　　// 以23分界，指数大于23则左移，否则右移
　　// 由于指数用偏移表示，所以23+127=150
　　int exponent = 150 - ((ival >> 23) & 0xff);
　　if (exponent < 0)
　　ival = (mantissa << -exponent);
　　else
　　ival = (mantissa >> exponent);
　　// 如果小于0，则将结果取反
　　if ((*(int *)&fval) & 0x80000000)
　　ival = -ival;
　　}
　　*******************************************************/
　　程序流程写的很清晰.兴奋之余马上在自己的机器上测试了一把,结果发现总是不正确,大体如下:
　　//主函数
　　...
　　int ival;
　　float fval=4.23;
　　TruncToInt32 (&ival, fval);
　　//这里的输出结果不正确
　　printf("ival=%d ",ival);
　　...
　　//函数定义
　　void TruncToInt32 (int &ival, float fval)
　　{
　　//skip
　　...
　　// 如果小于0，则将结果取反
　　if ((*(int *)&fval) & 0x80000000)
　　ival = -ival;
　　//发现在这里输出的结果是正确的
　　printf("ival=%d ",ival);
　　}
　　仔细看了一眼，才发现，原来是函数调用时参数传递有问题。主要上理解指针和指针的指针，大家对照我的思路考虑考虑吧。
　　修改主程序，如下所示：
　　...
　　int ival;
　　float fval=4.23;
　　RightTruncToInt32 (&ival, fval);//第n行,假定是经过修正后的函数
　　ival=(int)fval;//第n+1行
　　printf("ival=%d ",ival);
　　...
　　对n行和n+1行进行时间效率测试，显示，第n+1行所需的计算时间少于第n行。
　　对主程序进一步改善，大致如下：
　　...
　　int ival;
　　float fval=4.23;
　　//-------Begin more fast section?-----------
　　// 提取尾数
　　int ival= ((*(int *)(&fval)) & 0x07fffff) | 0x800000;
　　// 提取指数
　　int exponent = 150 - (((*(int *)(&fval)) >> 23) & 0xff);
　　if (exponent < 0)
　　ival = (ival<< -exponent);
　　else
　　ival = (ival >> exponent);
　　// 如果小于0，则将结果取反
　　if ((*(int *)&fval) & 0x80000000)
　　ival = -ival;
　　//-------End more fast section?-----------
　　ival=(int)fval;//第n+1行
　　...
　　结果显示，在RedHat 9.0平台上测试显示，如果编译不加优化选项，more fast section的执行时间效率优与第n+1行；加优化选项，第n+1行的执行时间效率优与more fast section。
　　分析结论：
　　由于more fast section还有一些代码行可以优化，我想，在一些特定场合下，还是可以借鉴这个思路的 ----------------------______________________________________ -----------————————————————————————
本文介绍在计算机体系结构中浮点floating-point和定点fixed-point数据的表示，浮点到定点转换的各种舍入和截断处理方式，定点数据的Q格式表示，根据信号处理算法来进行定点数据的定标方法，定点数据的加减乘除各种数学模拟操作方法，最后以FFT运算为例说明如何在定点DSP上提高定点运算的精度 计算机体系结构中浮点和定点数据的表示 1、定点数：
定点数指小数点在数中的位置是固定不变的，通常有定点整数和定点小数或者说是定点分数。在对小数点位置作出选择之后即定标定了Q值后，运算中的所有数均应统一为定点整数或定点小数，在运算中不再考虑小数点的位置问题。
（1）定义：数据中小数点位置固定不变的数
（2）种类：定点整数、定点小数(分数fraction)
（3）小数点在符号位与有效位之间，
需要注意的是 定点数受字长的限制，超出范围会有溢出。
2、浮点数：
浮点数中小数点的位置是不固定的，用指数和尾数来表示。通常尾数为纯小数，指数为整数，尾数和指数均为带符号数。尾数的符号表示数的正负；指数的符号则表明小数点的实际位置。
（1）形式：V＝MX2E
（2）M：尾数Mantissa
（3）E：指数Exponent
（4）在计算机中M和E表示形式为尾数符号、指数、尾数 需要注意的是，浮点数的精度由尾数决定，数的表示范围由指数决定。
IEEE754的浮点表示 a) 两种基本浮点格式：单精度和双精度。 IEEE单精度格式具有24位有效数字（1个符号位，23位尾数），8位的指数部分，总共占用32 位。IEEE双精度格式具有53位有效数字精度（1个符号位，52位尾数），11位的指数部分，并总共占用64位。 说明：基本浮点格式是固定格式，相对应的十进制有效数字分别为7位和17位。基本浮点格式对应的C/C++类型为float和double。 b) 浮点到定点的四种舍入(截断)方向： 向最接近的可表示的值；当有两个最接近的可表示的值时首选"偶数"值；向负无穷大（向下）；向正无穷大（向上）以及向0（截断）。 说明：舍入模式也是比较容易引起误解的地方之一。IEEE 754标准根本不支持四舍五入模式，它的默认模式是最近舍入（Round to Nearest），它与四舍五入只有一点不同，对.5的舍入上，采用取偶数的方式。如： 最近舍入模式：Round(0.5) = 0; Round(1.5) = 2; Round(2.5) = 2; 四舍五入模式：Round(0.5) = 1; Round(1.5) = 2; Round(2.5) = 3; 向0（截断）舍入：C/C++的类型转换。(int) 1.324 = 1，(int) -1.324 = -1; 向负无穷大（向下）舍入：C/C++函数floor()。例如：floor(1.324) = 1，floor(-1.324) = -2。 向正无穷大（向上）舍入：C/C++函数ceil()。ceil(1.324) = 2。Ceil(-1.324) = -1; 后两种舍入方法据说是为了数值计算中的区间算法，但很少被应用。 3、定点数与浮点数区别
定点表示法运算直观，但数的表示动态范围较小(对于32位定点整数而言，数据范围只能从-2147483648到2147483647)，不同的数运算时要考虑比例因子的选取，以防止溢出，对于某些应用场合，如高品质音频，要求高达24bit即144dB的动态范围，使用定点运算就需要不断的调整比例因子以满足精度要求。浮点表示法运算时可以不考虑溢出（即使对于32位单精度浮点而言，其数据范围也能从，但浮点运算在实际的平台实现复杂度较高，即运算复杂度较高，要掌握定、浮点数的转换方法及浮点数规格化方法。  单精度格式浮点数表示范围
             指数     尾数       数值(二进制)           数值(十进制)
最大数       0xfe     0x7fffff   1.(23个1)  * 2^(+127)  3,4028236692093846346337460743177e+38
最小数       0x1      0x0        1.(23个0)  * 2^(-126)  1,1754943508222875079687365372222e-38
最小弱规范数 0x0      0x1        0.(22个0)1 * 2^(-126)  1,4012984643248170709237295832899e-45 定点数据的定标Q 一般的定点处理器的寄存器字长都是有限的，一般为16位或者32位，还有24位的DSP处理器，当然字长越长，所能表示的数据范围越大，精度也越高。定点DSP芯片的数以2的补码形式表示，以16位字长为例。每个16位数用一个符号位来表示数的正负，0表示数值为正，l则表示数值为负。其余15位表示数值的大小。因此，
二进制数0010000000000000b=0x2000H=8192
二进制数1111111111111100b=0xfffcH=-4
对于定点处理器而言，它并不清楚数据的小数点位置，它只能根据指令把数据做统一的运算处理，此时程序员需要确定小数点的位置即定标，然后定标过程中需要考虑数据的动态范围和精度要求。定点数的定标一般采用Q表示法，Q表示小数点在字宽中的位置。同样一个16位数，若小数点设定的位置不同（即定标不同），它所表示的数也就不同。例如，
16进制数2000H=8192，用Q0表示
16进制数2000H=0.25，用Q15 显然，不同的Q所表示的数不仅范围不同，而且精度也不相同。Q越大，数值范围越小，但精度越高；相反，Q越小，数值范围越大，但精度就越低。例如，Q0 的数值范围是一32768到+32767，其精度为1，而Q15的数值范围为-1到0.9999695，精度为1/32768=0.00003051。因此，对定点数而言，数值范围与精度是一对矛盾，一个变量要想能够表示比较大的数值范围，必须以牺牲精度为代价；而想精度提高，则数的表示范围就相应地减 小。在实际的定点算法中，为了达到最佳的性能，必须充分考虑到这一点，即权衡动态范围和精度。
浮点数与定点数的转换关系可表示为：
浮点数(x)转换为定点数(xq)：xq=(int)x* 2Q，当然浮点转换为顶点时需要考虑以上提到的舍入方式，不同的舍入方式，其精度也不同。
定点数(xq)转换为浮点数(x)：x=(float)xq*2-Q 定点模拟的浮点操作 浮点的加减法需要先把小数点统一到一个相同的位置，然后再加减，统一小数点位置可以用左移到一个Q值高的定点位置或者右移到Q值低的定点位置，显然左移的实现需要保证加减的结果不会导致溢出，而右移的方法需要考虑数据精度是否满足需求。一般的DSP处理器都不会考虑加减运算导致的数值溢出问题，因而程序员需要用更长的字长来防止溢出，或者采用饱和的方式以防止数据的上溢或者下溢，有些DSP芯片会有专门的饱和加减的操作，有些DSP还会包含更长字长的累加器来存储加减的结果以防止溢出。 定点乘法运算的模拟，则需要知道乘数和被乘数的Q值以及乘积结果的Q值，
设浮点乘法运算的表达式为：
float x，y，z；
z=xy；
假设经过统计后x的定标值为Qx，y的定标值为Qy，乘积z的定标值为Qz，则
z=xy
zq*2-Qx=xq*yq*2-(Qx+Qy)
zq=(xqyq)2Qz-(Qx+Qy)
所以定点表示的乘法为：
int x，y，z；
long temp；
temp=(long)x；
z=(temp*y)>>(Qx+Qy-Qz)；
例 定点乘法。
设x=18.4，y=36.8，则浮点运算值为=18.4*36.8=677.12；
根据上节，得Qx=10，Qy=9，Qz=5，所以
x=18841；y=18841；
temp=18841L；
z=(18841L*18841)>>(10+9-5)=354983281L>>14=21666；
因为z的定标值为5，故定点z=21666，即为浮点的z=21666/32=677.08。 对于定点除法运算而言，和乘法类似，需要清楚数据的Q值，然后对结果做调整。
设浮点除法运算的表达式为：
float x，y，z；
z=x/y；
假设经过统计后被除数x的定标值为Qx，除数y的定标值为Qy，商z的定标值为Qz，则
z=x/y
zq*2-Qz=(xq*2-Qx)/(yq*2-Qy)
zq=(xq*2(Qz-Qx+Qy))/yq
所以定点表示的除法为：
int x，y，z；
long temp；
temp=(long)x；
z=(temp<<(Qz-Qx+Qy))/y；
例1.6定点除法。
设x=18.4，y=36.8，浮点运算值为z=x/y=18.4/36.8=0.5；
根据上节，得Qx=10，Qy=9，Qz=15；所以有
z=18841，y=18841；
temp=(long)18841；
z=(18841L<<(15-10+9)/18841=3O8690944L/18841=16384；
因为商z的定标值为15，所以定点z=16384，即为浮点z=16384/215=0.5。 因而为了适应各种定点运算操作，需要清楚操作数和结果的Q值，实际上Q值对应于变量的动态范围，只要清楚数据的动态范围，也就能确定Q值了。一般可以通过对数据的动态范围进行理论和统计的分析来得出合理的动态范围，而且在进行数据操作时，还需要考虑饱和运算以防止溢出。 如何优化定点的运算精度     定点计算需要注意的就是保证精度，数据动态范围很大的情况下，有限的字长表示不了那么大动态范围的情况下，有以下方法保证尽可能高的运算精度

1. 分级缩放数据，以FFT计算为例，可以分级scale，蝶形运算的理论表明每一级的增益不会大于2，所以每次蝶形运算后做一次右移要比先把输入数据缩放到不会溢出的水平上精度更高。
2. Normalization的归一化，同样的FFT运算，首先根据输入数据的幅度的最大情况，计算一个不会导致溢出的缩放因子，然后对整个输入做归一化的scale，FFT之后再de-normalization能显著的提高精度，同时又能适应输入数据较大的动态范围。
3. Pseudo floating-point：即把定点数据表示为指数和尾数的形式，用虚拟浮点数据的方式来满足高动态范围的需求，当然pseudo floating-point的表示还是需要做归一化来求解指数的。

Reference： http://www.blog.163.com/houh-1984/ 本文介绍在计算机体系结构中浮点floating-point和定点fixed-point数据的表示，浮点到定点转换的各种舍入和截断处理方式，定点数据的Q格式表示，根据信号处理算法来进行定点数据的定标方法，定点数据的加减乘除各种数学模拟操作方法，最后以FFT运算为例说明如何在定点DSP上提高定点运算的精度
===============================================

浮点数转换成定点整数 在计算机图形运算中，常常要将浮点数转换为整数，例如在图像的光栅化阶段，就要执行大量的类型转换，以便将浮点数表示的坐标转化为整数表示的屏幕坐标。Ok，it's so easy：
----------------------------------------------------------------------------------------
//
// 强制类型转换
// 小数部分将被裁剪掉
//
int_val = (int)float_val;
----------------------------------------------------------------------------------------
嘿嘿，很高兴你居然和我一样单纯！这个操作实在是太TINY了，以至于我从来没想过它是怎么实现的，直到某天某个家伙跟我说，不要使用标准C类型转换，因为那太慢了！我当时的震惊不下于倒霉的冒险者遇上了龙。

标准C类型转换最大的优点是，它是独立于平台的，无论是在X86上跑，还是在PowerPC上跑，你什么都不用担心，编译器会为你搞定一切。而这也恰恰是它最大的缺点——严重依赖于编译器的实现。而实际测试表明，编译器所生成的代码，其速度实在不尽人意。

一个替代的方法是直接对数据位进行操作。如果你对IEEE浮点数的表示法比较熟悉的话（如果你压根什么都不知道，请先查阅文末附录中的资料），这是显而易见的。它提取指数和尾数，然后对尾数执行移位操作。代码如下：
----------------------------------------------------------------------------------------
//
// 将32位浮点数fval转换为32位整数并存储在ival中
// 小数部分将被裁剪掉
//
void TruncToInt32 (int &ival, float fval)
{
ival = *(int *)&fval;

// 提取尾数
// 注意实际的尾数前面还有一个被省略掉的1
int mantissa = (ival & 0x07fffff) | 0x800000;

// 提取指数
// 以23分界，指数大于23则左移，否则右移
// 由于指数用偏移表示，所以23+127=150
int exponent = 150 - ((ival >> 23) & 0xff);

if (exponent < 0)
ival = (mantissa << -exponent);
else
ival = (mantissa >> exponent);

// 如果小于0，则将结果取反
if ((*(int *)&fval) & 0x80000000)
ival = -ival;
}
----------------------------------------------------------------------------------------
该函数有一个BUG，那就是当fval=0时，返回值是2。原因是对于语句mantissa>>exponent，编译器使用了循环移位指令。解决方法是要么对0作特殊处理，要么直接用汇编来实现。

这个函数比标准的C转换要快，而且由于整个过程只使用了整数运算，可以和FPU并行运行。但缺点是，(1)依赖于硬件平台。例如根据小尾和大尾顺序的不同，要相应地修改函数。(2)对于float和double要使用不同的实现，因为二者的数据位不同。(3)对于float，只能保留24位有效值，尽管int有32位。

更快的方法是使用FPU指令FISTP，它将栈中的浮点数弹出并保存为整数：
----------------------------------------------------------------------------------------
//
// 将64位浮点数fval转换为32位整数并存储在ival中
// 小数部分将四舍五入到偶数
//
inline void RoundToIntFPU (int &ival, double fval)
{
_asm
{
fld fval
mov edx, dword ptr [ival]
fistp dword ptr [edx]
}
}
----------------------------------------------------------------------------------------
很好，无论速度还是精度似乎都相当令人满意。但如果换一个角度来看的话，fistp指令需要6个cycle，而浮点数乘法才仅仅需要3个cycle！更糟的是，当fistp运行的时候，它必须占用FPU，也就是说，其他的浮点运算将不能执行。仅仅为了一次类型转换操作就要付出如此大的代价，光想想就觉得心疼。

当然，它也有很多优点：更快的速度，更精确的数值（四舍五入到偶数），更强的适用性（float和double均可)。要注意的是，FPU默认的四舍五入到偶数（round to even）和我们平常所说的四舍五入（round）是不同的。二者的区别在于对中间值的处理上。考虑十进制的3.5和4.5。四舍五入到偶数是使其趋向于邻近的偶数，所以舍入的结果是3.5->4，4.5->4；而传统的四舍五入则是3.5->4，4.5->5。四舍五入到偶数可以产生更精确，更稳定的数值。

除此之外，还有更好，更快的方法吗？有的，那就是华丽的 Magic Number ！

请看下面的函数：
----------------------------------------------------------------------------------------
//
// 将64位浮点数转换为32位整数
// 小数部分将四舍五入到偶数
//
inline void RoundToInt64 (int &val, double dval)
{
static double magic = 6755399441055744.0;
dval += magic;
val = *(int*)&dval;
}
----------------------------------------------------------------------------------------
快，这就是它最伟大的地方！你所需要的仅仅是一次浮点数加法，你还能再奢望些什么呢？

当然，绝对不要使用莫名其妙的代码，现在马上就让我们来看看它是怎么变的戏法。

首先，我们要搞清楚FPU是怎样执行浮点数加法的。考虑一下8.75加2^23。8.75的二进制表示是1000.11，转化为IEEE标准浮点数格式是1.00011*2^3。假设二者均是32位的float，它们在内存中的存储为：
----------------------------------------------------------------------------------------
符号位(31), 指数(30-23), 尾数(22-0)

8.75: 0, 10000010, 0001 1000 0000 0000 0000 000

2^23: 0, 10010110，0000 0000 0000 0000 0000 000
----------------------------------------------------------------------------------------
FPU所执行的操作是：(1)提升指数较小者的指数，使二者指数相同；(2)将二者的尾数相加；(3)将结果规整为标准格式。让我们具体来看看整个过程：
----------------------------------------------------------------------------------------
1，提升8.75的指数，尾数要相应地右移23-3=20位：

8.75 = 1.00011*2^3 = .0000000000000000000100011*2^23

2，将二者相加。必须注意的是，
   由于float只有23位尾数，所以8.75的低位被移除掉了：

8.75: 0, 10010110, 0000 0000 0000 0000 0001 000
+
2^23: 0, 10010110，0000 0000 0000 0000 0000 000
=
0, 10010110, 0000 0000 0000 0000 0001 000

3，将规整为标准格式：

别忘了2^23还有个前导1，所以结果是规整的，无需执行任何操作
----------------------------------------------------------------------------------------
你发现什么了吗？不错，将结果的尾数部分提取出来，正是int(8.75)！是的，magic number的奥妙就在这里，通过迫使FPU将尾数移位，从而获得我们需要的结果。

但是别高兴得太早，让我们来看看负数的情况。以-8.75为例，-8.75加2^23相当于2^23减去8.75，过程如下：
----------------------------------------------------------------------------------------
2^23: 0, 10010110，0000 0000 0000 0000 0000 000
-
8.75: 0, 10010110, 0000 0000 0000 0000 0001 000
=
0, 10010110, 1111 1111 1111 1111 1110 000
----------------------------------------------------------------------------------------
很好，尾数部分正是int(-8.75)=-8的补码。然而，2^23的前导1在减法过程中被借位了，所以结果的尾数前面并没有1，FPU将执行规整操作，最后我们得到的是：
----------------------------------------------------------------------------------------
0, 10010110, 1111 1111 1111 1111 1100 000
----------------------------------------------------------------------------------------
功亏一篑！等等，如果将2^23的尾数的最高位置1，那么在减法过程中不就可以保全前导1了吗？完全正确，这就是我们需要的。所以最后的magic number是0,10010110,1000 0000 0000 0000 0000 000，也就是1.5*2^23。

然而，我们要清楚这个方法的一些局限性。首先，在将结果的float值保存为整数时，我们需要使用某些mask值将22-31位屏蔽掉。其次，float只能保有最多23位的有效值，在将尾数最高位置1后，有效值更是降到了22位，这意味着我们对大于2^23-1的数值无能为力。

相比之下，如果我们只处理double，那么所有的问题就都迎刃而解了。首先，double的指数位，符号位和尾数的最高位全部都在高32位，无需屏蔽操作。其次，double有多达52位的尾数，对于32位的int来说，实在是太奢侈了！

用于double的magic number是1.5*2^52=6755399441055744.0，推导过程是一样的。

根据同样的原理，我们还可以计算出将float和double转化为定点数的magic number。例如对于16-16的定点数，尾数右移的位数比原先转化为整数时要少16，所以对于double来说，相应的magic number就是1.5*2^36。如果要转化为8-24的定点数，则移位次数要少24，magic number是1.5*2^28。对于其他格式的定点数，以此类推。比起int(float_val*65536)这样的表达式，无论是速度还是精度都要优胜得多。

另外，不得不再次重申的是，对于在最后的结果中被移除掉的数值，FPU会将其四舍五入到偶数。然而，有时我们确实需要像floor和ceil那样的操作，那该怎么办呢？很简单，将原数加上或减去一个修正值就可以了，如下所示：
----------------------------------------------------------------------------------------
// 修正值
static double magic_delta=0.499999999999;

// 截取到整数
inline void Floor64 (int &val, double dval)
{
RoundToInt64(val,dval-delta);
}

// 进位到整数
inline void Ceil64 (int &val, double dval)
{
RoundToInt64(val,dval+delta);
}
----------------------------------------------------------------------------------------
这个世界上没有免费的午餐。我们获得了速度，相对应地，也必须付出一些代价，那就是可移植性。上述方法全都基于以下几个假设：(1)在x86上跑；(2)符合IEEE的浮点数标准；(3)int为32位，float为32位，double为64位。局限性其实是蛮大的，相比之下，int_val=(int)float_val就要省事多了。当然，你也可以使用条件编译。

最后，让我们来看两组实际的测试数值。这份报告来自于Sree Kotay和他的朋友Herf：
----------------------------------------------------------------------------------------
平台：Pentium IV/1.2

64位浮点数到32位整数的转换：

int(f)：        2772.65 ms
fistp：         679.269 ms
magic number：  622.519 ms

64位浮点数到16-16定点数的转换：

int(f*65536)：  2974.57 ms
fistp：         3100.84 ms
magic number：  606.80 ms

floor函数：

ANSI floor：    7719.400 ms
magic number：  687.177 ms
----------------------------------------------------------------------------------------
平台：Pentium D/3.2

64位浮点数到32位整数的转换：

int(f)：        1948.470 ms
fistp：         523.792 ms
magic number：  333.033 ms

64位浮点数到16-16定点数的转换：

int(f*65536)：  2163.56 ms
fistp：         7103.66 ms
magic number：  335.03 ms

floor函数：

ANSI floor：    3771.55 ms
magic number：  380.25 ms
----------------------------------------------------------------------------------------
性能的差距实在令人惊讶，不是吗？我想说的是，写编译器的家伙绝对不是傻瓜（恐怕比你我都要聪明得多），他们当然知道所有这些优秀的算法。但为什么编译器的表现会如此糟糕呢？其中一个理由是，为了使浮点数运算尽可能精确和迅速，IEEE在算法的选择上有很多的限制。而另一方面，IEEE的舍入规则（四舍五入到偶数）尽管从维持浮点数的连贯性上来看非常棒，但却不符合ANSI C在浮点数到整数的类型转换上的说明（截尾）。于是，编译器不得不做一大堆的工作来保证行为的正确性（符合标准）。这在大部分情况下都不是什么问题——除非你在写图形/声音/多媒体之类的代码。

刚刚在我的赛扬1.1G上实际测试了一下。编译器是VC2003，使用PerformenceCounter来计算时间，在DEBUG模式下，C标准转换/FISTP/MagicNumber三种方法所耗费时间的比约为5/3/2。但在优化全开的RELEASE模式下，标准C类型转换的速度却是遥遥领先于所有其他的方法！也不知道是我的测试方法有问题，还是该赞VS厉害。
--------------------------------------------------------------------------------------
参考文献和相关资源可到鄙人的小店下载：http://greatsorcerer.go2.icpcn.com/info/float2int.html

1，What Every Computer Scientist Should Know About Floating-Point Arithmetic by David Goldberg（PDF）：
这篇论文囊括了关于浮点数的所有东西，正如其名字所示，任何想要了解浮点数的人都必读的文献。（其中关于精度的讨论实在令我受益非浅。）

2，Let's Get to The Floating Point by Chris Hecker（PDF）：
早期的关于浮点数的文章，写得非常棒，值得一看。

3，IEEE Standard 754 for Binary Floating Point Arithmetic by Prof.W.Kahan（PDF）：
关于IEEE754标准的说明。

4，IA64 Floating Point Operations and the IEEE Standard for Binary Floating Point Arithmetic（PDF）：
关于IA64的浮点数实现。

5，Know Your FPU: Fixing Floating Fast by Sree Kotay
=================================================== 有 一 二三
请看网址 ： http://blog.sina.com.cn/s/blog_5f853eb10100sjh3.html

算法移植－浮点转定点运算（一）

什么是定点数、浮点数？

首先我们要认清一个概念，定点数不一定是整数，浮点数不一定是小数。

    如其名，浮点数和定点数的区别就在于浮点和定点上，点就是指小数点。浮点数就是小数点是浮动的，定点数就是小数点是固定不动的。


   具体，什么是浮点数？

   浮点数是在计算机中用以近似表示任意某个实数。具体的说，这个实数由一个整数或定点数（即尾数）乘以某个基数（计算机中通常是2）的整数次幂得到，这种表示方法类似于基数为10的科学记数法。

   一个浮点数a由两个数m和e来表示：a = m ×b^e。在任意一个这样的系统中，我们选择一个基数b（记数系统的基）和精度p（即使用多少位来存储）。m（即尾数）是形如±d.ddd...ddd的p位数（每一位是一个介于0到b-1之间的整数，包括0和b-1）。如果m的第一位是非0整数，m称作规格化的。有一些描述使用一个单独的符号位（s代表+或者-）来表示正负，这样m必须是正的。e是指数。

  对于一些外文资料中，指数即：exponent，尾数即：mantissa。

 

 

     在IEEE754中，定义了两种浮点数，即我们熟悉的float和double型。如图所示。对于float型的浮点数来说，最高一位是符号位，不用说了，1为负号，0为正。紧跟着指数位是8位，尾数是23位。由于尾数是规格化的，最高一位肯定非零，并且最高一位隐藏。所以对于尾数来说，实际上可以有23+1=24位。

比如说，如果mantissa为010，exponent为3，s为0，则尾数实际上是1.010，因此这个数是+1.010*2^3=1010b,即为十进制10.0。

   Double的位数说明类似于float。


什么是定点数？

   我们在上述的浮点数中可以看到，浮点的小数位是可变的（随exponent变化），因此浮点数可表达的小数范围非常广。但浮点数运算量非常大（从它的定义上就知道了）。并且在目前市场占有量最大的定点DSP并不支持浮点运算。

   因此，定点数应运而生。定点数就是指在一个数中，整数部分和小数部分位数固定。比如，我们定点数总共32位，其中小数占低13位：

typedef int fix_t;

#define FIX_FRACBITS 13

如何将一个int型的整数转化为定点数？

  Int类型转定点数最简单。上面的定义，小数位占低13位。

#define INTTOFIX(fix_t, fracbits, x) fix_t（(x)<< (fracbits)）

   即我们只需要将x左移就行了，当然要注意，如果x本身超过了整数部分的最大位数，则会产生溢出。


如何将一个double型小数转化为定点数？

   浮点数转为定点数与int型转为定点数类似，比如说1.1101，转为有13位小数的定点数就是1 1 1010 00000000，也就是1.1101*2^13。

   当然，在int型转定点时，我们使用了逻辑运算，而不是直接乘以2^13，是因为逻辑运算的速度一般比算术运算快几十倍。但是逻辑左移并不能运用于小数，double转定点时于是只能用算术乘法了。

   下面是double转定点：

#define DBLTOFIX(fix_t, fracbits, x)

    fix_t（((x) *（double)(（fix_t）(1) << (fracbits))))）

 

定点数直接如何相加？

  很简单，如下。如下，下面的宏定义的定点数相加没有考虑到数据溢出的问题，如果当x和y同号并且两数足够大时，是会产生数据位溢出的。在具体编程实现时要特别注意。

#define FIX_ADD(fix_t, fracbits, x, y) ((x) + (y))


定点数如何相乘？

   乘法更容易产生数据位溢出的问题，因此，如果我们运算时就要特别小心。同样，也是没有考虑到数据溢出的问题。


Typedef _int64 fixbig_t;

#define FIX_MUL(fix_t, fracbits, bigfix_t, x, y)

   Fix_t((bigfix_t(x) * bigfix_t(y)) >>(fracbits))


定点数如何实现除法？

  如下式所示，我们把移位放到被除数那里，是因为这样的处理基本上不会产生下溢出，提高了精度。

#define   JAS_FIX_DIV_FAST(fix_t, fracbits, bigfix_t, x, y)

   JAS_CAST(fix_t, (JAS_CAST(bigfix_t, x)<< (fracbits)) / (y))
=========================================================
一  ＤＳＰ定点算数运算
1  数的定标
    在定点DSP芯片中,采用定点数进行数值运算,其操作数一般采用整型数来表示。一个整型数的最大表示范围取决于DSP芯片所给定的字长,一般为16位或24位。显然,字长越长,所能表示的数的范围越大,精度也越高。如无特别说明,本书均以16位字长为例。
DSP芯片的数以2的补码形式表示。每个16位数用一个符号位来表示数的正负,0表示数值为正,l则表示数值为负。其余15位表示数值的大小。因此,
      二进制数0010000000000011b=8195
      二进制数1111111111111100b= -4
    对DSP芯片而言,参与数值运算的数就是16位的整型数。但在许多情况下,数学运算过程中的数不一定都是整数。那么,DSP芯片是如何处理小数的呢？应该说,DSP芯片本身无能为力。那么是不是说DSP芯片就不能处理各种小数呢？当然不是。这其中的关键就是由程序员来确定一个数的小数点处于16位中的哪一位。这就是数的定标。
通过设定小数点在16位数中的不同位置,就可以表示不同大小和不同精度的小数了。数的定标有Q表示法和S表示法两种。表1.1列出了一个16位数的16种Q表示、S表示及它们所能表示的十进制数值范围。
    从表1.1可以看出,同样一个16位数,若小数点设定的位置不同,它所表示的数也就不同。例如,
         16进制数2000H=8192,用Q0表示
         16进制数2000H=0.25,用Q15表示
但对于DSP芯片来说,处理方法是完全相同的。
    从表1.1还可以看出,不同的Q所表示的数不仅范围不同,而且精度也不相同。Q越大,数值范围越小,但精度越高;相反,Q越小,数值范围越大,但精度就越低。例如,Q0 的数值范围是一32768到+32767,其精度为1,而Q15的数值范围为-1到0.9999695,精度为1/32768=0.00003051。因此,对定点数而言,数值范围与精度是一对矛盾,一个变量要想能够表示比较大的数值范围,必须以牺牲精度为代价;而想精度提高,则数的表示范围就相应地减小。在实际的定点算法中,为了达到最佳的性能,必须充分考虑到这一点。
浮点数与定点数的转换关系可表示为：
        浮点数(x)转换为定点数(xq)：xq=(int)x* 2Q
        定点数(xq)转换为浮点数(x)：x=(float)xq*2-Q
    例如,浮点数x=0.5,定标Q=15,则定点数xq=L0.5*32768J=16384,式中LJ表示下取整。反之,一个用Q=15表示的定点数16384,其浮点数为163幼*2-15=16384/32768=0.5。浮点数转换为定点数时,为了降低截尾误差,在取整前可以先加上0.5。





表1.1    Q表示、S表示及数值范围
Q表示    S表示    十进制数表示范围
Q15    S0.15    -1≤x≤0.9999695
Q14    S1.14    -2≤x≤1.9999390
Q13    S2.13    -4≤x≤3.9998779
Q12    S3.12    -8≤x≤7.9997559
Q11    S4.11    -16≤x≤15.9995117
Q10    S5.10    -32≤x≤31.9990234
Q9    S6.9    -64≤x≤63.9980469
Q8    S7.8    -128≤x≤127.9960938
Q7    S8.7    -256≤x≤255.9921875
Q6    S9.6    -512≤x≤511.9804375
Q5    S10.5    -1024≤x≤1023.96875
Q4    S11.4    -2048≤x≤2047.9375
Q3    S12.3    -4096≤x≤4095.875
Q2    S13.2    -8192≤x≤8191.75
Q1    S14.1    -16384≤x≤16383.5
Q0    S15.0    -32768≤x≤32767

2  高级语言：从浮点到定点
    我们在编写DSP模拟算法时,为了方便,一般都是采用高级语言(如C语言)来编写模拟程序。程序中所用的变量一般既有整型数,又有浮点数。如例1.1程序中的变量i是整型数,而pi是浮点数,hamwindow则是浮点数组。
例1.1  256点汉明窗计算
int i;+
float pi=3.14l59;
float hamwindow[256];
for(i=0;i<256;i++)  hamwindow[i]=0.54-0.46*cos(2.0*pi*i/255);
    如果我们要将上述程序用某种足点DSP芯片来实现,则需将上述程序改写为DSP芯片的汇编语言程序。为了DSP程序调试的方便及模拟定点DSP实现时的算法性能,在编写DSP汇编程序之前一般需将高级语言浮点算法改写为高级语言定点算法。下面我们讨论基本算术运算的定点实现方法。
2.1  加法/减法运算的C语言定点摸拟
设浮点加法运算的表达式为：
float x,y,z;
z=x+y;
将浮点加法/减法转化为定点加法/减法时最重要的一点就是必须保证两个操作数的定标
temp=x+temp; 
z=temp>>(Qx-Qz),若Qx>=Qz
z=temp<<(Qz-Qx),若Qx<=Qz
例1.4结果超过16位的定点加法
设x=l5000,y=20000,则浮点运算值为z=x+y=35000,显然z>32767,因此
Qx=1,Qy=0,Qz=0,则定点加法为：
x=30000;y=20000;
temp=20000<<1=40000; 
temp=temp+x=40000+30000=70000;
z=70000L>>1=35000;
    因为z的Q值为0,所以定点值z=35000就是浮点值,这里z是一个长整型数。当加法或加法的结果超过16位表示范围时,如果程序员事先能够了解到这种情况,并且需要保持运算精度时,则必须保持32位结果。如果程序中是按照16位数进行运算的,则超过16位实际上就是出现了溢出。如果不采取适当的措施,则数据溢出会导致运算精度的严重恶化。一般的定点DSP芯片都没有溢出保护功能,当溢出保护功能有效时,一旦出现溢出,则累加器ACC的结果为最大的饱和值(上溢为7FFFH,下溢为8001H),从而达到防止溢出引起精度严重恶化的目的。
2.2乘法运算的C语言定点模拟
设浮点乘法运算的表达式为：
float x,y,z;
z=xy; 
假设经过统计后x的定标值为Qx,y的定标值为Qy,乘积z的定标值为Qz,则
z=xy
zq*2-Qx=xq*yq*2-(Qx+Qy)
zq=(xqyq)2Qz-(Qx+Qy)
所以定点表示的乘法为：
int  x,y,z;
long temp;
temp=(long)x; 
z=(temp*y)>>(Qx+Qy-Qz);
例1.5定点乘法。
设x=18.4,y=36.8,则浮点运算值为=18.4*36.8=677.12;
根据上节,得Qx=10,Qy=9,Qz=5,所以
x=18841;y=18841;
temp=18841L;
z=(18841L*18841)>>(10+9-5)=354983281L>>14=21666;
因为z的定标值为5,故定点z=21666,即为浮点的z=21666/32=677.08。
2.3除法运算的C语言定点摸拟
设浮点除法运算的表达式为：
float x,y,z;
z=x/y;
假设经过统计后被除数x的定标值为Qx,除数y的定标值为Qy,商z的定标值为Qz,则
z=x/y
zq*2-Qz=(xq*2-Qx)/(yq*2-Qy)
zq=(xq*2(Qz-Qx+Qy))/yq
所以定点表示的除法为：
int x,y,z;
long temp;
temp=(long)x;
z=(temp<<(Qz-Qx+Qy))/y;
例1.6定点除法。
设x=18.4,y=36.8,浮点运算值为z=x/y=18.4/36.8=0.5;
根据上节,得Qx=10,Qy=9,Qz=15;所以有
z=18841,y=18841;
temp=(long)18841;
z=(18841L<<(15-10+9)/18841=3O8690944L/18841=16384;
因为商z的定标值为15,所以定点z=16384,即为浮点z=16384/215=0.5。
2.4程序变量的Q值确定
    在前面几节介绍的例子中,由于x,y,z的值都是已知的,因此从浮点变为定点时Q值很好确定。在实际的DSP应用中,程序中参与运算的都是变量,那么如何确定浮点程序中变量的Q值呢？从前面的分析可以知道,确定变量的Q值实际上就是确定变量的动态范围,动态范围确定了,则Q值也就确定了。
设变量的绝对值的最大值为 max ,注意 max 必须小于或等于32767。取一个整数n,使满足
2n-1< max <2n
则有
2-Q=2-15*2n=2-(15-n)
Q=15-n
例如,某变量的值在-1至+1之间,即 max <1,因此n=0,Q=15-n=15。
    既然确定了变量的 max 就可以确定其Q值,那么变量的 max 又是如何确定的呢？一般来说,确定变量的 max 有两种方法。一种是理论分析法,另一种是统计分析法。
  1.  理论分析法
    有些变量的动态范围通过理论分析是可以确定的。例如：
(1)三角函数。y=sin(x)或y=cos(x),由三角函数知识可知, y <=1。
(2)汉明窗。y(n)=0.54一0.46cos[nπn/(N-1)],0<=n<=N-1。因为-1<=cos[2πn/(N-1)]<=1,所以0.08<=y(n)<=1.0。
(3)FIR卷积。y(n)=∑h(k)x(n-k),设∑ h(k) =1.0,且x(n)是模拟信号12位量化值,即有 x(n) <=211,则 y(n) <=211。
(4)理论已经证明,在自相关线性预测编码(LPC)的程序设计中,反射系数ki满足下列不等式： ki <1.0,i=1,2,...,p,p为LPC的阶数。
  2.  统计分析法
    对于理论上无法确定范围的变量,一般采用统计分析的方法来确定其动态范围。所谓统计分析,就是用足够多的输入信号样值来确定程序中变量的动态范围,这里输入信号一方面要有一定的数量,另一方面必须尽可能地涉及各种情况。例如,在语音信号分析中,统计分析时就必须来集足够多的语音信号样值,并且在所采集的语音样值中,应尽可能地包含各种情况。如音量的大小,声音的种类(男声、女声等)。只有这样,统计出来的结果才能具有典型性。
    当然,统计分析毕竟不可能涉及所有可能发生的情况,因此,对统计得出的结果在程序设计时可采取一些保护措施,如适当牺牲一些精度,Q值取比统计值稍大些,使用DSP芯片提供的溢出保护功能等。
2.5浮点至定点变换的C程序举例
    本节我们通过一个例子来说明C程序从浮点变换至定点的方法。这是一个对语音信号(0.3~3.4kHz)进行低通滤波的C语言程序,低通滤波的截止频率为800Hz,滤波器采用19点的有限冲击响应FIR滤波。语音信号的采样频率为8kHz,每个语音样值按16位整型数存放在insp.dat文件中。
例1.7语音信号800Hz 19点FIR低通滤波C语言浮点程序。
#include  
const int length=180
void filter(int xin[],int xout[],int n,float h[]);

static float h[19]= 
{0.01218354,-0.009012882,-0.02881839,-0.04743239,-0.04584568,
-0.008692503,0.06446265,0.1544655,0.2289794,0.257883,
0.2289794,0.1544655,0.06446265,-0.008692503,-0.04584568,
-0.04743239,-0.02881839,-0.009012882,O.01218354};
static int xl[length+20];

void filter(int xin[],int xout[],int n,float h[])
{
int i,j;
float sum;
for(i=0;i for(i=0;i＜length;i++)
{
sum=0.0;
for(j=0;j＜n;j++)sum+=h[j]*x1[i-j+n-1];
xout[i]=(int)sum;
for(i=0;i＜(n-l);i++)x1[n-i-2]=xin[length-1-i];
} 

void main()
FILE *fp1,*fp2;
int frame,indata[length],outdata[length];
fp1=fopen(insp.dat,"rb");
fp2=fopen(Outsp.dat,"wb");
frame=0;
while(feof(fp1) ==0) 
{
frame++;
printf(“frame=％d＼n”,frame);
for(i=0;i＜length;i++)indata[i]=getw(fp1);  
filter(indata,outdata,19,h);
for(i=0;i＜length;i++)putw(outdata[i],fp2);
} 
fcloseall();
return(0);
}
例1.8语音信号800Hz l9点FIR低通滤波C语言定点程序。
#include 
const int length=180; 
void  filter (int xin[],int xout[],int n,int h[]);
static int h[19]={399,-296,-945,-1555,-1503,-285,2112,5061,7503,8450,
7503,5061,2112,-285,-1503,-1555,-945,-296,399};
static int x1[length+20];

void filter(int xin[],int xout[],int n,int h[])
int i,j;
long sum;
for(i=0;i＜length;i++)x1[n＋i-111=xin][i];
for(i=0;i＜1ength;i++)
sum=0;
for(j=0;j＜n;j++)sum+=(long)h[j]*x1[i-j＋n-1];
xout[i]=sum>>15;
for(i=0;i＜(n-1);i＋＋)x1[n-i-2]=xin[length-i-1];
}
主程序与浮点的完全一样。“
3  DSP定点算术运算
    定点DSP芯片的数值表示基于2的补码表示形式。每个16位数用l个符号位、i个整数位和15-i个小数位来表示。因此：
00000010.10100000
表示的值为：
21＋2-1＋2-3=2.625
     这个数可用Q8格式(8个小数位)来表示,其表示的数值范围为-128至＋l27.996,一个Q8定点数的小数精度为1/256=0.004。
    虽然特殊情况(如动态范围和精度要求)必须使用混合表示法。但是,更通常的是全部以Q15格式表示的小数或以Q0格式表示的整数来工作。这一点对于主要是乘法和累加的信号处理算法特别现实,小数乘以小数得小数,整数乘以整数得整数。当然,乘积累加时可能会出现溢出现象,在这种情况下,程序员应当了解数学里面的物理过程以注意可能的溢出情况。下面我们来讨论乘法、加法和除法的DSP定点运算,汇编程序以TMS320C25为例。
3.1定点乘法
    两个定点数相乘时可以分为下列三种情况：
1.  小数乘小数
例1.9  Q15*Q15=Q30
0.5*0.5=0.25
0.100000000000000;Q15
  *  0.100000000000000;Q15
--------------------------------------------
00.010000000000000000000000000000=0.25;Q30
    两个Q15的小数相乘后得到一个Q30的小数,即有两个符号位。一般情况下相乘后得到的满精度数不必全部保留,而只需保留16位单精度数。由于相乘后得到的高16位不满15位的小数据度,为了达到15位精度,可将乘积左移一位,下面是上述乘法的TMS320C25程序：
LT  OP1;OP1=4000H(0.5/Q15)
MPY OP2;oP2=4000H(0.5/Ql5)
PAC
SACH  ANS,1;ANS=2000H(0.25/Q15)