轉自:https://blog.csdn.net/renjie10/article/details/114933766
1.目的
針對光譜離散數據,尋峰完成后截取near峰值的數據,利用高斯擬合重繪單峰曲線,進而實現分峰功能
2.原理
3.代碼
import numpy as np from math import log, exp import matplotlib.pyplot as plt from scipy.optimize import curve_fit from scipy import asarray as ar,exp # 將txt文件讀入numpy數組 yOriginal = np.loadtxt('C:\\工作\數據.txt') #一維數據 #yOriginal = np.array([5.81528E-05, 0.000111682, 0.000271214, 0.000391546, 0.000786933, 0.002034528, 0.002968284, 0.005004177, 0.007329225, 0.011119662, 0.017025547, 0.02488255, 0.04219861, 0.040429801, 0.035320014, 0.05154864, 0.06894745, # 0.105841984, 0.083166325, 0.110311517, 0.055681743, 0.093540639, 0.066621081, 0.056688568, 0.045128754, 0.045911571, 0.028179728, 0.021262112, 0.018781554, 0.008240159, 0.008562607, 0.004372914, 0.002847578, 0.001717186, 0.001081616]) xOriginal = np.arange(len(yOriginal)) print("輸入數據n", yOriginal) average = sum(yOriginal)/len(yOriginal) #y = yOriginal #x = xOriginal #print("過濾數據n", y) y = np.log(y) zMatrix = np.matrix(y) print("取對數n", y) # 構造 X 矩陣 xMatrixT = np.matrix(np.reshape(np.concatenate( (np.ones((len(y))), x, x*x)), (3, len(y)))) xMatrix = np.matrix(xMatrixT.T) print("X 矩陣n", xMatrix) #np.matrix.__mul__ = np.dot # 重載運算符 print(xMatrixT*xMatrix) bMatrix = ((xMatrixT*xMatrix).I*xMatrixT)*zMatrix.T#矩陣運算 print("B 矩陣n", bMatrix) b2, b1, b0 = float(bMatrix[2][0]), float(bMatrix[1][0]), float(bMatrix[0][0]) print("b0={}b1={}b2={}".format(b0, b1, b2)) s = -1/b2 xMaxi = s*b1/2 yMaxi = exp(b0+xMaxi**2/s) print(yMaxi, xMaxi, s) def gaussian(x,*param): return param[0]*np.exp(-np.power(x - param[1], 2.) / (2 * np.power(param[2], 2.)))#高斯公式 popt,pcov = curve_fit(gaussian,xOriginal,yOriginal,p0=[yMaxi,xMaxi,s]) print(popt) print(pcov) plt.plot(xOriginal,yOriginal,'b+:',label='data') plt.plot(xOriginal,gaussian(xOriginal,*popt),'ro:',label='fit')#繪圖 plt.legend() plt.show()
其中傳參數時s應改為 np.sqrt(s/2)
4.效果
參考:
最小二乘法: https://www.zhihu.com/question/37031188
高斯曲線擬合公式推導: https://blog.csdn.net/Abaqus3_0/article/details/82011159
1