使用 fft 的二阶导数

问题描述

所有，我正在尝试采用以下函数的拉普拉斯算子：

g(x,y) = 1/2cx^2+1/2dy2

拉普拉斯算子是 c + d，这是一个常数。使用 FFT 我应该得到相同的结果（在我的 FFT 示例中，我填充函数以避免边缘效应）。

这是我的代码：

public class Alert{
    public Guid AlertId{get;set;}
    public List<Info> Infos {get;set;}

    public Alert() {
        Infos = new List<Info>();
    }
}

[

第一张图是原函数g(x,y)的绘图，第二张图是g的解析拉普拉斯算子，第三张图是里约热内卢的糖面包（笑），实际上是使用FFT的拉普拉斯算子。我在这里做错了什么？

编辑：评论涟漪效应。 Cris 你的意思是下图中 set_zlimit 引起的涟漪效应？只是记住你的结果应该是 8。

编辑 2：使用非对称的 x 和 y 值，生成两个图像。

解决方法

填充不会改变边界条件：您通过复制函数进行填充，镜像，四次。该函数是对称的，因此镜像不会改变它。因此，您的填充只是将函数重复四次。通过 DFT（您正在尝试实现）的卷积使用周期性边界条件，因此已经将输入函数视为周期性。复制函数不会改善边缘的卷积结果。

为了改善边缘的结果，您需要实现不同的边界条件，最有效的一个（因为输入无论如何都是解析的）是简单地扩展域，然后在应用卷积后裁剪它。这引入了边界扩展，通过查看原始域之外的更多数据来填充图像。它是一种理想的边界扩展，适用于我们不必处理现实世界数据的理想情况。

这通过大大简化的代码通过 DFT 实现了拉普拉斯，我们忽略任何边界扩展，以及样本间距（基本上设置 dx=1 和 dy=1）：

import numpy as np
import matplotlib.pyplot as pp

n = 30 # number of points
c = 4
d = 4
x = np.arange(-n//2,n//2)
y = np.arange(-n//2,n//2)
g = (1/2)*c*x[None,:]**2 + (1/2)*d*y[:,None]**2

kx = 2 * np.pi * np.fft.fftfreq(n)
ky = 2 * np.pi * np.fft.fftfreq(n)
lapg = np.real(np.fft.ifft2(np.fft.fft2(g) * (-kx[None,:]**2 - ky[:,None]**2)))

fig = pp.figure()
ax = fig.add_subplot(121,projection='3d')
ax.plot_surface(x[None,:],y[:,None],g)
ax = fig.add_subplot(122,lapg)
pp.show()

---

编辑：边界扩展的工作方式如下：

import numpy as np
import matplotlib.pyplot as pp

n_true = 30 # number of pixels we want to compute
n_boundary = 15 # number of pixels to extend the image in all directions
c = 4
d = 4

# First compute g and lapg including boundary extenstion
n = n_true + n_boundary * 2
x = np.arange(-n//2,None]**2
kx = 2 * np.pi * np.fft.fftfreq(n)
ky = 2 * np.pi * np.fft.fftfreq(n)
lapg = np.real(np.fft.ifft2(np.fft.fft2(g) * (-kx[None,None]**2)))

# Now crop the two images to our desired size
x = x[n_boundary:-n_boundary]
y = y[n_boundary:-n_boundary]
g = g[n_boundary:-n_boundary,n_boundary:-n_boundary]
lapg = lapg[n_boundary:-n_boundary,n_boundary:-n_boundary]

# Display
fig = pp.figure()
ax = fig.add_subplot(121,g)
ax.set_zlim(0,800)
ax = fig.add_subplot(122,lapg)
ax.set_zlim(0,800)
pp.show()

请注意，我以相同的方式缩放两个图的 z 轴，以免过多地增强边界的影响。像这样的傅立叶域滤波通常比空间域（或时域）滤波对边缘效应更敏感，因为滤波器具有无限长的脉冲响应。如果您省略 set_zlim 命令，您会在原本平坦的 lapg 图像中看到涟漪效应。涟漪非常小，但无论多小，在完全平坦的函数上它们看起来都会很大，因为它们会从图的底部延伸到顶部。两个图中相等的 set_zlim 只是将噪声按比例放置。