本帖最后由 cmz 于 2019-12-17 18:57 编辑
前言：因为在学校参加电赛选择的题目是仪器仪表类，经常要用到各种校准，之前的做法是采集几个样本点在Excel上计算拟合曲线，这样的方法在自己实验室里面数据还算准确的。但是到了电赛测评地点，因为仪器、环境的差异，当时迫切希望在测评现场进行校准。那样最后测量显示的数据将会非常准确，只要校准的阶数足够大，样本数量足够多。后来出来工作，在生产的时候做一致性校准也是将会用到。以下的算法得以成功，得益于在网上的学习，深知开源的强大。
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
最小二乘法应用 一、 模拟数据 1.模拟数据 模拟数据是用两个电阻分压测量出来的，主要用来验证算法的正确性，主要与Excel的拟合曲线做比较。
2.2.模拟数据在Excel的拟合公式 （1）.一元一次方程拟合结果
一元一次方程拟合结果
本算法计算结果： 得到拟合公式：y=7.585912x + 2.49931，（一阶还是非常准的） （2）.二元一次方程拟合结果
二元一次方程拟合结果
本算法计算结果： 得到拟合公式：y= 0.000003* x^2 + 7.584531x + 2.554939。（二阶计算结果的最后的常数项有点差别，未找到原因，，但是一点不影响拟合） 最小二乘法就是要确定这个一元一次方程和二元一次方程。当然三次，四次拟合也是可以的。 二、 最小二乘法实现 1.复习方差的计算方法 2.学习最小二乘法的几个网站 最小二乘法的中心思想：使总的误差的平方最小： 。 附录1：二元一次方程组通解： 附录2：三元一次方程组通解： 附录3：代码。 #include "stdio.h" #include "string.h" /* 程序版本 0001 */ #define OLS_ORDER_NUM 2 //阶数：1 y = a0 + a1*x //阶数：2 y = a0 + a1*x + a2*x^2 #define OLS_SAMPLE_NUM 7 //样本数量 #define OLS_POW2(x) (x*x) //计算平方 #define OLS_POW3(x) (x*x*x) //计算三次方 #define OLS_POW4(x) (x*x*x*x) //计算四次方 typedef struct { float leftBuf[OLS_ORDER_NUM + 1][OLS_ORDER_NUM + 1]; //左边的矩阵 float righeBuf[OLS_ORDER_NUM + 1]; //右边的矩阵 float resultBuf[OLS_ORDER_NUM + 1]; //结果矩阵，resultBuf[0]为a0，resultBuf[1]为a1 ... float sampleBuf[2][OLS_SAMPLE_NUM]; //样本Buf，sampleBuf[0]保存Xi，sampleBuf[1]保存Yi }OLS_t; OLS_t g_OLS; /******************************************************************************* * xxx@2019-12-13 * Function Name : OLS_Init * Description : 初始化 * Input : None * Output : None * Return : None *******************************************************************************/ void OLS_Init(void) { memset((unsigned char*)&g_OLS, 0, sizeof(OLS_t)); g_OLS.sampleBuf[0][0] = 0; //填充样本Xi g_OLS.sampleBuf[0][1] = 13; g_OLS.sampleBuf[0][2] = 145; g_OLS.sampleBuf[0][3] = 266; g_OLS.sampleBuf[0][4] = 357; g_OLS.sampleBuf[0][5] = 411; g_OLS.sampleBuf[0][6] = 490; g_OLS.sampleBuf[1][0] = 0; //填充样本Yi g_OLS.sampleBuf[1][1] = 103; g_OLS.sampleBuf[1][2] = 1106; g_OLS.sampleBuf[1][3] = 2014; g_OLS.sampleBuf[1][4] = 2716; g_OLS.sampleBuf[1][5] = 3118; g_OLS.sampleBuf[1][6] = 3720; } /******************************************************************************* * xxx@2019-12-13 * Function Name : OLS_CalcBufData * Description : 填充矩阵数据，根据最小二乘法拟合多项式的推导结果，来计算矩阵的数据 * Input : None * Output : None * Return : None *******************************************************************************/ void OLS_CalcBufData(void) { int i; #if (OLS_ORDER_NUM == 1) { g_OLS.leftBuf[0][0] = OLS_SAMPLE_NUM; //样本数量 for (i = 0; i < OLS_SAMPLE_NUM; i++) { g_OLS.leftBuf[0][1] += g_OLS.sampleBuf[0]; //Xi的累加和 g_OLS.leftBuf[1][1] += OLS_POW2(g_OLS.sampleBuf[0]); //Xi^2的累加和 g_OLS.righeBuf[0] += g_OLS.sampleBuf[1]; //Yi的累加和 g_OLS.righeBuf[1] += (g_OLS.sampleBuf[0] * g_OLS.sampleBuf[1]); //Xi * Yi的累加和 } g_OLS.leftBuf[1][0] = g_OLS.leftBuf[0][1]; //其中两个的 Xi的累加和 相等 } #elif (OLS_ORDER_NUM == 2) { g_OLS.leftBuf[0][0] = OLS_SAMPLE_NUM; //样本数量 for (i = 0; i < OLS_SAMPLE_NUM; i++) { g_OLS.leftBuf[0][1] += g_OLS.sampleBuf[0]; //Xi的累加和 g_OLS.leftBuf[0][2] += OLS_POW2(g_OLS.sampleBuf[0]); //Xi^2的累加和 g_OLS.leftBuf[1][2] += OLS_POW3(g_OLS.sampleBuf[0]); //Xi^3的累加和 g_OLS.leftBuf[2][2] += OLS_POW4(g_OLS.sampleBuf[0]); //Xi^4的累加和 g_OLS.righeBuf[0] += g_OLS.sampleBuf[1]; //Yi的累加和 g_OLS.righeBuf[1] += (g_OLS.sampleBuf[0] * g_OLS.sampleBuf[1]); //Xi * Yi的累加和 g_OLS.righeBuf[2] += (OLS_POW2(g_OLS.sampleBuf[0]) * g_OLS.sampleBuf[1]); //Xi^2 * Yi的累加和 } g_OLS.leftBuf[1][0] = g_OLS.leftBuf[0][1]; //其中两个的 Xi的累加和 相等 g_OLS.leftBuf[1][1] = g_OLS.leftBuf[0][2]; g_OLS.leftBuf[2][0] = g_OLS.leftBuf[0][2]; g_OLS.leftBuf[2][1] = g_OLS.leftBuf[1][2]; } #endif } /******************************************************************************* * xxx@2019-12-13 * Function Name : OLS_CalcResult * Description : 计算结果，计算最终的拟合公式的各个系数 * Input : None * Output : None * Return : None *******************************************************************************/ void OLS_CalcResult(void) { float A; //计算A的行列式 #if (OLS_ORDER_NUM == 1) { /* 计算二元一次方程组： ax + by = m cx + dy = n 通解： x = (md - bn)/(ad - bc) y = (mc - an)/(bc - ad) */ g_OLS.resultBuf[0] = (g_OLS.righeBuf[0] * g_OLS.leftBuf[1][1] - g_OLS.leftBuf[0][1] * g_OLS.righeBuf[1]) / \ (g_OLS.leftBuf[0][0] * g_OLS.leftBuf[1][1] - g_OLS.leftBuf[0][1] * g_OLS.leftBuf[1][0]); g_OLS.resultBuf[1] = (g_OLS.righeBuf[0] * g_OLS.leftBuf[1][0] - g_OLS.leftBuf[0][0] * g_OLS.righeBuf[1]) / \ (g_OLS.leftBuf[0][1] * g_OLS.leftBuf[1][0] - g_OLS.leftBuf[0][0] * g_OLS.leftBuf[1][1]); } #elif (OLS_ORDER_NUM == 2) { /* 计算三元一次方程组： ax + by + cz = α dx + ey + fz = β gx + hy + iz = γ 通解： A的行列式|A| = a(ei - hf) + b(fg - id) + c(dh - ge) x = (α(ei - hf) + β(hc - bi) + γ(bf - ec))/|A| y = (α(gf - di) + β(ai - gc) + γ(dc - af))/|A| z = (α(dh - ge) + β(gb - ah) + γ(ae - db))/|A| */ A = g_OLS.leftBuf[0][0] * (g_OLS.leftBuf[1][1] * g_OLS.leftBuf[2][2] - g_OLS.leftBuf[2][1] * g_OLS.leftBuf[1][2]) + \ g_OLS.leftBuf[0][1] * (g_OLS.leftBuf[1][2] * g_OLS.leftBuf[2][0] - g_OLS.leftBuf[2][2] * g_OLS.leftBuf[1][0]) + \ g_OLS.leftBuf[0][2] * (g_OLS.leftBuf[1][0] * g_OLS.leftBuf[2][1] - g_OLS.leftBuf[2][0] * g_OLS.leftBuf[1][1]); g_OLS.resultBuf[0] = (g_OLS.righeBuf[0] * (g_OLS.leftBuf[1][1] * g_OLS.leftBuf[2][2] - g_OLS.leftBuf[2][1] * g_OLS.leftBuf[1][2]) + \ g_OLS.righeBuf[1] * (g_OLS.leftBuf[2][1] * g_OLS.leftBuf[0][2] - g_OLS.leftBuf[0][1] * g_OLS.leftBuf[2][2]) + \ g_OLS.righeBuf[2] * (g_OLS.leftBuf[0][1] * g_OLS.leftBuf[1][2] - g_OLS.leftBuf[1][1] * g_OLS.leftBuf[0][2])) / A; g_OLS.resultBuf[1] = (g_OLS.righeBuf[0] * (g_OLS.leftBuf[2][0] * g_OLS.leftBuf[1][2] - g_OLS.leftBuf[1][0] * g_OLS.leftBuf[2][2]) + \ g_OLS.righeBuf[1] * (g_OLS.leftBuf[0][0] * g_OLS.leftBuf[2][2] - g_OLS.leftBuf[2][0] * g_OLS.leftBuf[0][2]) + \ g_OLS.righeBuf[2] * (g_OLS.leftBuf[1][0] * g_OLS.leftBuf[0][2] - g_OLS.leftBuf[0][0] * g_OLS.leftBuf[1][2])) / A; g_OLS.resultBuf[2] = (g_OLS.righeBuf[0] * (g_OLS.leftBuf[1][0] * g_OLS.leftBuf[2][1] - g_OLS.leftBuf[2][0] * g_OLS.leftBuf[1][1]) + \ g_OLS.righeBuf[1] * (g_OLS.leftBuf[2][0] * g_OLS.leftBuf[0][1] - g_OLS.leftBuf[0][0] * g_OLS.leftBuf[2][1]) + \ g_OLS.righeBuf[2] * (g_OLS.leftBuf[0][0] * g_OLS.leftBuf[1][1] - g_OLS.leftBuf[1][0] * g_OLS.leftBuf[0][1])) / A; } #endif } int main(void) { OLS_Init(); OLS_CalcBufData(); OLS_CalcResult(); #if (OLS_ORDER_NUM == 1) { printf("a0 = %f\n\r", g_OLS.resultBuf[0]); printf("a1 = %f\n\r", g_OLS.resultBuf[1]); } #elif (OLS_ORDER_NUM == 2) { printf("a0 = %f\n\r", g_OLS.resultBuf[0]); printf("a1 = %f\n\r", g_OLS.resultBuf[1]); printf("a2 = %f\n\r", g_OLS.resultBuf[2]); } #endif return 0; } 当然你们也可以直接用矩阵的运算，直接多少阶都可以直接计算。那样就不用找多元方程组的通解怎么解了。（注意float的精度问题，精度不够可以用double）
附上代码和Word文档： | 
/2 
|手机版|OpenEdv-开源电子网
( 粤ICP备12000418号-1 )
发表于 2019-12-13 19:20:49
楼主
