在实际应用中，我们可能遇到这样一些数据：大小随温度变化但不是直线变化，无法用简单的一次线性函数描述，所以经常要对某点数据进行近似求取，例如MEMS陀螺仪输出与温度的关系，通常生产商已经为大家提供了温度输出，为了进行较为准确的补偿，就需要进行曲线拟合。最小二乘拟合是一种数学上的近似和优化，利用已知的数据得出一条直线或者曲线，使之在坐标系上与已知数据之间的距离的平方和最小。

下面的程序是实现对10个抽样点进行拟合，生成一条用y=a0+a1*x+a2*x^2+a3*x^3+a4*x^4描述的曲线，最终求出factor,即为系数a的矩阵，大小为5行1列.factor[0][0]为常数项a0，factor[4][0]为4次项系数a4.

#include   /*注意头文件*/
#include
#include       /*注意头文件math必须有*/

#define RM 5       /*矩阵行数为5，宏常量*/
#define CM 5       /*矩阵列数为5，宏常量*/
   
 
 float mat[RM][CM]={0};   /*矩阵ATA大小为5行5列*/
 float imat[RM][CM]={0};  /*矩阵ATA的转置矩阵大小为5行5列*/
 float ATb[5][1]={0};      /*矩阵ATb大小为5行1列*/
 float factor[5][1]={0};  /*最终结果系数a的矩阵大小为5行1列*/

 /*软件动态给定初始参数，10个样本点，x0---x9,y0----y9，共20个浮点型数据*/
 float x0,x1,x2,x3,x4,x5,x6,x7,x8,x9,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9;

 float X[10]={0};  /*由初始给定数据x组成的矩阵*/
 float Y[10]={0};  /*由初始给定数据y组成的矩阵*/
 
 /*矩阵ATA的25个元素定义为a1---a5,b1---b5,c1---c5,d1----d5,e1---e5，其中a1为固定值10*/
 float a1=10;
 float a2=0,a3=0,a4=0,a5=0,b1=0,b2=0,b3=0,b4=0,b5=0,c1=0,c2=0,c3=0,c4=0,c5=0,d1=0,d2=0,d3=0,d4=0,d5=0,e1=0,e2=0,e3=0,e4=0,e5=0;

 /*矩阵ATb的5个元素定义为u,v,w,p,q*/
 float u=0,v=0,w=0,p=0,q=0;

/*********************************************************************************
*函数名称：Inv()
*函数功能：求解矩阵ATA的逆矩阵
*调用函数：无
*入口参数：无
*出口参数：无
**********************************************************************************/ 
 
 void Inv()
{
    int i=0,j=0,k=0;

    for(i=0;i<5;i++)
            imat[i][i]=1;
    for(i=0;i<5;i++)
    {
        for(j=0;j<5;j++)
        {
            if(i!=j)
            {
                float t=mat[j][i]/mat[i][i];
                for(k=0;k<5;k++)
                {
                    mat[j][k]-=mat[i][k]*t;
                    imat[j][k]-=imat[i][k]*t;
                }
            }
        }
    }

    for(i=0;i<5;i++)
        if(mat[i][i]!=1)
        {
            float t=mat[i][i];
            for(j=0;j<5;j++)
            {
                mat[i][j]=mat[i][j]/t;
                imat[i][j]=imat[i][j]/t;
            }
        }



    
}

/*********************************************************************************
*函数名称：main()
*函数功能：主函数
*调用函数：Inv()
*入口参数：无
*出口参数：无
**********************************************************************************/ 
void main()
{
 int m=0,n=0,i=0,j=0;

/*给定20个数据做例子测试，实际中可以动态给定，注意：实际应用时请自行赋值！！！！！！*/
 x0=1;x1=2;x2=3;x3=4;x4=5;x5=6;x6=7;x7=8;x8=9;x9=10;y0=6;y1=2;y2=5;y3=6;y4=7;y5=8;y6=3;y7=2;y8=9;y9=4;

 X[0]=x0;
 X[1]=x1;
 X[2]=x2;
 X[3]=x3;
 X[4]=x4;
 X[5]=x5;
 X[6]=x6;
 X[7]=x7;
 X[8]=x8;
 X[9]=x9;

 Y[0]=y0;
 Y[1]=y1;
 Y[2]=y2;
 Y[3]=y3;
 Y[4]=y4;
 Y[5]=y5;
 Y[6]=y6;
 Y[7]=y7;
 Y[8]=y8;
 Y[9]=y9;

/*根据x0---x9求矩阵ATA*/
 a1=10;

for(i=0;i<10;i++)
{
 a2 = a2 + X[i];
}

for(i=0;i<10;i++)
{
 a3 = a3 + pow(X[i],2);
}

for(i=0;i<10;i++)
{
 a4 = a4 + pow(X[i],3);
}
for(i=0;i<10;i++)
{
 a5 = a5 + pow(X[i],4);
}
b1=a2;
b2=a3;
b3=a4;
b4=a5;


for(i=0;i<10;i++)
{
 b5 = b5 + pow(X[i],5);
}

c1=b2;
c2=b3;
c3=b4;
c4=b5;

for(i=0;i<10;i++)
{
 c5 = c5 + pow(X[i],6);
}

d1=c2;
d2=c3;
d3=c4;
d4=c5;

for(i=0;i<10;i++)
{
 d5 = d5 + pow(X[i],7);
}

e1=d2;
e2=d3;
e3=d4;
e4=d5;

for(i=0;i<10;i++)
{
 e5 = e5 + pow(X[i],8);
}

 mat[0][0]=a1;
 mat[0][1]=a2;
 mat[0][2]=a3;
 mat[0][3]=a4;
 mat[0][4]=a5;

 mat[1][0]=b1;
 mat[1][1]=b2;
 mat[1][2]=b3;
 mat[1][3]=b4;
 mat[1][4]=b5;

 mat[2][0]=c1;
 mat[2][1]=c2;
 mat[2][2]=c3;
 mat[2][3]=c4;
 mat[2][4]=c5;

 mat[3][0]=d1;
 mat[3][1]=d2;
 mat[3][2]=d3;
 mat[3][3]=d4;
 mat[3][4]=d5;

 mat[4][0]=e1;
 mat[4][1]=e2;
 mat[4][2]=e3;
 mat[4][3]=e4;
 mat[4][4]=e5;


     Inv();    /*求矩阵ATA的转置矩阵imat*/


/*根据x0---x9,y0----y9求矩阵ATb*/
 for(j=0;j<10;j++)
{
 u=u+Y[j];
}


for(j=0;j<10;j++)
{
 v = v + X[j]*Y[j];
}

for(j=0;j<10;j++)
{
 w = w + pow(X[j],2)*Y[j];
} 

for(j=0;j<10;j++)
{
 p = p + pow(X[j],3)*Y[j];
} 
for(j=0;j<10;j++)
{
 q = q + pow(X[j],4)*Y[j];
} 

 ATb[0][0]=u;
 ATb[1][0]=v;
 ATb[2][0]=w;
 ATb[3][0]=p;
 ATb[4][0]=q;

   /*根据ATA,ATb求矩阵factor,即为系数a的矩阵，大小为5行1列.factor[0][0]为常数项，factor[4][0]为4次项系数，其余类推*/
   factor[0][0] = imat[0][0]*ATb[0][0]+imat[0][1]*ATb[1][0]+imat[0][2]*ATb[2][0]+imat[0][3]*ATb[3][0]+imat[0][4]*ATb[4][0];
   factor[1][0] = imat[1][0]*ATb[0][0]+imat[1][1]*ATb[1][0]+imat[1][2]*ATb[2][0]+imat[1][3]*ATb[3][0]+imat[1][4]*ATb[4][0];
   factor[2][0] = imat[2][0]*ATb[0][0]+imat[2][1]*ATb[1][0]+imat[2][2]*ATb[2][0]+imat[2][3]*ATb[3][0]+imat[2][4]*ATb[4][0];
   factor[3][0] = imat[3][0]*ATb[0][0]+imat[3][1]*ATb[1][0]+imat[3][2]*ATb[2][0]+imat[3][3]*ATb[3][0]+imat[3][4]*ATb[4][0];
   factor[4][0] = imat[4][0]*ATb[0][0]+imat[4][1]*ATb[1][0]+imat[4][2]*ATb[2][0]+imat[4][3]*ATb[3][0]+imat[4][4]*ATb[4][0];


}