長いシーケンスx(n)があり、最初の要素だけがゼロと異なるとします。たとえば、とと仮定しています。そのようなシーケンスのFFTWをFFTWで計算したいと思います。これは、 length のシーケンスを持ち、 にゼロをパディングすることと同じです。とは「大きい」可能性があるため、かなりのゼロ パディングがあります。明示的なゼロパディングを避けて、計算時間を節約できるかどうかを調査しています。K * NNN << KN = 10K = 100000NK * NNK
ケースK = 2
ケースを検討することから始めましょうK = 2。この場合、 の DFT は次のx(n)ように記述できます。
kが偶数、つまりの場合k = 2 * m、
つまり、このような DFT の値は、Nではなく長さ のシーケンスの FFT を介して計算できますK * N。
kが奇数、つまりの場合k = 2 * m + 1、
つまり、このような DFT の値は、長さ のシーケンスの FFT を介して再度計算できますが、 ではNありませんK * N。
したがって、結論として、長さの 1 つの FFT を長さの FFT と交換2 * Nでき2ますN。
恣意的なケースK
この場合、
書くk = m * K + tと、私たちは持っています
したがって、結論として、長さの 1 つの FFT を長さの FFT と交換K * NできKます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);
}
私が開発したアプローチは、次の 3 つのステップで構成されています。
- 入力シーケンスを「ひねる」複雑な指数で乗算します。
- 実行する
fftw_many; - 結果の再編成。
は、入力ポイントfftw_manyで単一の FFTW よりも高速です。K * Nただし、ステップ 1 と 3 は、そのようなゲインを完全に破壊します。ステップ 1 と 3 は、ステップ 2 よりもはるかに軽量であると予想されます。
私の質問は次のとおりです。
- ステップ 1 と 3 が、ステップ 2 よりも計算量が多くなる可能性があるのはなぜですか?
- ステップ 1 と 3 を改善して、「標準的な」アプローチに対して純利益を得るにはどうすればよいですか?
ヒントをありがとうございました。
編集
私は Visual Studio 2013 で作業しており、リリース モードでコンパイルしています。