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)，獲取詳細算法。