пересчет, интерполирующая матрица

Question

Я пытаюсь интерполировать некоторые данные с целью построения. Например, учитывая N точек данных, я хотел бы иметь возможность генерировать «гладкий» график, состоящий из 10 * N или около того интерполированных точек данных.

Мой подход заключается в создании матрицы N × 10 * N и вычислении внутреннего произведения исходного вектора и сгенерированной мной матрицы с получением вектора 1 × 10 * N. Я уже разработал математику, которую хотел бы использовать для интерполяции, но мой код довольно медленный. Я довольно новичок в Python, поэтому я надеюсь, что некоторые из присутствующих здесь экспертов могут дать мне несколько идей о том, как я могу попытаться ускорить мой код.

Я думаю, что часть проблемы заключается в том, что для генерации матрицы требуется 10 * N ^ 2 вызовов следующей функции:

def sinc(x):
    import math
    try:
        return math.sin(math.pi * x) / (math.pi * x)
    except ZeroDivisionError:
        return 1.0

(Это взято из теории дискретизации . По сути, я пытаюсь воссоздать сигнал из его выборок и повысить частоту до более высокой частоты.)

Матрица генерируется следующим образом:

def resampleMatrix(Tso, Tsf, o, f):
    from numpy import array as npar
    retval = []

    for i in range(f):
        retval.append([sinc((Tsf*i - Tso*j)/Tso) for j in range(o)])

    return npar(retval)

Я собираюсь разбить задачу на более мелкие части, потому что мне не нравится идея матрицы N ^ 2, сидящей в памяти. Я мог бы, вероятно, превратить resampleMatrix в функцию генератора и выполнять внутренний продукт построчно, но я не думаю, что это сильно ускорит мой код, пока я не начну разбирать содержимое в памяти и из памяти.

Заранее спасибо за ваши предложения!

Answer 1

Это повышение частоты дискретизации. См. Справка по повторной дискретизации / повышающей дискретизации для некоторых примеров решений.

Быстрый способ сделать это (для автономных данных, таких как ваше графическое приложение) - это использовать БПФ. Это то, что делает нативная функция SciPy resample() . Тем не менее, он принимает периодический сигнал , поэтому он не совсем такой же . См. эту ссылку :

Вот второй вопрос, касающийся интерполяции реального сигнала во временной области, и это действительно большая проблема. Этот точный алгоритм интерполяции дает правильные результаты только в том случае, если исходная последовательность x (n) является периодической в ​​течение полного интервала времени.

Ваша функция предполагает, что все выборки сигнала находятся за пределами заданного диапазона, поэтому оба метода будут отклоняться от центральной точки. Если сначала заполнить сигнал большим количеством нулей, это даст очень близкий результат. Еще несколько нулей за краем участка, не показанного здесь:

Кубическая интерполяция не будет корректной для целей передискретизации. Этот пример является крайним случаем (около частоты дискретизации), но, как вы можете видеть, кубическая интерполяция даже не близка. Для низких частот это должно быть довольно точно.

Answer 2

Если вы хотите интерполировать данные достаточно общим и быстрым способом, сплайны или полиномы очень полезны. У Scipy есть модуль scipy.interpolate, который очень полезен. Вы можете найти много примеров на официальных страницах.

Answer 3

Если ваш единственный интерес состоит в том, чтобы «сгенерировать« гладкий »график», я бы просто использовал простую полиномиальную кривую сплайна:

Для любых двух смежных точек данных коэффициенты полиномиальной функции третьей степени могут быть вычислены из координат этих точек данных и двух дополнительных точек слева и справа (без учета граничных точек). Это создаст точки на хорошем гладкая кривая с непрерывным первым диривативом. Есть прямая формула для преобразования 4 координат в 4 полиномиальных коэффициента, но я не хочу лишать вас удовольствия от ее поиска; o).

Answer 4

Я не совсем уверен, что вы пытаетесь сделать, но есть некоторые ускорения, которые вы можете сделать, чтобы создать матрицу. Предложение Braincore использовать numpy.sinc - это первый шаг, но второй заключается в том, чтобы понять, что функции numpy хотят работать с массивами numpy, где они могут делать циклы на скорости C и делать это быстрее, чем на отдельные элементы.

def resampleMatrix(Tso, Tsf, o, f):
    retval = numpy.sinc((Tsi*numpy.arange(i)[:,numpy.newaxis]
                         -Tso*numpy.arange(j)[numpy.newaxis,:])/Tso)
    return retval

Хитрость в том, что, индексируя апельсины с помощью numpy.newaxis, numpy преобразует массив с формой i в единицу с формой i x 1, а массив с формой j - в форму 1 x j. На этапе вычитания numpy будет «транслировать» каждый вход, чтобы действовать как массив в форме i x j, и выполнять вычитание. («Трансляция» - это термин «numpy», отражающий тот факт, что для увеличения размера i x 1 до i x j не создается дополнительная копия).

Теперь numpy.sinc может перебирать все элементы в скомпилированном коде, гораздо быстрее, чем любой цикл for, который вы можете написать.

(Существует дополнительное ускорение, если вы выполняете деление до вычитания, тем более что в последнем случае деление отменяется.)

Единственным недостатком является то, что теперь вы платите за дополнительный массив Nx10 * N, чтобы удержать разницу. Это может быть нарушителем, если N большое, а память - проблема.

В противном случае вы сможете написать это, используя numpy.convolve. Из того, что я только что узнал о sinc-интерполяции, я бы сказал, что вы хотите что-то вроде numpy.convolve(orig,numpy.sinc(numpy.arange(j)),mode="same"). Но я, вероятно, ошибаюсь по поводу специфики.

Answer 5

Вот минимальный пример 1d-интерполяции со scipy - не так весело, как переизобретать, но.
Сюжет выглядит как sinc, что не случайно: попробуйте Google Spline Resample "Approx Sinc".
(Предположительно меньше локальных / больше нажатий & rArr; лучшее приближение, но я понятия не имею, как локальные UnivariateSplines.)

""" interpolate with scipy.interpolate.UnivariateSpline """
from __future__ import division
import numpy as np
from scipy.interpolate import UnivariateSpline
import pylab as pl

N = 10 
H = 8
x = np.arange(N+1)
xup = np.arange( 0, N, 1/H )
y = np.zeros(N+1);  y[N//2] = 100

interpolator = UnivariateSpline( x, y, k=3, s=0 )  # s=0 interpolates
yup = interpolator( xup )
np.set_printoptions( 1, threshold=100, suppress=True )  # .1f
print "yup:", yup

pl.plot( x, y, "green",  xup, yup, "blue" )
pl.show()

Добавлен февраль 2010: см. Также basic-spline-interpolation-in-a-little-of-numpy

Answer 6

Ваш вопрос не совсем понятен; вы пытаетесь оптимизировать код, который вы разместили, верно?

Переписывание sinc, как это должно значительно ускорить его. Эта реализация избегает проверки того, что математический модуль импортируется при каждом вызове, не выполняет доступ к атрибутам три раза и заменяет обработку исключений условным выражением:

from math import sin, pi
def sinc(x):
    return (sin(pi * x) / (pi * x)) if x != 0 else 1.0

Можно также попытаться избежать создания матрицы дважды (и удерживать ее дважды в памяти параллельно), создавая массив numpy.array напрямую (не из списка списков):

def resampleMatrix(Tso, Tsf, o, f):
    retval = numpy.zeros((f, o))
    for i in xrange(f):
        for j in xrange(o):
            retval[i][j] = sinc((Tsf*i - Tso*j)/Tso)
    return retval

(заменить xrange диапазоном на Python 3.0 и выше)

Наконец, вы можете создавать строки с помощью numpy.arange, а также вызывать numpy.sinc для каждой строки или даже для всей матрицы:

def resampleMatrix(Tso, Tsf, o, f):
    retval = numpy.zeros((f, o))
    for i in xrange(f):
        retval[i] = numpy.arange(Tsf*i / Tso, Tsf*i / Tso - o, -1.0)
    return numpy.sinc(retval)

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

Answer 7

Я рекомендую вам проверить свой алгоритм, так как это нетривиальная проблема. В частности, я предлагаю вам получить доступ к статье Ху и Павлидиса (1991) «Разработка функций с использованием конических сплайнов» (IEEE Computer Graphics and Applications). Их реализация алгоритма позволяет адаптивную выборку функции, так что время рендеринга меньше, чем при подходах с регулярным интервалом.

Резюме следует:

Представлен метод, в соответствии с которым математическое описание функция, аппроксимирующая конический сплайн сюжет функции производится. Конические дуги были выбраны в качестве примитивные кривые, потому что есть простые алгоритмы инкрементного построения для коник уже включены в некоторые драйверы устройств, и есть простые алгоритмы для локальных приближений коники. Алгоритм разделения и слияния для адаптивного выбора узлов, согласно анализу формы оригинальная функция на основе его производные первого порядка, это введено.

Answer 8

Небольшое улучшение. Используйте встроенную функцию numpy.sinc (x), которая выполняется в скомпилированном коде C.

Возможное большее улучшение: Можете ли вы сделать интерполяцию на лету (как происходит построение графика)? Или вы привязаны к библиотеке печати, которая принимает только матрицу?