加速 FFTW 修剪以避免大量零填充

Question

假设我有一个长度为 K * N 的序列 x(n)，并且只有前 N 个元素不为零。我假设 N << K，例如 N = 10 和 K = 100000。我想通过 FFTW 计算这样一个序列的 FFT。这相当于拥有一个长度为 N 的序列，并有一个 K * N 的零填充。由于 N 和 K 可能是 "large"，我有一个重要的零填充。我正在探索是否可以节省一些计算时间避免显式零填充。

案例K = 2

让我们从案例K = 2开始。在这种情况下，x(n)的DFT可以写成

若k为偶数，即k = 2 * m，则

这意味着DFT的这些值可以通过长度为N的序列的FFT来计算，而不是K * N。

若k为奇数，即k = 2 * m + 1，则

这意味着 DFT 的这些值可以通过长度为 N 的序列的 FFT 再次计算，而不是 K * N.

因此，总而言之，我可以将长度为 2 * N 的单个 FFT 与长度为 N.

的 2 个 FFT 进行交换

大小写随意K

在这种情况下，我们有

在写作 k = m * K + t 时，我们有

因此，总而言之，我可以将长度为 K * N 的单个 FFT 与长度为 N 的 K 个 FFT 进行交换。由于 FFTW 有 fftw_plan_many_dft，我可以期望在单个 FFT 的情况下有所收获。

为了验证这一点，我设置了以下代码

#include <stdio.h>
#include <stdlib.h>     /* srand, rand */
#include <time.h>       /* time */
#include <math.h>
#include <fstream>

#include <fftw3.h>

#include "TimingCPU.h"

#define PI_d            3.141592653589793

void main() {

    const int N = 10;
    const int K = 100000;

    fftw_plan plan_zp;

    fftw_complex *h_x = (fftw_complex *)malloc(N     * sizeof(fftw_complex));
    fftw_complex *h_xzp = (fftw_complex *)calloc(N * K, sizeof(fftw_complex));
    fftw_complex *h_xpruning = (fftw_complex *)malloc(N * K * sizeof(fftw_complex));
    fftw_complex *h_xhatpruning = (fftw_complex *)malloc(N * K * sizeof(fftw_complex));
    fftw_complex *h_xhatpruning_temp = (fftw_complex *)malloc(N * K * sizeof(fftw_complex));
    fftw_complex *h_xhat = (fftw_complex *)malloc(N * K * sizeof(fftw_complex));

    // --- Random number generation of the data sequence
    srand(time(NULL));
    for (int k = 0; k < N; k++) {
        h_x[k][0] = (double)rand() / (double)RAND_MAX;
        h_x[k][1] = (double)rand() / (double)RAND_MAX;
    }

    memcpy(h_xzp, h_x, N * sizeof(fftw_complex));

    plan_zp = fftw_plan_dft_1d(N * K, h_xzp, h_xhat, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_plan plan_pruning = fftw_plan_many_dft(1, &N, K, h_xpruning, NULL, 1, N, h_xhatpruning_temp, NULL, 1, N, FFTW_FORWARD, FFTW_ESTIMATE);

    TimingCPU timerCPU;
    timerCPU.StartCounter();
    fftw_execute(plan_zp);
    printf("Stadard %f\n", timerCPU.GetCounter());

    timerCPU.StartCounter();
    double factor = -2. * PI_d / (K * N);
    for (int k = 0; k < K; k++) {
        double arg1 = factor * k;
        for (int n = 0; n < N; n++) {
            double arg = arg1 * n;
            double cosarg = cos(arg);
            double sinarg = sin(arg);
            h_xpruning[k * N + n][0] = h_x[n][0] * cosarg - h_x[n][1] * sinarg;
            h_xpruning[k * N + n][1] = h_x[n][0] * sinarg + h_x[n][1] * cosarg;
        }
    }
    printf("Optimized first step %f\n", timerCPU.GetCounter());

    timerCPU.StartCounter();
    fftw_execute(plan_pruning);
    printf("Optimized second step %f\n", timerCPU.GetCounter());
    timerCPU.StartCounter();
    for (int k = 0; k < K; k++) {
        for (int p = 0; p < N; p++) {
            h_xhatpruning[p * K + k][0] = h_xhatpruning_temp[p + k * N][0];
            h_xhatpruning[p * K + k][1] = h_xhatpruning_temp[p + k * N][1];
        }
    }
    printf("Optimized third step %f\n", timerCPU.GetCounter());

    double rmserror = 0., norm = 0.;
    for (int n = 0; n < N; n++) {
        rmserror = rmserror + (h_xhatpruning[n][0] - h_xhat[n][0]) * (h_xhatpruning[n][0] - h_xhat[n][0]) + (h_xhatpruning[n][1] - h_xhat[n][1]) * (h_xhatpruning[n][1] - h_xhat[n][1]);
        norm = norm + h_xhat[n][0] * h_xhat[n][0] + h_xhat[n][1] * h_xhat[n][1];
    }
    printf("rmserror %f\n", 100. * sqrt(rmserror / norm));

    fftw_destroy_plan(plan_zp);

}

我开发的方法包括三个步骤：

将输入序列乘以"twiddle"复指数；
执行 fftw_many;
重新组织结果。

fftw_many 在 K * N 个输入点上比单个 FFTW 更快。然而，步骤#1 和#3 完全破坏了这种增益。我希望步骤 #1 和 #3 在计算上比步骤 #2 轻得多。

我的问题是：

第 1 步和第 3 步怎么可能比第 2 步的计算要求更高？
我如何改进步骤 #1 和 #3 以获得相对于 "standard" 方法的净收益？

非常感谢您的任何提示。

编辑

我正在使用 Visual Studio 2013 并在发布模式下编译。

Answer 1

运行更快的几个选项：

运行多线程，如果你只是运行单线程并且有多个内核可用。
创建并保存 FFTW wisdom 文件，尤其是在 FFT 维度事先已知的情况下。使用FFTW_EXHAUSTIVE，重新加载FFTW智慧，而不是每次都重新计算。如果您希望结果一致，这一点也很重要。由于 FFTW 可能会根据不同的计算智慧计算出不同的 FFT，并且智慧结果不一定总是相同，因此当给定相同的输入数据时，您的过程的不同运行可能会产生不同的结果。
如果您使用的是 x86，运行 64 位。 FFTW 算法非常占用寄存器，64 位模式下的 x86 CPU 运行ning 比 32 位模式下的运行ning 有更多可用的通用寄存器。位模式.
由于 FFTW 算法占用大量寄存器，我通过使用防止使用预取的编译器选项编译 FFTW，在提高 FFTW 性能方面取得了很好的成功并防止函数的隐式内联。

Answer 2

对于第三步，您可能想尝试切换循环的顺序：

for (int p = 0; p < N; p++) {
    for (int k = 0; k < K; k++) {
        h_xhatpruning[p * K + k][0] = h_xhatpruning_temp[p + k * N][0];
        h_xhatpruning[p * K + k][1] = h_xhatpruning_temp[p + k * N][1];
    }
}

因为存储地址连续比加载地址连续更有利。

不管怎样，你的访问模式都是缓存不友好的。您可以尝试使用块来改进这一点，例如假设 N 是 4 的倍数：

for (int p = 0; p < N; p += 4) {
    for (int k = 0; k < K; k++) {
        for (int p0 = 0; p0 < 4; p0++) {
            h_xhatpruning[(p + p0) * K + k][0] = h_xhatpruning_temp[(p + p0) + k * N][0];
            h_xhatpruning[(p + p0) * K + k][1] = h_xhatpruning_temp[(p + p0) + k * N][1];
        }
    }
}

这应该有助于在一定程度上减少高速缓存行的流失。如果确实如此，那么也许还可以尝试使用 4 以外的块大小来查看是否存在 "sweet spot".

Answer 3

我还根据 Paul R 的评论改进了我的代码。现在，替代方法比标准（零填充）方法更快。下面是完整的 C++ 脚本。对于第 1 步和第 3 步，我已经评论了其他尝试过的解决方案，这些解决方案显示出与未评论的解决方案一样慢或一样快。我有特权非嵌套 for 循环，也是为了更简单的未来 CUDA 并行化。我还没有为 FFTW 使用多线程。

#include <stdio.h>
#include <stdlib.h>     /* srand, rand */
#include <time.h>       /* time */
#include <math.h>
#include <fstream>

#include <omp.h>

#include <fftw3.h>

#include "TimingCPU.h"

#define PI_d            3.141592653589793

/******************/
/* STEP #1 ON CPU */
/******************/
void step1CPU(fftw_complex * __restrict h_xpruning, const fftw_complex * __restrict h_x, const int N, const int K) {

//  double factor = -2. * PI_d / (K * N);
//  int n;
//  omp_set_nested(1);
//#pragma omp parallel for private(n) num_threads(4)
//  for (int k = 0; k < K; k++) {
//      double arg1 = factor * k;
//#pragma omp parallel for num_threads(4)
//      for (n = 0; n < N; n++) {
//          double arg = arg1 * n;
//          double cosarg = cos(arg);
//          double sinarg = sin(arg);
//          h_xpruning[k * N + n][0] = h_x[n][0] * cosarg - h_x[n][1] * sinarg;
//          h_xpruning[k * N + n][1] = h_x[n][0] * sinarg + h_x[n][1] * cosarg;
//      }
//  }

    //double factor = -2. * PI_d / (K * N);
    //int k;
    //omp_set_nested(1);
    //#pragma omp parallel for private(k) num_threads(4)
    //for (int n = 0; n < N; n++) {
    //  double arg1 = factor * n;
    //  #pragma omp parallel for num_threads(4)
    //  for (k = 0; k < K; k++) {
    //      double arg = arg1 * k;
    //      double cosarg = cos(arg);
    //      double sinarg = sin(arg);
    //      h_xpruning[k * N + n][0] = h_x[n][0] * cosarg - h_x[n][1] * sinarg;
    //      h_xpruning[k * N + n][1] = h_x[n][0] * sinarg + h_x[n][1] * cosarg;
    //  }
    //}

    //double factor = -2. * PI_d / (K * N);
    //for (int k = 0; k < K; k++) {
    //  double arg1 = factor * k;
    //  for (int n = 0; n < N; n++) {
    //      double arg = arg1 * n;
    //      double cosarg = cos(arg);
    //      double sinarg = sin(arg);
    //      h_xpruning[k * N + n][0] = h_x[n][0] * cosarg - h_x[n][1] * sinarg;
    //      h_xpruning[k * N + n][1] = h_x[n][0] * sinarg + h_x[n][1] * cosarg;
    //  }
    //}

    //double factor = -2. * PI_d / (K * N);
    //for (int n = 0; n < N; n++) {
    //  double arg1 = factor * n;
    //  for (int k = 0; k < K; k++) {
    //      double arg = arg1 * k;
    //      double cosarg = cos(arg);
    //      double sinarg = sin(arg);
    //      h_xpruning[k * N + n][0] = h_x[n][0] * cosarg - h_x[n][1] * sinarg;
    //      h_xpruning[k * N + n][1] = h_x[n][0] * sinarg + h_x[n][1] * cosarg;
    //  }
    //}

    double factor = -2. * PI_d / (K * N);
    #pragma omp parallel for num_threads(8)
    for (int n = 0; n < K * N; n++) {
        int row = n / N;
        int col = n % N;
        double arg = factor * row * col;
        double cosarg = cos(arg);
        double sinarg = sin(arg);
        h_xpruning[n][0] = h_x[col][0] * cosarg - h_x[col][1] * sinarg;
        h_xpruning[n][1] = h_x[col][0] * sinarg + h_x[col][1] * cosarg;
    }
}

/******************/
/* STEP #3 ON CPU */
/******************/
void step3CPU(fftw_complex * __restrict h_xhatpruning, const fftw_complex * __restrict h_xhatpruning_temp, const int N, const int K) {

    //int k;
    //omp_set_nested(1);
    //#pragma omp parallel for private(k) num_threads(4)
    //for (int p = 0; p < N; p++) {
    //  #pragma omp parallel for num_threads(4)
    //  for (k = 0; k < K; k++) {
    //      h_xhatpruning[p * K + k][0] = h_xhatpruning_temp[p + k * N][0];
    //      h_xhatpruning[p * K + k][1] = h_xhatpruning_temp[p + k * N][1];
    //  }
    //} 

    //int p;
    //omp_set_nested(1);
    //#pragma omp parallel for private(p) num_threads(4)
    //for (int k = 0; k < K; k++) {
    //  #pragma omp parallel for num_threads(4)
    //  for (p = 0; p < N; p++) {
    //      h_xhatpruning[p * K + k][0] = h_xhatpruning_temp[p + k * N][0];
    //      h_xhatpruning[p * K + k][1] = h_xhatpruning_temp[p + k * N][1];
    //  }
    //}

    //for (int p = 0; p < N; p++) {
    //  for (int k = 0; k < K; k++) {
    //      h_xhatpruning[p * K + k][0] = h_xhatpruning_temp[p + k * N][0];
    //      h_xhatpruning[p * K + k][1] = h_xhatpruning_temp[p + k * N][1];
    //  }
    //}

    //for (int k = 0; k < K; k++) {
    //  for (int p = 0; p < N; p++) {
    //      h_xhatpruning[p * K + k][0] = h_xhatpruning_temp[p + k * N][0];
    //      h_xhatpruning[p * K + k][1] = h_xhatpruning_temp[p + k * N][1];
    //  }
    //}

    #pragma omp parallel for num_threads(8)
    for (int p = 0; p < K * N; p++) {
        int col = p % N;
        int row = p / K;
        h_xhatpruning[col * K + row][0] = h_xhatpruning_temp[col + row * N][0];
        h_xhatpruning[col * K + row][1] = h_xhatpruning_temp[col + row * N][1];
    }

    //for (int p = 0; p < N; p += 2) {
    //  for (int k = 0; k < K; k++) {
    //      for (int p0 = 0; p0 < 2; p0++) {
    //          h_xhatpruning[(p + p0) * K + k][0] = h_xhatpruning_temp[(p + p0) + k * N][0];
    //          h_xhatpruning[(p + p0) * K + k][1] = h_xhatpruning_temp[(p + p0) + k * N][1];
    //      }
    //  }
    //}

}

/********/
/* MAIN */
/********/
void main() {

    int N = 10;
    int K = 100000;

    // --- CPU memory allocations
    fftw_complex *h_x = (fftw_complex *)malloc(N     * sizeof(fftw_complex));
    fftw_complex *h_xzp = (fftw_complex *)calloc(N * K, sizeof(fftw_complex));
    fftw_complex *h_xpruning = (fftw_complex *)malloc(N * K * sizeof(fftw_complex));
    fftw_complex *h_xhatpruning = (fftw_complex *)malloc(N * K * sizeof(fftw_complex));
    fftw_complex *h_xhatpruning_temp = (fftw_complex *)malloc(N * K * sizeof(fftw_complex));
    fftw_complex *h_xhat = (fftw_complex *)malloc(N * K * sizeof(fftw_complex));
    //double2        *h_xhatGPU = (double2 *)malloc(N * K * sizeof(double2));


    // --- Random number generation of the data sequence on the CPU - moving the data from CPU to GPU
    srand(time(NULL));
    for (int k = 0; k < N; k++) {
        h_x[k][0] = (double)rand() / (double)RAND_MAX;
        h_x[k][1] = (double)rand() / (double)RAND_MAX;
    }
    //gpuErrchk(cudaMemcpy(d_x, h_x, N * sizeof(double2), cudaMemcpyHostToDevice));

    memcpy(h_xzp, h_x, N * sizeof(fftw_complex));

    // --- FFTW and cuFFT plans
    fftw_plan h_plan_zp      = fftw_plan_dft_1d(N * K, h_xzp, h_xhat, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_plan h_plan_pruning = fftw_plan_many_dft(1, &N, K, h_xpruning, NULL, 1, N, h_xhatpruning_temp, NULL, 1, N, FFTW_FORWARD, FFTW_ESTIMATE);

    double totalTimeCPU = 0., totalTimeGPU = 0.;
    double partialTimeCPU, partialTimeGPU;

    /****************************/
    /* STANDARD APPROACH ON CPU */
    /****************************/
    printf("Number of processors available = %i\n", omp_get_num_procs());
    printf("Number of threads              = %i\n", omp_get_max_threads());

    TimingCPU timerCPU;
    timerCPU.StartCounter();
    fftw_execute(h_plan_zp);
    printf("\nStadard on CPU: \t \t %f\n", timerCPU.GetCounter());

    /******************/
    /* STEP #1 ON CPU */
    /******************/
    timerCPU.StartCounter();
    step1CPU(h_xpruning, h_x, N, K);
    partialTimeCPU = timerCPU.GetCounter();
    totalTimeCPU = totalTimeCPU + partialTimeCPU;
    printf("\nOptimized first step CPU: \t %f\n", totalTimeCPU);

    /******************/
    /* STEP #2 ON CPU */
    /******************/
    timerCPU.StartCounter();
    fftw_execute(h_plan_pruning);
    partialTimeCPU = timerCPU.GetCounter();
    totalTimeCPU = totalTimeCPU + partialTimeCPU;
    printf("Optimized second step CPU: \t %f\n", timerCPU.GetCounter());

    /******************/
    /* STEP #3 ON CPU */
    /******************/
    timerCPU.StartCounter();
    step3CPU(h_xhatpruning, h_xhatpruning_temp, N, K);
    partialTimeCPU = timerCPU.GetCounter();
    totalTimeCPU = totalTimeCPU + partialTimeCPU;
    printf("Optimized third step CPU: \t %f\n", partialTimeCPU);

    printf("Total time CPU: \t \t %f\n", totalTimeCPU);

    double rmserror = 0., norm = 0.;
    for (int n = 0; n < N; n++) {
        rmserror = rmserror + (h_xhatpruning[n][0] - h_xhat[n][0]) * (h_xhatpruning[n][0] - h_xhat[n][0]) + (h_xhatpruning[n][1] - h_xhat[n][1]) * (h_xhatpruning[n][1] - h_xhat[n][1]);
        norm = norm + h_xhat[n][0] * h_xhat[n][0] + h_xhat[n][1] * h_xhat[n][1];
    }
    printf("\nrmserror %f\n", 100. * sqrt(rmserror / norm));

    fftw_destroy_plan(h_plan_zp);

}

案例

N = 10
K = 100000

我的时间如下

Stadard on CPU:                  23.895417

Optimized first step CPU:        4.472087
Optimized second step CPU:       4.926603
Optimized third step CPU:        2.394958
Total time CPU:                  11.793648

加速 FFTW 修剪以避免大量零填充

Accelerating FFTW pruning to avoid massive zero padding

c++

c

fft

fftw

pruning