使用numpy.fft在傅里叶空间中进行数值积分

Question

我想在傅里叶空间中整合一个函数进行数值积分。

下面的代码展示了一个可以正常工作的例子：

import numpy as np
from pylab import *
from numpy.fft import fft, ifft, fftshift, ifftshift

N = 2**16
x = np.linspace(- np.pi , np.pi,N)
y = np.exp(-x**2)                          # function f(x)
ys = np.exp(-x**2) * (-2*x)                # derivative f'(x)
T = x[-1] - x[0] # the whole range
w = (np.arange(N) - N /2.) / T + 0.00000001 # slightly shifted
# integration
fourier =  ifft(ifftshift(1./ ( 2 * np.pi * 1j * w) * fftshift(fft(ys)) )  ) 
# differentiation 
fourier2 =  ifft(ifftshift(( 2 * np.pi * 1j * w) * fftshift(fft(fourier)) )  ) 

你可能会注意到在频率 w 的定义中有一个 + 0.00000001。我需要这个，因为不加它的话会出现 ZeroDivisionError（除零错误）或者 numpy 的警告。这是一个临时解决办法，对于上面的例子来说似乎没问题，但在我遇到的一个更复杂的问题中就不行了。我的同事告诉我，如果我使用偏移的频率值 (np.arange(N) - N /2. + 1./2) / T，就可以简单地避免这个问题。那么在 numpy 中怎么做呢？有没有办法指定 numpy fft 的输出网格？

谢谢！

Answer 1

问题在于，w里有0（这是正常的），而你却用w去做除法。w里的0代表的是“直流”（DC）频率，它对应的是傅里叶级数中的常数项。

如果你对一个有非零直流成分的函数进行积分，假设这个成分的系数是A0，那么得到的函数里会有一个形如A0*t的项，而这个项不属于傅里叶技术适用的周期函数空间。所以你必须假设输入的直流成分是0。在这种情况下，当你用w去除时，你会得到(0+0j)/0，这个结果是(nan+nanj)。如果输入的直流成分不为零，你会得到(inf+nanj)。无论哪种情况，解决办法就是忽略你得到的结果，并在用ifft反转之前把直流傅里叶系数设为0。

实现这个有几种方法。其中一种方法是把这些代码行：

w = (np.arange(N) - N /2.) / T + 0.00000001 # slightly shifted
# integration
fourier =  ifft(ifftshift(1./ ( 2 * np.pi * 1j * w) * fftshift(fft(ys)) )  ) 

改成这样（我加了一些中间变量）：

w = (np.arange(N) - N /2.) / T

# integration
Fys = fft(ys)
with np.errstate(divide="ignore", invalid="ignore"):
    modFys = ifftshift(1./ (2 * np.pi * 1j * w) * fftshift(Fys))

# modFys[0] will hold the result of dividing the DC component of y by 0, so it
# will be nan or inf.  Setting modFys[0] to 0 amounts to choosing a specific
# constant of integration.
modFys[0] = 0

fourier = ifft(modFys).real

我还取了ifft结果的实部。从理论上讲，虚部应该都是0；但实际上，由于浮点运算的不精确性，它们会非常小但不为零。

顺便说一下，如果你不想自己实现这个技术，可以使用scipy.fftpack.diff。