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

Posted 2023-03-24

【中文标题】有没有办法用 scipy.fft 在傅里叶空间中进行数值积分？

【英文标题】：Is there a way to numerically integrate in Fourier space with scipy.fft?

【发布时间】：2021-09-28 11:42:55

【问题描述】：

在使用 scipy 获取一些数据后，我有兴趣在傅立叶空间中进行积分。我一直在关注这个堆栈交换帖子numerical integration in Fourier space with 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 在傅立叶空间中正确集成函数？

【问题讨论】：

看起来您的集成和标准化常数已关闭...

可能不想将 DC 项设置为零

此外，您的“倾斜”实际上是在为所有内容添加线性偏差。

这些 cmets 没有帮助，只是重申我在代码中尝试过的问题/事情。如果您没有任何实质内容要添加，我建议您不要发表评论（在此帖子和将来的帖子中）。

【参考方案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 不需要填充，函数在采样域的边缘达到零。我们也可以将这种行为强加于其他情况。在离散傅里叶分析中，这称为windowing。这看起来像这样：

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()

（不显示输出：两个图完全重叠。）

【讨论】：

您能否详细说明您关于循环卷积的观点？将此 FFT 乘以 f 的任何函数都是卷积，对吗？这与循环卷积有何区别？我同意这归结为一个边界问题，但我不明白循环卷积如何与之相关。

@MSPaintPhysics 当卷积核到达采样域之外时，它会环绕并使用域另一侧的值。这称为循环卷积，相当于无限重复函数的采样部分。