Вычисления свертки в Numpy / Scipy

Профилирование некоторой вычислительной работы, которую я делаю, показало мне, что одним узким местом в моей программе была функция, которая в основном делала это ( npnumpy , sp scipy ):

 def mix1(signal1, signal2): spec1 = np.fft.fft(signal1, axis=1) spec2 = np.fft.fft(signal2, axis=1) return np.fft.ifft(spec1*spec2, axis=1) 

Оба сигнала имеют форму (C, N) где C – количество наборов данных (обычно меньше 20), а N – количество выборок в каждом наборе (около 5000). Вычисление для каждого набора (строки) полностью не зависит от какого-либо другого набора.

Я понял, что это была просто свертка, поэтому я попытался ее заменить:

 def mix2(signal1, signal2): outputs = np.empty_like(signal1) for idx, row in enumerate(outputs): outputs[idx] = sp.signal.convolve(signal1[idx], signal2[idx], mode='same') return outputs 

… просто чтобы узнать, получили ли я те же результаты. Но я этого не делал, и мои вопросы:

  1. Почему нет?
  2. Есть ли лучший способ вычислить эквивалент mix1() ?

(Я понимаю, что mix2 вероятно, не был бы быстрее, как есть, но, возможно, это была хорошая отправная точка для параллелизации.)

Вот полный сценарий, который я использовал, чтобы быстро проверить это:

 import numpy as np import scipy as sp import scipy.signal N = 4680 C = 6 def mix1(signal1, signal2): spec1 = np.fft.fft(signal1, axis=1) spec2 = np.fft.fft(signal2, axis=1) return np.fft.ifft(spec1*spec2, axis=1) def mix2(signal1, signal2): outputs = np.empty_like(signal1) for idx, row in enumerate(outputs): outputs[idx] = sp.signal.convolve(signal1[idx], signal2[idx], mode='same') return outputs def test(num, chans): sig1 = np.random.randn(chans, num) sig2 = np.random.randn(chans, num) res1 = mix1(sig1, sig2) res2 = mix2(sig1, sig2) np.testing.assert_almost_equal(res1, res2) if __name__ == "__main__": np.random.seed(0x1234ABCD) test(N, C) 

3 Solutions collect form web for “Вычисления свертки в Numpy / Scipy”

Поэтому я проверил это и теперь могу подтвердить несколько вещей:

1) numpy.convolve не является циркулярным, что дает вам код fft:

2) БПФ не внутренне прокладывает до мощности 2. Сравните значительно разные скорости следующих операций:

 x1 = np.random.uniform(size=2**17-1) x2 = np.random.uniform(size=2**17) np.fft.fft(x1) np.fft.fft(x2) 

3) Нормализация не является разницей – если вы делаете наивную круговую свертку, добавляя a (k) * b (ik), вы получите результат кода FFT.

Дело в том, что отступы до степени 2 будут менять ответ. Я слышал рассказы о том, что есть способы справиться с этим, умело используя простые факторы длины (упомянутые, но не закодированные в Numerical Recipes), но я никогда не видел, чтобы люди действительно это делали.

scipy.signal.fftconvolve выполняет свертку с помощью FFT, это код python. Вы можете изучить исходный код и исправить функцию mix1.

Как упоминалось ранее, функция scipy.signal.convolve не выполняет круговую свертку. Если вам нужна круговая свертка, выполняемая в реальном пространстве (в отличие от использования fft), я предлагаю использовать функцию scipy.ndimage.convolve. Он имеет параметр режима, который может быть установлен на «wrap», что делает его круговой сверткой.

 for idx, row in enumerate(outputs): outputs[idx] = sp.ndimage.convolve(signal1[idx], signal2[idx], mode='wrap') 
Interesting Posts

Питоновский способ написать цикл for, который не использует индекс цикла

Совместимость функций OpenCV для нескольких изображений

Значение пытается быть установлено на копии среза из DataFrame

Показать последовательные изображения / массивы с imshow как повторяющиеся анимации в python

Пропустить каждый n-й индекс массива numpy

python cgitb не работает через браузер

Почему не удается передать * args и ** kwargs в __init__ дочернего класса

Сколько элементов было очищено для start_url

Получить полное имя компьютера из буквы сетевого диска в python

Как мне вызвать один вид фляжки из другого?

Замените все элементы массива Python NumPy, которые больше некоторого значения

unicode и кодировка для персидского или арабского в python3

Автоматическое заполнение поля created_by с административным сайтом Django

Оптимальные методы импорта Python (и Django)

SelfReferenceProperty против ListProperty Google App Engine

Python - лучший язык программирования в мире.