c++ - 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()

(編集: メソッドは Radix-2 DIT です。明確でない場合は申し訳ありません。)

私はそれで逃げることができますか？私が理解している限り、finish() はブロッキング関数であり、最後の for ループは、すべてのカーネルがそのグローバル範囲全体で計算を完了するまで完了しません (ここでは fftSize / 2、基数 2 のバタフライ演算に関する文献を参照してください)。、おまけとして、残りのカーネルをキューに入れている間に、flush() が原因で一部のカーネルが既に実行されています。

全体として、この特定のソフトウェアで openCL/c++ を使用すると、ファンキー/ガベージの結果が得られます。私は完全なデータパイプラインを pythonに実装しました（アルゴリズムは「トポロジ的に同等」です。明らかに、ホスト<-->デバイスバッファ/命令またはPythonメソッドを使用したデバイス側操作はありません）、シミュレートしましたカーネルの実行方法と、scipy.fftpack モジュールを使用して信号データベクトルを操作する場合と同じ結果が生成されます。

いくつかの写真が役立つと思います。これがまさに両方のプログラムで起こっていることです。

1) ガウスベクトルを生成します。 2) ガウスベクトルを 2 の次の最大べき乗までゼロパディングします。

これは、 scipy.fftpack.fft( vector ) を使用した場合と比較した、私のカーネルの python シミュレーションです。

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

それらは同じです。それを次のいずれかと比較します。

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

(x 軸のインデックスは無視してください。これらはすべて自然順序の FFT 結果です)

それらはすべて同じタイプの開始データです (0 から N までのガシアンで、N/2 を中心とし、この場合は 2N までゼロが埋め込まれます)。そして、それらはすべて画像 1 の緑/青の線のように見えるはずですが、そうではありません。その 2 番目のプログラムのホスト/デバイスコードを長い間見つめていたので、私の目は曇っていますが、タイプミスや間違ったアルゴリズムは見当たりません。私が気付いていないデバイス側で何かが起こっているのではないかと強く思っているので、ここに投稿します。明らかに、アルゴリズムは正しく動作しているように見えます (赤/赤の一般的な形状は、開始データに関係なく、ほぼ青/緑の形状です。異なる開始セットでアルゴリズムを実行しましたが、一貫して青/緑のように見えますが、その無意味なノイズ/エラーで）、しかし何かがおかしい.

だから私はインターウェブに目を向けます。前もって感謝します。

編集: 以下のポスターの 1 つは、mem fencing に関する質問があるため、少なくともデバイス側のコードを見ないとコメントするのは難しいと示唆していたので、以下にカーネルコードを投稿しました。

//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」を実行します

score 1 · Accepted Answer

したがって、何も通知するコードがなく、コードが実際に「同じ」であるという仮定の下で実行されている場合、使用されているアルゴリズムが Python と OpenCLの間で実際に同じように見えると仮定すると、おそらく同期の問題であると言う傾向があります。 . グローバル同期の問題 (カーネル呼び出し間で作業を適切に分割する) ではない場合、グローバルメモリフェンスの欠落、またはカーネル呼び出しごとのローカルグループごとのローカルメモリフェンスの問題です。

どのように通話を整理していますか? 提供されている疑似コードは、再帰的な FFT を分割している正確さに関してあいまいに見えます。私はあなたが「正しいこと」をしていると推測しています...ええと、あなたがDITやDIF、またはFFTアルゴリズムで利用できる膨大な種類のデータフローのいずれかを行っているかどうかさえわかりません。私が知っている限りでは、適切にメンフェンシングせずにバタフライを実行している可能性があります。または、FFT ステップを厳密に同期しているため、バタフライは再帰的な FFT ステップとはまったく異なるカーネル呼び出しの一部であり、私の提案は完全に無効になっている可能性があります。そして無効。

（私はこれをコメントに骨抜きにしていたでしょうが、そうする評判がないので、これが「回答」として投稿されたことをお詫びします）

c++ - OpenCL FFT 実装 - 無意味な出力データ - おそらく正しいアルゴリズム

1 に答える 1

Related

Reference