最近在做线激光测距的工作, 需要设计算法查找线激光中心, 于是用到了一些拟合算法, 下面介绍一下用到的高斯函数拟合.
高斯函数拟合的目标函数为
$$y = y_{\max} \exp{\left(-\frac{(x - x_{\max})^2}{S}\right)}$$
式中待估计参数 $y_{\max}$, $x_{\max}$ 和 $S$ 分别为高斯曲线的峰值, 峰值位置和半宽度信息. 对上式两边取自然对数, 化为,
$$\begin{aligned} \ln y &= \ln y_{\max} - \frac{(x - x_{\max})^2}{S} \\ &= \left( \ln y_{\max} - \frac{x_{\max}^2}{S} \right) + \frac{2 x_{\max} x}{S} - \frac{x^2}{S} \end{aligned}$$
令 $z=\ln y$, $b_0 = \ln y_{\max} - \frac{x_{\max}^2}{S}$, $b_1 = \frac{2 x_{\max}}{S}$, $b_2=- \frac{1}{S}$, 所以目标函数转化为:
$$ z = b_2 x^2 + b_1 x + b_0$$
已知点集 $\{(x_1,y_1), (x_2,y_2), \dots, (x_n,y_n)\}$, $n\geq 3$, 先按上面的 $z=\ln y$ 将其转化为 $\{(x_1,z_1), (x_2,z_2), \dots, (x_n,z_n)\}$, 最小二乘拟合的目标函数为:
$$
\sum\nolimits_{i = 1}^{n} {{\left( {b_2{x_i}^2 + b_1{x_i} + b_0 - {z_i}} \right)}^2} = {\left\| {X\vec \theta - \vec z} \right\|^2}
$$
其中
$X = \begin{bmatrix}
1 & x_1 & x_1^2 \\
1 & x_2 & x_2^2 \\
\vdots & \vdots & \vdots \\
1 & x_n & x_n^2
\end{bmatrix}$ (称为 Vandermonde 矩阵), $\vec \theta = \begin{bmatrix} b_0 \\ b_1 \\ b_2 \end{bmatrix}$ 和 $\vec z = \begin{bmatrix} z_1 \\ z_2 \\ \vdots \\ z_n \end{bmatrix}$.
由 $b_2=- \frac{1}{S}$ 可得 $S=-\frac{1}{b_2}$.
再由 $b_1 = \frac{2 x_{\max}}{S}$ 可得 $x_{\max}=\frac{Sb_1}{2}=-\frac{b_1}{2b_2}$.
再由 $b_0 = \ln y_{\max} - \frac{x_{\max}^2}{S}$ 可得 $y_{\max} =\exp{\left( b_0 + \frac{x_{\max}^2 }{S}\right)} =\exp{\left( b_0 - \frac{b_1^2 }{4b_2}\right)}$.
测试代码:
import numpy as np
import matplotlib.pyplot as plt
def gaussian_curve_fitting(x, y):
"""高斯曲线拟合
"""
vander_mat = np.vander(x, 3, increasing=True)
rep_coeffs = np.linalg.lstsq(vander_mat, np.log(y), rcond=-1)[0]
coeffs = np.empty_like(rep_coeffs)
coeffs[2] = S = -1 / rep_coeffs[2]
coeffs[1] = x_max = S * rep_coeffs[1] / 2
coeffs[0] = y_max = np.exp(rep_coeffs[0] + x_max **2 / S)
return coeffs
if __name__ == '__main__':
x = np.linspace(-5, 5, 100)
y_max, x_max, S = 50, 1.1, 3.0
y = y_max * np.exp(-(x - x_max)**2 / S)
y += 0.001 * np.abs(np.random.randn(*y.shape))
coeffs = gaussian_curve_fitting(x, y)
print(coeffs)
y_hat = coeffs[0] * np.exp(-(x - coeffs[1] )**2 / coeffs[2])
plt.scatter(x, y, c='b', s=2)
plt.plot(x, y_hat, c='g')
plt.show()