Numpy/Scipy中的卷积计算

Question

分析我正在做的一些计算工作向我展示了我的程序中的一个瓶颈是基本上这样做的功能(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)

Answer 1

所以我测试了这个,现在可以确认一些事情:

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)归一化不是一个区别 - 如果你通过加上一个(k)*b(ik)进行一个简单的循环卷积,你将得到FFT代码的结果.

事情是填充2的力量将改变答案.我听说有很多方法可以通过巧妙地使用长度的素数因子(在数字食谱中提到但没有编码)来处理这个问题,但我从未见过人们真的这样做.