Python中的FFT窗口选择与优化

Question

我现在正在尝试计算总谐波失真（THD）、噪声底和其他音频测量（如互调失真IMD、频率响应），使用的是Python。为了做到这一点，我先把波形文件导入到numpy数组中，然后用scipy模块计算快速傅里叶变换（FFT）。为了避免混叠，我需要在进行FFT之前对数据进行窗函数处理。所以我尝试比较不同的窗函数，这里是一些结果（997 kHz的正弦波，32位，192 kHz是用Adobe Audition生成的）：

我在寻找精确度：噪声底应该尽可能低，峰值之外的响应应该尽量平坦。那么我的问题是：Rife-Vincent窗函数真的是我最好的选择吗？我是否错过了其他我不知道且没有测试的“秘密”窗函数？

如果我决定使用Rife-Vincent窗函数，问题就是计算时间！其他窗函数已经在scipy模块中实现，计算速度非常快。我计算Rife-Vincent系数的方式是：

w = np.empty(M,dtype=np.float64)
a = 2*np.pi/M
for i in np.arange(0, M):
    w[i] = (35 - 56*np.cos(a*i) + 28*np.cos(2*a*i) - 8*np.cos(3*a*i) + np.cos(4*a*i))/128

其中M是我的数据长度，可能会很长。这非常耗时，有人能帮我优化一下吗？

Answer 1

关于如何计算窗口的问题：

a = 2*np.pi/M
x = np.arange(0, M)
w = (35 - 56 * np.cos(a * x) + 28 * np.cos(2 * a * x) - 8 * np.cos(3 * a * x) + np.cos(4 * a * x))/128

这个过程应该挺快的。你可以通过这样做来省去一些临时变量：

w = 35 - 56 * np.cos(a * x)
w +=  28 * np.cos(2 * a * x)
w -= 8 * np.cos(3 * a * x)
w += np.cos(4 * a * x)
w /= 128.

不过这样做大约只能节省7秒中的1秒，针对3000万个点来说。

如果你想要更快的速度，可以使用numexpr：

w = ne.evaluate('(35 - 56 * cos(a * x) + 28 * cos(2 * a * x) - 8 * cos(3 * a * x) + cos(4 * a * x))/128')

它会并行编译和计算；在我的电脑上，时间从7秒减少到2秒。