Вычисления свертки в Numpy / Scipy
Профилирование некоторой вычислительной работы, которую я делаю, показало мне, что одним узким местом в моей программе была функция, которая в основном делала это ( 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). Вычисление для каждого набора (строки) полностью не зависит от какого-либо другого набора.
- Функция конденсированной матрицы для нахождения пар
- как я могу обновить scipy в winpython на окнах?
- Многопоточные вызовы целевой функции scipy.optimize.leastsq
- Matplotlib: укажите формат поплавков для тиков
- Matplotlib - Быстрый способ создания множества подзаголовков?
- Как быстро выполнить установку наименьших квадратов по множеству наборов данных?
Я понял, что это была просто свертка, поэтому я попытался ее заменить:
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 вероятно, не был бы быстрее, как есть, но, возможно, это была хорошая отправная точка для параллелизации.)
Вот полный сценарий, который я использовал, чтобы быстро проверить это:
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')