5

我想弄清楚我在这里不明白的是什么。

我正在关注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]我应该尝试不同的拟合方法吗?对哪方面的建议?

在此处输入图像描述

4

2 回答 2

2

下面是一些代码实现了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能够收敛到一个合理的答案。

于 2012-11-15T22:13:38.567 回答
2

From what I can see from playing a bit with leastsq (without fancy stuff from the cookbook, just plain direct calls to leastsq --- and by the way, full_output=True is your friend here), is that it's very hard to fit all three of the amplitude, frequency and phase in one go. On the other hand, if I fix the amplitude and fit the frequency and phase, it works; if I fix the frequency and fit the amplitude and phase, it works too.

There is more than one way out here. What might be the simplest one --- if you are sure you only have one sine wave (and this is easy to check with the Fourier transform), then you know the frequency from just the distance between consecutive maxima of your signal. Then fit the two remaining parameters.

If what you have is a mixture of several harmonics, well, again, Fourier transform will tell you that.

于 2012-11-15T21:19:19.410 回答