线性插值.如何在C中实现此算法？(给出Python版本)

Question

存在一种非常好的线性插值方法.它执行线性插值,每个输出样本最多需要一个乘法.我在里昂的第三版理解DSP中找到了它的描述.此方法涉及一个特殊的保持缓冲区.给定要在任意两个输入样本之间插入的多个样本,它使用线性插值产生输出点.在这里,我使用Python重写了这个算法:

temp1, temp2 = 0, 0
iL = 1.0 / L
for i in x:
   hold = [i-temp1] * L
   temp1 = i
   for j in hold:
      temp2 += j
      y.append(temp2 *iL)

其中x包含输入样本,L是要插入的点数,y将包含输出样本.

我的问题是如何以最有效的方式在ANSI C中实现这样的算法,例如是否可以避免第二个循环？

注意:提供的Python代码只是为了理解该算法的工作原理.

更新:这是一个如何在Python中工作的示例:

x=[]
y=[]
hold=[]
num_points=20
points_inbetween = 2

temp1,temp2=0,0

for i in range(num_points):
   x.append( sin(i*2.0*pi * 0.1) )

L = points_inbetween
iL = 1.0/L
for i in x:
   hold = [i-temp1] * L
   temp1 = i
   for j in hold:
      temp2 += j
      y.append(temp2 * iL)

比方说x = [.... 10,20,30 ......].然后,如果L = 1,它将产生[... 10,15,20,25,30 ......]

Answer 1

"信号采样率增加"意义上的插值

...或者我称之为"上采样"(错误的术语,可能是.免责声明:我还没读过里昂斯).我只需要了解代码的作用,然后重新编写它以便于阅读.鉴于它有几个问题:

a)它是低效的 - 两个循环是可以的,但它确实为每个输出项增殖; 它也使用中间列表(hold),生成结果append(小啤酒)

b)它在第一个间隔内插错误; 它在第一个元素前面生成假数据.假设我们有乘数= 5和seq = [20,30] - 它将生成[0,4,8,12,16,20,22,24,28,30]而不是[20,22,24,26, 28,30].

所以这里是生成器形式的算法:

def upsampler(seq, multiplier):
    if seq:
        step = 1.0 / multiplier
        y0 = seq[0];
        yield y0
        for y in seq[1:]:
            dY = (y-y0) * step
            for i in range(multiplier-1):
                y0 += dY;
                yield y0
            y0 = y;
            yield y0

好的,现在进行一些测试:

>>> list(upsampler([], 3))  # this is just the same as [Y for Y in upsampler([], 3)]
[]
>>> list(upsampler([1], 3))
[1]
>>> list(upsampler([1,2], 3))
[1, 1.3333333333333333, 1.6666666666666665, 2]
>>> from math import sin, pi
>>> seq = [sin(2.0*pi * i/10) for i in range(20)]
>>> seq
[0.0, 0.58778525229247314, 0.95105651629515353, 0.95105651629515364, 0.58778525229247325, 1.2246063538223773e-016, -0.58778525229247303, -0.95105651629515353, -0.95105651629515364, -0.58778525229247336, -2.4492127076447545e-016, 0.58778525229247214, 0.95105651629515353, 0.95105651629515364, 0.58778525229247336, 3.6738190614671318e-016, -0.5877852522924728, -0.95105651629515342, -0.95105651629515375, -0.58778525229247347]
>>> list(upsampler(seq, 2))
[0.0, 0.29389262614623657, 0.58778525229247314, 0.76942088429381328, 0.95105651629515353, 0.95105651629515364, 0.95105651629515364, 0.7694208842938135, 0.58778525229247325, 0.29389262614623668, 1.2246063538223773e-016, -0.29389262614623646, -0.58778525229247303, -0.76942088429381328, -0.95105651629515353, -0.95105651629515364, -0.95105651629515364, -0.7694208842938135, -0.58778525229247336, -0.29389262614623679, -2.4492127076447545e-016, 0.29389262614623596, 0.58778525229247214, 0.76942088429381283, 0.95105651629515353, 0.95105651629515364, 0.95105651629515364, 0.7694208842938135, 0.58778525229247336, 0.29389262614623685, 3.6738190614671318e-016, -0.29389262614623618, -0.5877852522924728, -0.76942088429381306, -0.95105651629515342, -0.95105651629515364, -0.95105651629515375, -0.76942088429381361, -0.58778525229247347]

这是我对C的翻译,适合Kratz的fn模板:

/**
 *
 * @param src caller supplied array with data
 * @param src_len len of src
 * @param steps to interpolate
 * @param dst output param will be filled with (src_len - 1) * steps + 1 samples
 */
float* linearInterpolation(float* src, int src_len, int steps, float* dst)
{
    float step, y0, dy;
    float *src_end;
    if (src_len > 0) {
        step = 1.0 / steps;
        for (src_end = src+src_len; *dst++ = y0 = *src++, src < src_end; ) {
            dY = (*src - y0) * step;
            for (int i=steps; i>0; i--) {
                *dst++ = y0 += dY;
            }
        }
    }
}

请注意C片段是"已键入但从未编译或运行",因此可能存在语法错误,1个错误等错误.但总体而言,这个想法就在那里.