python - 使用 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 警告。这是一种解决方法，对于上面的示例似乎没问题，但对于我遇到的更复杂的问题却失败了。一位同事告诉我，如果我得到偏移频率值的 fft (np.arange(N) - N/2. + 1./2)/T，我可以简单地避免它.如何在 numpy 中做到这一点？有没有办法指定 numpy fft 的输出网格？

谢谢!

Answer 1

问题是 w 包含 0(它应该如此)，并且您除以 w。 w 中的 0 是“DC”频率；对应傅里叶级数的常数项。

如果您对具有系数 A0 的非零 DC 分量的函数进行积分，则生成的函数包含 A0*t 形式的项，该项不在该傅立叶技术适用的周期函数空间中。因此，您必须假设输入的直流分量为 0。在这种情况下，当您除以 w 时，您将得到 (0+0j)/0，即 (nan +nanj)。如果输入的直流分量不为零，您将得到 (inf+nanj)。无论哪种方式，解决方案都是简单地忽略你得到的任何东西，并在使用 ifft 反转之前将 DC 傅里叶系数设置为 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 .