Реализация OpenCL FFT — бессмысленные выходные данные — предположительно правильный алгоритм

У меня возникли некоторые проблемы с модулем кросс-корреляции сигналов FFT, который я создаю (использует теорему о круговой свертке и т. Д. И т. Д.). Я хотел бы просто подтвердить, что следующая схема будет гарантировать, что определенные уровни рекурсии вычислений бабочки БПФ завершены до того, как начнется следующий уровень, и что буферы, содержащие данные, полностью записаны в / выполнены. Таким образом, круговая корреляция / свертка включает в себя FFT, векторное внутреннее произведение, а затем IFFT.

Из-за этой схемы у меня нет ядер, которые упорядочивают данные в обратном порядке индексов. Ядро прямого БПФ генерирует БПФ с битовым обратным порядком, и после внутренних продуктов IFFT просто использует этот результат для вычисления решения естественного порядка.

Я должен упомянуть, что у меня есть несколько графических процессоров.

В любом случае, вот псевдокодовое представление того, что происходит для каждого FFT / IFFT (алгоритмы доступа / работы эквивалентны, кроме сопряженных тиддл-факторов, ядро ​​нормализации приходит позже:

    for numDevices:
data -> buffers
buffers -> kernel arguments

for fftRecursionLevels:
for numDevices:
recursionLevel -> kernel arguments
deviceCommandQueue -> enqueueNDRangeKernel
deviceCommandQueue ->  flush()

for numDevices:
deviceCommandQueue -> finish()

(EDIT: метод Radix-2 DIT, если это не ясно, извините.)

Могу ли я сойти с рук? Насколько я понимаю, finish () является блокирующей функцией, и этот самый последний цикл for не завершится, пока каждое ядро ​​не завершит вычисления во всем своем глобальном диапазоне (здесь fftSize / 2, см. Любую литературу по операциям с бабочками в Radix-2), и Что касается бонусных баллов, некоторые ядра уже выполняются из-за flush (), а я ставлю в очередь остальные ядра.

В целом, я получаю некоторые результаты с использованием openCL / c ++ для этого конкретного программного обеспечения. Я реализовал полный конвейер данных в Python, (алгоритм «топологически эквивалентен», если хотите, очевидно, что хоста нет<-> буфер устройства / инструкции или операции на стороне устройства с методом python) и моделирование работы ядра и это дает идентичные результаты, когда я использую модули scipy.fftpack и просто оперировать векторами данных сигнала.

Я думаю, что некоторые фотографии помогут. Вот именно то, что происходит в обеих программах.

1) создать гауссов вектор
2) Гауссов вектор нулевой площадки к следующей следующей по величине степени длины 2
3) вперед БПФ, в результате чего в естественном порядке (в ш) результат
4) сюжет

Вот симуляция Python моего ядра, по сравнению с использованием scipy.fftpack.fft (vector):

http://i.imgur.com/pGcYTrL.png

Они одинаковые. Теперь сравните это с одним из:

http://i.imgur.com/pbiYGpR.png

(Игнорировать значения на осях x, все они — результаты БПФ естественного порядка)

Все они относятся к одному и тому же типу начальных данных (в этом случае — от гассианского от 0 до N, с центром в N / 2 и с добавлением нуля до 2N). И все они должны выглядеть как зеленая / синяя линия на изображении один, но это не так. Мои глаза застеклены тем, как долго я смотрю на код хоста / устройства для этой второй программы, и я не вижу никаких опечаток или неправильных алгоритмов. Я очень подозреваю, что что-то происходит на стороне устройства, о котором я не знаю, поэтому моя публикация здесь. Ясно, что алгоритм выглядит так, как будто он работает правильно (общая форма красного / красного приблизительно равна синему / зеленому независимо от начальных данных. Я запускал алгоритм на разных начальных наборах, и он последовательно выглядит как синий / зеленый, но с этим бессмысленным шумом / ошибкой), но что-то не так.

Итак, я перехожу к межсетям. Заранее спасибо.

РЕДАКТИРОВАТЬ: один из постеров ниже предположил, что трудно комментировать, не видя хотя бы код со стороны устройства, так как есть вопросы по поводу фехтования, поэтому я разместил код ядра ниже.

//fftCorr.cl
//
//OpenCL Kernels/Methods for FFT Cross Correlation of Signals
//
//Copyright (C) 2013  Steve Novakov
//
//This file is part of OCLSIGPACK.
//
//OCLSIGPACK is free software: you can redistribute it and/or modify
//it under the terms of the GNU General Public License as published by
//the Free Software Foundation, either version 3 of the License, or
//(at your option) any later version.
//
//OCLSIGPACK is distributed in the hope that it will be useful,
//but WITHOUT ANY WARRANTY; without even the implied warranty of
//MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//GNU General Public License for more details.
//
//You should have received a copy of the GNU General Public License
//along with OCLSIGPACK.  If not, see <http://www.gnu.org/licenses/>.
//#define PIE 3.14159265359fvoid DFT2(float2 * a,float2 * b, float2 *w){
float2 tmp;
float2 bmul = ( (*w).x*((*b).x) - (*w).y*((*b).y),  (*w).x*((*b).y) + (*w).y*((*b).x) );
tmp = (*a) - bmul;
(*a) += bmul;
(*b) = tmp;
}
//
//
// Spin Factor Calc
//
// Computes spin/twiddle factor for particular bit reversed index.
//
//

float2 spinFact(unsigned int N, unsigned int k)
{
float phi = -2.0 * PIE * (float) k / (float)  N;
// \bar{w}^offset_(groupDist)

float spinRe, spinIm;
spinIm = sincos( phi, &spinRe);

return (float2) (spinRe, spinIm);
}float2 spinFactR(unsigned int N, unsigned int k)
{
float phi = 2.0 * PIE * (float) k / (float)  N;
// w^offset_(groupDist)

float spinRe, spinIm;
spinIm = sincos( phi, &spinRe);

return (float2) (spinRe, spinIm);
}//
// Bit-Reversed Index Reversal, (that sounds confusing)
//
unsigned int BRIR( unsigned int index, unsigned int fftDepth)
{
unsigned int rev = index;

rev = (((rev & 0xaaaaaaaa) >> 1 ) | ((rev & 0x55555555) << 1 ));
rev = (((rev & 0xcccccccc) >> 2 ) | ((rev & 0x33333333) << 2 ));
rev = (((rev & 0xf0f0f0f0) >> 4 ) | ((rev & 0x0f0f0f0f) << 4 ));
rev = (((rev & 0xff00ff00) >> 8 ) | ((rev & 0x00ff00ff) << 8 ));
rev = (((rev & 0xffff0000) >> 16) | ((rev & 0x0000ffff) << 16));

rev >>= (32-fftDepth);

return rev;
}//
//
// Index Bit Reversal Kernel, if Necessary/for Testing.
//
// Maybe I should figure out an in-place swap algorithm later.
//
//
__kernel void bitRevKernel(     __global float2 * fftSetX,
__global float2 * fftSetY,
__global float2 * fftRevX,
__global float2 * fftRevY,
unsigned int fftDepth
)
{
unsigned int glID = get_global_id(0);

unsigned int revID = BRIR(glID, fftDepth);

fftRevX[revID] = fftSetX[glID];
fftRevY[revID] = fftSetY[glID];

}//
//
// FFT Radix-2 Butterfly Operation Kernel
//
// This is an IN-PLACE algorithm. It calculates both bit-reversed indeces and spin factors in the same thread and
// updates the original set of data with the "butterfly" results.
//
// recursionLevel is the level of recursion of the butterfly operation
// # of threads is half the vector size N/2, (glID is between 0 and this value, non inclusive)
//
// Assumes natural order data input. Produces bit-reversed order FFT output.
//
//
__kernel void fftForwKernel(    __global float2 * fftSetX,
__global float2 * fftSetY,
unsigned int recursionLevel,
unsigned int totalDepth
)
{

unsigned int glID = get_global_id(0);

unsigned int gapSize = 1 << (recursionLevel - 1);
unsigned int groupSize = 1 << recursionLevel;
unsigned int base = (glID >> (recursionLevel - 1)) * groupSize;
unsigned int offset = glID & (gapSize - 1 );

unsigned int bitRevIdA = (unsigned int) base + offset;
unsigned int bitRevIdB = (unsigned int) bitRevIdA + gapSize;

unsigned int actualIdA = BRIR(bitRevIdA, totalDepth);
unsigned int actualIdB = BRIR(bitRevIdB, totalDepth);float2 tempXA = fftSetX[actualIdA];
float2 tempXB = fftSetX[actualIdB];

float2 tempYA = fftSetY[actualIdA];
float2 tempYB = fftSetY[actualIdB];

float2 spinF = spinFact(groupSize, offset);

// size 2 DFT
DFT2(&tempXA, &tempXB, &spinF);
DFT2(&tempYA, &tempYB, &spinF);fftSetX[actualIdA] = tempXA;
fftSetX[actualIdB] = tempXB;

fftSetY[actualIdA] = tempYA;
fftSetY[actualIdB] = tempYB;

}

Для данных, представленных на картинке. Я запускаю «fftForwKernel», как описано в начале поста, а затем запускаю «bitRevKernel»