这个傅里叶变换有什么问题？（在蟒蛇中）答案

【问题标题】：What is wrong with this Fourier transform? (in python)这个傅里叶变换有什么问题？（在蟒蛇中）
【发布时间】：2019-07-02 19:10:34
【问题描述】：

正如阅读我之前在此站点上的帖子的人所知道的那样，我一直在尝试实现一个在 Python 中使用 FFT 的 PDE 求解器。编程部分基本搞定了，但是程序产生了（非常适合本站的）溢出错误（基本上会增长很多，直到变成NaN）。

在排除了所有其他可能性后，我将问题归结为 FFT 以及我尝试进行导数的方式，因此我决定测试两个不同的 FFT（numpy 的 fft 模块和 @ 987654332@包），代码如下：

import pyfftw
import numpy as np
import matplotlib.pyplot as plt


def fftw_(y: np.ndarray) -> np.ndarray:
    a = pyfftw.empty_aligned((N, N), dtype=np.float64)
    b = pyfftw.empty_aligned((N, N//2+1), dtype=np.complex128)
    fft_object = pyfftw.FFTW(a, b, axes=(0, 1), direction='FFTW_FORWARD', flags=('FFTW_MEASURE',), threads=12)
    y_hat = fft_object(y)
    return y_hat


def ifftw_(y_hat: np.ndarray) -> np.ndarray:
    a = pyfftw.empty_aligned((N, N//2+1), dtype=np.complex128)
    b = pyfftw.empty_aligned((N, N), dtype=np.float64)
    fft_object = pyfftw.FFTW(a, b, axes=(0, 1), direction='FFTW_BACKWARD', flags=('FFTW_MEASURE',), threads=12)
    y = fft_object(y_hat)
    return y


def func(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    return np.exp(x)*np.sin(y)

dx = 0.02

x = np.arange(-1, 1, dx)
y = np.arange(-1, 1, dx)
X, Y = np.meshgrid(x, y)

N = len(x)

kxw, kyw = np.meshgrid(np.fft.rfftfreq(N, dx), np.fft.fftfreq(N, dx))
Lapw = -4*np.pi**2*(kxw**2+kyw**2)

kxnp, kynp = np.meshgrid(np.fft.fftfreq(N, dx), np.fft.fftfreq(N, dx))
Lapnp = -4*np.pi**2*(kxnp**2+kynp**2)



z = func(X, Y)

lap_z_w = ifftw_(Lapw*fftw_(z))
lap_z_np = np.fft.ifft2(Lapnp*np.fft.fft2(z))

lap_z_np_mag = np.abs(lap_z_np)
lap_z_np_ang = np.angle(lap_z_np)

plt.imshow(z, cmap='plasma')
plt.colorbar()
plt.savefig("f.png", dpi=200)
plt.clf()

plt.imshow(lap_z_w, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_fftw.png", dpi=200)
plt.clf()

plt.imshow(lap_z_np_mag, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_np_mag.png", dpi=200)
plt.clf()

plt.imshow(lap_z_np_ang, cmap='plasma')
plt.colorbar()
plt.savefig("Lap_np_ang.png", dpi=200)
plt.clf()

这里np.ndarray 的命名为Lapw 和Lapnp 是我认为 应该做离散拉普拉斯算子。而我选择的函数eˣsin(y)是调和函数，所以它的拉普拉斯应该为零。

但是程序的结果与这个预期值相差甚远。特别是我得到：

原函数f

使用 pyFFTW 的 f 的“拉普拉斯算子”

使用 Numpy 的 f 的“拉普拉斯算子”的幅度和相位

查看这些图的值（请注意颜色栏中的范围以及 20000 不是 0 的任何体面近似值这一事实）清楚地说明了为什么我制作的程序会溢出，但我没有不知道如何纠正。任何帮助将不胜感激。

【问题讨论】：

情节对我来说似乎很好。你期望一个零拉普拉斯算子。你得到一个零拉普拉斯算子。角度由浮点舍入误差给出（复数值 + i* 的角度当然是任意的）。您还会有一些边缘效应，这在 DFT 中是预期的。您需要意识到您的数值计算会产生您希望得到的结果的近似值。
不为零。侧面的颜色条是自动生成的，如您所见，它会上升到 20000-25000，具体取决于所采用的变换。这是不正确的。
如果误差是 1/100 甚至 1/10，我会说这是一个很好的近似值。即使函数的一阶偏导数也应该限制在 2-4 左右，到 20000 意味着必须有一个错误
那些是边缘效应。获取输入数据的多个副本，假装它们是瓦片，然后将瓦片彼此相邻放置。您会在两个图块之间看到非常尖锐的过渡。基于 DFT 的滤波对这种转变作出反应。 (1) 忽略域边缘附近的输出，或 (2) 对输入数据应用窗口函数，然后也忽略域边缘附近的输出。你不能期望边缘附近有有意义的结果，那里没有数据可以计算导数。
是的，如果您通过傅立叶域进行过滤，则您的过滤器根据定义具有无限支持。您的所有输出都受到边缘效应的影响。只有靠近边缘的少数样本很明显。如果您以不同方式缩放数据，您可能会注意到更多样本受到影响。如果效果太强，请使用窗口功能。或者改为对导数使用有限差分逼近。

标签： python numpy fft pyfftw

【解决方案1】：

这里的错误是您对函数的假设不正确。 e^x sin(y) 可能看起来很谐波，但您只计算了 -1 fft 将隐式地定期继续它，即您在函数的所有边缘处都会出现不连续性。如果函数不连续，则它不是谐波，尤其是在傅里叶变换中会出现发散。这就是使您的 FFT 在边缘发散的原因。除此之外，“远离”边缘，结果看起来与预期一样。

【讨论】：