<p>下面是一些代码,实现了振亚的一些想法。
它使用</p>
<pre><code>yhat = fftpack.rfft(yReal)
idx = (yhat**2).argmax()
freqs = fftpack.rfftfreq(N, d = (xReal[1]-xReal[0])/(2*pi))
frequency = freqs[idx]
</code></pre>
<p>猜测数据的主频,以及</p>
<pre><code>amplitude = yReal.max()
</code></pre>
<p>去猜振幅。</p>
<hr/>
<pre><code>import numpy as np
import scipy.optimize as optimize
import scipy.fftpack as fftpack
import matplotlib.pyplot as plt
pi = np.pi
plt.figure(figsize = (15, 5))
# generate a perfect data set (my real data have tiny error)
def mysine(x, a1, a2, a3):
return a1 * np.sin(a2 * x + a3)
N = 5000
xmax = 10
xReal = np.linspace(0, xmax, N)
a1 = 200.
a2 = 2*pi/10.5 # omega, 10.5 is the period
a3 = np.deg2rad(10.) # 10 degree phase offset
print(a1, a2, a3)
yReal = mysine(xReal, a1, a2, a3) + 0.2*np.random.normal(size=len(xReal))
yhat = fftpack.rfft(yReal)
idx = (yhat**2).argmax()
freqs = fftpack.rfftfreq(N, d = (xReal[1]-xReal[0])/(2*pi))
frequency = freqs[idx]
amplitude = yReal.max()
guess = [amplitude, frequency, 0.]
print(guess)
(amplitude, frequency, phase), pcov = optimize.curve_fit(
mysine, xReal, yReal, guess)
period = 2*pi/frequency
print(amplitude, frequency, phase)
xx = xReal
yy = mysine(xx, amplitude, frequency, phase)
# plot the real data
plt.plot(xReal, yReal, 'r', label = 'Real Values')
plt.plot(xx, yy , label = 'fit')
plt.legend(shadow = True, fancybox = True)
plt.show()
</code></pre>
<p>收益率</p>
<pre><code>(200.0, 0.5983986006837702, 0.17453292519943295) # (a1, a2, a3)
[199.61981404516041, 0.61575216010359946, 0.0] # guess
(200.06145097308041, 0.59841420869261097, 0.17487141943703263) # fitted parameters
</code></pre>
<p>注意,通过使用fft,对频率的猜测已经非常接近最终的拟合参数。</p>
<p>似乎不需要修复任何参数。
通过使频率猜测更接近实际值,<code>optimize.curve_fit</code>能够收敛到一个合理的答案。</p>