Предположим, что у меня есть последовательность x(n) длиной K * N и что только первые N элемента отличны от нуля. Я предполагаю, что N << K, скажем, например, N = 10 и K = 100000. Я хочу рассчитать FFT с помощью FFTW такой последовательности. Это эквивалентно последовательности длиной N и дополнению нулями до K * N. Поскольку N и K могут быть "большими", у меня есть значительное заполнение нулями. Я изучаю, могу ли я сэкономить время вычислений, избегая явного заполнения нулями.
Дело K = 2
Начнем с рассмотрения случая K = 2. В этом случае ДПФ x(n) можно записать как
Если k четное, а именно k = 2 * m, то
что означает, что такие значения ДПФ могут быть вычислены через БПФ последовательности длиной N, а не K * N.
Если k нечетное, а именно k = 2 * m + 1, то
что означает, что такие значения ДПФ можно снова вычислить через БПФ последовательности длины N, а не K * N.
Итак, в заключение я могу обменять одно БПФ длиной 2 * N на 2 БПФ длиной N.
Случай произвольного K
В этом случае мы имеем
При записи k = m * K + t имеем
Итак, в заключение я могу обменять одно БПФ длиной K * N на K БПФ длиной N. Поскольку FFTW имеет fftw_plan_many_dft, я могу ожидать некоторого выигрыша по сравнению со случаем одного БПФ.
Чтобы убедиться в этом, я установил следующий код
#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 быстрее, чем одиночный FFTW на K * N входных точках. Однако шаги №1 и №3 полностью разрушают такой выигрыш. Я ожидаю, что шаги № 1 и № 3 будут намного легче в вычислительном отношении, чем шаг № 2.
Мои вопросы:
- Как возможно, что шаги № 1 и № 3 более требовательны к вычислительным ресурсам, чем шаг № 2?
- Как я могу улучшить шаги № 1 и № 3, чтобы получить чистую прибыль по сравнению со «стандартным» подходом?
Большое спасибо за любую подсказку.
ИЗМЕНИТЬ
Я работаю с Visual Studio 2013 и компилирую в режиме Release.
gcc -O3 ...? - person Paul R schedule 16.11.2016