4

生成された 2 次元のおもちゃのデータを使用して、正弦曲線の振幅、周波数、位相を当てはめようとしました。(最後にコード)

3 つのパラメーターの推定値を得るために、まず FFT を実行します。FFT からの値を実際の周波数と位相の初期推定として使用し、次にそれらを (行ごとに) 適合させます。周波数を入れたいFFTのビンを入力するようにコードを書いたので、フィッティングがうまく機能しているかどうかを確認できます。しかし、かなり奇妙な動作があります。入力ビンが 3.1 の場合 (非整数ビンであるため、FFT では適切な周波数が得られません)、フィットは素晴らしく機能します。しかし、入力ビンが 3 の場合 (FFT が正確な周波数を出力するため)、適合は失敗し、その理由を理解しようとしています。

入力ビン (X 方向と Y 方向) にそれぞれ 3.0 と 2.1 を指定したときの出力は次のとおりです。

(右側のプロットはデータ - 適合) 図1

入力ビンに 3.0 と 2.0 を指定した場合の出力は次のとおりです。 図2

質問:曲線の正確な度数を入力すると、非線形フィットが失敗するのはなぜですか?

コード：

#! /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
4

5 に答える 5

3

あなたの質問の重要な部分をこの回答にまとめようとしました。

  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()
于 2013-08-08T16:15:22.093 に答える
0

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 であると仮定すると、うまくいくはずです。しかし、繰り返しになりますが、非常に多くの相互リンクされた変数とコメントがまったくないこのような複雑なコードを理解することは非常に困難です。コードは常にボトムアップで構築する必要があります。つまり、1 つのフィットが機能しない場合は、フィットのループを実行せず、ノイズのない例に合わせてコードが機能するまでノイズを追加しません...良い幸運！

于 2013-08-05T14:46:49.517 に答える
0

「空想的なものは何もない」とは、フィットに関係のないものをすべて削除し、次のような単純化されたモックの例を実行するようなものを意味しました。

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]

この振る舞いを説明できますか？

于 2013-08-06T11:15:01.567 に答える