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

Question

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

案例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与 2长度为 N 的 FFT .

任意大小写K

在这种情况下，我们有

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

因此，总而言之，我可以交换长度为 K * N 的单个 FFT与 K长度为 N 的 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);

}

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

1. 将输入序列乘以“twiddle”复指数；
2. 执行 fftw_many ;
3. 重组结果。

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

我的问题是:

1. 第 1 步和第 3 步的计算要求怎么可能比第 2 步高？
2. 我如何改进第 1 步和第 3 步以获得相对于“标准”方法的净 yield ？

非常感谢您的任何提示。

编辑

我正在使用 Visual Studio 2013 并在 Release模式下进行编译。

Answer 1

几个运行速度更快的选项:

1. 如果您只运行单线程并且有多核可用，则运行多线程。

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

3. 如果您使用的是 x86，请运行 64 位。 FFTW 算法占用大量寄存器，运行在 64 位模式下的 x86 CPU 比运行在 32 位模式下的 x86 CPU 有更多可用的通用寄存器。

4. 由于 FFTW 算法的寄存器密集程度如此之高，我通过使用编译器选项编译 FFTW 来成功提高 FFTW 性能，这些选项防止使用预取并防止函数的隐式内联.