Плохие результаты с простым FFT 1D OpenCL пример

Question

Я пытаюсь сделать простой БПФ 1D с clAmdFft OpenCL libray. Для этого я использую косинус-сигнал с частотой, равной 10 Гц, и частотой дискретизации sizex * frequency_signal с sizex количеством точек выборки.

В мои результаты я не могу получить импульс дирака при f = 10 Гц. В течение первых частот (K * f_sampling/SizeX с к = 0 до SizeX), я следующие выходы (колонка 1: к * f_sampling/размер х и столбцов 2: X_k)

0.000000 -5.169492e-06 
10.000000 -5.169449e-06 
20.000000 2.498745e-02 
30.000000 2.499508e-02 
40.000000 2.832322e-06 
50.000000 9.587315e-07 
60.000000 4.648561e-07 
70.000000 9.241408e-08 
80.000000 1.400218e-07 
90.000000 3.699876e-07 
100.000000 2.663564e-07 
110.000000 2.523218e-07 
120.000000 1.631196e-07 
130.000000 7.783362e-08 
140.000000 7.932793e-08 
150.000000 9.296434e-08 
160.000000 1.039581e-07 
170.000000 7.396823e-08 
180.000000 6.584698e-08 
190.000000 1.468647e-08 
200.000000 1.694581e-08 
210.000000 3.367836e-08 
220.000000 3.538253e-08 
... 

Как видно , кажется, что есть дирак (я говорю, что из-за значения «e-02» по сравнению с другими значениями, но я не уверен) между f = 20 Гц и f = 30 Гц, а не для 10 Гц, как ожидалось.

Вот код:

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

///////////////////////////////////// 
// OpenCL FFT 1D function /////////// 
///////////////////////////////////// 

int FftOpenCL(float *tab, const char* direction, int Ng) 
{ 
cl_int err; 
    cl_platform_id platform = 0; 
    cl_device_id device = 0; 
    cl_context_properties props[3] = { CL_CONTEXT_PLATFORM, 0, 0 }; 
    cl_context ctx = 0; 
    cl_command_queue queue = 0; 
    cl_mem bufX; 
    int ret = 0; 
    size_t N = Ng/2; 

    /* FFT library realted declarations */ 
    clAmdFftPlanHandle planHandle; 
    clAmdFftDim dim = CLFFT_1D; 
    size_t clLengths[1] = {N}; 

    /* Setup OpenCL environment. */ 
    err = clGetPlatformIDs(1, &platform, NULL); 
    err = clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL); 

    props[1] = (cl_context_properties)platform; 
    ctx = clCreateContext(props, 1, &device, NULL, NULL, &err); 
    queue = clCreateCommandQueue(ctx, device, 0, &err); 

    /* Setup clFFT. */ 
    clAmdFftSetupData fftSetup; 
    err = clAmdFftInitSetupData(&fftSetup); 
    err = clAmdFftSetup(&fftSetup); 

    /* Prepare OpenCL memory objects and place data inside them. */ 
    bufX = clCreateBuffer(ctx, CL_MEM_READ_WRITE, N * 2 * sizeof(float), NULL, &err); 

    err = clEnqueueWriteBuffer(queue, bufX, CL_TRUE, 0, 
    N * 2 * sizeof(float), tab, 0, NULL, NULL); 

    /* Create a default plan for a complex FFT. */ 
    err = clAmdFftCreateDefaultPlan(&planHandle, ctx, dim, clLengths); 

    /* Set plan parameters. */ 
    err = clAmdFftSetPlanPrecision(planHandle, CLFFT_SINGLE); 
    err = clAmdFftSetLayout(planHandle, CLFFT_COMPLEX_INTERLEAVED, CLFFT_COMPLEX_INTERLEAVED); 
    err = clAmdFftSetResultLocation(planHandle, CLFFT_INPLACE); 

    /* Bake the plan. */ 
    err = clAmdFftBakePlan(planHandle, 1, &queue, NULL, NULL); 


    if(strcmp(direction,"forward") == 0) 
    { 
    /* Execute the plan. */ 
    err = clAmdFftEnqueueTransform(planHandle, CLFFT_FORWARD, 1, &queue, 0, NULL, NULL, &bufX, NULL, NULL); 
    } 

    /* Wait for calculations to be finished. */ 
    err = clFinish(queue); 

    /* Fetch results of calculations. */ 
    err = clEnqueueReadBuffer(queue, bufX, CL_TRUE, 0, N * 2 * sizeof(float), tab, 0, NULL, NULL); 

    /* Release OpenCL memory objects. */ 
    clReleaseMemObject(bufX); 

    /* Release the plan. */ 
    err = clAmdFftDestroyPlan(&planHandle); 

    /* Release clFFT library. */ 
    clAmdFftTeardown(); 

    /* Release OpenCL working objects. */ 
    clReleaseCommandQueue(queue); 
    clReleaseContext(ctx); 

    return ret; 

} 


int main(void) { 

    int i; 

    // temporal array 
    float *x_temp; 

    // signal array and FFT output array 
    float *Array; 

    // number of sampling points 
    int sizex = 10000; 

    float h = 0; 

    // signal frequency 
    float frequency_signal = 10; 

    // size x points between 0 and T_signal 
    float frequency_sampling = sizex*frequency_signal; 

    // step = T_sampling 
    float step = 1.0/frequency_sampling; 

    FILE *Fft1D; 

    // Allocation of Array 

    Array = (float*) malloc(sizex*sizeof(float)); 
    x_temp = (float*) malloc(sizex*sizeof(float)); 

    // Initialization of 1D ArrayInput 

    for(i=0; i<sizex; i++) 
    { 
     Array[i] = cos(2*M_PI*frequency_signal*h); 
     x_temp[i] = h; 
     h = h + step; 
    } 

    if (FftOpenCL(Array,"forward", sizex) == 0) 
     printf("Fft passed !\n");; 

    // Printf ArrayOutput - Multiply FFT Array result by period T_sampling 

    Fft1D = fopen("Fft1D.dat","w"); 

    for(i=0; i<sizex; i++) 
     fprintf(Fft1D,"%f %e\n", i*frequency_sampling/sizex, 1.0/frequency_sampling*Array[i]); 

    fclose(Fft1D); 

    return 0; 

} 

Любой мог видеть, что это не так?

Благодаря

Answer 1

Похоже, что существует несоответствие типа между реальными и сложными массивами. Предполагая, что COMPLEX_INTERLEAVED означает пары (п, им) значений, вам нужно будет выделить массиву 2 * SizeX и инициализации:

array[2*i ] = cos(...); 
array[2*i+1] = 0.0f; 

Или изменить тип входа к реальным.