2010-11-19 4 views
2

У меня есть следующая функция, в которой я хочу интерполировать из таблицы с заданным значением. Фокус в том, что таблица определена в логарифмическом смысле так, что прямые линии между точками в log-log действительно экспоненциальны. Таким образом, я не могу использовать ни одну из типичных подпрограмм с интернетомными scipy.numpy magic to clean up function

Итак, вот что у меня есть:

PSD = np.array([[5.0, 0.001], 
       [25.0, 0.03], 
       [30.0, 0.03], 
       [89.0, 0.321], 
       [90.0, 1.0], 
       [260.0, 1.0], 
       [261.0, 0.03], 
       [359.0, 0.03], 
       [360.0, 0.5], 
       [520.0, 0.5], 
       [540.0, 0.25], 
       [780.0, 0.25], 
       [781.0, 0.03], 
       [2000.0, 0.03]]) 

def W_F(freq): 
    ''' 
    A line connecting two points in a log-log plot are exponential 
    ''' 
    w_f = [] 
    for f in freq: 
     index = np.searchsorted(PSD[:,0], f) 
     if index <= 0: 
      w_f.append(PSD[:,1][0]) 
     elif index + 1>= PSD.shape[0]: 
      w_f.append(PSD[:,1][-1]) 
     x0 = PSD[:,0][index-1] 
     F0 = PSD[:,1][index-1] 
     x1 = PSD[:,0][index] 
     F1 = PSD[:,1][index] 
     w_f.append(F0*(f/x0)**(math.log(F1/F0)/math.log(x1/x0))) 
    return np.array(w_f) 

Я ищу лучше, чище, «Numpy-иш» способ реализации этого

+1

'for freq = in f:' - это просто опечатка транскрипции или ваш код не работает? Кроме того, разве вы не должны использовать freq во второй строке для последней строки вместо f? –

+0

Да, функция неверна. Должно быть прочитано: – user90855

+0

Я исправил ошибку – user90855

ответ

3

Самый простой путь, чтобы просто взять логарифм PSD, а затем использовать SciPy интерполяции функции:

logPSD = numpy.log(PSD) 
logW_F = scipy.interpolate.interp1d(logPSD[:,0], logPSD[:,1]) 
W_F = numpy.exp(logW_F(numpy.log(f))) 

Это будет сгенерировано сообщение об ошибке для недоступный значений. Чтобы избежать этой ошибки, вы можете

  • Pass bounds_error=False функции interp1d() см the documentation.

  • Добавить запись в начале и конце PSD с очень маленьким и очень большим значением x для захвата всех возможных значений.

В качестве альтернативы использованию interp1d(), можно vectorise код, но я бы только это по причине.

+0

Только то, что я искал. Я думаю, что журнал (f) на последней строке должен быть изменен на numpy.log (f) – user90855

+0

Да, вы правы. Не тестировал код :) –