作业需要,需要在单片机上显示距离与ad值大小的关系,理论推倒得到这两个的关系是抛物线,因此查书写了这段c语言程序。
具体的证明过程我也不会,书上提到了对于抛物线采用多项式拟合的方法,好像是将x,x^2,看成了两个变量。
但核心公式只有一个,下面这个正则方程组
把上面这个方程简写成 xa=y(x、a、y分别对应一个矩阵)可以得到a=x^(-1)y
其中x^(-1)表示x的逆矩阵
具体的算法步骤是
1.求x的伴随矩阵
2.对矩阵x进行行列式操作
3.x的逆矩阵=x的伴随矩阵/x的行列式的值
4计算x^(-1)y
在下面的代码中对求伴随矩阵进行了简化,减少了一点运算量。
#include "stdio.h"
#ifndef uint8
#define uint8 unsigned short int
#endif
char Least_Squarel_Linear_Fit_y_2a0_4a1x_4a2xx(float *x,float *y,uint8 n,float *a0,float *a1,float *a2)
{
float temp=0;
float x0=n,x1=0,x2=0,x3=0,x4=0,x0y1=0,x1y1=0,x2y1=0;//表示各项的求和
float x0x2,x0x3,x0x4,x1x1,x1x2,x1x3,x1x4,x2x2,x2x3,x2x4,x3x3;
float a_banshui[3][3],a_ni[3][3];//这个是伴随矩阵和逆矩阵
float a_abs;//矩阵a的行列式
for(uint8 i=0;i<n;i++)
{
x0y1+=*(y+i);
temp=*(x+i);
x1+=temp;
x1y1+=temp*(*(y+i));
temp*=*(x+i);
x2+=temp;
x2y1+=temp*(*(y+i));
temp*=*(x+i);
x3+=temp;
temp*=*(x+i);
x4+=temp;
}
//后面要用到的数据,先算出来
x0x2=x0*x2;
x0x3=x0*x3;
x0x4=x0*x4;
x1x1=x1*x1;
x1x2=x1*x2;
x1x3=x1*x3;
x1x4=x1*x4;
x2x2=x2*x2;
x2x3=x2*x3;
x2x4=x2*x4;
x3x3=x3*x3;
//计算伴随矩阵,行列式,求逆矩阵,其实可以利用对称性再减少运算
a_banshui[0][0]= (x2x4-x3x3); a_banshui[0][1]=-(x1x4-x2x3); a_banshui[0][2]= (x1x3-x2x2);
a_banshui[1][0]=-(x1x4-x2x3); a_banshui[1][1]= (x0x4-x2x2); a_banshui[1][2]=-(x0x3-x1x2);
a_banshui[2][0]= (x1x3-x2x2); a_banshui[2][1]=-(x0x3-x1x2); a_banshui[2][2]= (x0x2-x1x1);
//计算矩阵对应行列式的值
a_abs=(x0*a_banshui[0][0]+x1*a_banshui[0][1]+x2*a_banshui[0][2]);
//计算逆矩阵
for(uint8 i=0;i<3;i++)
{
for(uint8 j=0;j<3;j++)
{
a_ni[i][j]=a_banshui[i][j]/a_abs;
}
}
*a0=a_ni[0][0]*x0y1+a_ni[0][1]*x1y1+a_ni[0][2]*x2y1;
*a1=a_ni[1][0]*x0y1+a_ni[1][1]*x1y1+a_ni[1][2]*x2y1;
*a2=a_ni[2][0]*x0y1+a_ni[2][1]*x1y1+a_ni[2][2]*x2y1;
return 1;
}
#define N 10
float x[N]={1,2,3,4,6,7,8,9,11,12};
float y[N]={2,3,6,7,5,3,2,1,1,1};
float a0,a1,a2;
int main()
{
Least_Squarel_Linear_Fit_y_2a0_4a1x_4a2xx(x,y,N,&a0,&a1,&a2);
printf("a0=%f\na1=%f\na2=%f\n",a0,a1,a2);
while(1);
}
结果与excel比较