爱魔域

数据拟合 曲线求二阶导数归零 ， 求拐点
为啥 不能上传文件呢
我找到了一个 拟合程序和 拐点程序 怎么将 两个连起来 呢 ？
public partial class Form1 : Form
{
public Form1()
{
InitializeComponent();
}
public struct point// public struct PersonStruct
{
public double x;
public double y;
};


// point[] p; //PersonStruct p1, p2; //与类一样，但可以不new int[] arr1;

point[] p = new point[200];
//正向直线方程
double fun(point p1, point p2, point p3)
{

double s;
s = (p2.x - p1.x) * (p3.y - p1.y) + (p1.y - p2.y) * (p3.x - p1.x);
return s;
}


//获取测试点拐点
void getPoint()
{
double i;
int j = 0;
for (i = 0; i < 7; i += 0.035)
{
p[j].x = i;
p[j].y = Math.Sin(i); //这里是sin 要是别的公式怎么办呢 ？？？？？？

////函数到底是什么？？？？？系数是多少 ，多少个系数？？？
j++;
}


}
//第一步： 最小二乘法曲线拟合 matlab polyfit
//第2步： 二阶导数归零，求离散拐点 matlab diff
private void button1_Click(object sender, EventArgs e)
{
//double[] y = { 110, 160, 190, 235, 264, 274, 295 };
//double[] x = { 10, 50, 130, 188, 255, 316, 487 };
double[] x = { 60, 120, 180, 240, 300, 360, 420 };
// double[] y = { 30,31, 32, 33, 33.5, 34, 34.5 };
double[] y = { -30, -31, -32, -33, -33.5, -34, -34.5 };
double[] f = ecf.MultiLine(x, y, 7, 4);//拟合函数，返回值是函数的系数，拟合代码在ecf.cs
//例如：y=a0+a1*x 返回值则为a0 a1
//例如：y=a0+a1*x+a2*x*x 返回值则为a0 a1 a2
//将ecf.cs拷到新项目的文件夹内，然后用添加现有项将其添加到项目中即可，注意要修改ecf.cs的namespace
//函数到底是什么？？？？？系数是多少 ，多少个系数？？？
Bitmap b1 = new Bitmap(500, 300);
Graphics G2 = Graphics.FromImage(b1);
SolidBrush mysbrush1 = new SolidBrush(Color.DarkOrchid);
Pen p1 = new Pen(mysbrush1);
p1.Color = Color.Green;
p1.Width = 3;

int i;
for (i = 0; i < y.Length - 1; i++)
G2.DrawLine(p1, new PointF((float)x[i], 300 - (float)y[i]),
new PointF((float)x[i + 1], 300 - (float)y[i + 1]));
p1.Color = Color.Red;
for (i = 0; i < (int)x.Max(); i++)
{
float y1 = (float)(f[0] + f[1] * i + f[2] * Math.Pow(i, 2) + f[3] * Math.Pow(i, 3) + f[4] * Math.Pow(i, 4));
float y2 = (float)(f[0] + f[1] * (i + 1) + f[2] * Math.Pow(i + 1, 2) + f[3] * Math.Pow(i + 1, 3) + f[4] * Math.Pow(i + 1, 4));
//G2.DrawLine(p1, new PointF((float)i, 300 - y1),
// new PointF((float)(i + 1), 300 - y2));
G2.DrawLine(p1, new PointF((float)i, (0 - 7*y1)),
new PointF((float)(i + 1), (0 - 7 * y2)));

}
G2.Save();
pictureBox1.Image = b1;
int j;


getPoint();
double s1, s2;
for (i = 3; i < 199; i++)
{
s1 = fun(p[i - 2], p[i - 1], p[i]);
s2 = fun(p[i - 1], p[i], p[i + 1]);


if (s1 * s2 < 0)
{
//printf(" %lf %lf\n", p[i].x, p[i].y);
label1.Text = p[i].x.ToString();
label2.Text = p[i].y.ToString();
}


}
}
}
public static double[] MultiLine(double[] arrX, double[] arrY, int length, int dimension)//二元多次线性方程拟合曲线
{
int n = dimension + 1; //dimension次方程需要求 dimension+1个 系数
double[,] Guass = new double[n, n + 1]; //高斯矩阵 例如：y=a0+a1*x+a2*x*x
for (int i = 0; i < n; i++)
{
int j;
for (j = 0; j < n; j++)
{
Guass[i, j] = SumArr(arrX, j + i, length);
}
Guass[i, j] = SumArr(arrX, i, arrY, 1, length);
}
return ComputGauss(Guass, n);

}

public static double SumArr(double[] arr, int n, int length) //求数组的元素的n次方的和
{
double s = 0;
for (int i = 0; i < length; i++)
{
if (arr[i] != 0 || n != 0)
s = s + Math.Pow(arr[i], n);
else
s = s + 1;
}
return s;
}
public static double SumArr(double[] arr1, int n1, double[] arr2, int n2, int length)
{
double s = 0;
for (int i = 0; i < length; i++)
{
if ((arr1[i] != 0 || n1 != 0) && (arr2[i] != 0 || n2 != 0))
s = s + Math.Pow(arr1[i], n1) * Math.Pow(arr2[i], n2);
else
s = s + 1;
}
return s;

}

public static double[] ComputGauss(double[,] Guass, int n)
{
int i, j;
int k, m;
double temp;
double max;
double s;
double[] x = new double[n];

for (i = 0; i < n; i++) x[i] = 0.0;//初始化


for (j = 0; j < n; j++)
{
max = 0;

k = j;
for (i = j; i < n; i++)
{
if (Math.Abs(Guass[i, j]) > max)
{
max = Guass[i, j];
k = i;
}
}



if (k != j)
{
for (m = j; m < n + 1; m++)
{
temp = Guass[j, m];
Guass[j, m] = Guass[k, m];
Guass[k, m] = temp;

}
}

if (0 == max)
{
// "此线性方程为奇异线性方程"

return x;
}


for (i = j + 1; i < n; i++)
{

s = Guass[i, j];
for (m = j; m < n + 1; m++)
{
Guass[i, m] = Guass[i, m] - Guass[j, m] * s / (Guass[j, j]);

}
}


}//结束for (j=0;j<n;j++)


for (i = n - 1; i >= 0; i--)
{
s = 0;
for (j = i + 1; j < n; j++)
{
s = s + Guass[i, j] * x[j];
}

x[i] = (Guass[i, n] - s) / Guass[i, i];

}

return x;
}//返回值是函数的系数
//例如：y=a0+a1*x 返回值则为a0 a1
//例如：y=a0+a1*x+a2*x*x 返回值则为a0 a1 a2