2016-04-25 19 views
2
写真のように、私は現在、次スプラインを使用して測定を補間しています

曲線の交点を見つけるために、より神託の方法

Interpolation of a normalised harmonic signal using cubic splines.

アイデア私はその補間の半値全幅を見つけたいです。そのために私は私のラインy=0.5との交点のソートされたインデックスを返すコード

f = interpolate.interp1d(hhg_data, signal_array, kind='cubic') 
idx2 = np.argsort(np.abs(f(hhg_interp)-0.5)) 

の小片を使用しています。しかし、私は曲線の左端と右端に解を欲しく、時にはそれが私に連続する2点を返します。これを避けるためのエレガントな無邪気な方法はありますか?少なくとも私のハックなソリューションよりもはるかに優れています。

idx_sorted = [] 
counter = 0 
counter2 = 0 
while(counter <= 1): 
    if idx2[counter2] != idx2[counter2-1]+1 and idx2[counter2] != idx2[counter2-1]-1: 
     idx_sorted.append(idx2[counter2]) 
     counter+=1 
    counter2+=1 

ありがとうございました! hhg_interp(すなわち離散点であなたの関数を計算し、それらの値で動作)、ソート、唯一の最大値があるとgridsearchをしたいと思っている、私は次のことを行うだろうと仮定

答えて

0

# hhg_interp could be something like 
# hhg_interp = linspace(18.5, 19.5, 10000) 
y = f(hhg_interp) # hhg_interp is an array filled with finely 
        # spaced x-axis values 
# get the indices for all y larger than the half maximum (0.5 in this case) 
indices_above_fwhm = np.where(y>0.5)[0] 
# if you're interested in the indices 
# the first element of the "indices larger then the half maximum" array 
# corresponds to your left edge 
index_left_edge = indices_above_fwhm[0] 
# the last element of the "indices larger then the half maximum array" 
# corresponds to your right edge 
index_right_edge = indices_above_fwhm[-1] 


# then you can get the corresponding x-axis values for these indices: 
x_left = hhg_interp[index_left_edge] 
x_right = hhg_interp[index_right_edge] 

# or directly get the x-axis values of the left and right edges 
cond = y > 0.5 # select all points that are above your half maximum 
x_left, x_right = hhg_interp[cond][[0,-1]] # then get the first 
           # and last of the points above maximum 

ですそれはあなたが探している情報ですか?

+0

本当に分かりませんが、hhg_interpは基本的にx軸変数のnp.arrayです。その場合、hhg_interp = linspace(18.5、19.5、10000)とすると、fをプロットする細かいグリッドが得られます。 – Roland

+0

左端と右端のインデックスだけでなく、いくつかのコメントを計算するためのコードをいくつか追加しました。これは明らかにするのに役立ちますか? – cobaltfiftysix

+0

ああ、私は参照してください。ありがとうございました。おそらく私はこのアイデアを取り上げて、自分の興味のあるケースに最も適したものにするでしょう。 – Roland

関連する問題