回答此问题可获得 20 贡献值,回答如果被采纳可获得 50 分。
<p>我有一个<code>space-separated csv</code>文件,其中包含一个度量值。第一列为测量时间,第二列为相应的测量值,第三列为误差。<a href="https://raw.githubusercontent.com/zabop/SO/master/signal.data" rel="nofollow noreferrer">The file can be found here.</a>我想使用Python将函数g的参数<code>a_i, f, phi_n</code>调整到数据中:</p>
<p><a href="https://i.stack.imgur.com/UgYVq.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/UgYVq.png" alt="enter image description here"/></a></p>
<p>在以下位置读取数据:</p>
<pre><code>import numpy as np
data=np.genfromtxt('signal.data')
time=data[:,0]
signal=data[:,1]
signalerror=data[:,2]
</code></pre>
<p>绘制数据:</p>
^{pr2}$
<p>得到结果:</p>
<p><a href="https://i.stack.imgur.com/OhzPb.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/OhzPb.png" alt="enter image description here"/></a></p>
<p>现在让我们计算周期信号的初步频率猜测:</p>
<pre><code>from gatspy.periodic import LombScargleFast
dmag=0.000005
nyquist_factor=40
model = LombScargleFast().fit(time, signal, dmag)
periods, power = model.periodogram_auto(nyquist_factor)
model.optimizer.period_range=(0.2, 10)
period = model.best_period
</code></pre>
<p>我们得到的结果是:0.5467448186001437</p>
<p>我为<code>N=10</code>定义了适合的函数:</p>
<pre><code>def G(x, A_0,
A_1, phi_1,
A_2, phi_2,
A_3, phi_3,
A_4, phi_4,
A_5, phi_5,
A_6, phi_6,
A_7, phi_7,
A_8, phi_8,
A_9, phi_9,
A_10, phi_10,
freq):
return (A_0 + A_1 * np.sin(2 * np.pi * 1 * freq * x + phi_1) +
A_2 * np.sin(2 * np.pi * 2 * freq * x + phi_2) +
A_3 * np.sin(2 * np.pi * 3 * freq * x + phi_3) +
A_4 * np.sin(2 * np.pi * 4 * freq * x + phi_4) +
A_5 * np.sin(2 * np.pi * 5 * freq * x + phi_5) +
A_6 * np.sin(2 * np.pi * 6 * freq * x + phi_6) +
A_7 * np.sin(2 * np.pi * 7 * freq * x + phi_7) +
A_8 * np.sin(2 * np.pi * 8 * freq * x + phi_8) +
A_9 * np.sin(2 * np.pi * 9 * freq * x + phi_9) +
A_10 * np.sin(2 * np.pi * 10 * freq * x + phi_10))
</code></pre>
<p>现在我们需要一个适合<code>G</code>的函数:</p>
<pre><code>def fitter(time, signal, signalerror, LSPfreq):
from scipy import optimize
pfit, pcov = optimize.curve_fit(lambda x, _A_0,
_A_1, _phi_1,
_A_2, _phi_2,
_A_3, _phi_3,
_A_4, _phi_4,
_A_5, _phi_5,
_A_6, _phi_6,
_A_7, _phi_7,
_A_8, _phi_8,
_A_9, _phi_9,
_A_10, _phi_10,
_freqfit:
G(x, _A_0, _A_1, _phi_1,
_A_2, _phi_2,
_A_3, _phi_3,
_A_4, _phi_4,
_A_5, _phi_5,
_A_6, _phi_6,
_A_7, _phi_7,
_A_8, _phi_8,
_A_9, _phi_9,
_A_10, _phi_10,
_freqfit),
time, signal, p0=[11, 2, 0, #p0 is the initial guess for numerical fitting
1, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
LSPfreq],
sigma=signalerror, absolute_sigma=True)
error = [] # DEFINE LIST TO CALC ERROR
for i in range(len(pfit)):
try:
error.<a href="https://www.cnpython.com/list/append" class="inner-link">append</a>(np.absolute(pcov[i][i]) ** 0.5) # CALCULATE SQUARE ROOT OF TRACE OF COVARIANCE MATRIX
except:
error.append(0.00)
perr_curvefit = np.array(error)
return pfit, perr_curvefit
</code></pre>
<p>看看我们得到了什么:</p>
<pre><code>LSPfreq=1/period
pfit, perr_curvefit = fitter(time, signal, signalerror, LSPfreq)
plt.figure()
model=G(time,pfit[0],pfit[1],pfit[2],pfit[3],pfit[4],pfit[5],pfit[6],pfit[7],pfit[8],pfit[8],pfit[10],pfit[11],pfit[12],pfit[13],pfit[14],pfit[15],pfit[16],pfit[17],pfit[18],pfit[19],pfit[20],pfit[21])
plt.scatter(time,model,marker='+')
plt.plot(time,signal,c='r')
plt.show()
</code></pre>
<p>屈服:</p>
<p><a href="https://i.stack.imgur.com/La2c3.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/La2c3.png" alt="enter image description here"/></a></p>
<p>这显然是错误的。如果我使用函数<code>fitter</code>定义中的初始猜测<code>p0</code>,我可以得到稍微好一点的结果。设置</p>
<pre><code>p0=[11, 1, 0,
0.1, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
LSPfreq]
</code></pre>
<p>给我们(放大):</p>
<p><a href="https://i.stack.imgur.com/iPQUG.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/iPQUG.png" alt="enter image description here"/></a></p>
<p>哪个好一点。高频分量仍然存在,尽管高频分量的振幅被猜测为零。原始的<code>p0</code>似乎也比基于对数据的视觉检查的修改版本更合理。在</p>
<p>我对<code>p0</code>使用了不同的值,虽然改变<code>p0</code>确实会改变结果,但我没有得到一条与数据相当吻合的行。在</p>
<p><strong>为什么这种模型拟合方法失败?我怎样才能更好的适应?</strong></p>
<p><a href="https://github.com/zabop/DSP1/blob/master/DSPcosinefitting.py" rel="nofollow noreferrer">The whole code can be found here.</a></p>
<hr/>
<p>这个问题的原始版本已经发布了<a href="https://dsp.stackexchange.com/questions/55240/how-can-i-improve-my-fit-of-cosines-to-periodic-data-using-python">here</a>。在</p>
<hr/>
<p><strong>编辑:</strong></p>
<p>PyCharm对代码的<code>p0</code>部分发出警告:</p>
<p><a href="https://i.stack.imgur.com/diibf.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/diibf.png" alt="enter image description here"/></a></p>
<blockquote>
<p>Expected type 'Union[None,int,float,complex]', got 'List[Union[int,float],Any]]' instead</p>
</blockquote>
<p>我不知道怎么处理,但可能有关系。在</p>