Cython и скорость numpy

Изменить:
Теперь есть scipy.signal.fftconvolve, который был бы предпочтительным подходом к реализации подхода свертки на основе БПФ, который я описываю ниже. Я оставлю исходный ответ, чтобы объяснить проблему скорости, но на практике используйте scipy.signal.fftconvolve.

Исходный ответ:
Использование БПФ и теорема свертки даст вам резкое увеличение скорости за счет преобразования задачи из O (n ^ 2) в O (n log n). Это особенно полезно для длинных наборов данных, таких как ваш, и может дать прирост скорости на 1000 секунд или намного больше, в зависимости от длины. Это также легко сделать: просто БПФ обоих сигналов умножьте и обратное БПФ произведите. numpy.correlate не использует метод БПФ в программе взаимной корреляции, и его лучше использовать с очень маленькими ядрами.

Вот пример

from timeit import Timer
from numpy import *

times = arange(0, 100, .001)

xdata = 1.*sin(2*pi*1.*times) + .5*sin(2*pi*1.1*times + 1.)
ydata = .5*sin(2*pi*1.1*times)

def xcorr(x, y):
    return correlate(x, y, mode='same')

def fftxcorr(x, y):
    fx, fy = fft.fft(x), fft.fft(y[::-1])
    fxfy = fx*fy
    xy = fft.ifft(fxfy)
    return xy

if __name__ == "__main__":
    N = 10
    t = Timer("xcorr(xdata, ydata)", "from __main__ import xcorr, xdata, ydata")
    print 'xcorr', t.timeit(number=N)/N
    t = Timer("fftxcorr(xdata, ydata)", "from __main__ import fftxcorr, xdata, ydata")
    print 'fftxcorr', t.timeit(number=N)/N

Это дает время работы на цикл (в секундах для 10 000 длинных сигналов).

xcorr 34.3761689901
fftxcorr 0.0768054962158

Понятно, что метод fftxcorr намного быстрее.

Если вы нанесете на график результаты, вы увидите, что они очень похожи при почти нулевом временном сдвиге. Обратите внимание, однако, что по мере удаления xcorr будет уменьшаться, а fftxcorr - нет. Это потому, что немного неоднозначно, что делать с частями сигнала, которые не перекрываются при смещении сигналов. xcorr рассматривает его как ноль, а БПФ рассматривает формы волны как периодические, но если это проблема, ее можно исправить с помощью заполнения нулями.

Уловка с подобными вещами состоит в том, чтобы найти способ разделять и властвовать.

В настоящее время вы перемещаетесь в каждую позицию и проверяете каждую точку в каждой позиции - фактически операция O (n ^ 2).

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

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

У вас может быть «чек на каждую 8-ю позицию», который вы умножаете на 8. Если это слишком мало, пропустите его и двигайтесь дальше. Если это достаточно много, проверьте все значения, чтобы увидеть, нашли ли вы максимумы.

Проблема в том, сколько времени требуется для выполнения всех этих умножений - (f[<unsigned int>(i+j)] * g[j]) Фактически, вы заполняете большую матрицу всеми этими продуктами и выбираете строку с максимальной суммой. Вы не хотите подсчитывать «все» продукты. Достаточно товаров, чтобы убедиться, что вы нашли максимальную сумму.

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

(Думаю, это может сработать. Я не пробовал.)

Если бы вы использовали max(g)-g[j] для работы с отрицательными числами, вы бы искали самые маленькие, а не самые большие. Вы можете вычислить корреляцию для первой позиции. Все, что суммировалось с большим значением, можно было немедленно остановить - больше никаких умножений или сложений для этого смещения, переходите к другому.

• вы можете извлечь диапазон (размер2) из ​​внешнего цикла
• вы можете использовать sum () вместо цикла для вычисления current_correlation
• вы можете хранить корреляции и задержки в списке, а затем использовать max (), чтобы получить самый большой