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

2022-04-02 00:00:00 python numpy fft

问题描述

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

我们教授的指导:

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

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

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


解决方案

此部分错误:

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,则您的输入高斯完全适合输入信号,并且您不会因为过滤被切断的信号而产生丑陋的边缘效果。

相关文章