C# Interpolación lineal

Question

Tengo problemas para interpolar un archivo de datos, que he convertido de .csv en una matriz X y matriz Y donde X [0] corresponde al punto Y [0], por ejemplo.

Necesito interpolar entre los valores para darme un archivo sin problemas al final. Estoy usando un Picoscope para generar la función, lo que significa que cada línea está equidistantemente espaciada en el tiempo, por lo que solo uso valores Y, por lo que trato de hacer esto de forma extraña cuando ve mi código.

el tipo de valores que tiene que hacer frente son:

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 

Así que 2 valores Y junto a la otra son los mismos que es una línea recta entre ellos, pero cuando se varían como desde x = 10 en wards se convierte en una onda sinusoidal en este ejemplo.

He estado buscando las ecuaciones para interpolación, etc. y otras publicaciones aquí, pero hace años que no hago ese tipo de cálculos y, por desgracia, ya no puedo resolverlo, porque hay 2 incógnitas y puedo No pienses cómo programar eso en C#.

Lo que tengo es:

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++; 
} 

y esto crea una forma de onda de cerca, pero todavía es muy steppy.

Entonces, ¿alguien puede ver por qué esto no es muy preciso? ¿Cómo mejorar su precisión? ¿O una forma diferente de hacer esto usando matrices?

Gracias por mirar :)

Answer 1

Pruebe este método para mí:

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); 
} 

eficaz, se deberían poder tomar sus matrices y utilizar de esta manera:

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

Saqué el código de here, pero verificó que el algoritmo coincidía con la teoría here, y entonces I t Hink es correcto. Sin embargo, es probable que deba considerar el uso de interpolación polinomial si esto todavía es pequeño, tenga en cuenta el enlace de la teoría, muestra que la interpolación lineal produce ondas estelares.

Por lo tanto, el primer enlace que he dado, cuando me agarró el código de, también tiene un algoritmo polinomio:

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; 
} 

Para utilizar éste va a tener que decidir cómo desea intensificar sus x valores, por lo que vamos a decir que queríamos hacerlo buscando el punto medio entre la iteración actual y la siguiente:

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

Tenga en cuenta que habrá más de este bucle, al igual que asegurar que haya un i+1 y tal , pero quería ver si yo c debería comenzar.

Answer 2

Theoretical base at Wolfram

La solución a continuación calcula los promedios de los valores de Y para puntos dados con la misma X, al igual que la función polyfit Matlab hace

Linq y .NET Framework versión> 3.5 son obligatorios para esta rápida API. Comentarios dentro del código.

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; 
    } 
} 

usarlo como:

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} 
    }); 

O:

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 });