scipy.ndimage.filters.convolve 和乘法傅里叶变换给出不同的结果

Mya*_*ath 5 python numpy convolution scipy

这是我的代码:

from scipy.ndimage import filters
import numpy

a = numpy.array([[2,43,42,123,461],[453,12,111,123,55] ,[123,112,233,12,255]])
b = numpy.array([[0,2,2,3,0],[0,15,12,100,0],[0,45,32,22,0]])

ab = filters.convolve(a,b, mode='constant', cval=0)

af = numpy.fft.fftn(a)
bf = numpy.fft.fftn(b)

abf = af*bf

abif = numpy.fft.ifftn(abf)

print numpy.around(ab)
print numpy.around(abif)
Run Code Online (Sandbox Code Playgroud)

结果是:

[[ 1599  2951  7153 13280 18311]
 [ 8085 51478 13028 40239 30964]
 [18192 32484 23527 36122  8726]]

[[ 37416.+0.j  32251.+0.j  46375.+0.j  32660.+0.j  23986.+0.j]
 [ 30265.+0.j  33206.+0.j  62450.+0.j  19726.+0.j  17613.+0.j]
 [ 40239.+0.j  38095.+0.j  24492.+0.j  51478.+0.j  13028.+0.j]]
Run Code Online (Sandbox Code Playgroud)

如何修复使用 FFT 进行卷积的方式,以保证它给出与 相同的结果scipy.ndimage.filters.convolve

谢谢。

sen*_*rle 5

事实证明这是一个有趣的问题。使用离散傅里叶变换(由 实现numpy.fft.fftn)的卷积似乎相当于循环卷积。所以我们需要做的就是使用'wrap'卷积模式,并适当设置原点:

>>> filters.convolve(a, b, mode='wrap', origin=(-1, -2))
array([[37416, 32251, 46375, 32660, 23986],
       [30265, 33206, 62450, 19726, 17613],
       [40239, 38095, 24492, 51478, 13028]])
>>> numpy.fft.ifftn(numpy.fft.fftn(a) * numpy.fft.fftn(b))
array([[ 37416.+0.j,  32251.+0.j,  46375.+0.j,  32660.+0.j,  23986.+0.j],
       [ 30265.+0.j,  33206.+0.j,  62450.+0.j,  19726.+0.j,  17613.+0.j],
       [ 40239.+0.j,  38095.+0.j,  24492.+0.j,  51478.+0.j,  13028.+0.j]])
>>> (filters.convolve(a, b, mode='wrap', origin=(-1, -2)) ==
...  numpy.around(numpy.fft.ifftn(numpy.fft.fftn(a) * numpy.fft.fftn(b))))
array([[ True,  True,  True,  True,  True],
       [ True,  True,  True,  True,  True],
       [ True,  True,  True,  True,  True]], dtype=bool)
Run Code Online (Sandbox Code Playgroud)

差异仅与处理边缘的方式有关filters.convolve。在其他模式下执行卷积的方法fftn并没有立即引起我的注意;对于解决该问题的聪明(且事后显而易见)的方法,请参阅Warren Weckesser的出色答案。


War*_*ser 4

正如@senderle指出的,当你使用FFT来实现卷积时,你得到的是循环卷积。@senderle 的答案显示了如何调整参数filters.convolve来进行循环卷积。要修改 FFT 计算以生成与原始使用相同的结果filters.convolve,可以用 0 填充参数,然后提取结果的适当部分:

from scipy.ndimage import filters
import numpy

a = numpy.array([[2.0,43,42,123,461], [453,12,111,123,55], [123,112,233,12,255]])
b = numpy.array([[0.0,2,2,3,0], [0,15,12,100,0], [0,45,32,22,0]])

ab = filters.convolve(a,b, mode='constant', cval=0)

print numpy.around(ab)
print

nrows, ncols = a.shape
# Assume b has the same shape as a.
# Pad the bottom and right side of a and b with zeros.
pa = numpy.pad(a, ((0, nrows-1), (0, ncols-1)), mode='constant')
pb = numpy.pad(b, ((0, nrows-1), (0, ncols-1)), mode='constant')
paf = numpy.fft.fftn(pa)
pbf = numpy.fft.fftn(pb)
pabf = paf*pbf
p0 = nrows // 2
p1 = ncols // 2
pabif = numpy.fft.ifftn(pabf).real[p0:p0+nrows, p1:p1+ncols]

print pabif
Run Code Online (Sandbox Code Playgroud)

输出:

[[  1599.   2951.   7153.  13280.  18311.]
 [  8085.  51478.  13028.  40239.  30964.]
 [ 18192.  32484.  23527.  36122.   8726.]]

[[  1599.   2951.   7153.  13280.  18311.]
 [  8085.  51478.  13028.  40239.  30964.]
 [ 18192.  32484.  23527.  36122.   8726.]]
Run Code Online (Sandbox Code Playgroud)

  • 查看 [`scipy.signal.fftconvolve`](https://github.com/scipy/scipy/blob/v0.15.1/scipy/signal/signaltools.py#L350) 的源代码也很有启发性,其中显示基于原生 fft 的卷积会自动执行此填充(到最接近的 2 的幂),并且当要求输出的形状“相同”(与输入形状匹配)时,它会提取相关的子部分结果。 (2认同)