Interpolare liniară. Cum se implementează acest algoritm în C? (este dată versiunea Python) (Programare, Python, C, Algoritm, Math, Prelucrarea Semnalelor)

psihodelia a intrebat.
a intrebat.

Există o metodă foarte bună de interpolare liniară. Aceasta realizează o interpolare liniară care necesită cel mult o multiplicare pentru fiecare eșantion de ieșire. Am găsit descrierea ei într-o a treia ediție a cărții Understanding DSP de Lyons. Această metodă implică un buffer de reținere special. Având în vedere un număr de eșantioane care trebuie inserate între două eșantioane de intrare oarecare, aceasta produce puncte de ieșire folosind interpolarea liniară. Aici, am rescris acest algoritm folosind 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)

unde x conține eșantioane de intrare, L este un număr de puncte care trebuie inserate, y va conține eșantioane de ieșire.

Întrebarea mea este cum să implementez un astfel de algoritm în ANSI C într-un mod cât mai eficient, de exemplu, este posibil să se evite cea de-a doua buclă?

NOTĂ: codul Python prezentat este doar pentru a înțelege cum funcționează acest algoritm.

UPDATE: iată un exemplu despre cum funcționează în 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)

Să spunem că x=[…. 10, 20, 30 ….]. Atunci, dacă L=1, se va produce [… 10, 15, 20, 20, 25, 30 …]

Comentarii

  • Dacă doriți doar să o implementați în C pentru performanță, dar totuși să o folosiți din Python, vă recomand Cython. –  > Por Björn Pollex.
  • codul Python prezentat este doar pentru a înțelege cum funcționează acest algoritm –  > Por psihodelia.
  • Este mai ușor să înțelegi cum funcționează un algoritm dacă folosești nume de variabile semnificative. –  > Por Lennart Regebro.
  • @Lennart: Ați actualizat. Acum ar trebui să fie ușor de înțeles –  > Por psihodelia.
  • Există o referință online la algoritm? De asemenea, evitarea înmulțirilor sună ca o optimizare din anii ’80… – user180326
5 răspunsuri
Nas Banov

Interpolare în sensul de „creștere a ratei de eșantionare a semnalului”

… sau eu îi spun, „upsampling” (termen greșit, probabil. disclaimer: nu am citit lucrarea lui Lyons). Trebuia doar să înțeleg ce face codul și apoi să-l rescriu pentru lizibilitate. Așa cum a fost dat, are câteva probleme:

a) este ineficient – două bucle este în regulă, dar face multiplicare pentru fiecare element de ieșire; de asemenea, folosește liste intermediare(hold), generează un rezultat cu append (bere mică)

b) interpolează greșit primul interval; generează date false în fața primului element. Să zicem că avem multiplicator=5 și seq=[20,30] – va genera [0,4,8,8,12,16,20,22,24,28,30] în loc de [20,22,24,26,28,30].

Iată, așadar, algoritmul sub forma unui 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 și acum câteva teste:

>>> 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]

Și iată traducerea mea în C, adaptată la șablonul fn al lui 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;
            }
        }
    }
}

Vă rugăm să rețineți că fragmentul C este „tastat, dar niciodată compilat sau rulat”, așa că pot exista erori de sintaxă, erori off-by-1 etc. Dar, în general, ideea este acolo.

Comentarii

  • oh bummer, recompensa a expirat și doar 1/2 a fost acordată. celelalte 25pts – phewt! – au dispărut de pe tărâmul stackoverflow. –  > Por Nas Banov.
Ignas Mikalajūnas

În acest caz, cred că poți evita a doua buclă:

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)

pur și simplu calculezi direct membrul din poziția dorită. Deși – s-ar putea să nu fie cea mai eficientă modalitate de a face acest lucru. Singura modalitate de a fi sigur este să compilați și să vedeți care este mai rapidă.

Comentarii

  • Nu, ați copiat o propoziție greșită. Eu am vorbit despre ce face bufferul de reținere. Algoritmul complet produce o interpolare liniară, nu doar elemente care se repetă. Vă rog să vedeți întrebarea mea actualizată. –  > Por psihodelia.
  • Da, numărul de calcule nu este mai mic în acest fel. Bucla în sine nu necesită timp semnificativ, deci nu este nevoie să o evităm. Totuși, ar trebui să funcționeze. +1 pentru că ai arătat soluții alternative (chiar dacă inutile). 😉 –  > Por Lennart Regebro.
Lennart Regebro

Ei bine, în primul rând, codul tău este stricat. L nu este definit, și nici y sau x.

După ce acest lucru este rezolvat, rulați cython pe codul rezultat:

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)

Și se pare că a funcționat. Nu am încercat să îl compilez, totuși, și puteți, de asemenea, să îmbunătățiți mult viteza prin adăugarea diferitelor optimizări.

„De exemplu, este posibil să se evite a doua buclă?”

Dacă da, atunci este posibil și în Python. Și nu văd cum, deși nu înțeleg de ce ai face-o în felul în care o faci tu. În primul rând crearea unei liste cu lungimea L a lui i-temp este complet inutilă. Pur și simplu faceți o buclă de L ori:

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)

Totuși, totul pare prea complicat pentru ceea ce se obține. Ce încerci să faci, de fapt? Să interpolați ceva? (Duh, așa scrie în titlu. Îmi pare rău pentru asta).

Există cu siguranță modalități mai simple de interpolare.

Update, o funcție de interpolare mult simplificată:

# 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)

Nu văd o modalitate de a scăpa de bucla interioară, nici un motiv pentru care să vrei să faci asta. conversia în C o las pe seama altcuiva, sau chiar mai bine, Cython, deoarece C este un limbaj excelent dacă vrei să vorbești cu hardware, dar altfel este inutil de dificil.

Comentarii

  • @psihodelia: Dar i-temp1 nu face o listă, deci nu va fi [10, 20, 10]. [i-temp1] va fi întotdeauna o listă cu o singură valoare. –  > Por Lennart Regebro.
  • @Lennart, * L va face ca lista să se repete. [1]*2 este [1,1]. Dar da, am observat același lucru, lista repetată nu prea ajută la exprimarea algoritmului. –  > Por mtrw.
  • @Lennart – Am înțeles greșit comentariul tău anterior, nu mi-am dat seama că noi doi subliniam același lucru. Îmi cer scuze. –  > Por mtrw.
  • Noua ta interpolare simplificată este mult mai lentă! Aveți o operație de împărțire într-o buclă. Scopul algoritmului prezentat este de a avea doar o singură înmulțire pe eșantion. Numărul de operații pe eșantion este foarte important în prelucrarea semnalelor. –  > Por psihodelia.
  • @psihodelia – este o împărțire cu o constantă. Compilatorul ar trebui să vadă asta și să o transforme într-o multiplicare. Dacă profilarea sau verificarea ansamblului arată că nu o face, atunci doar declară un temp și precalculează-l. Oricum, acest algoritm oferă rezultatul corect, ceea ce pare mai important. –  > Por mtrw.
mtrw

Cred că aveți nevoie de cele două bucle. Trebuie să treci peste eșantioanele din x pentru a inițializa interpolatorul, ca să nu mai vorbim de copierea valorilor lor în y, și trebuie să treci peste eșantioanele de ieșire pentru a le completa valorile. Presupun că ai putea face o buclă pentru a copia x în locurile corespunzătoare din y, urmată de o altă buclă pentru a utiliza toate valorile din y, dar acest lucru va necesita totuși o logică de trecere. Este mai bine să folosiți abordarea cu bucle imbricate.

(Și, după cum subliniază Lennart Regebro) Ca o notă secundară, nu înțeleg de ce faceți hold = [i-temp1] * L. În schimb, de ce să nu faceți hold = i-temp, și apoi bucla for j in xrange(L): și temp2 += hold? Acest lucru va utiliza mai puțină memorie, dar în rest se va comporta exact la fel.

Comentarii

  • Vă rugăm să nu acordați o mare atenție codului Python. Eu doar l-am scris pe baza unei diagrame vizuale a sistemului LTI. Este doar pentru a înțelege cum funcționează. Da, nu este nevoie să folosiți un buffer în cod, dar bufferii sunt ușor de plasat pe diagramă. –  > Por psihodelia.
Kratz

Iată încercarea mea de a face o implementare în C pentru algoritmul tău. Înainte de a încerca să îl optimizați în continuare, vă sugerez să faceți un profil al performanței sale cu toate optimizările compilatorului activate.

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