线性插值。如何在 C 中实现这个算法？（给出Python版本）答案

【问题标题】：Linear Interpolation. How to implement this algorithm in C ? (Python version is given)线性插值。如何在 C 中实现这个算法？（给出Python版本）
【发布时间】：2010-12-16 10:49:53
【问题描述】：

存在一种非常好的线性插值方法。它执行线性插值，需要每个输出样本最多一次乘法。我在 Lyons 的《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 ...]

【问题讨论】：

如果你只是想在 C 中实现它以获得性能，但仍然在 Python 中使用它，我推荐Cython。
提供的 Python 代码只是为了了解这个算法是如何工作的
如果使用有意义的变量名，则更容易理解算法的工作原理。
@Lennart：已更新。现在应该很容易理解了
算法有在线参考吗？此外，避免乘法听起来像是对 80 年代的优化......

标签： python c algorithm math signal-processing

【解决方案1】：

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

...或者我称之为“上采样”（可能是错误的术语。免责声明：我没有读过里昂的文章）。我只需要了解代码的作用，然后重新编写它以提高可读性。鉴于它有几个问题：

a) 效率低下 - 两个循环是可以的，但它对每个输出项进行乘法运算；它还使用中间列表（hold），使用append（小啤酒）生成结果

b) 它在第一个区间内插错了；它在第一个元素前面生成假数据。假设我们有 multiplier=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 sn-p 是“已键入但从未编译或运行”，因此可能存在语法错误、off-by-1 错误等。但总体而言，这个想法是存在的。

【讨论】：

哦，真可惜，赏金到期了，只有 1/2 获得了奖励。剩下的 25 分——呸！ - 从 stackoverflow 领域消失。

【解决方案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)

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

【讨论】：

不，你复制粘贴了一个错误的句子。我已经谈到了保持缓冲区的作用。完整的算法产生线性插值，而不仅仅是重复项。请查看我更新的问题。
是的，这样计算的次数不会少。循环本身不会花费大量时间，因此没有必要避免它。不过，它应该可以工作。 +1 用于显示替代（如果无用）解决方案。 ;)

【解决方案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 是您想与硬件交谈的一门很棒的语言，但否则只是不必要的困难。

【讨论】：

@psihodelia: 但是 i-temp1 没有列出列表，所以不会是 [10, 20, 10]。 [i-temp1] 将始终是一个仅包含一个值的列表。
@Lennart, * L 将导致列表重复。 [1]*2 是 [1,1]。但是，是的，我注意到了同样的事情，重复的列表并不能真正帮助表达算法。
@Lennart - 我误解了你之前的评论，我没有意识到你和我指出了同样的事情。对不起。
您的新简化插值要慢得多！您在循环中有一个除法运算。所提出算法的目的是每个样本只有一次乘法。每个样本的操作数在信号处理中非常重要。
@psihodelia - 这是一个常数的除法。编译器应该看到并将其转换为乘法。如果分析或检查程序集显示它没有，那么只需声明一个临时并预先计算它。无论如何，这个算法提供了正确的结果，这似乎更重要。

【解决方案4】：

我认为您需要两个循环。您必须跨过x 中的样本来初始化插值器，更不用说将它们的值复制到y，并且您必须跨过输出样本来填充它们的值。我想您可以执行一个循环将x 复制到y 的适当位置，然后再执行另一个循环以使用y 中的所有值，但这仍然需要一些步进逻辑。最好使用嵌套循环方法。

（并且，正如Lennart Regebro 指出的那样）作为旁注，我不明白你为什么要这样做hold = [i-temp1] * L。相反，为什么不做hold = i-temp，然后循环for j in xrange(L): 和temp2 += hold？这将使用更少的内存，但其他方面的行为完全相同。

【讨论】：

请不要太关注 Python 代码。我刚刚根据可视化 LTI 系统图编写了它。这只是为了了解它是如何工作的。是的，代码中不需要使用缓冲区，但是缓冲区很容易放在图表上。

【解决方案5】：

这是我对您算法的 C 实现的尝试。在尝试进一步优化它之前，建议您在启用所有编译器优化的情况下分析其性能。

/**
 *
 * @param src caller supplied array with data
 * @param src_len len of src
 * @param steps to interpolate
 * @param dst output param needs to be of size src_len * steps
 */
float* linearInterpolation(float* src, size_t src_len, size_t steps, float* dst)
{
  float* dst_ptr = dst;
  float* src_ptr = src;
  float stepIncrement = 1.0f / steps;
  float temp1 = 0.0f;
  float temp2 = 0.0f;
  float hold;
  size_t idx_src, idx_steps;
  for(idx_src = 0; idx_src < src_len; ++idx_src)
  {
    hold = *src_ptr - temp1;
    temp1 = *src_ptr;
    ++src_ptr;
    for(idx_steps = 0; idx_steps < steps; ++idx_steps)
    {
      temp2 += hold;
      *dst_ptr = temp2 * stepIncrement;
      ++dst_ptr;
    }
  }
}

【讨论】：