scipy.signal.fftconvolve不能给出所需的结果。

Question

我有一个关于Python的fftconvolve的问题。在我的当前研究中，我需要计算两个函数之间的某些卷积。为此，我使用傅里叶变换进行计算（使用了numpy.fft并对其进行了归一化）。问题是，如果我想使用fftconvolve包进行比较，它无法给出正确的结果。以下是我的代码：

#!/usr/bin/python
import numpy as np
from scipy.signal import fftconvolve , convolve 

def FFT(array , sign):
  if sign==1:
    return np.fft.fftshift(np.fft.fft(np.fft.fftshift(array))) * dw / (2.0 * np.pi)
  elif sign==-1:
    return np.fft.fftshift(np.fft.ifft(np.fft.fftshift(array))) * dt * len(array)


def convolve_arrays(array1,array2,sign):
  sign = int(sign)
  temp1 = FFT(array1 , sign,)
  temp2 = FFT(array2 , sign,)
  temp3 = np.multiply(temp1 , temp2)
  return  FFT(temp3 , -1 * sign) / (2. * np.pi) 

""" EXAMPLE """ 

dt    = .1
N     = 2**17
t_max = N * dt / 2
time  = dt * np.arange(-N / 2 , N / 2 , 1)

dw    = 2. * np.pi / (N * dt)
w_max = N * dw / 2.
w     = dw * np.arange(-N / 2 , N / 2 , 1)

eta_fourier = 1e-10




Gamma   = 1.
epsilon = .5
omega   = .5


G    = zeros(N , complex)
G[:] = 1. / (w[:] - epsilon + 1j * eta_fourier)

D    = zeros(N , complex)
D[:] = 1. / (w[:] - omega + 1j * eta_fourier) - 1. / (w[:] + omega + 1j * eta_fourier)

H    = convolve_arrays(D , G , 1)     
J    = fftconvolve(D , G , mode = 'same') * np.pi  / (2. * N) 

如果您绘制 H 和 J 的实部/虚部，您会看到 w 轴发生了偏移，并且为了尽可能接近正确的结果，我不得不乘以 J 的结果。

有什么建议吗？

谢谢！

- user2857354

---

请看一下scipy.fftconvolve函数，注意到它的算法不包含你所使用的奇怪的傅里叶变换移位或缩放。你想通过使用FFT函数达到什么目的？ - Henry Gomersall

Answer 1

当计算卷积时，边界条件非常重要。

当你将两个信号进行卷积时，结果的边缘取决于你假设输入边缘之外的值是什么。fftconvolve 计算卷积时假设边界为零填充。

看一下 fftconvolve 的源代码。注意他们为了实现零填充边界条件所做的花招，特别是这些行：

size = s1 + s2 - 1

...

fsize = 2 ** np.ceil(np.log2(size)).astype(int) #For speed; often suboptimal!
fslice = tuple([slice(0, int(sz)) for sz in size])

...

ret = ifftn(fftn(in1, fsize) * fftn(in2, fsize))[fslice].copy()

...

return _centered(ret, s1) #strips off padding

这是很好的东西！如果想要理解基于傅里叶变换的卷积，建议仔细阅读fftconvolve的代码，并接受良好的教育。

简要概述：

正向FFT对每个信号进行零填充以避免周期性边界条件：

a = np.array([3, 4, 5])
b = np.fft.ifftn(np.fft.fftn(a, (5,))).real
print b #[ 3.  4.  5.  0.  0.]

正向FFT的乘积的反向FFT会给出一个填充结果：

a = np.array([3, 4, 5])
b = np.array([0., 0.9, 0.1])
b = np.fft.ifftn(np.fft.fftn(a, (5,))*
                 np.fft.fftn(b, (5,))
                 ).real
print b #[ 0.   2.7  3.9  4.9  0.5]

而_centered函数会在结尾剥离多余的填充像素（假设您使用mode='same'选项）。

- Andrew

---

看起来 fftconvolve 的新版本在选择填充量以获得快速FFT计算方面做得更好。 (https://github.com/scipy/scipy/blob/v0.15.1/scipy/signal/signaltools.py#L349) - Andrew