fft - 使用 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 的解析拉普拉斯算子，第三张图是里约热内卢的甜面包（lol），实际上是使用 FFT 的拉普拉斯算子。我在这里做错了什么？

编辑：评论涟漪效应。Cris 你的意思是由于下图中的 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只是按比例放置了这种噪声。