python - 非线性拟合正弦曲线失败

给定一些生成的二维玩具数据，我一直在尝试拟合正弦曲线的幅度、频率和相位。 (代码在最后)

为了获得三个参数的估计值，我首先执行 FFT。我使用来自 FFT 的值作为实际频率和相位的初始猜测，然后适合它们(逐行)。我写了我的代码，这样我就可以输入我想要频率所在的 FFT 的哪个 bin，这样我就可以检查拟合是否正常工作。但是有一些非常奇怪的行为。如果我的输入 bin 是 3.1(一个非整数 bin，所以 FFT 不会给我正确的频率)，那么拟合效果很好。但是如果输入 bin 是 3(所以 FFT 输出确切的频率)，那么我的拟合失败，我试图理解为什么。

这是我分别将输入 bin(在 X 和 Y 方向)设置为 3.0 和 2.1 时的输出:

(右图为数据-拟合)

这是我将输入 bin 设为 3.0 和 2.0 时的输出:

问题: 为什么当我输入曲线的确切频率时非线性拟合失败？

代码:

#! /usr/bin/python

# For the purposes of this code, it's easier to think of the X-Y axes as transposed,
# so the X axis is vertical and the Y axis is horizontal

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as optimize
import itertools
import sys

PI = np.pi

# Function which accepts paramters to define a sin curve
# Used for the non linear fit
def sineFit(t, a, f, p):
   return a * np.sin(2.0 * PI * f*t + p)

xSize    = 18
ySize    = 60
npt      = xSize * ySize

# Get frequency bin from user input
xFreq    = float(sys.argv[1])
yFreq    = float(sys.argv[2])

xPeriod  = xSize/xFreq
yPeriod  = ySize/yFreq

# arrays should be defined here

# Generate the 2D sine curve
for jj in range (0, xSize):
   for ii in range(0, ySize):
      sineGen[jj, ii] = np.cos(2.0*PI*(ii/xPeriod + jj/yPeriod))

# Compute 2dim FFT as well as freq bins along each axis
fftData  = np.fft.fft2(sineGen)
fftMean  = np.mean(fftData)
fftRMS   = np.std(fftData)
xFreqArr = np.fft.fftfreq(fftData.shape[1]) # Frequency bins along x
yFreqArr = np.fft.fftfreq(fftData.shape[0]) # Frequency bins along y

# Find peak of FFT, and position of peak
maxVal = np.amax(np.abs(fftData))
maxPos = np.where(np.abs(fftData) == maxVal)

# Iterate through peaks in the FFT
# For this example, number of loops will always be only one

prevPhase = -1000
for col, row in itertools.izip(maxPos[0], maxPos[1]):

   # Initial guesses for fit parameters from FFT
   init_phase  = np.angle(fftData[col,row])
   init_amp    = 2.0 * maxVal/npt
   init_freqY  = yFreqArr[col]
   init_freqX  = xFreqArr[row]

   cntr  = 0
   if prevPhase == -1000:
      prevPhase = init_phase

   guess = [init_amp, init_freqX, prevPhase]
   # Fit each row of the 2D sine curve independently
   for rr in sineGen:
      (amp, freq, phs), pcov = optimize.curve_fit(sineFit, xDat, rr, guess)
      # xDat is an linspace array, containing a list of numbers from 0 to xSize-1

      # Subtract fit from original data and plot
      fitData     = sineFit(xDat, amp, freq, phs)
      sub1        = rr - fitData

      # Plot
      fig1 = plt.figure()
      ax1  = fig1.add_subplot(121)
      p1,  = ax1.plot(rr, 'g')
      p2,  = ax1.plot(fitData, 'b')
      plt.legend([p1,p2], ["data", "fit"])

      ax2  = fig1.add_subplot(122)
      p3,  = ax2.plot(sub1)
      plt.legend([p3], ['residual1'])

      fig1.tight_layout()

      plt.show()
      cntr += 1
      prevPhase = phs # Update guess for phase of sine curve

最佳答案

我试图将您问题的重要部分提炼成这个答案。

首先，尝试拟合单个数据块，而不是数组。一旦你确信你的模型足够了，你就可以继续前进。

你的拟合只会和你的模型一样好，如果你转向不是“正弦”的东西——比如你需要相应地调整。

拟合是一门“艺术”，因为初始条件可以极大地改变 改变误差函数的收敛性。此外，您的拟合中可能有多个最小值，因此您经常不得不担心所提出的解决方案的独特性。

虽然您的 FFT 想法是正确的，但我认为您的实现并不完全正确。下面的代码应该是一个很棒的玩具系统。它生成 f(x) = a0*sin(a1*x+a2) 类型的随机数据。有时随机的初始猜测会起作用，有时它会失败。但是，使用频率的 FFT 猜测收敛应该始终适用于该系统。示例输出:

import numpy as np
import pylab as plt
import scipy.optimize as optimize

# This is your target function
def sineFit(t, (a, f, p)):
    return a * np.sin(2.0*np.pi*f*t + p)

# This is our "error" function
def err_func(p0, X, Y, target_function):
    err = ((Y - target_function(X, p0))**2).sum()
    return err


# Try out different parameters, sometimes the random guess works
# sometimes it fails. The FFT solution should always work for this problem
inital_args = np.random.random(3)

X = np.linspace(0, 10, 1000)
Y = sineFit(X, inital_args)

# Use a random inital guess
inital_guess = np.random.random(3)

# Fit
sol = optimize.fmin(err_func, inital_guess, args=(X,Y,sineFit))

# Plot the fit
Y2 = sineFit(X, sol)
plt.figure(figsize=(15,10))
plt.subplot(211)
plt.title("Random Inital Guess: Final Parameters: %s"%sol)
plt.plot(X,Y)
plt.plot(X,Y2,'r',alpha=.5,lw=10)

# Use an improved "fft" guess for the frequency
# this will be the max in k-space
timestep = X[1]-X[0]
guess_k = np.argmax( np.fft.rfft(Y) )
guess_f = np.fft.fftfreq(X.size, timestep)[guess_k]
inital_guess[1] = guess_f

# Guess the amplitiude by taking the max of the absolute values
inital_guess[0] = np.abs(Y).max()

sol = optimize.fmin(err_func, inital_guess, args=(X,Y,sineFit))
Y2 = sineFit(X, sol)

plt.subplot(212)
plt.title("FFT Guess          : Final Parameters: %s"%sol)
plt.plot(X,Y)
plt.plot(X,Y2,'r',alpha=.5,lw=10)
plt.show()

关于python - 非线性拟合正弦曲线失败，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/18059819/