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

Question

有一种非常好的线性插值方法。它在每个输出样本的计算中，最多只需要进行一次乘法运算。我在《Understanding 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

首先，你的代码有问题。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语言在和硬件打交道时很棒，但在其他情况下就显得不必要地复杂。

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

信号采样率增加的插值

...或者我称之为“上采样”（这个词可能用错了。声明一下：我没有读过Lyons的书）。我只是想搞明白这段代码在干嘛，然后把它重新写得更容易理解。原来的代码有几个问题：

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的函数模板：

/**
 *
 * @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代码是“写了但从未编译或运行过”，所以可能会有语法错误、越界错误等。但总体思路是对的。