[Math]Audio开发中定点的sin函数实现和评测

Audio开发中定点的sin函数实现和评测

-- By Water

在音频开发中经常会使用到sin这个函数，例如产生一个1KHz的纯音就会用到sin(2*pi*fc*n/fs)。但是在嵌入式系统或者linux kernel开发或者DSP开发中经常会有Fixed point开发需求，这样就不呢使用系统浮点函数库，就得自己来实现fixed point sin函数了。

最简单直白的方式就是创建一个sin table在使用的时候去查表获得相应的sin值。

产生sin表格的方式参考如下：

    // Step1. 产生需要的Sin Table
    {
        double phase = 0;
        double delta = 2*3.1415926535897932384626433832795/sin_tab_size;

        for (i=0; i<sin_tab_size; i++)
        {
            psintab[i] = (long)((0x7FFFFFFF * sin(phase)) + 0.5);
            phase += delta;
        }
    }

然后在使用的时候就可以直接查表获得sin值：

long lib_sin_lite(unsigned long phase)
{
    return tab_sine_1024_32bits_Q31[phase & 0x3FF];
}

这样做的好处是效率比较高，但精度较低。

在Table Size为1024的时候，最大的误差大约在-44dB左右。

如果要获得更高的精度就需要增大Table Size，每增大一倍就能减少6dB的误差。

在增大到4096的时候最大误差减少到-56dB左右。

而且在Table Size比较小的时候每两个点之间的Phase差会比较大。

例如：1024个点的时候，每两个点之间的phase差是0.006。

但在48KHz 采样率下要达到1Hz的精度就需要Phase差小于2*pi/48000 = 0.000131。

如果是Hi-Res audio sampling rate达到384000Hz的话要分辨1Hz的audio就需要Phase差小于0.0000164。

这样1024个点的Sin Table就远远不能满足要求了，就需要524288个点这么大的Sin Table了。

虽然目前嵌入式系统资源不像以前使用8032时候那么苛刻，但memory在运算量能保证情况下还是能省尽量省。

这个时候最简单用小Table获得更高精度的办法就是线性插值。在已知两个点之间的phase采用线性计算获得。

例如还是1024的Table，但是Phase可以用22bits表达，这样高10bits的值通过查表获得，低12bits就用于线性计算。

long lib_sin(long phase)
{
    int sflag=0;
    const long *pstab;
    long sa, sb, da, db;

    if (phase < 0)
    {
        phase = -phase;
        sflag = -1;
    }

    da = phase & 0xFFF;
    db = 0x1000 - da;
    phase >>= 12;
    pstab = tab_sine_1024_32bits_Q19 + phase;
    sa = *pstab++;
    sb = *pstab;

    if (sflag < 0)
    {
        return -(sa*db + sb*da);
    }

    return (sa*db + sb*da);
}

这样就能把1024 Table的最大误差缩小到-104.7dB。两个点之间的Phase差为0.0000015 < 0.0000164 完全满足384Khz以及一般音频的应用。

相关SourceCode可以参考：

http://git.oschina.net/webwater/AudioKits/tree/master/src/math

或者直接Git下载：https://git.oschina.net/webwater/AudioKits.git