x(n)がK * Nであり、最初のN要素のみがゼロと異なると仮定します。私はN << Kと言って、例えば、N = 10とK = 100000と言っています。私はそのようなシーケンスのFFTWでFFTを計算したいと思います。これは、長さがNで、ゼロパディングがK * Nになっていることに相当します。 NとKは「大」である可能性があるため、大きなゼロ詰めがあります。明示的なゼロパディングを避けるために計算時間を節約できるかどうか検討しています。大量のゼロパディングを避けるためにFFTWプルーニングを高速化する
場合
K = 2は、私たちはケースK = 2を考慮することから始めましょう。 kがDFTのような値を算出することができることを意味しても、すなわちk = 2 * m、次いで
ある場合この場合、x(n)のDFTは
のように書くことができます。長さがNであり、K * NではないFFTである。
k場合は、DFTのような値が再び長NのシーケンスのFFTにより算出されずK * Nすることができることを意味する、すなわちk = 2 * m + 1、次いで
奇数です。
結論として、長さ2 * Nの単一のFFTを2の長さのNのFFTと交換することができます。
この場合、任意のK
の場合、我々はk = m * K + tを書くことに
を持って、我々は結論で、だから、
を持っています、私はexchaすることができます長さがK * Nの単一のFFTをKの長さのFFTに設定します。長さはNです。 FFTWにはfftw_plan_many_dftがあるので、1つのFFTの場合に対していくらかの利益を得ることができます。それを確認する
、私が開発したアプローチは、3つのステップから構成され、次のコード
#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を実行する。- 結果を再編成する。
入力ポイントがK * Nの場合、fftw_manyは単一のFFTWよりも高速です。しかしながら、ステップ1及びステップ3は、このような利得を完全に破壊する。私はステップ＃1と＃3がステップ＃2より計算上非常に軽いと期待します。
私の質問は以下のとおりです。
- どのようにそれは＃1と＃3ステップということも可能であるので、計算ステップ＃2よりも厳しいですか？
- ステップ＃1と＃3を改善して、「標準」のアプローチに対して純利益を得るにはどうすればよいですか？
ありがとうございました。
EDIT
私は、Visual Studio 2013での作業とリリースモードでコンパイルしています。より高速に実行する
あなたは、例えば、有効な最適化を使用してテスト・コードをコンパイルしたのです'gcc -O3 ...'？ –
@PaulRリリースモードでVisual Studio 2013を使用してコンパイルしています。私はそれに応じて投稿を編集しました。 – JackOLantern