拟合算法

3.1 直线拟合

原理：已知一组数据Xi,Yi。Xi与Yi具有线性关系；利用最小二乘法，有如下关系

Φ(a,b)= (aX_i+b-Y_i)²    （1）

要使得Φ(a,b)取得最小值，那么

Φ(a,b)/ a=0, Φ(a,b)/ b=0

是极小值点的必要条件，于是(1)可化为

2X_i(aX_i+b-Y_i)=0

2_i(aX_i+b-Y_i)=0

对上式求解可得：

a=[(n+1) X_iY_i- X_i Y_i]/[(n+1) X_i²-( X_i)²]

b=[ X_i² Yi- X_i XiY_i]/[(n+1) X_i²-( X_i)²]

typedef struct tagLpoint

{

       double x;

       double y;

}Lpoint;

 

int LineMin2(HDC hDC,Lpoint *lp,int n)

{

       int           i;

       double     a,b,sumX,sumX2,sumY,sumXY;

 

       sumX=0.0;

       sumX2=0.0;

       sumY=0.0;

       sumXY=0.0;

 

       for(i=0;i<n;i++)

       {

              sumX+=lp[i].x;

              sumX2+=lp[i].x*lp[i].x;

              sumY+=lp[i].y;

              sumXY+=lp[i].x*lp[i].y;

       }

 

       a=(n*sumXY-sumX*sumY)/(n*sumX2-sumX*sumX);

       b=(sumX2*sumY-sumX*sumXY)/(n*sumX2-sumX*sumX);

 

       MoveToEx(hDC,0,-b/a,NULL);

       LineTo(hDC,600,(600-b)/a);  //600表示从高度600开始画

       SelectObject(hDC,CreatePen(0,1,RGB(255,0,0)));

       MoveToEx(hDC,0,b,NULL);

       LineTo(hDC,700,700*a+b);//0-700的x值

       return 1;

}

3.2 二次以上拟合

typedef struct tagLpoint

{

       double     x;

       double     y;

}Lpoint;

//检查方程组是否有解

int Check(double s[10][10],int R,int k)

{

       int           i;

 

       i=k;

       do i=i+1;

       while((i<=R)&&(abs(s[i][k])<=ER));

       if(i>=R)

       {

              MessageBox(NULL,"方程组无解或无惟一解，程序终止",NULL,NULL);

              return 0;

       }

       return 1;

}

///把最大元素作为对角线元数

int Pivot(double s[10][10],double c[10],int R,int k)

{

       int           l,j;

       double     swap;

 

       l=k;

       for(j=k;j<R;j++)

              if(abs(s[j][k])>abs(s[l][k]))

                     l=j;//即找出最大绝对值的下标

       if(l==k)

              return 0;

       swap=c[k];

       c[k]=c[l];

       c[l]=swap;

       for(j=k;j<R;j++)

       {

              swap=s[k][j];

              s[k][j]=s[l][j];

              s[l][j]=swap;

       }

       return 1;

}

 

//高斯消元法解线性方程组

int GussianSalve(double s[10][10],double c[10],int R,double a[10])

{

       int           i,j,k,ret;

       double     f;

 

       for(k=0;k<R-1;k++)

       {

              ret=Check(s,R,k);

              if(ret==0)

                     return 0;

              Pivot(s,c,R,k);

              for(i=k+1;i<R;i++)

              {

                     f=s[i][k]/s[k][k];

                     c[i]=c[i]-c[k]*f;

                     for(j=k+1;j<R;j++)

                            s[i][j]=s[i][j]-s[k][j]*f;

              }

       }

       a[R-1]=c[R-1]/s[R-1][R-1];

       for(i=R-2;i>=0;i--)

       {

              f=0.0;

              for(j=i+1;j<R;j++)

                     f=f+s[i][j]*a[j];

              a[i]=(c[i]-f)/s[i][i];

       }

       return 1;

}

 

//多项式最小二乘法算法

/*名称：PolyMin

 功能：对一组数据进行多项式拟合

 参数：hDC绘画句柄，tempasc时个包含x,y值的数组，R表示几次如2表示一次，意思是要求ab两个参数，n表示数组个数

*/

void PolyMin(HDC hDC,Lpoint *tempasc,int R,int n)

{

       int           i,j,k;

       double     s[10][10],c[10],a[10];

       double     **px,temp,te;

 

       px=new double*[n];

       for(i=0;i<n;i++)

              px[i]=new double[10];

 

       for(i=0;i<R;i++)

       {

              c[i]=0.0;

              for(j=0;j<R;j++)

                     s[i][j]=0.0;

       }

       for(i=0;i<n;i++)

       {

              px[i][0]=1.0;

              for(j=1;j<R;j++)

                     px[i][j]=px[i][j-1]*tempasc[i].x;

       }

       for(i=0;i<R;i++)

       {

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

                     c[i]=c[i]+tempasc[j].y*px[j][i];

              for(j=0;j<R;j++)

                     for(k=0;k<n;k++)

                            s[i][j]=s[i][j]+px[k][i]*px[k][j];

       }

       GussianSalve(s,c,R,a);

       SelectObject(hDC,CreatePen(0,1,RGB(255,0,0)));

       MoveToEx(hDC,40,400-50*a[0],NULL);

       for(te=0;te<31;te+=0.25)

       {

              temp=a[R-1];

              for(j=R-1;j>0;j--)

                     temp=temp*te+a[j-1];

              LineTo(hDC,40+te*50,400-temp*50);//从坐标点40,500处开始画

       }

       for(i=0;i<n;i++)

              delete []px[i];

       delete []px;

}

3.3 拉格朗日插值曲线

double     LagBasChg(Lpoint *p,int n,double ly,int k)

{

       int           i;

       double     lf;

 

       lf=1.0;

       for(i=0;i<n;i++)

              if(i!=k)

                     lf=lf*(ly-p[i].x)/(p[k].x-p[i].x);

       return lf;

}

void LagrangeInterpChg(HDC hDC,int color,Lpoint *p,int n)

{

       int           i;

       double     ly,te;

       HPEN     hPen;

 

       ly=0.0;

       for(i=0;i<n;i++)

              ly=ly+p[i].pasc*LagBasChg(p,n,0.0,i);

       hPen=CreatePen(0,1,color);

       SelectObject(hDC,hPen);

       MoveToEx(hDC,40,500-ly,NULL);

       for(te=0;te<=30;te+=0.25)

       {

              ly=0.0;

              for(i=0;i<N;i++)

                     ly=ly+p[i].y*LagBasChg(p,n,te,i);

              LineTo(hDC,40+te*20,500-ly);

       }

       DeleteObject(hPen);

}