基于 C 的高斯曲线拟合原理以及实现 [Transfer
https://blog.****.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.****.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),获取详细算法。
原文地址:https://www.cnblogs.com/sggggr/p/14247183.html
推荐阅读
-
基于 C 的高斯曲线拟合原理以及实现 [Transfer
-
windows下进程间通信的(13种方法)-摘 要 本文讨论了进程间通信与应用程序间通信的含义及相应的实现技术,并对这些技术的原理、特性等进行了深入的分析和比较。 ---- 关键词 信号 管道 消息队列 共享存储段 信号灯 远程过程调用 Socket套接字 MQSeries 1 引言 ---- 进程间通信的主要目的是实现同一计算机系统内部的相互协作的进程之间的数据共享与信息交换,由于这些进程处于同一软件和硬件环境下,利用操作系统提供的的编程接口,用户可以方便地在程序中实现这种通信;应用程序间通信的主要目的是实现不同计算机系统中的相互协作的应用程序之间的数据共享与信息交换,由于应用程序分别运行在不同计算机系统中,它们之间要通过网络之间的协议才能实现数据共享与信息交换。进程间通信和应用程序间通信及相应的实现技术有许多相同之处,也各有自己的特色。即使是同一类型的通信也有多种的实现方法,以适应不同情况的需要。 ---- 为了充分认识和掌握这两种通信及相应的实现技术,本文将就以下几个方面对这两种通信进行深入的讨论:问题的由来、解决问题的策略和方法、每种方法的工作原理和实现、每种实现方法的特点和适用的范围等。 2 进程间的通信及其实现技术 ---- 用户提交给计算机的任务最终都是通过一个个的进程来完成的。在一组并发进程中的任何两个进程之间,如果都不存在公共变量,则称该组进程为不相交的。在不相交的进程组中,每个进程都独立于其它进程,它的运行环境与顺序程序一样,而且它的运行环境也不为别的进程所改变。运行的结果是确定的,不会发生与时间相关的错误。 ---- 但是,在实际中,并发进程的各个进程之间并不是完全互相独立的,它们之间往往存在着相互制约的关系。进程之间的相互制约关系表现为两种方式: ---- (1) 间接相互制约:共享CPU ---- (2) 直接相互制约:竞争和协作 ---- 竞争——进程对共享资源的竞争。为保证进程互斥地访问共享资源,各进程必须互斥地进入各自的临界段。 ---- 协作——进程之间交换数据。为完成一个共同任务而同时运行的一组进程称为同组进程,它们之间必须交换数据,以达到协作完成任务的目的,交换数据可以通知对方可以做某事或者委托对方做某事。 ---- 共享CPU问题由操作系统的进程调度来实现,进程间的竞争和协作由进程间的通信来完成。进程间的通信一般由操作系统提供编程接口,由程序员在程序中实现。UNIX在这个方面可以说最具特色,它提供了一整套进程间的数据共享与信息交换的处理方法——进程通信机制(IPC)。因此,我们就以UNIX为例来分析进程间通信的各种实现技术。 ---- 在UNIX中,文件(File)、信号(Signal)、无名管道(Unnamed Pipes)、有名管道(FIFOs)是传统IPC功能;新的IPC功能包括消息队列(Message queues)、共享存储段(Shared memory segment)和信号灯(Semapores)。 ---- (1) 信号 ---- 信号机制是UNIX为进程中断处理而设置的。它只是一组预定义的值,因此不能用于信息交换,仅用于进程中断控制。例如在发生浮点错、非法内存访问、执行无效指令、某些按键(如ctrl-c、del等)等都会产生一个信号,操作系统就会调用有关的系统调用或用户定义的处理过程来处理。 ---- 信号处理的系统调用是signal,调用形式是: ---- signal(signalno,action) ---- 其中,signalno是规定信号编号的值,action指明当特定的信号发生时所执行的动作。 ---- (2) 无名管道和有名管道 ---- 无名管道实际上是内存中的一个临时存储区,它由系统安全控制,并且独立于创建它的进程的内存区。管道对数据采用先进先出方式管理,并严格按顺序操作,例如不能对管道进行搜索,管道中的信息只能读一次。 ---- 无名管道只能用于两个相互协作的进程之间的通信,并且访问无名管道的进程必须有共同的祖先。 ---- 系统提供了许多标准管道库函数,如: pipe——打开一个可以读写的管道; close——关闭相应的管道; read——从管道中读取字符; write——向管道中写入字符; ---- 有名管道的操作和无名管道类似,不同的地方在于使用有名管道的进程不需要具有共同的祖先,其它进程,只要知道该管道的名字,就可以访问它。管道非常适合进程之间快速交换信息。 ---- (3) 消息队列(MQ) ---- 消息队列是内存中独立于生成它的进程的一段存储区,一旦创建消息队列,任何进程,只要具有正确的的访问权限,都可以访问消息队列,消息队列非常适合于在进程间交换短信息。 ---- 消息队列的每条消息由类型编号来分类,这样接收进程可以选择读取特定的消息类型——这一点与管道不同。消息队列在创建后将一直存在,直到使用msgctl系统调用或iqcrm -q命令删除它为止。 ---- 系统提供了许多有关创建、使用和管理消息队列的系统调用,如: ---- int msgget(key,flag)——创建一个具有flag权限的MQ及其相应的结构,并返回一个唯一的正整数msqid(MQ的标识符); ---- int msgsnd(msqid,msgp,msgsz,msgtyp,flag)——向队列中发送信息; ---- int msgrcv(msqid,cmd,buf)——从队列中接收信息; ---- int msgctl(msqid,cmd,buf)——对MQ的控制操作; ---- (4) 共享存储段(SM) ---- 共享存储段是主存的一部分,它由一个或多个独立的进程共享。各进程的数据段与共享存储段相关联,对每个进程来说,共享存储段有不同的虚拟地址。系统提供的有关SM的系统调用有: ---- int shmget(key,size,flag)——创建大小为size的SM段,其相应的数据结构名为key,并返回共享内存区的标识符shmid; ---- char shmat(shmid,address,flag)——将当前进程数据段的地址赋给shmget所返回的名为shmid的SM段; ---- int shmdr(address)——从进程地址空间删除SM段; ---- int shmctl (shmid,cmd,buf)——对SM的控制操作; ---- SM的大小只受主存限制,SM段的访问及进程间的信息交换可以通过同步读写来完成。同步通常由信号灯来实现。SM非常适合进程之间大量数据的共享。 ---- (5) 信号灯 ---- 在UNIX中,信号灯是一组进程共享的数据结构,当几个进程竞争同一资源时(文件、共享内存或消息队列等),它们的操作便由信号灯来同步,以防止互相干扰。 ---- 信号灯保证了某一时刻只有一个进程访问某一临界资源,所有请求该资源的其它进程都将被挂起,一旦该资源得到释放,系统才允许其它进程访问该资源。信号灯通常配对使用,以便实现资源的加锁和解锁。 ---- 进程间通信的实现技术的特点是:操作系统提供实现机制和编程接口,由用户在程序中实现,保证进程间可以进行快速的信息交换和大量数据的共享。但是,上述方式主要适合在同一台计算机系统内部的进程之间的通信。 3 应用程序间的通信及其实现技术 ---- 同进程之间的相互制约一样,不同的应用程序之间也存在竞争和协作的关系。UNIX操作系统也提供一些可用于应用程序之间实现数据共享与信息交换的编程接口,程序员可以通过自己编程来实现。如远程过程调用和基于TCP/IP协议的套接字(Socket)编程。但是,相对普通程序员来说,它们涉及的技术比较深,编程也比较复杂,实现起来困难较大。 ---- 于是,一种新的技术应运而生——通过将有关通信的细节完全掩盖在某个独立软件内部,即底层的通讯工作和相应的维护管理工作由该软件内部来实现,用户只需要将通信任务提交给该软件去完成,而不必理会它的具体工作过程——这就是所谓的中间件技术。 ---- 我们在这里分别讨论这三种常用的应用程序间通信的实现技术——远程过程调用、会话编程技术和MQSeries消息队列技术。其中远程过程调用和会话编程属于比较低级的方式,程序员参与的程度较深,而MQSeries消息队列则属于比较高级的方式,即中间件方式,程序员参与的程度较浅。 ---- 4.1 远程过程调用(RPC)