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

Question

假设我有一个x(n)很K * N长的序列,只有第一个N元素与零不同.我假设,N << K例如,N = 10和K = 100000.我想通过FFTW计算这种序列的FFT.这相当于具有一系列长度N并且具有零填充K * N.由于N并且K可能是"大",我有一个显着的零填充.我正在探索是否可以节省一些计算时间,避免显式零填充.

案子 K = 2

让我们从考虑案例开始K = 2.在这种情况下,DFT x(n)可以写成

如果k是偶数k = 2 * m,那么

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

如果k是奇数k = 2 * m + 1,那么

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

因此,总之,我可以用长度2 * N的2FFT 交换长度的单个FFT N.

任意的情况 K

在这种情况下,我们有

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

因此,总之,我可以用长度K * N的KFFT 交换长度的单个FFT N.由于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);

}

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

1. 将输入序列乘以"旋转"复数指数;
2. 表演fftw_many;
3. 重组结果.

它fftw_many比K * N输入点上的单个FFTW更快.但是,步骤#1和#3完全破坏了这种增益.我希望步骤#1和#3在计算上比步骤#2轻得多.

我的问题是:

1. 步骤#1和#3如何比步骤#2的计算要求更高？
2. 如何改进步骤#1和#3以获得"标准"方法的净收益？

非常感谢您的任何提示.

编辑

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

Answer 1

几个选项可以更快地运行:

1. 如果您只运行单线程并且有多个可用核心,则运行多线程.

2. 创建并保存FFTW智能文件,尤其是在预先知道FFT尺寸的情况下.使用FFTW_EXHAUSTIVE并重新加载FFTW智慧,而不是每次重新计算它.如果您希望结果一致,这也很重要.由于FFTW可能以不同的计算智慧计算FFT,并且智慧结果不一定总是相同的,因此当给出相同的输入数据时,不同的运行过程可能会产生不同的结果.

3. 如果您使用的是x86,请运行64位.FFTW算法非常容易注册,并且运行在64位模式下的x86 CPU比32位模式下运行的通用寄存器要多得多.

4. 由于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 以外的块大小，看看是否存在“最佳点”。