为什么`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包，或者它不准确？你知道吗

---

更新

以下是八度音阶部分脚本的等效形式：

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]上的误差在八度中几乎为零。（我没有检查其他示例，因为我不知道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或质疑它的准确性（因为这显然是一个错误的警报）。你知道吗