OpenCL FFT实现.无意义输出数据.大概正确的算法

Question

我正在创建的FFT信号互相关模块(利用循环卷积定理等)有一些问题。我只想确认以下方案将确保在下一级别开始之前完成FFT蝶形计算的特定递归级别，并确保包含数据的缓冲区被完全写入/完成。因此，循环相关/卷积涉及到一个FFT，一个向量的内积，然后一个IFFT。

由于这种方案，我没有内核来按位反转索引顺序排列数据。前向FFT内核产生一个有点倒序的FFT，并且在内部积之后，IFFT只是使用这个结果来计算一个自然的顺序解。

我应该说我有多个GPU。

无论如何，下面是每个FFT/IFFT正在发生的事情的伪代码表示(访问/操作算法是等价的，除了共轭的旋转因子外，规范化内核随后出现：

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

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

        for numDevices:
            deviceCommandQueue -> finish()

(编辑:方法是基数-2，如果不清楚，对不起。)

我能逃脱吗？据我所知，finish()是一个阻塞函数，在每个内核完成其全局范围内的计算之前，最后一个for循环是不可能完成的，(在这里，fftSize / 2，请参阅关于基-2蝶形运算的任何文献)，而且，为了获得额外的分数，在我为剩下的内核排队时，一些内核已经执行了。

总的来说，我使用openCL/c++对这个特定的软件使用openCL/c++获得了一些有趣的/垃圾结果。I在python中实现了完整的数据管道(算法是“拓扑等效的”，如果需要的话，显然没有主机<->设备缓冲区/指令或设备端操作w/ python方法)，并模拟了内核应该如何运行，并且它产生的结果与我使用scipy.fftpack模块而只对信号数据向量操作时的结果相同。

我想有些照片会有帮助的。以下是在这两个项目中正在发生的事情。

1)生成高斯矢量2)零pad高斯向量到下一最高的2长3)正向FFT，得到自然阶(W)结果4)图。

下面是我的内核的python模拟，与仅仅使用scipy.fftpack.fft(向量)相比：

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

他们是一样的。现在，将其与以下任一项进行比较：

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

(忽略x轴上的识别器，它们都是自然的FFT结果)

它们都是相同类型的起始数据，(从0到N的guassian，以N/2为中心，在本例中为零填充到2N )。它们都应该看起来像图像中的绿/蓝线，但它们不是。我的眼睛已经盯着第二个程序的主机/设备代码多久了，我没有看到任何打字或不正确的算法。我高度怀疑正在发生的设备方面的事情，我不知道，因此我在这里张贴。显然，该算法操作正确(无论起始数据如何，红色/红色的一般形状大约是蓝色/绿色。我在不同的启动集上运行了这个算法，它看起来一直像蓝色/绿色，但是带有那种毫无意义的噪音/错误)，但是有些地方不对劲。

所以我把话题转到网上。提前谢谢。

编辑:下面的一张海报建议，至少很难评论w/o看到的设备端代码，因为有关于击剑的问题，所以我在下面发布了内核代码。

//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”

Answer 1

因此，如果没有代码可以提供任何信息，并且假设代码是“相同的”，我倾向于说，假设所使用的算法在Python和OpenCL之间确实是相同的，那么这可能是一个同步问题。如果不是全局同步问题(在内核调用之间适当地分割工作)，那么在每次内核调用时，每个本地组都会出现缺少全局内存围栏甚至本地内存栅栏的问题。

你是怎么组织电话的？所提供的伪代码对于划分递归FFT究竟是如何进行的似乎不太清楚。我猜你这么做是为了.嗯，我甚至不确定你是在做DIT还是DIF，或者是FFT算法可以使用的大量数据流。据我所知，这可能是因为你做蝴蝶时没有适当的记忆，或者你可能严格地同步你的FFT步骤，以至于蝴蝶是与递归FFT步骤完全不同的内核调用的一部分，我的建议是完全无效的。

(我本想把这句话写成评论，但我没有这样做的声誉，所以我很抱歉，这是作为一个‘答案’张贴)