Numpy/Scipy中的卷积计算 [英] Convolution computations in Numpy/Scipy

问题描述

分析我正在做的一些计算工作，向我表明程序中的一个瓶颈基本上是一个实现此功能的函数(np是numpy，sp是scipy):

Profiling some computational work I'm doing showed me that one bottleneck in my program was a function that basically did this (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个).每个集合(行)的计算完全独立于任何其他集合.

Both signals have shape (C, N) where C is the number of sets of data (usually less than 20) and N is the number of samples in each set (around 5000). The computation for each set (row) is completely independent of any other set.

我认为这只是一个简单的卷积，所以我尝试将其替换为:

I figured that this was just a simple convolution, so I tried to replace it with:

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

...只是看看我是否得到了相同的结果.但是我没有，我的问题是:

...just to see if I got the same results. But I didn't, and my questions are:

1. 为什么不呢?
2. 是否有更好的方法来计算mix1()的当量?

1. Why not?
2. Is there a better way to compute the equivalent of mix1()?

(我意识到mix2可能不会保持原样更快，但它可能是并行化的一个很好的起点.)

(I realise that mix2 probably wouldn't have been faster as-is, but it might have been a good starting point for parallelisation.)

这是我用来快速检查此内容的完整脚本:

Here's the full script I used to quickly check this:

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)

推荐答案

所以我对此进行了测试，现在可以确认一些事情:

So I tested this out and can now confirm a few things:

1)numpy.convolve不是循环的，这是fft代码为您提供的:

1) numpy.convolve is not circular, which is what the fft code is giving you:

2)FFT在内部没有填充2的幂.比较以下操作的速度差异很大:

2) FFT does not internally pad to a power of 2. Compare the vastly different speeds of the following operations:

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代码的结果.

3) Normalization is not a difference -- if you do a naive circular convolution by adding up a(k)*b(i-k), you will get the result of the FFT code.

填充2的幂将改变答案.我听说过一些故事，有一些方法可以通过巧妙地使用长度的素数来解决这一问题(提到但未在《数字食谱》中进行编码)，但我从未见过有人真正这样做过.

The thing is padding to a power of 2 is going to change the answer. I've heard tales that there are ways to deal with this by cleverly using prime factors of the length (mentioned but not coded in Numerical Recipes) but I've never seen people actually do that.

这篇关于Numpy/Scipy中的卷积计算的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！