NumPy 中 exp(-x^2) 的快速傅里叶变换答案

【问题标题】：Fast Fourier Transform of exp(-x^2) in NumPyNumPy 中 exp(-x^2) 的快速傅里叶变换
【发布时间】：2020-05-13 14:51:32
【问题描述】：

我必须用数值计算高斯函数的二阶导数：
我已经在这里阅读了有关此主题的所有问题，但无法得出一个好的结果。我选择 NumPy 作为我的首选工具。

我们教授的指示：

使用步骤dx = 1 获取大小为N = 128 的x 数组。所以，-64, -63, ..., 62, 63。计算f(x)
对f(x) 执行FFT 并接收转换后的数组f_m。
将f_m 乘以，其中是虚数单位，是导数和
执行逆 FFT 以接收导数。
在某些 FFT 实现中，您可能需要按 1/n 进行扩展（但这是目前最小的问题）

现在这是我的代码，尽可能简单。

import numpy as np

# Set some parameters
n = 128
dx = 1
a = 0.001

# Create x, calculate f(x) and its FFT
x = np.arange(-n/2, n/2) * dx
psi = np.exp(-a * x * x)
f_m = np.fft.fft(psi)

# k_m creation according to professor (point 3. in my instruction)
k_m = np.arange(-n/2, n/2, dtype=float)
k_m[:int(n / 2)] = (2 * np.pi * k_m[:int(n / 2)]) / (n * dx)
k_m[int(n / 2):] = (2 * np.pi * (k_m[int(n / 2):] - n)) / (n * dx)

# Multiply f_m by (j * k_m)^q. For q=2, this is -k_m^2
f_m *= -k_m * k_m
# Inverse FFT on the result to get the second derivative and scale by 1 / n
f_m = np.fft.ifft(f_m) / n

我无法得到的一件事是结果仍然有虚部，所以有些事情是不对的。有人可以帮忙吗？

编辑：Cris Luengo 的回答有效。

【问题讨论】：

标签： python numpy fft

【解决方案1】：

这部分错了：

k_m = np.arange(-n/2, n/2, dtype=float)

第 3 步中的说明谈到 m 从 0 到 n-1。代码应如下所示：

k_m = np.arange(0, n, dtype=float)
half = int(n / 2) + 1;  # notice the + 1 here!
k_m[:half] = (2 * np.pi * k_m[:half]) / (n * dx)
k_m[half:] = (2 * np.pi * (k_m[half:] - n)) / (n * dx)

FFT 产生一个输出，其中第一个元素（索引 0）是 0 频率，而不是频率 -n/2。

如果您使用 fftshift 将 0 频率区间移动到阵列的中间，您当前版本的 k_m 阵列可能是正确的，尽管我不完全确定（可能是 -n in后半部分应该去掉？）。

最后，这里不需要除以n：

f_m = np.fft.ifft(f_m) / n

NumPy IFFT 已经标准化。

并记得在验证虚部几乎为零之后绘制f_m.real（这些值应该与零不同，只是因为数字舍入错误）。

如果您将a 稍大一点，例如a=0.005，那么您的输入高斯将完全适合输入信号，并且您不会因过滤被截断的信号而产生难看的边缘效应。

【讨论】：

感谢您的回答。更改为f_m = np.fft.fftshift(np.fft.fft(psi)) 和k_m = np.arange(0, n, dtype=float) 后，它仍然无法正常工作......
@lkky7：是的，你应该做这两件事中的一个，而不是两者都做。
谢谢，这改变了一些东西，但结果仍然很糟糕。并且删除 -n 也不会改变这一点:( 添加结果进行编辑
@lkky7 当我使用我的更改运行您的代码时，我确实看到了正确的二阶导数。请注意，您正在将其与图表中的一阶导数进行比较。您还将结果按比例缩小了n，这是您不应该做的。结果确实有一个难看的边缘效果，如果你让a 更大，这个效果应该会减少。另请注意，您需要忽略（几乎为零）虚部，并绘制f_m.real。

【解决方案2】：

您可以使用更简单的 k，只要您在某个时间点执行正确的 FT 移位，其实现方式与您的讲师或 @CrisLuengo 明确编写的相同。

import numpy as np


# Set some parameters
n = 128
dx = 1
a = 0.001

# Create x, calculate f(x) and its FFT
x = np.arange(-n // 2, n // 2) * dx
f_x = np.exp(-a * x ** 2)
dd_f_x = 2 * a * np.exp(-a * x ** 2) * (2 * a * x ** 2 - 1)

f_k = np.fft.fft(f_x)
k = np.fft.ifftshift(np.arange(-n // 2, n // 2))
k = (2 * np.pi * k / (n * dx))
dd_f_k = -k ** 2 * f_k
dd_f_x_ = np.fft.ifft(dd_f_k)

按预期工作：

import matplotlib.pyplot as plt


fig, ax = plt.subplots(1, 1, squeeze=True)
ax.plot(x, dd_f_x_.real, label='∂²/∂x² f(x) with DFT')
ax.plot(x, dd_f_x, label='∂²/∂x² f(x)')
ax.legend()

【讨论】：