text1

最小二乘法拟合用在嵌入式系统的代码,单片机上可以测试过

/*拟合直线方程为y = ax +b*/

#define N_MAX  20    //拟合折线的总数据个数
#define N_TOTAL 15   //实际需拟合折线的总数据个数
#define N_SECTION 3  //一段拟合的数据个数
#define SECTION (unsigend char) (N_TOTAL/N_SECTION)//段数
unsigned int x[N_MAX],y[N_MAX]; //拟合折线的数据
unsigned int dx[N_SECTION],dy[N_SECTION]; //一段拟合的数据
long k[2][3],m1,m2,m0;
int i = 0,j = 0;
//求∑Ai*Bi
long fsum(unsigned int *a,unsigned int *b,int c)
{
        long sum=0;
        for(i=0;i < c;i++)
                sum+=a[i]*b[i];
        return sum; 
}

//求矩阵
long fmatrix(unsigned char m,unsiged char n)
{
        long matrix;
        matrix=k[0][m]*k[1][n]-k[0][n]*k[1][m];
        return matrix;
}
//求折线参数
unsigned char get_value(float *pa,float *pb)
{
	long n_section_temp;
        double mi[N_SECTION];                //大小为1的数列，矩阵求和时匹配使
        for(i=0;i < N_SECTION;i++)
                mi[i]=1;
	n_section_temp = (long) N_SECTION;
//求线性方程系数
        k[0][0]=fsum(dx,dx,n_section_temp);
        k[0][1]=fsum(dx,mi,n_section_temp);
        k[0][2]=-fsum(dx,dy,n_section_temp);
        k[1][0]=fsum(dx,mi,n_section_temp);
        k[1][1]=n_section_temp;
        k[1][2]=-fsum(mi,dy,n_section_temp);
//输出线性方程系数
        m0=fmatrix(0,1);
        m1=fmatrix(1,2);
        m2=fmatrix(2,0);

        if(m0==0) return 0;
	else
	{	 
	*pa = (float) (m2/m0);
	*pb = (float) (m1/m0);
	return 1;
	}	
}
//求最大值
float get_max(float * pp,int pp_size)
{
	float max_temp = *pp;
	for(i = 1;i<pp_size;i++)
	{
	if(max_temp<*(pp+i))
	{
	max_temp = *(pp+i);
	}
	}
	return max_temp;
}
main()
{
	float aBuffer[SECTION],bBuffer[SECTION];
	for(i = 0;i<SECTION;i++)
	{
		for(j = 0;j<N_SECTION;j++)
		{
		dx[j] = x[i*N_SECTION+j];
		dy[j] = y[i*N_SECTION+j];
		}
		get_value(aBuffer,bBuffer);
		aBuffer ++;
		bBuffer ++;
	}
}