fft - 使用 fft 的二阶导数
问题描述
所有,我正在尝试采用以下函数的拉普拉斯算子:
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。
解决方案
填充不会改变边界条件:您通过复制函数进行填充,镜像,四次。该函数是对称的,因此镜像不会改变它。因此,您的填充只是重复该功能四次。通过 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
只是按比例放置了这种噪声。
推荐阅读
- windows - 我应该如何更改项目的目标处理器架构?
- javascript - 如何在 javascript 中拆分一些文本正文行?
- regex - 使用 termux 和 regexp 在 git 上搜索历史日志和文件
- c# - 如何在没有正则表达式或任何外部方法的情况下摆脱双空格?
- python - 这段代码给了我一个 else: ^ IndentationError: unexpected unindent error
- ruby-on-rails - Rails 隐藏字段不保存参数
- java - java中的Apache Solr日期字段转换错误
- python - windows 10上的python无法升级虚拟环境中的模块
- rust - 如何设置货物构建的功能选项?
- javascript - ReactJs - 输入组件从函数返回时表现异常