我不喜欢 Matlab，所以我不知道 Matlab 拟合做了哪些额外的工作来估计非线性拟合的起始值。不过，我可以说，这curve_fit根本不是，即所有值都假定为1. x最简单的方法是将轴重新缩放到 range [0, 2 pi]。因此，OP 的问题再次是错误的起始值。然而，重新缩放需要知道要拟合的主波大约是数据集的宽度。此外，我们需要假设所有其他拟合参数也是顺序的1。幸运的是，情况就是这样，所以这会奏效：
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
xdat, ydat = np.loadtxt( "data.tsv", unpack=True, skiprows=1 )
def fourier(x, *as_bs):
sum_a = 0
sum_b = 0
j = 1
w = as_bs[0]
a0 = as_bs[1]
for i in range(2, len( as_bs ) - 1, 2 ):
sum_a += as_bs[i] * np.cos( j * w * x )
sum_b += as_bs[i+1] * np.sin( j * w * x )
j = j + 1
return a0 + sum_a + sum_b
"""
lets rescale the data to get the base frequency in the range of one
"""
xmin = min( xdat )
xmax = max( xdat )
xdat = ( xdat - xmin ) / (xmax - xmin ) * 2 * np.pi
popt, pcov = curve_fit(
fourier,
xdat, ydat,
p0 = np.ones(8)
)
### here I assume that higher order are similar to lower orders
### but slightly smaller. ... hoping that the fit correts errors in
### this assumption
print(popt)
### scale back w noting that it scales inverse to x
print( popt[0] * 2 * np.pi / (xmax - xmin ) )
data_fit = fourier( xdat, *popt )
如果我们不能做出上述假设，我们只能假设存在一个对信号有主要贡献的基频（请注意，这并不总是正确的）。在这种情况下，我们可以以非迭代方式预先计算起始猜测。
解决方案看起来有点复杂：
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from scipy.integrate import cumtrapz
xdat, ydat = np.loadtxt( "data.tsv", unpack=True, skiprows=1 )
def fourier(x, *as_bs):
sum_a = 0
sum_b = 0
j = 1
w = as_bs[0]
a0 = as_bs[1]
for i in range(2, len( as_bs ) - 1, 2 ):
sum_a += as_bs[i] * np.cos( j * w * x )
sum_b += as_bs[i+1] * np.sin( j * w * x )
j = j + 1
return a0 + sum_a + sum_b
#### initial guess
"""
This uses the fact that if y = a sin w t + b cos w t + c we have
int int y = -y/w^2 + c/2 t^2 + d t + e
i.e. we can get 1/w^2 as linear fit parameter without the danger of
a non-linear fit iterative process running into a local minimum
for details see:
https://scikit-guess.readthedocs.io/en/sine/_downloads/4b4ed1e691ff195be3ca73879a674234/Regressions-et-equations-integrales.pdf
"""
Sy = cumtrapz( ydat, xdat, initial=0 )
SSy = cumtrapz( Sy, xdat, initial=0 )
ST = np.array( [
ydat, xdat**2, xdat, np.ones( len( xdat ) )
] )
S = np.transpose( ST )
eta = np.dot( ST, SSy )
A = np.dot( ST, S )
sol = np.linalg.solve( A, eta )
wFit = np.sqrt( -1 / sol[0] )
### linear parameters
"""
Once we have a good guess for w we can get starting guesses for
a, b and c from a standard linear fit
"""
ST = np.array( [
np.sin( wFit * xdat ), np.cos( wFit * xdat ), np.ones( len( xdat ) )
])
S = np.transpose( ST )
eta = np.dot( ST, ydat )
A = np.dot( ST, S )
sol = np.linalg.solve( A, eta )
a1 = sol[0]
b1 = sol[1]
a0 = sol[2]
### final non-linear fit
"""
Now we can use the guesses from above as input for the final
non-linear fit. Hopefully, we are now close enough to the global minimum
and have the algorithm converge reasonably
"""
popt, pcov = curve_fit(
fourier,
xdat, ydat,
p0=[
wFit, a0, a1, b1,
a1 / 2, b1 / 2,
a1 / 4, b1 / 4
]
)
### here I assume that higher order are similar to lower orders
### but slightly smaller. ... hoping that the fit correts errors in
### this assumption
print(popt)
data_fit = fourier( xdat, *popt )
plt.plot( xdat, ydat, ls="", marker="o", ms=0.5, label="data" )
plt.plot( xdat, data_fit, label='fitting')
plt.legend()
plt.show()
两者都提供基本相同的解决方案，后一种代码适用于更多情况，假设更少。
