C# Linear Interpolation

Question

I'm having issues interpolating a data file, which i have converted from .csv into an X array and Y array where X[0] corresponds to point Y[0] for example.

I need to interpolate between the values to give me a smooth file at the end. I am using a Picoscope to output the function which means each line is equally spaced in time, so only using Y values, which is why I'm trying to do this in a strange way when you see my code.

The kind of values it has to deal with are:

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

So where 2 Y values next to each other are the same its a straight line between them but when they vary like from X = 10 on wards it becomes a sine wave in this example.

I have been looking at the equations for interpolation etc and other posts on here but i haven't done that sort of maths for years and sadly i can't figure it out any more, because there's 2 unknowns and i can't think how to program that into c#.

What i have is:

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

and this creates a close waveform but it is still very steppy.

So can anyone see why this isn't very accurate? How to improve its accuracy? Or a different way of doing this using arrays?

Thanks for looking :)

Answer 1

Try this method for me:

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

Effectively you should be able to take your arrays and use it like this:

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

I pulled the code from here, but verified that the algorithm matched the theory here, and so I think it's right. However, you probably should consider using polynomial interpolation if this is still steppy, please note the theory link, it shows that linear interpolation produces steppy waves.

So, the first link I gave, where I grabbed this code from, also has a polynomial algorithm:

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

To use this one you're going to have to decide how you want to step up your x values, so let's say we wanted to do it by finding the midpoint between the current iteration and the next:

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

Please note that there will be more to this loop, like ensuring there is an i+1 and such, but I wanted to see if I could give you a start.

Answer 2

Theoretical base at Wolfram

The solution below computes the averages of Y values for given points with same X, just as the Matlab polyfit function does

Linq and .NET framework version > 3.5 are mandatory for this swift API. Comments inside the code.

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

Use it like:

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

Or:

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