16

分析我正在做的一些计算工作向我表明,我的程序中的一个瓶颈是一个基本上可以做到这一点的函数(npnumpyspscipy):

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)
4

3 回答 3

11

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

1) numpy.convolve 不是循环的,这是 fft 代码给你的:

2) FFT 不会在内部填充 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 的幂会改变答案。我听说过有一些方法可以通过巧妙地使用长度的主要因素(在数字食谱中提到但未编码)来解决这个问题,但我从未见过人们真正这样做。

于 2011-07-28T07:27:04.520 回答
2

scipy.signal.fftconvolve 确实通过 FFT 进行卷积,它是 python 代码。您可以研究源代码,并更正您的 mix1 功能。

于 2011-07-29T00:40:31.547 回答
1

如前所述,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')
于 2013-08-20T15:21:25.750 回答