skcuda.fft не то же самое, что numpy.fft.rfft? - PullRequest
Новая
       67

skcuda.fft не то же самое, что numpy.fft.rfft?

1 голос
/ 28 июня 2019

Я пытался проверить вывод fft против numpy fft для модульного тестирования, я понял вскоре после того, как он потерпел неудачу, это было не потому, что я сделал что-то не так, но skcuda буквально не дает тот же ответ , Я знал, что они будут немного отличаться, но по крайней мере одно из чисел на несколько порядков отличается от того, что производит numpy, и оба allclose и almost_equal возвращают огромные ошибки (33% и 25% для rtol=1e-6, 16% для atol=1e-6). Что я здесь не так делаю? Можно это исправить?

Тестовый файл: 

import pycuda.autoinit
from skcuda import fft
from pycuda import gpuarray
import numpy as np

def test_skcuda():
    array_0 = np.array([[1, 2, 3, 4, 5, 4, 3, 2, 1, 0]], dtype=np.float32)
    array_1 = array_0 * 10
    time_domain_signal = np.array([array_0[0], array_1[0]], dtype=np.float32)
    fft_point_count = 10
    fft_plan = fft.Plan(fft_point_count, np.float32, np.complex64,
                        batch=2)
    fft_reserved = gpuarray.empty((2, fft_point_count // 2 + 1), dtype=np.complex64)
    fft.fft(gpuarray.to_gpu(time_domain_signal), fft_reserved, fft_plan)

    np.testing.assert_array_almost_equal(
        np.fft.rfft(time_domain_signal, fft_point_count), fft_reserved.get())

test_skcuda()

Ошибка подтверждения: 

AssertionError: 
Arrays are not almost equal to 6 decimals

(mismatch 25.0%)
 x: array([[ 2.500000e+01+0.000000e+00j, -8.472136e+00-6.155367e+00j,
        -1.193490e-15+2.331468e-15j,  4.721360e-01-1.453085e+00j,
         2.664535e-15+0.000000e+00j,  1.000000e+00+0.000000e+00j],...
 y: array([[ 2.500000e+01+0.000000e+00j, -8.472136e+00-6.155367e+00j,
         8.940697e-08+5.960464e-08j,  4.721359e-01-1.453085e+00j,
         0.000000e+00+0.000000e+00j,  1.000000e+00+0.000000e+00j],...

вывод на печать: 

#numpy
[[ 2.50000000e+01+0.00000000e+00j -8.47213595e+00-6.15536707e+00j
  -1.19348975e-15+2.33146835e-15j  4.72135955e-01-1.45308506e+00j
   2.66453526e-15+0.00000000e+00j  1.00000000e+00+0.00000000e+00j]
 [ 2.50000000e+02+0.00000000e+00j -8.47213595e+01-6.15536707e+01j
  -1.11022302e-14+2.39808173e-14j  4.72135955e+00-1.45308506e+01j
   3.55271368e-14+7.10542736e-15j  1.00000000e+01+0.00000000e+00j]]

#skcuda
[[ 2.5000000e+01+0.0000000e+00j -8.4721355e+00-6.1553669e+00j
   8.9406967e-08+5.9604645e-08j  4.7213593e-01-1.4530852e+00j
   0.0000000e+00+0.0000000e+00j  1.0000000e+00+0.0000000e+00j]
 [ 2.5000000e+02+0.0000000e+00j -8.4721359e+01-6.1553673e+01j
   1.4305115e-06-4.7683716e-07j  4.7213597e+00-1.4530851e+01j
   0.0000000e+00+1.9073486e-06j  1.0000000e+01+0.0000000e+00j]]

Ответы [ 3 ]

1 голос
/ 28 июня 2019

это очень похоже на ошибки округления, числа с плавающей точкой одинарной точности имеют ~ 8 десятичных цифр точности (двойные числа имеют ~ 16)

вместо использования numpy.fft альтернативой будет использование fftpack от scipy , который поддерживает поплавки одинарной точности, например: 

from scipy import fftpack

x = np.array([1, 2, 3, 4, 5, 4, 3, 2, 1, 0])

y = fftpack.fft(
    np.array([x, x * 10], dtype=np.float32)
)
print(y[:,:6])

вывод: 

[[ 2.5000000e+01+0.0000000e+00j -8.4721355e+00-6.1553669e+00j
   8.9406967e-08+5.9604645e-08j  4.7213593e-01-1.4530852e+00j
   0.0000000e+00+0.0000000e+00j  1.0000000e+00+0.0000000e+00j]
 [ 2.5000000e+02+0.0000000e+00j -8.4721359e+01-6.1553673e+01j
   1.1920929e-06+1.9073486e-06j  4.7213583e+00-1.4530851e+01j
   0.0000000e+00+1.9073486e-06j  1.0000000e+01+0.0000000e+00j]]

, что выглядит намного ближе

1 голос
/ 28 июня 2019

На выходе БПФ имеется ошибка, которая связана с величиной входных значений. Каждый выходной элемент вычисляется на основе объединения всех входных элементов, и, следовательно, именно их величины определяют точность результата.

Вы вычисляете два одномерных БПФ в одном массиве. Каждый из них имеет разные входные величины и, следовательно, должен иметь разные допуски.

Следующий быстрый код демонстрирует, как вы могли бы реализовать это. Я не знаю, как настроить любую из функций в numpy.testing, чтобы сделать это. 

import numpy as np

array_0 = np.array([[1, 2, 3, 4, 5, 4, 3, 2, 1, 0]], dtype=np.float32)
array_1 = array_0 * 10
time_domain_signal = np.array([array_0[0], array_1[0]], dtype=np.float32)

# numpy result
a=np.array([[ 2.50000000e+01+0.00000000e+00j, -8.47213595e+00-6.15536707e+00j,
  -1.19348975e-15+2.33146835e-15j,  4.72135955e-01-1.45308506e+00j,
   2.66453526e-15+0.00000000e+00j,  1.00000000e+00+0.00000000e+00j],
 [ 2.50000000e+02+0.00000000e+00j, -8.47213595e+01-6.15536707e+01j,
  -1.11022302e-14+2.39808173e-14j,  4.72135955e+00-1.45308506e+01j,
   3.55271368e-14+7.10542736e-15j,  1.00000000e+01+0.00000000e+00j]])

# skcuda result
b=np.array([[ 2.5000000e+01+0.0000000e+00j, -8.4721355e+00-6.1553669e+00j,
   8.9406967e-08+5.9604645e-08j,  4.7213593e-01-1.4530852e+00j,
   0.0000000e+00+0.0000000e+00j,  1.0000000e+00+0.0000000e+00j],
 [ 2.5000000e+02+0.0000000e+00j, -8.4721359e+01-6.1553673e+01j,
   1.4305115e-06-4.7683716e-07j,  4.7213597e+00-1.4530851e+01j,
   0.0000000e+00+1.9073486e-06j,  1.0000000e+01+0.0000000e+00j]])

# Tolerance for result array row relative to the mean absolute input values
# 1e-6 because we're using single-precision floats
tol = np.mean(np.abs(time_domain_signal), axis=1) * 1e-6

# Compute absolute difference and compare that to our tolearances
diff = np.abs(a-b)
if np.any(diff > tol[:,None]):
   print('ERROR!!!')
0 голосов
/ 29 июня 2019

Крошечные результирующие значения e-8 (с плавающей запятой) и e-15 (двойные) из БПФ (при полномасштабном вводе или выводе где-либо около 1,0) по существу равны нулю (плюс числовой шум округления).

и

ноль + шум == ноль + шум

так что ваши результаты могут фактически совпадать.

...