为什么不使用Scipy的FFT代码的结果与Scipy FFT不相似？

w9apscun 于 2022-11-23 发布在其他

关注(0)|答案(1)|浏览(197)

下面的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()

这给出：

scipy

来源：https://stackoverflow.com/questions/58977386/why-is-the-result-in-fft-code-without-using-scipy-not-similar-to-scipy-fft

1条答案

按热度按时间

omhiaaxx1#

以上所有评论，* 即 * 舍入误差和实现正确性，都是正确的，但您忽略了一个重要问题... 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

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

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

或者甚至使用按位and运算符（真正优雅的测试IMO）：

if N & N-1:
    raise ValueError('size of input y must be a power of 2')

正如评论中所建议的，如果信号的大小不能很容易地修改，零填充绝对是解决这种采样问题的方法。
希望这对你有帮助。

赞(0）回复(0）举报 2022-11-23

我来回答

为什么不使用Scipy的FFT代码的结果与Scipy FFT不相似？

1条答案

相关问题

热门标签

最新问答