使用 fft 的二阶导数

Question

全部，我正在尝试取以下函数的拉普拉斯：

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

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

这是我的代码：

Define a 2D function
n = 30 # number of points
Lx = 30 # extension in x
Ly = 30 # extension in x
dx = n/Lx # Step in x 
dy = n/Ly # Step in x 
c=4
d=4
x=np.arange(-Lx/2,Lx/2)
y=np.arange(-Ly/2,Ly/2)
g = np.zeros((Lx,Ly))
lapg = np.zeros((Lx,Ly))

for j in range(Ly):
    for i in range(Lx):
        g[i,j] = (1/2)*c*x[i]**2 + (1/2)*d*y[j]**2
        lapg[i,j] = c + d

kxpad = 2*np.pi*np.fft.fftfreq(2*Lx,d=dx)
#kxpad = (2*np.pi/(2*Lx))*np.arange(-2*Lx/2,2*Lx/2)
#kxpad = np.fft.fftshift(kxpad)

#kypad = (2*np.pi/(2*Ly))*np.arange(-2*Ly/2,2*Ly/2)
#kypad = np.fft.fftshift(kypad)

kypad = 2*np.pi*np.fft.fftfreq(2*Ly,d=dy)
kpad = np.zeros((2*Lx,2*Ly))


for j in range(2*Ly):
    for i in range(2*Lx):
        kpad[i,j] =  math.sqrt(kxpad[i]**2+kypad[j]**2)

kpad = np.fft.fftshift(kpad)

gpad = np.zeros((2*Lx,2*Ly))


gpad[:Lx,:Ly] = g # Filling main part of g in gpad

gpad[:Lx,Ly:] = g[:,-1::-1] # Filling the last 3 columns of gpad with g flipped

gpad[Lx:,:Ly] = g[-1::-1,:]# Filling the last 3 lines of gpad with g flipped

gpad[Lx:,Ly:] = g[-1::-1, -1::-1]# Filling the last 3 lines and last 3 columns of gpad with g flipped in line and column

rdFFT2D = np.zeros((Lx,Ly))

gpadhat = np.fft.fft2(gpad)
dgpadhat = -(kpad**2)*gpadhat #taking the derivative iwFFT(f)
rdpadFFT2D = np.real(np.fft.ifft2(dgpadhat))
rdFFT2D = rdpadFFT2D[:Lx,:Ly]

[

第一张图是原函数g(x,y)的图，第二张图是g的解析拉普拉斯算子，第三张图是里约热内卢的甜面包（笑），其实是拉普拉斯算子用快速傅立叶变换。我在这里做错了什么？

编辑 : 评论涟漪效应。克里斯，你的意思是下图中 set_zlimit 引起的涟漪效应？请记住，结果应该是 8。

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

Answer 1

填充不会改变边界条件：您通过复制函数进行填充，镜像，四次。该函数是对称的，因此镜像不会改变它。因此，您的填充只是将函数重复四次。通过 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, projection='3d')
ax.plot_surface(x[None,:], y[:,None], 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,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)))

# 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, projection='3d')
ax.plot_surface(x[None,:], y[:,None], g)
ax.set_zlim(0, 800)
ax = fig.add_subplot(122, projection='3d')
ax.plot_surface(x[None,:], y[:,None], lapg)
ax.set_zlim(0, 800)
pp.show()

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

使用 fft 的二阶导数

Second derivative using fft

fft

derivative

laplacian