Ошибка нелинейной подгонки к синусоидальной кривой

Question

Я старался соответствовать амплитуде, частоте и фазе кривой синуса, учитывая некоторые генерируемые двумерные данные игрушек. (Код в конце)

Чтобы получить оценки для трех параметров, я сначала выполняю БПФ. Я использую значения из БПФ в качестве исходных догадок для фактической частоты и фазы, а затем для них (строки за строкой). Я написал свой код таким образом, что я ввожу, в каком бункере БПФ мне нужна частота, поэтому я могу проверить, работает ли фитинг. Но есть странное поведение. Если мой входной бит равен 3.1 (не встроенный бит, поэтому БПФ не даст мне правильную частоту), то пригонка отлично работает. Но если входной бит равен 3 (поэтому FFT выводит точную частоту), то моя подгонка не удалась, и я пытаюсь понять, почему.

Вот вывод, когда я дам входных бункеров (в направлении X и Y), как 3,0 и 2,1 соответственно:

(Участок справа есть данные - подходят)

Вот выход, когда я даю входные лотки как 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 

Answer 1

Я попытался перевести важные части вашего вопроса в этот ответ.

1. Прежде всего, попытайтесь подгонки одного блок данных, а не массив. Как только вы уверены, что ваша модель достаточна, вы можете двигаться дальше.
2. Ваша пригонка будет только такой же хорошей, как ваша модель, если вы перейдете к чему-то не «синусоидальному», как вам нужно будет соответствующим образом отрегулировать.
3. Установка является «искусством», поскольку начальные условия могут значительно изменить сходимость функции ошибки. Кроме того, в ваших припадках может быть более одного минимума, поэтому вам часто приходится беспокоиться о уникальности вашего предлагаемого решения.

Пока вы были на правильном пути с идеей 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() 

Answer 2

Без понимания большей части кода, в соответствии с http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html:

(amp2, freq2, phs2), pcov = optimize.curve_fit(sineFit, tDat, 
                sub1, guess2) 

должны стать:

(amp2, freq2, phs2), pcov = optimize.curve_fit(sineFit, tDat, 
                 sub1, p0=guess2) 

Предполагая, что tDat и sub1 являются x и y, что Я делаю трюк. Но, опять же, довольно сложно понять такой сложный код с таким количеством взаимосвязанных переменных и комментариев вообще. Код всегда должен строиться снизу вверх, а это значит, что вы не выполняете цикл пригонок, когда один из них не работает, вы не добавляете шум, пока код не работает, чтобы соответствовать нешумным примерам ... Хорошо удачи!

Answer 3

К «Ничего особенного» я ​​имел в виду что-то вроде удаления все, что не связано с подгонкой и делать упрощенный макет, например, такие как:

import numpy as np 
import scipy.optimize as optimize 

def sineFit(t, a, f, p): 
     return a * np.sin(2.0 * np.pi * f*t + p) 


# Create array of x and y with given parameters 
x = np.asarray(range(100)) 
y = sineFit(x, 1, 0.05, 0) 

# Give a guess and fit, printing result of the fitted values 
guess = [1., 0.05, 0.] 
print optimize.curve_fit(sineFit, x, y, guess)[0] 

Результатом этого является именно ответ:

[1. 0.05 0.] 

Но если изменить догадываюсь не слишком много, достаточно просто:

# Give a guess and fit, printing result of the fitted values 
guess = [1., 0.06, 0.] 
print optimize.curve_fit(sineFit, x, y, guess)[0] 

т он приводит к абсурдно неправильным цифрам:

[ 0.00823701 0.06391323 -1.20382787] 

Можете ли вы объяснить это поведение?

Answer 4

Проблема заключается в том, чтобы из-за плохое начальное приближение фазы, а не частоты. При циклическом перемещении по строкам genSine (внутренний цикл) вы используете результат соответствия предыдущей строки как начальное предположение для следующей строки, которая не работает всегда. Если вы определяете фазу из fft текущей строки и используете ее как исходную догадку, она будет успешной. Вы можете изменить внутренний цикл следующим образом:

for n,rr in enumerate(sineGen): 
    fftx = np.fft.fft(rr) 
    fftx = fftx[:len(fftx)/2] 
    idx = np.argmax(np.abs(fftx)) 
    init_phase = np.angle(fftx[idx]) 
    print fftx[idx], init_phase 
    ... 

Кроме того, необходимо изменить

def sineFit(t, a, f, p): 
    return a * np.sin(2.0 * np.pi * f*t + p) 

в

def sineFit(t, a, f, p): 
    return a * np.cos(2.0 * np.pi * f*t + p) 

поскольку фаза = 0 означает, что мнимая часть FFT является нуль и, следовательно, функция является косинусоидальной.

КПП. в вашем примере выше все еще отсутствуют определения sineGen и xDat.

Answer 5

Вы можете использовать curve_fit с серией тригонометрических функций, как правило, очень надежные и складываемые в точности, что вам нужно только за счет увеличения числа членов ... вот пример:

from scipy import sin, cos, linspace 
def f(x, a0,s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12, 
      c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12): 
    return a0 + s1*sin(1*x) + c1*cos(1*x) \ 
       + s2*sin(2*x) + c2*cos(2*x) \ 
       + s3*sin(3*x) + c3*cos(3*x) \ 
       + s4*sin(4*x) + c4*cos(4*x) \ 
       + s5*sin(5*x) + c5*cos(5*x) \ 
       + s6*sin(6*x) + c6*cos(6*x) \ 
       + s7*sin(7*x) + c7*cos(7*x) \ 
       + s8*sin(8*x) + c8*cos(8*x) \ 
       + s9*sin(9*x) + c9*cos(9*x) \ 
      + s10*sin(9*x) + c10*cos(9*x) \ 
      + s11*sin(9*x) + c11*cos(9*x) \ 
      + s12*sin(9*x) + c12*cos(9*x) 

from scipy.optimize import curve_fit 
pi/2./(x.max() - x.min()) 
x_norm *= norm_factor 
popt, pcov = curve_fit(f, x_norm, y) 
x_fit = linspace(x_norm.min(), x_norm.max(), 1000) 
y_fit = f(x_fit, *popt) 
plt.plot(x_fit/x_norm, y_fit)