sin定点化

sin定点化 默认分类 2010-01-19 15:26 阅读4 评论0 字号： 大 中 小 sincospowtan 1. sin定点化 转自某处：用递推数列计算sin ---;month=04 本文介绍快速计算正弦波的方法，并且就速度和精度问题进行一些讨论。 在DSP运用中，经常需要产生正弦波。如果 直接用c的数学函数sin，当然可以产生正弦波，但是由于sin函数本身的效率很低，产生正弦波所需要的MIPS就会占去DSP处理能力的相当大的一部 分。本章介用递推数列算正弦波的 方法，先介绍原理，推导出递推，然后用浮点小数实现计算，再用定点小数 进一步优化算法，最后进行误差，并提出更精 确的定点小数算法。 先来看看如何推导出递推数列的公式。 我们所要产生的正弦波，其实是一系列的整数，把这些整数按照一定的取样频率 发送给数模转换器，就可以变成真正的正弦波了。假设取样周期是Ts，产生的正弦波的圆频率为w，那么我们需要产生的数列就是： sin(0), sin(w*Ts), sin(2*w*Ts), ... sin(n*w*Ts) 假设f(n)= sin(n*w*Ts)，则问题就变成，从f(n-1), f(n-2), f(n-3),..., 如何计算f(n) 了。解决了这个问题，也就找到了递推公式。 下面是这个递推公式的求解过程，假设x=w*Ts： 公式：sin( a + b) = sin(a)*cos(b) + cos(a)*sin(b) sin( x + (n-1)x ) = sin(x)*cos( (n-1)x ) + cos(x)*sin( (n-1)x ) 公式：Sin(a)*cos(b) = 1/2 * [ sin( a+b ) + sin( a-b ) ] sin(x)*cos( (n-1)x ) = 1/2 * [ sin(nx) - sin( (n-2)x ) ] sin(nx) = 1/2 * [ sin(nx) - sin( (n-2)x ) ] + cos(x)*sin( (n-1)x ) sin(nx) = 2*cos(x)*sin( (n-1)x ) - sin( (n-2)x ) 我们看到这个递推公式是: f(n)=2*cos(w*Ts)*f(n-1) - f(n-2) 也就是说只要知道最初始的两项f(0)和f(1)，就可以计算出整个正弦波了。 根据上面的递推公式，很容易写出下面的正弦波计算程序。只要事先计算一次sin(w*Ts)和cos(w*Ts)，以后的值就可以通过递推公式得到，所以计 算一个值所需要的工作就是一次乘法，一次加法，两次变量复制而已了。 float y[3] = {0, sin(w*Ts),0}; // y(n), y(n-1), y(n-2) float a1=2*cos(w*Ts); float a2=-1; float singen(){ y[0]=a1*y[1]+a2*y[2]; y[2]=y[1]; y[1]=y[0]; return y[0]; } 假如我们需要产生取样频率为8KHz的440Hz的正弦波，那么a1=2*cos(2*pi*440/8000)=1.8817615，而y[1]=sin(2*pi*440/8000)=0.33873792。 现在看如何用定点小数来更快的计算正弦波。我们使用16bit也就是short型的整数来表示定点小数。首先需要决定的是小数的Q值，虽然我们最 后计算的正弦波的值都是小于1的，但是在计算过程中需要用2*cos(w*Ts)，而这个值最大为2，所以我们选择的Q值必须至少最大能表示2。这 里我们选择Q14，Q14的定点小数能表示-2到2的取值范围，对于本例的正弦波计算正好合适。1.8817615的Q14值是1.8817615*2^14=5550= 0x786F，同样0.33873792的Q14值为0x15AE。 下面就是完整的计算8KHz取样频率的400Hz的定点小数的正弦波的程序。 short y[3] = {0, 0x15AE,0}; // y(n), y(n-1), y(n-2) short a1=0x786F; short a2=0xC000; short singen(){ y[0]=( (long)a1*(long)y[1]+(long)a2*(long)y[2] )>>14; y[2]=y[1]; y[1]=y[0]; return y[0]; } 使用定点小数计算不但速度比浮点更快，而且计算得出来的值是整数，这个数值 可以直接传递给DAC(数模转换器)转换为模拟的声音信号，如果使用浮点小数计 算的话，还必须把浮点数转换为整数才能传递给DAC。 使用定点小数计算必须仔细分析误差，下面来看看我们产生的正弦波的误差是多 少。定点小数计算中的误差就是由定点小数表达精度决定的。在上面的例子中我 们用0x786F表示1.8817615，这存在一定的误差，把Q14的0x786F再转换为浮点数就是0x786F/2^14=1.8817749，我们可以看到相对误差非常小，也就是 说最终得到的正弦波在频率上的误差也是非常小的。 但是，定点小数并不是什么时候都这么精确。例如如果用CD音质的取样频率44100Hz来产生100Hz的正弦波，那么a1=2*cos(2*pi*440/44100)= 1.996071 3，这个数转换为16比特的Q14的值是0x7fc0。我们可以看到这时定点小数已 经十分接近0x7fff了，最终产生的正弦波的频率也会有很大的误差。为了能够精 确地计算这样的正弦波，必须使用32bit的Q30定点小数。关于32bit定点小数的计算方法将在别的章节介绍。 另外上面的singen函数每调用一次只产生一个值，如果要产生实时的正弦波的 话，函数的调用频率和取样频率相同，DSP的负担相对比较大。一般DSP计算 都采取块计算方式，一次计算n个（例如64）个取样值，这样不但减少了函数 的调用负担，也可以减少中间的内存移动的次数（y[2]=y[1];y[1]=y[0];）。 2. 开方定点化程序---未验证 , #include , #include , #define max(a, b) (((a) > (b)) ? (a) : (b)) , #define min(a, b) (((a) < (b)) ? (a) : (b)) , #define maxnum 0x7fffffff , #define testbits 0x80000000 , , , int ComputeSqrt(int Data,int DataLeftShift,int *OutData, int *OutDataShift) , { , int i,temp1=0,temp1shift=0,temp2=0,divshift,quotient,n=0; , if (Data==0) , { , *OutData=0;*OutDataShift=0; , return 0; , } , temp1 = testbits; , for (i=0;i<31;i++) , { , if(Data&temp1) , break; , Data=Data<<1; , DataLeftShift++; , } , temp1 = 1<<30; , temp1shift = 15; , *OutDataShift = DataLeftShift>>1; , Data = Data>>(DataLeftShift-(*OutDataShift<<1)); , *OutDataShift = *OutDataShift+15; , while (1)//循环迭代时，将精度定为Q15.开方程序精度可以达到8位十进 制有效数字 , { , , quotient=compdivide(Data,temp1,0,temp1shift,&divshift); , if (divshift>=15) , quotient = quotient>>(divshift-15); , , temp2 = (temp1>>1)+(quotient>>1); , if ((abs(temp1-temp2)<=1)||(temp2==0)||(++n==20)) , break; , temp1=temp2; , } , *OutData =temp2; , return 1; , } , , , main() , { , int i=8,qrt,shift,sh=0; , , computesqrt(i,sh,&qrt,&shift); , printf("%d %f %d %d\n",i,((double)qrt/pow(2,shift)),qrt,shift); , }