16

ここに画像の説明を入力

プログラムでスペクトル バンドを判断する方法を尋ねたところ、 @detlyは FWHM (半値全幅) を使用してピークの幅を決定することを提案しました。私は周りを検索し、モデルのフィッティングにFWHMを使用できることを発見しました（私は実際にはこれに素人です！）、特にGuassianモデル。具体的に2.354 * sigmaは、Guassian モデルの幅です。

悪いピークが存在するため、ガウス フィットを見ています。この写真には 4 つの整形式のピークがあります。次に、「ダブル」ピーク (両方が重要な場合があります) と 2 つの広がりのあるピークがあります。それらは、単純な FWHM に対する不可能な挑戦であることが証明されます。

x 座標のおおよその位置を考えると、Scip/Numpy のピークのガウス フィッティング (FWHM 計算の目的で) を生成するのを手伝ってもらえますか? Guassian が悪い選択である場合は、他のスキームを使用してください。

4

1 に答える 1

21

ガウス分布をフィッティングすることは良いアプローチです。また、初期パラメータ値を正しく推測できる場合は、それらすべてを一度に推測することができます。大きな問題はノイズです。実際には、各ピークを個別に適合させるか（つまり、特定のピークが一度に存在する範囲のみを考慮する）、データ全体のベースラインノイズ曲線を取得して減算する必要があります。

これは、複数のガウス分布に適合しようとするコードです。私は、8つの最も顕著なピークであると考えたものに対して、かなり緩いパラメーターをいくつか入力しました。さらに、バックグラウンドノイズをキャプチャするために、非常に広いグアシアンを1つ入力しました。処理する前に、投稿した画像を少しクリーンアップして、データを取得できるようにしました（マウスカーソルと軸の端を削除し、画像を反転しました）。

ここに画像の説明を入力してください

コード：

import Image
from scipy import *
from scipy.optimize import leastsq
import numpy as np

im = Image.open("proc.png")
im = im.convert('L')
h, w = im.size
#Extract data from the processed image:
im = np.asarray(im)
y_vals, junk  = np.mgrid[w:0:-1, h:0:-1]
y_vals[im < 255] = 0
y_vals = np.amax(y_vals,axis=0)

def gaussian(x, A, x0, sig):
    return A*exp(-(x-x0)**2/(2.0*sig**2))

def fit(p,x):
    return np.sum([gaussian(x, p[i*3],p[i*3+1],p[i*3+2]) 
                   for i in xrange(len(p)/3)],axis=0)

err = lambda p, x, y: fit(p,x)-y

#params are our intitial guesses for fitting gaussians, 
#(Amplitude, x value, sigma):
params = [[50,40,5],[50,110,5],[100,160,5],[100,220,5],
          [50,250,5],[100,260,5],[100,320,5], [100,400,5],   
          [30,300,150]]  # this last one is our noise estimate
params = np.asarray(params).flatten()

x  = xrange(len(y_vals))
results, value = leastsq(err, params, args=(x, y_vals))

for res in results.reshape(-1,3):
    print "amplitude, position, sigma", res

import pylab
pylab.subplot(211, title='original data')
pylab.plot(y_vals)
pylab.subplot(212, title='guassians fit')
y = fit(results, x)
pylab.plot(x, y)
pylab.savefig('fig2.png')
pylab.show()

これらは、フィッティングされた出力グアシアンパラメータです。

#Output:
amplitude, position, sigma [ 23.47418352  41.27086097   5.91012897]
amplitude, position, sigma [  16.26370263  106.39664543    3.67827824]
amplitude, position, sigma [  59.74500239  163.11210316    2.89866545]
amplitude, position, sigma [  57.91752456  220.24444939    2.91145375]
amplitude, position, sigma [  39.78579841  251.25955921    2.74646334]
amplitude, position, sigma [  86.50061756  260.62004831    3.08367483]
amplitude, position, sigma [  79.26849867  319.64343319    3.57224402]
amplitude, position, sigma [ 127.97009966  399.27833126    3.14331212]
amplitude, position, sigma [  20.21224956  379.41066063  195.47122312]
于 2012-06-04T16:09:36.127 に答える