C# Interpolacja liniowa

Question

Mam problemy z interpolacją pliku danych, który przekonwertowałem z .csv na tablicę X i tablicę Y, gdzie X [0] odpowiada na przykład Y [0].

Potrzebuję interpolacji między wartościami, aby uzyskać gładki plik na końcu. Używam Picoscope do wyprowadzania funkcji, co oznacza, że ​​każda linia jest równo rozłożona w czasie, więc tylko przy użyciu wartości Y, dlatego próbuję to zrobić w dziwny sposób, gdy widzisz mój kod.

Rodzaj wartości ma do czynienia to:

X  Y 
0  0 
2.5 0 
2.5 12000 
7.5 12000 
7.5 3000 
10 3000 
10 6000 
11 6625.254 
12 7095.154 

Więc gdzie 2 wartości Y obok siebie są takie same, jest to linia prosta pomiędzy nimi, ale gdy różnią się one niczym od X = 10 na Oddziały stają się falą sinusoidalną w tym przykładzie.

Patrzyłem na równania interpolacji itp. I inne posty tutaj, ale nie robiłem tego rodzaju matematyki od lat i niestety nie mogę tego rozgryźć, ponieważ są 2 niewiadome i mogę nie myśl, jak to zaprogramować na C#.

co mam to:

int a = 0, d = 0, q = 0; 
bool up = false; 
double times = 0, change = 0, points = 0, pointchange = 0; 
double[] newy = new double[8192]; 
while (a < sizeoffile-1 && d < 8192) 
{ 
    Console.Write("..."); 
    if (Y[a] == Y[a+1])//if the 2 values are the same add correct number of lines in to make it last the correct time 
    { 
     times = (X[a] - X[a + 1]);//number of repetitions 
     if ((X[a + 1] - X[a]) > times)//if that was a negative number try the other way around 
      times = (X[a + 1] - X[a]);//number of repetitions 
     do 
     { 
      newy[d] = Y[a];//add the values to the newy array which replaces y later on in the program 
      d++;//increment newy position 
      q++;//reduce number of reps in this loop 
     } 
     while (q < times + 1 && d < 8192); 
     q = 0;//reset reps 
    } 
    else if (Y[a] != Y[a + 1])//if the 2 values are not the same interpolate between them 
    { 
     change = (Y[a] - Y[a + 1]);//work out difference between the values 
     up = true;//the waveform is increasing 
     if ((Y[a + 1] - Y[a]) > change)//if that number was a negative try the other way around 
     { 
      change = (Y[a + 1] - Y[a]);//work out difference bwteen values 
      up = false;//the waveform is decreasing 
     } 
     points = (X[a] - X[a + 1]);//work out amount of time between given points 
     if ((X[a + 1] - X[a]) > points)//if that was a negative try other way around 
      points = (X[a + 1] - X[a]);///work out amount of time between given points 
     pointchange = (change/points);//calculate the amount per point in time the value changes by 
     if (points > 1)//any lower and the values cause errors 
     { 
      newy[d] = Y[a];//load first point 
      d++; 
      do 
      { 
       if (up == true)// 
        newy[d] = ((newy[d - 1]) + pointchange); 
       else if (up == false) 
        newy[d] = ((newy[d - 1]) - pointchange); 
       d++; 
       q++; 
      } 
      while (q < points + 1 && d < 8192); 
      q = 0; 
     } 
     else if (points != 0 && points > 0) 
     { 
      newy[d] = Y[a];//load first point 
      d++; 
     } 
    } 
    a++; 
} 

i tworzy ścisły fali, ale nadal jest bardzo steppy.

Ktoś może zobaczyć, dlaczego nie jest to bardzo dokładne? Jak poprawić dokładność? Albo inny sposób robienia tego za pomocą tablic?

Dzięki za szuka :)

Answer 1

Spróbuj tej metody dla mnie:

static public double linear(double x, double x0, double x1, double y0, double y1) 
{ 
    if ((x1 - x0) == 0) 
    { 
     return (y0 + y1)/2; 
    } 
    return y0 + (x - x0) * (y1 - y0)/(x1 - x0); 
} 

Skutecznie powinieneś być w stanie podjąć swoich tablice i używać go tak:

var newY = linear(X[0], X[0], X[1], Y[0], Y[1]); 

Wyciągnąłem kod z here, ale zweryfikowano, że algorytm pasuje do teorii here, a więc I t hink W porządku. Prawdopodobnie jednak powinieneś rozważyć zastosowanie interpolacji wielomianowej, jeśli nadal jest to steppy, proszę zwrócić uwagę na teorię, pokazuje ona, że ​​interpolacja liniowa wytwarza fale steppy.

Tak, pierwszy link dałem, gdzie złapałem ten kod, posiada również wielomian algorytmu:

static public double lagrange(double x, double[] xd, double[] yd) 
{ 
    if (xd.Length != yd.Length) 
    { 
     throw new ArgumentException("Arrays must be of equal length."); //$NON-NLS-1$ 
    } 
    double sum = 0; 
    for (int i = 0, n = xd.Length; i < n; i++) 
    { 
     if (x - xd[i] == 0) 
     { 
      return yd[i]; 
     } 
     double product = yd[i]; 
     for (int j = 0; j < n; j++) 
     { 
      if ((i == j) || (xd[i] - xd[j] == 0)) 
      { 
       continue; 
      } 
      product *= (x - xd[i])/(xd[i] - xd[j]); 
     } 
     sum += product; 
    } 
    return sum; 
} 

Aby korzystać z tego jeden pan będzie musiał zdecydować, jak chcesz, aby przyspieszyć twoje wartości x, więc powiedzmy, że chcemy to zrobić poprzez znalezienie środkowy między bieżącej iteracji i następny:

for (int i = 0; i < X.Length; i++) 
{ 
    var newY = lagrange(new double[] { X[i]d, X[i+1]d }.Average(), X, Y); 
} 

Należy pamiętać, że nie będzie więcej do tej pętli, jak zapewnienie istnieje i+1 i takie , ale chciałem zobaczyć, czy c Mógłbyś zacząć.

Answer 2

Theoretical base at Wolfram

Rozwiązanie poniżej oblicza średnie wartości Y dla poszczególnych punktów z tego samego X, podobnie jak funkcja polyfit Matlab robi

LINQ i .NET framework wersja> 3,5 są obowiązkowe dla ten szybki API. Komentarze wewnątrz kodu.

using System; 
using System.Collections.Generic; 
using System.Linq; 


/// <summary> 
/// Linear Interpolation using the least squares method 
/// <remarks>http://mathworld.wolfram.com/LeastSquaresFitting.html</remarks> 
/// </summary> 
public class LinearLeastSquaresInterpolation 
{ 
    /// <summary> 
    /// point list constructor 
    /// </summary> 
    /// <param name="points">points list</param> 
    public LinearLeastSquaresInterpolation(IEnumerable<Point> points) 
    { 
     Points = points; 
    } 
    /// <summary> 
    /// abscissae/ordinates constructor 
    /// </summary> 
    /// <param name="x">abscissae</param> 
    /// <param name="y">ordinates</param> 
    public LinearLeastSquaresInterpolation(IEnumerable<float> x, IEnumerable<float> y) 
    { 
     if (x.Empty() || y.Empty()) 
      throw new ArgumentNullException("null-x"); 
     if (y.Empty()) 
      throw new ArgumentNullException("null-y"); 
     if (x.Count() != y.Count()) 
      throw new ArgumentException("diff-count"); 

     Points = x.Zip(y, (unx, uny) => new Point { x = unx, y = uny }); 
    } 

    private IEnumerable<Point> Points; 
    /// <summary> 
    /// original points count 
    /// </summary> 
    public int Count { get { return Points.Count(); } } 

    /// <summary> 
    /// group points with equal x value, average group y value 
    /// </summary> 
                public IEnumerable<Point> UniquePoints 
{ 
    get 
    { 
     var grp = Points.GroupBy((p) => { return p.x; }); 
     foreach (IGrouping<float,Point> g in grp) 
     { 
      float currentX = g.Key; 
      float averageYforX = g.Select(p => p.y).Average(); 
      yield return new Point() { x = currentX, y = averageYforX }; 
     } 
    } 
} 
    /// <summary> 
    /// count of point set used for interpolation 
    /// </summary> 
    public int CountUnique { get { return UniquePoints.Count(); } } 
    /// <summary> 
    /// abscissae 
    /// </summary> 
    public IEnumerable<float> X { get { return UniquePoints.Select(p => p.x); } } 
    /// <summary> 
    /// ordinates 
    /// </summary> 
    public IEnumerable<float> Y { get { return UniquePoints.Select(p => p.y); } } 
    /// <summary> 
    /// x mean 
    /// </summary> 
    public float AverageX { get { return X.Average(); } } 
    /// <summary> 
    /// y mean 
    /// </summary> 
    public float AverageY { get { return Y.Average(); } } 

    /// <summary> 
    /// the computed slope, aka regression coefficient 
    /// </summary> 
    public float Slope { get { return ssxy/ssxx; } } 

    // dotvector(x,y)-n*avgx*avgy 
    float ssxy { get { return X.DotProduct(Y) - CountUnique * AverageX * AverageY; } } 
    //sum squares x - n * square avgx 
    float ssxx { get { return X.DotProduct(X) - CountUnique * AverageX * AverageX; } } 

    /// <summary> 
    /// computed intercept 
    /// </summary> 
    public float Intercept { get { return AverageY - Slope * AverageX; } } 

    public override string ToString() 
    { 
     return string.Format("slope:{0:F02} intercept:{1:F02}", Slope, Intercept); 
    } 
} 

/// <summary> 
/// any given point 
/// </summary> 
public class Point 
{ 
    public float x { get; set; } 
    public float y { get; set; } 
} 

/// <summary> 
/// Linq extensions 
/// </summary> 
public static class Extensions 
{ 
    /// <summary> 
    /// dot vector product 
    /// </summary> 
    /// <param name="a">input</param> 
    /// <param name="b">input</param> 
    /// <returns>dot product of 2 inputs</returns> 
    public static float DotProduct(this IEnumerable<float> a, IEnumerable<float> b) 
    { 
     return a.Zip(b, (d1, d2) => d1 * d2).Sum(); 
    } 
    /// <summary> 
    /// is empty enumerable 
    /// </summary> 
    /// <typeparam name="T"></typeparam> 
    /// <param name="a"></param> 
    /// <returns></returns> 
    public static bool Empty<T>(this IEnumerable<T> a) 
    { 
     return a == null || a.Count() == 0; 
    } 
} 

używać go jak:

var llsi = new LinearLeastSquaresInterpolation(new Point[] 
    { 
     new Point {x=1, y=1}, new Point {x=1, y=1.1f}, new Point {x=1, y=0.9f}, 
     new Point {x=2, y=2}, new Point {x=2, y=2.1f}, new Point {x=2, y=1.9f}, 
     new Point {x=3, y=3}, new Point {x=3, y=3.1f}, new Point {x=3, y=2.9f}, 
     new Point {x=10, y=10}, new Point {x=10, y=10.1f}, new Point {x=10, y=9.9f}, 
     new Point {x=100, y=100}, new Point{x=100, y=100.1f}, new Point {x=100, y=99.9f} 
    }); 

czyli

var llsi = new LinearLeastSquaresInterpolation(
    new float[] { 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9 }, 
    new float[] { 1, 2, 3, 4, 5, 6, 7, 8, 9, 1, 2, 3, 4, 5, 6, 7, 8, 9 });