分析我正在做的一些计算工作表明,我的程序中的一个瓶颈是一个基本上做到了这一点的函数(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
...just查看我是否得到了相同的结果。但我没有,我的问题是:
(我意识到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)
发布于 2011-07-28 15:27:05
所以我对此进行了测试,现在可以确认几件事:
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(i-k)相加来进行简单的循环卷积,您将得到FFT代码的结果。
填充到2的幂将会改变答案。我听说过有一些方法可以通过巧妙地使用长度的质数因子(提到但没有用数字配方编码)来解决这个问题,但我从来没有见过人们真正这样做。
发布于 2011-07-29 08:40:32
scipy.signal.fftconvolve确实通过快速傅立叶变换进行卷积,这是python代码。你可以研究一下源代码,改正你的mix1函数。
发布于 2013-08-20 23:21:26
如前所述,scipy.signal.convolve函数不执行循环卷积。如果你想在实数空间中执行循环卷积(与使用快速傅立叶变换相反),我建议使用scipy.ndimage.convolve函数。它有一个模式参数,可以设置为“wrap”,使其成为循环卷积。
for idx, row in enumerate(outputs):
outputs[idx] = sp.ndimage.convolve(signal1[idx], signal2[idx], mode='wrap')
https://stackoverflow.com/questions/6855169
复制相似问题