Фурье... разное в CUFFT, Mathematica, fftw3 Опции

[AlexCorvis]

Здравствуйте!
работал в пакете Mathematica, отлаживал свои алгоритмы... все работает отлично!
Начал кодить реализацию... для рассчета Фурье взял либу CUFFT, которая считает с использованием CUDA
Вобщем не сходятся спектры математики и CUFFT(
Попробовал fftw3 еще... спектр такой же как у CUFFT...

Фурье берется от массива действительных чисел размер которого 1024.
CUFFT имеет функцию которая из массива действительных чисел и выдает комплексные... я мнимую часть у них обнуляю и не учитываю ее. Т.к. мне нужно только реальная часть.
Математика после Abs[] отбрасывает мнимую часть.

Объясните, пожалуйста в чем разница и как считает фурье Mathematica... алгоритм работает только в математике из-за неправильных фурье
Вот так Mathematica

А вот так считает CUFFT и fftw3


Визуально они очень похожи, но значения разные... я думал что значения отличаются на какой-нибудь множитель, но это не так... некоторые отсчеты в два-три раза отличаются по значению, некоторые в 10-15!

[AlexCorvis]

В документации математики накопал что Фурье считается вот так:

Реализовал эту формулу в математике, но убрал мнимую часть по формуле Эйлера:


В коде на скорую руку получилось нечто такое:
В Математике

Код

inputPtr;

n = Length[inputPtr]

fftNew = ConstantArray[0, n];
For[s = 1, s <= n, s++, {
  For[r = 1, r <= n, r++, {
    fftNew[[s]] =
      fftNew[[s]] + inputPtr[[r]]*Cos[2*Pi*(r - 1)*(s - 1)/n];
    }]
  }]

Функция CUDA... вызывается для каждого элементов fftPtr параллельно. потоков 512, поэтому в каждом считаю два отсчета Фурье

Код

  __global__ void fft(float* inputPtr, float* fftPtr)
{
    //Получаем id текущей нити.
    int idx = threadIdx.x;
    int id1 = idx*2+0;
    int id2 = idx*2+1;


    fftPtr[id1] = 0;
    //Расчитываем результат.
    for(int i = 0; i < 1024; i++)
    {
        fftPtr[id1] += inputPtr[i]*cos((3.141592653*i*id1)/512);
        fftPtr[id2] += inputPtr[i]*cos((3.141592653*i*id2)/512);
    }

    fftPtr[id1] = abs(fftPtr[id1]/32);
    fftPtr[id2] = abs(fftPtr[id2]/32);

}

Вобщем моя реализация что в математике, что в С++ выдает почти одинаковый результат довольно похожий на CUFFT, но он опять же отличается от того что выдает стандартное Fourier[] Математики!

Сообщение отредактировал AlexCorvis - Jun 2 2010, 04:23

[AndriAno]

Возникла парочка вопросов?
А вы уверены что вам надо брать реальную часть а не модуль?(ABS если не ошибаюсь должен модуль возвращать?)

CUFFT считает с одинарной точностью если память не изменяет, Возможно ошибки округления.

[AlexCorvis]

AndriAno, Спасибо Вам большое!!!!
Я и забыл вообще что такое модуль комплексного... Всю ночь не спал, бился над этой проблемой!
После Вашей подсказки за две минуты накодил код и все сошлось!!!!
Счастью нет предела! Можно идти спать)))
Спасибо еще раз!!

[AndriAno]

Незачто.