2012-05-26 4 views
7

enter image description hereこれらのスペクトル帯域は、目で判断するのに使用されました。

演算子位置と各ピークのを知って、スペクトルを調べるために使用され、スペクトルが属するピースを判断します。新しい方法では、画像がカメラによってスクリーンにキャプチャされます。そして、各バンドの幅はプログラム的に計算されなければならない。

旧システム:分光器 - >人間の目 新システム:分光器 - >カメラ - >プログラム

に良い方法は各バンドの幅を計算している何が、彼らのおおよそのX軸位置を与えられました。このタスクは目で完全に実行されていたため、プログラムで実行する必要があります。

申し訳ありませんが詳細が不足していますが、希少です。


前のグラフを生成したプログラムリスト。私はそれが関連している願っています:おおよその出発点を考えると

import Image 
from scipy import * 
from scipy.optimize import leastsq 

# Load the picture with PIL, process if needed 
pic   = asarray(Image.open("spectrum.jpg")) 

# Average the pixel values along vertical axis 
pic_avg  = pic.mean(axis=2) 
projection = pic_avg.sum(axis=0) 

# Set the min value to zero for a nice fit 
projection /= projection.mean() 
projection -= projection.min() 

#print projection 

# Fit function, two gaussians, adjust as needed 
def fitfunc(p,x): 
    return p[0]*exp(-(x-p[1])**2/(2.0*p[2]**2)) + \ 
     p[3]*exp(-(x-p[4])**2/(2.0*p[5]**2)) 
errfunc = lambda p, x, y: fitfunc(p,x)-y 

# Use scipy to fit, p0 is inital guess 
p0 = array([0,20,1,0,75,10]) 
X = xrange(len(projection)) 
p1, success = leastsq(errfunc, p0, args=(X,projection)) 
Y = fitfunc(p1,X) 

# Output the result 
print "Mean values at: ", p1[1], p1[4] 

# Plot the result 
from pylab import * 
#subplot(211) 
#imshow(pic) 
#subplot(223) 
#plot(projection) 
#subplot(224) 
#plot(X,Y,'r',lw=5) 
#show() 

subplot(311) 
imshow(pic) 
subplot(312) 
plot(projection) 
subplot(313) 
plot(X,Y,'r',lw=5) 
show() 
+5

人間の目が背景と色を区別できないしきい値を見つけます。それはかなり一定でなければなりません。次に、しきい値を上回るものが人間の目で「見られる」ようにデータポイントをしきい値化し、それらのデータポイントをクラスタ化して各クラスタの幅を見つけます。 – Blender

+0

山と坂を簡単に見つけるためにカーブの派生を取るかもしれませんか?バンドの幅はかなり一定ではありませんか?私がそれを理解するにつれて変化する振幅と位置です。 –

+0

@Blender、あなたはあなたの要点を詳述できますか? – aitchnyu

答えて

3

、あなたはこの点に最も近い極大値を見つける単純なアルゴリズムを使用することができます。あなたのフィッティングコードはすでにそれを行っているかもしれません(私があなたがそれをうまく使っているかどうかはわかりませんでした)。

ここでは、ユーザから与えられた出発点から見つける簡単なピークを示しているいくつかのコードです:

#!/usr/bin/env python 
from __future__ import division 
import numpy as np 
from matplotlib import pyplot as plt 

# Sample data with two peaks: small one at t=0.4, large one at t=0.8 
ts = np.arange(0, 1, 0.01) 
xs = np.exp(-((ts-0.4)/0.1)**2) + 2*np.exp(-((ts-0.8)/0.1)**2) 

# Say we have an approximate starting point of 0.35 
start_point = 0.35 

# Nearest index in "ts" to this starting point is... 
start_index = np.argmin(np.abs(ts - start_point)) 

# Find the local maxima in our data by looking for a sign change in 
# the first difference 
# From http://stackoverflow.com/a/9667121/188535 
maxes = (np.diff(np.sign(np.diff(xs))) < 0).nonzero()[0] + 1 

# Find which of these peaks is closest to our starting point 
index_of_peak = maxes[np.argmin(np.abs(maxes - start_index))] 

print "Peak centre at: %.3f" % ts[index_of_peak] 

# Quick plot showing the results: blue line is data, green dot is 
# starting point, red dot is peak location 
plt.plot(ts, xs, '-b') 
plt.plot(ts[start_index], xs[start_index], 'og') 
plt.plot(ts[index_of_peak], xs[index_of_peak], 'or') 
plt.show() 

上昇アップピークはあなたの出発点から完全に平滑である場合、このメソッドはのみ機能しますが。これが雑音に対してより復元力が必要な場合は、私はそれを使用していませんが、PyDSToolが役に立ったようです。このSciPy postは、ノイズの多いデータセットの1Dピークを検出するためにそれを使用する方法を詳しく説明しています。

この時点で、ピークの中心を見つけたとします。幅については、いくつかの方法がありますが、最も簡単なのはおそらく「半値全幅」(FWHM)です。繰り返しますが、これは簡単で脆弱です。近くのダブルピーク、またはノイズの多いデータの場合には破損します。

FWHMは、その名前が示す通りです。ピークの幅が最大値の半分であることがわかります。なお、ここではないいくつかのコードは、(それだけで上から上続けて)です:

# FWHM... 
half_max = xs[index_of_peak]/2 

# This finds where in the data we cross over the halfway point to our peak. Note 
# that this is global, so we need an extra step to refine these results to find 
# the closest crossovers to our peak. 

# Same sign-change-in-first-diff technique as above 
hm_left_indices = (np.diff(np.sign(np.diff(np.abs(xs[:index_of_peak] - half_max)))) > 0).nonzero()[0] + 1 
# Add "index_of_peak" to result because we cut off the left side of the data! 
hm_right_indices = (np.diff(np.sign(np.diff(np.abs(xs[index_of_peak:] - half_max)))) > 0).nonzero()[0] + 1 + index_of_peak 

# Find closest half-max index to peak 
hm_left_index = hm_left_indices[np.argmin(np.abs(hm_left_indices - index_of_peak))] 
hm_right_index = hm_right_indices[np.argmin(np.abs(hm_right_indices - index_of_peak))] 

# And the width is...  
fwhm = ts[hm_right_index] - ts[hm_left_index] 

print "Width: %.3f" % fwhm 

# Plot to illustrate FWHM: blue line is data, red circle is peak, red line 
# shows FWHM 
plt.plot(ts, xs, '-b') 
plt.plot(ts[index_of_peak], xs[index_of_peak], 'or') 
plt.plot(
    [ts[hm_left_index], ts[hm_right_index]], 
    [xs[hm_left_index], xs[hm_right_index]], '-r') 
plt.show() 

それは、半分値幅である必要はありません - 1つのコメンターが指摘するように、あなたは把握しようとすることができますピーク検出のためのオペレータの通常のしきい値がどこにあるかを決定し、プロセスのこのステップのアルゴリズムに変換します。

ピークを中心とするデータのサブセット、たとえば一方のローカルミニマムから他方のローカルミニマムにガウス曲線(または独自のモデル)をフィットさせるのがより堅牢な方法かもしれません。その曲線のパラメータの1つ(例えば、シグマ)を使用して幅を計算します。

これはたくさんのコードだとわかりましたが、私は意図的に "私の作業を表示する"ためのインデックス検索機能を考慮していませんでした。もちろんプロット関数があります。

これは、少なくとも、あなたの特定のセットに適したものを思い付くための出発点となることを望みます。

+0

「ベース」がしきい値(この場合はピークの0.5)よりも低い場合に例外をスローする変更を含めることができますか? – aitchnyu

+0

@aitchnyuここで私は少し単純化しました:私はちょうど全高の0.5と言ったが、FWHMはベース*(またはノイズフロア)より上のピーク*の高さの半分として測定されています。したがって、ベースがピークの0.5よりも低いことは意味をなさない。 – detly

+0

はい、アルゴリズムが異常なピークに達した場合に例外をスローするようにします。最大ピーク値0.5まで単調に低下するピークはありません(グラフの二重項を参照)。そして、私は**幅だけ**が必要です、振幅は、機器の特性のために変更することができます。 「ノイズフロア」は理想的には無音で、黒の強さに対応していなければなりません。 – aitchnyu

0

最良の方法は、多数のメソッドを人間の結果と統計的に比較することです。

多種多様なデータとさまざまな測定値(さまざまなしきい値、さまざまなしきい値を上回るエリア、さまざまなしきい値の選択方法、2番目の瞬間、さまざまな程度の多項式曲線適合、パターンマッチングなど)を取ります。これらの推定値を同じデータセットの人間測定値と比較する。熟練した人間の結果と最も相関する推定方法を選ぶ。それともなど

2

後期パーティーに、将来的にこの質問を越えてくる人のために...他のピークから様々な分離のために、様々な高さのそれぞれに最適なものをいくつかの方法を選択し、

目の動きのデータは、これと非常によく似ています。私はNystrom + Holmqvist, 2010で使用されているアプローチに基づいています。 Savitsky-Golayフィルタ(scipy v0.14 +のscipy.signal.savgol_filter)を使用してデータをスムージングして、大きなピークを損なわずに低レベルのノイズを除去してください。オーダー2とウィンドウサイズ約を使用することをお勧めします検出できるようにしたい最小ピークの幅の2倍です。特定のy値より上のすべての値を任意に削除することで、バンドの位置を知ることができます(numpy.nanに設定します)。次に残りの(平均)標準偏差と(標準偏差)標準偏差を取って、平均+ [パラメータ] *標準より大きな値をすべて削除します(私は彼らが紙で6を使うと思います)。データポイントを削除しない限り繰り返しますが、データによっては[パラメータ]の値が安定しないことがあります。次に、numpy.isnan()を使用してイベントと非イベントを見つけ、numpy.diff()を使用して、それぞれのイベントの開始と終了を見つけます(値はそれぞれ-1と1)。さらに正確な開始点と終了点を取得するには、各開始点から後方に向かってデータをスキャンし、各端点から前方に移動して、平均+ [別のパラメータ] * stdより小さい値を持つ最も近いローカル最小点を見つけ出すことができます紙の中で)。次に、開始点と終了点の間のデータポイントを数えるだけです。

これは、ダブルピークでは機能しません。あなたはそれについていくつか外挿をしなければならないでしょう。

関連する問題