在使用边界时，我如何优化我的鼻窦拟合?

Question

我正试图在曲线的给定部分上做一个正弦拟合。这里有两个限制。首先，我拟合的正弦曲线的偏移量应该是0，其次，我拟合的正弦曲线的振幅应该与原始数据的最小值相同。

当我使用下面的代码时，拟合看起来像我添加(1)的图片。在我看来，正弦函数的周期应该更高。拟合曲线与原始数据的最小匹配，拟合曲线宽度不够。当我不使用c和A的边界时，我的拟合看起来很好(2)，我做错了什么?是否有一种方法可以修改拟合，使正弦曲线在使用a和c的边界时更适合?

无界拟合

用边界拟合

编辑:我提到的一件事是，拟合非常依赖于振幅的起始值(ff_guess)。当我手动将其更改为X时(例如。10或30)，而拟合的窦性曲线总是显示在X附近的振幅(10.3或32.5)。是否有任何我还没有考虑的设置，可以阻止优化器改变振幅值?

import matplotlib.pyplot as plt
from matplotlib.ticker import FormatStrFormatter
from matplotlib.pyplot import figure, rcParams
import numpy as np
from scipy.optimize import curve_fit

#Time
t = [313.544, 313.545, 313.546, 313.547, 313.548, 313.549, 313.55,  313.551, 313.552, 313.553, 313.554, 313.555, 313.556, 313.557, 313.558, 313.559, 313.56,  313.561, 313.562, 313.563, 313.564, 313.565, 313.566, 313.567,]
#y-Values
s = [0.911188, -0.43135, -1.80997, -3.27816, -4.85784, -6.59428, -8.2214, -9.53617, -10.6892, -11.6003, -12.0844, -12.0524, -11.9749, -11.4891, -10.6131, -9.49873, -8.1154, -6.41442, -5.09357, -3.99165, -2.72991, -1.71446, -0.56306, 0.440741]


#fourier frequency
ff = np.fft.fftfreq(len(t), (t[1]-t[0]))                               
#fourier amplitude
fa = abs(np.fft.fft(s, len(t)))                                            
#Position of maximum Amplitude 
pos_amax = np.argmax(fa[1:])+1                                        
#Frequency at maximum Amplitude (w/2pi)
ff_max = abs(ff[pos_amax])                                             
ff_guess = ff_max   
T_guess = 1000/ff_max 
#A_guess = np.std(s) *2. **0.5                                          
A_guess = min(s)
#c_guess = np.mean(s)
c_guess = 0
#First Guess for all paramters                           
f_guess = np.array([A_guess, 2*np.pi*ff_guess, 0., c_guess])  
#Sinus_Curve
def sin_func(t, A, w, phi, c):
return A * np.sin(w*t + phi) + c    
#Defining Bounds for A and c
c_bound = 0.1
A_bound = min(s)
#Bounds Array for curve_fit
param_bounds=([1.01*A_bound, -np.inf, -np.inf, -1*c_bound],[0.99*A_bound, np.inf, np.inf, c_bound])
popt, pcov = curve_fit(sin_func, t, s, p0=f_guess, bounds=param_bounds, maxfev=10000000)
#popt, pcov = curve_fit(sin_func, t, s, p0=f_guess, maxfev=10000000)
#
A, w, phi, c = popt 
f = w/(2.*np.pi)
T = 1000/f
t = np.array(t)  
s = np.array(s)
plt.figure(1)
#Generate Sinus Function
s_fit = A * np.sin(w*t + phi) + c
#Plotting
rcParams['figure.figsize'] =10, 5
fig, ax = plt.subplots()  
plt.plot(t, s, "b", label="Original")     
plt.plot(t, s_fit, 'k--', label="Fitting")  
ytitle='ytitle'
xtitle='xtitle'
ax.set(xlabel=xtitle, ylabel=ytitle)
ax.xaxis.set_major_formatter(FormatStrFormatter('%.2f'))               
ax.yaxis.set_major_formatter(FormatStrFormatter('%.2f'))  
ax.grid()
#Sidetext
ausgabe = ("Sinus-Fit nAmplitude = {:.2f} m/s^2 nPeriode = {:.2f} ms nOffset = {:.2f} m/s^2".format(A, abs(T), c))
plt.text(0.795, 0.7, ausgabe, family="sans-serif", fontsize=10, ha='left', va='top', transform=fig.transFigure)  
box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width * 0.85, box.height]) 
plt.show()

Answer 1

如果振幅是固定的，偏移量应该是零，为什么要把它们放在首位呢?此外，不需要fft来估计参数，因为有一个简单的线性方法。

看起来像这样:

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from scipy.integrate import cumtrapz

#Time
t = np.array([
313.544, 313.545, 313.546, 313.547, 313.548, 313.549, 313.55,
313.551, 313.552, 313.553, 313.554, 313.555, 313.556, 313.557,
313.558, 313.559, 313.56,  313.561, 313.562, 313.563, 313.564,
313.565, 313.566, 313.567
])
# ~t -= 313.544
#y-Values
s = np.array([
0.911188, -0.43135, -1.80997, -3.27816, -4.85784, -6.59428,
-8.2214, -9.53617, -10.6892, -11.6003, -12.0844, -12.0524,
-11.9749, -11.4891, -10.6131, -9.49873, -8.1154, -6.41442,
-5.09357, -3.99165, -2.72991, -1.71446, -0.56306, 0.440741
])
###fixing the amplitude the the maximum value of the data
amp = np.max( np.abs( s ) )
### sine function with the fixed amplitude from above and no offset
def sine(t, w, f ):
return amp * np.sin( w * t + f )

### getting the nonlinear parameters by using the fact that
### int int y = -y/w**2 + const1 * t + const2
### so the integro-equation as w hidden in a linear factor
### and we can use linear optimization to get it
Sy = cumtrapz( s, x=t, initial = 0 )    ### single integration
SSy = cumtrapz( Sy, x=t, initial = 0 )  ### double integratiom
VMXT = np.array( [ s, t, np.ones( len( t ) ) ] ) ### matrix describing the linear relationship
VMX = np.transpose( VMXT )
A = np.dot( VMXT, VMX )
SV = np.dot( VMXT, SSy )
AI = np.linalg.inv( A )
alpha = np.dot( AI , SV ) ### solution
wstart = np.sqrt( -1 / alpha[0] ) ### as mentioned above, the linear factor ist -1/w**2...so inverse function
### now getting phi by actiully fitting A sin + B cos
### and using tan phi = B/A
VMXT = np.array( [ np.sin( wstart * t ), np.cos( wstart * t ) ] )
VMX = np.transpose( VMXT )
A = np.dot( VMXT, VMX )
SV = np.dot( VMXT, s )
AI = np.linalg.inv( A )
alpha = np.dot( AI , SV ) ### solution
phistart = np.arctan2( alpha[1], alpha[0] )

op, cv = curve_fit( sine, t, s, p0=( wstart, phistart ) )
print( op )
yth = sine( t, *op )
fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.plot( t, s )
ax.plot( t, yth )
plt.show()

可以正常工作。