Дискретное преобразование Фурье в fftw3 и в CUDA

Если вам вдруг понадобилось вычислить дискретное преобразование Фурье, то для этой цели прекрасно подойдут два инструмента: http://www.fftw.org/ и CUDA. К достоинствам первого можно отнести его открытость, специализация именно на преобразовании Фурье, предельная простота установки. 

Рассмотрим для начала пример с использованием fftw.

Создадим файл fftw3.c следующего содержания:

#include <stdlib.h>
#include <stdio.h>
#include <fftw3.h>

int
main ( void )
{
    int i;
    fftw_complex *in;
    fftw_complex *in2;
    int n = 10;
    fftw_complex *out;
    fftw_plan plan_backward;
    fftw_plan plan_forward;

    //Создаем входные данные
    in = fftw_malloc ( sizeof ( fftw_complex ) * n );
    printf ( "Input Data:\n" );
    for ( i = 0; i < n; i++ ){
        in[i][0] = i;in[i][1] = 0;
        printf ( "  %3d  %12f  %12f\n", i, in[i][0], in[i][1] );
    }

    //Выполняем прямое преобразование Фурье
    out = fftw_malloc ( sizeof ( fftw_complex ) * n );
    plan_forward = fftw_plan_dft_1d ( n, in, out, FFTW_FORWARD, FFTW_ESTIMATE );
    fftw_execute ( plan_forward );

    printf ( "Output FFT Coefficients:\n" );
    for ( i = 0; i < n; i++ )
        printf ( "  %3d  %12f  %12f\n", i, out[i][0], out[i][1] );

    //Выполняем обратное преобразование Фурье
    in2 = fftw_malloc ( sizeof ( fftw_complex ) * n );
    plan_backward = fftw_plan_dft_1d ( n, out, in2, FFTW_BACKWARD, FFTW_ESTIMATE );
    fftw_execute ( plan_backward );

    printf ( "Recovered input data:\n" );
    for ( i = 0; i < n; i++ )
        printf ( "  %3d  %12f  %12f\n", i, in2[i][0]/(double)(n), in2[i][1]/(double)(n) );

    //Прибираем за собой
    fftw_destroy_plan ( plan_forward );
    fftw_destroy_plan ( plan_backward );
    fftw_free ( in );
    fftw_free ( in2 );
    fftw_free ( out );

  return 0;
}

Компилируем его с помощью следующей команды:

#gcc fftw3.c -o fftw3 -lfftw3 -lm

 Как видите, если убрать код, выводящий данные на стандартный вывод, то остальной код, выполняющий главную функциональность, будет смотреться предельно просто и ясно. Всего несколько строчек. Для завершения картины запускаем программу:

#./fftw3

и изучаем вывод. Должны увидеть, что обратное преобразование восстанавливает входные данные для прямого преобразования.

Далее рассмотрим решение этой же задачи средствами CUDA. К недостаткам этого метода можно отнести не всегда простую установку CUDA  и обязательное наличие видеокарты, поддерживающую эту технологию. К достоинствам отнесем всю мощь данной технологии. Код взят из примеров, поставляемых вместе с тулкитом с сайта NVIDIA. Называется он (код, с помощью которого была написана программа, исходники которой далее приводятся) иронично simpleCUFFT. Хотя простоты там, на мой взгляд нет совсем, потому что пример, показывающий как пользоваться преобразованием Фурье, завален кодом, делающим конволюцию сигнала с фильтром (зачем это людям, изучающим CUDA?), и выковырять пару строчек, которые и делают прямую задачу примера, было для меня не просто.  За деревьями леса не видно. Поэтому, думаю, этот пост будет полезен.

Итак, создаем файл truelySimpleFFT.cu следующего содержания:

#include<stdlib.h>
#include<stdio.h>
#include<string.h>
#include<math.h>

#include <cuda_runtime.h>
#include <cufft.h>

typedef float2 Complex;
#define N_ELEM 10
int
main(int argc, char **argv)
{
    Complex *h_signal;
    Complex *d_signal;
    int mem_size;
    cufftHandle plan;
   
    printf("[truely simpleCUFFT] is starting...\n");
    //Выделяем память для входных данных на хосте
    h_signal = (Complex *)malloc(sizeof(Complex) * N_ELEM);

    //Инициализируем входные данные
    for (unsigned int i = 0; i < N_ELEM; ++i){
        h_signal[i].x = (float)i;
        h_signal[i].y = 0;
        printf("%f\t%f\n", h_signal[i].x, h_signal[i].y);
    }  
   
    //Выделяем память для входных данных на видеокарте и копируем их туда
    mem_size = sizeof(Complex) * N_ELEM;     
    cudaMalloc((void **)&d_signal, mem_size);                                  
    cudaMemcpy(d_signal, h_signal, mem_size, cudaMemcpyHostToDevice);                 
     
    //Собственно выполняем преобразование Фурье
    printf("Transforming signal cufftExecC2C\n");
    cufftPlan1d(&plan, N_ELEM, CUFFT_C2C, 1);                        
    cufftExecC2C(plan, (cufftComplex *)d_signal, (cufftComplex *)d_signal, CUFFT_FORWARD);
   
    //Копируем результат вычисления с видеокарты на хост
    cudaMemcpy(h_signal, d_signal, mem_size, cudaMemcpyDeviceToHost);
    for (unsigned int i = 0; i < N_ELEM; ++i)
        printf("%f\t%f\n", h_signal[i].x, h_signal[i].y);
   
    //Прибираем за собой
    cufftDestroy(plan);
    free(h_signal);
    cudaFree(d_signal);
    cudaDeviceReset();
   
    exit(0);
}

 

Компилируем его с помощью следующих двух команд:

#nvcc  -o simpleCUFFT.o -c truelySimpleFFT.cu

#g++ -o simpleCUFFT simpleCUFFT.o -lcudart -lcudart -lcufft
Необходимо отметить следующее. Если вы устанавливали SDK с сайта NVIDIA "руками", то вам надо вставить в выше приведенные команды ключи -L и -I с соответствующими путями к заголовочным файлам и библиотекам, которые вы указывали при установке (а также nvcc с полным путем). Я устанавливал с помощью репозитария (пользуюсь Дебиан) и в моем случае надо вводить команды как написано выше.

Также обратите внимание на подобие кода, выполняющего преобразование Фурье в обоих примерах.

Запускаем на выполнение командой:

 #./simpleCUFFT

Обе программы компиллировались и ваполнялись на Debian GNU/Linux 7.0 (wheezy).