Линейная интерполяция. Как реализовать этот алгоритм в C? (Предоставляется версия Python) - PullRequest
Новая
       12

Линейная интерполяция. Как реализовать этот алгоритм в C? (Предоставляется версия Python)

4 голосов
/ 16 декабря 2010

Существует один очень хороший метод линейной интерполяции. Он выполняет линейную интерполяцию, требующую не более одного умножения на выходную выборку . Я нашел его описание в третьем издании Understanding DSP by Lyons. Этот метод включает специальный буфер удержания. Учитывая количество выборок, которые будут вставлены между любыми двумя входными выборками, это производит выходные точки, используя линейную интерполяцию. Здесь я переписал этот алгоритм с использованием 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)

Скажем, х = [.... 10, 20, 30 ....]. Тогда, если L = 1, будет получено [... 10, 15, 20, 25, 30 ...]

Ответы [ 5 ]

8 голосов
/ 27 декабря 2010

Интерполяция в смысле «увеличения частоты дискретизации сигнала»

... или я называю это "повышающей дискретизацией" (неправильный термин, вероятно. Отказ от ответственности: я не читал "Лиона"). Мне просто нужно было понять, что делает код, а затем переписать его для удобства чтения. Как дано, у него есть пара проблем:

а) это неэффективно - два цикла в порядке, но они умножают для каждого элемента вывода; также он использует промежуточные списки (hold), генерирует результат с append (мелочь)

б) неправильно интерполирует первый интервал; он генерирует поддельные данные перед первым элементом. Скажем, у нас есть множитель = 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, вписывающийся в шаблон 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-1 и т. Д. Но в целом идея существует.

2 голосов
/ 16 декабря 2010

В этом случае, я думаю, вы можете избежать второго цикла: 

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 голос
/ 16 декабря 2010

Ну, во-первых, ваш код не работает.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 - отличный язык для вас, желающих поговорить с аппаратным обеспечением, но в остальном просто неоправданно сложно.

0 голосов
/ 26 декабря 2010

Вот моя попытка реализации 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;
    }
  }
}
0 голосов
/ 16 декабря 2010

Я думаю, вам нужны две петли.Чтобы инициализировать интерполятор, вам нужно перешагнуть через сэмплы в x, не говоря уже о копировании их значений в y, и вам нужно перешагнуть через выходные сэмплы, чтобы заполнить их значения.Я полагаю, вы могли бы сделать один цикл, чтобы скопировать x в соответствующие места в y, а затем еще один цикл, чтобы использовать все значения из y, но это все равно потребует некоторой логики шага.Лучше использовать подход с вложенными циклами.

(И, как указывает Lennart Regebro ) В качестве примечания, я не понимаю, почему вы делаете hold = [i-temp1] * L.Вместо этого, почему бы не сделать hold = i-temp, а затем выполнить цикл for j in xrange(L): и temp2 += hold?Это будет использовать меньше памяти, но в остальном будет вести себя точно так же.

...