Linear interpolation. How to implement this algorithm in C? (Python version provided)

Question

There is one very good linear interpolation method. It performs linear interpolation, requiring no more than one multiplication by the output sample . I found his description in the third edition of "Understanding DSP" by Lyons. This method includes a special hold buffer. Given the number of samples to be inserted between any two input samples, it creates output points using linear interpolation. Here I rewrote this algorithm using 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 the number of points to be inserted, y will contain output samples.

My question is how to implement such an algorithm in ANSI C in the most efficient way , for example. Is it possible to avoid the second cycle?

NOTE: the presented Python code is just to understand how this algorithm works.

UPDATE: here is an example of how it works in 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 ...]

Answer 1

Interpolation in the sense of increasing the sampling rate of a signal

... or I call it "upsampling" (the wrong term, perhaps a disclaimer: I did not read "Lyon"). I just needed to understand what the code was doing and then rewrite it for readability. Since this has a couple of problems:

a) it is inefficient - two loops in order, but this is a multiplication for each individual output element; also it uses intermediate lists ( hold ), generates the result with append (small beer)

b) it interpolates an incorrect first interval; it generates fake data before the first element. Let's say we have a factor = 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].

So, here is the algorithm in the 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] , 0.58778525229247325, 1.2246063538223773e- >>> 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] , -2.4492127076447545e- >>> 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] , 3.6738190614671318e- >>> 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 into C, fits into the Kratz fn template:

 /** * * @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; } } } } 

Note that C snippet is “typed, but never compiled or launched,” so syntax errors, one-by-one errors, etc. may occur. But in general, there is an idea.

Answer 2

In this 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 compute the member in the place you want directly. Although - this may not be the most effective for this. The only way to make sure that you need to compile it and see which one is faster.

Answer 3

Well, firstly, your code is broken. L is undefined, and neither y nor x.

Once this is fixed, I ran cython for the resulting code:

 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 it seemed to work. However, I did not try to compile it, and you can also significantly increase the speed by adding various optimizations.

"for example, is it possible to avoid the second cycle?"

If so, then this is also possible in Python. And I do not understand how, although I do not understand why you will do it the way you do. At first, creating an L-length list of i-temp is completely pointless. Only cycle L times:

 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 too complicated for you to go out. What are you really trying to do? Interpolate something? (Spirit says so in the title. Sorry.)

There are, of course, simpler interpolation methods.

Update, greatly simplified interpolation function:

 # A simple list, so it 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 do not see a way to get rid of the inner cycle, and I do not want to do this. Converting it to C, I will leave it only to someone, or even better, Keaton, since C is a great tool that you want to talk to equipment, but otherwise it is simply useless difficult.

Answer 4

I think you need two loops. You must step over the patterns in x to initialize the interpolator, not to mention copy their values ​​into y , and you need to go to the output patterns to fill in their values. I suppose you could do one loop to copy x to the appropriate places in y , and then another loop to use all the values ​​from y , but it still needs some step logic. Better use a nested loop approach.

(And, as Lennart Regebro points out). As a side note, I don't understand 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 will behave exactly the same.

Answer 5

Here is my attempt to implement C for your algorithm. Before trying to optimize it, id will offer you to profile its performance with all compiler optimizations.

 /** * * @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; } } }