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

Question

我知道大多数FFT/IFFT例程都有误差。我期望NumPy的FFT的误差范围与FFTW相同（例如1e-15），但以下实验显示误差的数量级为1e-5。

考虑计算一个矩形函数的IDFT。众所周知，结果是类似于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);

在Octave中，x [0]的误差几乎为零。（我没有检查其他样本，因为我不知道Octave中diric函数的等效物。）

- MikeL

---

你的NumPy和Octave代码中的X数组似乎不相等；NumPy数组中没有类似于Octave中的X(N-K:N) = 1;。 - user2357112

你的代码中似乎有一个参数在dirichlet函数中偏差了一位。 - Mad Physicist

@user2357112 这是因为Octave没有特定的“实数”FFT/IFFT例程。因此，我必须指定完整的输入到IFFT，因此我将对应于负频率的元素设置为1。 - MikeL

你确定不想要 X = np.ones(K+1, dtype=np.complex) 吗？（或者等价地，在期望结果中将 2*K + 1 替换为 2*K-1。） - Mark Dickinson

@MarkDickinson 谢谢！我的例子确实是错的。应用你提到的更正，我得到了正确的结果。 - MikeL

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的准确性或挑战它的准确性（因为这显然是一个虚假警报）。

- MikeL