计算科学 - 在 numpy 中为 DFT 指定网格间距 - 吾爱随笔录

在 numpy 中为 DFT 指定网格间距

计算科学 Python 麻木的傅立叶分析傅里叶变换

2021-12-13 02:01:01

我在 Python 3.7.2 中测试 numpy 1.16.1 的.fft 包。特别是，我试图验证转换是否类似于分析转换：

f (x) = e x p [- {(\frac{x - 5}{2})}^{2}]

$f(x) = \mathrm{exp}\left[-\left(\frac{x-5}{2}\right)^{2}\right]$

我从Wolfram Alpha那里得到 $\hat{f} = \mathcal{F}[f]$ 看起来像这样：

然后我尝试用 numpy 和 matplotlib 复制这个图，代码如下：

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(0, 10, 1/1000)
y = np.exp(-((x-5)**2)/4)

y_hat = np.fft.fftshift(np.fft.fft(y))
re_y_hat = np.real(y_hat)
im_y_hat = np.imag(y_hat)

fig, ax = plt.subplots()
ax.plot(x, re_y_hat, "b-", x, im_y_hat, "r-")
plt.show()
plt.close()

但我获得的图像与 Wolfram 给出的图像大不相同：

在最后一张图像中，零频率通过使用偏移到中心，np.fft.fftshift()因此尖峰对应于频率零。

我已经发现问题np.fft.fft()在于 $\Delta x$ 被指定，所以 numpy 解释的是我的数据变化非常缓慢，几乎恒定 $^{1}$ ，因此变换接近于常数函数的变换。

我查看了 numpy 文档和其他 SE 帖子以了解如何解决此问题，但一无所获。有谁知道如何解决这一问题？

$^{1}$ 我们可以计算 numpy看到的函数的平均斜率 $\frac{\mathrm{max}\{f\}-\mathrm{min}\{f\}}{x_{f_{\mathrm{max}}}-x_{f_{\mathrm{min}}}} = \frac{f(5)-f(0)}{n\Delta x} \approx \frac{1}{n\Delta x}$ 在哪里 $n$ 是将最大值与最小值分开的节点数。在这种情况下，由于 numpy 需要 $\Delta x = 1$ 默认情况下，斜率约为 1/5000=0.0002

1个回答

我很想说我完全理解我在下面所做的所有比例因子和变化，但我主要是玩因子直到事情匹配:) 我会等到我想通了，但我想把这个解决方案给你，所以请原谅我不完整的答案。

import numpy as np
import matplotlib.pyplot as plt

fs=1e2
t0=-100
t1=100.0
x = np.arange(t0,t1, 1./fs)
y = np.exp(-((x-5)**2)/4)

y_hat = np.fft.fftshift(np.fft.fft(np.fft.fftshift(y))) / x.shape[0] *(t1-t0)/np.sqrt(2.0*np.pi)
freq=np.arange(-fs/2,fs/2,1.0/(t1-t0))
omega=2*np.pi*freq
y_hat_exact=np.sqrt(2.0)*np.exp(-omega**2-5j*omega)

plt.ion()

ff=plt.figure(1)
ff.clf()
ax=ff.gca()
ax.plot(x,y,'b-')
ax.set_xlabel('x')
ax.set_ylabel('y(x)')

ff=plt.figure(2)
ff.clf()
ax=ff.gca()
ax.plot(omega, y_hat.real, "b.-",label='re fft')
ax.plot(omega, y_hat.imag, "r.-",label='im fft')
ax.plot(omega,y_hat_exact.real,'k--',label='re exact')
ax.plot(omega,y_hat_exact.imag,'k:',label='im exact')
ax.set_xlabel('Radial frequency $\omega$ (rad/s)')
ax.set_ylabel('F[y]')
ax.set_xlim([-1.0,1.0])

ax.legend()

plt.show()

一些评论：

您不能将x其用作绘制 FFT 的自变量，您需要频率（或径向频率）。该代码显示了如何计算它。
您必须将结果除以样本数。这是 FFT 的约定。
除以 $\sqrt{2\pi}$ 是匹配 Wolfram Alpha 使用的约定傅里叶变换约定
总时间的比例是通过t1-t0实验计算出来的。如果我考虑得够久，我确信这与解析傅里叶变换在这段时间内是一个积分这一事实有关。
FFT 的输入需要从 x=0 开始。如果负 x 中有非零分量，则需要将它们缠绕到最后，fftshift因此y
您的 Wolfram Alpha 链接和相应的图像适用于 $f(x) = e^{-(x-2)^2/4}$ 当您的代码用于 $f(x) = e^{-(x-5)^2/4}$ . 我的示例使用后者。
Wolfram Alpha 使用 $+i\omega t$ 约定，而此 Pythonfft函数使用 $-i\omega t$ 约定，这就是为什么与 Wolfram 相比，我的情节中的虚构部分被翻转了。

其它你可能感兴趣的问题

上一篇最终边界条件未知，Backward Euler 是否仍然有效执行？下一篇使用 fmincon 进行参数估计