Почему `numpy.fft.irfft` так неточно? - PullRequest
Новая
       10

Почему `numpy.fft.irfft` так неточно?

3 голосов
/ 31 октября 2019

Я понимаю, что большинство подпрограмм FFT / IFFT имеют минимальный уровень ошибки. Я ожидал, что FFT NumPy будет иметь минимальный уровень ошибок в тех же порядках, что и FFTW (скажем, 1e-15), но следующий эксперимент показывает ошибки в порядке 1e-5.

Рассмотрите возможность вычисления IDFT для блока, хорошо известно , что в результате получается ядро ​​Dirichlet, похожее на sinc. Но это не то, что я получаю от numpy.fft.irfft. На самом деле, даже первый образец, который должен просто равняться ширине прямоугольника, деленной на количество точек БПФ, отключается на величину около 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)

Это выглядит какошибка имеет синусоидальную форму: enter image description here

Кто-нибудь испытывал такую ​​же проблему? Я что-то неправильно понимаю о пакете FFT NumPy или он просто не точен?

Обновление

Вот эквивалент части скрипта в 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] практически равна нулю в октаве. (Я не проверял другие сэмплы, потому что мне неизвестен эквивалент diric функции в октаве.)

1 Ответ

1 голос
/ 31 октября 2019

Благодаря 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 является точным. Вот график ошибки: enter image description here

Numpy правильно, моя математика была неверна. Прошу прощения за размещение этого вводящего в заблуждение вопроса. Я не знаю, какова стандартная процедура в этих случаях. Должен ли я удалить свой вопрос или оставить его здесь с этим ответом? Я просто не хочу, чтобы это подрывало NumPy или бросало вызов его точности (поскольку это явно было ложной тревогой).

...