2

バックグラウンド

OpenCL を使用して基数 2 の FFT (Stockham 自動並べ替え) 用に Microsoft Research からこのアルゴリズムを実装しました。

非整数ポイントでサンプリングする必要があり、それをテクスチャ サンプリング ハードウェアに委任する方がよいと考えたため、カーネルの入力と出力に浮動小数点テクスチャ (256 列 XN 行) を使用します。FFT は常に 256 ポイントのシーケンスであることに注意してください (テクスチャのすべての行)。この時点で、使用している GPU と許可されている最大 2D テクスチャ サイズに応じて、私の N は 16384 または 32768 です。

また、4 つの実数値シーケンスの FFT を一度に実行する必要があるため、カーネルは FFT(a、b、c、d) を FFT(a + ib、c + id) として実行し、そこから 4 つの複素数を抽出できます。 O（n）アルゴリズムを使用して後でシーケンスアウトします。誰かが望むなら、これについて詳しく説明できますが、この質問の範囲に収まるとは思いません。

カーネルソース

const sampler_t fftSampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST;

__kernel void FFT_Stockham(read_only image2d_t input, write_only image2d_t output, int fftSize, int size)
{
    int x = get_global_id(0);
    int y = get_global_id(1);
    int b = floor(x / convert_float(fftSize)) * (fftSize / 2);
    int offset = x % (fftSize / 2);
    int x0 = b + offset;
    int x1 = x0 + (size / 2);

    float4 val0 = read_imagef(input, fftSampler, (int2)(x0, y));
    float4 val1 = read_imagef(input, fftSampler, (int2)(x1, y));

    float angle = -6.283185f * (convert_float(x) / convert_float(fftSize));

    // TODO: Convert the two calculations below into lookups from a __constant buffer
    float tA = native_cos(angle);
    float tB = native_sin(angle);

    float4 coeffs1 = (float4)(tA, tB, tA, tB);
    float4 coeffs2 = (float4)(-tB, tA, -tB, tA);
    float4 result = val0 + coeffs1 * val1.xxzz + coeffs2 * val1.yyww;

    write_imagef(output, (int2)(x, y), result);
}

ホスト コードは、単純にこのカーネルを log2(256) 回呼び出し、入力テクスチャと出力テクスチャをピンポンします。

注: native_cosとnative_sinを削除して、タイミングに影響があるかどうかを確認しようとしましたが、状況はあまり変わらないようです。いずれにせよ、私が探している要因ではありません。

アクセス パターン 私はおそらくメモリ帯域幅の制約を受けているので、基数 2 の FFT のメモリ アクセス パターン (行ごと) を次に示します。

FFTアクセスパターン

  • X0 - 組み合わせる (読み取る) 要素 1
  • X1 - 組み合わせる (読み取る) 要素 2
  • X - 書き込む要素 (書き込み)

質問

だから私の質問は - 誰かが私を助けてくれたり、このアルゴリズムのより高い基数の定式化に向けて私を指摘したりできますか? ほとんどの FFT は、大規模なケースと単一の実数/複素数値シーケンス用に最適化されているためです。彼らのカーネルジェネレーターも大文字と小文字に大きく依存しており、内部をいじろうとするとすぐに壊れてしまいます。

単純に基数 8 または 16 のカーネルに移行するよりも優れたオプションはありますか?

私の制約のいくつかは次のとおりです-OpenCLを使用する必要があります（cuFFTなし）。この目的で ACML のclAmdFftを使用することもできません。CPU の最適化 (このカーネルは CPU で大きな時間を費やします) についても話すとよいでしょうが、GPU でより少ない反復で実行することが私の主なユースケースです。

これをすべて読んで、助けようとしてくれてありがとう！

4

1 に答える 1

1

いくつかのバージョンを試しましたが、CPU と GPU で最高のパフォーマンスを発揮したのは、私の特定のケースでは基数 16 のカーネルでした。

参照用のカーネルは次のとおりです。これは、Eric Bainville の (最も優れた) Web サイトから取得され、完全な帰属とともに使用されています。

// #define M_PI 3.14159265358979f

//Global size is x.Length/2, Scale = 1 for direct, 1/N to inverse (iFFT)
__kernel void ConjugateAndScale(__global float4* x, const float Scale)
{
   int i = get_global_id(0);

   float temp = Scale;
   float4 t = (float4)(temp, -temp, temp, -temp);

   x[i] *= t;
}


// Return a*EXP(-I*PI*1/2) = a*(-I)
float2 mul_p1q2(float2 a) { return (float2)(a.y,-a.x); }

// Return a^2
float2 sqr_1(float2 a)
{ return (float2)(a.x*a.x-a.y*a.y,2.0f*a.x*a.y); }

// Return the 2x DFT2 of the four complex numbers in A
// If A=(a,b,c,d) then return (a',b',c',d') where (a',c')=DFT2(a,c)
// and (b',d')=DFT2(b,d).
float8 dft2_4(float8 a) { return (float8)(a.lo+a.hi,a.lo-a.hi); }

// Return the DFT of 4 complex numbers in A
float8 dft4_4(float8 a)
{
  // 2x DFT2
  float8 x = dft2_4(a);
  // Shuffle, twiddle, and 2x DFT2
  return dft2_4((float8)(x.lo.lo,x.hi.lo,x.lo.hi,mul_p1q2(x.hi.hi)));
}

// Complex product, multiply vectors of complex numbers

#define MUL_RE(a,b) (a.even*b.even - a.odd*b.odd)
#define MUL_IM(a,b) (a.even*b.odd + a.odd*b.even)

float2 mul_1(float2 a, float2 b)
{ float2 x; x.even = MUL_RE(a,b); x.odd = MUL_IM(a,b); return x; }
float4 mul_1_F4(float4 a, float4 b)
{ float4 x; x.even = MUL_RE(a,b); x.odd = MUL_IM(a,b); return x; }


float4 mul_2(float4 a, float4 b)
{ float4 x; x.even = MUL_RE(a,b); x.odd = MUL_IM(a,b); return x; }

// Return the DFT2 of the two complex numbers in vector A
float4 dft2_2(float4 a) { return (float4)(a.lo+a.hi,a.lo-a.hi); }

// Return cos(alpha)+I*sin(alpha)  (3 variants)
float2 exp_alpha_1(float alpha)
{
  float cs,sn;
  // sn = sincos(alpha,&cs);  // sincos
  //cs = native_cos(alpha); sn = native_sin(alpha);  // native sin+cos
  cs = cos(alpha); sn = sin(alpha); // sin+cos
  return (float2)(cs,sn);
}
// Return cos(alpha)+I*sin(alpha)  (3 variants)
float4 exp_alpha_1_F4(float alpha)
{
  float cs,sn;
  // sn = sincos(alpha,&cs);  // sincos
  // cs = native_cos(alpha); sn = native_sin(alpha);  // native sin+cos
  cs = cos(alpha); sn = sin(alpha); // sin+cos
  return (float4)(cs,sn,cs,sn);
}


// mul_p*q*(a) returns a*EXP(-I*PI*P/Q)
#define mul_p0q1(a) (a)

#define mul_p0q2 mul_p0q1
//float2  mul_p1q2(float2 a) { return (float2)(a.y,-a.x); }

__constant float SQRT_1_2 = 0.707106781186548; // cos(Pi/4)
#define mul_p0q4 mul_p0q2
float2  mul_p1q4(float2 a) { return (float2)(SQRT_1_2)*(float2)(a.x+a.y,-a.x+a.y); }
#define mul_p2q4 mul_p1q2
float2  mul_p3q4(float2 a) { return (float2)(SQRT_1_2)*(float2)(-a.x+a.y,-a.x-a.y); }

__constant float COS_8 = 0.923879532511287; // cos(Pi/8)
__constant float SIN_8 = 0.382683432365089; // sin(Pi/8)
#define mul_p0q8 mul_p0q4
float2  mul_p1q8(float2 a) { return mul_1((float2)(COS_8,-SIN_8),a); }
#define mul_p2q8 mul_p1q4
float2  mul_p3q8(float2 a) { return mul_1((float2)(SIN_8,-COS_8),a); }
#define mul_p4q8 mul_p2q4
float2  mul_p5q8(float2 a) { return mul_1((float2)(-SIN_8,-COS_8),a); }
#define mul_p6q8 mul_p3q4
float2  mul_p7q8(float2 a) { return mul_1((float2)(-COS_8,-SIN_8),a); }

// Compute in-place DFT2 and twiddle
#define DFT2_TWIDDLE(a,b,t) { float2 tmp = t(a-b); a += b; b = tmp; }

// T = N/16 = number of threads.
// P is the length of input sub-sequences, 1,16,256,...,N/16.
__kernel void FFT_Radix16(__global const float4 * x, __global float4 * y, int pp)
{
  int p = pp;
  int t = get_global_size(0); // number of threads
  int i = get_global_id(0); // current thread


//////  y[i] = 2*x[i];
//////  return;

  int k = i & (p-1); // index in input sequence, in 0..P-1
  // Inputs indices are I+{0,..,15}*T
  x += i;
  // Output indices are J+{0,..,15}*P, where
  // J is I with four 0 bits inserted at bit log2(P)
  y += ((i-k)<<4) + k;

  // Load
  float4 u[16];
  for (int m=0;m<16;m++) u[m] = x[m*t];

  // Twiddle, twiddling factors are exp(_I*PI*{0,..,15}*K/4P)
  float alpha = -M_PI*(float)k/(float)(8*p);
  for (int m=1;m<16;m++) u[m] = mul_1_F4(exp_alpha_1_F4(m * alpha), u[m]);

  // 8x in-place DFT2 and twiddle (1)
  DFT2_TWIDDLE(u[0].lo,u[8].lo,mul_p0q8);
  DFT2_TWIDDLE(u[0].hi,u[8].hi,mul_p0q8);

  DFT2_TWIDDLE(u[1].lo,u[9].lo,mul_p1q8);
  DFT2_TWIDDLE(u[1].hi,u[9].hi,mul_p1q8);

  DFT2_TWIDDLE(u[2].lo,u[10].lo,mul_p2q8);
  DFT2_TWIDDLE(u[2].hi,u[10].hi,mul_p2q8);

  DFT2_TWIDDLE(u[3].lo,u[11].lo,mul_p3q8);
  DFT2_TWIDDLE(u[3].hi,u[11].hi,mul_p3q8);

  DFT2_TWIDDLE(u[4].lo,u[12].lo,mul_p4q8);
  DFT2_TWIDDLE(u[4].hi,u[12].hi,mul_p4q8);

  DFT2_TWIDDLE(u[5].lo,u[13].lo,mul_p5q8);
  DFT2_TWIDDLE(u[5].hi,u[13].hi,mul_p5q8);

  DFT2_TWIDDLE(u[6].lo,u[14].lo,mul_p6q8);
  DFT2_TWIDDLE(u[6].hi,u[14].hi,mul_p6q8);

  DFT2_TWIDDLE(u[7].lo,u[15].lo,mul_p7q8);
  DFT2_TWIDDLE(u[7].hi,u[15].hi,mul_p7q8);


  // 8x in-place DFT2 and twiddle (2)
  DFT2_TWIDDLE(u[0].lo,u[4].lo,mul_p0q4);
  DFT2_TWIDDLE(u[0].hi,u[4].hi,mul_p0q4);

  DFT2_TWIDDLE(u[1].lo,u[5].lo,mul_p1q4);
  DFT2_TWIDDLE(u[1].hi,u[5].hi,mul_p1q4);

  DFT2_TWIDDLE(u[2].lo,u[6].lo,mul_p2q4);
  DFT2_TWIDDLE(u[2].hi,u[6].hi,mul_p2q4);

  DFT2_TWIDDLE(u[3].lo,u[7].lo,mul_p3q4);
  DFT2_TWIDDLE(u[3].hi,u[7].hi,mul_p3q4);

  DFT2_TWIDDLE(u[8].lo,u[12].lo,mul_p0q4);
  DFT2_TWIDDLE(u[8].hi,u[12].hi,mul_p0q4);

  DFT2_TWIDDLE(u[9].lo,u[13].lo,mul_p1q4);
  DFT2_TWIDDLE(u[9].hi,u[13].hi,mul_p1q4);

  DFT2_TWIDDLE(u[10].lo,u[14].lo,mul_p2q4);
  DFT2_TWIDDLE(u[10].hi,u[14].hi,mul_p2q4);

  DFT2_TWIDDLE(u[11].lo,u[15].lo,mul_p3q4);
  DFT2_TWIDDLE(u[11].hi,u[15].hi,mul_p3q4);

  // 8x in-place DFT2 and twiddle (3)
  DFT2_TWIDDLE(u[0].lo,u[2].lo,mul_p0q2);
  DFT2_TWIDDLE(u[0].hi,u[2].hi,mul_p0q2);

  DFT2_TWIDDLE(u[1].lo,u[3].lo,mul_p1q2);
  DFT2_TWIDDLE(u[1].hi,u[3].hi,mul_p1q2);

  DFT2_TWIDDLE(u[4].lo,u[6].lo,mul_p0q2);
  DFT2_TWIDDLE(u[4].hi,u[6].hi,mul_p0q2);

  DFT2_TWIDDLE(u[5].lo,u[7].lo,mul_p1q2);
  DFT2_TWIDDLE(u[5].hi,u[7].hi,mul_p1q2);

  DFT2_TWIDDLE(u[8].lo,u[10].lo,mul_p0q2);
  DFT2_TWIDDLE(u[8].hi,u[10].hi,mul_p0q2);

  DFT2_TWIDDLE(u[9].lo,u[11].lo,mul_p1q2);
  DFT2_TWIDDLE(u[9].hi,u[11].hi,mul_p1q2);

  DFT2_TWIDDLE(u[12].lo,u[14].lo,mul_p0q2);
  DFT2_TWIDDLE(u[12].hi,u[14].hi,mul_p0q2);

  DFT2_TWIDDLE(u[13].lo,u[15].lo,mul_p1q2);
  DFT2_TWIDDLE(u[13].hi,u[15].hi,mul_p1q2);

  // 8x DFT2 and store (reverse binary permutation)
  y[0]    = u[0]  + u[1];
  y[p]    = u[8]  + u[9];
  y[2*p]  = u[4]  + u[5];
  y[3*p]  = u[12] + u[13];
  y[4*p]  = u[2]  + u[3];
  y[5*p]  = u[10] + u[11];
  y[6*p]  = u[6]  + u[7];
  y[7*p]  = u[14] + u[15];
  y[8*p]  = u[0]  - u[1];
  y[9*p]  = u[8]  - u[9];
  y[10*p] = u[4]  - u[5];
  y[11*p] = u[12] - u[13];
  y[12*p] = u[2]  - u[3];
  y[13*p] = u[10] - u[11];
  y[14*p] = u[6]  - u[7];
  y[15*p] = u[14] - u[15];
}

1 つではなく 2 つの複素数シーケンスの FFT を一度に実行するようにカーネルを変更したことに注意してください。また、はるかに大きなシーケンスで一度に 256 要素の FFT しか必要ないため、このカーネルを 2 回実行するだけで、より大きな配列に 256 の長さの DFT が残ります。

関連するホスト コードの一部を次に示します。

var ev = new[] { new Cl.Event() };
var pEv = new[] { new Cl.Event() };

int fftSize = 1;
int iter = 0;
int n = distributionSize >> 5;
while (fftSize <= n)
{
    Cl.SetKernelArg(fftKernel, 0, memA);
    Cl.SetKernelArg(fftKernel, 1, memB);
    Cl.SetKernelArg(fftKernel, 2, fftSize);

    Cl.EnqueueNDRangeKernel(commandQueue, fftKernel, 1, null, globalWorkgroupSize, localWorkgroupSize,
        (uint)(iter == 0 ? 0 : 1),
        iter == 0 ? null : pEv,
        out ev[0]).Check();
    if (iter > 0)
        pEv[0].Dispose();
    Swap(ref ev, ref pEv);

    Swap(ref memA, ref memB); // ping-pong

    fftSize = fftSize << 4;
    iter++;

    Cl.Finish(commandQueue);
}

Swap(ref memA, ref memB);

これが誰かを助けることを願っています!

于 2013-05-16T14:35:21.410 に答える