基于C语言的高斯曲线拟合原理以及实现【转】
https://blog.csdn.net/dingzj2000/article/details/103719368?utm_medium=distribute.pc_relevant.none-task-blog-OPENSEARCH-3.control&depth_1-utm_source=distribute.pc_relevant.none-task-blog-OPENSEARCH-3.control
1.意义
高斯曲线 ，又叫做gaussian curve，是正态分布中的一条标准曲线。具有以下特征：
1.1 正态曲线在横轴上方均数处最高；
1.2 正在分布以均数为中心，左右对称；
1.3 正态分布有两个参数，即均数和标准差；标准正态分布用N(0，1)表示；
1.4 正态曲线下的面积分布有一定的规律。
在分析仪器的测量中，有许多具有明确的物理意义的二维图谱，如光谱图、色谱图等，许多测量图谱都可以用高斯曲线予以描述。高斯曲线虽然也是非线性函数，但它的各个参数具有明确的物理意义，因为高斯拟合在分析仪器的测量中具有广泛的应用前景。利用它来描述或拟合求出一些实验数据的分析，往往能起到常规方法不能达到的作用。
2.结果展示
最近在做一个医用仪器项目，需要用到高斯拟合，之前也没有搞过相关内容，于是先在PC端实现高斯拟合函数，能够实现要求后，再移植到仪器上使用。
先展示一下结果：
图2.1 高斯曲线拟合
先随机选6个测试点（蓝色点），根据这6个测试点，进行高斯拟合，红色曲线就是拟合出来的曲线。拟合出来的曲线基本在选取的6个测试点附近。通过这6个点，找出了互相之间的关系。达到了设计目的。
3.原理
高斯拟合即使用形如：Gi(x) = Ai*exp((x-Bi)^2/Ci^2)的高斯函数对数据点集进行函数逼近的拟合方法，高斯拟合跟多项式拟合类似，不同的是多项式拟合是用幂函数系，而高斯拟合用的是高斯函数系。使用高斯函数拟合来进行拟合，优点在于计算积分十分简单快捷。
3.1 高斯函数
高斯函数：
a表示得到曲线的高度，c是指曲线在x轴的中心，b指width(与半峰全宽有关),图形如下：图形如下：
3.2 高斯拟合原理
设有一组实验数据(i = 1,2,3,...N),可用高斯函数描述：
式（3.1）中待估参数a,c和b，分别代表的物理意义为高斯曲线的峰高、峰位置和半宽度信息。将式(3.1)两边取自然对数，化为：
令
则式（3.2）化为二次多项式拟合函数
考虑全部数据和测量误差，并以矩阵形式表示如下
简记为：
在不考虑总量程误差E的影响情况下，根据最小二乘原理，可求得拟合常数b0，b1，b2构成的矩阵B的广义最小二乘解为：
进而根据式（3.3），可求得待估参数a，b，c：
，
，
(3.7)
将所得的a，b和c带入式（3.1），就能够得到由实验数据(i = 1,2,3,...N)拟合的高斯函数。
3.3 最小二乘法
最小二乘法（又称最小平方法）是一种数学优化技术。它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。
1801年，意大利天文学家朱赛普·皮亚齐发现了第一颗小行星谷神星。经过40天的跟踪观测后，由于谷神星运行至太阳背后，使得皮亚齐失去了谷神星的位置。随后全世界的科学家利用皮亚齐的观测数据开始寻找谷神星，但是根据大多数人计算的结果来寻找谷神星都没有结果。时年24岁的高斯也计算了谷神星的轨道。奥地利天文学家海因里希·奥尔伯斯根据高斯计算出来的轨道重新发现了谷神星。
高斯使用的最小二乘法的方法发表于1809年他的著作《天体运动论》中。法国科学家勒让德于1806年独立发明“最小二乘法”，但因不为世人所知而默默无闻。勒让德曾与高斯为谁最早创立最小二乘法原理发生争执。1829年，高斯提供了最小二乘法的优化效果强于其他方法的证明，因此被称为高斯-马尔可夫定理。
4.实现
4.1 源码
这里先把C代码展示出来：
/* 2019.12.24 验证高斯拟合代码 环境:Ubuntu32 16.04.1 cairo库实现 */ #include <stdlib.h> #include <stdio.h> #include <string.h> #include <cairo.h> #include <math.h> #include "GuassFitting.h" /* 采集样本数量 */ #define N 6 /* 原点坐标 */ #define OriginX 20 #define OriginY 20 /* 坐标转换 */ #define X(n) ((n)*1.0) #define Y(n) ((400-(n))*1.0) /* 坐标系内位置 */ #define PosX(n) (X(OriginX+(n))) #define PosY(n) (Y(OriginY+(n))) #define ANGLE(ang) (ang * 3.1415926 / 180.0) struct POS{ int x; int y; }; cairo_surface_t *image_surface_create_from_png(const char *filename) { cairo_status_t cst; cairo_surface_t *image_sf=cairo_image_surface_create_from_png(filename); cst = cairo_surface_status (image_sf); if (cst!=CAIRO_STATUS_SUCCESS) { printf( "failed to cairo_image_surface_create_from_png cairo_status_t is:%d file: %s",cst, filename); image_sf = NULL; //if (cst == CAIRO_STATUS_NO_MEMORY) { //image_sf = cairo_image_surface_create_from_jpeg(filename); //} } return image_sf; } void draw_png2surface(cairo_t *cr, double x, double y, cairo_surface_t *surface){ if(surface != NULL){ cairo_set_source_surface(cr, surface, x, y); cairo_paint(cr); } } /* 创建背景 */ void createBackground(cairo_t *cr,char *file) { cairo_surface_t *g_background; if(cr == NULL || file == NULL) { return; } /* 背景图 */ g_background = image_surface_create_from_png(file); draw_png2surface(cr, 0, 0, g_background); } /* 构建坐标系 */ void createCoordinate(cairo_t *cr) { if(cr == NULL) { return; } /* 划线 */ cairo_set_source_rgb (cr, 0.0, 0.0, 0.0);/* 设置颜色 -黑色 */ cairo_set_line_cap (cr, CAIRO_LINE_CAP_ROUND); cairo_set_line_width (cr, 1.0); cairo_move_to (cr, X(OriginX), Y(OriginY));//X轴 cairo_line_to (cr, X(OriginX+650), Y(OriginY)); cairo_stroke (cr); cairo_move_to (cr, X(OriginX), Y(OriginY));//Y轴 cairo_line_to (cr, X(OriginX), Y(OriginY+350)); cairo_stroke (cr); /* 量程 */ cairo_select_font_face (cr, "serif", CAIRO_FONT_SLANT_NORMAL, CAIRO_FONT_WEIGHT_BOLD); cairo_set_font_size (cr, 15.0); //原点 cairo_move_to (cr, X(OriginX-5), Y(OriginY-15)); cairo_show_text (cr, "0"); //X轴 cairo_move_to (cr, X(OriginX+650-10), Y(OriginY-15)); cairo_show_text (cr, "650"); //Y轴 cairo_move_to (cr, X(OriginX-10), Y(OriginY+350+5)); cairo_show_text (cr, "350"); cairo_stroke (cr); } /* 构建坐标点 */ void createPoints(cairo_t *cr,int x,int y) { char buf[64]; if(cr == NULL) { return; } cairo_set_source_rgb(cr, 0.0, 0.0, 1.0);/* 设置颜色 -蓝色 */ cairo_set_line_width(cr, 4); cairo_arc(cr, PosX(x), PosY(y), 2, ANGLE(0), ANGLE(360)); cairo_stroke (cr); /* 显示坐标 */ memset(buf,0x00,sizeof(buf)); sprintf(buf,"(%d,%d)",x,y); cairo_set_source_rgb(cr, 0.0, 0.0, 0.0);/* 设置颜色 -黑色 */ cairo_select_font_face (cr, "serif", CAIRO_FONT_SLANT_NORMAL, CAIRO_FONT_WEIGHT_BOLD); cairo_set_font_size (cr, 10.0); cairo_move_to(cr,PosX(x-8), PosY(y+8)); cairo_show_text (cr, buf); cairo_stroke (cr); } /* 创建高斯曲线 */ void createGuassCurve(cairo_t *cr,double a,double b,double c) { int x,y; if(cr == NULL) { return; } cairo_set_source_rgb(cr, 1.0, 0.0, 0.0);/* 设置颜色 -红色 */ cairo_set_line_width (cr, 1.0); for(x = 0; x < 630; x++) { y = a*exp(-(pow(x-c,2.0)/b)); cairo_move_to(cr, PosX(x), PosY(y)); cairo_line_to(cr, PosX(x), PosY(y)); cairo_stroke (cr); } } /* 指数和对数函数测试 */ int mathTest(void) { printf("pow(x,y) x= 10,y=2 value=%lf\n",pow(10.0,2.0)); printf("powl(x,y) log x=10,y=100 value=%lf\n",powf(100.0,2.0)); printf("exp(x) e x=1 value=%f\n",exp(1)); printf("exp(x) e x=2 value=%f\n",exp(2)); printf("loge=%f\n",log(10)); //以e为底的对数函数 printf("loge=%f\n",log(2.718282)); //以e为底的对数函数 printf("log10=%f\n",log10(100)); //以10为底的对数函数 printf("sqrt=%f\n",sqrt(16));//平方根 } int main(int argc,char *argv[]) { int i; double b[N*4],a[N*3*4], q[N*4*4];//对应缓存扩大4倍，防止core dumped double b0,b1,b2; double ga,gb,gc; int m,n; /* 采样点 */ #if 0 struct POS pos[N] = { {200,100}, {300,220}, {400,300}, {500,220} }; #else struct POS pos[N] = { {200,200}, {300,280}, {400,320}, {440,310}, {500,280}, {600,220} }; #endif /************************ 创建cairo ****************************/ cairo_surface_t *surface = cairo_image_surface_create (CAIRO_FORMAT_ARGB32, 700, 400); cairo_t *cr = cairo_create (surface); /*********************** 创建背景 *****************************/ createBackground(cr,"background.png"); /*********************** 构建坐标系 *****************************/ createCoordinate(cr); /************************ 绘制采样点 *****************************/ for(i = 0; i < N; i++) { createPoints(cr,pos[i].x,pos[i].y); } /********************* 对采样点进行数据预处理 *******************/ for(i = 0;i < N;i++) { b[i] = log(pos[i].y);//Z } for(i = 0;i < N;i++) { a[i*3+0] = 1.0; a[i*3+1] = pos[i].x*1.0; a[i*3+2] = pow(pos[i].x,2)*1.0; } m = N;n = 3; /*********************** 高斯拟合 最小二乘解 ********************/ if(GuassFitting_Gmqr(a,m,n,b,q)==0) { printf("GuassFitting Fail!\n"); } else { for (i=0; i<n; i++) { printf("b%d=%13.7f\n",i,b[i]); } b0 = b[0];b1 = b[1];b2 = b[2];//最小二乘解系数 /*********************** 计算高斯参数 ********************/ GuassFitting_GetCurvePara(b0,b1,b2,&ga,&gb,&gc); printf("a = %f\n",ga); printf("b = %f\n",gb); printf("c = %f\n",gc); /************************ 绘制高斯曲线 ****************************/ createGuassCurve(cr,ga,gb,gc); } /************************ 创建图片 ****************************/ cairo_surface_write_to_png (surface, "guassfitting.png"); /************************ 销毁cairo ****************************/ cairo_destroy (cr); cairo_surface_destroy (surface); return 0; }
4.2 cairo库
基于Linux C编程，图像显示使用cairo库。至于如何使用不在本文论述的重点。详见：
https://blog.csdn.net/dingzj2000/article/details/103719104
4.3 数据预处理
采样的坐标分别为(200,200),(300,280),(400,320),(440,310),(500,280),(600,220)
基于A=XB形式，由式（3.3），对于A来说.
存在将y位置转换为对数形式。所以：b[i] = log(pos[i].y);//Z
对于X的形式是[1 x x^2]，所以：
a[i*3+0] = 1.0;
a[i*3+1] = pos[i].x*1.0;
a[i*3+2] = pow(pos[i].x,2)*1.0;
4.4 函数参数定义
int GuassFitting_Gmqr(double a[],int m,int n,double b[],double q[])
作用：
函数语句与形参说明最小二乘问题的豪斯荷尔德变化法，输入矩阵相关参数；获取最小二乘解
参数:
double a[]:存放超定方程组的系数矩阵A。返回时存放QR分解式中的R矩阵
理解为a[m][n],以一位数组展示出来,组合成一位数组
int m: 系数矩阵A的行数,m>=n，获取的采样点数
int n: 系数矩阵A的列数,n<=m ，固定为3
doube b[]:存放方程组右端B的常数向量。返回时前n个分量存放方程组的最小二乘解
doube q[]:返回时存放QR的分解式中的正交矩阵Q。理解为q[m][n],以一维数组表现出来。
返回:
返回0，则表示程序工作失败(如A列线性相关);
若返回的标志值不为0，则表示正常返回。
b0 = b[0];b1 = b[1];b2 = b[2];//最小二乘解系数
int GuassFitting_GetCurvePara(double b0,double b1,double b2,double *a,double *b,double *c)
此函数就是实现式（3.7），通过b0，b1和b2获取a，b，c最终获取高斯曲线函数
4.5 参数打印
打印参数如下：
通过上述打印清楚的看到B的系数，以及高斯曲线参数。通过值发现，在x=409.360814，是高斯曲线的峰值。
5.核心算法
核心高斯拟合算法采用C语言实现,见下图：
6.获取算法
加微信(微信号：dingzj2000)，获取详细算法。
