SciPy leastsq適合正弦波失敗

Question

我想弄清楚這是什麼，我不明白這裏。

我正在關注http://www.scipy.org/Cookbook/FittingData並嘗試適合正弦波。真正的問題是衛星磁強計數據在旋轉的航天器上產生很好的正弦波。我創建了一個數據集，然後試圖使其適合於恢復輸入。

這裏是我的代碼：

import numpy as np 
from scipy import optimize 

from scipy.optimize import curve_fit, leastsq 

import matplotlib.pyplot as plt 


class Parameter: 
    def __init__(self, value): 
      self.value = value 

    def set(self, value): 
      self.value = value 

    def __call__(self): 
      return self.value 

def fit(function, parameters, y, x = None): 
    def f(params): 
     i = 0 
     for p in parameters: 
      p.set(params[i]) 
      i += 1 
     return y - function(x) 

    if x is None: x = np.arange(y.shape[0]) 
    p = [param() for param in parameters] 
    return optimize.leastsq(f, p, full_output=True, ftol=1e-6, xtol=1e-6) 

# generate a perfect data set (my real data have tiny error) 
def mysine(x, a1, a2, a3): 
    return a1 * np.sin(a2 * x + a3) 

xReal = np.arange(500)/10. 
a1 = 200. 
a2 = 2*np.pi/10.5 # omega, 10.5 is the period 
a3 = np.deg2rad(10.) # 10 degree phase offset 
yReal = mysine(xReal, a1, a2, a3) 

# plot the real data 
plt.figure(figsize=(15,5)) 
plt.plot(xReal, yReal, 'r', label='Real Values') 

# giving initial parameters 
amplitude = Parameter(175.) 
frequency = Parameter(2*np.pi/8.) 
phase = Parameter(0.0) 

# define your function: 
def f(x): return amplitude() * np.sin(frequency() * x + phase()) 

# fit! (given that data is an array with the data to fit) 
out = fit(f, [amplitude, frequency, phase], yReal, xReal) 
period = 2*np.pi/frequency() 
print amplitude(), period, np.rad2deg(phase()) 

xx = np.linspace(0, np.max(xReal), 50) 
plt.plot(xx, f(xx) , label='fit') 
plt.legend(shadow=True, fancybox=True) 

這使得這個陰謀：

回收的[44.2434221897 8.094832581 -61.6204033699]擬合參數都沒有相似之處什麼我開始用。

對於我不理解或做錯的任何想法？

scipy.__version__ 
'0.10.1' 

---

編輯： 固定一個參數建議。在上例中，將幅度固定爲np.histogram(yReal)[1][-1]仍然會產生不可接受的輸出。適合：[175.0 8.31681375217 6.0]我應該嘗試一種不同的裝配方法嗎？有哪些建議？

Answer 1

下面是一些實現Zhenya想法的代碼。 它使用

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() 

猜測振幅。

---

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() 

產量

(200.0, 0.5983986006837702, 0.17453292519943295) # (a1, a2, a3) 
[199.61981404516041, 0.61575216010359946, 0.0]  # guess 
(200.06145097308041, 0.59841420869261097, 0.17487141943703263) # fitted parameters 

注意，通過使用FFT，對於頻率的猜測是已經非常接近最終擬合參數。

看來你不需要修復任何參數。 通過使頻率猜測更接近實際值，optimize.curve_fit能夠收斂到合理的答案。

Answer 2

從我可以從打了一下，用leastsq看（從菜譜沒有花哨的東西，只是簡單的直接調用leastsq ---順便說一句，full_output=True是你的朋友在這裏），是難以一次性適應全部三個幅度，頻率和相位。另一方面，如果我確定幅度並適合頻率和相位，它就可以工作;如果我確定頻率並適合幅度和相位，它也可以工作。

這裏有不止一條路。最簡單的一個---如果你確信你只有一個正弦波（這很容易通過傅里葉變換來檢查），那麼你就可以知道信號連續最大值之間的距離。然後再擬合剩下的兩個參數。

如果你有什麼是幾個諧波的混合，那麼傅立葉變換會告訴你。