Numpy fft 函数给出的输出不同于使用公式计算的 dft答案

【问题标题】：Numpy fft function giving output different from the dft calculated using formulaNumpy fft 函数给出的输出不同于使用公式计算的 dft
【发布时间】：2019-10-31 12:40:49
【问题描述】：

我正在尝试在 python 中实现 dft。我正在使用标准公式：

这是我的代码：

k = np.array([np.arange(-50, 50)])
fs, xn = wavfile.read('voice_recording.wav')
nbits = 16
max_nbits = float(2**(nbits-1))
xn = xn / (max_nbits + 1.0)
xn = np.expand_dims(xn[:,0], axis=1)
N = len(xn)
n = np.array([np.arange(0, N)])
Xk = np.sum(xn*np.exp(((-1j*2*math.pi)/N)*np.matmul(n.T, k)), axis=0)

这里，xn 是从 .wav 文件 (voice_recording.wav) 读取的音频信号。 FFT的代码是：

Xk1 = np.fft.fftshift(np.fft.fft(xn, n=100, axis=0))

但两个结果完全不同，即使它们应该相同。 DFT 图：

和 FFT 图：

我做错了什么？

【问题讨论】：

什么是xn？它是在哪里定义的？
@SeverinPappadeux 它是一个 numpy 数组。它是从 .wav 文件中读取的
你在一行代码中做了太多的事情，很难看到它做了什么，也很难调试。您应该编写代码，使其易于阅读。接下来，查看中间结果，确保每个结果都符合您的预期。这样你很快就会发现你的错误。
@CrisLuengo 并没有使用比 Python 循环更有效的矩阵计算。它使用 C 循环
我不建议添加循环。我建议写更短的陈述。一行中有太多的括号，无法轻松地分辨出对 exp 的调用中的内容是什么，什么不是！而且您无法检查np.matmul 的结果是否符合您的预期。 np.exp 的结果是您所期望的吗？现在你无法判断你是否逐行运行你的代码。将其拆分，以便您查看这些中间结果。

标签： python numpy signal-processing fft dft

【解决方案1】：

如果不下载您的数据文件，我假设它有超过 100 个样本。如果是的话，那么

np.fft.fft(xn, n=100, axis=0)

截断前 100 个样本并计算这些样本的 FFT。也就是说，它计算的内容与您的代码不同。

当我使用 xn = np.random.randn(100) 并运行您的代码时，Xk 和 Xk1 与 1e-13 左右相同。这表明您的代码是正确的。

要使用 FFT 算法仅计算 k 值的子集，请先计算完整的变换，然后丢弃不需要的值。例如：

Xk1 = np.fft.fft(xn, axis=0)
Xk1 = np.fft.fftshift(Xk1)
Xk1 = Xk1[(N//2 - 50):(N//2 + 50)]

【讨论】：

输出的维度应该等于k吗？当我们对 n 求和时？
是函数x[n]还是x[k]的输出？
@ShantanuShinde：您的代码可以计算任何k 的 DFT，它独立计算每个。 FFT 算法必须一次计算所有这些。它无法计算子集。您必须计算整个事物 (fft(xn, axis=0))，然后选择您想要查看的 100 个元素。
当我取 k 和 n 相等时，输出与 fft 输出匹配。那么您能否建议自己如何实现 fftshift？