2012-10-11 97 views
7

我遇到了内插数据文件的问题,例如,我已经将.csv转换为X数组和Y数组,其中X [0]对应于点Y [0]。C#线性插值

我需要在这些值之间进行插值以在最后给我一个平滑的文件。 我正在使用一个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++; 
} 

,这将创建一个紧密的波形,但它仍然是非常steppy。

所以任何人都可以看到为什么这不是很准确? 如何提高其准确性? 或使用数组做这个不同的方式?

为寻找谢谢:)

+0

可能重复http://stackoverflow.com/questions/ 350852 /最小二乘法c-sharp-library) –

回答

11

试试这个方法对我来说:

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

我拉来自here的代码,但验证了该算法与理论here匹配,所以我t hink是的。但是,如果这仍然很简单,那么可能应该考虑使用多项式插值,请注意理论上的联系,它表明线性插值会产生粗糙波。

所以,我给第一个链接,在那里我抓住这个代码,还有一个多项式算法:

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和这样的,但我想看看我是否应该给你一个开始。

+0

那么我把什么值放入x? ,因为它的一个未知数是不是? – user1548411

+0

@ user1548411,请参阅我的编辑。 –

+0

如果y0,y1很大,则返回(y0 + y1)/ 2可能会引发溢出异常。您可以使用y0 +(y1 - y0)/ 2来修复此购买。 – openshac

1

Theoretical base at Wolfram

下面的溶液计算Y值的平均值用于与相同的X给定的点,正如Matlab的polyfit函数做

LINQ和.NET框架版本> 3.5是强制性的这个快速的API。 评论里面的代码。

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 }); 
的[最小二乘C#库(
+0

你写这个吗?如果不是,它是从哪里来的,代码的授权是什么? – ozziwald

+0

是的。您可以免费使用它。这是一个基于wolfram链接中引用的最小平方的朴素实现 –