有没有办法用 scipy.fft 在傅立叶空间中进行数值积分？

Question

我有兴趣在使用 scipy 对一些数据进行快速傅里叶变换后集成到傅立叶空间中。我一直在关注这个堆栈交换帖子numpy.fft在傅立叶空间中的数值集成，但它没有正确集成我一直在使用的一些测试用例。我添加了几行来解决这个问题，但仍然没有恢复正确的积分。下面是我用来集成测试用例的代码。代码顶部是我一直在使用的 3 个测试用例。

import numpy as np
import scipy.special as sp
from scipy.fft import fft, ifft, fftfreq
import matplotlib.pyplot as plt

#set number of points in array
Ns = 2**16 
#create array in space
x = np.linspace(-np.pi, np.pi, Ns)

#test case 1 from stack exchange post
# y = np.exp(-x**2)                          # function f(x)
# ys = np.exp(-x**2) * (-2 *x)               # derivative f'(x)

#test case 2 
# y = np.exp(-x**2) * x - 1/2 *np.sqrt(np.pi)*sp.erf(x)
# ys =  np.exp(-x**2) * -2*x**2

#test case 3
y = np.sin(x**2) + (1/4)*np.exp(x)
ys = 1/4*(np.exp(x) + 8*x*np.cos(x**2))


#find spacing in space array
ss = x[1]-x[0]

#definte fft integration function
def fft_int(N,s,dydt):
    #create frequency array
    f = fftfreq(N,s)

    # integration step ignoring divide by 0 errors
    Fys = fft(dydt)
    with np.errstate(divide="ignore", invalid="ignore"):
        modFys = Fys / (2*np.pi*1j*f) 

    #set DC term to 0, was a nan since we divided by 0
    modFys[0] = 0

    #take inverse fft and subtract by integration constant
    fourier = ifft(modFys)
    fourier = fourier-fourier[0]
    
    #tilt correction if function doesn't approach 0 at its ends
    tilt = np.sum(dydt)*s*(np.arange(0,N)/(N-1) - 1/2)
    fourier = fourier + tilt
    
    return fourier

测试用例 1 来自上面的堆栈交换帖子。如果您复制粘贴顶部答案中的代码并绘制，您将得到如下所示的结果：

蓝色实线为fft积分方法，橙色虚线为解析解。我用下面的代码行来解释这个偏移量：

fourier = fourier-fourier[0]

因为我不相信代码设置了积分常数。

接下来，对于测试用例 2，我得到如下图：

同样，蓝色实线是 fft 积分方法，橙色虚线是解析解。我使用以下代码行来解释解决方案中的这种倾斜

tilt = np.sum(dydt)*s*(np.arange(0,N)/(N-1) - 1/2)
fourier = fourier + tilt

最后我们到达测试用例 3。结果如下图所示：

同样，蓝色实线是 fft 积分方法，橙色虚线是解析解。这就是我陷入困境的地方，这个偏移量再次出现，我不知道为什么。

TLDR：如何使用 scipy.fft 在傅立叶空间中正确积分函数？

Answer 1

tilt

组件没有任何意义。它修复了一个功能，但它不是问题的通用解决方案。

问题在于 FFT 会引起信号的周期性，这意味着您计算不同函数的积分。将信号的 FFT 乘以

1/(2*np.pi*1j*f)

相当于信号与

ifft(1/(2*np.pi*1j*f))

的循环卷积。 “循环”是这里的关键。这只是一个边界问题。

用零填充函数是尝试解决此问题的一种方法：

import numpy as np
import scipy.special as sp
from scipy.fft import fft, ifft, fftfreq
import matplotlib.pyplot as plt

def fft_int(s, dydt, N=0):
    dydt_padded = np.pad(dydt, (0, N))
    f = fftfreq(dydt_padded.shape[0], s)
    F = fft(dydt_padded)
    with np.errstate(divide="ignore", invalid="ignore"):
        F = F / (2*np.pi*1j*f) 
    F[0] = 0
    y_padded = np.real(ifft(F))
    y = y_padded[0:dydt.shape[0]]
    return y - np.mean(y)

N = 2**16
x = np.linspace(-np.pi, np.pi, N)
s = x[1] - x[0]

# Test case 3
y = np.sin(x**2) + (1/4)*np.exp(x)
dy = 1/4*(np.exp(x) + 8*x*np.cos(x**2))
plt.plot(y - np.mean(y))
plt.plot(fft_int(s, dy))
plt.plot(fft_int(s, dy, N))
plt.plot(fft_int(s, dy, 10*N))
plt.show()

（蓝色是预期输出，没有填充的计算解决方案是橙色，随着填充量的增加，计算解决方案是绿色和红色。）

在这里，我通过绘制所有删除均值的函数来解决“偏移”问题。将 DC 分量设置为 0 等于减去平均值。但在裁剪掉填充后，平均值会发生变化，因此

fft_int

在裁剪后再次减去平均值。

无论如何，请注意随着填充的增加，我们如何获得越来越好的近似值。为了获得准确的结果，需要无限量的填充，这当然是不现实的。

测试用例 #1 不需要填充，函数在采样域的边缘达到零。我们也可以将这种行为强加于其他情况。在离散傅立叶分析中，这称为“加窗”。这看起来像这样：

def fft_int(s, dydt):
    dydt_windowed = dydt * np.hanning(dydt.shape[0])
    f = fftfreq(dydt.shape[0], s)
    F = fft(dydt_windowed)
    with np.errstate(divide="ignore", invalid="ignore"):
        F = F / (2*np.pi*1j*f) 
    F[0] = 0
    y = np.real(ifft(F))
    return y

但是，这里我们仅在域的中间获得正确的积分结果，并且向末端逐渐抑制值。所以这也不是一个实用的解决方案。

我的结论是，不，这是不可能的。使用

np.cumsum

:

 计算积分要容易得多

yp = np.cumsum(dy) * s

plt.plot(y - np.mean(y))
plt.plot(yp - np.mean(yp))
plt.show()

（未显示输出：两个图完美重叠。）