使用numpy进行FFT归一化

Question

我刚开始使用numpy包并从简单的任务开始计算输入信号的FFT。以下是代码：

import numpy as np
import matplotlib.pyplot as plt

#Some constants
L = 128
p = 2
X = 20
x = np.arange(-X/2,X/2,X/L)
fft_x = np.linspace(0,128,128, True)

fwhl = 1

fwhl_y = (2/fwhl) \
*(np.log([2])/np.pi)**0.5*np.e**(-(4*np.log([2]) \
*x**2)/fwhl**2)

fft_fwhl = np.fft.fft(fwhl_y, norm='ortho')

ampl_fft_fwhl = np.abs(fft_fwhl)

plt.bar(fft_x, ampl_fft_fwhl, width=.7, color='b')

plt.show()

因为我使用的是带有一些常数除以 pi 的指数函数，所以我期望在傅里叶空间中得到指数函数，其中 FFT 的常数部分始终等于 1（零频率）。但是，我使用 numpy 得到的该组分值较大（约为 1.13）。这里我有一个幅度谱，它被 1 /（number_of_counts）**0.5 归一化（这是我在 numpy 文档中读到的）。我不明白有什么问题... 能有人帮帮我吗？

谢谢！

[编辑] 看起来问题已经解决，你只需要将 FFT 乘以步长（在我的例子中是 X/L），就可以得到 Fourier 积分和 FFT 相同的结果。至于 numpy.fft.fft（...，norm ='ortho'）的标准化选项，它仅用于保存变换的比例，否则您将需要将反向 FFT 的结果除以样本数量。感谢每个人的帮助！

- GuntherOnFire

---

在进行除法运算时要小心 - 如果你在Python中对两个整数进行除法运算，它们不会返回浮点数。如果你将15除以10，你会得到1而不是1.5。同样地，20/128=0而不是0.15。 - flyingmeatball

@flyingmeatball 在 Python 2.x 上是正确的。在 Python 3.x 上，它们被显式地转换为浮点数 ;) - BPL

零频率对应于输入的平均值：fft_fwhl[0]＃(0.56568542494923801+0j)，它与np.mean(fwhl_y)* np.sqrt(len(fwhl_y))＃0.56568542494923812相匹配（这是因为您正在使用norm =“ortho”）。 - Ahmed Fasih

我觉得以DFT定义来思考FFT比连续傅里叶变换更有帮助：http://docs.scipy.org/doc/numpy/reference/routines.fft.html#implementation-details，有时让DFT/FFT得到与连续傅里叶变换相同的答案还是有一点难度。 - Ahmed Fasih

对于那些寻求更深入了解FFT缩放（和窗口）的人，我推荐Heinzel等人的这篇论文，这篇论文帮助我很多地理解了FFT缩放。它比大多数入门文本都更加应用导向，但仍具有足够的理论，可以超越简单的公式：离散傅立叶变换(DFT)的频谱与谱密度估计，包括全面的窗口函数列表和一些新的平顶窗，G. Heinzel等，2002 - cx05

Answer 1

我终于解决了我的问题。你需要将FFT与傅里叶积分结合起来，只需将变换（FFT）的结果乘以步长（在我的情况下为X/L，FFT X/L），它通常有效。在我的情况下，由于我有一个额外的函数转换规则，所以情况要复杂一些。我必须确保曲线下面积等于1，因为这是δ函数的模型，所以由于步长不可改变，我必须满足step sum（fwhl_y）= 1条件，即X/L = 1/sum（fwhl_y）。因此，为了获得正确的结果，我必须做以下几件事：

1. 计算FFT fft_fwhl = np.fft.fft(fwhl_y)
2. 消除相位组件，这是由于fwhl_y函数的对称性导致的，该函数定义在 [-T/2,T/2] 间隔中，其中T是周期，np.fft.fft操作认为我的函数定义在 [0，T] 间隔中。因此，为了仅获得幅度谱（这就是我需要的），我只需使用 np.abs（FFT）
3. 为了获得所期望的值，我应该将我在上一步获得的结果乘以X / L，即 np.abs（FFT）* X/L
4. 我有一个关于曲线下面积的额外条件，因此为 X/L * sum（fwhl_y）= 1 ，最终我得出了 np.abs（FFT）* X/L = np.abs（FFT）/sum（fwhl_y）

希望至少能对任何人有所帮助。

- GuntherOnFire

Answer 2

这是一个可能的解决方案：

import numpy as np
import matplotlib.pyplot as plt
from scipy import fft
from numpy import log, pi, e

# Signal setup
Fs = 150
Ts = 1.0 / Fs
t = np.arange(0, 1, Ts)
ff = 50
fwhl = 1
y = (2 / fwhl) * (log([2]) / pi)**0.5 * e**(-(4 * log([2]) * t**2) / fwhl**2)

# Plot original signal
plt.subplot(2, 1, 1)
plt.plot(t, y, 'k-')
plt.xlabel('time')
plt.ylabel('amplitude')

# Normalized FFT
plt.subplot(2, 1, 2)
n = len(y)
k = np.arange(n)
T = n / Fs
frq = k / T
freq = frq[range(n / 2)]

Y = np.fft.fft(y) / n
Y = Y[range(n / 2)]

plt.plot(freq, abs(Y), 'r-')
plt.xlabel('freq (Hz)')
plt.ylabel('|Y(freq)|')

plt.show()

当fwhl=1时：

当fwhl=0.1时：

从上面的图表中可以看出，当fwhl接近0时，指数和FFT图的变化情况。

- BPL

---

但是FFT的零分量仍然不同，尽管它应该始终保持相同。我的意思是，如果我解决积分以找到这个函数的FFT，我得到这个结果。所以无论参数的值是什么，我在零点上始终有1。而在你的情况下，FFT的零分量总是有不同的值。 - GuntherOnFire

@GuntherOnFire 哦，是的，DFT / FFT 定义在周期离散序列上。所有时间上的傅里叶积分与离散函数的FFT没有直接关系。 - Ahmed Fasih

@Ahmed Fasih，看起来我终于找到了如何将傅里叶积分与FFT / DFT连接起来的解决方案，你只需要将FFT的结果乘以您选择的步骤即可。在我的情况下，它是X / L。 - GuntherOnFire

@GuntherOnFire 不错！你能将它写成一个答案并接受吗？ - Ahmed Fasih

@Ahmed Fasih，终于做好了。希望它会有用。 - GuntherOnFire