天天看點

單片機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

繼續閱讀