作業需要,需要在單片機上顯示距離與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比較