为啥不使用 Scipy 的 FFT 代码中的结果与 Scipy FFT 不相似？

Posted 2023-03-12

【中文标题】为啥不使用 Scipy 的 FFT 代码中的结果与 Scipy FFT 不相似？

【英文标题】：Why is the result in FFT code without using Scipy not similar to Scipy FFT?为什么不使用 Scipy 的 FFT 代码中的结果与 Scipy FFT 不相似？

【发布时间】：2020-03-17 12:44:36

【问题描述】：

下面的 FFT 代码没有给出类似于 Python 的scipy 库的结果。但我不知道这段代码有什么问题。

import numpy as np
import matplotlib.pyplot as plt
#from scipy.fftpack import fft

def omega(p, q):
   return np.exp((-2j * np.pi * p) / q)

def fft(x):
   N = len(x)
   if N <= 1: return x
   even = fft(x[0::2])
   odd = fft(x[1::2])
   combined = [0] * N
   for k in range(N//2):
     combined[k] = even[k] + omega(k,N) * odd[k]
     combined[k + N//2] = even[k] - omega(k,N) * odd[k]
   return combined

 N = 600
 T = 1.0 / 800.0
 x = np.linspace(0, N*T, N)
 #y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
 y = np.sin(50.0 * 2.0*np.pi*x)
 xf = np.linspace(0.0, 1.0/(2.0*T), N//2)
 yf = fft(y)
 yfa = 2.0/N * np.abs(yf[0:N//2])
 plt.plot(xf, yfa)
 plt.show()

这给出了：

【问题讨论】：

您已正确实施，但尚未完全优化。 FFT 基本上是函数的能量。所以你可能会得到一点点的失真。如果想看实际实现请参考github.com/scipy/scipy/blob/master/scipy/fftpack/basic.py

绘制 scipy 版本。差异应该很小。根据您对计算的排序和优化方式，舍入误差会在不同的地方蔓延。

【参考方案1】：

上述所有 cmets，ie 舍入误差和实现正确性都是正确的，但您错过了一个重要的事情...... FFT Cooley 和 Tukey 原始算法仅在样本数 N 时才有效是 2 的幂。你确实注意到了

np.allclose(yfa,yfa_sp)
>>> False

对于您当前的输入 N = 600，您的输出和 numpy/scipy 之间的差异很大。但是现在，让我们使用 2 的幂，在本例中为 N = 2**9 = 512，它给出了

np.allclose(yfa,yfa_sp)
>>> True

太棒了！这次的输出是相同的，并且可以验证输入信号y 的其他 2 次方（除了奈奎斯特准则）大小。对于深入的解释，您可以阅读this question 的已接受答案，以了解为什么 numpy/scipy fft 函数可能允许所有 N（当N 是 2 的幂时效率最高，而当 N 是prime) 而不是像你应该处理的那样处理这个错误，像

if np.log2(N) % 1 > 0:
    raise ValueError('size of input y must be a power of 2') 

正如 cmets 中所建议的，如果信号的大小不能那么容易地修改，零填充绝对是解决此类采样问题的方法。

希望这会有所帮助。

【讨论】：

@Yacoba：哦，非常感谢您的详细解释。这对我很有用。

...如果您的信号长度不正确，您可以简单地附加零，直到它是 2 的幂：一个域中的零填充会导致双域中的 sinc 插值。