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

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

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

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

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

    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/С++ для этого конкретного программного обеспечения. Я реализовал полный конвейер данных в python (алгоритм «топологически эквивалентен», если хотите, очевидно, что нет хоста --> буфера/инструкций устройства или операций на стороне устройства с/ метод python) и смоделировал, как должно работать ядро, и оно дает ИДЕНТИЧНЫЕ результаты, когда я использую модули scipy.fftpack и просто оперирую векторами данных сигнала.

Я думаю, некоторые фотографии не помешали бы. Именно это и происходит в обеих программах.

1) сгенерировать гауссовский вектор 2) нулевой гауссовский вектор до следующей следующей наивысшей степени длины 2 3) прямое БПФ, что приводит к результату естественного порядка (в w ) 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.14159265359f


void 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».