sinx/x插值

将离散的抽样信通过一个理想的低通滤波器就能恢复原始波形，

理想低通滤波器的冲阶响应为

sin(PI * t / T) / (PI *t / T)

其中PI指圆周率，t指时间，T则是抽样时间间隔

滤波器的输出为

∞

∑ƒ(mT)*sin(PI * (t - mT) / T) / (PI * (t - mT) / T) 

m=-∞

其中PI是圆周率，t是时间，m是抽样数组的索引，这里是一个无穷数组，T是抽样的时间间隔

但我们实际上的信号不能无穷，所以我们只使用过去M个抽样作为参考，虽然这样会产生一定的误差

M-1

∑ƒ(mT)*sin(PI * (t - mT) / T) / (PI * (t - mT) / T) = ƒ(t)

m=0

 ƒ(t)指的是插值之后的数据。ƒ(mT)  也就是原始数据

按照公式，使用所有原始数据作为参考进行插值， 编写代码如下：

#include <stdio.h>
#include <stdint.h>
#include <math.h>
#define PI 3.1415926                       //圆周率

#define BUF_LEN  64                        //原始数据长度
#define INTER_LEN 1024                     //插值以后的数据长度
#define INTER    (INTER_LEN / BUF_LEN)     //插值比例
int main()
{
    double buf[BUF_LEN];                   //原始数据
    //初始化为一个周期正弦波形，并打印
    for(int32_t i = 0; i < BUF_LEN; i++)
    {
        buf[i] = sin(PI * 2 * i / BUF_LEN);
        printf("%f,",buf[i]);
        for(int j = 0; j < (buf[i] + 1) * 64; j++)
            printf(" ");
        printf("*\n");
    }

    printf("\n##############################\n");
    double buf_inter[INTER_LEN];          //插值结果
    for(int t = 0; t < INTER_LEN; t++)
    {
        if(t % INTER == 0)                //原始数据不需要运算
            buf_inter[t] = buf[t / INTER];
        else
        {
            double sum = 0;
            for(int m = 0; m < BUF_LEN; m++)
            {
                double x = PI * ((double)t - m * INTER) / INTER;
                sum += (buf[m] * sin(x) / x);
            }
            buf_inter[t] = sum;
        }
    }
    for(int t = 0; t < INTER_LEN; t++)
    {
        printf("%d,%f,",t,buf_inter[t]);
        for(int i = 0; i< (buf_inter[t] + 1)* 64; i++)
            printf(" ");
        printf("*\n");
    }
    return 0;
}

 使用整数运算，并将系数保存

#include <stdio.h>
#include <stdint.h>
#include <math.h>
#define PI 3.1415926                            //圆周率

#define REF_LEN        64                        //参考数据长度
#define INTER        16                        //插值比例
int32_t sinx[INTER * REF_LEN][REF_LEN];
void sinx_init()
{
    for(int t = 0; t < INTER * REF_LEN; t++)
    {
        for(int m = 0; m < REF_LEN; m++)
        {
            int temp = t - m * INTER;
            if(temp == 0)
            {
                sinx[t][m] = 0x7fff;
            }
            else
            {
                double x = PI * ((double)temp) / INTER;
                sinx[t][m] = 0x7fff * sin(x) / x;
            }
        }
    }
}
void sinx_do(int32_t ref[], int32_t result[])
{
    int off = 0;
    for(int t = 0; t < INTER * REF_LEN; t++)
    {
        int32_t sum = 0;
        for(int m = 0; m < REF_LEN; m++)
            sum += ref[m] * sinx[t][m];
        result[off] = (sum / 0x7fff);
        off++;
    }
}

int main()
{
    int32_t buf[REF_LEN];                        //原始数据
    for(int32_t i = 0; i < REF_LEN; i++)                //初始化为一个周期正弦波形，并打印
    {
        buf[i] = 0x7fff * sin(PI * 2 * i / REF_LEN);
        printf("%d,",buf[i]);
        for(int j = 0; j < (buf[i] >> 9) + 64; j++)
            printf(" ");
        printf("*\n");
    }
    printf("\n##############################\n");

    sinx_init();

    int32_t buf_inter[REF_LEN * INTER ];
    sinx_do(buf, buf_inter);

    for(int t = 0; t < REF_LEN * INTER; t++)
    {
        printf("%d,%d,",t,buf_inter[t]);
        for(int i = 0; i< (buf_inter[t] >> 9) + 64; i++)
            printf(" ");
        printf("*\n");
    }
    return 0;
}