7
假设我有一个序列x(n)，它的长度为K * N，并且只有第一个N元素与零不同。我假设N << K，例如，N = 10和K = 100000。我想通过FFTW计算这样一个序列的FFT。这相当于具有长度为N的序列并具有零填充为K * N。由于N和K可能是“大”，我有一个重要的零填充。我在探索是否可以节省一些计算时间，避免显式零填充。加速FFTW修剪以避免大量零填充
的情况下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。
因此，总而言之，我可以用长度为N的长度为2 * N的长度为2的FFT交换单个FFT。
的任意K
的情况下，在这种情况下，我们有
上书写k = m * K + t，我们有
所以，综上所述，我可以交换nge长度为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);
}
我已经开发的方法包括三个步骤：
- 通过“玩弄”复指数相乘的输入序列;
- 执行
fftw_many; - 重组结果。
fftw_many比在K * N输入点上的单个FFTW快。然而，步骤＃1和＃3完全摧毁了这样的收益。我希望步骤＃1和＃3在计算上比步骤＃2轻得多。
我的问题是：
- 这怎么可能，步骤＃1和＃3这样计算量比第2步的要求更高？
- 如何改善步骤＃1和＃3以获得与“标准”方法相比的净收益？
非常感谢您的任何提示。
编辑
我与Visual Studio 2013的工作，并在Release模式下进行编译。
你编译优化测试代码启用，例如'gcc -O3 ...'？ –
@PaulR我在发布模式下使用Visual Studio 2013进行编译。我已经相应地编辑了这篇文章。 – JackOLantern