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

Question

我必须在数值上计算一个高斯函数的第2次导数。 我读了这里所有关于这个主题的问题， 但不能得出一个好的结果。我选择了NumPy作为我的工具。

我们教授的说明 :

1. 获取一个 x 大小阵列 N = 128 有步骤 dx = 1. 所以.., -64, -63, ..., 62, 63. 计算 f(x)
2. 进行FFT f(x) 并接收转换后的数组 f_m.
3. 乘以 f_m 由 ，其中 是虚单位。 是派生程度和
4. 执行逆FFT来接收导数。
5. 在某些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

有一点我无法理解，就是结果仍然有虚部，所以有些地方不对。谁能帮帮我？

EDIT：Cris Luengo的答案有效。

Answer 1

这个部分是错误的。

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

第三步的说明中说到 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.

您当前的版本是 k_m 数组可能是正确的，如果你使用 fftshift 将0频点移到数组的中间，虽然我并不完全确定（可能是将0频点移到数组的中间）。-n 在下半场应该删除）。)

---

最后，除以 n 这里没有必要。

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

NumPy的IFFT已经归一化了。

记住要绘制 f_m.real在验证了虚分量几乎为零之后（这些值应该与零相差无几，只是因为数值上的四舍五入误差）。

如果您使 a 大一点 a=0.005，那么你的输入高斯完全适合于输入信号，你不会因为过滤一个被切断的信号而产生丑陋的边缘效应。

Answer 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()