线性插值。如何在 C 中实现这个算法？（给出Python版本）
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/4459930/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Linear Interpolation. How to implement this algorithm in C ? (Python version is given)
提问by psihodelia
There exists one very good linear interpolation method. It performs linear interpolation requiring at most one multiply per output sample. I found its description in a third edition of Understanding DSP by Lyons. This method involves a special hold buffer. Given a number of samples to be inserted between any two input samples, it produces output points using linear interpolation. Here, I have rewritten this algorithm using Python:
存在一种非常好的线性插值方法。它执行线性插值，每个输出样本最多需要一次乘法。我在 Lyons 的《理解 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)
where x contains input samples, L is a number of points to be inserted, y will contain output samples.
其中 x 包含输入样本，L 是要插入的点数，y 将包含输出样本。
My question is how to implement such algorithm in ANSI C in a most effective way, e.g. is it possible to avoid the second loop?
我的问题是如何以最有效的方式在 ANSI C 中实现这种算法，例如是否可以避免第二个循环？
NOTE: presented Python code is just to understand how this algorithm works.
注意：提供的 Python 代码只是为了了解该算法的工作原理。
UPDATE: here is an example how it works in 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)
Let's say x=[.... 10, 20, 30 ....]. Then, if L=1, it will produce [... 10, 15, 20, 25, 30 ...]
假设 x=[.... 10, 20, 30 ....]。然后，如果 L=1，它将产生 [... 10, 15, 20, 25, 30 ...]
采纳答案by Nas Banov
Interpolation in the sense of "signal sample rate increase"
“信号采样率增加”意义上的插值
... or i call it, "upsampling" (wrong term, probably. disclaimer: i have not read Lyons'). I just had to understand what the code does and then re-write it for readability. As given it has couple of problems:
......或者我称之为“上采样”（可能是错误的术语。免责声明：我没有读过里昂的书）。我只需要了解代码的作用，然后重新编写它以提高可读性。鉴于它有几个问题：
a) it is inefficient - two loops is ok but it does multiplication for every single output item; also it uses intermediary lists(hold), generates result with append(small beer)
a）效率低下 - 两个循环是可以的，但它对每个输出项进行乘法运算；它也使用中间列表（hold），用append（小啤酒）生成结果
b) it interpolates wrong the first interval; it generates fake data in front of the first element. Say we have multiplier=5 and seq=[20,30] - it will generate [0,4,8,12,16,20,22,24,28,30] instead of [20,22,24,26,28,30].
b）它在第一个间隔内插错了；它在第一个元素前面生成假数据。假设我们有 multiplier=5 和 seq=[20,30] - 它会生成 [0,4,8,12,16,20,22,24,28,30] 而不是 [20,22,24,26, 28,30]。
So here is the algorithm in form of a generator:
所以这里是生成器形式的算法：
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
Ok and now for some tests:
好的，现在进行一些测试：
>>> 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]
And here is my translation to C, fit into Kratz's fn template:
这是我对 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;
}
}
}
}
Please note the C snippet is "typed but never compiled or run", so there might be syntax errors, off-by-1 errors etc. But overall the idea is there.
请注意，C 代码片段是“已输入但从未编译或运行”，因此可能存在语法错误、off-by-1 错误等。但总体思路是存在的。
回答by Lennart Regebro
Well, first of all, your code is broken. L is not defined, and neither is y or x.
好吧，首先，您的代码已损坏。L 未定义，y 或 x 也未定义。
Once that is fixed, I run cython on the resulting code:
一旦修复，我在生成的代码上运行 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)
And that seemed to work. I haven't tried to compile it, though, and you can also improve the speed a lot by adding different optimizations.
这似乎奏效了。不过我没有尝试编译它，你也可以通过添加不同的优化来提高速度。
"e.g. is it possible to avoid the second loop?"
“例如，是否可以避免第二个循环？”
If it is, then it's possible in Python too. And I don't see how, although I don't see why you would do it the way you do. First creating a list of L length of i-temp is completely pointless. Just loop L times:
如果是，那么在 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)
It all seems overcomplicated for what you get out though. What are you trying to do, actually? Interpolate something? (Duh it says so in the title. Sorry about that.)
不过，对于你得到的东西来说，这一切似乎都过于复杂了。你到底想做什么？插值什么？（呃，标题上是这么说的。抱歉。）
There are surely easier ways of interpolating.
肯定有更简单的插值方法。
Update, a much simplified interpolation function:
更新，一个大大简化的插值函数：
# 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)
I don't see a way to get rid of the inner loop, nor a reason for wanting to do so. Converting it to C I'll leave up to someone else, or even better, Cython, as C is a great langauge of you want to talk to hardware, but otherwise just needlessly difficult.
我没有看到摆脱内部循环的方法，也没有想要这样做的理由。将它转换为 C 我会留给其他人，甚至更好，Cython，因为 C 是一种很好的语言，你想与硬件交谈，但否则就不必要地困难。
回答by mtrw
I think you need the two loops. You have to step over the samples in xto initialize the interpolator, not to mention copy their values into y, and you have to step over the output samples to fill in their values. I suppose you could do one loop to copy xinto the appropriate places in y, followed by another loop to use all the values from y, but that will still require some stepping logic. Better to use the nested loop approach.
我认为你需要两个循环。您必须跳过样本x以初始化插值器，更不用说将它们的值复制到 中y，并且必须跳过输出样本以填充它们的值。我想你可以做一个循环来复制x到 中的适当位置y，然后是另一个循环来使用 中的所有值y，但这仍然需要一些步进逻辑。最好使用嵌套循环方法。
(And, as Lennart Regebropoints out) As a side note, I don't see why you do hold = [i-temp1] * L. Instead, why not do hold = i-temp, and then loop for j in xrange(L):and temp2 += hold? This will use less memory but otherwise behave exactly the same.
（而且，正如Lennart Regebro指出的那样）作为旁注，我不明白你为什么这样做hold = [i-temp1] * L。相反，为什么不先做hold = i-temp，然后循环for j in xrange(L):和temp2 += hold？这将使用更少的内存，但其他方面的行为完全相同。
回答by Ignas Mikalajūnas
In that case I think you can avoid the second loop:
在这种情况下，我认为您可以避免第二个循环：
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)
you just calculate the member in the position you want directly. Though - that might not be the most efficient to do it. The only way to be sure is to compile it and see which one is faster.
您只需直接计算您想要的位置的成员。虽然 - 这可能不是最有效的。唯一确定的方法是编译它并查看哪个更快。
回答by Kratz
Heres my try at a C implementation for your algorithm. Before trying to further optimize it id suggest you profile its performance with all compiler optimizations enabled.
这是我为您的算法尝试的 C 实现。在尝试进一步优化它之前，id 建议您在启用所有编译器优化的情况下分析其性能。
/**
*
* @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;
}
}
}
