Сверточные вычисления в Numpy / Scipy - PullRequest
Новая
       16

Сверточные вычисления в Numpy / Scipy

16 голосов
/ 28 июля 2011

Профилирование некоторой вычислительной работы, которую я выполняю, показало, что одним узким местом в моей программе была функция, которая в основном выполняла это (np - numpy, 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 ]

11 голосов
/ 28 июля 2011

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

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 (i-k), вы получите результат кода FFT.

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

2 голосов
/ 29 июля 2011

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

1 голос
/ 20 августа 2013

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

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