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

问题描述

我必须数值计算高斯函数的二阶导数：
我读了这里关于这个话题的每一个问题，但都没有得到好的结果。我已选择NumPy作为我的选择工具。

我们教授的指导：

1. 通过步骤dx = 1获取大小N = 128的x数组。所以，-64, -63, ..., 62, 63。计算f(x)
2. 对f(x)执行FFT，收到转换后的数组f_m。
乘以，其中是虚单位，是导数的程度，
3. 执行逆FFT以获得导数。
4. 在某些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

我不能得到的一件事是，结果仍然有虚部，所以有些东西是不正确的。有人能帮忙吗？

编辑：克里斯·卢恩戈的答案奏效。

---

解决方案

此部分错误：

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

步骤3中的说明介绍了从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？)。

---

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

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

NumPy IFFT已正常化。

并记住在验证虚分量几乎为零(这些值应仅因数字舍入误差而与零不同)后绘制f_m.real。

如果您将a放大一点，例如a=0.005，则您的输入高斯完全适合输入信号，并且您不会因为过滤被切断的信号而产生丑陋的边缘效果。