用递推数列计算正弦波
在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 ) 1式
我们看到这个递推公式是:
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);由1式可知y(n-2)=0,则
float a1=2*cos(w*Ts); //y(n-1)必须等于sin(w*Ts),这里的y(n)初值随意
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=30831=0x786E同样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.9960713，这个数转换为16比特的Q14的值是0x7fc0。我们可以看到这时定点小数已经十分接近0x7fff了，最终产生的正弦波的频率也会有很大的误差。为了能够精确地计算这样的正弦波，必须使用32bit的Q30定点小数。关于32bit定点小数的计算方法将在别的章节介绍。
另外上面的singen函数每调用一次只产生一个值，如果要产生实时的正弦波的话，函数的调用频率和取样频率相同，DSP的负担相对比较大。一般DSP计算都采取块计算方式，一次计算n个（例如64）个取样值，这样不但减少了函数的调用负担，也可以减少中间的内存移动的次数（y[2]=y[1];y[1]=y[0];）。
DSP基础--定点小数运算
许多DSP芯片只支持整数运算，如果现在这些芯片上进行小数运算的话，定点小数运算应该是最佳选择了，此外即使芯片支持浮点数，定点小数运算也是最佳的速度选择。
在DSP世界中，由于DSP芯片的限制,经常使用定点小数运算。所谓定点小数，实际上就是用整数来进行小数运算。下面先介绍定点小数的一些理论知识，然后以C语言为例，介绍一下定点小数运算的方法。在TI C5000 DSP系列中使用16比特为最小的储存单位，所以我们就用16比特的整数来进行定点小数运算。
先从整数开始，16比特的储存单位最多可以表示0x0000到0xffff，65536种状态，如果它表示C语言中的无符号整数的话，就是从0到65535。如果需要表示负数的话，那么最高位就是符号位，而剩下的15位可以表示32768种状态。这里可以看出，对于计算机或者DSP芯片来说，符号并没有什么特殊的储存方式，其实是和数字一起储存的。为了使得无论是无符号数还是符号数，都可以使用同样的加法减法规则，符号数中的负数用正数的补码表示。
我们都知道-1 + 1 =0，而0x0001表示1，那么-1用什么来表示才能使得-1 + 1 =0呢？答案很简单：0xffff。现在就可以打开Windows的计算器，用16进制计算一下0xffff+0x0001，结果是0x10000。那么0x10000和0x0000等价麽，我们刚才说过用16比特来表达整数，最高位的1是第17位，这一位是溢出位，在运算寄存器中没有储存这一位，所以结果是低16位，也就是0x0000。现在我们知道负数的表达方式了。举个例子：-100。首先我们需要知道100的16进制，用计算器转换一下，可以知道是0x0064，那么-100就是0x10000 - 0x0064，用计算器算一下得0xff9c。
还有一种简单的转换符号的方法，就是取反加一：把数x写成二进制格式，每位0变1，1变0，最后把结果加1就是-x了。
好，复习了整数的相关知识之后，我们进入定点小数运算环节。所谓定点小数，就是小数点的位置是固定的。我们是要用整数来表示定点小数，由于小数点的位置是固定的，所以就没有必要储存它（如果储存了小数点的位置，那就是浮点数了）。既然没有储存小数点的位置，那么计算机当然就不知道小数点的位置，所以这个小数点的位置是我们写程序的人自己需要牢记的。
先以10进制为例。如果我们能够计算12+34=46的话，当然也就能够计算1.2+3.4 或者 0.12+0.34了。所以定点小数的加减法和整数的相同，并且和小数点的位置无关。乘法就不同了。 12*34=408，而1.2*3.4=4.08。这里1.2的小数点在第1位之前，而4.08的小数点在第2位之前，小数点发生了移动。所以在做乘法的时候，需要对小数点的位置进行调整？！可是既然我们是做定点小数运算，那就说小数点的位置不能动！！怎么解决这个矛盾呢，那就是舍弃最低位。 也就说1.2*3.4=4.1，这样我们就得到正确的定点运算的结果了。所以在做定点小数运算的时候不仅需要牢记小数点的位置，还需要记住表达定点小数的有效位数。上面这个例子中，有效位数为2，小数点之后有一位。
现在进入二进制。我们的定点小数用16位二进制表达，最高位是符号位，那么有效位就是15位。小数点之后可以有0 - 15位。我们把小数点之后有n位叫做Qn，例如小数点之后有12位叫做Q12格式的定点小数，而Q0就是我们所说的整数。
Q12的正数的最大值是 0 111 . 111111111111，第一个0是符号位，后面的数都是1，那么这个数是十进制的多少呢，很好运算，就是 0x7fff / 2^12 = 7.999755859375。对于Qn格式的定点小数的表达的数值就它的整数值除以2^n。在计算机中还是以整数来运算，我们把它想象成实际所表达的值的时候，进行这个运算。反过来把一个实际所要表达的值x转换Qn型的定点小数的时候，就是x*2^n了。例如 0.2的Q12型定点小数为：0.2*2^12 = 819.2，由于这个数要用整数储存， 所以是819 即 0x0333。因为舍弃了小数部分，所以0x0333不是精确的0.2，实际上819/2^12=0.199951171875
我们用数学表达式做一下总结：
x表示实际的数（*一个浮点数）， q表示它的Qn型定点小数（一个整数）。
q = (int) (x * 2^n)
x = (float)q/2^n
由以上公式我们可以很快得出定点小数的+-*/算法：
假设q1，q2，q3表达的值分别为x1，x2，x3
q3 = q1 + q2 若 x3 = x1 + x2
q3 = q1 - q2 若 x3 = x1 - x2
q3 = q1 * q2 / 2^n若 x3 = x1 * x2
q3 = q1 * 2^n / q2若 x3 = x1 / x2
我们看到加减法和一般的整数运算相同，而乘除法的时候，为了使得结果的小数点位不移动，对数值进行了移动。
用c语言来写定点小数的乘法就是：
short q1,q2,q3;
....
q3=((long q1) * (long q2)) >> n;
由于/ 2^n和* 2^n可以简单的用移位来计算，所以定点小数的运算比浮点小数要快得多。下面我们用一个例子来验证一下上面的公式：
用Q12来计算2.1 * 2.2，先把2.1 2.2转换为Q12定点小数：
2.1 * 2^12 = 8601.6 = 8602
2.2 * 2^12 = 9011.2 = 9011
(8602 * 9011) >> 12 = 18923
18923的实际值是18923/2^12 = 4.619873046875 和实际的结果 4.62相差0.000126953125，对于一般的计算已经足够精确了。
