SciPy最小二乘拟合正弦波失败

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]我应该尝试不同的拟合方法吗？有什么建议吗？

- Brian Larsen

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能够收敛到一个合理的答案。

- unutbu

---

更好的猜测在这种情况下似乎是一个非常聪明的事情。感谢您提出这个想法。 - Brian Larsen

Answer 2

从我对leastsq进行的一些简单操作中可以看出，同时拟合振幅、频率和相位非常困难。然而，如果我固定振幅并拟合频率和相位，它就有效；如果我固定频率并拟合振幅和相位，也是有效的。

这里有多种方法可选。可能最简单的方法是，如果您确定只有一个正弦波（傅里叶变换很容易检查），则只需通过信号连续极大值之间的距离即可确定频率。然后拟合剩余的两个参数。

如果您的数据是几个谐波的混合物，同样可以通过傅里叶变换来确定。

- ev-br

---

谢谢，我确定我只有一个正弦波，但振幅会随时间变化很大，而频率只会略微变化。我将尝试通过仅获取数据中的高值来修复振幅问题。 - Brian Larsen

我会修复频率，而不是振幅。或者将所有内容转换为傅里叶空间，并在频率域中进行所有拟合。 - ev-br

是的，这里的问题是太空船在进入日食时会改变频率，而我只关心频率和相位，而不关心幅度，所以这似乎更自然。 - Brian Larsen

如果频率在时间上不是恒定的，那么它就不是单一的正弦波，你无法用一个正弦波来拟合它。因此，在频域中工作看起来更有前途。 - ev-br

我会看一下。由于周期变化缓慢，因此我想我可以分段处理。我将研究频域方法，但我对此不太熟悉，所以一直试图避免使用它。 - Brian Larsen