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是你的朋友在这里），是难以一次性适应全部三个幅度，频率和相位。另一方面，如果我确定幅度并适合频率和相位，它就可以工作;如果我确定频率并适合幅度和相位，它也可以工作。

这里有不止一条路。最简单的一个---如果你确信你只有一个正弦波（这很容易通过傅里叶变换来检查），那么你就可以知道信号连续最大值之间的距离。然后再拟合剩下的两个参数。

如果你有什么是几个谐波的混合，那么傅立叶变换会告诉你。