时间戳

首页 > 解决方案 > 使用 fftconvolve 评估两个连续函数的卷积

python - 使用 fftconvolve 评估两个连续函数的卷积

问题描述

我正在尝试使用 评估两个连续函数的卷积scipy.signal.fftconvolve。代码的场景如下：

我正在尝试近似以下双重积分：

在此处输入图像描述

，即在区域 C_1(x',y') 中，表示以 (x', y') 为中心的半径为 1 的圆。这可以通过以下积分来近似：

在此处输入图像描述

其中函数 K 被选为连续可积函数，例如 ，exp(-x^2-y^2)其形状近似为半径为 1 的圆。如果我采用函数 K'(x,y)=K(-x,-y) ，那么积分就是两个函数的卷积：

在此处输入图像描述

所以我尝试将这两个函数离散化成数组，然后进行卷积。

以下代码将用 Julia 编写，fftconvolve函数将使用PyCall.jl.

using PyCall
using Interpolations

r = 1
xc = -10:0.05:10
yc = -10:0.05:10

K(x, y) = exp(-(x^2+y^2)/r^2)
rho(x, y) = x^2+y^3    # Try some arbitrary function

ss = pyimport("scipy.signal")  # Import scipy.signal module from Python
a = [rho(x,y) for x in xc, y in yc]
b = [K(-x,-y) for x in xc, y in yc]
c = ss.fftconvolve(a,b,mode="same")  # zero-paddings beyond boundary, unimportant since rho is near zero beyond the boundary anyway
c_unscaled = interpolate(c', BSpline(Cubic(Line(OnCell())))) 
# Adjoint because the array comprehension switched x and y, then interpolate the array
c_scaled = Interpolations.scale(c_unscaled, xc, yc) # Scale the interpolated function w.r.t. xc, yc
print(c_scaled(0.0,0.0)) # The result of the integral for (x', y') = (0, 0)

结果是628.3185307178969，而数值积分的结果是0.785398。这里有什么问题？

标签： pythonscipyjuliafftconvolution

解决方案

您可能会尝试使用 scipy.signal.convolve 来对两个 N 维数组进行卷积，但不能使用快速傅立叶变换。它使用直接方法来计算卷积。在这里，我的意思是卷积是直接从总和中确定的。

因此，您可以尝试将计算 c 的行替换为以下行：

c = ss.convolve(a,b,mode="same", method='direct')

推荐阅读