I have to calculate numerically the 2nd derivative of a Gaussian function:
I've read every question on this topic here, but can't come to a good result. I've chosen NumPy as my tool of choice.
我们教授的指示：
x大小数组N = 128有步骤dx = 1. So, -64, -63, ..., 62, 63。计算f(x)
f(x)并接收转换后的数组f_m.f_m by 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 中的说明讲述了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 频率仓移到阵列的中间，尽管我不完全确定这一点（也许-n后半部分应该删除吗？）。
最后，除以n这里没有必要：
f_m = np.fft.ifft(f_m) / n
NumPy IFFT 已经标准化。
并记住情节f_m.real，在验证虚部几乎为零之后（这些值应仅由于数字舍入误差而不同于零）。
如果你做a稍微大一点，例如a=0.005，那么您的输入高斯完全适合输入信号，并且您不会因过滤被切断的信号而产生丑陋的边缘效应。