c# - C# 線形補間

Question

たとえば、.csv から X[0] がポイント Y[0] に対応する X 配列と Y 配列に変換したデータ ファイルの補間に問題があります。

最後に滑らかなファイルを得るために、値の間を補間する必要があります。私は Picoscope を使用して関数を出力しています。これは、各行が時間的に等間隔であることを意味するため、Y 値のみを使用しているため、コードを見るときに奇妙な方法でこれを実行しようとしています。

処理する値の種類は次のとおりです。

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

したがって、隣り合う 2 つの Y 値は同じで、それらの間の直線ですが、X = 10 のように変化すると、この例では正弦波になります。

私は補間などの方程式やここの他の投稿を見てきましたが、そのような数学を何年も行っていませんでした.2つの未知数があり、どのように考えることができないので、残念ながらもう理解できません.それをC＃にプログラムします。

私が持っているものは次のとおりです。

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

これにより、近い波形が作成されますが、それでも非常に階段状です。

では、なぜこれがあまり正確ではないのか、誰にもわかりますか? その精度を向上させる方法は？または、配列を使用してこれを行う別の方法はありますか?

見てくれてありがとう:)

Answer 1

私のためにこの方法を試してください：

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

事実上、配列を取得して次のように使用できるはずです。

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

ここからコードを引っ張ってきましたが、アルゴリズムがここの理論と一致することを確認したので、正しいと思います。ただし、これがまだ階段状である場合は、おそらく多項式補間の使用を検討する必要があります。理論のリンクに注意してください。線形補間が階段状の波を生成することが示されています。

したがって、このコードを取得した最初のリンクには、多項式アルゴリズムもあります。

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

これを使用するには、値をステップアップする方法を決定する必要がxあるため、現在の反復と次の反復の間の中間点を見つけることでそれを実行したいとしましょう:

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

などがあることを確認するなど、このループにはさらに多くのことがあることに注意してくださいi+1。

Answer 2

Wolframの理論的基盤

以下のソリューションは、Matlab の polyfit 関数と同じように、同じ X を持つ特定のポイントの Y 値の平均を計算します。

この迅速な API には、Linq および .NET Framework バージョン > 3.5 が必須です。 コード内のコメント。

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

次のように使用します。

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

または：

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