python - 为什么 `numpy.fft.irfft` 如此不精确？

Question

据我所知，大多数 FFT/IFFT 例程都有错误层。我原以为 NumPy 的 FFT 与 FFTW 的顺序相同(例如 1e-15)，但下面的实验显示错误的顺序为 1e-5 .

考虑计算一个盒子的 IDFT。是well-known结果是类似 sinc 的 Dirichlet 内核。但这不是我从 numpy.fft.irfft 中得到的。事实上，即使第一个样本应该简单地等于框的宽度除以 FFT 点数，也会偏离 4e-5 左右的数量，如下例所示:

import numpy as np
import matplotlib.pyplot as plt

from scipy.special import diric

N = 40960
K = 513

X = np.ones(K, dtype=np.complex)
x = np.fft.irfft(X, N)

print("x[0] = %g: expected %g - error = %g" % (x[0], (2*K+1)/N, x[0]-(2*K+1)/N))

# expected IDFT of a box is Dirichlet function (see
# https://en.wikipedia.org/wiki/Discrete_Fourier_transform#Some_discrete_Fourier_transform_pairs)

y = diric(2*np.pi*np.arange(N)/N, 2*K+1) * (2*K+1) / N

plt.figure()
plt.plot(x[:1024] - y[:1024])
plt.title('error')

plt.show(block=True)

看起来误差是正弦形式的:

有人遇到同样的问题吗？我是不是对 NumPy 的 FFT 包有什么误解，或者它只是不准确？

---

更新

以下是 Octave 中脚本的一部分:

N = 40960;
K = 513;

X = zeros(1, N);


X(1:K) = 1;
X(N-K:N) = 1;

x = ifft(X);

fprintf("x[0] = %g, expected = %g - error = %g\n", x(1), (2*K+1)/N, x(1)-(2*K+1)/N);

x[0] 上的错误在 Octave 中几乎为零。 (我没有检查其他示例，因为我不知道 Octave 中的 diric 函数的等价物。)

Answer 1

感谢MarkDickinson ，我意识到我的数学是错误的。正确的比较将通过以下方式进行:

import numpy as np
import matplotlib.pyplot as plt

from scipy.special import diric

N = 40960
K = 513

X = np.ones(K+1, dtype=np.complex)
x = np.fft.irfft(X, N)

print("x[0] = %g: expected %g - error = %g" % (x[0], (2*K+1)/N, x[0]-(2*K+1)/N))

# expected IDFT of a box is Dirichlet function (see
# https://en.wikipedia.org/wiki/Discrete_Fourier_transform#Some_discrete_Fourier_transform_pairs)

y = diric(2*np.pi*np.arange(N)/N, 2*K+1) * (2*K+1) / N

plt.figure()
plt.plot(x[:1024] - y[:1024])
plt.title('error')

plt.show(block=True)

这表明 irfft 是准确的。这是错误图:

Numpy 是正确的，我的数学不正确。我很抱歉发布这个误导性的问题。我不知道在这些情况下的标准程序是什么。我应该删除我的问题还是将其与此答案一起留在此处？我只是不希望它破坏 NumPy 或挑战其准确性(因为这显然是误报)。