天天看点

单片机c语言拟合二次曲线y=a0+a1x+a2x^2

作业需要,需要在单片机上显示距离与ad值大小的关系,理论推倒得到这两个的关系是抛物线,因此查书写了这段c语言程序。

具体的证明过程我也不会,书上提到了对于抛物线采用多项式拟合的方法,好像是将x,x^2,看成了两个变量。

但核心公式只有一个,下面这个正则方程组

单片机c语言拟合二次曲线y=a0+a1x+a2x^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比较

单片机c语言拟合二次曲线y=a0+a1x+a2x^2

继续阅读