线性插值。如何在C语言中实现该算法？（给出Python版本）

Question

有一种很好的线性插值方法。它执行线性插值，每个输出样本最多需要一个乘法。我在里昂的《理解数字信号处理器》第三版中找到了它的描述。此方法涉及一个特殊的保持缓冲区。给定要在任意两个输入样本之间插入的样本数，它使用线性插值生成输出点。在这里，我用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，序列为[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代码片段是“类型化的，但从未编译或运行”，因此可能会有语法错误、off-by-1错误等，但总的来说是这样的。

Answer 2

在这种情况下，我认为您可以避免第二个循环：

def interpolate2(x, L):
    new_list = []
    new_len = (len(x) - 1) * (L + 1)
    for i in range(0, new_len):
        step = i / (L + 1)
        substep = i % (L + 1)
        fr = x[step]
        to = x[step + 1]
        dy = float(to - fr) / float(L + 1)
        y = fr + (dy * substep)
        new_list.append(y)
    new_list.append(x[-1])
    return new_list

print interpolate2([10, 20, 30], 3)

你只需直接计算你想要的位置上的成员。不过，这可能不是最有效的方法。唯一确定的方法是编译它，看看哪个更快。

Answer 3

首先，你的代码坏了。L没有定义，y和x也没有

一旦解决了这个问题，我就对生成的代码运行cython：

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

这似乎奏效了。不过，我还没有尝试编译它，您还可以通过添加不同的优化来大大提高速度。

“例如，是否可以避免第二个循环？”

如果是，那么在Python中也是可能的。我不知道怎么做，虽然我不知道你为什么会这样做。首先创建一个长度为L的i-temp列表是完全没有意义的。循环L次：

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

但对于你所得到的一切来说，似乎都太复杂了。你到底想干什么？插入什么？（但它在标题中是这么说的。很抱歉。）

当然有更简单的插值方法。

更新，一个非常简化的插值函数：

# A simple list, so it's easy to see that you interpolate.
indata = [float(x) for x in range(0, 110, 10)]
points_inbetween = 3

outdata = [indata[0]]

for point in indata[1:]: # All except the first
    step = (point - outdata[-1]) / (points_inbetween + 1)
    for i in range(points_inbetween):
        outdata.append(outdata[-1] + step)

我看不到摆脱内环的方法，也找不到想要这样做的理由。 把它转换成C，我会留给其他人，或者更好的，Cython，因为C是一个很好的语言，你们想和硬件交流，但在其他方面只是不必要的困难。