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

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

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

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

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

Здесь полный script я использовал для быстрого проверки:

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)

Ответы

Ответ 1

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

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 изменит ответ. Я слышал рассказы о том, что есть способы справиться с этим, умело используя простые факторы длины (упомянутые, но не закодированные в Numerical Recipes), но я никогда не видел, чтобы люди действительно это делали.

Ответ 2

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

Ответ 3

Как упоминалось ранее, функция 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')