2016-04-25 2 views
2

В настоящее время я интерполирование измерения с помощью кубических сплайнов, как показан на рисунке:Более вещий способ найти кривое пересечение

Interpolation of a normalised harmonic signal using cubic splines.

Идея заключается в том, я хочу найти всю ширину полувысота этой интерполяции , Для этого я использую небольшой фрагмент кода

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

, который возвращает мне отсортированные индексы пересечения с линией y=0.5. Тем не менее, я хочу решения на левом краю и правый край кривой, а иногда он возвращает мне две точки подряд. Есть ли элегантный питонический способ избежать этого? По крайней мере, намного лучше, чем мое хакерское решение:

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 

Благодарим за ответы!

ответ

0

Предполагая hhg_interp сортируется, и есть только один максимум, и хотел бы сделать gridsearch (т.е. рассчитать функцию в дискретных точках и работать с этими значениями), я хотел бы сделать следующее:

# 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 в основном представляет собой np.array переменной оси x. В этом случае вы можете взять его как hhg_interp = linspace (18.5, 19.5, 10000), чтобы иметь мелкую сетку для построения f против. – Roland

+0

Добавлен код для вычисления индексов левого и правого краев, а также несколько комментариев. Помогает ли это разъяснить? – cobaltfiftysix

+0

А, ок, я вижу. Спасибо, он работает красиво. Возможно, я обдумаю эту идею, чтобы лучше адаптировать ее к моим интересам. – Roland

Смежные вопросы