打印

快速计算正弦波

[复制链接]
553|19
手机看帖
扫描二维码
随时随地手机跟帖
跳转到指定楼层
楼主
FCCdsp|  楼主 | 2017-7-10 18:23 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
快速计算正弦波


在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.9960713,这个数转换为16比特的Q14的值是0x7fc0。我们可以看到这时定点小数已经十分接近0x7fff了,最终产生的正弦波的频率也会有很大的误差。为了能够精确地计算这样的正弦波,必须使用32bit的Q30定点小数。关于32bit定点小数的计算方法将在别的章节介绍。

另外上面的singen函数每调用一次只产生一个值,如果要产生实时的正弦波的话,函数的调用频率和取样频率相同,DSP的负担相对比较大。一般DSP计算都采取块计算方式,一次计算n个(例如64)个取样值,这样不但减少了函数的调用负担,也可以减少中间的内存移动的次数(y[2]=y[1];y[1]= y[0];)。


相关帖子

沙发
quickman| | 2017-7-14 09:27 | 只看该作者
最快的就是查表输出数据

使用特权

评论回复
板凳
jstgotodo| | 2017-7-14 09:28 | 只看该作者
计算的速度能够提高多少?

使用特权

评论回复
地板
iamaiqiyi| | 2017-7-14 09:28 | 只看该作者

使用特权

评论回复
5
dzfansman| | 2017-7-14 09:29 | 只看该作者
以后的程序编写也是需要靠创新的。

使用特权

评论回复
6
sanxingnote7| | 2017-7-14 09:29 | 只看该作者
主要是数学没有学习好。

使用特权

评论回复
7
backlugin| | 2017-7-14 09:30 | 只看该作者
有没有快速实现FFT的方法?

使用特权

评论回复
8
sdCAD| | 2017-7-14 09:30 | 只看该作者
这个数据计算量能够减少到多少

使用特权

评论回复
9
fengm| | 2017-7-14 09:32 | 只看该作者
定点小数

使用特权

评论回复
10
uiint| | 2017-7-14 09:32 | 只看该作者
DSP计算浮点型还是比较快的吧。

使用特权

评论回复
11
quickman| | 2017-7-14 09:32 | 只看该作者
做一个正弦波的表,快速查表。

使用特权

评论回复
12
jstgotodo| | 2017-7-14 09:33 | 只看该作者
以前都是使用sin函数自己计算得到的。

使用特权

评论回复
13
iamaiqiyi| | 2017-7-14 09:33 | 只看该作者
原来还可以这样计算正弦波函数。

使用特权

评论回复
14
dzfansman| | 2017-7-14 09:33 | 只看该作者
确实学习到了新技术,这个也可以。

使用特权

评论回复
15
sanxingnote7| | 2017-7-14 09:33 | 只看该作者
这个根本想不到怎么计算的。

使用特权

评论回复
16
backlugin| | 2017-7-14 09:33 | 只看该作者
最近做1024的FFT性能确实有待提高。

使用特权

评论回复
17
sdCAD| | 2017-7-14 09:33 | 只看该作者
使用TMS320F28335能够实现吗?

使用特权

评论回复
18
uiint| | 2017-7-14 09:33 | 只看该作者
这样的代码工作量也会提高。

使用特权

评论回复
19
fengm| | 2017-7-14 09:33 | 只看该作者
这个的思路想不很明白,是什么原理?

使用特权

评论回复
20
jcx0324| | 2017-7-18 08:47 | 只看该作者
看看,内容没了

使用特权

评论回复
发新帖 我要提问
您需要登录后才可以回帖 登录 | 注册

本版积分规则

967

主题

1447

帖子

9

粉丝