2013-04-06 3 views
2

Я пытаюсь реализовать фильтр FIR (конечный импульсный отклик) в CUDA. Мой подход достаточно прост и выглядит примерно так:FIR-фильтр в CUDA (как 1D свертка)

#include <cuda.h> 

__global__ void filterData(const float *d_data, 
          const float *d_numerator, 
          float *d_filteredData, 
          const int numeratorLength, 
          const int filteredDataLength) 
{ 
    int i = blockDim.x * blockIdx.x + threadIdx.x; 

    float sum = 0.0f; 

    if (i < filteredDataLength) 
    { 
     for (int j = 0; j < numeratorLength; j++) 
     { 
      // The first (numeratorLength-1) elements contain the filter state 
      sum += d_numerator[j] * d_data[i + numeratorLength - j - 1]; 
     } 
    } 

    d_filteredData[i] = sum; 
} 

int main(void) 
{ 
    // (Skipping error checks to make code more readable) 

    int dataLength = 18042; 
    int filteredDataLength = 16384; 
    int numeratorLength= 1659; 

    // Pointers to data, filtered data and filter coefficients 
    // (Skipping how these are read into the arrays) 
    float *h_data = new float[dataLength]; 
    float *h_filteredData = new float[filteredDataLength]; 
    float *h_filter = new float[numeratorLength]; 


    // Create device pointers 
    float *d_data = nullptr; 
    cudaMalloc((void **)&d_data, dataLength * sizeof(float)); 

    float *d_numerator = nullptr; 
    cudaMalloc((void **)&d_numerator, numeratorLength * sizeof(float)); 

    float *d_filteredData = nullptr; 
    cudaMalloc((void **)&d_filteredData, filteredDataLength * sizeof(float)); 


    // Copy data to device 
    cudaMemcpy(d_data, h_data, dataLength * sizeof(float), cudaMemcpyHostToDevice); 
    cudaMemcpy(d_numerator, h_numerator, numeratorLength * sizeof(float), cudaMemcpyHostToDevice); 

    // Launch the kernel 
    int threadsPerBlock = 256; 
    int blocksPerGrid = (filteredDataLength + threadsPerBlock - 1)/threadsPerBlock; 
    filterData<<<blocksPerGrid,threadsPerBlock>>>(d_data, d_numerator, d_filteredData, numeratorLength, filteredDataLength); 

    // Copy results to host 
    cudaMemcpy(h_filteredData, d_filteredData, filteredDataLength * sizeof(float), cudaMemcpyDeviceToHost); 

    // Clean up 
    cudaFree(d_data); 
    cudaFree(d_numerator); 
    cudaFree(d_filteredData); 

    // Do stuff with h_filteredData... 

    // Clean up some more 
    delete [] h_data; 
    delete [] h_filteredData; 
    delete [] h_filter; 
} 

Фильтр работает, но, как я новичок в программировании CUDA, и я не знаю, как оптимизировать его.

Небольшая проблема, которую я вижу в том, что dataLength, filteredDataLength и numeratorLength не известны заранее в заявке я намерен использовать фильтр. Кроме того, даже если dataLength кратно 32 в приведенном выше коде, это не гарантируется, что в конечном приложении.

Когда я сравниваю свой код выше с ArrayFire, мой код занимает примерно в три раза больше времени для выполнения.

Есть ли у кого-нибудь идеи о том, как ускорить процесс?

EDIT: Сменили все filterLength на numeratorLength.

+1

ли 'numeratorLength' так же, как' filterLength'? Я не вижу определения 'numatorLength' где-либо в том, что вы разместили. Эта проблема, по сути, является проблемой одномерного трафарета. Стандартная оптимизация для проблем с трафаретом заключается в том, чтобы довести часть входных данных в разделяемую память, достаточную для потоковых потоков блоков для вычисления их выходов, а затем позволить этим потокам работать из копии общей памяти. –

+1

Если вы в конечном итоге избиваете ArrayFire, дайте нам знать! Если нет, вы всегда можете просто использовать ArrayFire, так как это быстрее :) – arrayfire

+0

@RobertCrovella Да, numatorLength такая же, как filterLength. Я решил изменить имя, но, видимо, пропустил несколько мест. Мой плохой, извините. Я изменил исходный пост так, что есть только numatorLength. Спасибо за предложение использовать общую память. Я читал, что они значительно быстрее, чем глобальная память, но я немного уверен в том, как наилучшим образом реализовать это, поскольку общая память ограничена по размеру, а длины фильтров могут быть довольно длинными. Я буду играть с ним и посмотреть, как это происходит. – Elfendahl

ответ

1

я могу предложить следующее, чтобы ускорить ваш код:

  1. Использование разделяемой памяти: это маленький кэш-памяти, как, но чрезвычайно быстрее, чем мировой карты памяти. Вы можете найти больше об этом , ища ключевое слово __shared__ в документации CUDA. Например, для вы можете предварительно получить числители фильтра и большие фрагменты данных в общей памяти, что значительно повысит производительность вашего . Вы должны обратить особое внимание на данные выравнивание в этом случае, так как это действительно имеет значение, и это может замедлить ваш код .
  2. Подумайте о разворачивании петли числителя сумма. Вы можете проверить пример сокращения вектора в документации CUDA .
  3. Вы также можете подумать о распараллеливании цикла числителя сам по себе. Это можно сделать, добавив дополнительный размер (например, 'y') к вашему блоку потоков. Вам нужно будет сделать также общий вектор, который имеет размерность numberatorLength. Вы также можете проверить пример вектора уменьшения на том, как быстро взять сумму этого вектора в конце.
+0

Спасибо за предложения! Я рассмотрю использование разделяемой памяти, и пример сокращения вектора кажется очень интересным. Но мне особенно нравится ваше третье предложение, потому что цикл цикла в моем нынешнем коде выглядит как потенциальная шея бутылки при увеличении длины фильтра. – Elfendahl

1

Вы пытаетесь вычислить выход фильтра, непосредственно оценивая 1D свертку через ядро ​​CUDA.

В случае, когда длительность отклика фильтра составляет long, вы можете сделать что-то, чтобы оценить отфильтрованный вход, выполняет вычисления непосредственно в сопряженном домене с использованием БПФ. Ниже я сообщаю пример кода с использованием CUDA Thrust и библиотеки cuFFT. Это прямой перевод, например Matlab на основе сообщили в

Low-Pass Filtering by FFT Convolution

Позволь мне отрицает, что некоторые оптимизации возможны с этим кодом, но я предпочел оставить его, как это так, что она может быть легко по сравнению с его аналогом Matlab.

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

#include <cufft.h> 

#include <thrust\device_vector.h> 
#include <thrust\sequence.h> 

#define pi_f 3.14159265358979f     // Greek pi in single precision 

/****************/ 
/* SIN OPERATOR */ 
/****************/ 
class sin_op { 

    float fk_, Fs_; 

    public: 

     sin_op(float fk, float Fs) { fk_ = fk; Fs_ = Fs; } 

     __host__ __device__ float operator()(float x) const { return sin(2.f*pi_f*x*fk_/Fs_); } 
}; 

/*****************/ 
/* SINC OPERATOR */ 
/*****************/ 
class sinc_op { 

    float fc_, Fs_; 

    public: 

     sinc_op(float fc, float Fs) { fc_ = fc; Fs_ = Fs; } 

     __host__ __device__ float operator()(float x) const 
     { 
      if (x==0) return (2.f*fc_/Fs_); 
      else   return (2.f*fc_/Fs_)*sin(2.f*pi_f*fc_*x/Fs_)/(2.f*pi_f*fc_*x/Fs_); 
     } 
}; 

/********************/ 
/* HAMMING OPERATOR */ 
/********************/ 
class hamming_op { 

    int L_; 

    public: 

     hamming_op(int L) { L_ = L; } 

     __host__ __device__ float operator()(int x) const 
     { 
      return 0.54-0.46*cos(2.f*pi_f*x/(L_-1)); 
     } 
}; 


/*********************************/ 
/* MULTIPLY CUFFTCOMPLEX NUMBERS */ 
/*********************************/ 
struct multiply_cufftComplex { 
    __device__ cufftComplex operator()(const cufftComplex& a, const cufftComplex& b) const { 
     cufftComplex r; 
     r.x = a.x * b.x - a.y * b.y; 
     r.y = a.x * b.y + a.y * b.x; 
     return r; 
    } 
}; 

/********/ 
/* MAIN */ 
/********/ 
void main(){ 

    // Signal parameters: 
    int M = 256;       // signal length 
    const int N = 4; 
    float f[N] = { 440, 880, 1000, 2000 };    // frequencies 
    float Fs = 5000.;      // sampling rate 

    // Generate a signal by adding up sinusoids: 
    thrust::device_vector<float> d_x(M,0.f);   // pre-allocate 'accumulator' 
    thrust::device_vector<float> d_n(M);    // discrete-time grid 
    thrust::sequence(d_n.begin(), d_n.end(), 0, 1); 

    thrust::device_vector<float> d_temp(M); 
    for (int i=0; i<N; i++) { 
     float fk = f[i]; 
     thrust::transform(d_n.begin(), d_n.end(), d_temp.begin(), sin_op(fk,Fs)); 
     thrust::transform(d_temp.begin(), d_temp.end(), d_x.begin(), d_x.begin(), thrust::plus<float>()); 
    } 

    // Filter parameters: 
    int L = 257;      // filter length 
    float fc = 600.f;     // cutoff frequency 

    // Design the filter using the window method: 
    thrust::device_vector<float> d_hsupp(L);    
    thrust::sequence(d_hsupp.begin(), d_hsupp.end(), -(L-1)/2, 1); 
    thrust::device_vector<float> d_hideal(L);   
    thrust::transform(d_hsupp.begin(), d_hsupp.end(), d_hideal.begin(), sinc_op(fc,Fs)); 
    thrust::device_vector<float> d_l(L);     
    thrust::sequence(d_l.begin(), d_l.end(), 0, 1); 
    thrust::device_vector<float> d_h(L);     
    thrust::transform(d_l.begin(), d_l.end(), d_h.begin(), hamming_op(L)); 
    // h is our filter 
    thrust::transform(d_hideal.begin(), d_hideal.end(), d_h.begin(), d_h.begin(), thrust::multiplies<float>()); 

    // --- Choose the next power of 2 greater than L+M-1 
    int Nfft = pow(2,(ceil(log2((float)(L+M-1))))); // or 2^nextpow2(L+M-1) 

    // Zero pad the signal and impulse response: 
    thrust::device_vector<float> d_xzp(Nfft,0.f); 
    thrust::device_vector<float> d_hzp(Nfft,0.f); 
    thrust::copy(d_x.begin(), d_x.end(), d_xzp.begin()); 
    thrust::copy(d_h.begin(), d_h.end(), d_hzp.begin()); 

    // Transform the signal and the filter: 
    cufftHandle plan; 
    cufftPlan1d(&plan, Nfft, CUFFT_R2C, 1); 
    thrust::device_vector<cufftComplex> d_X(Nfft/2+1); 
    thrust::device_vector<cufftComplex> d_H(Nfft/2+1); 
    cufftExecR2C(plan, (cufftReal*)thrust::raw_pointer_cast(d_xzp.data()), (cufftComplex*)thrust::raw_pointer_cast(d_X.data())); 
    cufftExecR2C(plan, (cufftReal*)thrust::raw_pointer_cast(d_hzp.data()), (cufftComplex*)thrust::raw_pointer_cast(d_H.data())); 

    thrust::device_vector<cufftComplex> d_Y(Nfft/2+1); 
    thrust::transform(d_X.begin(), d_X.end(), d_H.begin(), d_Y.begin(), multiply_cufftComplex()); 

    cufftPlan1d(&plan, Nfft, CUFFT_C2R, 1); 
    thrust::device_vector<float> d_y(Nfft); 
    cufftExecC2R(plan, (cufftComplex*)thrust::raw_pointer_cast(d_Y.data()), (cufftReal*)thrust::raw_pointer_cast(d_y.data())); 

    getchar(); 

} 
1

Кроме того, мой другой ответ, который я ожидаю будет более удобно для ядер свертки с длинной длительности, ниже я хочу сообщить другую реализацию, которая более соответствует первоначальной попытки OP и я ожидаю, будет более удобный для ядер свертки с короткий продолжительность. Такая реализация основана на рукописном ядре, использующем кеширование в общей памяти. Более подробную информацию можно найти в книге D.B. Кирк и В.-М. В. Hwu

Programming Massively Parallel Processors, Second Edition: A Hands-on Approach

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

#include "TimingGPU.cuh" 
#include "Utilities.cuh" 

#define RG   10 
#define BLOCKSIZE 8 

/****************/ 
/* CPU FUNCTION */ 
/****************/ 
void h_convolution_1D(const float * __restrict__ h_Signal, const float * __restrict__ h_ConvKernel, float * __restrict__ h_Result_CPU, 
         const int N, const int K) { 

    for (int i = 0; i < N; i++) { 

     float temp = 0.f; 

     int N_start_point = i - (K/2); 

     for (int j = 0; j < K; j++) if (N_start_point + j >= 0 && N_start_point + j < N) { 
      temp += h_Signal[N_start_point+ j] * h_ConvKernel[j]; 
     } 

     h_Result_CPU[i] = temp; 
    } 
} 

/********************/ 
/* BASIC GPU KERNEL */ 
/********************/ 
__global__ void d_convolution_1D_basic(const float * __restrict__ d_Signal, const float * __restrict__ d_ConvKernel, float * __restrict__ d_Result_GPU, 
             const int N, const int K) { 

    int i = blockIdx.x * blockDim.x + threadIdx.x; 

    float temp = 0.f; 

    int N_start_point = i - (K/2); 

    for (int j = 0; j < K; j++) if (N_start_point + j >= 0 && N_start_point + j < N) { 
     temp += d_Signal[N_start_point+ j] * d_ConvKernel[j]; 
    } 

    d_Result_GPU[i] = temp; 
} 

/***************************/ 
/* GPU KERNEL WITH CACHING */ 
/***************************/ 
__global__ void d_convolution_1D_caching(const float * __restrict__ d_Signal, const float * __restrict__ d_ConvKernel, float * __restrict__ d_Result_GPU, 
             const int N, const int K) { 

    int i = blockIdx.x * blockDim.x + threadIdx.x; 

    __shared__ float d_Tile[BLOCKSIZE]; 

    d_Tile[threadIdx.x] = d_Signal[i]; 
    __syncthreads(); 

    float temp = 0.f; 

    int N_start_point = i - (K/2); 

    for (int j = 0; j < K; j++) if (N_start_point + j >= 0 && N_start_point + j < N) { 

      if ((N_start_point + j >= blockIdx.x * blockDim.x) && (N_start_point + j < (blockIdx.x + 1) * blockDim.x)) 

       // --- The signal element is in the tile loaded in the shared memory 
       temp += d_Tile[threadIdx.x + j - (K/2)] * d_ConvKernel[j]; 

      else 

       // --- The signal element is not in the tile loaded in the shared memory 
       temp += d_Signal[N_start_point + j] * d_ConvKernel[j]; 

    } 

    d_Result_GPU[i] = temp; 
} 

/********/ 
/* MAIN */ 
/********/ 
int main(){ 

    const int N = 15;   // --- Signal length 
    const int K = 5;   // --- Convolution kernel length 

    float *h_Signal   = (float *)malloc(N * sizeof(float)); 
    float *h_Result_CPU  = (float *)malloc(N * sizeof(float)); 
    float *h_Result_GPU  = (float *)malloc(N * sizeof(float)); 
    float *h_ConvKernel  = (float *)malloc(K * sizeof(float)); 

    float *d_Signal;  gpuErrchk(cudaMalloc(&d_Signal,  N * sizeof(float))); 
    float *d_Result_GPU; gpuErrchk(cudaMalloc(&d_Result_GPU, N * sizeof(float))); 
    float *d_ConvKernel; gpuErrchk(cudaMalloc(&d_ConvKernel, K * sizeof(float))); 

    for (int i=0; i < N; i++) { h_Signal[i] = (float)(rand() % RG); } 

    for (int i=0; i < K; i++) { h_ConvKernel[i] = (float)(rand() % RG); } 

    gpuErrchk(cudaMemcpy(d_Signal,  h_Signal,  N * sizeof(float), cudaMemcpyHostToDevice)); 
    gpuErrchk(cudaMemcpy(d_ConvKernel, h_ConvKernel, K * sizeof(float), cudaMemcpyHostToDevice)); 

    h_convolution_1D(h_Signal, h_ConvKernel, h_Result_CPU, N, K); 

    d_convolution_1D_basic<<<iDivUp(N, BLOCKSIZE), BLOCKSIZE>>>(d_Signal, d_ConvKernel, d_Result_GPU, N, K); 
    gpuErrchk(cudaPeekAtLastError()); 
    gpuErrchk(cudaDeviceSynchronize()); 

    gpuErrchk(cudaMemcpy(h_Result_GPU, d_Result_GPU, N * sizeof(float), cudaMemcpyDeviceToHost)); 

    for (int i = 0; i < N; i++) if (h_Result_CPU[i] != h_Result_GPU[i]) {printf("mismatch2 at %d, cpu: %d, gpu %d\n", i, h_Result_CPU[i], h_Result_GPU[i]); return 1;} 

    printf("Test basic passed\n"); 

    d_convolution_1D_caching<<<iDivUp(N, BLOCKSIZE), BLOCKSIZE>>>(d_Signal, d_ConvKernel, d_Result_GPU, N, K); 
    gpuErrchk(cudaPeekAtLastError()); 
    gpuErrchk(cudaDeviceSynchronize()); 

    gpuErrchk(cudaMemcpy(h_Result_GPU, d_Result_GPU, N * sizeof(float), cudaMemcpyDeviceToHost)); 

    for (int i = 0; i < N; i++) if (h_Result_CPU[i] != h_Result_GPU[i]) {printf("mismatch2 at %d, cpu: %d, gpu %d\n", i, h_Result_CPU[i], h_Result_GPU[i]); return 1;} 

    printf("Test caching passed\n"); 

    return 0; 
}